From a14cb8afb8fe340812c28a1eddbba57d6dafac9e Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Wed, 28 Mar 2018 15:57:09 +0200 Subject: Add some TII testing scripts --- tii/README.md | 1 + tii/analyse_tii_from_file.py | 453 +++++++++++++++++++++++++++++++++++++++++++ tii/merge_iq.py | 53 +++++ tii/tii.ini.template | 51 +++++ tii/tii.py | 219 +++++++++++++++++++++ 5 files changed, 777 insertions(+) create mode 100644 tii/README.md create mode 100755 tii/analyse_tii_from_file.py create mode 100755 tii/merge_iq.py create mode 100644 tii/tii.ini.template create mode 100755 tii/tii.py (limited to 'tii') 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. -- cgit v1.2.3