1 files changed, 152 insertions, 0 deletions
diff --git a/align/GenerateExampleTxRxIQ.py b/align/GenerateExampleTxRxIQ.py
new file mode 100755
index 0000000..dc8b572
--- /dev/null
+++ b/align/GenerateExampleTxRxIQ.py
@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+#
+# Generate an example RX and TX dataset, with a subsample delay and try to resolve it afterwards
+#
+#
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Matthias P. Braendli
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import matplotlib.pyplot as plt
+import scipy.signal
+import numpy as np
+
+## Configuration
+# whether to correct for delays larger than one sample
+# Not necessary unless you have delay larger than oversample/2
+do_integer_compensation = 0
+
+# by how much to oversample the signal before applying the delay
+oversample = 16
+
+# Add a delay of delay/oversample samples to the input signal
+delay = 7
+
+print("Expecting a delay of {} samples".format(delay/oversample))
+
+## Generate signal
+
+iq_file = "/home/bram/dab/aux/odr-dab-cir/phasereference.2048000.fc64.iq"
+
+iq_data = np.fromfile(iq_file, np.complex64)
+
+# oversampling the input signal doesn't make much of a difference
+phase_ref_iq = scipy.signal.resample(iq_data, 2 * len(iq_data))
+
+# make the signal periodic by duplicating the signal
+phase_ref_iq = np.concatenate((phase_ref_iq, phase_ref_iq))
+
+noise_iq = np.random.normal(scale = np.max(np.abs(phase_ref_iq)) * 0.02,
+                           size=len(phase_ref_iq))
+
+phase_ref_iq = phase_ref_iq + noise_iq
+
+# exp(-2i pi f) is the Fourier transform of a unity delay.
+# exp(2i pi f) is a negative delay.
+bin_frequencies = np.concatenate(
+    (np.linspace(0, 0.5, len(phase_ref_iq)/2, endpoint=False),
+    np.linspace(-0.5, 0, len(phase_ref_iq)/2, endpoint=False)))
+
+phase_ref_uc = scipy.signal.resample(phase_ref_iq, oversample * len(phase_ref_iq))
+
+
+phase_ref_uc_delayed = np.roll(phase_ref_uc, delay)
+
+phase_ref_delayed = scipy.signal.resample(phase_ref_uc_delayed, len(phase_ref_iq))
+
+## Integer delay
+corr_begin_ix = -32
+corr_end_ix = 32
+
+corr = [np.abs(np.corrcoef(phase_ref_delayed, np.roll(phase_ref_iq, i))[0,1]) 
+        for i in range(corr_begin_ix, corr_end_ix)]
+# TODO check for negative real correlation peak
+if do_integer_compensation:
+    delay_in = np.argmax(corr) + corr_begin_ix
+else:
+    delay_in = 0
+
+phase_ref_int_delay_removed = np.roll(phase_ref_delayed, -delay_in)
+
+print("Integer delay corrected: {}".format(delay_in))
+
+## Fractional delay
+signal_fft = np.fft.fft(phase_ref_int_delay_removed)
+reference_fft = np.fft.fft(phase_ref_iq)
+
+# rotate each bin backwards with the phase of the reference. As we have already resolved the
+# integer delay, we should find at most one 2*pi wrapping.
+u = signal_fft * np.conj(reference_fft)
+
+
+# the phase signal will still wrap around, and will have values between -pi/4 and pi/4
+phase_wrapping = np.angle(u)
+
+unwrap_with_deriv_integrate = True
+if unwrap_with_deriv_integrate:
+    # to unwrap, take the derivative, remove peaks, integrate
+    phase_deriv = phase_wrapping - np.roll(phase_wrapping, 1)
+
+    def filter_phase_deriv(p):
+        if np.abs(p) < 0.5:
+            return p
+        else:
+            return 0
+
+    phase_deriv_nopeaks = [filter_phase_deriv(p) for p in phase_deriv]
+    phase_unwrapped = np.cumsum(phase_deriv_nopeaks)
+else:
+    # doesn't always work, sometimes there are smaller jumps in phase
+    phase_unwrapped = np.mod(phase_wrapping + np.pi/2, np.ones(len(phase_wrapping)) * np.pi / 2)
+
+
+# Find the slope using a linear regression
+slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(phase_unwrapped,range(len(phase_unwrapped)))
+
+if p_value < 0.05:
+    frac_delay = slope / len(phase_unwrapped)
+    frac_rotate = intercept / len(phase_unwrapped)
+    print("Applying subsample correction: {} {}".format(frac_delay, frac_rotate))
+    print(slope, intercept, r_value, p_value, std_err)    
+else:
+    print("Skipping subsample correction")
+    print(slope, intercept, r_value, p_value, std_err)
+    frac_delay = None
+
+plt.figure()
+plt.plot(phase_wrapping)
+plt.plot(phase_unwrapped)
+
+
+if frac_delay:
+    fine_shift_fft = np.exp((0+2j * np.pi * frac_delay) * bin_frequencies) * np.exp(0+2j * np.pi * frac_rotate)
+    sig_delay_removed_fft = signal_fft * fine_shift_fft
+
+    sig_delay_removed = np.fft.ifft(sig_delay_removed_fft)
+
+    plt.figure()
+    plt.plot(np.angle(np.fft.fftshift(fine_shift_fft)))
+
+    plt.figure()
+    plt.plot(np.abs(np.fft.fftshift(np.fft.fft(phase_ref_iq))))
+    plt.plot(np.abs(np.fft.fftshift(np.fft.fft(sig_delay_removed-phase_ref_iq))))
+
+plt.show()