Rework GUI and DPDCE

5 files changed, 894 insertions, 0 deletions

diff --git a/python/dpd/old/apply_adapt_dumps.py b/python/dpd/old/apply_adapt_dumps.py
new file mode 100755
index 0000000..20bc013
--- /dev/null
+++ b/python/dpd/old/apply_adapt_dumps.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# DPD Computation Engine, apply stored configuration.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import glob
+import logging
+
+dt = datetime.datetime.now().isoformat()
+logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s',
+                    datefmt='%Y-%m-%d %H:%M:%S',
+                    level=logging.DEBUG)
+
+import src.Adapt as Adapt
+import argparse
+
+parser = argparse.ArgumentParser(
+        description="Load pkl dumps DPD settings into ODR-DabMod")
+parser.add_argument('--port', default=50055, type=int,
+                    help='port of DPD server to connect to (default: 50055)',
+                    required=False)
+parser.add_argument('--rc-port', default=9400, type=int,
+                    help='port of ODR-DabMod ZMQ Remote Control to connect to (default: 9400)',
+                    required=False)
+parser.add_argument('--coefs', default='poly.coef',
+                    help='File with DPD coefficients, which will be read by ODR-DabMod',
+                    required=False)
+parser.add_argument('file', help='File to read the DPD settings from')
+
+cli_args = parser.parse_args()
+
+port = cli_args.port
+port_rc = cli_args.rc_port
+coef_path = cli_args.coefs
+filename = cli_args.file
+
+# No need to initialise a GlobalConfig since adapt only needs this one field
+class DummyConfig:
+    def __init__(self):
+        self.plot_location = None
+
+c = DummyConfig()
+
+adapt = Adapt.Adapt(c, port_rc, coef_path)
+
+print("Loading and applying DPD settings from {}".format(filename))
+adapt.load(filename)
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+# Copyright (c) 2017 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/python/dpd/old/iq_file_server.py b/python/dpd/old/iq_file_server.py
new file mode 100755
index 0000000..7a4e570
--- /dev/null
+++ b/python/dpd/old/iq_file_server.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This example server simulates the ODR-DabMod's
+# DPD server, taking samples from an IQ file
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import argparse
+import numpy as np
+from datetime import datetime
+
+SIZEOF_SAMPLE = 8 # complex floats
+# Constants for TM 1
+NbSymbols = 76
+NbCarriers = 1536
+Spacing = 2048
+NullSize = 2656
+SymSize = 2552
+FicSizeOut = 288
+FrameSize = NullSize + NbSymbols*SymSize
+
+def main():
+    parser = argparse.ArgumentParser(description="Simulate ODR-DabMod DPD server")
+    parser.add_argument('--port', default='50055',
+            help='port to listen on (default: 50055)',
+            required=False)
+    parser.add_argument('--file', help='I/Q File to read from (complex floats)',
+            required=True)
+    parser.add_argument('--samplerate', default='8192000', help='Sample rate',
+            required=False)
+
+    cli_args = parser.parse_args()
+
+    serve_file(cli_args)
+
+def recv_exact(sock, num_bytes):
+    bufs = []
+    while num_bytes > 0:
+        b = sock.recv(num_bytes)
+        if len(b) == 0:
+            break
+        num_bytes -= len(b)
+        bufs.append(b)
+    return b''.join(bufs)
+
+def serve_file(options):
+    oversampling = int(int(options.samplerate) / 2048000)
+    consumesamples = 8*FrameSize * oversampling
+    iq_data = np.fromfile(options.file, count=consumesamples, dtype=np.complex64)
+
+    print("Loaded {} samples of IQ data".format(len(iq_data)))
+
+    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+    s.bind(('localhost', int(options.port)))
+    s.listen()
+
+    try:
+        while True:
+            sock, addr_info = s.accept()
+            print("Got a connection from {}".format(addr_info))
+
+            ver = recv_exact(sock, 1)
+            (num_samps,) = struct.unpack("=I", recv_exact(sock, 4))
+            num_bytes = num_samps * SIZEOF_SAMPLE
+
+            if num_bytes > len(iq_data):
+                print("Truncating length to {} samples".format(len(iq_data)))
+                num_samps = len(iq_data)
+                num_bytes = num_samps * 4
+
+            tx_sec = datetime.now().timestamp()
+            tx_pps = int(16384000 * (tx_sec - int(tx_sec)))
+            tx_second = int(tx_sec)
+
+            # TX metadata and data
+            sock.sendall(struct.pack("=III", num_samps, tx_second, tx_pps))
+            sock.sendall(iq_data[-num_samps:].tobytes())
+
+            # RX metadata and data
+            rx_second = tx_second + 1
+            rx_pps = tx_pps
+            sock.sendall(struct.pack("=III", num_samps, rx_second, rx_pps))
+            sock.sendall(iq_data[-num_samps:].tobytes())
+
+            print("Sent {} samples".format(num_samps))
+
+            sock.close()
+    finally:
+        s.close()
+        raise
+
+main()
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 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/python/dpd/old/main.py b/python/dpd/old/main.py
new file mode 100755
index 0000000..9ea5a39
--- /dev/null
+++ b/python/dpd/old/main.py
@@ -0,0 +1,338 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# DPD Computation Engine standalone main file.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+"""This Python script is the main file for ODR-DabMod's DPD Computation Engine running
+in stand-alone mode.
+
+This engine calculates and updates the parameter of the digital
+predistortion module of ODR-DabMod."""
+
+import sys
+import datetime
+import os
+import argparse
+import matplotlib
+
+matplotlib.use('Agg')
+
+parser = argparse.ArgumentParser(
+    description="DPD Computation Engine for ODR-DabMod")
+parser.add_argument('--port', default=50055, type=int,
+                    help='port of DPD server to connect to (default: 50055)',
+                    required=False)
+parser.add_argument('--rc-port', default=9400, type=int,
+                    help='port of ODR-DabMod ZMQ Remote Control to connect to (default: 9400)',
+                    required=False)
+parser.add_argument('--samplerate', default=8192000, type=int,
+                    help='Sample rate',
+                    required=False)
+parser.add_argument('--coefs', default='poly.coef',
+                    help='File with DPD coefficients, which will be read by ODR-DabMod',
+                    required=False)
+parser.add_argument('--txgain', default=-1,
+                    help='TX Gain, -1 to leave unchanged',
+                    required=False,
+                    type=int)
+parser.add_argument('--rxgain', default=30,
+                    help='TX Gain, -1 to leave unchanged',
+                    required=False,
+                    type=int)
+parser.add_argument('--digital_gain', default=0.01,
+                    help='Digital Gain',
+                    required=False,
+                    type=float)
+parser.add_argument('--target_median', default=0.05,
+                    help='The target median for the RX and TX AGC',
+                    required=False,
+                    type=float)
+parser.add_argument('--samps', default='81920', type=int,
+                    help='Number of samples to request from ODR-DabMod',
+                    required=False)
+parser.add_argument('-i', '--iterations', default=10, type=int,
+                    help='Number of iterations to run',
+                    required=False)
+parser.add_argument('-L', '--lut',
+                    help='Use lookup table instead of polynomial predistorter',
+                    action="store_true")
+parser.add_argument('--enable-txgain-agc',
+                    help='Enable the TX gain AGC',
+                    action="store_true")
+parser.add_argument('--plot',
+                    help='Enable all plots, to be more selective choose plots in GlobalConfig.py',
+                    action="store_true")
+parser.add_argument('--name', default="", type=str,
+                    help='Name of the logging directory')
+parser.add_argument('-r', '--reset', action="store_true",
+                    help='Reset the DPD settings to the defaults.')
+parser.add_argument('-s', '--status', action="store_true",
+                    help='Display the currently running DPD settings.')
+parser.add_argument('--measure', action="store_true",
+                    help='Only measure metrics once')
+
+cli_args = parser.parse_args()
+
+port = cli_args.port
+port_rc = cli_args.rc_port
+coef_path = cli_args.coefs
+digital_gain = cli_args.digital_gain
+num_iter = cli_args.iterations
+rxgain = cli_args.rxgain
+txgain = cli_args.txgain
+name = cli_args.name
+plot = cli_args.plot
+
+# Logging
+import logging
+
+# Simple usage scenarios don't need to clutter /tmp
+if not (cli_args.status or cli_args.reset or cli_args.measure):
+    dt = datetime.datetime.now().isoformat()
+    logging_path = '/tmp/dpd_{}'.format(dt).replace('.', '_').replace(':', '-')
+    if name:
+        logging_path += '_' + name
+    print("Logs and plots written to {}".format(logging_path))
+    os.makedirs(logging_path)
+    logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s',
+                        datefmt='%Y-%m-%d %H:%M:%S',
+                        filename='{}/dpd.log'.format(logging_path),
+                        filemode='w',
+                        level=logging.DEBUG)
+    # also log up to INFO to console
+    console = logging.StreamHandler()
+    console.setLevel(logging.INFO)
+    # set a format which is simpler for console use
+    formatter = logging.Formatter('%(asctime)s - %(module)s - %(levelname)s - %(message)s')
+    # tell the handler to use this format
+    console.setFormatter(formatter)
+    # add the handler to the root logger
+    logging.getLogger('').addHandler(console)
+else:
+    dt = datetime.datetime.now().isoformat()
+    logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s',
+                        datefmt='%Y-%m-%d %H:%M:%S',
+                        level=logging.INFO)
+    logging_path = None
+
+logging.info("DPDCE starting up with options: {}".format(cli_args))
+
+import numpy as np
+import traceback
+from src.Model import Lut, Poly
+import src.Heuristics as Heuristics
+from src.Measure import Measure
+from src.ExtractStatistic import ExtractStatistic
+from src.Adapt import Adapt, dpddata_to_str
+from src.RX_Agc import Agc
+from src.TX_Agc import TX_Agc
+from src.Symbol_align import Symbol_align
+from src.GlobalConfig import GlobalConfig
+from src.MER import MER
+from src.Measure_Shoulders import Measure_Shoulders
+
+c = GlobalConfig(cli_args, logging_path)
+SA = Symbol_align(c)
+MER = MER(c)
+MS = Measure_Shoulders(c)
+meas = Measure(c, cli_args.samplerate, port, cli_args.samps)
+extStat = ExtractStatistic(c)
+adapt = Adapt(c, port_rc, coef_path)
+
+if cli_args.status:
+    txgain = adapt.get_txgain()
+    rxgain = adapt.get_rxgain()
+    digital_gain = adapt.get_digital_gain()
+    dpddata = dpddata_to_str(adapt.get_predistorter())
+
+    logging.info("ODR-DabMod currently running with TXGain {}, RXGain {}, digital gain {} and {}".format(
+        txgain, rxgain, digital_gain, dpddata))
+    sys.exit(0)
+
+if cli_args.lut:
+    model = Lut(c)
+else:
+    model = Poly(c)
+
+# Models have the default settings on startup
+adapt.set_predistorter(model.get_dpd_data())
+adapt.set_digital_gain(digital_gain)
+
+# Set RX Gain
+if rxgain == -1:
+    rxgain = adapt.get_rxgain()
+else:
+    adapt.set_rxgain(rxgain)
+
+# Set TX Gain
+if txgain == -1:
+    txgain = adapt.get_txgain()
+else:
+    adapt.set_txgain(txgain)
+
+tx_gain = adapt.get_txgain()
+rx_gain = adapt.get_rxgain()
+digital_gain = adapt.get_digital_gain()
+
+dpddata = adapt.get_predistorter()
+
+logging.info("TX gain {}, RX gain {}, digital_gain {}, {!s}".format(
+        tx_gain, rx_gain, digital_gain, dpddata_to_str(dpddata)))
+
+if cli_args.reset:
+    logging.info("DPD Settings were reset to default values.")
+    sys.exit(0)
+
+# Automatic Gain Control
+agc = Agc(meas, adapt, c)
+agc.run()
+
+if cli_args.measure:
+    txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples()
+
+    print("TX signal median {}".format(np.median(np.abs(txframe_aligned))))
+    print("RX signal median {}".format(rx_median))
+
+    tx, rx, phase_diff, n_per_bin = extStat.extract(txframe_aligned, rxframe_aligned)
+
+    off = SA.calc_offset(txframe_aligned)
+    print("off {}".format(off))
+    tx_mer = MER.calc_mer(txframe_aligned[off:off + c.T_U], debug_name='TX')
+    print("tx_mer {}".format(tx_mer))
+    rx_mer = MER.calc_mer(rxframe_aligned[off:off + c.T_U], debug_name='RX')
+    print("rx_mer {}".format(rx_mer))
+
+    mse = np.mean(np.abs((txframe_aligned - rxframe_aligned) ** 2))
+    print("mse {}".format(mse))
+
+    digital_gain = adapt.get_digital_gain()
+    print("digital_gain {}".format(digital_gain))
+
+    #rx_shoulder_tuple = MS.average_shoulders(rxframe_aligned)
+    #tx_shoulder_tuple = MS.average_shoulders(txframe_aligned)
+    sys.exit(0)
+
+# Disable TXGain AGC by default, so that the experiments are controlled
+# better.
+tx_agc = None
+if cli_args.enable_txgain_agc:
+    tx_agc = TX_Agc(adapt, c)
+
+state = 'report'
+i = 0
+lr = None
+n_meas = None
+while i < num_iter:
+    try:
+        # Measure
+        if state == 'measure':
+            # Get Samples and check gain
+            txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples()
+            if tx_agc and tx_agc.adapt_if_necessary(txframe_aligned):
+                continue
+
+            # Extract usable data from measurement
+            tx, rx, phase_diff, n_per_bin = extStat.extract(txframe_aligned, rxframe_aligned)
+
+            n_meas = Heuristics.get_n_meas(i)
+            if extStat.n_meas >= n_meas:  # Use as many measurements nr of runs
+                state = 'model'
+            else:
+                state = 'measure'
+
+        # Model
+        elif state == 'model':
+            # Calculate new model parameters and delete old measurements
+            if any([x is None for x in [tx, rx, phase_diff]]):
+                logging.error("No data to calculate model")
+                state = 'measure'
+                continue
+
+            lr = Heuristics.get_learning_rate(i)
+            model.train(tx, rx, phase_diff, lr=lr)
+            dpddata = model.get_dpd_data()
+            extStat = ExtractStatistic(c)
+            state = 'adapt'
+
+        # Adapt
+        elif state == 'adapt':
+            adapt.set_predistorter(dpddata)
+            state = 'report'
+
+        # Report
+        elif state == 'report':
+            try:
+                txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples()
+
+                # Store all settings for pre-distortion, tx and rx
+                adapt.dump()
+
+                # Collect logging data
+                off = SA.calc_offset(txframe_aligned)
+                tx_mer = MER.calc_mer(txframe_aligned[off:off + c.T_U], debug_name='TX')
+                rx_mer = MER.calc_mer(rxframe_aligned[off:off + c.T_U], debug_name='RX')
+                mse = np.mean(np.abs((txframe_aligned - rxframe_aligned) ** 2))
+                tx_gain = adapt.get_txgain()
+                rx_gain = adapt.get_rxgain()
+                digital_gain = adapt.get_digital_gain()
+                tx_median = np.median(np.abs(txframe_aligned))
+                rx_shoulder_tuple = MS.average_shoulders(rxframe_aligned)
+                tx_shoulder_tuple = MS.average_shoulders(txframe_aligned)
+
+                # Generic logging
+                logging.info(list((name, eval(name)) for name in
+                                  ['i', 'tx_mer', 'tx_shoulder_tuple', 'rx_mer',
+                                   'rx_shoulder_tuple', 'mse', 'tx_gain',
+                                   'digital_gain', 'rx_gain', 'rx_median',
+                                   'tx_median', 'lr', 'n_meas']))
+
+                # Model specific logging
+                if dpddata[0] == 'poly':
+                    coefs_am = dpddata[1]
+                    coefs_pm = dpddata[2]
+                    logging.info('It {}: coefs_am {}'.
+                                 format(i, coefs_am))
+                    logging.info('It {}: coefs_pm {}'.
+                                 format(i, coefs_pm))
+                elif dpddata[0] == 'lut':
+                    scalefactor = dpddata[1]
+                    lut = dpddata[2]
+                    logging.info('It {}: LUT scalefactor {}, LUT {}'.
+                                 format(i, scalefactor, lut))
+            except:
+                logging.error('Iteration {}: Report failed.'.format(i))
+                logging.error(traceback.format_exc())
+            i += 1
+            state = 'measure'
+
+    except:
+        logging.error('Iteration {} failed.'.format(i))
+        logging.error(traceback.format_exc())
+        i += 1
+        state = 'measure'
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+# Copyright (c) 2017 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/python/dpd/old/show_spectrum.py b/python/dpd/old/show_spectrum.py
new file mode 100755
index 0000000..f23dba2
--- /dev/null
+++ b/python/dpd/old/show_spectrum.py
@@ -0,0 +1,276 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This is an example tool that shows how to connect to ODR-DabMod's dpd TCP
+# server and get samples from there.
+#
+# Since the TX and RX samples are not perfectly aligned, the tool has to align
+# them properly, which is done in two steps: First on sample-level using a
+# correlation, then with subsample accuracy using a FFT approach.
+#
+# It requires SciPy and matplotlib.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import numpy as np
+import matplotlib.pyplot as pp
+from matplotlib.animation import FuncAnimation
+import argparse
+from scipy.misc import imsave
+
+SIZEOF_SAMPLE = 8 # complex floats
+
+# Constants for TM 1
+NbSymbols = 76
+NbCarriers = 1536
+Spacing = 2048
+NullSize = 2656
+SymSize = 2552
+FicSizeOut = 288
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot the spectrum of ODR-DabMod's DPD feedback")
+    parser.add_argument('--samps', default='10240', help='Number of samples to request at once',
+            required=False)
+    parser.add_argument('--port', default='50055',
+            help='port to connect to ODR-DabMod DPD (default: 50055)',
+            required=False)
+    parser.add_argument('--animated', action='store_true', help='Enable real-time animation')
+    parser.add_argument('--constellation', action='store_true', help='Draw constellaton plot')
+    parser.add_argument('--samplerate', default='8192000', help='Sample rate',
+            required=False)
+
+    cli_args = parser.parse_args()
+
+    if cli_args.constellation:
+        plot_constellation_once(cli_args)
+    elif cli_args.animated:
+        plot_spectrum_animated(cli_args)
+    else:
+        plot_spectrum_once(cli_args)
+
+def recv_exact(sock, num_bytes):
+    bufs = []
+    while num_bytes > 0:
+        b = sock.recv(num_bytes)
+        if len(b) == 0:
+            break
+        num_bytes -= len(b)
+        bufs.append(b)
+    return b''.join(bufs)
+
+def get_samples(port, num_samps_to_request):
+    """Connect to ODR-DabMod, retrieve TX and RX samples, load
+    into numpy arrays, and return a tuple
+    (tx_timestamp, tx_samples, rx_timestamp, rx_samples)
+    where the timestamps are doubles, and the samples are numpy
+    arrays of complex floats, both having the same size
+    """
+
+    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+    s.connect(('localhost', port))
+
+    print("Send version");
+    s.sendall(b"\x01")
+
+    print("Send request for {} samples".format(num_samps_to_request))
+    s.sendall(struct.pack("=I", num_samps_to_request))
+
+    print("Wait for TX metadata")
+    num_samps, tx_second, tx_pps = struct.unpack("=III", recv_exact(s, 12))
+    tx_ts = tx_second + tx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} TX samples".format(num_samps))
+        txframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        txframe = np.fromstring(txframe_bytes, dtype=np.complex64)
+    else:
+        txframe = np.array([], dtype=np.complex64)
+
+
+    print("Wait for RX metadata")
+    rx_second, rx_pps = struct.unpack("=II", recv_exact(s, 8))
+    rx_ts = rx_second + rx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} RX samples".format(num_samps))
+        rxframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64)
+    else:
+        rxframe = np.array([], dtype=np.complex64)
+
+    print("Disconnecting")
+    s.close()
+
+    return (tx_ts, txframe, rx_ts, rxframe)
+
+def recv_rxtx(port, num_samps_to_request):
+    tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request)
+
+    # convert to complex doubles for more dynamic range
+    txframe = txframe.astype(np.complex128)
+    rxframe = rxframe.astype(np.complex128)
+
+    print("Received {} & {} frames at {} and {}".format(
+        len(txframe), len(rxframe), tx_ts, rx_ts))
+    return tx_ts, txframe, rx_ts, rxframe
+
+def get_spectrum(port, num_samps_to_request):
+    tx_ts, txframe, rx_ts, rxframe = recv_rxtx(port, num_samps_to_request)
+
+    print("Calculate TX and RX spectrum assuming 8192000 samples per second")
+    tx_spectrum = np.fft.fftshift(np.fft.fft(txframe, fft_size))
+    tx_power = 20*np.log10(np.abs(tx_spectrum))
+
+    rx_spectrum = np.fft.fftshift(np.fft.fft(rxframe, fft_size))
+    rx_power = 20*np.log10(np.abs(rx_spectrum))
+    return tx_power, rx_power
+
+def remove_guard_intervals(frame, options):
+    """Remove the cyclic prefix. The frame needs to be aligned to the
+    end of the transmission frame. Transmission Mode 1 is assumed"""
+    oversample = int(int(options.samplerate) / 2048000)
+
+    # From the end, take 2048 samples, then skip 504 samples
+    frame = frame[::-1]
+
+    stride_len = Spacing * oversample
+    stride_advance = SymSize * oversample
+
+    # Truncate the frame to an integer length of strides
+    newlen = len(frame) - (len(frame) % stride_advance)
+    print("Truncating frame from {} to {}".format(len(frame), newlen))
+    frame = frame[:newlen]
+
+    # Remove the cyclic prefix
+    frame = frame.reshape(-1, stride_advance)[:,:stride_len].reshape(-1)
+
+    # Reverse again
+    return frame[::-1]
+
+
+def plot_constellation_once(options):
+    port = int(options.port)
+    num_samps_to_request = int(options.samps)
+
+    tx_ts, txframe, rx_ts, rxframe = recv_rxtx(port, num_samps_to_request)
+
+    frame = remove_guard_intervals(txframe, options)
+
+    oversample = int(int(options.samplerate) / 2048000)
+
+    n = Spacing * oversample # is also number of samples per symbol
+    if len(frame) % n != 0:
+        raise ValueError("Frame length doesn't contain exact number of symbols")
+    num_syms = int(len(frame) / n)
+    print("frame {} has {} symbols".format(len(frame), num_syms))
+    spectrums = np.array([np.fft.fftshift(np.fft.fft(frame[n*i:n*(i+1)], n)) for i in range(num_syms)])
+
+    def normalise(x):
+        """Normalise a real-valued array x to the range [0,1]"""
+        y = x + np.min(x)
+        return x / np.max(x)
+
+    imsave("spectrums.png", np.concatenate([
+        normalise(np.abs(spectrums)),
+        normalise(np.angle(spectrums))]))
+
+    # Only take bins that are supposed to contain energy
+    # i.e. the middle 1536 bins, excluding the bin at n/2
+    assert(n % 2 == 0)
+    n_half = int(n/2)
+    spectrums = np.concatenate(
+            [spectrums[...,n_half-768:n_half],
+             spectrums[...,n_half + 1:n_half + 769]], axis=1)
+
+    sym_indices = (np.tile(np.arange(num_syms-1).reshape(num_syms-1,1), (1,NbCarriers)) +
+                   np.tile(np.linspace(-0.4, 0.4, NbCarriers), (num_syms-1, 1) ) )
+    sym_indices = sym_indices.reshape(-1)
+    diff_angles = np.mod(np.diff(np.angle(spectrums, deg=1), axis=0), 360)
+    #sym_points = spectrums[:-1].reshape(-1)
+    # Set amplitude and phase of low power points to zero, avoid cluttering diagram
+    #sym_points[np.abs(sym_points) < np.mean(np.abs(sym_points)) * 0.1] = 0
+
+    print("ix {}  spec {} da {}".format(
+        sym_indices.shape, spectrums.shape, diff_angles.shape))
+
+    fig = pp.figure()
+
+    fig.suptitle("Constellation")
+    ax1 = fig.add_subplot(111)
+    ax1.set_title("TX")
+    ax1.scatter(sym_indices, diff_angles.reshape(-1), alpha=0.1)
+
+    pp.show()
+
+fft_size = 4096
+
+def plot_spectrum_once(options):
+    port = int(options.port)
+    num_samps_to_request = int(options.samps)
+    freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./int(options.samplerate)))
+
+    tx_power, rx_power = get_spectrum(port, num_samps_to_request)
+    fig = pp.figure()
+
+    fig.suptitle("TX and RX spectrum")
+    ax1 = fig.add_subplot(211)
+    ax1.set_title("TX")
+    ax1.plot(freqs, tx_power, 'r')
+    ax2 = fig.add_subplot(212)
+    ax2.set_title("RX")
+    ax2.plot(freqs, rx_power, 'b')
+    pp.show()
+
+def plot_spectrum_animated(options):
+    port = int(options.port)
+    num_samps_to_request = int(options.samps)
+    freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./int(options.samplerate)))
+
+    fig, axes = pp.subplots(2, sharex=True)
+    line1, = axes[0].plot(freqs, np.ones(len(freqs)), 'r', animated=True)
+    axes[0].set_title("TX")
+    line2, = axes[1].plot(freqs, np.ones(len(freqs)), 'b', animated=True)
+    axes[1].set_title("RX")
+    lines = [line1, line2]
+
+    axes[0].set_ylim(-30, 50)
+    axes[1].set_ylim(-60, 40)
+
+    def update(frame):
+        tx_power, rx_power = get_spectrum(port, num_samps_to_request)
+
+        lines[0].set_ydata(tx_power)
+        lines[1].set_ydata(rx_power)
+        return lines
+
+    ani = FuncAnimation(fig, update, blit=True)
+    pp.show()
+
+main()
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 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/python/dpd/old/store_received.py b/python/dpd/old/store_received.py
new file mode 100755
index 0000000..19b735e
--- /dev/null
+++ b/python/dpd/old/store_received.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This is an example tool that shows how to connect to ODR-DabMod's dpd TCP server
+# and get samples from there.
+#
+# Since the TX and RX samples are not perfectly aligned, the tool has to align them properly,
+# which is done in two steps: First on sample-level using a correlation, then with subsample
+# accuracy using a FFT approach.
+#
+# It requires SciPy and matplotlib.
+#
+# Copyright (C) 2017 Matthias P. Braendli
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import argparse
+import os
+import time
+from src.GlobalConfig import GlobalConfig
+from src.Measure import Measure
+
+SIZEOF_SAMPLE = 8 # complex floats
+
+parser = argparse.ArgumentParser(description="Plot the spectrum of ODR-DabMod's DPD feedback")
+parser.add_argument('--samps', default='10240', type=int,
+        help='Number of samples to request at once',
+        required=False)
+parser.add_argument('--port', default='50055', type=int,
+        help='port to connect to ODR-DabMod DPD (default: 50055)',
+        required=False)
+parser.add_argument('--count', default='1', type=int,
+        help='Number of recordings',
+        required=False)
+parser.add_argument('--verbose', type=int, default=0,
+        help='Level of verbosity',
+        required=False)
+parser.add_argument('--plot',
+                    help='Enable all plots, to be more selective choose plots in GlobalConfig.py',
+                    action="store_true")
+parser.add_argument('--samplerate', default=8192000, type=int,
+                    help='Sample rate',
+                    required=False)
+
+cli_args = parser.parse_args()
+
+cli_args.target_median = 0.05
+
+c = GlobalConfig(cli_args, None)
+
+meas = Measure(c, cli_args.samplerate, cli_args.port, cli_args.samps)
+
+for i in range(int(cli_args.count)):
+    if i>0:
+        time.sleep(0.1)
+
+    txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples()
+
+    txframe_aligned.tofile("%d_tx_record.iq" % i)
+    rxframe_aligned.tofile("%d_rx_record.iq" % i)
+
+# 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.