Add some TII testing scripts

6 files changed, 778 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
index fce88b8..ec33b70 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,4 @@ rdsparse/Makefile
 rdsparse/cmake_install.cmake
 rdsparse/cmake_uninstall.cmake
 rdsparse/rdsparse
+ofdm*iq
diff --git a/tii/README.md b/tii/README.md
new file mode 100644
index 0000000..5887ee4
--- /dev/null
+++ b/tii/README.md
@@ -0,0 +1 @@
+A few helper scripts to develop and debug TII
diff --git a/tii/analyse_tii_from_file.py b/tii/analyse_tii_from_file.py
new file mode 100755
index 0000000..9d834ea
--- /dev/null
+++ b/tii/analyse_tii_from_file.py
@@ -0,0 +1,453 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Analyse TII of iq file. complexf 2048000
+#
+# It requires SciPy and matplotlib.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import functools
+import numpy as np
+import matplotlib.pyplot as pp
+import argparse
+import subprocess
+from scipy.misc import imsave
+from tii import calculate_carrier_indices, calculate_reduced_carrier_indices
+
+all_tii_patterns = {}
+cp_per_carrier = {}
+
+def prepare_tii_patterns():
+    print("Building TII table...")
+    global all_tii_patterns
+    global cp_per_carrier
+
+    for c in range(24):
+        for p in range(70):
+            carriers = calculate_carrier_indices(c, p)
+            all_tii_patterns[tuple(carriers)] = (c, p)
+
+            carriers = calculate_reduced_carrier_indices(c, p)
+            for k in carriers:
+                if k not in cp_per_carrier:
+                    cp_per_carrier[k] = set()
+                cp_per_carrier[k].add((c, p))
+
+SIZEOF_SAMPLE = 8 # complex floats
+
+# Constants for TM 1
+NbSymbols = 76
+NbCarriers = 1536
+Spacing = 2048
+NullSize = 2656
+SymSize = 2552
+FicSizeOut = 288
+
+CyclicPrefix = SymSize - Spacing
+
+prs_phases = np.array([0, -1, 1, -1, -1, -1, -1, 0, 0, 1, -1, -1, -1, 1, 1, 0, 0, -1, 1,
+        -1, -1, -1, -1, 0, 0, 1, -1, -1, -1, 1, 1, 0, 0, 1, 2, -1, 2, 1, 0, 0,
+        -1, -1, 0, -1, 2, -1, 2, 0, -1, 1, 2, -1, 2, 1, 0, 0, -1, -1, 0, -1, 2,
+        -1, 2, 0, -1, 1, 1, 1, -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 1,
+        1, 1, -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 1, 0, -1, 0, 1, 2,
+        0, 1, -1, 2, -1, 0, -1, 0, 0, 1, 1, 0, -1, 0, 1, 2, 0, 1, -1, 2, -1, 0,
+        -1, 0, 0, 1, 2, 0, 2, 2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, 0,
+        2, 2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, -1, 0, -1, 2, 1, 1, 0,
+        0, 1, 0, -1, 0, -1, 1, 0, 2, -1, 0, -1, 2, 1, 1, 0, 0, 1, 0, -1, 0, -1,
+        1, 0, 1, 1, 1, -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 1, 1, 1,
+        -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 0, -1, 2, -1, 0, 1, -1, 0,
+        2, 1, 2, -1, 2, -1, -1, 0, 0, -1, 2, -1, 0, 1, -1, 0, 2, 1, 2, -1, 2,
+        -1, -1, 0, 2, 0, 2, 2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, 0, 2,
+        2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, -1, 0, -1, 2, 1, 1, 0, 0,
+        1, 0, -1, 0, -1, 1, 0, 2, -1, 0, -1, 2, 1, 1, 0, 0, 1, 0, -1, 0, -1, 1,
+        0, -1, -1, -1, 1, -1, 1, 0, 2, 1, 1, -1, 1, 1, -1, 0, 2, -1, -1, -1, 1,
+        -1, 1, 0, 2, 1, 1, -1, 1, 1, -1, 0, 2, -1, 2, 1, 2, -1, 0, 2, -1, 1, 0,
+        1, 2, 1, 2, 2, -1, -1, 2, 1, 2, -1, 0, 2, -1, 1, 0, 1, 2, 1, 2, 2, -1,
+        0, 2, 0, 0, 0, 0, 1, 1, 2, 0, 0, 0, 2, 2, 1, 1, 0, 2, 0, 0, 0, 0, 1, 1,
+        2, 0, 0, 0, 2, 2, 1, 1, 2, -1, 0, -1, 2, 1, 1, 0, 0, 1, 0, -1, 0, -1,
+        1, 0, 2, -1, 0, -1, 2, 1, 1, 0, 0, 1, 0, -1, 0, -1, 1, 0, 1, 1, 1, -1,
+        1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 1, 1, 1, -1, 1, -1, 2, 0, -1,
+        -1, 1, -1, -1, 1, 2, 0, -1, 2, 1, 2, -1, 0, 2, -1, 1, 0, 1, 2, 1, 2, 2,
+        -1, -1, 2, 1, 2, -1, 0, 2, -1, 1, 0, 1, 2, 1, 2, 2, -1, -1, 1, -1, -1,
+        -1, -1, 0, 0, 1, -1, -1, -1, 1, 1, 0, 0, -1, 1, -1, -1, -1, -1, 0, 0,
+        1, -1, -1, -1, 1, 1, 0, 0, -1, 0, 1, 0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0,
+        2, 1, -1, 0, 1, 0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0, 2, 1, -1, -1, -1, 1,
+        -1, 1, 0, 2, 1, 1, -1, 1, 1, -1, 0, 2, -1, -1, -1, 1, -1, 1, 0, 2, 1,
+        1, -1, 1, 1, -1, 0, 2, 0, -1, 2, -1, 0, 1, -1, 0, 2, 1, 2, -1, 2, -1,
+        -1, 0, 0, -1, 2, -1, 0, 1, -1, 0, 2, 1, 2, -1, 2, -1, -1, 0, -1, 1, -1,
+        -1, -1, -1, 0, 0, 1, -1, -1, -1, 1, 1, 0, 0, -1, 1, -1, -1, -1, -1, 0,
+        0, 1, -1, -1, -1, 1, 1, 0, 0, 0, 1, 2, 1, 0, -1, -1, 2, 2, -1, 2, 1, 2,
+        1, -1, 2, 0, 1, 2, 1, 0, -1, -1, 2, 2, -1, 2, 1, 2, 1, -1, 2, 1, 1, 1,
+        -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 1, 1, 1, -1, 1, -1, 2, 0,
+        -1, -1, 1, -1, -1, 1, 2, 0, 1, 0, -1, 0, 1, 2, 0, 1, -1, 2, -1, 0, -1,
+        0, 0, 1, 1, 0, -1, 0, 1, 2, 0, 1, -1, 2, -1, 0, -1, 0, 0, 1, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 1, 1, 1, 1, 2, 2,
+-1, 1, 1, 1, -1, -1, 2, 2, 1, -1, 1, 1, 1, 1, 2, 2, -1, 1, 1, 1, -1, -1, 2, 2,
+2, 1, 0, 1, 2, -1, 1, 2, 0, -1, 0, 1, 0, 1, 1, 2, 2, 1, 0, 1, 2, -1, 1, 2, 0,
+-1, 0, 1, 0, 1, 1, 2, 0, 0, 0, 2, 0, 2, 1, -1, 2, 2, 0, 2, 2, 0, 1, -1, 0, 0,
+0, 2, 0, 2, 1, -1, 2, 2, 0, 2, 2, 0, 1, -1, 1, 2, -1, 2, 1, 0, 0, -1, -1, 0,
+-1, 2, -1, 2, 0, -1, 1, 2, -1, 2, 1, 0, 0, -1, -1, 0, -1, 2, -1, 2, 0, -1, -1,
+1, -1, -1, -1, -1, 0, 0, 1, -1, -1, -1, 1, 1, 0, 0, -1, 1, -1, -1, -1, -1, 0,
+0, 1, -1, -1, -1, 1, 1, 0, 0, 2, 1, 0, 1, 2, -1, 1, 2, 0, -1, 0, 1, 0, 1, 1, 2,
+2, 1, 0, 1, 2, -1, 1, 2, 0, -1, 0, 1, 0, 1, 1, 2, 2, 2, 2, 0, 2, 0, -1, 1, 0,
+0, 2, 0, 0, 2, -1, 1, 2, 2, 2, 0, 2, 0, -1, 1, 0, 0, 2, 0, 0, 2, -1, 1, -1, 0,
+1, 0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0, 2, 1, -1, 0, 1, 0, -1, 2, 2, 1, 1, 2, 1,
+0, 1, 0, 2, 1, 2, 0, 2, 2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, 0, 2, 2,
+2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 1, 0, -1, 0, 1, 2, 0, 1, -1, 2, -1, 0,
+-1, 0, 0, 1, 1, 0, -1, 0, 1, 2, 0, 1, -1, 2, -1, 0, -1, 0, 0, 1, 2, 2, 2, 0, 2,
+0, -1, 1, 0, 0, 2, 0, 0, 2, -1, 1, 2, 2, 2, 0, 2, 0, -1, 1, 0, 0, 2, 0, 0, 2,
+-1, 1, -1, 0, 1, 0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0, 2, 1, -1, 0, 1, 0, -1, 2, 2,
+1, 1, 2, 1, 0, 1, 0, 2, 1, 1, -1, 1, 1, 1, 1, 2, 2, -1, 1, 1, 1, -1, -1, 2, 2,
+1, -1, 1, 1, 1, 1, 2, 2, -1, 1, 1, 1, -1, -1, 2, 2, 2, 1, 0, 1, 2, -1, 1, 2, 0,
+-1, 0, 1, 0, 1, 1, 2, 2, 1, 0, 1, 2, -1, 1, 2, 0, -1, 0, 1, 0, 1, 1, 2, -1, -1,
+-1, 1, -1, 1, 0, 2, 1, 1, -1, 1, 1, -1, 0, 2, -1, -1, -1, 1, -1, 1, 0, 2, 1, 1,
+-1, 1, 1, -1, 0, 2, -1, 0, 1, 0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0, 2, 1, -1, 0, 1,
+0, -1, 2, 2, 1, 1, 2, 1, 0, 1, 0, 2, 1, 2, 0, 2, 2, 2, 2, -1, -1, 0, 2, 2, 2,
+0, 0, -1, -1, 2, 0, 2, 2, 2, 2, -1, -1, 0, 2, 2, 2, 0, 0, -1, -1, 2, 1, 0, 1,
+2, -1, 1, 2, 0, -1, 0, 1, 0, 1, 1, 2, 2, 1, 0, 1, 2, -1, 1, 2, 0, -1, 0, 1, 0,
+1, 1, 2, 2, 2, 2, 0, 2, 0, -1, 1, 0, 0, 2, 0, 0, 2, -1, 1, 2, 2, 2, 0, 2, 0,
+-1, 1, 0, 0, 2, 0, 0, 2, -1, 1, 1, 2, -1, 2, 1, 0, 0, -1, -1, 0, -1, 2, -1, 2,
+0, -1, 1, 2, -1, 2, 1, 0, 0, -1, -1, 0, -1, 2, -1, 2, 0, -1, 1, -1, 1, 1, 1, 1,
+2, 2, -1, 1, 1, 1, -1, -1, 2, 2, 1, -1, 1, 1, 1, 1, 2, 2, -1, 1, 1, 1, -1, -1,
+2, 2, -1, 2, 1, 2, -1, 0, 2, -1, 1, 0, 1, 2, 1, 2, 2, -1, -1, 2, 1, 2, -1, 0,
+2, -1, 1, 0, 1, 2, 1, 2, 2, -1, 1, 1, 1, -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1,
+2, 0, 1, 1, 1, -1, 1, -1, 2, 0, -1, -1, 1, -1, -1, 1, 2, 0, 2, -1, 0, -1, 2, 1,
+1, 0, 0, 1, 0, -1, 0, -1, 1, 0, 2, -1, 0, -1, 2, 1, 1, 0, 0, 1, 0, -1, 0, -1,
+1, 0])
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot TII")
+    parser.add_argument('--frame', default='0', help='Which transmission frame to analyse',
+            required=False)
+    parser.add_argument('--samplerate', default='2048000', help='Sample rate',
+            required=False)
+    parser.add_argument('--iq-file', default='ofdm.iq', help='File to read',
+            required=False)
+    parser.add_argument('--test-combs', action='store_const', const=True, help='Test all comb pattern values',
+            required=False)
+    parser.add_argument('--old-algo', action='store_const', const=True, help='Do an analysis with the old algorithm',
+            required=False)
+
+    cli_args = parser.parse_args()
+
+    if cli_args.test_combs:
+        prepare_tii_patterns()
+        prepare_comb_pattern_run(cli_args, 1, 11)
+        prepare_comb_pattern_run(cli_args, 1, 12)
+    elif cli_args.old_algo:
+        prepare_tii_patterns()
+        plot_tii_once(cli_args)
+    else:
+        prepare_tii_patterns()
+        algo1(cli_args)
+
+tii_ini_template = open("tii.ini.template", "r").read()
+
+def prepare_comb_pattern_run(options, comb, pattern):
+    fd = open("tii.ini", "w")
+    fd.write(tii_ini_template.format(comb=comb, pattern=pattern, oldvariant=0))
+    fd.close()
+
+    subprocess.call(['./odr-dabmod', 'tii.ini'])
+
+    ret = plot_tii_once(options)
+    if ret[0] == "ok":
+        _, c, p, variant = ret
+        if c == comb and p == pattern:
+            print("== OK    {}".format(ret))
+        else:
+            print("== WRONG {}".format(ret))
+    else:
+        print("== FAIL  {}".format(ret))
+
+
+def plot_tii_once(options):
+    oversample = int(int(options.samplerate) / 2048000)
+    num_samples_per_tf = oversample * (NullSize + NbSymbols * SymSize)
+    print("{} samples per TF".format(num_samples_per_tf))
+
+    N = 8
+    num_samples = num_samples_per_tf * (int(options.frame) + N)
+
+    frames = np.fromfile(options.iq_file, np.complex64, count=num_samples)
+
+    #skip to frame
+    frames = frames[num_samples_per_tf * int(options.frame):]
+    print("{} samples loaded".format(len(frames)))
+
+    #fig = pp.figure()
+    #fig.suptitle("TII analysis")
+    #ax1 = fig.add_subplot(N, 1, i+1)
+    #ax1.set_title("Frame {}".format(i))
+    #ax1.plot(np.abs(null_fft))
+
+    # analyse N successive frames
+    for i in range(N):
+        frame = frames[num_samples_per_tf * i:num_samples_per_tf * (i + 1)]
+
+        if np.max(np.abs(frame[:NullSize])) != 0:
+            print("Frame {} has nonzero power in null symbol".format(i))
+
+        # Take the NULL symbol from that frame, but skip the cyclic prefix and truncate
+        null_skip = NullSize - Spacing
+        frame_null = frame[oversample*null_skip:oversample*(null_skip + Spacing)]
+
+        null_fft = np.fft.fft(frame_null, oversample*Spacing)
+        null_fft_abs = np.abs(null_fft)
+
+        # The phase reference symbol, skip the cyclic prefix
+        frame_prs = frame[oversample*(NullSize + CyclicPrefix):oversample * (NullSize + CyclicPrefix + Spacing)]
+        prs_fft = np.fft.fft(frame_prs, oversample*Spacing)
+
+        threshold = 0.01
+        carrier_indices = np.nonzero(null_fft_abs > threshold)[0]
+        print("Frame {} has {} carriers".format(i, len(carrier_indices)))
+
+        if len(carrier_indices) != 32 and len(carrier_indices) != 0:
+            print("Invalid number of TII carriers")
+            return ('fail', )
+
+        def ix_to_k(ix):
+            "Convert a FFT index to a carrier index k"
+            # ix goes from 0 to 2047
+            # 0..1024 positive, 1024..2047 negative, 0 is DC, 1..768 and 1280..2047 have power.
+            carrieroffset = int(oversample*Spacing)
+            if ix <= 1024:
+                return ix
+            else:
+                return ix - carrieroffset
+
+        analyse_prs = False
+        if len(carrier_indices) and not analyse_prs:
+            print("k              =" + " ".join("{:>4}".format(ix_to_k(ix)) for ix in carrier_indices))
+            print("phases         =" + " ".join("{:>4}".format(int(round(p))) for p in np.angle(prs_fft[null_fft_abs > threshold], deg=True)))
+
+            phases_tii = np.angle(null_fft[null_fft_abs > threshold], deg=True)
+            phases_prs = np.angle(prs_fft[null_fft_abs > threshold], deg=True)
+            phase_error_wrong_impl = np.mod(np.asarray(np.around(phases_tii - phases_prs), dtype=np.int), 360)
+
+            print("phases wrong   =" + " ".join("{:>4}".format(p) for p in phase_error_wrong_impl))
+
+            pair_indices = carrier_indices.reshape(-1,2)[...,0]
+            correct_phase_indices = np.stack([pair_indices, pair_indices]).T.flatten()
+            phases_prs2 = np.angle(prs_fft[correct_phase_indices], deg=True)
+            phase_error_correct_impl = np.mod(np.asarray(np.around(phases_tii - phases_prs2), dtype=np.int), 360)
+            print("phases correct =" + " ".join("{:>4}".format(p) for p in phase_error_correct_impl))
+
+            variant = "unknown"
+            if np.all(phase_error_correct_impl == 0):
+                print("TII uses correct implementation variant")
+                variant = "correct"
+            elif np.all(phase_error_wrong_impl == 0):
+                print("TII uses old wrong implementation variant")
+                variant = "wrong"
+            else:
+                print("TII is wrong")
+
+            carrier_indices_t = tuple(sorted(ix_to_k(c) for c in carrier_indices))
+            print(carrier_indices_t)
+            if carrier_indices_t in all_tii_patterns:
+                c, p = all_tii_patterns[carrier_indices_t]
+                print("TII is using comb {} and pattern {}".format(c, p))
+                return ("ok", c, p, variant)
+            else:
+                print("Unrecognised TII comb and pattern")
+
+            return ("fail", variant)
+
+        if analyse_prs:
+            delta = ((prs_phases[null_fft_abs > threshold] * 90) - np.angle(prs_fft[null_fft_abs > threshold], deg=True))
+            print("prs=" + " ".join("{:>4}".format(int(round(p)) % 360) for p in delta))
+
+    #pp.show()
+
+def algo1(options):
+    oversample = 1
+    if int(options.samplerate) != 2048000:
+        print("oversampling not supported for new algo")
+        sys.exit(1)
+
+    num_samples_per_tf = NullSize + NbSymbols * SymSize
+    print("{} samples per TF".format(num_samples_per_tf))
+
+    def ix_to_k(ix):
+        "Convert a FFT index to a carrier index k"
+        # ix goes from 0 to 2047
+        # 0..1024 positive, 1024..2047 negative, 0 is DC, 1..768 and 1280..2047 have power.
+        carrieroffset = int(Spacing)
+        if ix <= 1024:
+            return ix
+        else:
+            return ix - carrieroffset
+
+    N = 8
+    num_samples = num_samples_per_tf * (int(options.frame) + N)
+
+    frames = np.fromfile(options.iq_file, np.complex64, count=num_samples)
+
+    #skip to frame
+    frames = frames[num_samples_per_tf * int(options.frame):]
+    print("{} samples loaded".format(len(frames)))
+
+    # analyse N successive frames
+    for i in range(N):
+        frame = frames[num_samples_per_tf * i:num_samples_per_tf * (i + 1)]
+
+        # Take the NULL symbol from that frame, but skip the cyclic prefix and truncate
+        null_skip = NullSize - Spacing
+        frame_null = frame[null_skip:(null_skip + Spacing)]
+
+        null_fft = np.fft.fft(frame_null, Spacing)
+
+        # The phase reference symbol, skip the cyclic prefix
+        frame_prs = frame[NullSize + CyclicPrefix:NullSize + CyclicPrefix + Spacing]
+        prs_fft = np.fft.fft(frame_prs, Spacing)
+
+        # In TM1, the carriers repeat four times:
+        # [-768, -384[
+        # [-384, 0[
+        # ]0, 384]
+        # ]384, 768]
+        # A consequence of the fact that the 0 bin is never used is that the
+        # first carrier of each pair is even for negative k, odd for positive k
+
+        blocks = [null_fft[-769:-385], null_fft[-385:-1], null_fft[1:385], null_fft[385:769]]
+        blocks_multiplied = np.zeros(384//2, dtype=np.complex128)
+        for block in blocks:
+            even_odd = block.reshape(-1, 2)
+            b = even_odd[...,0] * np.conj(even_odd[...,1])
+            blocks_multiplied += b
+
+        threshold = 0.01
+        indices_above_threshold = np.nonzero(np.abs(blocks_multiplied) > threshold)[0]
+        carrier_indices = indices_above_threshold * 2 + 1
+
+        if np.max(np.abs(frame[:NullSize])) != 0:
+            print("Frame {} has nonzero power in null symbol and {} TII carriers".format(i, len(carrier_indices)))
+
+        if len(carrier_indices) >= 4:
+            carriers = [ix_to_k(ix) for ix in carrier_indices]
+
+            cp_counts = {}
+            for k in carriers:
+                if k in cp_per_carrier:
+                    cps = cp_per_carrier[k]
+                    for cp in cps:
+                        if cp not in cp_counts:
+                            cp_counts[cp] = 0
+                        cp_counts[cp] += 1
+
+            if len(cp_counts) > 0:
+                for cp in cp_counts:
+                    if cp_counts[cp] >= 4:
+                        c, p = cp
+                        valid, delay, err = analyse_phase(c, p, null_fft, prs_fft)
+                        if valid:
+                            print("TII likelihood {}: comb {} and pattern {}, delay {} samples, err {}".format(
+                                cp_counts[cp], c, p,
+                                delay, err))
+                        else:
+                            print("TII likelihood {}: comb {} and pattern {}, invalid delay measurement, err {}".format(
+                                cp_counts[cp], c, p, err))
+            else:
+                print("Unrecognised TII comb and pattern")
+
+def convert_angles(angle):
+    """Convert 0 to 360 degree angles to 0 to 180, -179 to -1"""
+    if angle > 180:
+        return angle - 360
+    elif angle < -180:
+        return angle + 360
+    else:
+        return angle
+phase_corrector = np.vectorize(convert_angles)
+
+def analyse_phase(c, p, null_fft, prs_fft):
+    carriers = np.array(calculate_carrier_indices(c, p))
+
+    #print("k              =" + " ".join("{:>4}".format(k) for k in carriers))
+    #print("phases         =" + " ".join("{:>4}".format(int(round(p))) for p in np.angle(prs_fft[carriers], deg=True)))
+
+    phases_tii = np.angle(null_fft[carriers], deg=True)
+
+    pair_indices = carriers.reshape(-1,2)[...,0]
+    correct_phase_indices = np.stack([pair_indices, pair_indices]).T.flatten()
+    phases_prs2 = np.angle(prs_fft[correct_phase_indices], deg=True)
+
+    fig = pp.figure()
+    fig.suptitle("Phase error C {} P {}".format(c, p))
+
+    phase_errors = {}
+    search_range = (-40, 40)
+    for err4 in range(*search_range):
+        err = err4 / 4
+        rotate_vec = np.exp(2j * np.pi * err * carriers / 2048)
+        rotated = np.angle(null_fft[carriers] * rotate_vec, deg=True)
+        delta = np.mod(np.asarray(np.around(rotated - phases_prs2), dtype=np.int), 360)
+        delta = phase_corrector(delta)
+        #print("phases {:>3} =".format(err) + " ".join("{:>4}".format(p) for p in delta))
+        phase_errors[err] = delta
+
+        #ax1 = fig.add_subplot(7, 1, i+1)
+        #ax1.set_title("err {}".format(err))
+        #ax1.plot(correct_phase_indices, delta, '.')
+        #ax1.plot(correct_phase_indices, slopepoly(correct_phase_indices))
+    #pp.show()
+
+    best_err = min(phase_errors.items(), key=(lambda e: np.sum(np.abs(e[1]))))
+    #print("Best err: {}".format(best_err[0]))
+    #print("phases =" + " ".join("{:>4}".format(p) for p in best_err[1]))
+
+    meas_err = np.sum(np.abs(best_err[1]))
+
+    if best_err[0] != search_range[0] and best_err[0] != search_range[1]:
+        return (True, best_err[0], meas_err)
+    else:
+        return (False, best_err[0], meas_err)
+
+main()
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2018 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
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
diff --git a/tii/merge_iq.py b/tii/merge_iq.py
new file mode 100755
index 0000000..e8e3119
--- /dev/null
+++ b/tii/merge_iq.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Read two IQ files, complexf 2048000, and sum them
+# together with separate amplitudes and add noise.
+#
+# It requires SciPy and matplotlib.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import numpy as np
+
+import argparse
+
+def main():
+    parser = argparse.ArgumentParser(description="MergeIQ")
+    #parser.add_argument('--noise', default='0', help='Noise amplitude to add', required=False)
+    parser.add_argument('--delay', default='20', help='Number of samples delay of second file', required=False)
+    #parser.add_argument('--noise', action='store_const', const=True, help='Test all comb pattern values', required=False)
+
+    cli_args = parser.parse_args()
+
+    num_samples = 1024
+    delay = int(cli_args.delay)
+
+    fd1 = open("ofdm-c1p12.iq", "rb")
+    fd2 = open("ofdm-c2p12.iq", "rb")
+    fd_out = open("ofdm.iq", "wb")
+
+    frame1 = np.fromfile(fd1, np.complex64, count=num_samples)
+    frame2 = np.fromfile(fd2, np.complex64, count=num_samples-delay)
+    frame2 = np.concatenate((np.zeros(delay, dtype=np.complex64), frame2))
+
+    totalsize = 0
+    while True:
+        n = 0.000001
+        noise_r = np.random.normal(0, n, num_samples)
+        noise_i = np.random.normal(0, n, num_samples)
+        noise = noise_r + 1j * noise_i
+        sum_frame = 0.8 * frame1 + 0.22 * frame2 + noise.astype(np.complex64)
+        totalsize += len(sum_frame)
+        fd_out.write(sum_frame.tobytes())
+
+        frame1 = np.fromfile(fd1, np.complex64, count=num_samples)
+        frame2 = np.fromfile(fd2, np.complex64, count=num_samples)
+
+        if len(frame1) != num_samples or len(frame2) != num_samples:
+            print("Stop {}".format(totalsize))
+            break
+
+main()
diff --git a/tii/tii.ini.template b/tii/tii.ini.template
new file mode 100644
index 0000000..51b9bee
--- /dev/null
+++ b/tii/tii.ini.template
@@ -0,0 +1,51 @@
+[remotecontrol]
+telnet=0
+telnetport=2121
+zmqctrl=0
+zmqctrlendpoint=tcp://127.0.0.1:9400
+
+[log]
+syslog=0
+filelog=0
+filename=odr-dabmod.log
+
+[input]
+transport=file
+source=/home/bram/dab/mmbtools-aux/eti/funk.eti
+loop=0
+
+[modulator]
+digital_gain=0.8
+rate=2048000
+mode=1
+
+[cfr]
+enable=0
+clip=70.0
+error_clip=0.05
+
+[firfilter]
+enabled=0
+
+[poly]
+enabled=0
+polycoeffile=polyCoefs
+
+[output]
+output=file
+
+[fileoutput]
+format=complexf_normalised
+filename=ofdm.iq
+#format=u8
+#filename=20180207-ofdm.u8.iq
+
+[zmqoutput]
+listen=tcp://*:54001
+socket_type=pub
+
+[tii]
+enable=1
+comb={comb}
+pattern={pattern}
+old_variant={oldvariant}
diff --git a/tii/tii.py b/tii/tii.py
new file mode 100755
index 0000000..0ec9429
--- /dev/null
+++ b/tii/tii.py
@@ -0,0 +1,219 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Does some TII pattern calculations
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import argparse
+import time
+
+TII_PATTERN = [ # {{{
+    [0,0,0,0,1,1,1,1],
+    [0,0,0,1,0,1,1,1],
+    [0,0,0,1,1,0,1,1],
+    [0,0,0,1,1,1,0,1],
+    [0,0,0,1,1,1,1,0],
+    [0,0,1,0,0,1,1,1],
+    [0,0,1,0,1,0,1,1],
+    [0,0,1,0,1,1,0,1],
+    [0,0,1,0,1,1,1,0],
+    [0,0,1,1,0,0,1,1],
+    [0,0,1,1,0,1,0,1],
+    [0,0,1,1,0,1,1,0],
+    [0,0,1,1,1,0,0,1],
+    [0,0,1,1,1,0,1,0],
+    [0,0,1,1,1,1,0,0],
+    [0,1,0,0,0,1,1,1],
+    [0,1,0,0,1,0,1,1],
+    [0,1,0,0,1,1,0,1],
+    [0,1,0,0,1,1,1,0],
+    [0,1,0,1,0,0,1,1],
+    [0,1,0,1,0,1,0,1],
+    [0,1,0,1,0,1,1,0],
+    [0,1,0,1,1,0,0,1],
+    [0,1,0,1,1,0,1,0],
+    [0,1,0,1,1,1,0,0],
+    [0,1,1,0,0,0,1,1],
+    [0,1,1,0,0,1,0,1],
+    [0,1,1,0,0,1,1,0],
+    [0,1,1,0,1,0,0,1],
+    [0,1,1,0,1,0,1,0],
+    [0,1,1,0,1,1,0,0],
+    [0,1,1,1,0,0,0,1],
+    [0,1,1,1,0,0,1,0],
+    [0,1,1,1,0,1,0,0],
+    [0,1,1,1,1,0,0,0],
+    [1,0,0,0,0,1,1,1],
+    [1,0,0,0,1,0,1,1],
+    [1,0,0,0,1,1,0,1],
+    [1,0,0,0,1,1,1,0],
+    [1,0,0,1,0,0,1,1],
+    [1,0,0,1,0,1,0,1],
+    [1,0,0,1,0,1,1,0],
+    [1,0,0,1,1,0,0,1],
+    [1,0,0,1,1,0,1,0],
+    [1,0,0,1,1,1,0,0],
+    [1,0,1,0,0,0,1,1],
+    [1,0,1,0,0,1,0,1],
+    [1,0,1,0,0,1,1,0],
+    [1,0,1,0,1,0,0,1],
+    [1,0,1,0,1,0,1,0],
+    [1,0,1,0,1,1,0,0],
+    [1,0,1,1,0,0,0,1],
+    [1,0,1,1,0,0,1,0],
+    [1,0,1,1,0,1,0,0],
+    [1,0,1,1,1,0,0,0],
+    [1,1,0,0,0,0,1,1],
+    [1,1,0,0,0,1,0,1],
+    [1,1,0,0,0,1,1,0],
+    [1,1,0,0,1,0,0,1],
+    [1,1,0,0,1,0,1,0],
+    [1,1,0,0,1,1,0,0],
+    [1,1,0,1,0,0,0,1],
+    [1,1,0,1,0,0,1,0],
+    [1,1,0,1,0,1,0,0],
+    [1,1,0,1,1,0,0,0],
+    [1,1,1,0,0,0,0,1],
+    [1,1,1,0,0,0,1,0],
+    [1,1,1,0,0,1,0,0],
+    [1,1,1,0,1,0,0,0],
+    [1,1,1,1,0,0,0,0] ] # }}}
+
+def calculate_A_c_p(c, p, k):
+    """ETSI EN 300 401 14.8.1"""
+    r = 0
+
+    if k == 0 or k == -769:
+        return 0
+    elif -768 <= k and k < -384:
+        for b in range(8):
+            if k == -768 + 2*c + 48*b and TII_PATTERN[p][b]:
+                r += 1
+    elif -384 <= k and k < 0:
+        for b in range(8):
+            if k == -384 + 2*c + 48*b and TII_PATTERN[p][b]:
+                r += 1
+    elif 0 < k and k <= 384:
+        for b in range(8):
+            if k == 1 + 2*c + 48*b and TII_PATTERN[p][b]:
+                r += 1
+    elif 384 < k and k <= 768:
+        for b in range(8):
+            if k == 385 + 2*c + 48*b and TII_PATTERN[p][b]:
+                r += 1
+    else:
+        raise ValueError("Invalid k={}".format(k))
+
+    # I'm expecting that r means "enable or disable", nothing else
+    assert(r < 2);
+    return r
+
+def calculate_reduced_carrier_indices(c, p):
+    carriers = []
+    for k in range(384):
+        for b in range(8):
+            if k == 1 + 2*c + 48*b and TII_PATTERN[p][b]:
+                carriers.append(k)
+    return carriers
+
+def calculate_carrier_indices(c, p):
+    carriers = calculate_reduced_carrier_indices(c, p)
+
+    # because of z_m,0,k, we enable carrier k if A_c,p(k) or A_c,p(k-1) are
+    # true
+    carriers += [k + 1 for k in carriers]
+    carriers += [k - 769 for k in carriers] + \
+            [k - 385 for k in carriers] + \
+            [k + 384 for k in carriers]
+    carriers.sort()
+    return carriers
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Calculate TII")
+    parser.add_argument('--all', action='store_const', const=True,
+            help='Calulate all carrier indices', required=False)
+    parser.add_argument('--check-zero', action='store_const', const=True,
+            help='Check if any pattern sets zero freq', required=False)
+    parser.add_argument('--test1', action='store_const', const=True,
+            help='Do equivalence test 1', required=False)
+    cli_args = parser.parse_args()
+
+
+    if cli_args.all:
+        for c in range(24):
+            for p in range(70):
+                carriers = calculate_carrier_indices(c, p)
+                print(c, p, carriers, len(carriers))
+
+    elif cli_args.check_zero:
+        for c in range(24):
+            for p in range(70):
+                carriers = calculate_carrier_indices(c, p)
+                if 0 in carriers:
+                    print("Found 0 in")
+                elif 1 in carriers:
+                    print("Found 1 in")
+                elif -1 in carriers:
+                    print("Found -1 in")
+                else:
+                    continue
+                print("Car: ", c, p, carriers, len(carriers))
+
+    elif cli_args.test1:
+        print("Running compare test table")
+        duration1 = 0
+        duration2 = 0
+        for c in range(24):
+            for p in range(70):
+                start1 = time.time()
+                carriers0 = [k
+                        for k in range(-768, 769)
+                        if calculate_A_c_p(c, p, k) or calculate_A_c_p(c, p, k-1)]
+                carriers0.sort()
+                duration1 += (time.time() - start1)
+
+                start2 = time.time()
+                carriers1 = calculate_carrier_indices(c, p)
+                duration2 += (time.time() - start2)
+
+                if carriers0 != carriers1:
+                    print("0: {}".format(carriers0))
+                    print("1: {}".format(carriers1))
+                    sys.exit(1)
+        print("Comparison successful: {} vs {}".format(duration1, duration2))
+
+    else:
+        print("Example given in standard:")
+        c = 1
+        p = 11
+        acps = [k
+                for k in range(-768, 769)
+                if calculate_A_c_p(c, p, k)]
+        print("ACP: ", c, p, acps, len(acps))
+        carriers = calculate_carrier_indices(c, p)
+        print("Car: ", c, p, carriers, len(carriers))
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2018 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
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.