Fix lag calculator and add test

3 files changed, 71 insertions, 10 deletions

diff --git a/src/dab_util.py b/src/dab_util.py
index 2b23812..843f8a5 100644
--- a/src/dab_util.py
+++ b/src/dab_util.py
@@ -3,6 +3,7 @@ import scipy
 import matplotlib.pyplot as plt
 import fftconvolve
 import src.dabconst as dabconst
+from scipy import signal
 
 c = {}
 c["bw"]=1536000
@@ -47,27 +48,45 @@ def crop_signal(signal, n_window = 1000, n_zeros = 120000, debug = False):
     signal = signal[max(0,idx_start - n_zeros): min(idx_end + n_zeros, signal.shape[0] -1)]
     return signal
 
-#def fftlag(signal_original, signal_rec):
+#def fftlag(sig_orig, sig_rec):
 #    """
 #    Efficient way to find lag between two signals
 #    Args:
-#        signal_original: The signal that has been sent
-#        signal_rec: The signal that has been recored
+#        sig_orig: The signal that has been sent
+#        sig_rec: The signal that has been recored
 #    """
-#    c = np.flipud(scipy.signal.fftconvolve(signal_original,np.flipud(signal_rec)))
+#    c = np.flipud(scipy.signal.fftconvolve(sig_orig,np.flipud(sig_rec)))
 #    #plt.plot(c)
-#    return np.argmax(c) - signal_original.shape[0] + 1
+#    return np.argmax(c) - sig_orig.shape[0] + 1
 
-def fftlag(signal_original, signal_rec, n_upsampling = 1):
+def lag(sig_orig, sig_rec):
+    """
+    Find lag between two signals
+    Args:
+        sig_orig: The signal that has been sent
+        sig_rec: The signal that has been recored
+    """
+    off = sig_rec.shape[0]
+    c = signal.correlate(sig_orig, sig_rec)
+    return np.argmax(c) - off + 1
+
+def lag_upsampling(sig_orig, sig_rec, n_up):
+    sig_orig_up = signal.resample(sig_orig, sig_orig.shape[0] * n_up)
+    sig_rec_up  = signal.resample(sig_rec, sig_rec.shape[0] * n_up)
+    l = lag(sig_orig_up, sig_rec_up)
+    l_orig = float(l) / n_up
+    return l_orig
+
+def fftlag(sig_orig, sig_rec, n_upsampling = 1):
     """
     Efficient way to find lag between two signals
     Args:
-        signal_original: The signal that has been sent
-        signal_rec: The signal that has been recored
+        sig_orig: The signal that has been sent
+        sig_rec: The signal that has been recored
     """
-    c = np.flipud(fftconvolve.fftconvolve(signal_original,np.flipud(signal_rec), n_upsampling))
+    c = np.flipud(fftconvolve.fftconvolve(sig_orig,np.flipud(sig_rec), n_upsampling))
     #plt.plot(c)
-    return (np.argmax(c) - signal_original.shape[0] + 1)
+    return (np.argmax(c) - sig_orig.shape[0] + 1)
 
 def get_amp_ratio(ampl_1, ampl_2, a_out_abs, a_in_abs):
     idxs = (a_in_abs > ampl_1) & (a_in_abs < ampl_2)
diff --git a/src/dab_util_test.py b/src/dab_util_test.py
new file mode 100644
index 0000000..83813a9
--- /dev/null
+++ b/src/dab_util_test.py
@@ -0,0 +1,36 @@
+from scipy import signal
+import numpy as np
+import src.gen_source as gs
+reload(gs)
+import src.dab_util as du
+reload(du)
+
+def gen_test_signals(oversampling=4, sample_offset_float=0):
+    off = int(sample_offset_float)
+    phi_samples = sample_offset_float - off
+    phi = phi_samples*360/oversampling
+
+    s1 = np.zeros((1024))
+    s1[256:768] = gs.gen_sin(512, oversampling, 0)
+    s2 = np.zeros((1024))
+    s2[256+off:768+off] = gs.gen_sin(512, oversampling, phi)
+
+    return s1, s2
+
+def test_phase_offset(lag_function, tol):
+    def r():
+        return np.random.rand(1)*100-50
+    res = []
+    for i in range(100):
+        off = r()
+        s1, s2 = gen_test_signals(
+            oversampling=4, sample_offset_float=off)
+
+        off_meas = lag_function(s2, s1)
+        res.append(np.abs(off-off_meas)<tol)
+    return np.mean(res)
+
+for n_up in [1, 2, 3, 4, 7, 8, 16]:
+   correct_ratio = test_phase_offset(lambda x,y: du.lag_upsampling(x,y,n_up), tol=1./n_up)
+   print("%.1f%% of the tested offsets were measured within tolerance %.4f for n_up = %d" % (correct_ratio * 100, 1./n_up, n_up))
+
diff --git a/src/gen_source.py b/src/gen_source.py
index 7620bc3..c8509ed 100644
--- a/src/gen_source.py
+++ b/src/gen_source.py
@@ -97,3 +97,9 @@ def gen_file_i(frequency_0, frequency_1, x1 = 0, x2 = 0, x3 = 0, x4 = 0, samp_ra
     assert(np.isclose(a_load, two_tone).all()), "Inconsistent stored file"
 
     return path
+
+def gen_sin(samples, oversampling, phi):
+    t = np.arange(samples, dtype=np.float)
+    sig = np.sin(((2*np.pi)/oversampling) * t - np.pi*phi/180.)
+    return sig
+