Add first code experiments

6 files changed, 163 insertions, 0 deletions

diff --git a/README.md b/README.md
new file mode 100644
index 0000000..783a224
--- /dev/null
+++ b/README.md
@@ -0,0 +1,41 @@
+Autocorrelation scripts
+=======================
+
+This is a set of script to do autocorrelation measurements
+of a DAB signal. The goal is to one day do a channel impulse
+response measurement.
+
+Right now there are three scripts:
+
+* simulate_channel.py: Reads an I/Q file generated by ODR-DabMod and
+  adds some additional components, somehow equivalent to signal reflexions.
+  Each reflexion has a delay and an amplitude factor.
+
+* autocorrelate.py: Reads the I/Q file generated by the previous script
+  and calculate a set of autocorrelations on the whole signal.
+
+* autocorrelate_window.py: Same as above, but cut the file in pieces,
+  calculate the autocorrelation for each piece, and display as a 2D image.
+
+Example image: ![Example autocorrelation](./example_autocorr.png)
+
+The example image shows the autocorrelation (horizontal: time delay, vertical:
+time offset of the different pieces) of the simulated channel.
+The peaks at delay indices 0, 14 and 25 are clearly visible.
+
+Next steps
+----------
+Read in a recording using a RTL-SDR receiver.
+
+
+Requirements
+------------
+Python with NumPy and matplotlib
+
+The iq files must be complex float.
+
+
+Licence
+-------
+MIT. See LICENCE for details.
+
diff --git a/autocorrelate.py b/autocorrelate.py
new file mode 100755
index 0000000..9a66917
--- /dev/null
+++ b/autocorrelate.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+#
+# Do a single autocorrelation over the whole test.16.14.25.iq file
+#
+# Licence: see LICENCE file
+
+import numpy as np
+import matplotlib.pyplot as pp
+
+file_in = "test.16.14.25.iq"
+
+channel_out = np.fromfile(file_in, np.complex64)
+
+print("Autocorrelating")
+
+correlationlength = 50
+
+def autocorrelate(x, length):
+    return np.array([1] + [np.abs(np.corrcoef(x[:-i], x[i:])[0,1]) for i in range(1, length)])
+
+autocorr = autocorrelate(channel_out, correlationlength)
+
+print("Done")
+
+numpeaks = 6
+print("The first {} highest peaks are at".format(numpeaks))
+print("  index: amplitude")
+for ind in autocorr.argsort()[-numpeaks:][::-1]:
+    print("   {:4}: {}".format(ind, autocorr[ind]))
+
+fig = pp.figure()
+ax = fig.add_subplot(111)
+hi = ax.semilogy(autocorr)
+
+
+pp.show()
+
+
diff --git a/autocorrelate_window.py b/autocorrelate_window.py
new file mode 100755
index 0000000..c97f7fb
--- /dev/null
+++ b/autocorrelate_window.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+#
+# Do a set of autocorrelations over the test.16.14.25.iq file
+#
+# Licence: see LICENCE file
+
+import numpy as np
+import matplotlib.pyplot as pp
+
+file_in = "test.16.14.25.iq"
+
+channel_out = np.fromfile(file_in, np.complex64)
+
+print("Autocorrelating")
+
+correlationlength = 50
+
+def autocorrelate(x, length):
+    return np.array([1] + [np.abs(np.corrcoef(x[:-i], x[i:])[0,1]) for i in range(1, length)])
+
+
+reshape_width = correlationlength * 4
+
+channel_out_truncated = channel_out[:-(channel_out.size % reshape_width)]
+
+channel_out_reshaped = channel_out_truncated.reshape(channel_out_truncated.size / reshape_width, reshape_width)
+
+channel_autocorr_image = np.zeros((channel_out_reshaped.shape[0], correlationlength))
+for i, window in enumerate(channel_out_reshaped):
+    if i % 100 == 0:
+        print("Window {}".format(i))
+    channel_autocorr_image[i] = autocorrelate(window, correlationlength)
+
+rows, cols = channel_autocorr_image.shape
+
+aspect_ratio = 1.0
+
+fig = pp.figure()
+ax = fig.add_subplot(111)
+hi = ax.imshow(channel_autocorr_image, cmap='hot', aspect=aspect_ratio*(cols/rows))
+
+pp.show()
+
diff --git a/example_autocorr.png b/example_autocorr.png
new file mode 100644
index 0000000..19cad54
--- /dev/null
+++ b/example_autocorr.png
diff --git a/simulate_channel.py b/simulate_channel.py
new file mode 100755
index 0000000..893d74b
--- /dev/null
+++ b/simulate_channel.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python
+#
+# Adds additional components to the test.16.iq file, which is
+# an array of complex floats, and save it as
+# test.16.14.25.iq
+#
+# Licence: see LICENCE file
+
+import numpy as np
+import matplotlib.pyplot as pp
+
+file_in = "test.16.iq"
+
+
+# Simulate the channel with several components,
+# shifted by some number of samples and multiplied by
+# the amplitude
+#
+# The (0, 1.0) component is always present
+CIR = [(14, 0.4), (25, 0.3)]
+
+print("Simulate channel {}".format(CIR))
+
+maxdelay = max(delay for delay, ampl in CIR)
+
+
+original_iq = np.fromfile(file_in, np.complex64)
+
+original_iq_extended = np.append(original_iq, np.zeros(maxdelay, dtype=np.complex64))
+
+channel_out = original_iq_extended
+
+for delay, ampl in CIR:
+    print("Add component {}".format(delay))
+    channel_out += np.append(np.zeros(delay, dtype=np.complex64), ampl * original_iq_extended[:-delay])
+
+
+# The simulated channel output is in channel_out
+
+channel_out.tofile("test.16.14.25.iq")
+
diff --git a/test.16.iq b/test.16.iq
new file mode 100644
index 0000000..e3aeb3f
--- /dev/null
+++ b/test.16.iq