Add Amp characterization in sync-measurement.ipynb

2 files changed, 135 insertions, 2 deletions

diff --git a/src/dab_util.py b/src/dab_util.py
index 2ad2af6..aa2030f 100644
--- a/src/dab_util.py
+++ b/src/dab_util.py
@@ -25,7 +25,7 @@ def freq_to_fft_sample(freq, fft_size, sampling_rate):
     freq_ratio = 1.0 * fft_size / sampling_rate
     return int(freq * freq_ratio + fft_size / 2)
 
-def crop_signal(signal, n_window = 1000, n_zeros = 1000, debug = False):
+def crop_signal(signal, n_window = 1000, n_zeros = 120000, debug = False):
     #signal = signal[-10:-1]
     mag = abs(signal)
     window = np.ones(n_window) / float(n_window)
@@ -38,5 +38,70 @@ def crop_signal(signal, n_window = 1000, n_zeros = 1000, debug = False):
         plt.plot(mag < thr)
         plt.plot((idx_start,idx_start), (0,0.1), color='g', linewidth=2)
         plt.plot((idx_end,idx_end), (0,0.1), color='r', linewidth=2)
-    signal = signal[idx_start - n_zeros: idx_end + n_zeros]
+    signal = signal[max(0,idx_start - n_zeros): min(idx_end + n_zeros, signal.shape[0] -1)]
     return signal
+
+def lagcorr(x,y,lag=None,verbose=False):
+    '''Compute lead-lag correlations between 2 time series.
+
+    <x>,<y>: 1-D time series.
+    <lag>: lag option, could take different forms of <lag>:
+          if 0 or None, compute ordinary correlation and p-value;
+          if positive integer, compute lagged correlation with lag
+          upto <lag>;
+          if negative integer, compute lead correlation with lead
+          upto <-lag>;
+          if pass in an list or tuple or array of integers, compute
+          lead/lag correlations at different leads/lags.
+
+    Note: when talking about lead/lag, uses <y> as a reference.
+    Therefore positive lag means <x> lags <y> by <lag>, computation is
+    done by shifting <x> to the left hand side by <lag> with respect to
+    <y>.
+    Similarly negative lag means <x> leads <y> by <lag>, computation is
+    done by shifting <x> to the right hand side by <lag> with respect to
+    <y>.
+
+    Return <result>: a (n*2) array, with 1st column the correlation
+    coefficients, 2nd column correpsonding p values.
+
+    Currently only works for 1-D arrays.
+    '''
+
+    import numpy
+    from scipy.stats import pearsonr
+
+    if len(x)!=len(y):
+        raise Exception('Input variables of different lengths.')
+
+    #--------Unify types of <lag>-------------
+    if numpy.isscalar(lag):
+        if abs(lag)>=len(x):
+            raise Exception('Maximum lag equal or larger than array.')
+        if lag<0:
+            lag=-numpy.arange(abs(lag)+1)
+        elif lag==0:
+            lag=[0,]
+        else:
+            lag=numpy.arange(lag+1)
+    elif lag is None:
+        lag=[0,]
+    else:
+        lag=numpy.asarray(lag)
+
+    #-------Loop over lags---------------------
+    result=[]
+    if verbose:
+        print '\n#<lagcorr>: Computing lagged-correlations at lags:',lag
+
+    for ii in lag:
+        if ii<0:
+            result.append(pearsonr(x[:ii],y[-ii:]))
+        elif ii==0:
+            result.append(pearsonr(x,y))
+        elif ii>0:
+            result.append(pearsonr(x[ii:],y[:-ii]))
+
+    result=numpy.asarray(result)
+
+    return result
diff --git a/src/signal_gen.py b/src/signal_gen.py
new file mode 100644
index 0000000..66306f3
--- /dev/null
+++ b/src/signal_gen.py
@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def gen_single_tone_with_silence(path = "./input.dat", n_periods=100, predist = None, par = None, debug = False):
+    period = 128
+    length = period * n_periods * 3
+
+    t = np.arange(0,length / 3)
+    sig = np.zeros(length, dtype=np.complex64)
+    sig[length/3:length*2/3] = np.exp(t * 2j * np.pi * 1./period)
+    #sin1 = np.cos(t * np.pi * 1./period)
+
+    if predist is None:
+        res = sig
+    else:
+        res = predist(sig, par)
+
+    res = res / np.max(res) * 0.95
+
+    res = res.astype(np.complex64)
+    res.tofile(path)
+
+    a_load = np.fromfile(path, dtype=np.complex64)
+    assert(np.isclose(a_load, res).all()), "Inconsistent stored file"
+
+    if debug == True:
+        plt.plot(np.abs(np.concatenate((a_load, a_load))))
+        plt.savefig(path + ".png")
+        plt.clf()
+
+        plt.plot(np.abs(np.fft.fftshift(np.fft.fft(np.concatenate((a_load, a_load))))), 'ro')
+        plt.savefig(path + "_fft.png")
+        plt.clf()
+
+    return path
+
+def gen_ramps(path = "./input.dat", n_periods=128, pause = 64, amplitudes = [1], predist = None, par = None, debug = False):
+    period = 128
+    l_pause = period * pause
+    l_sig = period * n_periods
+    length = l_pause + l_sig
+
+    t = np.arange(0, l_sig)
+    sig = np.zeros(length * len(amplitudes), dtype=np.complex64)
+    for i, amplitude in enumerate(amplitudes):
+        sig[length * i + l_pause : length * (i + 1)] = complex(amplitude) * np.exp(t * 2j * np.pi * 1./period)
+
+    if predist is None:
+        res = sig
+    else:
+        res = predist(sig, par)
+
+    res = res.astype(np.complex64)
+    res.tofile(path)
+
+    a_load = np.fromfile(path, dtype=np.complex64)
+    assert(np.isclose(a_load, res).all()), "Inconsistent stored file"
+
+    if debug == True:
+        plt.plot(np.abs(np.concatenate((a_load, a_load))))
+        plt.savefig(path + ".png")
+        plt.clf()
+
+        plt.plot(np.abs(np.fft.fftshift(np.fft.fft(np.concatenate((a_load, a_load))))), 'ro')
+        plt.savefig(path + "_fft.png")
+        plt.clf()
+
+    return path