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 --- dpd/README.md | 264 --------------------- dpd/apply_adapt_dumps.py | 75 ------ dpd/dpd.ini | 60 ----- dpd/img/setup_diagram.svg | 395 -------------------------------- dpd/img/setup_photo.svg | 184 --------------- dpd/img/shoulder_measurement_after.png | Bin 368619 -> 0 bytes dpd/img/shoulder_measurement_before.png | Bin 381174 -> 0 bytes dpd/iq_file_server.py | 120 ---------- dpd/lut.coef | 64 ------ dpd/main.py | 336 --------------------------- dpd/poly.coef | 12 - dpd/show_spectrum.py | 276 ---------------------- dpd/src/Adapt.py | 286 ----------------------- dpd/src/Dab_Util.py | 246 -------------------- dpd/src/ExtractStatistic.py | 196 ---------------- dpd/src/GlobalConfig.py | 108 --------- dpd/src/Heuristics.py | 56 ----- dpd/src/MER.py | 132 ----------- dpd/src/Measure.py | 136 ----------- dpd/src/Measure_Shoulders.py | 158 ------------- dpd/src/Model.py | 32 --- dpd/src/Model_AM.py | 122 ---------- dpd/src/Model_Lut.py | 60 ----- dpd/src/Model_PM.py | 124 ---------- dpd/src/Model_Poly.py | 101 -------- dpd/src/RX_Agc.py | 166 -------------- dpd/src/Symbol_align.py | 193 ---------------- dpd/src/TX_Agc.py | 131 ----------- dpd/src/__init__.py | 0 dpd/src/phase_align.py | 98 -------- dpd/src/subsample_align.py | 111 --------- dpd/store_received.py | 85 ------- 32 files changed, 4327 deletions(-) delete mode 100644 dpd/README.md delete mode 100755 dpd/apply_adapt_dumps.py delete mode 100644 dpd/dpd.ini delete mode 100644 dpd/img/setup_diagram.svg delete mode 100644 dpd/img/setup_photo.svg delete mode 100644 dpd/img/shoulder_measurement_after.png delete mode 100644 dpd/img/shoulder_measurement_before.png delete mode 100755 dpd/iq_file_server.py delete mode 100644 dpd/lut.coef delete mode 100755 dpd/main.py delete mode 100644 dpd/poly.coef delete mode 100755 dpd/show_spectrum.py delete mode 100644 dpd/src/Adapt.py delete mode 100644 dpd/src/Dab_Util.py delete mode 100644 dpd/src/ExtractStatistic.py delete mode 100644 dpd/src/GlobalConfig.py delete mode 100644 dpd/src/Heuristics.py delete mode 100644 dpd/src/MER.py delete mode 100644 dpd/src/Measure.py delete mode 100644 dpd/src/Measure_Shoulders.py delete mode 100644 dpd/src/Model.py delete mode 100644 dpd/src/Model_AM.py delete mode 100644 dpd/src/Model_Lut.py delete mode 100644 dpd/src/Model_PM.py delete mode 100644 dpd/src/Model_Poly.py delete mode 100644 dpd/src/RX_Agc.py delete mode 100644 dpd/src/Symbol_align.py delete mode 100644 dpd/src/TX_Agc.py delete mode 100644 dpd/src/__init__.py delete mode 100644 dpd/src/phase_align.py delete mode 100755 dpd/src/subsample_align.py delete mode 100755 dpd/store_received.py (limited to 'dpd') diff --git a/dpd/README.md b/dpd/README.md deleted file mode 100644 index 307a2f5..0000000 --- a/dpd/README.md +++ /dev/null @@ -1,264 +0,0 @@ -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/dpd/apply_adapt_dumps.py b/dpd/apply_adapt_dumps.py deleted file mode 100755 index 20bc013..0000000 --- a/dpd/apply_adapt_dumps.py +++ /dev/null @@ -1,75 +0,0 @@ -#!/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/dpd/dpd.ini b/dpd/dpd.ini deleted file mode 100644 index 31d6140..0000000 --- a/dpd/dpd.ini +++ /dev/null @@ -1,60 +0,0 @@ -[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/dpd/img/setup_diagram.svg b/dpd/img/setup_diagram.svg deleted file mode 100644 index 423ee1e..0000000 --- a/dpd/img/setup_diagram.svg +++ /dev/null @@ -1,395 +0,0 @@ - - - - - - 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/dpd/img/setup_photo.svg b/dpd/img/setup_photo.svg deleted file mode 100644 index a401fda..0000000 --- a/dpd/img/setup_photo.svg +++ /dev/null @@ -1,184 +0,0 @@ - - - - - - - - image/svg+xml - - - - - - - - - ① - ② - - - ③ - ④ - ⑤⑥ - ⑦ - ⑧ - diff --git a/dpd/img/shoulder_measurement_after.png b/dpd/img/shoulder_measurement_after.png deleted file mode 100644 index 2631256..0000000 Binary files a/dpd/img/shoulder_measurement_after.png and /dev/null differ diff --git a/dpd/img/shoulder_measurement_before.png b/dpd/img/shoulder_measurement_before.png deleted file mode 100644 index 3581a72..0000000 Binary files a/dpd/img/shoulder_measurement_before.png and /dev/null differ diff --git a/dpd/iq_file_server.py b/dpd/iq_file_server.py deleted file mode 100755 index 7a4e570..0000000 --- a/dpd/iq_file_server.py +++ /dev/null @@ -1,120 +0,0 @@ -#!/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/dpd/lut.coef b/dpd/lut.coef deleted file mode 100644 index a198d56..0000000 --- a/dpd/lut.coef +++ /dev/null @@ -1,64 +0,0 @@ -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/dpd/main.py b/dpd/main.py deleted file mode 100755 index 10a56fc..0000000 --- a/dpd/main.py +++ /dev/null @@ -1,336 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# -# DPD Computation Engine 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. -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/dpd/poly.coef b/dpd/poly.coef deleted file mode 100644 index 248d316..0000000 --- a/dpd/poly.coef +++ /dev/null @@ -1,12 +0,0 @@ -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/dpd/show_spectrum.py b/dpd/show_spectrum.py deleted file mode 100755 index f23dba2..0000000 --- a/dpd/show_spectrum.py +++ /dev/null @@ -1,276 +0,0 @@ -#!/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/dpd/src/Adapt.py b/dpd/src/Adapt.py deleted file mode 100644 index a57602f..0000000 --- a/dpd/src/Adapt.py +++ /dev/null @@ -1,286 +0,0 @@ -# -*- 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/dpd/src/Dab_Util.py b/dpd/src/Dab_Util.py deleted file mode 100644 index bc89a39..0000000 --- a/dpd/src/Dab_Util.py +++ /dev/null @@ -1,246 +0,0 @@ -# -*- 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/dpd/src/ExtractStatistic.py b/dpd/src/ExtractStatistic.py deleted file mode 100644 index 639513a..0000000 --- a/dpd/src/ExtractStatistic.py +++ /dev/null @@ -1,196 +0,0 @@ -# -*- 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/dpd/src/GlobalConfig.py b/dpd/src/GlobalConfig.py deleted file mode 100644 index 56839fc..0000000 --- a/dpd/src/GlobalConfig.py +++ /dev/null @@ -1,108 +0,0 @@ -# -*- 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/dpd/src/Heuristics.py b/dpd/src/Heuristics.py deleted file mode 100644 index 21d400b..0000000 --- a/dpd/src/Heuristics.py +++ /dev/null @@ -1,56 +0,0 @@ -# -*- 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/dpd/src/MER.py b/dpd/src/MER.py deleted file mode 100644 index 693058d..0000000 --- a/dpd/src/MER.py +++ /dev/null @@ -1,132 +0,0 @@ -# -*- 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/dpd/src/Measure.py b/dpd/src/Measure.py deleted file mode 100644 index 6d8007d..0000000 --- a/dpd/src/Measure.py +++ /dev/null @@ -1,136 +0,0 @@ -# -*- 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/dpd/src/Measure_Shoulders.py b/dpd/src/Measure_Shoulders.py deleted file mode 100644 index fd90050..0000000 --- a/dpd/src/Measure_Shoulders.py +++ /dev/null @@ -1,158 +0,0 @@ -# -*- 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/dpd/src/Model.py b/dpd/src/Model.py deleted file mode 100644 index b2c303f..0000000 --- a/dpd/src/Model.py +++ /dev/null @@ -1,32 +0,0 @@ -# -*- 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/dpd/src/Model_AM.py b/dpd/src/Model_AM.py deleted file mode 100644 index 75b226f..0000000 --- a/dpd/src/Model_AM.py +++ /dev/null @@ -1,122 +0,0 @@ -# -*- 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/dpd/src/Model_Lut.py b/dpd/src/Model_Lut.py deleted file mode 100644 index e70fdb0..0000000 --- a/dpd/src/Model_Lut.py +++ /dev/null @@ -1,60 +0,0 @@ -# -*- 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/dpd/src/Model_PM.py b/dpd/src/Model_PM.py deleted file mode 100644 index 7b80bf3..0000000 --- a/dpd/src/Model_PM.py +++ /dev/null @@ -1,124 +0,0 @@ -# -*- 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/dpd/src/Model_Poly.py b/dpd/src/Model_Poly.py deleted file mode 100644 index cdfd319..0000000 --- a/dpd/src/Model_Poly.py +++ /dev/null @@ -1,101 +0,0 @@ -# -*- 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/dpd/src/RX_Agc.py b/dpd/src/RX_Agc.py deleted file mode 100644 index f778dee..0000000 --- a/dpd/src/RX_Agc.py +++ /dev/null @@ -1,166 +0,0 @@ -# -*- 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/dpd/src/Symbol_align.py b/dpd/src/Symbol_align.py deleted file mode 100644 index 2a17a65..0000000 --- a/dpd/src/Symbol_align.py +++ /dev/null @@ -1,193 +0,0 @@ -# -*- 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/dpd/src/TX_Agc.py b/dpd/src/TX_Agc.py deleted file mode 100644 index 309193d..0000000 --- a/dpd/src/TX_Agc.py +++ /dev/null @@ -1,131 +0,0 @@ -# -*- 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/dpd/src/__init__.py b/dpd/src/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/dpd/src/phase_align.py b/dpd/src/phase_align.py deleted file mode 100644 index 8654333..0000000 --- a/dpd/src/phase_align.py +++ /dev/null @@ -1,98 +0,0 @@ -# -*- 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/dpd/src/subsample_align.py b/dpd/src/subsample_align.py deleted file mode 100755 index 20ae56b..0000000 --- a/dpd/src/subsample_align.py +++ /dev/null @@ -1,111 +0,0 @@ -# -*- 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/dpd/store_received.py b/dpd/store_received.py deleted file mode 100755 index 19b735e..0000000 --- a/dpd/store_received.py +++ /dev/null @@ -1,85 +0,0 @@ -#!/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