From d5cbe10c0e2298b0e40161607a3da158249bdb82 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Tue, 4 Dec 2018 10:18:33 +0100 Subject: Move python stuff to folder --- python/dpd/README.md | 264 +++++++++++++++++ python/dpd/apply_adapt_dumps.py | 75 +++++ python/dpd/dpd.ini | 60 ++++ python/dpd/img/setup_diagram.svg | 395 +++++++++++++++++++++++++ python/dpd/img/setup_photo.svg | 184 ++++++++++++ python/dpd/img/shoulder_measurement_after.png | Bin 0 -> 368619 bytes python/dpd/img/shoulder_measurement_before.png | Bin 0 -> 381174 bytes python/dpd/iq_file_server.py | 120 ++++++++ python/dpd/lut.coef | 64 ++++ python/dpd/main.py | 338 +++++++++++++++++++++ python/dpd/poly.coef | 12 + python/dpd/show_spectrum.py | 276 +++++++++++++++++ python/dpd/src/Adapt.py | 286 ++++++++++++++++++ python/dpd/src/Dab_Util.py | 246 +++++++++++++++ python/dpd/src/ExtractStatistic.py | 196 ++++++++++++ python/dpd/src/GlobalConfig.py | 108 +++++++ python/dpd/src/Heuristics.py | 56 ++++ python/dpd/src/MER.py | 132 +++++++++ python/dpd/src/Measure.py | 136 +++++++++ python/dpd/src/Measure_Shoulders.py | 158 ++++++++++ python/dpd/src/Model.py | 32 ++ python/dpd/src/Model_AM.py | 122 ++++++++ python/dpd/src/Model_Lut.py | 60 ++++ python/dpd/src/Model_PM.py | 124 ++++++++ python/dpd/src/Model_Poly.py | 101 +++++++ python/dpd/src/RX_Agc.py | 166 +++++++++++ python/dpd/src/Symbol_align.py | 193 ++++++++++++ python/dpd/src/TX_Agc.py | 131 ++++++++ python/dpd/src/__init__.py | 0 python/dpd/src/phase_align.py | 98 ++++++ python/dpd/src/subsample_align.py | 111 +++++++ python/dpd/store_received.py | 85 ++++++ 32 files changed, 4329 insertions(+) create mode 100644 python/dpd/README.md create mode 100755 python/dpd/apply_adapt_dumps.py create mode 100644 python/dpd/dpd.ini create mode 100644 python/dpd/img/setup_diagram.svg create mode 100644 python/dpd/img/setup_photo.svg create mode 100644 python/dpd/img/shoulder_measurement_after.png create mode 100644 python/dpd/img/shoulder_measurement_before.png create mode 100755 python/dpd/iq_file_server.py create mode 100644 python/dpd/lut.coef create mode 100755 python/dpd/main.py create mode 100644 python/dpd/poly.coef create mode 100755 python/dpd/show_spectrum.py create mode 100644 python/dpd/src/Adapt.py create mode 100644 python/dpd/src/Dab_Util.py create mode 100644 python/dpd/src/ExtractStatistic.py create mode 100644 python/dpd/src/GlobalConfig.py create mode 100644 python/dpd/src/Heuristics.py create mode 100644 python/dpd/src/MER.py create mode 100644 python/dpd/src/Measure.py create mode 100644 python/dpd/src/Measure_Shoulders.py create mode 100644 python/dpd/src/Model.py create mode 100644 python/dpd/src/Model_AM.py create mode 100644 python/dpd/src/Model_Lut.py create mode 100644 python/dpd/src/Model_PM.py create mode 100644 python/dpd/src/Model_Poly.py create mode 100644 python/dpd/src/RX_Agc.py create mode 100644 python/dpd/src/Symbol_align.py create mode 100644 python/dpd/src/TX_Agc.py create mode 100644 python/dpd/src/__init__.py create mode 100644 python/dpd/src/phase_align.py create mode 100755 python/dpd/src/subsample_align.py create mode 100755 python/dpd/store_received.py (limited to 'python/dpd') diff --git a/python/dpd/README.md b/python/dpd/README.md new file mode 100644 index 0000000..307a2f5 --- /dev/null +++ b/python/dpd/README.md @@ -0,0 +1,264 @@ +Digital Predistortion Computation Engine for ODR-DabMod +======================================================= + +This folder contains a digital predistortion prototype. +It was only tested in a laboratory system, and is not ready +for production usage. + +Concept +------- + +ODR-DabMod makes outgoing TX samples and feedback RX samples available to an +external tool. This external tool can request a buffer of samples for analysis, +can calculate coefficients for the predistorter in ODR-DabMod and load the new +coefficients using the remote control. + +The external tool is called the Digital Predistortion Computation Engine (DPDCE). +The DPDCE is written in python, and makes use of the numpy library for +efficient computation. Its sources reside in the *dpd* folder. + +The predistorter in ODR-DabMod supports two modes: polynomial and lookup table. +In the DPDCE, only the polynomial model is implemented at the moment. + +The *dpd/main.py* script is the entry point for the *DPD Computation Engine* +into which these features will be implemented. The tool uses modules from the +*dpd/src/* folder: + +- Sample transfer and time alignment with subsample accuracy is done by *Measure.py* +- Estimating the effects of the PA using some model and calculation of the updated + polynomial coefficients is done in *Model.py* and other specific *Model_XXX.py* files +- Finally, *Adapt.py* updates the ODR-DabMod predistortion setting and digital gain + +These modules themselves use additional helper scripts in the *dpd/src/* folder. + +Requirements +------------ + +- USRP B200. +- Power amplifier. +- A feedback connection from the power amplifier output, such that the average power level at + the USRP RX port is at -45dBm or lower. + Usually this is done with a directional coupler and additional attenuators. +- ODR-DabMod with enabled *dpd_port*, and with a samplerate of 8192000 samples per second. +- Synchronous=1 so that the USRP has the timestamping set properly, internal refclk and pps + are sufficient (not GPSDO necessary). +- A live mux source with TIST enabled. + +See dpd/dpd.ini for an example. + +The DPD server port can be tested with the *dpd/show_spectrum.py* helper tool, which can also display +a constellation diagram. + +Hardware Setup +-------------- + +![setup diagram](img/setup_diagram.svg) +![setup photo](img/setup_photo.svg) + +Our setup is depicted in the Figure above. We used components with the following properties: + 1. USRP TX (max +20dBm) + 2. Band III Filter (Mini-Circuits RBP-220W+, 190-250MHz, -3.5dB) + 3. Power amplifier (Mini-Circuits, max +15dBm output, +10 dB gain at 200MHz) + 4. Directional coupler (approx. -25dB @ 223MHz) + 5. Attenuator (-20 dB) + 6. Attenuator (-30 dB) + 7. USRP RX (max -15dBm input allowed, max -45dBm desired) + 8. Spectrum analyzer (max +30dBm allowed) + +It is important to make sure that the USRP RX port does not receive too much +power. Otherwise the USRP will break. Here is an example of how we calculated +the maximal USRP RX input power for our case. As this is only a rough +calculation to protect the port, the predistortion software will later +automatically apply a normalization for the RX input by adapting the USRP RX +gain. + + TX Power + PA Gain - Coupling Factor - Attenuation = 20dBm + 10dB -25dB -50dB = -45dBm + +Thus we have a margin of about 30dB for the input power of the USRP RX port. +Keep in mind we need to calculate using peak power, not average power, and it is +essential that there is no nonlinearity in the RX path! + +Software Setup +-------------- + +We assume that you already installed *ODR-DabMux* and *ODR-DabMod*. +You should install the required python dependencies for the DPDCE using +distribution packages. You will need at least scipy, matplotlib and +python-zeromq. + +Use the predistortion +---------------------- + +Make sure you have a ODR-DabMux running with a TCP output on port 9200. + +Then run the modulator, with the example dpd configuration file. + +``` +./odr-dabmod dpd/dpd.ini +``` + +This configuration file is different from usual defaults in several respects: + + * logging to /tmp/dabmod.log + * 4x oversampling: 8192000 sample rate + * a very small digital gain, which will be overridden by the DPDCE + * predistorter enabled + +The TX gain should be chosen so that you can drive your amplifier into +saturation with a digital gain of 0.1, so that there is margin for the DPD to +operate. + +You should *not modify txgain, rxgain, digital gain or coefficient settings manually!* +When the DPDCE is used, it controls these settings, and there are command line +options for you to define initial values. + +To verify that the communication between the DPDCE and ODR-DabMod is ok, +you can use the status and reset options: + +``` +cd dpd +python main.py --status +python main.py --reset +``` + +The reset option sets all DPD-related settings to the defaults (as shown in the +`--help` usage screen) and stops. + +When neither `--status` nor `--reset` is given, the DPDCE will run the predistortion +algorithm. As a first test you should run the DPDCE with the `--plot` +parameter. It preserves the output power and generates all available +visualisation plots in the newly created logging directory +`/tmp/dpd_`. As the predistortion should increase the peak to +shoulder ratio, you should select a *txgain* in the ODR-DabMod configuration +file such that the initial peak-to-soulder ratio visible on your spectrum +analyser. This way, you will be able to see a the +change. + +``` +cd dpd +python main.py --plot +``` + +The DPDCE now does 10 iterations, and tries to improve the predistortion effectiveness. +In each step the learning rate is decreased. The learning rate is the factor +with which new coefficients are weighted in a weighted mean with the old +coefficients. Moreover the nuber of measurements increases in each iteration. +You find more information about that in *Heuristic.py*. + +Each plot is stored to the logging directory under a filename containing its +time stamp and its label. Following plots are generated chronologically: + + - ExtractStatistic: Extracted information from one or multiple measurements. + - Model\_AM: Fitted function for the amplitudes of the power amplifier against the TX amplitude. + - Model\_PM: Fitted function for the phase difference of the power amplifier against the TX amplitude. + - adapt.pkl: Contains all settings for the predistortion. + You can load them again without doing measurements with the `apply_adapt_dumps.py` script. + - MER: Constellation diagram used to calculate the modulation error rate. + +After the run you should be able to observe that the peak-to-shoulder +difference has increased on your spectrum analyzer, similar to the figure below. + +Without digital predistortion: + +![shoulder_measurement_before](img/shoulder_measurement_before.png) + +With digital predistortion, computed by the DPDCE: + +![shoulder_measurement_after](img/shoulder_measurement_after.png) + +Now see what happens if you apply the predistortions for different TX gains. +You can either set the TX gain before you start the predistortion or using the +command line option `--txgain gain`. You can also try to adjust other +parameters. To see their documentation run `python main.py --help`. + +File format for coefficients +---------------------------- +The coef file contains the polynomial coefficients used in the predistorter. +The file format is very similar to the filtertaps file used in the FIR filter. +It is a text-based format that can easily be inspected and edited in a text +editor. + +The first line contains an integer that defines the predistorter to be used: +1 for polynomial, 2 for lookup table. + +For the polynomial, the subsequent line contains the number of coefficients +as an integer. The second and third lines contain the real, respectively the +imaginary parts of the first coefficient. Fourth and fifth lines give the +second coefficient, and so on. The file therefore contains 1 + 1 + 2xN lines if +it contains N coefficients. + +For the lookup table, the subsequent line contains a float scalefactor that is +applied to the samples in order to bring them into the range of 32-bit unsigned +integer. Then, the next pair of lines contains real and imaginary part of the first +lookup-table entry, which is multiplied to samples in first range. Then it's +followed by 31 other pairs. The entries are complex values close to 1 + 0j. +The file therefore contains 1 + 1 + 2xN lines if it contains N coefficients. + +TODO +---- + + - Understand and fix occasional ODR-DabMod crashes when using DPDCE. + - Make the predistortion more robust. At the moment the shoulders sometimes + increase instead of decrease after applying newly calculated predistortion + parameters. Can this behaviour be predicted from the measurement? This would + make it possible to filter out bad predistortion settings. + - Find a better measurement for the quality of the predistortion. The USRP + might not be good enough to measure large peak-to-shoulder ratios, because + the ADC has 12 bits and DAB signals have a large crest factor. + - Implement a Volterra polynomial to model the PA. Compared to the current + model this would also capture the time dependent behaviour of the PA (memory + effects). + - Continuously observe DAB signal in frequency domain and make sure the power + stays the same. At the moment only the power in the time domain is kept the + same. + - At the moment we assume that the USRP RX gain has to be larger than 30dB and + the received signal should have a median absolute value of 0.05 in order to + have a high quality quantization. Do measurements to support or improve + this heuristic. + - Check if we need to measure MER differently (average over more symbols?) + - Is -45dBm the best RX feedback power level? + +REFERENCES +---------- + +Some papers: + +The paper Raich, Qian, Zhou, "Orthogonal Polynomials for Power Amplifier +Modeling and Predistorter Design" proposes other base polynomials that have +less numerical instability. + +Aladrén, Garcia, Carro, de Mingo, and Sanchez-Perez, "Digital Predistortion +Based on Zernike Polynomial Functions for RF Nonlinear Power Amplifiers". + +Jiang and Wilford, "Digital predistortion for power amplifiers using separable functions" + +Changsoo Eun and Edward J. Powers, "A New Volterra Predistorter Based on the Indirect Learning Architecture" + +Raviv Raich, Hua Qian, and G. Tong Zhou, "Orthogonal Polynomials for Power Amplifier Modeling and Predistorter Design" + + +Models without memory: + +Complex polynomial: y[i] = a1 x[i] + a2 x[i]^2 + a3 x[i]^3 + ... + +The complex polynomial corresponds to the input/output relationship that +applies to the PA in passband (real-valued signal). According to several +sources, this gets transformed to another representation if we consider complex +baseband instead. In the following, all variables are complex. + +Odd-order baseband: y[i] = (b1 + b2 abs(x[i])^2 + b3 abs(x[i])^4) + ...) x[i] + +Complete baseband: y[i] = (b1 + b2 abs(x[i]) + b3 abs(x[i])^2) + ...) x[i] + +with + b_k = 2^{1-k} \binom{k}{(k-1)/2} a_k + + +Models with memory: + + - Hammerstein model: Nonlinearity followed by LTI filter + - Wiener model: LTI filter followed by NL + - Parallel Wiener: input goes to N delays, each delay goes to a NL, all NL outputs summed. + +Taken from slide 36 of [ECE218C Lecture 15](http://www.ece.ucsb.edu/Faculty/rodwell/Classes/ece218c/notes/Lecture15_Digital%20Predistortion_and_Future%20Challenges.pdf) + diff --git a/python/dpd/apply_adapt_dumps.py b/python/dpd/apply_adapt_dumps.py new file mode 100755 index 0000000..20bc013 --- /dev/null +++ b/python/dpd/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/dpd.ini b/python/dpd/dpd.ini new file mode 100644 index 0000000..31d6140 --- /dev/null +++ b/python/dpd/dpd.ini @@ -0,0 +1,60 @@ +[remotecontrol] +telnet=1 +telnetport=2121 +zmqctrl=1 +zmqctrlendpoint=tcp://127.0.0.1:9400 + +[log] +syslog=0 +filelog=1 +filename=/tmp/dabmod.log + +[input] +transport=tcp +source=localhost:9200 + +[modulator] +gainmode=var +rate=8192000 + +# keep in mind that the DPDCE will set the digital gain through the RC! +digital_gain=0.001 + +[firfilter] +enabled=1 + +[poly] +enabled=1 +polycoeffile=dpd/poly.coef + +# How many threads to use for the predistorter. +# If not set, detect automatically. +#num_threads=2 + +[output] +# to prepare a file for the dpd/iq_file_server.py script, +# use output=file +output=uhd + +[fileoutput] +filename=dpd.iq + +[uhdoutput] +device= +master_clock_rate=32768000 +type=b200 +txgain=80 +channel=13C +refclk_source=internal +pps_source=none +behaviour_refclk_lock_lost=ignore +max_gps_holdover_time=600 +dpd_port=50055 +rxgain=15 + +[delaymanagement] +; Use synchronous=1 so that the USRP time is set. This works +; even in the absence of a reference clk and PPS +synchronous=1 +mutenotimestamps=1 +offset=4.0 diff --git a/python/dpd/img/setup_diagram.svg b/python/dpd/img/setup_diagram.svg new file mode 100644 index 0000000..423ee1e --- /dev/null +++ b/python/dpd/img/setup_diagram.svg @@ -0,0 +1,395 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + ① USRP TX + + + ① USRP TX + + + + + + + + + + ② Filter + + + ② Filter + + + + + + + + + + ③ PA + + + ③ PA + + + + + + + ④ Dir. Coupler + + + ④ Dir. Coupler + + + + + + + + + + + + + + + + ⑤,⑥ Attenuators + + + ⑤,⑥Attenuators + + + + + + + + ⑦ USRP RX + + + ⑦ USRP RX + + + + + + + + + ⑧ Spectrum Analyzer + + + ⑧ Spectrum Analyzer + + + + diff --git a/python/dpd/img/setup_photo.svg b/python/dpd/img/setup_photo.svg new file mode 100644 index 0000000..a401fda --- /dev/null +++ b/python/dpd/img/setup_photo.svg @@ -0,0 +1,184 @@ + + + + + + + + image/svg+xml + + + + + + + + + ① + ② + + + ③ + ④ + ⑤⑥ + ⑦ + ⑧ + diff --git a/python/dpd/img/shoulder_measurement_after.png b/python/dpd/img/shoulder_measurement_after.png new file mode 100644 index 0000000..2631256 Binary files /dev/null and b/python/dpd/img/shoulder_measurement_after.png differ diff --git a/python/dpd/img/shoulder_measurement_before.png b/python/dpd/img/shoulder_measurement_before.png new file mode 100644 index 0000000..3581a72 Binary files /dev/null and b/python/dpd/img/shoulder_measurement_before.png differ diff --git a/python/dpd/iq_file_server.py b/python/dpd/iq_file_server.py new file mode 100755 index 0000000..7a4e570 --- /dev/null +++ b/python/dpd/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/lut.coef b/python/dpd/lut.coef new file mode 100644 index 0000000..a198d56 --- /dev/null +++ b/python/dpd/lut.coef @@ -0,0 +1,64 @@ +2 +4294967296 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 +1 +0 diff --git a/python/dpd/main.py b/python/dpd/main.py new file mode 100755 index 0000000..9ea5a39 --- /dev/null +++ b/python/dpd/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/poly.coef b/python/dpd/poly.coef new file mode 100644 index 0000000..248d316 --- /dev/null +++ b/python/dpd/poly.coef @@ -0,0 +1,12 @@ +1 +5 +1.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 diff --git a/python/dpd/show_spectrum.py b/python/dpd/show_spectrum.py new file mode 100755 index 0000000..f23dba2 --- /dev/null +++ b/python/dpd/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/src/Adapt.py b/python/dpd/src/Adapt.py new file mode 100644 index 0000000..a57602f --- /dev/null +++ b/python/dpd/src/Adapt.py @@ -0,0 +1,286 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine: updates ODR-DabMod's settings +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file +""" +This module is used to change settings of ODR-DabMod using +the ZMQ remote control socket. +""" + +import zmq +import logging +import numpy as np +import os +import datetime +import pickle + +LUT_LEN = 32 +FORMAT_POLY = 1 +FORMAT_LUT = 2 + + +def _write_poly_coef_file(coefs_am, coefs_pm, path): + assert (len(coefs_am) == len(coefs_pm)) + + f = open(path, 'w') + f.write("{}\n{}\n".format(FORMAT_POLY, len(coefs_am))) + for coef in coefs_am: + f.write("{}\n".format(coef)) + for coef in coefs_pm: + f.write("{}\n".format(coef)) + f.close() + + +def _write_lut_file(scalefactor, lut, path): + assert (len(lut) == LUT_LEN) + + f = open(path, 'w') + f.write("{}\n{}\n".format(FORMAT_LUT, scalefactor)) + for coef in lut: + f.write("{}\n{}\n".format(coef.real, coef.imag)) + f.close() + +def dpddata_to_str(dpddata): + if dpddata[0] == "poly": + coefs_am = dpddata[1] + coefs_pm = dpddata[2] + return "dpd_coefs_am {}, dpd_coefs_pm {}".format( + coefs_am, coefs_pm) + elif dpddata[0] == "lut": + scalefactor = dpddata[1] + lut = dpddata[2] + return "LUT scalefactor {}, LUT {}".format( + scalefactor, lut) + else: + raise ValueError("Unknown dpddata type {}".format(dpddata[0])) + +class Adapt: + """Uses the ZMQ remote control to change parameters of the DabMod + + Parameters + ---------- + port : int + Port at which the ODR-DabMod is listening to connect the + ZMQ remote control. + """ + + def __init__(self, config, port, coef_path): + logging.debug("Instantiate Adapt object") + self.c = config + self.port = port + self.coef_path = coef_path + self.host = "localhost" + self._context = zmq.Context() + + def _connect(self): + """Establish the connection to ODR-DabMod using + a ZMQ socket that is in request mode (Client). + Returns a socket""" + sock = self._context.socket(zmq.REQ) + poller = zmq.Poller() + poller.register(sock, zmq.POLLIN) + + sock.connect("tcp://%s:%d" % (self.host, self.port)) + + sock.send(b"ping") + + socks = dict(poller.poll(1000)) + if socks: + if socks.get(sock) == zmq.POLLIN: + data = [el.decode() for el in sock.recv_multipart()] + + if data != ['ok']: + raise RuntimeError( + "Invalid ZMQ RC answer to 'ping' at %s %d: %s" % + (self.host, self.port, data)) + else: + sock.close(linger=10) + raise RuntimeError( + "ZMQ RC does not respond to 'ping' at %s %d" % + (self.host, self.port)) + + return sock + + def send_receive(self, message): + """Send a message to ODR-DabMod. It always + returns the answer ODR-DabMod sends back. + + An example message could be + "get sdr txgain" or "set sdr txgain 50" + + Parameter + --------- + message : str + The message string that will be sent to the receiver. + """ + sock = self._connect() + logging.debug("Send message: %s" % message) + msg_parts = message.split(" ") + for i, part in enumerate(msg_parts): + if i == len(msg_parts) - 1: + f = 0 + else: + f = zmq.SNDMORE + + sock.send(part.encode(), flags=f) + + data = [el.decode() for el in sock.recv_multipart()] + logging.debug("Received message: %s" % message) + return data + + def set_txgain(self, gain): + """Set a new txgain for the ODR-DabMod. + + Parameters + ---------- + gain : int + new TX gain, in the same format as ODR-DabMod's config file + """ + # TODO this is specific to the B200 + if gain < 0 or gain > 89: + raise ValueError("Gain has to be in [0,89]") + return self.send_receive("set sdr txgain %.4f" % float(gain)) + + def get_txgain(self): + """Get the txgain value in dB for the ODR-DabMod.""" + # TODO handle failure + return float(self.send_receive("get sdr txgain")[0]) + + def set_rxgain(self, gain): + """Set a new rxgain for the ODR-DabMod. + + Parameters + ---------- + gain : int + new RX gain, in the same format as ODR-DabMod's config file + """ + # TODO this is specific to the B200 + if gain < 0 or gain > 89: + raise ValueError("Gain has to be in [0,89]") + return self.send_receive("set sdr rxgain %.4f" % float(gain)) + + def get_rxgain(self): + """Get the rxgain value in dB for the ODR-DabMod.""" + # TODO handle failure + return float(self.send_receive("get sdr rxgain")[0]) + + def set_digital_gain(self, gain): + """Set a new rxgain for the ODR-DabMod. + + Parameters + ---------- + gain : int + new RX gain, in the same format as ODR-DabMod's config file + """ + msg = "set gain digital %.5f" % gain + return self.send_receive(msg) + + def get_digital_gain(self): + """Get the rxgain value in dB for the ODR-DabMod.""" + # TODO handle failure + return float(self.send_receive("get gain digital")[0]) + + def get_predistorter(self): + """Load the coefficients from the file in the format given in the README, + return ("poly", [AM coef], [PM coef]) or ("lut", scalefactor, [LUT entries]) + """ + f = open(self.coef_path, 'r') + lines = f.readlines() + predistorter_format = int(lines[0]) + if predistorter_format == FORMAT_POLY: + coefs_am_out = [] + coefs_pm_out = [] + n_coefs = int(lines[1]) + coefs = [float(l) for l in lines[2:]] + i = 0 + for c in coefs: + if i < n_coefs: + coefs_am_out.append(c) + elif i < 2 * n_coefs: + coefs_pm_out.append(c) + else: + raise ValueError( + 'Incorrect coef file format: too many' + ' coefficients in {}, should be {}, coefs are {}' + .format(self.coef_path, n_coefs, coefs)) + i += 1 + f.close() + return 'poly', coefs_am_out, coefs_pm_out + elif predistorter_format == FORMAT_LUT: + scalefactor = int(lines[1]) + coefs = np.array([float(l) for l in lines[2:]], dtype=np.float32) + coefs = coefs.reshape((-1, 2)) + lut = coefs[..., 0] + 1j * coefs[..., 1] + if len(lut) != LUT_LEN: + raise ValueError("Incorrect number of LUT entries ({} expected {})".format(len(lut), LUT_LEN)) + return 'lut', scalefactor, lut + else: + raise ValueError("Unknown predistorter format {}".format(predistorter_format)) + + def set_predistorter(self, dpddata): + """Update the predistorter data in the modulator. Takes the same + tuple format as argument than the one returned get_predistorter()""" + if dpddata[0] == "poly": + coefs_am = dpddata[1] + coefs_pm = dpddata[2] + _write_poly_coef_file(coefs_am, coefs_pm, self.coef_path) + elif dpddata[0] == "lut": + scalefactor = dpddata[1] + lut = dpddata[2] + _write_lut_file(scalefactor, lut, self.coef_path) + else: + raise ValueError("Unknown predistorter '{}'".format(dpddata[0])) + self.send_receive("set memlesspoly coeffile {}".format(self.coef_path)) + + def dump(self, path=None): + """Backup current settings to a file""" + dt = datetime.datetime.now().isoformat() + if path is None: + if self.c.plot_location is not None: + path = self.c.plot_location + "/" + dt + "_adapt.pkl" + else: + raise Exception("Cannot dump Adapt without either plot_location or path set") + d = { + "txgain": self.get_txgain(), + "rxgain": self.get_rxgain(), + "digital_gain": self.get_digital_gain(), + "predistorter": self.get_predistorter() + } + with open(path, "wb") as f: + pickle.dump(d, f) + + return path + + def load(self, path): + """Restore settings from a file""" + with open(path, "rb") as f: + d = pickle.load(f) + + self.set_txgain(d["txgain"]) + self.set_digital_gain(d["digital_gain"]) + self.set_rxgain(d["rxgain"]) + self.set_predistorter(d["predistorter"]) + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger, 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/src/Dab_Util.py b/python/dpd/src/Dab_Util.py new file mode 100644 index 0000000..bc89a39 --- /dev/null +++ b/python/dpd/src/Dab_Util.py @@ -0,0 +1,246 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, utilities for working with DAB signals. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import numpy as np +import matplotlib + +matplotlib.use('agg') +import matplotlib.pyplot as plt +import src.subsample_align as sa +import src.phase_align as pa +from scipy import signal + + +def fromfile(filename, offset=0, length=None): + if length is None: + return np.memmap(filename, dtype=np.complex64, mode='r', offset=64 / 8 * offset) + else: + return np.memmap(filename, dtype=np.complex64, mode='r', offset=64 / 8 * offset, shape=length) + + +class Dab_Util: + """Collection of methods that can be applied to an array + complex IQ samples of a DAB signal + """ + + def __init__(self, config, sample_rate, plot=False): + """ + :param sample_rate: sample rate [sample/sec] to use for calculations + """ + self.c = config + self.sample_rate = sample_rate + self.dab_bandwidth = 1536000 # Bandwidth of a dab signal + self.frame_ms = 96 # Duration of a Dab frame + + self.plot = plot + + def lag(self, sig_orig, sig_rec): + """ + Find lag between two signals + Args: + sig_orig: The signal that has been sent + sig_rec: The signal that has been recored + """ + off = sig_rec.shape[0] + c = np.abs(signal.correlate(sig_orig, sig_rec)) + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + corr_path = self.c.plot_location + "/" + dt + "_tx_rx_corr.png" + plt.plot(c, label="corr") + plt.legend() + plt.savefig(corr_path) + plt.close() + + return np.argmax(c) - off + 1 + + def lag_upsampling(self, sig_orig, sig_rec, n_up): + if n_up != 1: + sig_orig_up = signal.resample(sig_orig, sig_orig.shape[0] * n_up) + sig_rec_up = signal.resample(sig_rec, sig_rec.shape[0] * n_up) + else: + sig_orig_up = sig_orig + sig_rec_up = sig_rec + l = self.lag(sig_orig_up, sig_rec_up) + l_orig = float(l) / n_up + return l_orig + + def subsample_align_upsampling(self, sig_tx, sig_rx, n_up=32): + """ + Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment + Using upsampling + """ + assert (sig_tx.shape[0] == sig_rx.shape[0]) + + if sig_tx.shape[0] % 2 == 1: + sig_tx = sig_tx[:-1] + sig_rx = sig_rx[:-1] + + sig1_up = signal.resample(sig_tx, sig_tx.shape[0] * n_up) + sig2_up = signal.resample(sig_rx, sig_rx.shape[0] * n_up) + + off_meas = self.lag_upsampling(sig2_up, sig1_up, n_up=1) + off = int(abs(off_meas)) + + if off_meas > 0: + sig1_up = sig1_up[:-off] + sig2_up = sig2_up[off:] + elif off_meas < 0: + sig1_up = sig1_up[off:] + sig2_up = sig2_up[:-off] + + sig_tx = signal.resample(sig1_up, sig1_up.shape[0] / n_up).astype(np.complex64) + sig_rx = signal.resample(sig2_up, sig2_up.shape[0] / n_up).astype(np.complex64) + return sig_tx, sig_rx + + def subsample_align(self, sig_tx, sig_rx): + """ + Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment + """ + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_sync_raw.png" + + fig, axs = plt.subplots(2) + axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") + axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") + axs[0].set_title("Raw Data") + axs[0].set_ylabel("Amplitude") + axs[0].set_xlabel("Samples") + axs[0].legend(loc=4) + + axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") + axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") + axs[1].set_title("Raw Data") + axs[1].set_ylabel("Real Part") + axs[1].set_xlabel("Samples") + axs[1].legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + off_meas = self.lag_upsampling(sig_rx, sig_tx, n_up=1) + off = int(abs(off_meas)) + + logging.debug("sig_tx_orig: {} {}, sig_rx_orig: {} {}, offset {}".format( + len(sig_tx), + sig_tx.dtype, + len(sig_rx), + sig_rx.dtype, + off_meas)) + + if off_meas > 0: + sig_tx = sig_tx[:-off] + sig_rx = sig_rx[off:] + elif off_meas < 0: + sig_tx = sig_tx[off:] + sig_rx = sig_rx[:-off] + + if off % 2 == 1: + sig_tx = sig_tx[:-1] + sig_rx = sig_rx[:-1] + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_sync_sample_aligned.png" + + fig, axs = plt.subplots(2) + axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") + axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") + axs[0].set_title("Sample Aligned Data") + axs[0].set_ylabel("Amplitude") + axs[0].set_xlabel("Samples") + axs[0].legend(loc=4) + + axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") + axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") + axs[1].set_ylabel("Real Part") + axs[1].set_xlabel("Samples") + axs[1].legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + sig_rx = sa.subsample_align(sig_rx, sig_tx) + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_sync_subsample_aligned.png" + + fig, axs = plt.subplots(2) + axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") + axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") + axs[0].set_title("Subsample Aligned") + axs[0].set_ylabel("Amplitude") + axs[0].set_xlabel("Samples") + axs[0].legend(loc=4) + + axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") + axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") + axs[1].set_ylabel("Real Part") + axs[1].set_xlabel("Samples") + axs[1].legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + sig_rx = pa.phase_align(sig_rx, sig_tx) + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_sync_phase_aligned.png" + + fig, axs = plt.subplots(2) + axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") + axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") + axs[0].set_title("Phase Aligned") + axs[0].set_ylabel("Amplitude") + axs[0].set_xlabel("Samples") + axs[0].legend(loc=4) + + axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") + axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") + axs[1].set_ylabel("Real Part") + axs[1].set_xlabel("Samples") + axs[1].legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + logging.debug( + "Sig1_cut: %d %s, Sig2_cut: %d %s, off: %d" % (len(sig_tx), sig_tx.dtype, len(sig_rx), sig_rx.dtype, off)) + return sig_tx, sig_rx + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/ExtractStatistic.py b/python/dpd/src/ExtractStatistic.py new file mode 100644 index 0000000..639513a --- /dev/null +++ b/python/dpd/src/ExtractStatistic.py @@ -0,0 +1,196 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, +# Extract statistic from received TX and RX data to use in Model +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import numpy as np +import matplotlib.pyplot as plt +import datetime +import os +import logging + + +def _check_input_extract(tx_dpd, rx_received): + # Check data type + assert tx_dpd[0].dtype == np.complex64, \ + "tx_dpd is not complex64 but {}".format(tx_dpd[0].dtype) + assert rx_received[0].dtype == np.complex64, \ + "rx_received is not complex64 but {}".format(rx_received[0].dtype) + # Check if signals have same normalization + normalization_error = np.abs(np.median(np.abs(tx_dpd)) - + np.median(np.abs(rx_received))) / ( + np.median(np.abs(tx_dpd)) + np.median(np.abs(rx_received))) + assert normalization_error < 0.01, "Non normalized signals" + + +def _phase_diff_value_per_bin(phase_diffs_values_lists): + phase_list = [] + for values in phase_diffs_values_lists: + mean = np.mean(values) if len(values) > 0 else np.nan + phase_list.append(mean) + return phase_list + + +class ExtractStatistic: + """Calculate a low variance RX value for equally spaced tx values + of a predefined range""" + + def __init__(self, c): + self.c = c + + # Number of measurements used to extract the statistic + self.n_meas = 0 + + # Boundaries for the bins + self.tx_boundaries = np.linspace(c.ES_start, c.ES_end, c.ES_n_bins + 1) + self.n_per_bin = c.ES_n_per_bin + + # List of rx values for each bin + self.rx_values_lists = [] + for i in range(c.ES_n_bins): + self.rx_values_lists.append([]) + + # List of tx values for each bin + self.tx_values_lists = [] + for i in range(c.ES_n_bins): + self.tx_values_lists.append([]) + + self.plot = c.ES_plot + + def _plot_and_log(self, tx_values, rx_values, phase_diffs_values, phase_diffs_values_lists): + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_ExtractStatistic.png" + sub_rows = 3 + sub_cols = 1 + fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) + i_sub = 0 + + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + ax.plot(tx_values, rx_values, + label="Estimated Values", + color="red") + for i, tx_value in enumerate(tx_values): + rx_values_list = self.rx_values_lists[i] + ax.scatter(np.ones(len(rx_values_list)) * tx_value, + np.abs(rx_values_list), + s=0.1, + color="black") + ax.set_title("Extracted Statistic") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("RX Amplitude") + ax.set_ylim(0, 0.8) + ax.set_xlim(0, 1.1) + ax.legend(loc=4) + + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + ax.plot(tx_values, np.rad2deg(phase_diffs_values), + label="Estimated Values", + color="red") + for i, tx_value in enumerate(tx_values): + phase_diff = phase_diffs_values_lists[i] + ax.scatter(np.ones(len(phase_diff)) * tx_value, + np.rad2deg(phase_diff), + s=0.1, + color="black") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("Phase Difference") + ax.set_ylim(-60, 60) + ax.set_xlim(0, 1.1) + ax.legend(loc=4) + + num = [] + for i, tx_value in enumerate(tx_values): + rx_values_list = self.rx_values_lists[i] + num.append(len(rx_values_list)) + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + ax.plot(num) + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("Number of Samples") + ax.set_ylim(0, self.n_per_bin * 1.2) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + def _rx_value_per_bin(self): + rx_values = [] + for values in self.rx_values_lists: + mean = np.mean(np.abs(values)) if len(values) > 0 else np.nan + rx_values.append(mean) + return rx_values + + def _tx_value_per_bin(self): + tx_values = [] + for start, end in zip(self.tx_boundaries, self.tx_boundaries[1:]): + tx_values.append(np.mean((start, end))) + return tx_values + + def _phase_diff_list_per_bin(self): + phase_values_lists = [] + for tx_list, rx_list in zip(self.tx_values_lists, self.rx_values_lists): + phase_diffs = [] + for tx, rx in zip(tx_list, rx_list): + phase_diffs.append(np.angle(rx * tx.conjugate())) + phase_values_lists.append(phase_diffs) + return phase_values_lists + + def extract(self, tx_dpd, rx): + """Extract information from a new measurement and store them + in member variables.""" + _check_input_extract(tx_dpd, rx) + self.n_meas += 1 + + tx_abs = np.abs(tx_dpd) + for i, (tx_start, tx_end) in enumerate(zip(self.tx_boundaries, self.tx_boundaries[1:])): + mask = (tx_abs > tx_start) & (tx_abs < tx_end) + n_add = max(0, self.n_per_bin - len(self.rx_values_lists[i])) + self.rx_values_lists[i] += \ + list(rx[mask][:n_add]) + self.tx_values_lists[i] += \ + list(tx_dpd[mask][:n_add]) + + rx_values = self._rx_value_per_bin() + tx_values = self._tx_value_per_bin() + + n_per_bin = np.array([len(values) for values in self.rx_values_lists]) + # Index of first not filled bin, assumes that never all bins are filled + idx_end = np.argmin(n_per_bin == self.c.ES_n_per_bin) + + phase_diffs_values_lists = self._phase_diff_list_per_bin() + phase_diffs_values = _phase_diff_value_per_bin(phase_diffs_values_lists) + + self._plot_and_log(tx_values, rx_values, phase_diffs_values, phase_diffs_values_lists) + + tx_values_crop = np.array(tx_values, dtype=np.float32)[:idx_end] + rx_values_crop = np.array(rx_values, dtype=np.float32)[:idx_end] + phase_diffs_values_crop = np.array(phase_diffs_values, dtype=np.float32)[:idx_end] + return tx_values_crop, rx_values_crop, phase_diffs_values_crop, n_per_bin + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/GlobalConfig.py b/python/dpd/src/GlobalConfig.py new file mode 100644 index 0000000..56839fc --- /dev/null +++ b/python/dpd/src/GlobalConfig.py @@ -0,0 +1,108 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, constants and global configuration +# +# Source for DAB standard: etsi_EN_300_401_v010401p p145 +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import numpy as np + +class GlobalConfig: + def __init__(self, cli_args, plot_location): + self.sample_rate = cli_args.samplerate + assert self.sample_rate == 8192000 # By now only constants for 8192000 + + self.plot_location = plot_location + + # DAB frame + # Time domain + oversample = int(self.sample_rate / 2048000) + self.T_F = oversample * 196608 # Transmission frame duration + self.T_NULL = oversample * 2656 # Null symbol duration + self.T_S = oversample * 2552 # Duration of OFDM symbols of indices l = 1, 2, 3,... L; + self.T_U = oversample * 2048 # Inverse of carrier spacing + self.T_C = oversample * 504 # Duration of cyclic prefix + + # Frequency Domain + # example: np.delete(fft[3328:4865], 768) + self.FFT_delta = 1536 # Number of carrier frequencies + self.FFT_delete = 768 + self.FFT_start = 3328 + self.FFT_end = 4865 + + # Calculate sample offset from phase rotation + # time per sample = 1 / sample_rate + # frequency per bin = 1kHz + # phase difference per sample offset = delta_t * 2 * pi * delta_freq + self.phase_offset_per_sample = 1. / self.sample_rate * 2 * np.pi * 1000 + + # Constants for ExtractStatistic + self.ES_plot = cli_args.plot + self.ES_start = 0.0 + self.ES_end = 1.0 + self.ES_n_bins = 64 # Number of bins between ES_start and ES_end + self.ES_n_per_bin = 128 # Number of measurements pre bin + + # Constants for Measure_Shoulder + self.MS_enable = False + self.MS_plot = cli_args.plot + + meas_offset = 976 # Offset from center frequency to measure shoulder [kHz] + meas_width = 100 # Size of frequency delta to measure shoulder [kHz] + shoulder_offset_edge = np.abs(meas_offset - self.FFT_delta) + self.MS_shoulder_left_start = self.FFT_start - shoulder_offset_edge - meas_width / 2 + self.MS_shoulder_left_end = self.FFT_start - shoulder_offset_edge + meas_width / 2 + self.MS_shoulder_right_start = self.FFT_end + shoulder_offset_edge - meas_width / 2 + self.MS_shoulder_right_end = self.FFT_end + shoulder_offset_edge + meas_width / 2 + self.MS_peak_start = self.FFT_start + 100 # Ignore region near edges + self.MS_peak_end = self.FFT_end - 100 + + self.MS_FFT_size = 8192 + self.MS_averaging_size = 4 * self.MS_FFT_size + self.MS_n_averaging = 40 + self.MS_n_proc = 4 + + # Constants for MER + self.MER_plot = cli_args.plot + + # Constants for Model + self.MDL_plot = cli_args.plot + + # Constants for Model_PM + # Set all phase offsets to zero for TX amplitude < MPM_tx_min + self.MPM_tx_min = 0.1 + + # Constants for TX_Agc + self.TAGC_max_txgain = 89 # USRP B200 specific + self.TAGC_tx_median_target = cli_args.target_median + self.TAGC_tx_median_max = self.TAGC_tx_median_target * 1.4 + self.TAGC_tx_median_min = self.TAGC_tx_median_target / 1.4 + + # Constants for RX_AGC + self.RAGC_min_rxgain = 25 # USRP B200 specific + self.RAGC_rx_median_target = cli_args.target_median + +# 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/src/Heuristics.py b/python/dpd/src/Heuristics.py new file mode 100644 index 0000000..21d400b --- /dev/null +++ b/python/dpd/src/Heuristics.py @@ -0,0 +1,56 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, heuristics we use to tune the parameters. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import numpy as np + + +def get_learning_rate(idx_run): + """Gradually reduce learning rate from lr_max to lr_min within + idx_max steps, then keep the learning rate at lr_min""" + idx_max = 10.0 + lr_min = 0.05 + lr_max = 0.4 + lr_delta = lr_max - lr_min + idx_run = min(idx_run, idx_max) + learning_rate = lr_max - lr_delta * idx_run / idx_max + return learning_rate + + +def get_n_meas(idx_run): + """Gradually increase number of measurements used to extract + a statistic from n_meas_min to n_meas_max within idx_max steps, + then keep number of measurements at n_meas_max""" + idx_max = 10.0 + n_meas_min = 10 + n_meas_max = 20 + n_meas_delta = n_meas_max - n_meas_min + idx_run = min(idx_run, idx_max) + learning_rate = n_meas_delta * idx_run / idx_max + n_meas_min + return int(np.round(learning_rate)) + +# 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/src/MER.py b/python/dpd/src/MER.py new file mode 100644 index 0000000..693058d --- /dev/null +++ b/python/dpd/src/MER.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, Modulation Error Rate. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import numpy as np +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + +class MER: + def __init__(self, c): + self.c = c + + self.plot = c.MER_plot + + def _calc_spectrum(self, tx): + fft = np.fft.fftshift(np.fft.fft(tx)) + return np.delete(fft[self.c.FFT_start:self.c.FFT_end], + self.c.FFT_delete) + + def _split_in_carrier(self, x, y): + if 0.5 < np.mean((np.abs(np.abs(x) - np.abs(y)) / + np.abs(np.abs(x) + np.abs(y)))): + # Constellation points are axis aligned on the Im/Re plane + x1 = x[(y < x) & (y > -x)] + y1 = y[(y < x) & (y > -x)] + + x2 = x[(y > x) & (y > -x)] + y2 = y[(y > x) & (y > -x)] + + x3 = x[(y > x) & (y < -x)] + y3 = y[(y > x) & (y < -x)] + + x4 = x[(y < x) & (y < -x)] + y4 = y[(y < x) & (y < -x)] + else: + # Constellation points are on the diagonal or Im/Re plane + x1 = x[(+x > 0) & (+y > 0)] + y1 = y[(+x > 0) & (+y > 0)] + + x2 = x[(-x > 0) & (+y > 0)] + y2 = y[(-x > 0) & (+y > 0)] + + x3 = x[(-x > 0) & (-y > 0)] + y3 = y[(-x > 0) & (-y > 0)] + + x4 = x[(+x > 0) & (-y > 0)] + y4 = y[(+x > 0) & (-y > 0)] + return (x1, y1), (x2, y2), (x3, y3), (x4, y4) + + def _calc_mer_for_isolated_constellation_point(self, x, y): + """Calculate MER for one constellation point""" + x_mean = np.mean(x) + y_mean = np.mean(y) + + U_RMS = np.sqrt(x_mean ** 2 + y_mean ** 2) + U_ERR = np.mean(np.sqrt((x - x_mean) ** 2 + (y - y_mean) ** 2)) + MER = 20 * np.log10(U_ERR / U_RMS) + + return x_mean, y_mean, U_RMS, U_ERR, MER + + def calc_mer(self, tx, debug_name=""): + """Calculate MER for input signal from a symbol aligned signal.""" + assert tx.shape[0] == self.c.T_U, "Wrong input length" + + spectrum = self._calc_spectrum(tx) + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_MER" + debug_name + ".png" + else: + fig_path = None + + MERs = [] + for i, (x, y) in enumerate(self._split_in_carrier( + np.real(spectrum), + np.imag(spectrum))): + x_mean, y_mean, U_RMS, U_ERR, MER =\ + self._calc_mer_for_isolated_constellation_point(x, y) + MERs.append(MER) + + tau = np.linspace(0, 2 * np.pi, num=100) + x_err = U_ERR * np.sin(tau) + x_mean + y_err = U_ERR * np.cos(tau) + y_mean + + if self.plot: + ax = plt.subplot(221 + i) + ax.scatter(x, y, s=0.2, color='black') + ax.plot(x_mean, y_mean, 'p', color='red') + ax.plot(x_err, y_err, linewidth=2, color='blue') + ax.text(0.1, 0.1, "MER {:.1f}dB".format(MER), transform=ax.transAxes) + ax.set_xlabel("Real") + ax.set_ylabel("Imag") + ylim = ax.get_ylim() + ax.set_ylim(ylim[0] - (ylim[1] - ylim[0]) * 0.1, ylim[1]) + + if fig_path is not None: + plt.tight_layout() + plt.savefig(fig_path) + plt.show() + plt.close() + + MER_res = 20 * np.log10(np.mean([10 ** (MER / 20) for MER in MERs])) + return MER_res + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Measure.py b/python/dpd/src/Measure.py new file mode 100644 index 0000000..6d8007d --- /dev/null +++ b/python/dpd/src/Measure.py @@ -0,0 +1,136 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, Measure signal using ODR-DabMod's +# DPD Server. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import socket +import struct +import numpy as np +import src.Dab_Util as DU +import os +import logging + +class Measure: + """Collect Measurement from DabMod""" + def __init__(self, config, samplerate, port, num_samples_to_request): + logging.info("Instantiate Measure object") + self.c = config + self.samplerate = samplerate + self.sizeof_sample = 8 # complex floats + self.port = port + self.num_samples_to_request = num_samples_to_request + + def _recv_exact(self, sock, num_bytes): + """Receive an exact number of bytes from a socket. This is + a wrapper around sock.recv() that can return less than the number + of requested bytes. + + Args: + sock (socket): Socket to receive data from. + num_bytes (int): Number of bytes that will be returned. + """ + 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 receive_tcp(self): + s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + s.settimeout(4) + s.connect(('localhost', self.port)) + + logging.debug("Send version") + s.sendall(b"\x01") + + logging.debug("Send request for {} samples".format(self.num_samples_to_request)) + s.sendall(struct.pack("=I", self.num_samples_to_request)) + + logging.debug("Wait for TX metadata") + num_samps, tx_second, tx_pps = struct.unpack("=III", self._recv_exact(s, 12)) + tx_ts = tx_second + tx_pps / 16384000.0 + + if num_samps > 0: + logging.debug("Receiving {} TX samples".format(num_samps)) + txframe_bytes = self._recv_exact(s, num_samps * self.sizeof_sample) + txframe = np.fromstring(txframe_bytes, dtype=np.complex64) + else: + txframe = np.array([], dtype=np.complex64) + + logging.debug("Wait for RX metadata") + rx_second, rx_pps = struct.unpack("=II", self._recv_exact(s, 8)) + rx_ts = rx_second + rx_pps / 16384000.0 + + if num_samps > 0: + logging.debug("Receiving {} RX samples".format(num_samps)) + rxframe_bytes = self._recv_exact(s, num_samps * self.sizeof_sample) + rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64) + else: + rxframe = np.array([], dtype=np.complex64) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG: + logging.debug('txframe: min {}, max {}, median {}'.format( + np.min(np.abs(txframe)), + np.max(np.abs(txframe)), + np.median(np.abs(txframe)))) + + logging.debug('rxframe: min {}, max {}, median {}'.format( + np.min(np.abs(rxframe)), + np.max(np.abs(rxframe)), + np.median(np.abs(rxframe)))) + + logging.debug("Disconnecting") + s.close() + + return txframe, tx_ts, rxframe, rx_ts + + + def get_samples(self): + """Connect to ODR-DabMod, retrieve TX and RX samples, load + into numpy arrays, and return a tuple + (txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median) + """ + + txframe, tx_ts, rxframe, rx_ts = self.receive_tcp() + + # Normalize received signal with sent signal + rx_median = np.median(np.abs(rxframe)) + rxframe = rxframe / rx_median * np.median(np.abs(txframe)) + + du = DU.Dab_Util(self.c, self.samplerate) + txframe_aligned, rxframe_aligned = du.subsample_align(txframe, rxframe) + + logging.info( + "Measurement done, tx %d %s, rx %d %s, tx aligned %d %s, rx aligned %d %s" + % (len(txframe), txframe.dtype, len(rxframe), rxframe.dtype, + len(txframe_aligned), txframe_aligned.dtype, len(rxframe_aligned), rxframe_aligned.dtype) ) + + return txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Measure_Shoulders.py b/python/dpd/src/Measure_Shoulders.py new file mode 100644 index 0000000..fd90050 --- /dev/null +++ b/python/dpd/src/Measure_Shoulders.py @@ -0,0 +1,158 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, calculate peak to shoulder difference. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import multiprocessing +import numpy as np +import matplotlib.pyplot as plt + + +def plt_next_axis(sub_rows, sub_cols, i_sub): + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + return i_sub, ax + + +def plt_annotate(ax, x, y, title=None, legend_loc=None): + ax.set_xlabel(x) + ax.set_ylabel(y) + if title is not None: + ax.set_title(title) + if legend_loc is not None: + ax.legend(loc=legend_loc) + + +def calc_fft_db(signal, offset, c): + fft = np.fft.fftshift(np.fft.fft(signal[offset:offset + c.MS_FFT_size])) + fft_db = 20 * np.log10(np.abs(fft)) + return fft_db + + +def _calc_peak(fft, c): + assert fft.shape == (c.MS_FFT_size,), fft.shape + idxs = (c.MS_peak_start, c.MS_peak_end) + peak = np.mean(fft[idxs[0]:idxs[1]]) + return peak, idxs + + +def _calc_shoulder_hight(fft_db, c): + assert fft_db.shape == (c.MS_FFT_size,), fft_db.shape + idxs_left = (c.MS_shoulder_left_start, c.MS_shoulder_left_end) + idxs_right = (c.MS_shoulder_right_start, c.MS_shoulder_right_end) + + shoulder_left = np.mean(fft_db[idxs_left[0]:idxs_left[1]]) + shoulder_right = np.mean(fft_db[idxs_right[0]:idxs_right[1]]) + + shoulder = np.mean((shoulder_left, shoulder_right)) + return shoulder, (idxs_left, idxs_right) + + +def calc_shoulder(fft, c): + peak = _calc_peak(fft, c)[0] + shoulder = _calc_shoulder_hight(fft, c)[0] + assert (peak >= shoulder), (peak, shoulder) + return peak, shoulder + + +def shoulder_from_sig_offset(arg): + signal, offset, c = arg + fft_db = calc_fft_db(signal, offset, c) + peak, shoulder = calc_shoulder(fft_db, c) + return peak - shoulder, peak, shoulder + + +class Measure_Shoulders: + """Calculate difference between the DAB signal and the shoulder hight in the + power spectrum""" + + def __init__(self, c): + self.c = c + self.plot = c.MS_plot + + def _plot(self, signal): + if self.c.plot_location is None: + return + + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_sync_subsample_aligned.png" + + fft = calc_fft_db(signal, 100, self.c) + peak, idxs_peak = _calc_peak(fft, self.c) + shoulder, idxs_sh = _calc_shoulder_hight(fft, self.c) + + sub_rows = 1 + sub_cols = 1 + fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) + i_sub = 0 + + i_sub, ax = plt_next_axis(sub_rows, sub_cols, i_sub) + ax.scatter(np.arange(fft.shape[0]), fft, s=0.1, + label="FFT", + color="red") + ax.plot(idxs_peak, (peak, peak)) + ax.plot(idxs_sh[0], (shoulder, shoulder), color='blue') + ax.plot(idxs_sh[1], (shoulder, shoulder), color='blue') + plt_annotate(ax, "Frequency", "Magnitude [dB]", None, 4) + + ax.text(100, -17, str(calc_shoulder(fft, self.c))) + + ax.set_ylim(-20, 30) + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + def average_shoulders(self, signal, n_avg=None): + if not self.c.MS_enable: + logging.info("Shoulder Measurement disabled via Const.py") + return None + + assert signal.shape[0] > 4 * self.c.MS_FFT_size + if n_avg is None: + n_avg = self.c.MS_averaging_size + + off_min = 0 + off_max = signal.shape[0] - self.c.MS_FFT_size + offsets = np.linspace(off_min, off_max, num=n_avg, dtype=int) + + args = zip( + [signal, ] * offsets.shape[0], + offsets, + [self.c, ] * offsets.shape[0] + ) + + pool = multiprocessing.Pool(self.c.MS_n_proc) + res = pool.map(shoulder_from_sig_offset, args) + shoulders_diff, shoulders, peaks = zip(*res) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + self._plot(signal) + + return np.mean(shoulders_diff), np.mean(shoulders), np.mean(peaks) + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Model.py b/python/dpd/src/Model.py new file mode 100644 index 0000000..b2c303f --- /dev/null +++ b/python/dpd/src/Model.py @@ -0,0 +1,32 @@ +# -*- coding: utf-8 -*- +from src.Model_Poly import Poly +from src.Model_Lut import Lut + +def select_model_from_dpddata(dpddata): + if dpddata[0] == 'lut': + return Lut + elif dpddata[0] == 'poly': + return Poly + +# 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/src/Model_AM.py b/python/dpd/src/Model_AM.py new file mode 100644 index 0000000..75b226f --- /dev/null +++ b/python/dpd/src/Model_AM.py @@ -0,0 +1,122 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, model implementation for Amplitude and not Phase +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import numpy as np +import matplotlib.pyplot as plt + + +def is_npfloat32(array): + assert isinstance(array, np.ndarray), type(array) + assert array.dtype == np.float32, array.dtype + assert array.flags.contiguous + assert not any(np.isnan(array)) + + +def check_input_get_next_coefs(tx_dpd, rx_received): + is_npfloat32(tx_dpd) + is_npfloat32(rx_received) + + +def poly(sig): + return np.array([sig ** i for i in range(1, 6)]).T + + +def fit_poly(tx_abs, rx_abs): + return np.linalg.lstsq(poly(rx_abs), tx_abs, rcond=None)[0] + + +def calc_line(coefs, min_amp, max_amp): + rx_range = np.linspace(min_amp, max_amp) + tx_est = np.sum(poly(rx_range) * coefs, axis=1) + return tx_est, rx_range + + +class Model_AM: + """Calculates new coefficients using the measurement and the previous + coefficients""" + + def __init__(self, + c, + learning_rate_am=1, + plot=False): + self.c = c + + self.learning_rate_am = learning_rate_am + self.plot = plot + + def _plot(self, tx_dpd, rx_received, coefs_am, coefs_am_new): + if self.plot and self.c.plot_location is not None: + tx_range, rx_est = calc_line(coefs_am, 0, 0.6) + tx_range_new, rx_est_new = calc_line(coefs_am_new, 0, 0.6) + + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_Model_AM.png" + sub_rows = 1 + sub_cols = 1 + fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) + i_sub = 0 + + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + ax.plot(tx_range, rx_est, + label="Estimated TX", + alpha=0.3, + color="gray") + ax.plot(tx_range_new, rx_est_new, + label="New Estimated TX", + color="red") + ax.scatter(tx_dpd, rx_received, + label="Binned Data", + color="blue", + s=1) + ax.set_title("Model_AM") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("RX Amplitude") + ax.set_xlim(-0.5, 1.5) + ax.legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + def get_next_coefs(self, tx_dpd, rx_received, coefs_am): + """Calculate the next AM/AM coefficients using the extracted + statistic of TX and RX amplitude""" + check_input_get_next_coefs(tx_dpd, rx_received) + + coefs_am_new = fit_poly(tx_dpd, rx_received) + coefs_am_new = coefs_am + \ + self.learning_rate_am * (coefs_am_new - coefs_am) + + self._plot(tx_dpd, rx_received, coefs_am, coefs_am_new) + + return coefs_am_new + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Model_Lut.py b/python/dpd/src/Model_Lut.py new file mode 100644 index 0000000..e70fdb0 --- /dev/null +++ b/python/dpd/src/Model_Lut.py @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, model implementation using polynomial +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import os +import logging +import numpy as np + +class Lut: + """Implements a model that calculates lookup table coefficients""" + + def __init__(self, + c, + learning_rate=1., + plot=False): + """ + + :rtype: + """ + logging.debug("Initialising LUT Model") + self.c = c + self.learning_rate = learning_rate + self.plot = plot + self.reset_coefs() + + def reset_coefs(self): + self.scalefactor = 0xFFFFFFFF # max uint32_t value + self.lut = np.ones(32, dtype=np.complex64) + + def train(self, tx_abs, rx_abs, phase_diff): + pass + + def get_dpd_data(self): + return "lut", self.scalefactor, self.lut + +# 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/src/Model_PM.py b/python/dpd/src/Model_PM.py new file mode 100644 index 0000000..7b80bf3 --- /dev/null +++ b/python/dpd/src/Model_PM.py @@ -0,0 +1,124 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, model implementation for Amplitude and not Phase +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import numpy as np +import matplotlib.pyplot as plt + + +def is_npfloat32(array): + assert isinstance(array, np.ndarray), type(array) + assert array.dtype == np.float32, array.dtype + assert array.flags.contiguous + assert not any(np.isnan(array)) + + +def check_input_get_next_coefs(tx_dpd, phase_diff): + is_npfloat32(tx_dpd) + is_npfloat32(phase_diff) + + +class Model_PM: + """Calculates new coefficients using the measurement and the previous + coefficients""" + + def __init__(self, + c, + learning_rate_pm=1, + plot=False): + self.c = c + + self.learning_rate_pm = learning_rate_pm + self.plot = plot + + def _plot(self, tx_dpd, phase_diff, coefs_pm, coefs_pm_new): + if self.plot and self.c.plot_location is not None: + tx_range, phase_diff_est = self.calc_line(coefs_pm, 0, 0.6) + tx_range_new, phase_diff_est_new = self.calc_line(coefs_pm_new, 0, 0.6) + + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_Model_PM.png" + sub_rows = 1 + sub_cols = 1 + fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) + i_sub = 0 + + i_sub += 1 + ax = plt.subplot(sub_rows, sub_cols, i_sub) + ax.plot(tx_range, phase_diff_est, + label="Estimated Phase Diff", + alpha=0.3, + color="gray") + ax.plot(tx_range_new, phase_diff_est_new, + label="New Estimated Phase Diff", + color="red") + ax.scatter(tx_dpd, phase_diff, + label="Binned Data", + color="blue", + s=1) + ax.set_title("Model_PM") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("Phase DIff") + ax.legend(loc=4) + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + def _discard_small_values(self, tx_dpd, phase_diff): + """ Assumes that the phase for small tx amplitudes is zero""" + mask = tx_dpd < self.c.MPM_tx_min + phase_diff[mask] = 0 + return tx_dpd, phase_diff + + def poly(self, sig): + return np.array([sig ** i for i in range(0, 5)]).T + + def fit_poly(self, tx_abs, phase_diff): + return np.linalg.lstsq(self.poly(tx_abs), phase_diff, rcond=None)[0] + + def calc_line(self, coefs, min_amp, max_amp): + tx_range = np.linspace(min_amp, max_amp) + phase_diff = np.sum(self.poly(tx_range) * coefs, axis=1) + return tx_range, phase_diff + + def get_next_coefs(self, tx_dpd, phase_diff, coefs_pm): + """Calculate the next AM/PM coefficients using the extracted + statistic of TX amplitude and phase difference""" + tx_dpd, phase_diff = self._discard_small_values(tx_dpd, phase_diff) + check_input_get_next_coefs(tx_dpd, phase_diff) + + coefs_pm_new = self.fit_poly(tx_dpd, phase_diff) + + coefs_pm_new = coefs_pm + self.learning_rate_pm * (coefs_pm_new - coefs_pm) + self._plot(tx_dpd, phase_diff, coefs_pm, coefs_pm_new) + + return coefs_pm_new + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Model_Poly.py b/python/dpd/src/Model_Poly.py new file mode 100644 index 0000000..cdfd319 --- /dev/null +++ b/python/dpd/src/Model_Poly.py @@ -0,0 +1,101 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, model implementation using polynomial +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import os +import logging +import numpy as np + +import src.Model_AM as Model_AM +import src.Model_PM as Model_PM + + +def assert_np_float32(x): + assert isinstance(x, np.ndarray) + assert x.dtype == np.float32 + assert x.flags.contiguous + + +def _check_input_get_next_coefs(tx_abs, rx_abs, phase_diff): + assert_np_float32(tx_abs) + assert_np_float32(rx_abs) + assert_np_float32(phase_diff) + + assert tx_abs.shape == rx_abs.shape, \ + "tx_abs.shape {}, rx_abs.shape {}".format( + tx_abs.shape, rx_abs.shape) + assert tx_abs.shape == phase_diff.shape, \ + "tx_abs.shape {}, phase_diff.shape {}".format( + tx_abs.shape, phase_diff.shape) + + +class Poly: + """Calculates new coefficients using the measurement and the previous + coefficients""" + + def __init__(self, + c, + learning_rate_am=1.0, + learning_rate_pm=1.0): + self.c = c + self.plot = c.MDL_plot + + self.learning_rate_am = learning_rate_am + self.learning_rate_pm = learning_rate_pm + + self.reset_coefs() + + self.model_am = Model_AM.Model_AM(c, plot=self.plot) + self.model_pm = Model_PM.Model_PM(c, plot=self.plot) + + def reset_coefs(self): + self.coefs_am = np.zeros(5, dtype=np.float32) + self.coefs_am[0] = 1 + self.coefs_pm = np.zeros(5, dtype=np.float32) + + def train(self, tx_abs, rx_abs, phase_diff, lr=None): + """ + :type tx_abs: np.ndarray + :type rx_abs: np.ndarray + :type phase_diff: np.ndarray + :type lr: float + """ + _check_input_get_next_coefs(tx_abs, rx_abs, phase_diff) + + if not lr is None: + self.model_am.learning_rate_am = lr + self.model_pm.learning_rate_pm = lr + + coefs_am_new = self.model_am.get_next_coefs(tx_abs, rx_abs, self.coefs_am) + coefs_pm_new = self.model_pm.get_next_coefs(tx_abs, phase_diff, self.coefs_pm) + + self.coefs_am = self.coefs_am + (coefs_am_new - self.coefs_am) * self.learning_rate_am + self.coefs_pm = self.coefs_pm + (coefs_pm_new - self.coefs_pm) * self.learning_rate_pm + + def get_dpd_data(self): + return "poly", self.coefs_am, self.coefs_pm + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/RX_Agc.py b/python/dpd/src/RX_Agc.py new file mode 100644 index 0000000..f778dee --- /dev/null +++ b/python/dpd/src/RX_Agc.py @@ -0,0 +1,166 @@ +# -*- coding: utf-8 -*- +# +# Automatic Gain Control +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import time +import numpy as np +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + +import src.Adapt as Adapt +import src.Measure as Measure + +class Agc: + """The goal of the automatic gain control is to set the + RX gain to a value at which all received amplitudes can + be detected. This means that the maximum possible amplitude + should be quantized at the highest possible digital value. + + A problem we have to face, is that the estimation of the + maximum amplitude by applying the max() function is very + unstable. This is due to the maximum’s rareness. Therefore + we estimate a far more robust value, such as the median, + and then approximate the maximum amplitude from it. + + Given this, we tune the RX gain in such a way, that the + received signal fulfills our desired property, of having + all amplitudes properly quantized.""" + + def __init__(self, measure, adapt, c): + assert isinstance(measure, Measure.Measure) + assert isinstance(adapt, Adapt.Adapt) + self.measure = measure + self.adapt = adapt + self.min_rxgain = c.RAGC_min_rxgain + self.rxgain = self.min_rxgain + self.peak_to_median = 1./c.RAGC_rx_median_target + + def run(self): + self.adapt.set_rxgain(self.rxgain) + + for i in range(2): + # Measure + txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median= \ + self.measure.get_samples() + + # Estimate Maximum + rx_peak = self.peak_to_median * rx_median + correction_factor = 20*np.log10(1/rx_peak) + self.rxgain = self.rxgain + correction_factor + + assert self.min_rxgain <= self.rxgain, ("Desired RX Gain is {} which is smaller than the minimum of {}".format( + self.rxgain, self.min_rxgain)) + + logging.info("RX Median {:1.4f}, estimated peak {:1.4f}, correction factor {:1.4f}, new RX gain {:1.4f}".format( + rx_median, rx_peak, correction_factor, self.rxgain + )) + + self.adapt.set_rxgain(self.rxgain) + time.sleep(0.5) + + def plot_estimates(self): + """Plots the estimate of for Max, Median, Mean for different + number of samples.""" + if self.c.plot_location is None: + return + + self.adapt.set_rxgain(self.min_rxgain) + time.sleep(1) + + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_agc.png" + fig, axs = plt.subplots(2, 2, figsize=(3*6,1*6)) + axs = axs.ravel() + + for j in range(5): + txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median =\ + self.measure.get_samples() + + rxframe_aligned_abs = np.abs(rxframe_aligned) + + x = np.arange(100, rxframe_aligned_abs.shape[0], dtype=int) + rx_max_until = [] + rx_median_until = [] + rx_mean_until = [] + for i in x: + rx_max_until.append(np.max(rxframe_aligned_abs[:i])) + rx_median_until.append(np.median(rxframe_aligned_abs[:i])) + rx_mean_until.append(np.mean(rxframe_aligned_abs[:i])) + + axs[0].plot(x, + rx_max_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.8,0.7)), + linestyle="-", linewidth=0.25) + axs[0].set_xlabel("Samples") + axs[0].set_ylabel("Amplitude") + axs[0].set_title("Estimation for Maximum RX Amplitude") + axs[0].legend() + + axs[1].plot(x, + rx_median_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), + linestyle="-", linewidth=0.25) + axs[1].set_xlabel("Samples") + axs[1].set_ylabel("Amplitude") + axs[1].set_title("Estimation for Median RX Amplitude") + axs[1].legend() + ylim_1 = axs[1].get_ylim() + + axs[2].plot(x, + rx_mean_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), + linestyle="-", linewidth=0.25) + axs[2].set_xlabel("Samples") + axs[2].set_ylabel("Amplitude") + axs[2].set_title("Estimation for Mean RX Amplitude") + ylim_2 = axs[2].get_ylim() + axs[2].legend() + + axs[1].set_ylim(min(ylim_1[0], ylim_2[0]), + max(ylim_1[1], ylim_2[1])) + + fig.tight_layout() + fig.savefig(fig_path) + + axs[3].hist(rxframe_aligned_abs, bins=60) + axs[3].set_xlabel("Amplitude") + axs[3].set_ylabel("Frequency") + axs[3].set_title("Histogram of Amplitudes") + axs[3].legend() + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/Symbol_align.py b/python/dpd/src/Symbol_align.py new file mode 100644 index 0000000..2a17a65 --- /dev/null +++ b/python/dpd/src/Symbol_align.py @@ -0,0 +1,193 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, Modulation Error Rate. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import numpy as np +import scipy +import matplotlib + +matplotlib.use('agg') +import matplotlib.pyplot as plt + + +def _remove_outliers(x, stds=5): + deviation_from_mean = np.abs(x - np.mean(x)) + inlier_idxs = deviation_from_mean < stds * np.std(x) + x = x[inlier_idxs] + return x + + +def _calc_delta_angle(fft): + # Introduce invariance against carrier + angles = np.angle(fft) % (np.pi / 2.) + + # Calculate Angle difference and compensate jumps + deltas_angle = np.diff(angles) + deltas_angle[deltas_angle > np.pi / 4.] = \ + deltas_angle[deltas_angle > np.pi / 4.] - np.pi / 2. + deltas_angle[-deltas_angle > np.pi / 4.] = \ + deltas_angle[-deltas_angle > np.pi / 4.] + np.pi / 2. + deltas_angle = _remove_outliers(deltas_angle) + + delta_angle = np.mean(deltas_angle) + + return delta_angle + + +class Symbol_align: + """ + Find the phase offset to the start of the DAB symbols in an + unaligned dab signal. + """ + + def __init__(self, c, plot=False): + self.c = c + self.plot = plot + pass + + def _calc_offset_to_first_symbol_without_prefix(self, tx): + tx_orig = tx[0:-self.c.T_U] + tx_cut_prefix = tx[self.c.T_U:] + + tx_product = np.abs(tx_orig - tx_cut_prefix) + tx_product_avg = np.correlate( + tx_product, + np.ones(self.c.T_C), + mode='valid') + tx_product_avg_min_filt = \ + scipy.ndimage.filters.minimum_filter1d( + tx_product_avg, + int(1.5 * self.c.T_S) + ) + peaks = np.ravel(np.where(tx_product_avg == tx_product_avg_min_filt)) + + offset = peaks[np.argmin([tx_product_avg[peak] for peak in peaks])] + + if self.plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_Symbol_align.png" + + fig = plt.figure(figsize=(9, 9)) + + ax = fig.add_subplot(4, 1, 1) + ax.plot(tx_product) + ylim = ax.get_ylim() + for peak in peaks: + ax.plot((peak, peak), (ylim[0], ylim[1])) + if peak == offset: + ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) + else: + ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) + ax.set_xlabel("Sample") + ax.set_ylabel("Conj. Product") + ax.set_title("Difference with shifted self") + + ax = fig.add_subplot(4, 1, 2) + ax.plot(tx_product_avg) + ylim = ax.get_ylim() + for peak in peaks: + ax.plot((peak, peak), (ylim[0], ylim[1])) + if peak == offset: + ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) + else: + ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) + ax.set_xlabel("Sample") + ax.set_ylabel("Conj. Product") + ax.set_title("Moving Average") + + ax = fig.add_subplot(4, 1, 3) + ax.plot(tx_product_avg_min_filt) + ylim = ax.get_ylim() + for peak in peaks: + ax.plot((peak, peak), (ylim[0], ylim[1])) + if peak == offset: + ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) + else: + ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) + ax.set_xlabel("Sample") + ax.set_ylabel("Conj. Product") + ax.set_title("Min Filter") + + ax = fig.add_subplot(4, 1, 4) + tx_product_crop = tx_product[peaks[0] - 50:peaks[0] + 50] + x = range(tx_product.shape[0])[peaks[0] - 50:peaks[0] + 50] + ax.plot(x, tx_product_crop) + ylim = ax.get_ylim() + ax.plot((peaks[0], peaks[0]), (ylim[0], ylim[1])) + ax.set_xlabel("Sample") + ax.set_ylabel("Conj. Product") + ax.set_title("Difference with shifted self") + + fig.tight_layout() + fig.savefig(fig_path) + plt.close(fig) + + # "offset" measures where the shifted signal matches the + # original signal. Therefore we have to subtract the size + # of the shift to find the offset of the symbol start. + return (offset + self.c.T_C) % self.c.T_S + + def _delta_angle_to_samples(self, angle): + return - angle / self.c.phase_offset_per_sample + + def _calc_sample_offset(self, sig): + assert sig.shape[0] == self.c.T_U, \ + "Input length is not a Symbol without cyclic prefix" + + fft = np.fft.fftshift(np.fft.fft(sig)) + fft_crop = np.delete(fft[self.c.FFT_start:self.c.FFT_end], self.c.FFT_delete) + delta_angle = _calc_delta_angle(fft_crop) + delta_sample = self._delta_angle_to_samples(delta_angle) + delta_sample_int = np.round(delta_sample).astype(int) + error = np.abs(delta_sample_int - delta_sample) + if error > 0.1: + raise RuntimeError("Could not calculate " \ + "sample offset. Error {}".format(error)) + return delta_sample_int + + def calc_offset(self, tx): + """Calculate the offset the first symbol""" + off_sym = self._calc_offset_to_first_symbol_without_prefix( + tx) + off_sam = self._calc_sample_offset( + tx[off_sym:off_sym + self.c.T_U]) + off = (off_sym + off_sam) % self.c.T_S + + assert self._calc_sample_offset(tx[off:off + self.c.T_U]) == 0, \ + "Failed to calculate offset" + return off + + def crop_symbol_without_cyclic_prefix(self, tx): + off = self.calc_offset(tx) + return tx[ + off: + off + self.c.T_U + ] + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/TX_Agc.py b/python/dpd/src/TX_Agc.py new file mode 100644 index 0000000..309193d --- /dev/null +++ b/python/dpd/src/TX_Agc.py @@ -0,0 +1,131 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, Automatic Gain Control. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import time +import numpy as np +import matplotlib + +matplotlib.use('agg') +import matplotlib.pyplot as plt + +import src.Adapt as Adapt + + +# TODO fix for float tx_gain +class TX_Agc: + def __init__(self, + adapt, + c): + """ + In order to avoid digital clipping, this class increases the + TX gain and reduces the digital gain. Digital clipping happens + when the digital analog converter receives values greater than + it's maximal output. This class solves that problem by adapting + the TX gain in a way that the peaks of the TX signal are in a + specified range. The TX gain is adapted accordingly. The TX peaks + are approximated by estimating it based on the signal median. + + :param adapt: Instance of Adapt Class to update + txgain and coefficients + :param max_txgain: limit for TX gain + :param tx_median_threshold_max: if the median of TX is larger + than this value, then the digital gain is reduced + :param tx_median_threshold_min: if the median of TX is smaller + than this value, then the digital gain is increased + :param tx_median_target: The digital gain is reduced in a way that + the median TX value is expected to be lower than this value. + """ + + assert isinstance(adapt, Adapt.Adapt) + self.adapt = adapt + self.max_txgain = c.TAGC_max_txgain + self.txgain = self.max_txgain + + self.tx_median_threshold_tolerate_max = c.TAGC_tx_median_max + self.tx_median_threshold_tolerate_min = c.TAGC_tx_median_min + self.tx_median_target = c.TAGC_tx_median_target + + def _calc_new_tx_gain(self, tx_median): + delta_db = 20 * np.log10(self.tx_median_target / tx_median) + new_txgain = self.adapt.get_txgain() - delta_db + assert new_txgain < self.max_txgain, \ + "TX_Agc failed. New TX gain of {} is too large.".format( + new_txgain + ) + return new_txgain, delta_db + + def _calc_digital_gain(self, delta_db): + digital_gain_factor = 10 ** (delta_db / 20.) + digital_gain = self.adapt.get_digital_gain() * digital_gain_factor + return digital_gain, digital_gain_factor + + def _set_tx_gain(self, new_txgain): + self.adapt.set_txgain(new_txgain) + txgain = self.adapt.get_txgain() + return txgain + + def _have_to_adapt(self, tx_median): + too_large = tx_median > self.tx_median_threshold_tolerate_max + too_small = tx_median < self.tx_median_threshold_tolerate_min + return too_large or too_small + + def adapt_if_necessary(self, tx): + tx_median = np.median(np.abs(tx)) + + if self._have_to_adapt(tx_median): + # Calculate new values + new_txgain, delta_db = self._calc_new_tx_gain(tx_median) + digital_gain, digital_gain_factor = \ + self._calc_digital_gain(delta_db) + + # Set new values. + # Avoid temorary increase of output power with correct order + if digital_gain_factor < 1: + self.adapt.set_digital_gain(digital_gain) + time.sleep(0.5) + txgain = self._set_tx_gain(new_txgain) + time.sleep(1) + else: + txgain = self._set_tx_gain(new_txgain) + time.sleep(1) + self.adapt.set_digital_gain(digital_gain) + time.sleep(0.5) + + logging.info( + "digital_gain = {}, txgain_new = {}, " \ + "delta_db = {}, tx_median {}, " \ + "digital_gain_factor = {}". + format(digital_gain, txgain, delta_db, + tx_median, digital_gain_factor)) + + return True + return False + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/__init__.py b/python/dpd/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/dpd/src/phase_align.py b/python/dpd/src/phase_align.py new file mode 100644 index 0000000..8654333 --- /dev/null +++ b/python/dpd/src/phase_align.py @@ -0,0 +1,98 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, phase-align a signal against a reference. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file +import datetime +import os +import logging +import numpy as np +import matplotlib.pyplot as plt + + +def phase_align(sig, ref_sig, plot=False): + """Do phase alignment for sig relative to the reference signal + ref_sig. + + Returns the aligned signal""" + + angle_diff = (np.angle(sig) - np.angle(ref_sig)) % (2. * np.pi) + + real_diffs = np.cos(angle_diff) + imag_diffs = np.sin(angle_diff) + + if plot and self.c.plot_location is not None: + dt = datetime.datetime.now().isoformat() + fig_path = self.c.plot_location + "/" + dt + "_phase_align.png" + + plt.subplot(511) + plt.hist(angle_diff, bins=60, label="Angle Diff") + plt.xlabel("Angle") + plt.ylabel("Count") + plt.legend(loc=4) + + plt.subplot(512) + plt.hist(real_diffs, bins=60, label="Real Diff") + plt.xlabel("Real Part") + plt.ylabel("Count") + plt.legend(loc=4) + + plt.subplot(513) + plt.hist(imag_diffs, bins=60, label="Imaginary Diff") + plt.xlabel("Imaginary Part") + plt.ylabel("Count") + plt.legend(loc=4) + + plt.subplot(514) + plt.plot(np.angle(ref_sig[:128]), label="ref_sig") + plt.plot(np.angle(sig[:128]), label="sig") + plt.xlabel("Angle") + plt.ylabel("Sample") + plt.legend(loc=4) + + real_diff = np.median(real_diffs) + imag_diff = np.median(imag_diffs) + + angle = np.angle(real_diff + 1j * imag_diff) + + logging.debug( + "Compensating phase by {} rad, {} degree. real median {}, imag median {}".format( + angle, angle*180./np.pi, real_diff, imag_diff + )) + sig = sig * np.exp(1j * -angle) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and plot: + plt.subplot(515) + plt.plot(np.angle(ref_sig[:128]), label="ref_sig") + plt.plot(np.angle(sig[:128]), label="sig") + plt.xlabel("Angle") + plt.ylabel("Sample") + plt.legend(loc=4) + plt.tight_layout() + plt.savefig(fig_path) + plt.close() + + return sig + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/src/subsample_align.py b/python/dpd/src/subsample_align.py new file mode 100755 index 0000000..20ae56b --- /dev/null +++ b/python/dpd/src/subsample_align.py @@ -0,0 +1,111 @@ +# -*- coding: utf-8 -*- +# +# DPD Computation Engine, utility to do subsample alignment. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file +import datetime +import logging +import os +import numpy as np +from scipy import optimize +import matplotlib.pyplot as plt + +def gen_omega(length): + if (length % 2) == 1: + raise ValueError("Needs an even length array.") + + halflength = int(length / 2) + factor = 2.0 * np.pi / length + + omega = np.zeros(length, dtype=np.float) + for i in range(halflength): + omega[i] = factor * i + + for i in range(halflength, length): + omega[i] = factor * (i - length) + + return omega + + +def subsample_align(sig, ref_sig, plot_location=None): + """Do subsample alignment for sig relative to the reference signal + ref_sig. The delay between the two must be less than sample + + Returns the aligned signal""" + + n = len(sig) + if (n % 2) == 1: + raise ValueError("Needs an even length signal.") + halflen = int(n / 2) + + fft_sig = np.fft.fft(sig) + + omega = gen_omega(n) + + def correlate_for_delay(tau): + # A subsample offset between two signals corresponds, in the frequency + # domain, to a linearly increasing phase shift, whose slope + # corresponds to the delay. + # + # Here, we build this phase shift in rotate_vec, and multiply it with + # our signal. + + rotate_vec = np.exp(1j * tau * omega) + # zero-frequency is rotate_vec[0], so rotate_vec[N/2] is the + # bin corresponding to the [-1, 1, -1, 1, ...] time signal, which + # is both the maximum positive and negative frequency. + # I don't remember why we handle it differently. + rotate_vec[halflen] = np.cos(np.pi * tau) + + corr_sig = np.fft.ifft(rotate_vec * fft_sig) + + return -np.abs(np.sum(np.conj(corr_sig) * ref_sig)) + + optim_result = optimize.minimize_scalar(correlate_for_delay, bounds=(-1, 1), method='bounded', + options={'disp': True}) + + if optim_result.success: + best_tau = optim_result.x + + if plot_location is not None: + corr = np.vectorize(correlate_for_delay) + ixs = np.linspace(-1, 1, 100) + taus = corr(ixs) + + dt = datetime.datetime.now().isoformat() + tau_path = (plot_location + "/" + dt + "_tau.png") + plt.plot(ixs, taus) + plt.title("Subsample correlation, minimum is best: {}".format(best_tau)) + plt.savefig(tau_path) + plt.close() + + # Prepare rotate_vec = fft_sig with rotated phase + rotate_vec = np.exp(1j * best_tau * omega) + rotate_vec[halflen] = np.cos(np.pi * best_tau) + return np.fft.ifft(rotate_vec * fft_sig).astype(np.complex64) + else: + # print("Could not optimize: " + optim_result.message) + return np.zeros(0, dtype=np.complex64) + +# The MIT License (MIT) +# +# Copyright (c) 2017 Andreas Steger +# +# 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/store_received.py b/python/dpd/store_received.py new file mode 100755 index 0000000..19b735e --- /dev/null +++ b/python/dpd/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. -- cgit v1.2.3