diff --git a/src/noise_module.py b/src/noise_module.py index d038020c..e0fb4a9e 100644 --- a/src/noise_module.py +++ b/src/noise_module.py @@ -1791,6 +1791,93 @@ def stretching(ref, cur, dv_range, nbtrial, para): return dv, error, cc, cdp +def stretching_vect(ref, cur, dv_range, nbtrial, para): + + """ + This function compares the Reference waveform to stretched/compressed current waveforms to get the relative seismic velocity variation (and associated error). + It also computes the correlation coefficient between the Reference waveform and the current waveform. + + PARAMETERS: + ---------------- + ref: Reference waveform (np.ndarray, size N) + cur: Current waveform (np.ndarray, size N) + dv_range: absolute bound for the velocity variation; example: dv=0.03 for [-3,3]% of relative velocity change ('float') + nbtrial: number of stretching coefficient between dvmin and dvmax, no need to be higher than 100 ('float') + para: vector of the indices of the cur and ref windows on wich you want to do the measurements (np.ndarray, size tmin*delta:tmax*delta) + For error computation, we need parameters: + fmin: minimum frequency of the data + fmax: maximum frequency of the data + tmin: minimum time window where the dv/v is computed + tmax: maximum time window where the dv/v is computed + RETURNS: + ---------------- + dv: Relative velocity change dv/v (in %) + cc: correlation coefficient between the reference waveform and the best stretched/compressed current waveform + cdp: correlation coefficient between the reference waveform and the initial current waveform + error: Errors in the dv/v measurements based on Weaver et al (2011), On the precision of noise-correlation interferometry, Geophys. J. Int., 185(3) + + Note: The code first finds the best correlation coefficient between the Reference waveform and the stretched/compressed current waveform among the "nbtrial" values. + A refined analysis is then performed around this value to obtain a more precise dv/v measurement . + + Originally by L. Viens 04/26/2018 (Viens et al., 2018 JGR) + modified by Chengxin Jiang + modified by Laura Ermert: vectorized version + """ + # load common variables from dictionary + twin = para['twin'] + freq = para['freq'] + dt = para['dt'] + tmin = np.min(twin) + tmax = np.max(twin) + fmin = np.min(freq) + fmax = np.max(freq) + tvec = np.arange(tmin, tmax, dt) + + # make useful one for measurements + dvmin = -np.abs(dv_range) + dvmax = np.abs(dv_range) + Eps = 1 + (np.linspace(dvmin, dvmax, nbtrial)) + cdp = np.corrcoef(cur, ref)[0, 1] # correlation coefficient between the reference and initial current waveforms + waveforms = np.zeros((nbtrial + 1, len(ref))) + waveforms[0, :] = ref + + # Set of stretched/compressed current waveforms + for ii in range(nbtrial): + nt = tvec * Eps[ii] + s = np.interp(x=tvec, xp=nt, fp=cur) + waveforms[ii + 1, :] = s + cof = np.corrcoef(waveforms)[0][1:] + + # find the maximum correlation coefficient + imax = np.nanargmax(cof) + if imax >= len(Eps)-2: + imax = imax - 2 + if imax < 2: + imax = imax + 2 + + # Proceed to the second step to get a more precise dv/v measurement + dtfiner = np.linspace(Eps[imax-2], Eps[imax+2], nbtrial) + #ncof = np.zeros(dtfiner.shape,dtype=np.float32) + waveforms = np.zeros((nbtrial + 1, len(ref))) + waveforms[0, :] = ref + for ii in range(len(dtfiner)): + nt = tvec * dtfiner[ii] + s = np.interp(x=tvec, xp=nt, fp=cur) + waveforms[ii + 1, :] = s + ncof = np.corrcoef(waveforms)[0][1: ] + cc = np.max(ncof) # Find maximum correlation coefficient of the refined analysis + dv = 100. * dtfiner[np.argmax(ncof)] - 100 # Multiply by 100 to convert to percentage (Epsilon = -dt/t = dv/v) + + # Error computation based on Weaver et al (2011), On the precision of noise-correlation interferometry, Geophys. J. Int., 185(3) + T = 1 / (fmax - fmin) + X = cc + wc = np.pi * (fmin + fmax) + t1 = np.min([tmin, tmax]) + t2 = np.max([tmin, tmax]) + error = 100*(np.sqrt(1-X**2)/(2*X)*np.sqrt((6* np.sqrt(np.pi/2)*T)/(wc**2*(t2**3-t1**3)))) + + return dv, error, cc, cdp + def dtw_dvv(ref, cur, para, maxLag, b, direction): """ Dynamic time warping for dv/v estimation. diff --git a/test/test_routines/test_stretching.py b/test/test_routines/test_stretching.py new file mode 100644 index 00000000..f36e2057 --- /dev/null +++ b/test/test_routines/test_stretching.py @@ -0,0 +1,58 @@ +import numpy as np +from src.noise_module import stretching_vect, stretching +from obspy.signal.invsim import cosine_taper +import pytest + +# This short script is intended as a test for the stretching routine +# it takes a generic sine curve with known stretching factor and ensures +# that the stretching routines in noise_module always recover this factor + +def test_stretching(): + t = np.linspace(0., 9.95, 2500) # 0.5 % perturbation + original_signal = np.sin(t * 10.) * cosine_taper(2500, p=0.75) + + t_stretch = np.linspace(0., 10.0, 2500) + stretched_signal = np.interp(t, t_stretch, original_signal) + + + para = {} + para["dt"] = 1. / 250. + para["twin"] = [0., 10.] + para["freq"] = [9.9, 10.1] + + dvv, error, cc, cdp = stretching(ref=original_signal, cur=stretched_signal, + dv_range=0.05, nbtrial=100, para=para) + + assert pytest.approx(cc) == 1.0 + assert dvv + 0.5 < para["dt"] # assert result is -0.5% + +def test_stretching_vect(): + t = np.linspace(0., 9.95, 2500) # 0.5 % perturbation + original_signal = np.sin(t * 10.) * cosine_taper(2500, p=0.75) + + t_stretch = np.linspace(0., 10.0, 2500) + stretched_signal = np.interp(t, t_stretch, original_signal) + + + para = {} + para["dt"] = 1. / 250. + para["twin"] = [0., 10.] + para["freq"] = [9.9, 10.1] + + dvv, error, cc, cdp = stretching_vect(ref=original_signal, cur=stretched_signal, + dv_range=0.05, nbtrial=100, para=para) + + assert pytest.approx(cc) == 1.0 + assert dvv + 0.5 < para["dt"] # assert result is -0.5% + +if __name__ == "__main__": + import sys + if sys.argv[1] == "normal": + for i in range(300): + test_stretching() + elif sys.argv[1] == "vect": + for i in range(300): + test_stretching_vect() + else: + raise ValueError("call with python test_stretching.py ,\ +where choice of stretching is \"normal\" or \"vect\".")