diff options
Diffstat (limited to 'dpd')
| -rw-r--r-- | dpd/README.md | 48 | ||||
| -rw-r--r-- | dpd/dpd.ini | 25 | ||||
| -rwxr-xr-x | dpd/iq_file_server.py | 120 | ||||
| -rwxr-xr-x | dpd/main.py | 176 | ||||
| -rw-r--r-- | dpd/poly.coef | 11 | ||||
| -rwxr-xr-x | dpd/show_spectrum.py | 132 | ||||
| -rw-r--r-- | dpd/src/Adapt.py | 190 | ||||
| -rw-r--r-- | dpd/src/Agc.py | 165 | ||||
| -rw-r--r-- | dpd/src/Dab_Util.py | 244 | ||||
| -rw-r--r-- | dpd/src/MER.py | 134 | ||||
| -rw-r--r-- | dpd/src/Measure.py | 145 | ||||
| -rw-r--r-- | dpd/src/Model.py | 355 | ||||
| -rw-r--r-- | dpd/src/Symbol_align.py | 191 | ||||
| -rw-r--r-- | dpd/src/TX_Agc.py | 109 | ||||
| -rw-r--r-- | dpd/src/Test_data.py | 136 | ||||
| -rw-r--r-- | dpd/src/__init__.py | 0 | ||||
| -rw-r--r-- | dpd/src/const.py | 38 | ||||
| -rw-r--r-- | dpd/src/phase_align.py | 102 | ||||
| -rwxr-xr-x | dpd/src/subsample_align.py | 112 | ||||
| -rw-r--r-- | dpd/src/test_dab_Util.py | 92 | ||||
| -rw-r--r-- | dpd/src/test_measure.py | 62 |
21 files changed, 2557 insertions, 30 deletions
diff --git a/dpd/README.md b/dpd/README.md index ec7cec2..b5a6b81 100644 --- a/dpd/README.md +++ b/dpd/README.md @@ -1,20 +1,56 @@ -Digital Predistortion for ODR-DabMod -==================================== +Digital Predistortion Calculation Engine for ODR-DabMod +======================================================= -This folder contains work in progress for digital predistortion. It requires: +This folder contains work in progress for digital predistortion. + +Concept +------- + +ODR-DabMod makes outgoing TX samples and feedback RX samples available for an external tool. This +external tool can request a buffer of samples for analysis, can calculate coefficients for the +polynomial predistorter in ODR-DabMod and load the new coefficients using the remote control. + +The *dpd/main.py* script is the entry point for the *DPD Calculation 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* +- Finally, *Adapt.py* loads them into ODR-DabMod. + +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, at an appropriate power level for the B200. - Usually this is done with a directional coupler. -- ODR-DabMod with enabled dpd_port, and with a samplerate of 8192000 samples per second. + 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 for this example. - 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. + +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 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 2xN + 1 lines if it contains N +coefficients. + TODO ---- -Implement a PA model that updates the predistorter. +Implement a PA model. +Implement cases for different oversampling for FFT bin choice. +Fix loads of missing and buggy aspects of the implementation. diff --git a/dpd/dpd.ini b/dpd/dpd.ini index 5e809e5..1bc51de 100644 --- a/dpd/dpd.ini +++ b/dpd/dpd.ini @@ -6,35 +6,48 @@ zmqctrlendpoint=tcp://127.0.0.1:9400 [log] syslog=0 -filelog=0 -filename=/dev/stderr +filelog=1 +filename=/tmp/dabmod.log [input] transport=tcp source=localhost:9200 [modulator] -digital_gain=0.8 +gainmode=var rate=8192000 [firfilter] -enabled=0 +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=75 +txgain=15 channel=13C refclk_source=internal pps_source=none behaviour_refclk_lock_lost=ignore max_gps_holdover_time=600 dpd_port=50055 -rxgain=0 +rxgain=15 [delaymanagement] ; Use synchronous=1 so that the USRP time is set. This works diff --git a/dpd/iq_file_server.py b/dpd/iq_file_server.py new file mode 100755 index 0000000..7a4e570 --- /dev/null +++ b/dpd/iq_file_server.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# This example server simulates the ODR-DabMod's +# DPD server, taking samples from an IQ file +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import sys +import socket +import struct +import argparse +import numpy as np +from datetime import datetime + +SIZEOF_SAMPLE = 8 # complex floats +# Constants for TM 1 +NbSymbols = 76 +NbCarriers = 1536 +Spacing = 2048 +NullSize = 2656 +SymSize = 2552 +FicSizeOut = 288 +FrameSize = NullSize + NbSymbols*SymSize + +def main(): + parser = argparse.ArgumentParser(description="Simulate ODR-DabMod DPD server") + parser.add_argument('--port', default='50055', + help='port to listen on (default: 50055)', + required=False) + parser.add_argument('--file', help='I/Q File to read from (complex floats)', + required=True) + parser.add_argument('--samplerate', default='8192000', help='Sample rate', + required=False) + + cli_args = parser.parse_args() + + serve_file(cli_args) + +def recv_exact(sock, num_bytes): + bufs = [] + while num_bytes > 0: + b = sock.recv(num_bytes) + if len(b) == 0: + break + num_bytes -= len(b) + bufs.append(b) + return b''.join(bufs) + +def serve_file(options): + oversampling = int(int(options.samplerate) / 2048000) + consumesamples = 8*FrameSize * oversampling + iq_data = np.fromfile(options.file, count=consumesamples, dtype=np.complex64) + + print("Loaded {} samples of IQ data".format(len(iq_data))) + + s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + s.bind(('localhost', int(options.port))) + s.listen() + + try: + while True: + sock, addr_info = s.accept() + print("Got a connection from {}".format(addr_info)) + + ver = recv_exact(sock, 1) + (num_samps,) = struct.unpack("=I", recv_exact(sock, 4)) + num_bytes = num_samps * SIZEOF_SAMPLE + + if num_bytes > len(iq_data): + print("Truncating length to {} samples".format(len(iq_data))) + num_samps = len(iq_data) + num_bytes = num_samps * 4 + + tx_sec = datetime.now().timestamp() + tx_pps = int(16384000 * (tx_sec - int(tx_sec))) + tx_second = int(tx_sec) + + # TX metadata and data + sock.sendall(struct.pack("=III", num_samps, tx_second, tx_pps)) + sock.sendall(iq_data[-num_samps:].tobytes()) + + # RX metadata and data + rx_second = tx_second + 1 + rx_pps = tx_pps + sock.sendall(struct.pack("=III", num_samps, rx_second, rx_pps)) + sock.sendall(iq_data[-num_samps:].tobytes()) + + print("Sent {} samples".format(num_samps)) + + sock.close() + finally: + s.close() + raise + +main() + + +# The MIT License (MIT) +# +# Copyright (c) 2017 Matthias P. Braendli +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. diff --git a/dpd/main.py b/dpd/main.py new file mode 100755 index 0000000..5e67c90 --- /dev/null +++ b/dpd/main.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# DPD Calculation 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 datetime +import os +import time + +import matplotlib +matplotlib.use('GTKAgg') + +import logging + +dt = datetime.datetime.now().isoformat() +logging_path = "/tmp/dpd_{}".format(dt).replace(".", "_").replace(":", "-") +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) + +import numpy as np +import traceback +import src.Measure as Measure +import src.Model as Model +import src.Adapt as Adapt +import src.Agc as Agc +import src.TX_Agc as TX_Agc +import src.Symbol_align +import src.const +import src.MER +import argparse + +parser = argparse.ArgumentParser( + description="DPD Computation Engine for ODR-DabMod") +parser.add_argument('--port', default='50055', + help='port of DPD server to connect to (default: 50055)', + required=False) +parser.add_argument('--rc-port', default='9400', + help='port of ODR-DabMod ZMQ Remote Control to connect to (default: 9400)', + required=False) +parser.add_argument('--samplerate', default='8192000', + 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=71, + help='TX Gain', + required=False, + type=int) +parser.add_argument('--rxgain', default=30, + help='TX Gain', + required=False, + type=int) +parser.add_argument('--digital_gain', default=1, + help='Digital Gain', + required=False, + type=float) +parser.add_argument('--samps', default='81920', + help='Number of samples to request from ODR-DabMod', + required=False) +parser.add_argument('-i', '--iterations', default='1', + help='Number of iterations to run', + required=False) +parser.add_argument('-l', '--load-poly', + help='Load existing polynomial', + action="store_true") + +cli_args = parser.parse_args() + +port = int(cli_args.port) +port_rc = int(cli_args.rc_port) +coef_path = cli_args.coefs +digital_gain = cli_args.digital_gain +txgain = cli_args.txgain +rxgain = cli_args.rxgain +num_req = int(cli_args.samps) +samplerate = int(cli_args.samplerate) +num_iter = int(cli_args.iterations) + +SA = src.Symbol_align.Symbol_align(samplerate) +MER = src.MER.MER(samplerate) +c = src.const.const(samplerate) + +meas = Measure.Measure(samplerate, port, num_req) + +adapt = Adapt.Adapt(port_rc, coef_path) +coefs_am, coefs_pm = adapt.get_coefs() +if cli_args.load_poly: + model = Model.Model(c, SA, MER, coefs_am, coefs_pm, plot=True) +else: + model = Model.Model(c, SA, MER, [1.0, 0, 0, 0, 0], [0, 0, 0, 0, 0], plot=True) +adapt.set_coefs(model.coefs_am, model.coefs_pm) +adapt.set_digital_gain(digital_gain) +adapt.set_txgain(txgain) +adapt.set_rxgain(rxgain) + +tx_gain = adapt.get_txgain() +rx_gain = adapt.get_rxgain() +digital_gain = adapt.get_digital_gain() +dpd_coefs_am, dpd_coefs_pm = adapt.get_coefs() +logging.info( + "TX gain {}, RX gain {}, dpd_coefs_am {}," + " dpd_coefs_pm {}, digital_gain {}".format( + tx_gain, rx_gain, dpd_coefs_am, dpd_coefs_pm, digital_gain + ) +) + +tx_agc = TX_Agc.TX_Agc(adapt) + +# Automatic Gain Control +agc = Agc.Agc(meas, adapt) +agc.run() + +for i in range(num_iter): + try: + txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples() + logging.debug("tx_ts {}, rx_ts {}".format(tx_ts, rx_ts)) + assert tx_ts - rx_ts < 1e-5, "Time stamps do not match." + + if tx_agc.adapt_if_necessary(txframe_aligned): + continue + + coefs_am, coefs_pm = model.get_next_coefs(txframe_aligned, rxframe_aligned) + adapt.set_coefs(coefs_am, coefs_pm) + + off = SA.calc_offset(txframe_aligned) + tx_mer = MER.calc_mer(txframe_aligned[off:off + c.T_U]) + rx_mer = MER.calc_mer(rxframe_aligned[off:off + c.T_U], debug=True) + logging.info("MER with lag in it. {}: TX {}, RX {}". + format(i, tx_mer, rx_mer)) + except Exception as e: + logging.warning("Iteration {} failed.".format(i)) + logging.warning(traceback.format_exc()) + +# 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/poly.coef b/dpd/poly.coef new file mode 100644 index 0000000..f65e369 --- /dev/null +++ b/dpd/poly.coef @@ -0,0 +1,11 @@ +5 +1 +0 +0 +0 +0 +0 +0 +0 +0 +0 diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index e92c1d0..f23dba2 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -1,16 +1,15 @@ #!/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. +# 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. +# 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 @@ -21,9 +20,18 @@ 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', @@ -31,18 +39,19 @@ def main(): 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() - port = int(cli_args.port) - num_samps_to_request = int(cli_args.samps) - - if cli_args.animated: - plot_spectrum_animated(port, num_samps_to_request) + if cli_args.constellation: + plot_constellation_once(cli_args) + elif cli_args.animated: + plot_spectrum_animated(cli_args) else: - plot_spectrum_once(port, num_samps_to_request) + plot_spectrum_once(cli_args) def recv_exact(sock, num_bytes): bufs = [] @@ -82,6 +91,7 @@ def get_samples(port, num_samps_to_request): 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 @@ -98,7 +108,7 @@ def get_samples(port, num_samps_to_request): return (tx_ts, txframe, rx_ts, rxframe) -def get_spectrum(port, num_samps_to_request): +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 @@ -107,6 +117,10 @@ def get_spectrum(port, num_samps_to_request): 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)) @@ -116,12 +130,90 @@ def get_spectrum(port, num_samps_to_request): 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() -sampling_rate = 8192000 fft_size = 4096 -freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate)) -def plot_spectrum_once(port, num_samps_to_request): +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() @@ -134,7 +226,11 @@ def plot_spectrum_once(port, num_samps_to_request): ax2.plot(freqs, rx_power, 'b') pp.show() -def plot_spectrum_animated(port, num_samps_to_request): +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") diff --git a/dpd/src/Adapt.py b/dpd/src/Adapt.py new file mode 100644 index 0000000..b4042d6 --- /dev/null +++ b/dpd/src/Adapt.py @@ -0,0 +1,190 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation Engine: updates ODR-DabMod's predistortion block. +# +# 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 + +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, port, coef_path): + logging.debug("Instantiate Adapt object") + 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) + sock.connect("tcp://%s:%d" % (self.host, self.port)) + + sock.send(b"ping") + data = [el.decode() for el in sock.recv_multipart()] + + if data != ['ok']: + raise RuntimeError( + "Could not ping server at %s %d: %s" % + (self.host, self.port, data)) + + 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 uhd txgain" or "set uhd 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 uhd 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 uhd 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 uhd 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 uhd 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 _read_coef_file(self, path): + """Load the coefficients from the file in the format given in the README, + return ([AM coef], [PM coef])""" + coefs_am_out = [] + coefs_pm_out = [] + f = open(path, 'r') + lines = f.readlines() + n_coefs = int(lines[0]) + coefs = [float(l) for l in lines[1:]] + 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(path, n_coefs, coefs)) + i += 1 + f.close() + return (coefs_am_out, coefs_pm_out) + + def get_coefs(self): + return self._read_coef_file(self.coef_path) + + def _write_coef_file(self, coefs_am, coefs_pm, path): + assert(len(coefs_am) == len(coefs_pm)) + + f = open(path, 'w') + f.write("{}\n".format(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 set_coefs(self, coefs_am, coefs_pm): + self._write_coef_file(coefs_am, coefs_pm, self.coef_path) + self.send_receive("set memlesspoly coeffile {}".format(self.coef_path)) + +# 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/Agc.py b/dpd/src/Agc.py new file mode 100644 index 0000000..978b607 --- /dev/null +++ b/dpd/src/Agc.py @@ -0,0 +1,165 @@ +# -*- 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 +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +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, min_rxgain=25, peak_to_median=20): + assert isinstance(measure, Measure.Measure) + assert isinstance(adapt, Adapt.Adapt) + self.measure = measure + self.adapt = adapt + self.min_rxgain = min_rxgain + self.rxgain = self.min_rxgain + self.peak_to_median = peak_to_median + + 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.""" + self.adapt.set_rxgain(self.min_rxgain) + time.sleep(1) + + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_agc.svg" + 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) + fig.clf() + + +# 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/Dab_Util.py b/dpd/src/Dab_Util.py new file mode 100644 index 0000000..e3dbfe3 --- /dev/null +++ b/dpd/src/Dab_Util.py @@ -0,0 +1,244 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation 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 +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +import numpy as np +import scipy +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 + +class Dab_Util: + """Collection of methods that can be applied to an array + complex IQ samples of a DAB signal + """ + def __init__(self, sample_rate, plot=False): + """ + :param sample_rate: sample rate [sample/sec] to use for calculations + """ + 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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + corr_path = (logging_path + "/" + dt + "_tx_rx_corr.svg") + plt.plot(c, label="corr") + plt.legend() + plt.savefig(corr_path) + plt.clf() + + 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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_sync_raw.svg" + + 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) + fig.clf() + + 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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_sync_sample_aligned.svg" + + 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) + fig.clf() + + sig_rx = sa.subsample_align(sig_rx, sig_tx) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_sync_subsample_aligned.svg" + + 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) + fig.clf() + + sig_rx = pa.phase_align(sig_rx, sig_tx) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_sync_phase_aligned.svg" + + 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) + fig.clf() + + 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 + + def fromfile(self, 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) + + +# 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/MER.py b/dpd/src/MER.py new file mode 100644 index 0000000..00fcc23 --- /dev/null +++ b/dpd/src/MER.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*- +# +# 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 +try: + logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) +except: + logging_path = "/tmp/" + +import src.const +import numpy as np +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + +class MER: + def __init__(self, sample_rate): + self.c = src.const.const(sample_rate) + + 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=False): + assert tx.shape[0] == self.c.T_U,\ + "Wrong input length" + + spectrum = self._calc_spectrum(tx) + + if debug: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_MER.svg" + + 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 debug: + 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 debug: + plt.tight_layout() + plt.savefig(fig_path) + plt.show() + plt.clf() + + 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 new file mode 100644 index 0000000..2450b8a --- /dev/null +++ b/dpd/src/Measure.py @@ -0,0 +1,145 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation 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 sys +import socket +import struct +import numpy as np +import logging +import src.Dab_Util as DU +import os +import datetime + +import logging +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +class Measure: + """Collect Measurement from DabMod""" + def __init__(self, samplerate, port, num_samples_to_request): + logging.info("Instantiate Measure object") + 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.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 %f, max %f, median %f" % + (np.min(np.abs(txframe)), + np.max(np.abs(txframe)), + np.median(np.abs(txframe)))) + + logging.debug("rxframe: min %f, max %f, median %f" % + (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 + (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 + """ + + txframe, tx_ts, rxframe, rx_ts = self.receive_tcp() + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG: + dt = datetime.datetime.now().isoformat() + txframe.tofile(logging_path + "/txframe_" + dt + ".iq") + + # 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.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/Model.py b/dpd/src/Model.py new file mode 100644 index 0000000..827027a --- /dev/null +++ b/dpd/src/Model.py @@ -0,0 +1,355 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation Engine, model implementation. +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging + +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +import numpy as np +import matplotlib.pyplot as plt +from sklearn import linear_model + +class Model: + """Calculates new coefficients using the measurement and the old + coefficients""" + + def __init__(self, + c, + SA, + MER, + coefs_am, + coefs_pm, + learning_rate_am=1., + learning_rate_pm=1., + plot=False): + self.c = c + self.SA = SA + self.MER = MER + + self.learning_rate_am = learning_rate_am + self.learning_rate_pm = learning_rate_pm + + self.coefs_am = coefs_am + self.coefs_am_history = [coefs_am, ] + self.mses_am = [] + self.errs_am = [] + + self.tx_mers = [] + self.rx_mers = [] + + self.coefs_pm = coefs_pm + self.coefs_pm_history = [coefs_pm, ] + self.errs_pm = [] + + self.plot = plot + + def sample_uniformly(self, tx_dpd, rx_received, n_bins=5): + """This function returns tx and rx samples in a way + that the tx amplitudes have an approximate uniform + distribution with respect to the tx_dpd amplitudes""" + mask = np.logical_and((np.abs(tx_dpd) > 0.01), (np.abs(rx_received) > 0.01)) + tx_dpd = tx_dpd[mask] + rx_received = rx_received[mask] + + txframe_aligned_abs = np.abs(tx_dpd) + ccdf_min = 0 + ccdf_max = np.max(txframe_aligned_abs) + tx_hist, ccdf_edges = np.histogram(txframe_aligned_abs, + bins=n_bins, + range=(ccdf_min, ccdf_max)) + n_choise = np.min(tx_hist) + tx_choice = np.zeros(n_choise * n_bins, dtype=np.complex64) + rx_choice = np.zeros(n_choise * n_bins, dtype=np.complex64) + + for idx, bin in enumerate(tx_hist): + indices = np.where((txframe_aligned_abs >= ccdf_edges[idx]) & + (txframe_aligned_abs <= ccdf_edges[idx + 1]))[0] + indices_choise = np.random.choice(indices, + n_choise, + replace=False) + rx_choice[idx * n_choise:(idx + 1) * n_choise] = \ + rx_received[indices_choise] + tx_choice[idx * n_choise:(idx + 1) * n_choise] = \ + tx_dpd[indices_choise] + + assert isinstance(rx_choice[0], np.complex64), \ + "rx_choice is not complex64 but {}".format(rx_choice[0].dtype) + assert isinstance(tx_choice[0], np.complex64), \ + "tx_choice is not complex64 but {}".format(tx_choice[0].dtype) + + return tx_choice, rx_choice + + def dpd_amplitude(self, sig, coefs=None): + if coefs is None: + coefs = self.coefs_am + assert isinstance(sig[0], np.complex64), "Sig is not complex64 but {}".format(sig[0].dtype) + sig_abs = np.abs(sig) + A_sig = np.vstack([np.ones(sig_abs.shape), + sig_abs ** 1, + sig_abs ** 2, + sig_abs ** 3, + sig_abs ** 4, + ]).T + sig_dpd = sig * np.sum(A_sig * coefs, axis=1) + return sig_dpd, A_sig + + def dpd_phase(self, sig, coefs=None): + if coefs is None: + coefs = self.coefs_pm + assert isinstance(sig[0], np.complex64), "Sig is not complex64 but {}".format(sig[0].dtype) + sig_abs = np.abs(sig) + A_phase = np.vstack([np.ones(sig_abs.shape), + sig_abs ** 1, + sig_abs ** 2, + sig_abs ** 3, + sig_abs ** 4, + ]).T + phase_diff_est = np.sum(A_phase * coefs, axis=1) + return phase_diff_est, A_phase + + def _next_am_coefficent(self, tx_choice, rx_choice): + """Calculate new coefficients for AM/AM correction""" + rx_dpd, rx_A = self.dpd_amplitude(rx_choice) + rx_dpd = rx_dpd * ( + np.median(np.abs(tx_choice)) / + np.median(np.abs(rx_dpd))) + err = np.abs(rx_dpd) - np.abs(tx_choice) + mse = np.mean(np.abs((rx_dpd - tx_choice) ** 2)) + self.mses_am.append(mse) + self.errs_am.append(np.mean(err**2)) + + reg = linear_model.Ridge(alpha=0.00001) + reg.fit(rx_A, err) + a_delta = reg.coef_ + new_coefs_am = self.coefs_am - self.learning_rate_am * a_delta + new_coefs_am = new_coefs_am * (self.coefs_am[0] / new_coefs_am[0]) + return new_coefs_am + + def _next_pm_coefficent(self, tx_choice, rx_choice): + """Calculate new coefficients for AM/PM correction + Assuming deviations smaller than pi/2""" + phase_diff_choice = np.angle( + (rx_choice * tx_choice.conjugate()) / + (np.abs(rx_choice) * np.abs(tx_choice)) + ) + plt.hist(phase_diff_choice) + plt.savefig('/tmp/hist_' + str(np.random.randint(0,1000)) + '.svg') + plt.clf() + phase_diff_est, phase_A = self.dpd_phase(rx_choice) + err_phase = phase_diff_est - phase_diff_choice + self.errs_pm.append(np.mean(np.abs(err_phase ** 2))) + + reg = linear_model.Ridge(alpha=0.00001) + reg.fit(phase_A, err_phase) + p_delta = reg.coef_ + new_coefs_pm = self.coefs_pm - self.learning_rate_pm * p_delta + + return new_coefs_pm, phase_diff_choice + + def get_next_coefs(self, 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" + + tx_choice, rx_choice = self.sample_uniformly(tx_dpd, rx_received) + new_coefs_am = self._next_am_coefficent(tx_choice, rx_choice) + new_coefs_pm, phase_diff_choice = self._next_pm_coefficent(tx_choice, rx_choice) + + logging.debug('txframe: min {:.2f}, max {:.2f}, ' \ + 'median {:.2f}; rxframe: min {:.2f}, max {:.2f}, ' \ + 'median {:.2f}; new coefs_am {};' \ + 'new_coefs_pm {}'.format( + np.min(np.abs(tx_dpd)), + np.max(np.abs(tx_dpd)), + np.median(np.abs(tx_dpd)), + np.min(np.abs(rx_choice)), + np.max(np.abs(rx_choice)), + np.median(np.abs(rx_choice)), + new_coefs_am, + new_coefs_pm)) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + off = self.SA.calc_offset(tx_dpd) + tx_mer = self.MER.calc_mer(tx_dpd[off:off + self.c.T_U]) + rx_mer = self.MER.calc_mer(rx_received[off:off + self.c.T_U], debug=True) + self.tx_mers.append(tx_mer) + self.rx_mers.append(rx_mer) + + if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_Model.svg" + + fig = plt.figure(figsize=(2 * 6, 2 * 6)) + + i_sub = 1 + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + ax.plot(np.abs(tx_dpd[:128]), + label="TX sent", + linestyle=":") + ax.plot(np.abs(rx_received[:128]), + label="RX received", + color="red") + ax.set_title("Synchronized Signals of Iteration {}" + .format(len(self.coefs_am_history))) + ax.set_xlabel("Samples") + ax.set_ylabel("Amplitude") + ax.text(0, 0, "TX (max {:01.3f}, mean {:01.3f}, " \ + "median {:01.3f})".format( + np.max(np.abs(tx_dpd)), + np.mean(np.abs(tx_dpd)), + np.median(np.abs(tx_dpd)) + ), size=8) + ax.legend(loc=4) + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + ccdf_min, ccdf_max = 0, 1 + tx_hist, ccdf_edges = np.histogram(np.abs(tx_dpd), + bins=60, + range=(ccdf_min, ccdf_max)) + tx_hist_normalized = tx_hist.astype(float) / np.sum(tx_hist) + ccdf = 1.0 - np.cumsum(tx_hist_normalized) + ax.semilogy(ccdf_edges[:-1], ccdf, label="CCDF") + ax.semilogy(ccdf_edges[:-1], + tx_hist_normalized, + label="Histogram", + drawstyle='steps') + ax.legend(loc=4) + ax.set_ylim(1e-5, 2) + ax.set_title("Complementary Cumulative Distribution Function") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("Ratio of Samples larger than x") + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + ax.semilogy(np.array(self.mses_am) + 1e-10, label="log(MSE)") + ax.semilogy(np.array(self.errs_am) + 1e-10, label="log(ERR)") + ax.legend(loc=4) + ax.set_title("MSE History") + ax.set_xlabel("Iterations") + ax.set_ylabel("MSE") + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + ax.plot(self.tx_mers, label="TX MER") + ax.plot(self.rx_mers, label="RX MER") + ax.legend(loc=4) + ax.set_title("MER History") + ax.set_xlabel("Iterations") + ax.set_ylabel("MER") + + ax = plt.subplot(4, 2, i_sub) + rx_range = np.linspace(0, 1, num=100, dtype=np.complex64) + rx_range_dpd = self.dpd_amplitude(rx_range)[0] + rx_range_dpd_new = self.dpd_amplitude(rx_range, new_coefs_am)[0] + i_sub += 1 + ax.scatter( + np.abs(tx_choice), + np.abs(rx_choice), + s=0.1) + ax.plot(rx_range_dpd / self.coefs_am[0], rx_range, linewidth=0.25, label="current") + ax.plot(rx_range_dpd_new / self.coefs_am[0], rx_range, linewidth=0.25, label="next") + ax.set_ylim(0, 1) + ax.set_xlim(0, 1) + ax.legend() + ax.set_title("Amplifier Characteristic") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("RX Amplitude") + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + coefs_am_history = np.array(self.coefs_am_history) + for idx, coef_hist in enumerate(coefs_am_history.T): + ax.plot(coef_hist, + label="Coef {}".format(idx), + linewidth=0.5) + ax.legend(loc=4) + ax.set_title("AM/AM Coefficient History") + ax.set_xlabel("Iterations") + ax.set_ylabel("Coefficient Value") + + phase_range = np.linspace(0, 1, num=100, dtype=np.complex64) + phase_range_dpd = self.dpd_phase(phase_range)[0] + phase_range_dpd_new = self.dpd_phase(phase_range, + coefs=new_coefs_pm)[0] + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + ax.scatter( + np.abs(tx_choice), + np.rad2deg(phase_diff_choice), + s=0.1) + ax.plot( + np.abs(phase_range), + np.rad2deg(phase_range_dpd), + linewidth=0.25, + label="current") + ax.plot( + np.abs(phase_range), + np.rad2deg(phase_range_dpd_new), + linewidth=0.25, + label="next") + ax.set_ylim(-60, 60) + ax.set_xlim(0, 1) + ax.legend() + ax.set_title("Amplifier Characteristic") + ax.set_xlabel("TX Amplitude") + ax.set_ylabel("Phase Difference") + + ax = plt.subplot(4, 2, i_sub) + i_sub += 1 + coefs_pm_history = np.array(self.coefs_pm_history) + for idx, coef_phase_hist in enumerate(coefs_pm_history.T): + ax.plot(coef_phase_hist, + label="Coef {}".format(idx), + linewidth=0.5) + ax.legend(loc=4) + ax.set_title("AM/PM Coefficient History") + ax.set_xlabel("Iterations") + ax.set_ylabel("Coefficient Value") + + fig.tight_layout() + fig.savefig(fig_path) + fig.clf() + + self.coefs_am = new_coefs_am + self.coefs_am_history.append(self.coefs_am) + self.coefs_pm = new_coefs_pm + self.coefs_pm_history.append(self.coefs_pm) + return 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/Symbol_align.py b/dpd/src/Symbol_align.py new file mode 100644 index 0000000..05a9049 --- /dev/null +++ b/dpd/src/Symbol_align.py @@ -0,0 +1,191 @@ +# -*- coding: utf-8 -*- +# +# 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 time +try: + logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) +except: + logging_path = "/tmp/" + +import numpy as np +import src.const +import scipy +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + +class Symbol_align: + """ + Find the phase offset to the start of the DAB symbols in an + unaligned dab signal. + """ + def __init__(self, sample_rate): + self.c = src.const.const(sample_rate) + pass + + def _calc_offset_to_first_symbol_without_prefix(self, tx, debug=False): + 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 debug: + 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() + + # "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 _remove_outliers(self, 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(self, fft, debug=False): + # 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 = self._remove_outliers(deltas_angle) + + delta_angle = np.mean(deltas_angle) + delta_angle_std = np.std(deltas_angle) + if debug: + plt.subplot(211) + plt.plot(angles, 'p') + plt.subplot(212) + plt.plot(deltas_angle, 'p') + return delta_angle + + def _delta_angle_to_samples(self, angle): + return - angle / self.c.phase_offset_per_sample + + def _calc_sample_offset(self, sig, debug=False): + 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 = self._calc_delta_angle(fft_crop, debug=debug) + 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") + return delta_sample_int + + def calc_offset(self, tx): + off_sym = self._calc_offset_to_first_symbol_without_prefix( + tx, debug=False) + 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 new file mode 100644 index 0000000..9274612 --- /dev/null +++ b/dpd/src/TX_Agc.py @@ -0,0 +1,109 @@ +# -*- 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 + +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +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, + max_txgain=89, + tx_median_target=0.1, + tx_median_threshold_max=0.12, + tx_median_threshold_min=0.08): + """ + 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 = max_txgain + self.txgain = self.max_txgain + + assert tx_median_threshold_max > tx_median_target,\ + "The tolerated tx_median has to be larger then the goal tx_median" + self.tx_median_threshold_tolerate_max = tx_median_threshold_max + self.tx_median_threshold_tolerate_min = tx_median_threshold_min + self.tx_median_target = tx_median_target + + def adapt_if_necessary(self, tx): + tx_median = np.median(np.abs(tx)) + if tx_median > self.tx_median_threshold_tolerate_max or\ + tx_median < self.tx_median_threshold_tolerate_min: + 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 + ) + self.adapt.set_txgain(new_txgain) + txgain = self.adapt.get_txgain() + + digital_gain_factor = 10 ** (delta_db / 20.) + digital_gain = self.adapt.get_digital_gain() * digital_gain_factor + self.adapt.set_digital_gain(digital_gain) + + logging.info( + "digital_gain = {}, txgain_new = {}, "\ + "delta_db = {}, tx_median {}, "\ + "digital_gain_factor = {}". + format(digital_gain, txgain, delta_db, + tx_median, digital_gain_factor)) + + time.sleep(1) + 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/Test_data.py b/dpd/src/Test_data.py new file mode 100644 index 0000000..9dd0913 --- /dev/null +++ b/dpd/src/Test_data.py @@ -0,0 +1,136 @@ +# -*- coding: utf-8 -*- +# +# 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 time +try: + logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) +except: + logging_path = "/tmp/" + +import src.const +import src.Dab_Util +import numpy as np +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + + +class Test_data: + def __init__(self, sample_rate, type): + """ + Standardized access to complex64 test data files. + + :param sample_rate: + :param type: + + Testing: + TD = src.Test_data.Test_data(8192000, 'file') + tx_orig = TD.get_symbol(0,0) + fig = plt.figure(figsize=(9,6)) + ax = fig.add_subplot(2,1,1) + ax.plot(tx_orig) + ax = fig.add_subplot(2,1,2) + plt.plot(np.angle(np.fft.fftshift(np.fft.fft(tx_orig))), 'p') + """ + + self.c = src.const.const(sample_rate) + self.du = src.Dab_Util.Dab_Util(sample_rate) + + self.file_paths = { + (2048000, 'file'): + ("./test_data/odr-dabmod_to_file_2048_NoFir_noDPD.iq", + ( + self.c.T_F + # Pipelineing + self.c.T_NULL + # NULL Symbol + self.c.T_S + # Synchronization Symbol + self.c.T_C # Cyclic Prefix + )), + (8192000, 'file'): + ("./test_data/odr-dabmod_to_file_8192_NoFir_noDPD.iq", + ( + self.c.T_F + # Pipelining + self.c.T_U + # Pipelining Resampler TODO(?) + self.c.T_NULL + # NULL Symbol + self.c.T_S + # Synchronization Symbol + self.c.T_C # Cyclic Prefix + )), + (8192000, 'rec_noFir'): + ("./test_data/odr-dabmod_reconding_8192_NoFir_DPD_2104.iq", + ( 64 )), + (8192000, 'rec_fir'): + ("./test_data/odr-dabmod_reconding_8192_Fir_DPD_2104.iq", + ( 232 )), + } + + config = (sample_rate, type) + if not config in self.file_paths.keys(): + raise RuntimeError("Configuration not found, possible are:\n {}". + format(self.file_paths)) + + self.path, self.file_offset = self.file_paths[(sample_rate, type)] + + def _load_from_file(self, offset, length): + print(offset, length, self.file_offset) + return self.du.fromfile( + self.path, + length=length, + offset=offset + self.file_offset) + + def get_symbol_without_prefix(self, + frame_idx=0, + symbol_idx=0, + off=0): + return self._load_from_file( + frame_idx*self.c.T_F + + symbol_idx*self.c.T_S + + off, + self.c.T_U) + + def get_symbol_with_prefix(self, + frame_idx=0, + symbol_idx=0, + n_symbols=1, + off=0): + offset = ( + frame_idx*self.c.T_F + + symbol_idx*self.c.T_S - + self.c.T_C + + off) + return self._load_from_file( offset, self.c.T_S * n_symbols) + + def get_file_length_in_symbols(self): + symbol_size = float( + 64/8 * # complex 64 + self.c.T_S # Symbol Size + ) + return os.path.getsize(self.path) / symbol_size + + +# 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 new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/dpd/src/__init__.py diff --git a/dpd/src/const.py b/dpd/src/const.py new file mode 100644 index 0000000..1011c6c --- /dev/null +++ b/dpd/src/const.py @@ -0,0 +1,38 @@ +# DAB Frame constants +# Sources: +# - etsi_EN_300_401_v010401p p145 +# - Measured with USRP B200 + +import numpy as np + +class const: + def __init__(self, sample_rate): + self.sample_rate = sample_rate + + # Time domain + self.T_F = sample_rate / 2048000 * 196608 # Transmission frame duration + self.T_NULL = sample_rate / 2048000 * 2656 # Null symbol duration + self.T_S = sample_rate / 2048000 * 2552 # Duration of OFDM symbols of indices l = 1, 2, 3,... L; + self.T_U = sample_rate / 2048000 * 2048 # Inverse of carrier spacing + self.T_C = sample_rate / 2048000 * 504 # Duration of cyclic prefix + + # Frequency Domain + # example: np.delete(fft[3328:4865], 768) + self.FFT_delete = 768 + self.FFT_delta = 1536 # Number of carrier frequencies + if sample_rate == 2048000: + self.FFT_start = 256 + self.FFT_end = 1793 + elif sample_rate == 8192000: + self.FFT_start = 3328 + self.FFT_end = 4865 + else: + raise RuntimeError("Sample Rate '{}' not supported".format( + sample_rate + )) + + # 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. / sample_rate * 2 * np.pi * 1000 diff --git a/dpd/src/phase_align.py b/dpd/src/phase_align.py new file mode 100644 index 0000000..7f82392 --- /dev/null +++ b/dpd/src/phase_align.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation 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 +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +import numpy as np +from scipy import signal, optimize +import sys +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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and plot: + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_phase_align.svg" + + 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.clf() + + 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 new file mode 100755 index 0000000..6d1cd2a --- /dev/null +++ b/dpd/src/subsample_align.py @@ -0,0 +1,112 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation 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 os +import logging +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +import numpy as np +from scipy import signal, 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=False): + """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: + corr = np.vectorize(correlate_for_delay) + ixs = np.linspace(-1, 1, 100) + taus = corr(ixs) + + dt = datetime.datetime.now().isoformat() + tau_path = (logging_path + "/" + dt + "_tau.svg") + plt.plot(ixs, taus) + plt.title("Subsample correlation, minimum is best: {}".format(best_tau)) + plt.savefig(tau_path) + plt.clf() + + # 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/src/test_dab_Util.py b/dpd/src/test_dab_Util.py new file mode 100644 index 0000000..ec15586 --- /dev/null +++ b/dpd/src/test_dab_Util.py @@ -0,0 +1,92 @@ +# -*- coding: utf-8 -*- +# +# Test code for DAB util +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +from unittest import TestCase + +import numpy as np +import pandas as pd +import src.Dab_Util as DU + +class TestDab_Util(TestCase): + + def test_subsample_align(self, sample_orig=r'../test_data/orig_rough_aligned.dat', + sample_rec =r'../test_data/recored_rough_aligned.dat', + length = 10240, max_size = 1000000): + du = DU.Dab_Util(8196000) + res1 = [] + res2 = [] + for i in range(10): + start = np.random.randint(50, max_size) + r = np.random.randint(-50, 50) + + s1 = du.fromfile(sample_orig, offset=start+r, length=length) + s2 = du.fromfile(sample_rec, offset=start, length=length) + + res1.append(du.lag_upsampling(s2, s1, 32)) + + s1_aligned, s2_aligned = du.subsample_align(s1, s2) + + res2.append(du.lag_upsampling(s2_aligned, s1_aligned, 32)) + + error_rate = np.mean(np.array(res2) != 0) + self.assertEqual(error_rate, 0.0, "The error rate for aligning was %.2f%%" + % error_rate * 100) + +#def test_using_aligned_pair(sample_orig=r'../data/orig_rough_aligned.dat', sample_rec =r'../data/recored_rough_aligned.dat', length = 10240, max_size = 1000000): +# res = [] +# for i in tqdm(range(100)): +# start = np.random.randint(50, max_size) +# r = np.random.randint(-50, 50) +# +# s1 = du.fromfile(sample_orig, offset=start+r, length=length) +# s2 = du.fromfile(sample_rec, offset=start, length=length) +# +# res.append({'offset':r, +# '1':r - du.lag_upsampling(s2, s1, n_up=1), +# '2':r - du.lag_upsampling(s2, s1, n_up=2), +# '3':r - du.lag_upsampling(s2, s1, n_up=3), +# '4':r - du.lag_upsampling(s2, s1, n_up=4), +# '8':r - du.lag_upsampling(s2, s1, n_up=8), +# '16':r - du.lag_upsampling(s2, s1, n_up=16), +# '32':r - du.lag_upsampling(s2, s1, n_up=32), +# }) +# df = pd.DataFrame(res) +# df = df.reindex_axis(sorted(df.columns), axis=1) +# print(df.describe()) +# +# +#print("Align using upsampling") +#for n_up in [1, 2, 3, 4, 7, 8, 16]: +# correct_ratio = test_phase_offset(lambda x,y: du.lag_upsampling(x,y,n_up), tol=1./n_up) +# print("%.1f%% of the tested offsets were measured within tolerance %.4f for n_up = %d" % (correct_ratio * 100, 1./n_up, n_up)) +#test_using_aligned_pair() +# +#print("Phase alignment") +#test_subsample_alignment() + + +# 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/test_measure.py b/dpd/src/test_measure.py new file mode 100644 index 0000000..ee3cfdb --- /dev/null +++ b/dpd/src/test_measure.py @@ -0,0 +1,62 @@ +# -*- coding: utf-8 -*- +# +# DPD Calculation Engine, test case for measure +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file +from unittest import TestCase +from Measure import Measure +import socket + + +class TestMeasure(TestCase): + + def _open_socks(self): + sock_server = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + sock_server.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1) + sock_server.bind(('localhost', 1234)) + sock_server.listen(1) + + sock_client = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + sock_client.connect(('localhost', 1234)) + + conn_server, addr_server = sock_server.accept() + return conn_server, sock_client + + def test__recv_exact(self): + m = Measure(1234, 1) + payload = b"test payload" + + conn_server, sock_client = self._open_socks() + conn_server.send(payload) + rec = m._recv_exact(sock_client, len(payload)) + + self.assertEqual(rec, payload, + "Did not receive the same message as sended. (%s, %s)" % + (rec, payload)) + + def test_get_samples(self): + self.fail() + + +# 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. |