Fix correlation against phase ref

Also remove the old autocorrelation scripts, not pursuing this approach

9 files changed, 75 insertions, 147 deletions

diff --git a/LICENSE b/LICENSE
index 0acb70d..186a6f0 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2015 Matthias P. Braendli
+Copyright (c) 2016 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
diff --git a/README.md b/README.md
index 783a224..0d24a56 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
-Autocorrelation scripts
-=======================
+CIR scripts
+===========
 
-This is a set of script to do autocorrelation measurements
+This is a set of script to do correlation measurement measurements
 of a DAB signal. The goal is to one day do a channel impulse
 response measurement.
 
@@ -11,16 +11,14 @@ Right now there are three scripts:
   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.
+* correlate_with_ref.py: Finds the NULL symbol of an IQ file, and runs
+  correlations against the known phase reference to find the components
 
-* 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 correlation](./example_corr.png)
 
-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 example image shows the correlation (horizontal: time delay, vertical:
+time offset of the different pieces) between a simulated channel
+and the phase reference.
 The peaks at delay indices 0, 14 and 25 are clearly visible.
 
 Next steps
@@ -32,7 +30,7 @@ Requirements
 ------------
 Python with NumPy and matplotlib
 
-The iq files must be complex float.
+The iq files must be complex float or interleaved unsigned 8-bit.
 
 
 Licence
diff --git a/autocorrelate.py b/autocorrelate.py
deleted file mode 100755
index 9a66917..0000000
--- a/autocorrelate.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#!/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
deleted file mode 100755
index cf69a77..0000000
--- a/autocorrelate_window.py
+++ /dev/null
@@ -1,77 +0,0 @@
-#!/usr/bin/env python
-#
-# Do a set of autocorrelations over the test.16.14.25.iq file
-#
-# Licence: see LICENCE file
-#
-# hackrf_transfer example:
-# hackrf_transfer -r hackrf_dab.iq -f 211648000 -s 8 -n 768000
-
-import numpy as np
-import matplotlib.pyplot as pp
-import sys
-
-if len(sys.argv) != 3:
-    print("Usage")
-    print(" script [f64|u8] <filename>")
-    print(" fc64: file is 32-bit float I + 32-bit float Q")
-    print(" u8:   file is 8-bit unsigned I + 8-bit unsigned Q")
-    sys.exit(1)
-
-file_in = sys.argv[2]
-
-if sys.argv[1] == "u8":
-    channel_out1 = np.fromfile(file_in, np.uint8)
-    print("Convert u8 IQ to fc64 IQ")
-    channel_out2 = channel_out1.reshape(2, len(channel_out1)/2)
-    channel_out3 = channel_out2[0,...] + 1j * channel_out2[1,...]
-    channel_out = channel_out3.astype(np.complex64) / 256.0 - (0.5+0.5j)
-elif sys.argv[1] == "fc64":
-    channel_out = np.fromfile(file_in, np.complex64)
-
-channel_out = channel_out[:channel_out.size / 4]
-
-print("Autocorrelating")
-
-correlationlength = 500
-
-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 - (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))
-
-num_windows = len(channel_out_reshaped)
-
-for i, window in enumerate(channel_out_reshaped):
-    if i % 100 == 0:
-        print("Window {}/{}".format(i, num_windows))
-    channel_autocorr_image[i] = autocorrelate(window, correlationlength)
-
-rows, cols = channel_autocorr_image.shape
-
-print("Shape: {}x{}".format(rows, cols))
-
-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))
-
-fig2 = pp.figure()
-ax = fig2.add_subplot(211)
-accumulated0 = channel_autocorr_image.sum(axis=0)
-ac1 = ax.plot(accumulated0)
-
-ax = fig2.add_subplot(212)
-accumulated1 = channel_autocorr_image.sum(axis=1)
-ac1 = ax.plot(accumulated1)
-
-pp.show()
-
diff --git a/correlate_with_ref.py b/correlate_with_ref.py
index 508206a..de68c9d 100755
--- a/correlate_with_ref.py
+++ b/correlate_with_ref.py
@@ -10,7 +10,7 @@ import sys
 
 if len(sys.argv) != 3:
     print("Usage")
-    print(" script [f64|u8] <filename>")
+    print(" script [fc64|u8] <filename>")
     print(" fc64: file is 32-bit float I + 32-bit float Q")
     print(" u8:   file is 8-bit unsigned I + 8-bit unsigned Q")
     sys.exit(1)
@@ -22,37 +22,75 @@ file_in = sys.argv[2]
 if sys.argv[1] == "u8":
     channel_out1 = np.fromfile(file_in, np.uint8)
     print("Convert u8 IQ to fc64 IQ")
