am-am am-pm diagram in syncmeasurement.ipynb

3 files changed, 96 insertions, 18 deletions

diff --git a/src/dab_tuning_lib.py b/src/dab_tuning_lib.py
new file mode 100644
index 0000000..8faafef
--- /dev/null
+++ b/src/dab_tuning_lib.py
@@ -0,0 +1,31 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import src.dab_util as du
+
+
+def calc_signal_sholder_ratio(fft, sampling_rate, debug = False, debug_path="", suffix=""):
+    fft_size = fft.shape[0]
+    n_sig = (du.freq_to_fft_sample(-du.c["bw"]/2., fft_size, sampling_rate),
+             du.freq_to_fft_sample( du.c["bw"]/2., fft_size, sampling_rate))
+    sig = np.mean(fft[n_sig[0]:n_sig[1]])
+
+    n_noise = (du.freq_to_fft_sample(-3000000., fft_size, sampling_rate),
+               du.freq_to_fft_sample(-2500000, fft_size, sampling_rate))
+    noise = np.mean(fft[n_noise[0]:n_noise[1]])
+
+    n_sholder = (du.freq_to_fft_sample(-1500000, fft_size, sampling_rate),
+               du.freq_to_fft_sample(-du.c["bw"]/2, fft_size, sampling_rate))
+    sholder = np.mean(fft[n_sholder[0]:n_sholder[1]])
+
+    score = -sholder
+
+    if debug == True:
+        print(n_sig, n_sholder, n_noise)
+        plt.plot(fft)
+        plt.plot((n_sig[0], n_sig[1]), (sig, sig), linewidth=5, color='g')
+        plt.plot((n_noise[0], n_noise[1]), (noise, noise), linewidth=5, color='r')
+        plt.plot((n_sholder[0], n_sholder[1]), (sholder, sholder), linewidth=5, color='y')
+        plt.savefig(debug_path + "/" + str(score) + suffix + ".png")
+        plt.clf()
+
+    return score
diff --git a/src/dab_util.py b/src/dab_util.py
new file mode 100644
index 0000000..2ad2af6
--- /dev/null
+++ b/src/dab_util.py
@@ -0,0 +1,42 @@
+import numpy as np
+import scipy
+import matplotlib.pyplot as plt
+
+c = {
+        "bw":1536000
+        }
+
+def calc_fft(signal, fft_size = 1024, sampling_rate = 1, plot = False):
+    """returns one numpy array for the frequencies and one for the corresponding fft"""
+    signal_spectrum = np.fft.fftshift(np.fft.fft(signal, fft_size))
+    freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate))
+    if plot == True:
+        plot_freq_spec(freqs, signal_spectrum)
+    return freqs, signal_spectrum
+
+def plot_freq_spec(freq, spec = None):
+    plt.figure(figsize=(10,5))
+    if spec == None:
+        plt.plot(freq)
+    else:
+        plt.plot(freq, np.abs(spec))
+
+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):
+    #signal = signal[-10:-1]
+    mag = abs(signal)
+    window = np.ones(n_window) / float(n_window)
+    mag = scipy.signal.convolve(window, mag)
+    mag = scipy.signal.convolve(window, mag)
+    thr = 0.05 * np.max(mag)
+    idx_start = np.argmax(mag > thr)
+    idx_end = mag.shape[0] - np.argmax(np.flipud(mag > thr))
+    if debug:
+        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]
+    return signal
diff --git a/src/two_tone_lib.py b/src/two_tone_lib.py
index a3c9675..beeb332 100644
--- a/src/two_tone_lib.py
+++ b/src/two_tone_lib.py
@@ -4,13 +4,15 @@ import matplotlib.pyplot as plt
 def gen_two_tone(path = "./input.dat", predist = None, par = None, debug = False):
     period1 = 3.875
     period2 = 4
-    t_both = 124
+    t_both = 124 * 10
     assert(t_both / period1 % 1 == 0)
     assert(t_both / period2 % 1 == 0)
 
     t = np.arange(0,t_both)
     sin1 = np.exp(t * 2j * np.pi * 1./period1)
     sin2 = np.exp(t * 2j * np.pi * 1./period2)
+    #sin1 = np.cos(t * np.pi * 1./period1)
+    #sin2 = np.cos(t * np.pi * 1./period2)
     sig = sin1 + sin2
 
     if predist is None:
@@ -18,7 +20,7 @@ def gen_two_tone(path = "./input.dat", predist = None, par = None, debug = False
     else:
         res = predist(sig, par)
 
-    res = res / np.max(res)
+    res = res / np.max(res) * 0.99
 
     res = res.astype(np.complex64)
     res.tofile(path)
@@ -43,33 +45,36 @@ def predist_poly(sig, coefs = []):
         res += sig * np.abs(sig)**(idx+1) * coef #+1 because first correction term is squared
     return res
 
-def analyse_power_spec(spec, debug = False, debug_path="", suffix=""):
+def analyse_power_spec(spec, threshold=40, debug = False, debug_path="", suffix=""):
     peak_1 = None
     peak_2 = None
-    spec_start = 4096
-    spec_end = 8192
-    first_peak = spec_start + 2048
-    second_peak = spec_start + 2114
     delta_freq = 66
+    first_peak = 4096 + 2048
+    second_peak = 4096 + 2114
+    spec_start = first_peak - delta_freq * 10
+    spec_end = second_peak + delta_freq * 10
     peak_other = []
     if debug: plt.plot(spec)
-    for x in [c * delta_freq + delta_freq//2 for c in range(spec_start//delta_freq)]:
-        start = spec_start + x 
+    #find peaks
+    for x in [c * delta_freq + delta_freq//2 for c in range((spec_end - spec_start)//delta_freq)]:
+        start = spec_start + x
         end = spec_start + x + delta_freq
         peak = spec[start:end].max()
-        if debug: plt.plot((start,end), (peak, peak))
         if start < first_peak and end > first_peak:
-            peak_1 = peak
-            if debug: plt.plot((start,end), (peak+1, peak+1))
+            peak_1 = (start, end, peak)
         elif start < second_peak and end > second_peak:
-            peak_2 = peak
-            if debug: plt.plot((start,end), (peak+1, peak+1))
+            peak_2 = (start, end, peak)
         else:
-            peak_other.append(peak)
-    mean_signal = (peak_1 + peak_2) / 2
-    mean_others = np.mean(peak_other)
-    score = mean_signal - mean_others
+            peak_other.append((start, end, peak))
+    mean_signal = (peak_1[2] + peak_2[2]) / 2
+    #peak_other = [[s,e,p] for s, e, p in peak_other if mean_signal - p < threshold]
+    meas = [mean_signal - p for s, e, p in peak_other]
+    score = np.min(meas)
     if debug:
+        plt.plot((peak_1[0],peak_1[1]), (peak_1[2], peak_1[2]), color='g', linewidth=2)
+        plt.plot((peak_2[0],peak_2[1]), (peak_2[2], peak_2[2]), color='g', linewidth=2)
+        for start, end, peak in peak_other:
+            plt.plot((start, end), (peak, peak), color='r', linewidth=2)
         plt.savefig(debug_path + "/" + str(score) + suffix + ".png")
         plt.clf()
     return score