Merge pull request #2 from lermert/master

numpy version for stretching

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.

test/test_routines/test_stretching.py

@@ -0,0 +1,66 @@
+import numpy as np
+import sys
+import os
+sys.path.append(os.getcwd())
+from noise_module import stretching_vect, stretching
+from obspy.signal.invsim import cosine_taper
+import pytest
+import time
+
+# 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
+# Note: The script has to be called from src/ directory, like
+# (in directory noisepy/src:)
+# python ../test/test_routines/test_stretching.py
+
+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__":
+    print("Running stretching...")
+    t = time.time()
+    for i in range(300):
+        test_stretching()
+    print("Done stretching, no errors, %4.2fs." %(time.time()-t))
+
+    print("Running stretching using numpy...")
+    t = time.time()
+    for i in range(300):
+        test_stretching_vect()
+    print("Done stretching, no errors, %4.2fs." %(time.time()-t))