From d81d2a564752c8cf085d2cdff12c0d0fdbaba102 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 1 May 2015 17:19:35 +0200 Subject: Add first code experiments --- README.md | 41 +++++++++++++++++++++++++++++++++++++++++ autocorrelate.py | 38 ++++++++++++++++++++++++++++++++++++++ autocorrelate_window.py | 43 +++++++++++++++++++++++++++++++++++++++++++ example_autocorr.png | Bin 0 -> 168763 bytes simulate_channel.py | 41 +++++++++++++++++++++++++++++++++++++++++ test.16.iq | Bin 0 -> 6286144 bytes 6 files changed, 163 insertions(+) create mode 100644 README.md create mode 100755 autocorrelate.py create mode 100755 autocorrelate_window.py create mode 100644 example_autocorr.png create mode 100755 simulate_channel.py create mode 100644 test.16.iq 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 Binary files /dev/null and b/example_autocorr.png differ 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 Binary files /dev/null and b/test.16.iq differ -- cgit v1.2.3