-    channel_out2 = channel_out1.reshape(2, len(channel_out1)/2)
+    channel_out2 = channel_out1.reshape(2, int(len(channel_out1)/2))
     channel_out3 = channel_out2[0,...] + 1j * channel_out2[1,...]
     channel_out = channel_out3.astype(np.complex64) / 256.0 - (0.5+0.5j)
 elif sys.argv[1] == "fc64":
     channel_out = np.fromfile(file_in, np.complex64)
 
-channel_out = channel_out[0:channel_out.size/2]
+print("  File contains {} samples ({}ms)".format(
+    len(channel_out), len(channel_out) / 2048000.0))
 
-print("Reading phase reference")
+# T = 1/2048000 s
+# NULL symbol is 2656 T (about 1.3ms) long.
+T_NULL = 2656
+# Full transmission frame in TM1 is 96ms = 196608 T.
+T_TF = 196608
 
+print("Reading phase reference")
 phase_ref = np.fromfile("phasereference.2048000.fc64.iq", np.complex64)
 
+# As we do not want to correlate of the whole recording that might be
+# containing several transmission frames, we first look for the null symbol in the
+# first 96ms
+
+print("Searching for NULL symbol")
+
+# Calculate power on blocks of length 2656 over the first 96ms. To gain speed,
+# we move the blocks by 10 samples.
+
+channel_out_power = np.array([np.abs(channel_out[t:t+T_NULL]).sum() for t in range(0, T_TF-T_NULL, 10)])
+
+# Look where the power is smallest, this gives the index where the NULL starts.
+# Because if the subsampling, we need to multiply the index.
+
+t_null = 10 * channel_out_power.argmin()
+
+print("  NULL symbol starts at ix={}".format(t_null))
+
+# The synchronisation channel occupies 5208 T and contains NULL symbol and
+# phase reference symbol. The phase reference symbol is 5208 - 2656 = 2552 T
+# long.
+
+if len(phase_ref) != 2552:
+    print("Warning: phase ref len is {} != 2552".format(len(phase_ref)))
+
+
+# We want to correlate our known phase reference symbol against the received
+# signal, and give us some more margin about the exact position of the NULL
+# symbol.
 print("Correlating")
 
-num_correlations = channel_out.size - phase_ref.size
-print("{} correlations to do...".format(num_correlations))
+# We start a bit earlier than the end of the null symbol
+corr_start_ix = t_null + T_NULL - 50
 
-correlations = np.array([np.abs(np.corrcoef(channel_out[i:phase_ref.size + i], phase_ref)[0,1]) for i in range(num_correlations)])
+# In TM1, the longest spacing between carrier components one can allow is
+# around 504 T (246us, or74km at speed of light). This gives us a limit
+# on the number of correlations it makes sense to do.
 
-print("Done")
+max_component_delay = 1000 # T
 
-numpeaks = 6
-print("The first {} highest peaks are at".format(numpeaks))
-print("  index: amplitude")
-for ind in correlations.argsort()[-numpeaks:][::-1]:
-    print("   {:4}: {}".format(ind, correlations[ind]))
+cir = np.array([np.abs(np.corrcoef(channel_out[corr_start_ix + i:corr_start_ix + phase_ref.size + i], phase_ref)[0,1]) for i in range(max_component_delay)])
+
+print("Done")
 
 fig = pp.figure()
 ax = fig.add_subplot(111)
-hi = ax.plot(correlations)
+hi = ax.plot(channel_out_power)
 
+fig = pp.figure()
+ax = fig.add_subplot(111)
+hi = ax.plot(cir)
 
 pp.show()
 
diff --git a/example_autocorr.png b/example_autocorr.png
deleted file mode 100644
index 19cad54..0000000
--- a/example_autocorr.png
+++ /dev/null
diff --git a/example_corr.png b/example_corr.png
new file mode 100644
index 0000000..810bb26
--- /dev/null
+++ b/example_corr.png
diff --git a/simulate_channel.py b/simulate_channel.py
index 893d74b..dc0ffb2 100755
--- a/simulate_channel.py
+++ b/simulate_channel.py
@@ -1,15 +1,22 @@
 #!/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
+# Adds additional components to the input file, which is
+# an array of complex floats, and save it as a new file
 #
 # Licence: see LICENCE file
 
 import numpy as np
 import matplotlib.pyplot as pp
+import sys
 
-file_in = "test.16.iq"
+if len(sys.argv) != 3:
+    print("Usage")
+    print(" script <filename_input> <filename_output>")
+    print(" input and output files are fc64 format")
+    sys.exit(1)
+
+file_in = sys.argv[1]
+file_out = sys.argv[2]
 
 
 # Simulate the channel with several components,
@@ -37,5 +44,5 @@ for delay, ampl in CIR:
 
 # The simulated channel output is in channel_out
 
-channel_out.tofile("test.16.14.25.iq")
+channel_out.tofile(file_out)
 
diff --git a/test.16.iq b/test.16.iq
deleted file mode 100644
index e3aeb3f..0000000
--- a/test.16.iq
+++ /dev/null