path: root/tii/analyse_tii_from_file.py

#!/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 load_iq(filename, u8, count):
    if u8:
        u8_interleaved = np.fromfile(filename, np.uint8, count=count)
        u8_iq = u8_interleaved.reshape(int(len(u8_interleaved)/2), 2)
        # This directly converts to fc64
        fc64_unscaled = u8_iq[...,0] + np.complex64(1j) * u8_iq[...,1]
        fc64_scaled = (fc64_unscaled - 127.0) / 128.0
        #fc64_dc_comp = fc64_scaled - np.average(fc64_scaled)
        return fc64_scaled
    else:
        return np.fromfile(filename, np.complex64, count=count)

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)
    parser.add_argument('--u8', action='store_const', const=True,
            help='File is in u8 and not cf32 format', required=False)
    parser.add_argument('--align', action='store_const', const=True,
            help='Find NULL and align to it', 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)


def prepare_comb_pattern_run(options, comb, pattern):
    tii_ini_template = open("tii.ini.template", "r").read()
    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 advance_to_null(frames):
    for i in range(0, len(frames), NullSize//2):
        if np.all(np.abs(frames[i:i+NullSize//2]) < 1e-2):
            print("NULL at?", i)
            # First non-zero sample is first PRS sample
            for j in range(NullSize):
                if np.abs(frames[i+j]) > 1e-2:
                    print("NULL at ", i+j - NullSize)
                    return frames[i+j - NullSize:]
    print("NULL not found")
    sys.exit(1)

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 = load_iq(options.iq_file, options.u8, 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 = load_iq(options.iq_file, options.u8, count=num_samples)

    if options.align:
        frames = advance_to_null(frames)

    #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[-768:-384], null_fft[-384:], null_fft[1:385], null_fft[385:769]]
        blocks_multiplied = np.zeros(384//2, dtype=np.complex128)
        #print("blocks {}".format(i))
        for block in blocks:
            even_odd = block.reshape(-1, 2)
            b = even_odd[...,0] * np.conj(even_odd[...,1])
            #print("".join('#' if i else '.' for i in np.abs(b) > 1e-7))
            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 len(frame) < NullSize:
            break

        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:
                print(carriers)
                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 2pi degree angles to -pi to pi"""
    if angle * 2.0 > np.pi:
        return angle - 2*np.pi
    elif angle * 2.0 < -np.pi:
        return angle + 2*np.pi
    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])

    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])

    fig = pp.figure()
    fig.suptitle("Phase error C {} P {}".format(c, p))

    phase_errors = {}
    search_range = (-20, 500)
    for err in range(*search_range):
        rotate_vec = np.exp(2j * np.pi * err * carriers / 2048)
        rotated = null_fft[carriers] * rotate_vec
        delta = np.asarray(np.around(np.angle(rotated) - phases_prs2), dtype=np.int)
        delta = phase_corrector(delta)
        #print("phases {:>3} =".format(err) + " ".join("{:>4}".format(int(180*p/np.pi)) 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.