diff options
| -rw-r--r-- | LICENSE | 2 | ||||
| -rw-r--r-- | README.md | 22 | ||||
| -rwxr-xr-x | autocorrelate.py | 38 | ||||
| -rwxr-xr-x | autocorrelate_window.py | 77 | ||||
| -rwxr-xr-x | correlate_with_ref.py | 66 | ||||
| -rw-r--r-- | example_autocorr.png | bin | 168763 -> 0 bytes | |||
| -rw-r--r-- | example_corr.png | bin | 0 -> 25850 bytes | |||
| -rwxr-xr-x | simulate_channel.py | 17 | ||||
| -rw-r--r-- | test.16.iq | bin | 6286144 -> 0 bytes |
9 files changed, 75 insertions, 147 deletions
@@ -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 @@ -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 image:  - -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 Binary files differdeleted file mode 100644 index 19cad54..0000000 --- a/example_autocorr.png +++ /dev/null diff --git a/example_corr.png b/example_corr.png Binary files differnew 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 Binary files differdeleted file mode 100644 index e3aeb3f..0000000 --- a/test.16.iq +++ /dev/null |
