Cumulate measurements over all TF in file

1 files changed, 44 insertions, 34 deletions

diff --git a/correlate_with_ref.py b/correlate_with_ref.py
index de68c9d..193fa8c 100755
--- a/correlate_with_ref.py
+++ b/correlate_with_ref.py
@@ -40,58 +40,68 @@ 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
+def calc_cir(channel, start_ix):
+    """Calculate correlation with phase reference"""
 
-print("Searching for NULL symbol")
+    channel_out
 
-# 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)])
+    # 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")
 
-# Look where the power is smallest, this gives the index where the NULL starts.
-# Because if the subsampling, we need to multiply the index.
+    # Calculate power on blocks of length 2656 over the first 96ms. To gain speed,
+    # we move the blocks by N samples.
+    N = 20
+    channel_out_power = np.array([np.abs(channel[start_ix+t:start_ix+t+T_NULL]).sum() for t in range(0, T_TF-T_NULL, N)])
 
-t_null = 10 * channel_out_power.argmin()
+    # 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 = N * channel_out_power.argmin()
 
-print("  NULL symbol starts at ix={}".format(t_null))
+    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.
+    # 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)))
 
-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")
 
-# 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")
+    # We start a bit earlier than the end of the null symbol
+    corr_start_ix = t_null + T_NULL - 50
 
-# We start a bit earlier than the end of the null symbol
-corr_start_ix = t_null + T_NULL - 50
+    # 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.
+    max_component_delay = 500 # T
 
-# 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.
+    cir = np.array([np.abs(np.corrcoef(channel[start_ix + corr_start_ix + i:start_ix + corr_start_ix + phase_ref.size + i], phase_ref)[0,1]) for i in range(max_component_delay)])
+    return cir
 
-max_component_delay = 1000 # T
+num_correlations = int(len(channel_out) / T_TF)
+print("Doing {} correlations".format(num_correlations))
 
-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)])
+cirs = np.array([
+    calc_cir(channel_out, i * T_TF)
+    for i in range(num_correlations) ])
 
-print("Done")
+print("Plotting")
 
-fig = pp.figure()
-ax = fig.add_subplot(111)
-hi = ax.plot(channel_out_power)
+pp.subplot(211)
+pp.plot(cirs.sum(axis=0))
+pp.subplot(212)
+pp.imshow(cirs)
 
-fig = pp.figure()
-ax = fig.add_subplot(111)
-hi = ax.plot(cir)
+print("Done")
 
 pp.show()
 
 
+