diff options
Diffstat (limited to 'dpd/src')
| -rw-r--r-- | dpd/src/Adapt.py | 286 | ||||
| -rw-r--r-- | dpd/src/Dab_Util.py | 246 | ||||
| -rw-r--r-- | dpd/src/ExtractStatistic.py | 196 | ||||
| -rw-r--r-- | dpd/src/GlobalConfig.py | 108 | ||||
| -rw-r--r-- | dpd/src/Heuristics.py | 56 | ||||
| -rw-r--r-- | dpd/src/MER.py | 132 | ||||
| -rw-r--r-- | dpd/src/Measure.py | 136 | ||||
| -rw-r--r-- | dpd/src/Measure_Shoulders.py | 158 | ||||
| -rw-r--r-- | dpd/src/Model.py | 32 | ||||
| -rw-r--r-- | dpd/src/Model_AM.py | 122 | ||||
| -rw-r--r-- | dpd/src/Model_Lut.py | 60 | ||||
| -rw-r--r-- | dpd/src/Model_PM.py | 124 | ||||
| -rw-r--r-- | dpd/src/Model_Poly.py | 101 | ||||
| -rw-r--r-- | dpd/src/RX_Agc.py | 166 | ||||
| -rw-r--r-- | dpd/src/Symbol_align.py | 193 | ||||
| -rw-r--r-- | dpd/src/TX_Agc.py | 131 | ||||
| -rw-r--r-- | dpd/src/__init__.py | 0 | ||||
| -rw-r--r-- | dpd/src/phase_align.py | 98 | ||||
| -rwxr-xr-x | dpd/src/subsample_align.py | 111 |
19 files changed, 0 insertions, 2456 deletions
diff --git a/dpd/src/Adapt.py b/dpd/src/Adapt.py deleted file mode 100644 index a57602f..0000000 --- a/dpd/src/Adapt.py +++ /dev/null @@ -1,286 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine: updates ODR-DabMod's settings -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file -""" -This module is used to change settings of ODR-DabMod using -the ZMQ remote control socket. -""" - -import zmq -import logging -import numpy as np -import os -import datetime -import pickle - -LUT_LEN = 32 -FORMAT_POLY = 1 -FORMAT_LUT = 2 - - -def _write_poly_coef_file(coefs_am, coefs_pm, path): - assert (len(coefs_am) == len(coefs_pm)) - - f = open(path, 'w') - f.write("{}\n{}\n".format(FORMAT_POLY, len(coefs_am))) - for coef in coefs_am: - f.write("{}\n".format(coef)) - for coef in coefs_pm: - f.write("{}\n".format(coef)) - f.close() - - -def _write_lut_file(scalefactor, lut, path): - assert (len(lut) == LUT_LEN) - - f = open(path, 'w') - f.write("{}\n{}\n".format(FORMAT_LUT, scalefactor)) - for coef in lut: - f.write("{}\n{}\n".format(coef.real, coef.imag)) - f.close() - -def dpddata_to_str(dpddata): - if dpddata[0] == "poly": - coefs_am = dpddata[1] - coefs_pm = dpddata[2] - return "dpd_coefs_am {}, dpd_coefs_pm {}".format( - coefs_am, coefs_pm) - elif dpddata[0] == "lut": - scalefactor = dpddata[1] - lut = dpddata[2] - return "LUT scalefactor {}, LUT {}".format( - scalefactor, lut) - else: - raise ValueError("Unknown dpddata type {}".format(dpddata[0])) - -class Adapt: - """Uses the ZMQ remote control to change parameters of the DabMod - - Parameters - ---------- - port : int - Port at which the ODR-DabMod is listening to connect the - ZMQ remote control. - """ - - def __init__(self, config, port, coef_path): - logging.debug("Instantiate Adapt object") - self.c = config - self.port = port - self.coef_path = coef_path - self.host = "localhost" - self._context = zmq.Context() - - def _connect(self): - """Establish the connection to ODR-DabMod using - a ZMQ socket that is in request mode (Client). - Returns a socket""" - sock = self._context.socket(zmq.REQ) - poller = zmq.Poller() - poller.register(sock, zmq.POLLIN) - - sock.connect("tcp://%s:%d" % (self.host, self.port)) - - sock.send(b"ping") - - socks = dict(poller.poll(1000)) - if socks: - if socks.get(sock) == zmq.POLLIN: - data = [el.decode() for el in sock.recv_multipart()] - - if data != ['ok']: - raise RuntimeError( - "Invalid ZMQ RC answer to 'ping' at %s %d: %s" % - (self.host, self.port, data)) - else: - sock.close(linger=10) - raise RuntimeError( - "ZMQ RC does not respond to 'ping' at %s %d" % - (self.host, self.port)) - - return sock - - def send_receive(self, message): - """Send a message to ODR-DabMod. It always - returns the answer ODR-DabMod sends back. - - An example message could be - "get sdr txgain" or "set sdr txgain 50" - - Parameter - --------- - message : str - The message string that will be sent to the receiver. - """ - sock = self._connect() - logging.debug("Send message: %s" % message) - msg_parts = message.split(" ") - for i, part in enumerate(msg_parts): - if i == len(msg_parts) - 1: - f = 0 - else: - f = zmq.SNDMORE - - sock.send(part.encode(), flags=f) - - data = [el.decode() for el in sock.recv_multipart()] - logging.debug("Received message: %s" % message) - return data - - def set_txgain(self, gain): - """Set a new txgain for the ODR-DabMod. - - Parameters - ---------- - gain : int - new TX gain, in the same format as ODR-DabMod's config file - """ - # TODO this is specific to the B200 - if gain < 0 or gain > 89: - raise ValueError("Gain has to be in [0,89]") - return self.send_receive("set sdr txgain %.4f" % float(gain)) - - def get_txgain(self): - """Get the txgain value in dB for the ODR-DabMod.""" - # TODO handle failure - return float(self.send_receive("get sdr txgain")[0]) - - def set_rxgain(self, gain): - """Set a new rxgain for the ODR-DabMod. - - Parameters - ---------- - gain : int - new RX gain, in the same format as ODR-DabMod's config file - """ - # TODO this is specific to the B200 - if gain < 0 or gain > 89: - raise ValueError("Gain has to be in [0,89]") - return self.send_receive("set sdr rxgain %.4f" % float(gain)) - - def get_rxgain(self): - """Get the rxgain value in dB for the ODR-DabMod.""" - # TODO handle failure - return float(self.send_receive("get sdr rxgain")[0]) - - def set_digital_gain(self, gain): - """Set a new rxgain for the ODR-DabMod. - - Parameters - ---------- - gain : int - new RX gain, in the same format as ODR-DabMod's config file - """ - msg = "set gain digital %.5f" % gain - return self.send_receive(msg) - - def get_digital_gain(self): - """Get the rxgain value in dB for the ODR-DabMod.""" - # TODO handle failure - return float(self.send_receive("get gain digital")[0]) - - def get_predistorter(self): - """Load the coefficients from the file in the format given in the README, - return ("poly", [AM coef], [PM coef]) or ("lut", scalefactor, [LUT entries]) - """ - f = open(self.coef_path, 'r') - lines = f.readlines() - predistorter_format = int(lines[0]) - if predistorter_format == FORMAT_POLY: - coefs_am_out = [] - coefs_pm_out = [] - n_coefs = int(lines[1]) - coefs = [float(l) for l in lines[2:]] - i = 0 - for c in coefs: - if i < n_coefs: - coefs_am_out.append(c) - elif i < 2 * n_coefs: - coefs_pm_out.append(c) - else: - raise ValueError( - 'Incorrect coef file format: too many' - ' coefficients in {}, should be {}, coefs are {}' - .format(self.coef_path, n_coefs, coefs)) - i += 1 - f.close() - return 'poly', coefs_am_out, coefs_pm_out - elif predistorter_format == FORMAT_LUT: - scalefactor = int(lines[1]) - coefs = np.array([float(l) for l in lines[2:]], dtype=np.float32) - coefs = coefs.reshape((-1, 2)) - lut = coefs[..., 0] + 1j * coefs[..., 1] - if len(lut) != LUT_LEN: - raise ValueError("Incorrect number of LUT entries ({} expected {})".format(len(lut), LUT_LEN)) - return 'lut', scalefactor, lut - else: - raise ValueError("Unknown predistorter format {}".format(predistorter_format)) - - def set_predistorter(self, dpddata): - """Update the predistorter data in the modulator. Takes the same - tuple format as argument than the one returned get_predistorter()""" - if dpddata[0] == "poly": - coefs_am = dpddata[1] - coefs_pm = dpddata[2] - _write_poly_coef_file(coefs_am, coefs_pm, self.coef_path) - elif dpddata[0] == "lut": - scalefactor = dpddata[1] - lut = dpddata[2] - _write_lut_file(scalefactor, lut, self.coef_path) - else: - raise ValueError("Unknown predistorter '{}'".format(dpddata[0])) - self.send_receive("set memlesspoly coeffile {}".format(self.coef_path)) - - def dump(self, path=None): - """Backup current settings to a file""" - dt = datetime.datetime.now().isoformat() - if path is None: - if self.c.plot_location is not None: - path = self.c.plot_location + "/" + dt + "_adapt.pkl" - else: - raise Exception("Cannot dump Adapt without either plot_location or path set") - d = { - "txgain": self.get_txgain(), - "rxgain": self.get_rxgain(), - "digital_gain": self.get_digital_gain(), - "predistorter": self.get_predistorter() - } - with open(path, "wb") as f: - pickle.dump(d, f) - - return path - - def load(self, path): - """Restore settings from a file""" - with open(path, "rb") as f: - d = pickle.load(f) - - self.set_txgain(d["txgain"]) - self.set_digital_gain(d["digital_gain"]) - self.set_rxgain(d["rxgain"]) - self.set_predistorter(d["predistorter"]) - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger, Matthias P. Braendli -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Dab_Util.py b/dpd/src/Dab_Util.py deleted file mode 100644 index bc89a39..0000000 --- a/dpd/src/Dab_Util.py +++ /dev/null @@ -1,246 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, utilities for working with DAB signals. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import numpy as np -import matplotlib - -matplotlib.use('agg') -import matplotlib.pyplot as plt -import src.subsample_align as sa -import src.phase_align as pa -from scipy import signal - - -def fromfile(filename, offset=0, length=None): - if length is None: - return np.memmap(filename, dtype=np.complex64, mode='r', offset=64 / 8 * offset) - else: - return np.memmap(filename, dtype=np.complex64, mode='r', offset=64 / 8 * offset, shape=length) - - -class Dab_Util: - """Collection of methods that can be applied to an array - complex IQ samples of a DAB signal - """ - - def __init__(self, config, sample_rate, plot=False): - """ - :param sample_rate: sample rate [sample/sec] to use for calculations - """ - self.c = config - self.sample_rate = sample_rate - self.dab_bandwidth = 1536000 # Bandwidth of a dab signal - self.frame_ms = 96 # Duration of a Dab frame - - self.plot = plot - - def lag(self, sig_orig, sig_rec): - """ - Find lag between two signals - Args: - sig_orig: The signal that has been sent - sig_rec: The signal that has been recored - """ - off = sig_rec.shape[0] - c = np.abs(signal.correlate(sig_orig, sig_rec)) - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - corr_path = self.c.plot_location + "/" + dt + "_tx_rx_corr.png" - plt.plot(c, label="corr") - plt.legend() - plt.savefig(corr_path) - plt.close() - - return np.argmax(c) - off + 1 - - def lag_upsampling(self, sig_orig, sig_rec, n_up): - if n_up != 1: - sig_orig_up = signal.resample(sig_orig, sig_orig.shape[0] * n_up) - sig_rec_up = signal.resample(sig_rec, sig_rec.shape[0] * n_up) - else: - sig_orig_up = sig_orig - sig_rec_up = sig_rec - l = self.lag(sig_orig_up, sig_rec_up) - l_orig = float(l) / n_up - return l_orig - - def subsample_align_upsampling(self, sig_tx, sig_rx, n_up=32): - """ - Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment - Using upsampling - """ - assert (sig_tx.shape[0] == sig_rx.shape[0]) - - if sig_tx.shape[0] % 2 == 1: - sig_tx = sig_tx[:-1] - sig_rx = sig_rx[:-1] - - sig1_up = signal.resample(sig_tx, sig_tx.shape[0] * n_up) - sig2_up = signal.resample(sig_rx, sig_rx.shape[0] * n_up) - - off_meas = self.lag_upsampling(sig2_up, sig1_up, n_up=1) - off = int(abs(off_meas)) - - if off_meas > 0: - sig1_up = sig1_up[:-off] - sig2_up = sig2_up[off:] - elif off_meas < 0: - sig1_up = sig1_up[off:] - sig2_up = sig2_up[:-off] - - sig_tx = signal.resample(sig1_up, sig1_up.shape[0] / n_up).astype(np.complex64) - sig_rx = signal.resample(sig2_up, sig2_up.shape[0] / n_up).astype(np.complex64) - return sig_tx, sig_rx - - def subsample_align(self, sig_tx, sig_rx): - """ - Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment - """ - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_sync_raw.png" - - fig, axs = plt.subplots(2) - axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") - axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") - axs[0].set_title("Raw Data") - axs[0].set_ylabel("Amplitude") - axs[0].set_xlabel("Samples") - axs[0].legend(loc=4) - - axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") - axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") - axs[1].set_title("Raw Data") - axs[1].set_ylabel("Real Part") - axs[1].set_xlabel("Samples") - axs[1].legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - off_meas = self.lag_upsampling(sig_rx, sig_tx, n_up=1) - off = int(abs(off_meas)) - - logging.debug("sig_tx_orig: {} {}, sig_rx_orig: {} {}, offset {}".format( - len(sig_tx), - sig_tx.dtype, - len(sig_rx), - sig_rx.dtype, - off_meas)) - - if off_meas > 0: - sig_tx = sig_tx[:-off] - sig_rx = sig_rx[off:] - elif off_meas < 0: - sig_tx = sig_tx[off:] - sig_rx = sig_rx[:-off] - - if off % 2 == 1: - sig_tx = sig_tx[:-1] - sig_rx = sig_rx[:-1] - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_sync_sample_aligned.png" - - fig, axs = plt.subplots(2) - axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") - axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") - axs[0].set_title("Sample Aligned Data") - axs[0].set_ylabel("Amplitude") - axs[0].set_xlabel("Samples") - axs[0].legend(loc=4) - - axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") - axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") - axs[1].set_ylabel("Real Part") - axs[1].set_xlabel("Samples") - axs[1].legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - sig_rx = sa.subsample_align(sig_rx, sig_tx) - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_sync_subsample_aligned.png" - - fig, axs = plt.subplots(2) - axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") - axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") - axs[0].set_title("Subsample Aligned") - axs[0].set_ylabel("Amplitude") - axs[0].set_xlabel("Samples") - axs[0].legend(loc=4) - - axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") - axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") - axs[1].set_ylabel("Real Part") - axs[1].set_xlabel("Samples") - axs[1].legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - sig_rx = pa.phase_align(sig_rx, sig_tx) - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_sync_phase_aligned.png" - - fig, axs = plt.subplots(2) - axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame") - axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame") - axs[0].set_title("Phase Aligned") - axs[0].set_ylabel("Amplitude") - axs[0].set_xlabel("Samples") - axs[0].legend(loc=4) - - axs[1].plot(np.real(sig_tx[:128]), label="TX Frame") - axs[1].plot(np.real(sig_rx[:128]), label="RX Frame") - axs[1].set_ylabel("Real Part") - axs[1].set_xlabel("Samples") - axs[1].legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - logging.debug( - "Sig1_cut: %d %s, Sig2_cut: %d %s, off: %d" % (len(sig_tx), sig_tx.dtype, len(sig_rx), sig_rx.dtype, off)) - return sig_tx, sig_rx - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/ExtractStatistic.py b/dpd/src/ExtractStatistic.py deleted file mode 100644 index 639513a..0000000 --- a/dpd/src/ExtractStatistic.py +++ /dev/null @@ -1,196 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, -# Extract statistic from received TX and RX data to use in Model -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import numpy as np -import matplotlib.pyplot as plt -import datetime -import os -import logging - - -def _check_input_extract(tx_dpd, rx_received): - # Check data type - assert tx_dpd[0].dtype == np.complex64, \ - "tx_dpd is not complex64 but {}".format(tx_dpd[0].dtype) - assert rx_received[0].dtype == np.complex64, \ - "rx_received is not complex64 but {}".format(rx_received[0].dtype) - # Check if signals have same normalization - normalization_error = np.abs(np.median(np.abs(tx_dpd)) - - np.median(np.abs(rx_received))) / ( - np.median(np.abs(tx_dpd)) + np.median(np.abs(rx_received))) - assert normalization_error < 0.01, "Non normalized signals" - - -def _phase_diff_value_per_bin(phase_diffs_values_lists): - phase_list = [] - for values in phase_diffs_values_lists: - mean = np.mean(values) if len(values) > 0 else np.nan - phase_list.append(mean) - return phase_list - - -class ExtractStatistic: - """Calculate a low variance RX value for equally spaced tx values - of a predefined range""" - - def __init__(self, c): - self.c = c - - # Number of measurements used to extract the statistic - self.n_meas = 0 - - # Boundaries for the bins - self.tx_boundaries = np.linspace(c.ES_start, c.ES_end, c.ES_n_bins + 1) - self.n_per_bin = c.ES_n_per_bin - - # List of rx values for each bin - self.rx_values_lists = [] - for i in range(c.ES_n_bins): - self.rx_values_lists.append([]) - - # List of tx values for each bin - self.tx_values_lists = [] - for i in range(c.ES_n_bins): - self.tx_values_lists.append([]) - - self.plot = c.ES_plot - - def _plot_and_log(self, tx_values, rx_values, phase_diffs_values, phase_diffs_values_lists): - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_ExtractStatistic.png" - sub_rows = 3 - sub_cols = 1 - fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) - i_sub = 0 - - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(tx_values, rx_values, - label="Estimated Values", - color="red") - for i, tx_value in enumerate(tx_values): - rx_values_list = self.rx_values_lists[i] - ax.scatter(np.ones(len(rx_values_list)) * tx_value, - np.abs(rx_values_list), - s=0.1, - color="black") - ax.set_title("Extracted Statistic") - ax.set_xlabel("TX Amplitude") - ax.set_ylabel("RX Amplitude") - ax.set_ylim(0, 0.8) - ax.set_xlim(0, 1.1) - ax.legend(loc=4) - - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(tx_values, np.rad2deg(phase_diffs_values), - label="Estimated Values", - color="red") - for i, tx_value in enumerate(tx_values): - phase_diff = phase_diffs_values_lists[i] - ax.scatter(np.ones(len(phase_diff)) * tx_value, - np.rad2deg(phase_diff), - s=0.1, - color="black") - ax.set_xlabel("TX Amplitude") - ax.set_ylabel("Phase Difference") - ax.set_ylim(-60, 60) - ax.set_xlim(0, 1.1) - ax.legend(loc=4) - - num = [] - for i, tx_value in enumerate(tx_values): - rx_values_list = self.rx_values_lists[i] - num.append(len(rx_values_list)) - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(num) - ax.set_xlabel("TX Amplitude") - ax.set_ylabel("Number of Samples") - ax.set_ylim(0, self.n_per_bin * 1.2) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - def _rx_value_per_bin(self): - rx_values = [] - for values in self.rx_values_lists: - mean = np.mean(np.abs(values)) if len(values) > 0 else np.nan - rx_values.append(mean) - return rx_values - - def _tx_value_per_bin(self): - tx_values = [] - for start, end in zip(self.tx_boundaries, self.tx_boundaries[1:]): - tx_values.append(np.mean((start, end))) - return tx_values - - def _phase_diff_list_per_bin(self): - phase_values_lists = [] - for tx_list, rx_list in zip(self.tx_values_lists, self.rx_values_lists): - phase_diffs = [] - for tx, rx in zip(tx_list, rx_list): - phase_diffs.append(np.angle(rx * tx.conjugate())) - phase_values_lists.append(phase_diffs) - return phase_values_lists - - def extract(self, tx_dpd, rx): - """Extract information from a new measurement and store them - in member variables.""" - _check_input_extract(tx_dpd, rx) - self.n_meas += 1 - - tx_abs = np.abs(tx_dpd) - for i, (tx_start, tx_end) in enumerate(zip(self.tx_boundaries, self.tx_boundaries[1:])): - mask = (tx_abs > tx_start) & (tx_abs < tx_end) - n_add = max(0, self.n_per_bin - len(self.rx_values_lists[i])) - self.rx_values_lists[i] += \ - list(rx[mask][:n_add]) - self.tx_values_lists[i] += \ - list(tx_dpd[mask][:n_add]) - - rx_values = self._rx_value_per_bin() - tx_values = self._tx_value_per_bin() - - n_per_bin = np.array([len(values) for values in self.rx_values_lists]) - # Index of first not filled bin, assumes that never all bins are filled - idx_end = np.argmin(n_per_bin == self.c.ES_n_per_bin) - - phase_diffs_values_lists = self._phase_diff_list_per_bin() - phase_diffs_values = _phase_diff_value_per_bin(phase_diffs_values_lists) - - self._plot_and_log(tx_values, rx_values, phase_diffs_values, phase_diffs_values_lists) - - tx_values_crop = np.array(tx_values, dtype=np.float32)[:idx_end] - rx_values_crop = np.array(rx_values, dtype=np.float32)[:idx_end] - phase_diffs_values_crop = np.array(phase_diffs_values, dtype=np.float32)[:idx_end] - return tx_values_crop, rx_values_crop, phase_diffs_values_crop, n_per_bin - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/GlobalConfig.py b/dpd/src/GlobalConfig.py deleted file mode 100644 index 56839fc..0000000 --- a/dpd/src/GlobalConfig.py +++ /dev/null @@ -1,108 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, constants and global configuration -# -# Source for DAB standard: etsi_EN_300_401_v010401p p145 -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import numpy as np - -class GlobalConfig: - def __init__(self, cli_args, plot_location): - self.sample_rate = cli_args.samplerate - assert self.sample_rate == 8192000 # By now only constants for 8192000 - - self.plot_location = plot_location - - # DAB frame - # Time domain - oversample = int(self.sample_rate / 2048000) - self.T_F = oversample * 196608 # Transmission frame duration - self.T_NULL = oversample * 2656 # Null symbol duration - self.T_S = oversample * 2552 # Duration of OFDM symbols of indices l = 1, 2, 3,... L; - self.T_U = oversample * 2048 # Inverse of carrier spacing - self.T_C = oversample * 504 # Duration of cyclic prefix - - # Frequency Domain - # example: np.delete(fft[3328:4865], 768) - self.FFT_delta = 1536 # Number of carrier frequencies - self.FFT_delete = 768 - self.FFT_start = 3328 - self.FFT_end = 4865 - - # Calculate sample offset from phase rotation - # time per sample = 1 / sample_rate - # frequency per bin = 1kHz - # phase difference per sample offset = delta_t * 2 * pi * delta_freq - self.phase_offset_per_sample = 1. / self.sample_rate * 2 * np.pi * 1000 - - # Constants for ExtractStatistic - self.ES_plot = cli_args.plot - self.ES_start = 0.0 - self.ES_end = 1.0 - self.ES_n_bins = 64 # Number of bins between ES_start and ES_end - self.ES_n_per_bin = 128 # Number of measurements pre bin - - # Constants for Measure_Shoulder - self.MS_enable = False - self.MS_plot = cli_args.plot - - meas_offset = 976 # Offset from center frequency to measure shoulder [kHz] - meas_width = 100 # Size of frequency delta to measure shoulder [kHz] - shoulder_offset_edge = np.abs(meas_offset - self.FFT_delta) - self.MS_shoulder_left_start = self.FFT_start - shoulder_offset_edge - meas_width / 2 - self.MS_shoulder_left_end = self.FFT_start - shoulder_offset_edge + meas_width / 2 - self.MS_shoulder_right_start = self.FFT_end + shoulder_offset_edge - meas_width / 2 - self.MS_shoulder_right_end = self.FFT_end + shoulder_offset_edge + meas_width / 2 - self.MS_peak_start = self.FFT_start + 100 # Ignore region near edges - self.MS_peak_end = self.FFT_end - 100 - - self.MS_FFT_size = 8192 - self.MS_averaging_size = 4 * self.MS_FFT_size - self.MS_n_averaging = 40 - self.MS_n_proc = 4 - - # Constants for MER - self.MER_plot = cli_args.plot - - # Constants for Model - self.MDL_plot = cli_args.plot - - # Constants for Model_PM - # Set all phase offsets to zero for TX amplitude < MPM_tx_min - self.MPM_tx_min = 0.1 - - # Constants for TX_Agc - self.TAGC_max_txgain = 89 # USRP B200 specific - self.TAGC_tx_median_target = cli_args.target_median - self.TAGC_tx_median_max = self.TAGC_tx_median_target * 1.4 - self.TAGC_tx_median_min = self.TAGC_tx_median_target / 1.4 - - # Constants for RX_AGC - self.RAGC_min_rxgain = 25 # USRP B200 specific - self.RAGC_rx_median_target = cli_args.target_median - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2017 Matthias P. Braendli -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Heuristics.py b/dpd/src/Heuristics.py deleted file mode 100644 index 21d400b..0000000 --- a/dpd/src/Heuristics.py +++ /dev/null @@ -1,56 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, heuristics we use to tune the parameters. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import numpy as np - - -def get_learning_rate(idx_run): - """Gradually reduce learning rate from lr_max to lr_min within - idx_max steps, then keep the learning rate at lr_min""" - idx_max = 10.0 - lr_min = 0.05 - lr_max = 0.4 - lr_delta = lr_max - lr_min - idx_run = min(idx_run, idx_max) - learning_rate = lr_max - lr_delta * idx_run / idx_max - return learning_rate - - -def get_n_meas(idx_run): - """Gradually increase number of measurements used to extract - a statistic from n_meas_min to n_meas_max within idx_max steps, - then keep number of measurements at n_meas_max""" - idx_max = 10.0 - n_meas_min = 10 - n_meas_max = 20 - n_meas_delta = n_meas_max - n_meas_min - idx_run = min(idx_run, idx_max) - learning_rate = n_meas_delta * idx_run / idx_max + n_meas_min - return int(np.round(learning_rate)) - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2017 Matthias P. Braendli -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/MER.py b/dpd/src/MER.py deleted file mode 100644 index 693058d..0000000 --- a/dpd/src/MER.py +++ /dev/null @@ -1,132 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, Modulation Error Rate. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import numpy as np -import matplotlib -matplotlib.use('agg') -import matplotlib.pyplot as plt - -class MER: - def __init__(self, c): - self.c = c - - self.plot = c.MER_plot - - def _calc_spectrum(self, tx): - fft = np.fft.fftshift(np.fft.fft(tx)) - return np.delete(fft[self.c.FFT_start:self.c.FFT_end], - self.c.FFT_delete) - - def _split_in_carrier(self, x, y): - if 0.5 < np.mean((np.abs(np.abs(x) - np.abs(y)) / - np.abs(np.abs(x) + np.abs(y)))): - # Constellation points are axis aligned on the Im/Re plane - x1 = x[(y < x) & (y > -x)] - y1 = y[(y < x) & (y > -x)] - - x2 = x[(y > x) & (y > -x)] - y2 = y[(y > x) & (y > -x)] - - x3 = x[(y > x) & (y < -x)] - y3 = y[(y > x) & (y < -x)] - - x4 = x[(y < x) & (y < -x)] - y4 = y[(y < x) & (y < -x)] - else: - # Constellation points are on the diagonal or Im/Re plane - x1 = x[(+x > 0) & (+y > 0)] - y1 = y[(+x > 0) & (+y > 0)] - - x2 = x[(-x > 0) & (+y > 0)] - y2 = y[(-x > 0) & (+y > 0)] - - x3 = x[(-x > 0) & (-y > 0)] - y3 = y[(-x > 0) & (-y > 0)] - - x4 = x[(+x > 0) & (-y > 0)] - y4 = y[(+x > 0) & (-y > 0)] - return (x1, y1), (x2, y2), (x3, y3), (x4, y4) - - def _calc_mer_for_isolated_constellation_point(self, x, y): - """Calculate MER for one constellation point""" - x_mean = np.mean(x) - y_mean = np.mean(y) - - U_RMS = np.sqrt(x_mean ** 2 + y_mean ** 2) - U_ERR = np.mean(np.sqrt((x - x_mean) ** 2 + (y - y_mean) ** 2)) - MER = 20 * np.log10(U_ERR / U_RMS) - - return x_mean, y_mean, U_RMS, U_ERR, MER - - def calc_mer(self, tx, debug_name=""): - """Calculate MER for input signal from a symbol aligned signal.""" - assert tx.shape[0] == self.c.T_U, "Wrong input length" - - spectrum = self._calc_spectrum(tx) - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_MER" + debug_name + ".png" - else: - fig_path = None - - MERs = [] - for i, (x, y) in enumerate(self._split_in_carrier( - np.real(spectrum), - np.imag(spectrum))): - x_mean, y_mean, U_RMS, U_ERR, MER =\ - self._calc_mer_for_isolated_constellation_point(x, y) - MERs.append(MER) - - tau = np.linspace(0, 2 * np.pi, num=100) - x_err = U_ERR * np.sin(tau) + x_mean - y_err = U_ERR * np.cos(tau) + y_mean - - if self.plot: - ax = plt.subplot(221 + i) - ax.scatter(x, y, s=0.2, color='black') - ax.plot(x_mean, y_mean, 'p', color='red') - ax.plot(x_err, y_err, linewidth=2, color='blue') - ax.text(0.1, 0.1, "MER {:.1f}dB".format(MER), transform=ax.transAxes) - ax.set_xlabel("Real") - ax.set_ylabel("Imag") - ylim = ax.get_ylim() - ax.set_ylim(ylim[0] - (ylim[1] - ylim[0]) * 0.1, ylim[1]) - - if fig_path is not None: - plt.tight_layout() - plt.savefig(fig_path) - plt.show() - plt.close() - - MER_res = 20 * np.log10(np.mean([10 ** (MER / 20) for MER in MERs])) - return MER_res - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Measure.py b/dpd/src/Measure.py deleted file mode 100644 index 6d8007d..0000000 --- a/dpd/src/Measure.py +++ /dev/null @@ -1,136 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, Measure signal using ODR-DabMod's -# DPD Server. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import socket -import struct -import numpy as np -import src.Dab_Util as DU -import os -import logging - -class Measure: - """Collect Measurement from DabMod""" - def __init__(self, config, samplerate, port, num_samples_to_request): - logging.info("Instantiate Measure object") - self.c = config - self.samplerate = samplerate - self.sizeof_sample = 8 # complex floats - self.port = port - self.num_samples_to_request = num_samples_to_request - - def _recv_exact(self, sock, num_bytes): - """Receive an exact number of bytes from a socket. This is - a wrapper around sock.recv() that can return less than the number - of requested bytes. - - Args: - sock (socket): Socket to receive data from. - num_bytes (int): Number of bytes that will be returned. - """ - bufs = [] - while num_bytes > 0: - b = sock.recv(num_bytes) - if len(b) == 0: - break - num_bytes -= len(b) - bufs.append(b) - return b''.join(bufs) - - def receive_tcp(self): - s = socket.socket(socket.AF_INET, socket.SOCK_STREAM) - s.settimeout(4) - s.connect(('localhost', self.port)) - - logging.debug("Send version") - s.sendall(b"\x01") - - logging.debug("Send request for {} samples".format(self.num_samples_to_request)) - s.sendall(struct.pack("=I", self.num_samples_to_request)) - - logging.debug("Wait for TX metadata") - num_samps, tx_second, tx_pps = struct.unpack("=III", self._recv_exact(s, 12)) - tx_ts = tx_second + tx_pps / 16384000.0 - - if num_samps > 0: - logging.debug("Receiving {} TX samples".format(num_samps)) - txframe_bytes = self._recv_exact(s, num_samps * self.sizeof_sample) - txframe = np.fromstring(txframe_bytes, dtype=np.complex64) - else: - txframe = np.array([], dtype=np.complex64) - - logging.debug("Wait for RX metadata") - rx_second, rx_pps = struct.unpack("=II", self._recv_exact(s, 8)) - rx_ts = rx_second + rx_pps / 16384000.0 - - if num_samps > 0: - logging.debug("Receiving {} RX samples".format(num_samps)) - rxframe_bytes = self._recv_exact(s, num_samps * self.sizeof_sample) - rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64) - else: - rxframe = np.array([], dtype=np.complex64) - - if logging.getLogger().getEffectiveLevel() == logging.DEBUG: - logging.debug('txframe: min {}, max {}, median {}'.format( - np.min(np.abs(txframe)), - np.max(np.abs(txframe)), - np.median(np.abs(txframe)))) - - logging.debug('rxframe: min {}, max {}, median {}'.format( - np.min(np.abs(rxframe)), - np.max(np.abs(rxframe)), - np.median(np.abs(rxframe)))) - - logging.debug("Disconnecting") - s.close() - - return txframe, tx_ts, rxframe, rx_ts - - - def get_samples(self): - """Connect to ODR-DabMod, retrieve TX and RX samples, load - into numpy arrays, and return a tuple - (txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median) - """ - - txframe, tx_ts, rxframe, rx_ts = self.receive_tcp() - - # Normalize received signal with sent signal - rx_median = np.median(np.abs(rxframe)) - rxframe = rxframe / rx_median * np.median(np.abs(txframe)) - - du = DU.Dab_Util(self.c, self.samplerate) - txframe_aligned, rxframe_aligned = du.subsample_align(txframe, rxframe) - - logging.info( - "Measurement done, tx %d %s, rx %d %s, tx aligned %d %s, rx aligned %d %s" - % (len(txframe), txframe.dtype, len(rxframe), rxframe.dtype, - len(txframe_aligned), txframe_aligned.dtype, len(rxframe_aligned), rxframe_aligned.dtype) ) - - return txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Measure_Shoulders.py b/dpd/src/Measure_Shoulders.py deleted file mode 100644 index fd90050..0000000 --- a/dpd/src/Measure_Shoulders.py +++ /dev/null @@ -1,158 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, calculate peak to shoulder difference. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import multiprocessing -import numpy as np -import matplotlib.pyplot as plt - - -def plt_next_axis(sub_rows, sub_cols, i_sub): - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - return i_sub, ax - - -def plt_annotate(ax, x, y, title=None, legend_loc=None): - ax.set_xlabel(x) - ax.set_ylabel(y) - if title is not None: - ax.set_title(title) - if legend_loc is not None: - ax.legend(loc=legend_loc) - - -def calc_fft_db(signal, offset, c): - fft = np.fft.fftshift(np.fft.fft(signal[offset:offset + c.MS_FFT_size])) - fft_db = 20 * np.log10(np.abs(fft)) - return fft_db - - -def _calc_peak(fft, c): - assert fft.shape == (c.MS_FFT_size,), fft.shape - idxs = (c.MS_peak_start, c.MS_peak_end) - peak = np.mean(fft[idxs[0]:idxs[1]]) - return peak, idxs - - -def _calc_shoulder_hight(fft_db, c): - assert fft_db.shape == (c.MS_FFT_size,), fft_db.shape - idxs_left = (c.MS_shoulder_left_start, c.MS_shoulder_left_end) - idxs_right = (c.MS_shoulder_right_start, c.MS_shoulder_right_end) - - shoulder_left = np.mean(fft_db[idxs_left[0]:idxs_left[1]]) - shoulder_right = np.mean(fft_db[idxs_right[0]:idxs_right[1]]) - - shoulder = np.mean((shoulder_left, shoulder_right)) - return shoulder, (idxs_left, idxs_right) - - -def calc_shoulder(fft, c): - peak = _calc_peak(fft, c)[0] - shoulder = _calc_shoulder_hight(fft, c)[0] - assert (peak >= shoulder), (peak, shoulder) - return peak, shoulder - - -def shoulder_from_sig_offset(arg): - signal, offset, c = arg - fft_db = calc_fft_db(signal, offset, c) - peak, shoulder = calc_shoulder(fft_db, c) - return peak - shoulder, peak, shoulder - - -class Measure_Shoulders: - """Calculate difference between the DAB signal and the shoulder hight in the - power spectrum""" - - def __init__(self, c): - self.c = c - self.plot = c.MS_plot - - def _plot(self, signal): - if self.c.plot_location is None: - return - - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_sync_subsample_aligned.png" - - fft = calc_fft_db(signal, 100, self.c) - peak, idxs_peak = _calc_peak(fft, self.c) - shoulder, idxs_sh = _calc_shoulder_hight(fft, self.c) - - sub_rows = 1 - sub_cols = 1 - fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) - i_sub = 0 - - i_sub, ax = plt_next_axis(sub_rows, sub_cols, i_sub) - ax.scatter(np.arange(fft.shape[0]), fft, s=0.1, - label="FFT", - color="red") - ax.plot(idxs_peak, (peak, peak)) - ax.plot(idxs_sh[0], (shoulder, shoulder), color='blue') - ax.plot(idxs_sh[1], (shoulder, shoulder), color='blue') - plt_annotate(ax, "Frequency", "Magnitude [dB]", None, 4) - - ax.text(100, -17, str(calc_shoulder(fft, self.c))) - - ax.set_ylim(-20, 30) - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - def average_shoulders(self, signal, n_avg=None): - if not self.c.MS_enable: - logging.info("Shoulder Measurement disabled via Const.py") - return None - - assert signal.shape[0] > 4 * self.c.MS_FFT_size - if n_avg is None: - n_avg = self.c.MS_averaging_size - - off_min = 0 - off_max = signal.shape[0] - self.c.MS_FFT_size - offsets = np.linspace(off_min, off_max, num=n_avg, dtype=int) - - args = zip( - [signal, ] * offsets.shape[0], - offsets, - [self.c, ] * offsets.shape[0] - ) - - pool = multiprocessing.Pool(self.c.MS_n_proc) - res = pool.map(shoulder_from_sig_offset, args) - shoulders_diff, shoulders, peaks = zip(*res) - - if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: - self._plot(signal) - - return np.mean(shoulders_diff), np.mean(shoulders), np.mean(peaks) - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Model.py b/dpd/src/Model.py deleted file mode 100644 index b2c303f..0000000 --- a/dpd/src/Model.py +++ /dev/null @@ -1,32 +0,0 @@ -# -*- coding: utf-8 -*- -from src.Model_Poly import Poly -from src.Model_Lut import Lut - -def select_model_from_dpddata(dpddata): - if dpddata[0] == 'lut': - return Lut - elif dpddata[0] == 'poly': - return Poly - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2017 Matthias P. Braendli -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Model_AM.py b/dpd/src/Model_AM.py deleted file mode 100644 index 75b226f..0000000 --- a/dpd/src/Model_AM.py +++ /dev/null @@ -1,122 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, model implementation for Amplitude and not Phase -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import numpy as np -import matplotlib.pyplot as plt - - -def is_npfloat32(array): - assert isinstance(array, np.ndarray), type(array) - assert array.dtype == np.float32, array.dtype - assert array.flags.contiguous - assert not any(np.isnan(array)) - - -def check_input_get_next_coefs(tx_dpd, rx_received): - is_npfloat32(tx_dpd) - is_npfloat32(rx_received) - - -def poly(sig): - return np.array([sig ** i for i in range(1, 6)]).T - - -def fit_poly(tx_abs, rx_abs): - return np.linalg.lstsq(poly(rx_abs), tx_abs, rcond=None)[0] - - -def calc_line(coefs, min_amp, max_amp): - rx_range = np.linspace(min_amp, max_amp) - tx_est = np.sum(poly(rx_range) * coefs, axis=1) - return tx_est, rx_range - - -class Model_AM: - """Calculates new coefficients using the measurement and the previous - coefficients""" - - def __init__(self, - c, - learning_rate_am=1, - plot=False): - self.c = c - - self.learning_rate_am = learning_rate_am - self.plot = plot - - def _plot(self, tx_dpd, rx_received, coefs_am, coefs_am_new): - if self.plot and self.c.plot_location is not None: - tx_range, rx_est = calc_line(coefs_am, 0, 0.6) - tx_range_new, rx_est_new = calc_line(coefs_am_new, 0, 0.6) - - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_Model_AM.png" - sub_rows = 1 - sub_cols = 1 - fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) - i_sub = 0 - - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(tx_range, rx_est, - label="Estimated TX", - alpha=0.3, - color="gray") - ax.plot(tx_range_new, rx_est_new, - label="New Estimated TX", - color="red") - ax.scatter(tx_dpd, rx_received, - label="Binned Data", - color="blue", - s=1) - ax.set_title("Model_AM") - ax.set_xlabel("TX Amplitude") - ax.set_ylabel("RX Amplitude") - ax.set_xlim(-0.5, 1.5) - ax.legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - def get_next_coefs(self, tx_dpd, rx_received, coefs_am): - """Calculate the next AM/AM coefficients using the extracted - statistic of TX and RX amplitude""" - check_input_get_next_coefs(tx_dpd, rx_received) - - coefs_am_new = fit_poly(tx_dpd, rx_received) - coefs_am_new = coefs_am + \ - self.learning_rate_am * (coefs_am_new - coefs_am) - - self._plot(tx_dpd, rx_received, coefs_am, coefs_am_new) - - return coefs_am_new - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Model_Lut.py b/dpd/src/Model_Lut.py deleted file mode 100644 index e70fdb0..0000000 --- a/dpd/src/Model_Lut.py +++ /dev/null @@ -1,60 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, model implementation using polynomial -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import os -import logging -import numpy as np - -class Lut: - """Implements a model that calculates lookup table coefficients""" - - def __init__(self, - c, - learning_rate=1., - plot=False): - """ - - :rtype: - """ - logging.debug("Initialising LUT Model") - self.c = c - self.learning_rate = learning_rate - self.plot = plot - self.reset_coefs() - - def reset_coefs(self): - self.scalefactor = 0xFFFFFFFF # max uint32_t value - self.lut = np.ones(32, dtype=np.complex64) - - def train(self, tx_abs, rx_abs, phase_diff): - pass - - def get_dpd_data(self): - return "lut", self.scalefactor, self.lut - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2017 Matthias P. Braendli -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Model_PM.py b/dpd/src/Model_PM.py deleted file mode 100644 index 7b80bf3..0000000 --- a/dpd/src/Model_PM.py +++ /dev/null @@ -1,124 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, model implementation for Amplitude and not Phase -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import numpy as np -import matplotlib.pyplot as plt - - -def is_npfloat32(array): - assert isinstance(array, np.ndarray), type(array) - assert array.dtype == np.float32, array.dtype - assert array.flags.contiguous - assert not any(np.isnan(array)) - - -def check_input_get_next_coefs(tx_dpd, phase_diff): - is_npfloat32(tx_dpd) - is_npfloat32(phase_diff) - - -class Model_PM: - """Calculates new coefficients using the measurement and the previous - coefficients""" - - def __init__(self, - c, - learning_rate_pm=1, - plot=False): - self.c = c - - self.learning_rate_pm = learning_rate_pm - self.plot = plot - - def _plot(self, tx_dpd, phase_diff, coefs_pm, coefs_pm_new): - if self.plot and self.c.plot_location is not None: - tx_range, phase_diff_est = self.calc_line(coefs_pm, 0, 0.6) - tx_range_new, phase_diff_est_new = self.calc_line(coefs_pm_new, 0, 0.6) - - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_Model_PM.png" - sub_rows = 1 - sub_cols = 1 - fig = plt.figure(figsize=(sub_cols * 6, sub_rows / 2. * 6)) - i_sub = 0 - - i_sub += 1 - ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(tx_range, phase_diff_est, - label="Estimated Phase Diff", - alpha=0.3, - color="gray") - ax.plot(tx_range_new, phase_diff_est_new, - label="New Estimated Phase Diff", - color="red") - ax.scatter(tx_dpd, phase_diff, - label="Binned Data", - color="blue", - s=1) - ax.set_title("Model_PM") - ax.set_xlabel("TX Amplitude") - ax.set_ylabel("Phase DIff") - ax.legend(loc=4) - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - def _discard_small_values(self, tx_dpd, phase_diff): - """ Assumes that the phase for small tx amplitudes is zero""" - mask = tx_dpd < self.c.MPM_tx_min - phase_diff[mask] = 0 - return tx_dpd, phase_diff - - def poly(self, sig): - return np.array([sig ** i for i in range(0, 5)]).T - - def fit_poly(self, tx_abs, phase_diff): - return np.linalg.lstsq(self.poly(tx_abs), phase_diff, rcond=None)[0] - - def calc_line(self, coefs, min_amp, max_amp): - tx_range = np.linspace(min_amp, max_amp) - phase_diff = np.sum(self.poly(tx_range) * coefs, axis=1) - return tx_range, phase_diff - - def get_next_coefs(self, tx_dpd, phase_diff, coefs_pm): - """Calculate the next AM/PM coefficients using the extracted - statistic of TX amplitude and phase difference""" - tx_dpd, phase_diff = self._discard_small_values(tx_dpd, phase_diff) - check_input_get_next_coefs(tx_dpd, phase_diff) - - coefs_pm_new = self.fit_poly(tx_dpd, phase_diff) - - coefs_pm_new = coefs_pm + self.learning_rate_pm * (coefs_pm_new - coefs_pm) - self._plot(tx_dpd, phase_diff, coefs_pm, coefs_pm_new) - - return coefs_pm_new - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Model_Poly.py b/dpd/src/Model_Poly.py deleted file mode 100644 index cdfd319..0000000 --- a/dpd/src/Model_Poly.py +++ /dev/null @@ -1,101 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, model implementation using polynomial -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import os -import logging -import numpy as np - -import src.Model_AM as Model_AM -import src.Model_PM as Model_PM - - -def assert_np_float32(x): - assert isinstance(x, np.ndarray) - assert x.dtype == np.float32 - assert x.flags.contiguous - - -def _check_input_get_next_coefs(tx_abs, rx_abs, phase_diff): - assert_np_float32(tx_abs) - assert_np_float32(rx_abs) - assert_np_float32(phase_diff) - - assert tx_abs.shape == rx_abs.shape, \ - "tx_abs.shape {}, rx_abs.shape {}".format( - tx_abs.shape, rx_abs.shape) - assert tx_abs.shape == phase_diff.shape, \ - "tx_abs.shape {}, phase_diff.shape {}".format( - tx_abs.shape, phase_diff.shape) - - -class Poly: - """Calculates new coefficients using the measurement and the previous - coefficients""" - - def __init__(self, - c, - learning_rate_am=1.0, - learning_rate_pm=1.0): - self.c = c - self.plot = c.MDL_plot - - self.learning_rate_am = learning_rate_am - self.learning_rate_pm = learning_rate_pm - - self.reset_coefs() - - self.model_am = Model_AM.Model_AM(c, plot=self.plot) - self.model_pm = Model_PM.Model_PM(c, plot=self.plot) - - def reset_coefs(self): - self.coefs_am = np.zeros(5, dtype=np.float32) - self.coefs_am[0] = 1 - self.coefs_pm = np.zeros(5, dtype=np.float32) - - def train(self, tx_abs, rx_abs, phase_diff, lr=None): - """ - :type tx_abs: np.ndarray - :type rx_abs: np.ndarray - :type phase_diff: np.ndarray - :type lr: float - """ - _check_input_get_next_coefs(tx_abs, rx_abs, phase_diff) - - if not lr is None: - self.model_am.learning_rate_am = lr - self.model_pm.learning_rate_pm = lr - - coefs_am_new = self.model_am.get_next_coefs(tx_abs, rx_abs, self.coefs_am) - coefs_pm_new = self.model_pm.get_next_coefs(tx_abs, phase_diff, self.coefs_pm) - - self.coefs_am = self.coefs_am + (coefs_am_new - self.coefs_am) * self.learning_rate_am - self.coefs_pm = self.coefs_pm + (coefs_pm_new - self.coefs_pm) * self.learning_rate_pm - - def get_dpd_data(self): - return "poly", self.coefs_am, self.coefs_pm - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/RX_Agc.py b/dpd/src/RX_Agc.py deleted file mode 100644 index f778dee..0000000 --- a/dpd/src/RX_Agc.py +++ /dev/null @@ -1,166 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Automatic Gain Control -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import time -import numpy as np -import matplotlib -matplotlib.use('agg') -import matplotlib.pyplot as plt - -import src.Adapt as Adapt -import src.Measure as Measure - -class Agc: - """The goal of the automatic gain control is to set the - RX gain to a value at which all received amplitudes can - be detected. This means that the maximum possible amplitude - should be quantized at the highest possible digital value. - - A problem we have to face, is that the estimation of the - maximum amplitude by applying the max() function is very - unstable. This is due to the maximum’s rareness. Therefore - we estimate a far more robust value, such as the median, - and then approximate the maximum amplitude from it. - - Given this, we tune the RX gain in such a way, that the - received signal fulfills our desired property, of having - all amplitudes properly quantized.""" - - def __init__(self, measure, adapt, c): - assert isinstance(measure, Measure.Measure) - assert isinstance(adapt, Adapt.Adapt) - self.measure = measure - self.adapt = adapt - self.min_rxgain = c.RAGC_min_rxgain - self.rxgain = self.min_rxgain - self.peak_to_median = 1./c.RAGC_rx_median_target - - def run(self): - self.adapt.set_rxgain(self.rxgain) - - for i in range(2): - # Measure - txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median= \ - self.measure.get_samples() - - # Estimate Maximum - rx_peak = self.peak_to_median * rx_median - correction_factor = 20*np.log10(1/rx_peak) - self.rxgain = self.rxgain + correction_factor - - assert self.min_rxgain <= self.rxgain, ("Desired RX Gain is {} which is smaller than the minimum of {}".format( - self.rxgain, self.min_rxgain)) - - logging.info("RX Median {:1.4f}, estimated peak {:1.4f}, correction factor {:1.4f}, new RX gain {:1.4f}".format( - rx_median, rx_peak, correction_factor, self.rxgain - )) - - self.adapt.set_rxgain(self.rxgain) - time.sleep(0.5) - - def plot_estimates(self): - """Plots the estimate of for Max, Median, Mean for different - number of samples.""" - if self.c.plot_location is None: - return - - self.adapt.set_rxgain(self.min_rxgain) - time.sleep(1) - - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_agc.png" - fig, axs = plt.subplots(2, 2, figsize=(3*6,1*6)) - axs = axs.ravel() - - for j in range(5): - txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median =\ - self.measure.get_samples() - - rxframe_aligned_abs = np.abs(rxframe_aligned) - - x = np.arange(100, rxframe_aligned_abs.shape[0], dtype=int) - rx_max_until = [] - rx_median_until = [] - rx_mean_until = [] - for i in x: - rx_max_until.append(np.max(rxframe_aligned_abs[:i])) - rx_median_until.append(np.median(rxframe_aligned_abs[:i])) - rx_mean_until.append(np.mean(rxframe_aligned_abs[:i])) - - axs[0].plot(x, - rx_max_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.8,0.7)), - linestyle="-", linewidth=0.25) - axs[0].set_xlabel("Samples") - axs[0].set_ylabel("Amplitude") - axs[0].set_title("Estimation for Maximum RX Amplitude") - axs[0].legend() - - axs[1].plot(x, - rx_median_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), - linestyle="-", linewidth=0.25) - axs[1].set_xlabel("Samples") - axs[1].set_ylabel("Amplitude") - axs[1].set_title("Estimation for Median RX Amplitude") - axs[1].legend() - ylim_1 = axs[1].get_ylim() - - axs[2].plot(x, - rx_mean_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), - linestyle="-", linewidth=0.25) - axs[2].set_xlabel("Samples") - axs[2].set_ylabel("Amplitude") - axs[2].set_title("Estimation for Mean RX Amplitude") - ylim_2 = axs[2].get_ylim() - axs[2].legend() - - axs[1].set_ylim(min(ylim_1[0], ylim_2[0]), - max(ylim_1[1], ylim_2[1])) - - fig.tight_layout() - fig.savefig(fig_path) - - axs[3].hist(rxframe_aligned_abs, bins=60) - axs[3].set_xlabel("Amplitude") - axs[3].set_ylabel("Frequency") - axs[3].set_title("Histogram of Amplitudes") - axs[3].legend() - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/Symbol_align.py b/dpd/src/Symbol_align.py deleted file mode 100644 index 2a17a65..0000000 --- a/dpd/src/Symbol_align.py +++ /dev/null @@ -1,193 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, Modulation Error Rate. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import numpy as np -import scipy -import matplotlib - -matplotlib.use('agg') -import matplotlib.pyplot as plt - - -def _remove_outliers(x, stds=5): - deviation_from_mean = np.abs(x - np.mean(x)) - inlier_idxs = deviation_from_mean < stds * np.std(x) - x = x[inlier_idxs] - return x - - -def _calc_delta_angle(fft): - # Introduce invariance against carrier - angles = np.angle(fft) % (np.pi / 2.) - - # Calculate Angle difference and compensate jumps - deltas_angle = np.diff(angles) - deltas_angle[deltas_angle > np.pi / 4.] = \ - deltas_angle[deltas_angle > np.pi / 4.] - np.pi / 2. - deltas_angle[-deltas_angle > np.pi / 4.] = \ - deltas_angle[-deltas_angle > np.pi / 4.] + np.pi / 2. - deltas_angle = _remove_outliers(deltas_angle) - - delta_angle = np.mean(deltas_angle) - - return delta_angle - - -class Symbol_align: - """ - Find the phase offset to the start of the DAB symbols in an - unaligned dab signal. - """ - - def __init__(self, c, plot=False): - self.c = c - self.plot = plot - pass - - def _calc_offset_to_first_symbol_without_prefix(self, tx): - tx_orig = tx[0:-self.c.T_U] - tx_cut_prefix = tx[self.c.T_U:] - - tx_product = np.abs(tx_orig - tx_cut_prefix) - tx_product_avg = np.correlate( - tx_product, - np.ones(self.c.T_C), - mode='valid') - tx_product_avg_min_filt = \ - scipy.ndimage.filters.minimum_filter1d( - tx_product_avg, - int(1.5 * self.c.T_S) - ) - peaks = np.ravel(np.where(tx_product_avg == tx_product_avg_min_filt)) - - offset = peaks[np.argmin([tx_product_avg[peak] for peak in peaks])] - - if self.plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_Symbol_align.png" - - fig = plt.figure(figsize=(9, 9)) - - ax = fig.add_subplot(4, 1, 1) - ax.plot(tx_product) - ylim = ax.get_ylim() - for peak in peaks: - ax.plot((peak, peak), (ylim[0], ylim[1])) - if peak == offset: - ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) - else: - ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) - ax.set_xlabel("Sample") - ax.set_ylabel("Conj. Product") - ax.set_title("Difference with shifted self") - - ax = fig.add_subplot(4, 1, 2) - ax.plot(tx_product_avg) - ylim = ax.get_ylim() - for peak in peaks: - ax.plot((peak, peak), (ylim[0], ylim[1])) - if peak == offset: - ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) - else: - ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) - ax.set_xlabel("Sample") - ax.set_ylabel("Conj. Product") - ax.set_title("Moving Average") - - ax = fig.add_subplot(4, 1, 3) - ax.plot(tx_product_avg_min_filt) - ylim = ax.get_ylim() - for peak in peaks: - ax.plot((peak, peak), (ylim[0], ylim[1])) - if peak == offset: - ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90) - else: - ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90) - ax.set_xlabel("Sample") - ax.set_ylabel("Conj. Product") - ax.set_title("Min Filter") - - ax = fig.add_subplot(4, 1, 4) - tx_product_crop = tx_product[peaks[0] - 50:peaks[0] + 50] - x = range(tx_product.shape[0])[peaks[0] - 50:peaks[0] + 50] - ax.plot(x, tx_product_crop) - ylim = ax.get_ylim() - ax.plot((peaks[0], peaks[0]), (ylim[0], ylim[1])) - ax.set_xlabel("Sample") - ax.set_ylabel("Conj. Product") - ax.set_title("Difference with shifted self") - - fig.tight_layout() - fig.savefig(fig_path) - plt.close(fig) - - # "offset" measures where the shifted signal matches the - # original signal. Therefore we have to subtract the size - # of the shift to find the offset of the symbol start. - return (offset + self.c.T_C) % self.c.T_S - - def _delta_angle_to_samples(self, angle): - return - angle / self.c.phase_offset_per_sample - - def _calc_sample_offset(self, sig): - assert sig.shape[0] == self.c.T_U, \ - "Input length is not a Symbol without cyclic prefix" - - fft = np.fft.fftshift(np.fft.fft(sig)) - fft_crop = np.delete(fft[self.c.FFT_start:self.c.FFT_end], self.c.FFT_delete) - delta_angle = _calc_delta_angle(fft_crop) - delta_sample = self._delta_angle_to_samples(delta_angle) - delta_sample_int = np.round(delta_sample).astype(int) - error = np.abs(delta_sample_int - delta_sample) - if error > 0.1: - raise RuntimeError("Could not calculate " \ - "sample offset. Error {}".format(error)) - return delta_sample_int - - def calc_offset(self, tx): - """Calculate the offset the first symbol""" - off_sym = self._calc_offset_to_first_symbol_without_prefix( - tx) - off_sam = self._calc_sample_offset( - tx[off_sym:off_sym + self.c.T_U]) - off = (off_sym + off_sam) % self.c.T_S - - assert self._calc_sample_offset(tx[off:off + self.c.T_U]) == 0, \ - "Failed to calculate offset" - return off - - def crop_symbol_without_cyclic_prefix(self, tx): - off = self.calc_offset(tx) - return tx[ - off: - off + self.c.T_U - ] - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/TX_Agc.py b/dpd/src/TX_Agc.py deleted file mode 100644 index 309193d..0000000 --- a/dpd/src/TX_Agc.py +++ /dev/null @@ -1,131 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, Automatic Gain Control. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file - -import datetime -import os -import logging -import time -import numpy as np -import matplotlib - -matplotlib.use('agg') -import matplotlib.pyplot as plt - -import src.Adapt as Adapt - - -# TODO fix for float tx_gain -class TX_Agc: - def __init__(self, - adapt, - c): - """ - In order to avoid digital clipping, this class increases the - TX gain and reduces the digital gain. Digital clipping happens - when the digital analog converter receives values greater than - it's maximal output. This class solves that problem by adapting - the TX gain in a way that the peaks of the TX signal are in a - specified range. The TX gain is adapted accordingly. The TX peaks - are approximated by estimating it based on the signal median. - - :param adapt: Instance of Adapt Class to update - txgain and coefficients - :param max_txgain: limit for TX gain - :param tx_median_threshold_max: if the median of TX is larger - than this value, then the digital gain is reduced - :param tx_median_threshold_min: if the median of TX is smaller - than this value, then the digital gain is increased - :param tx_median_target: The digital gain is reduced in a way that - the median TX value is expected to be lower than this value. - """ - - assert isinstance(adapt, Adapt.Adapt) - self.adapt = adapt - self.max_txgain = c.TAGC_max_txgain - self.txgain = self.max_txgain - - self.tx_median_threshold_tolerate_max = c.TAGC_tx_median_max - self.tx_median_threshold_tolerate_min = c.TAGC_tx_median_min - self.tx_median_target = c.TAGC_tx_median_target - - def _calc_new_tx_gain(self, tx_median): - delta_db = 20 * np.log10(self.tx_median_target / tx_median) - new_txgain = self.adapt.get_txgain() - delta_db - assert new_txgain < self.max_txgain, \ - "TX_Agc failed. New TX gain of {} is too large.".format( - new_txgain - ) - return new_txgain, delta_db - - def _calc_digital_gain(self, delta_db): - digital_gain_factor = 10 ** (delta_db / 20.) - digital_gain = self.adapt.get_digital_gain() * digital_gain_factor - return digital_gain, digital_gain_factor - - def _set_tx_gain(self, new_txgain): - self.adapt.set_txgain(new_txgain) - txgain = self.adapt.get_txgain() - return txgain - - def _have_to_adapt(self, tx_median): - too_large = tx_median > self.tx_median_threshold_tolerate_max - too_small = tx_median < self.tx_median_threshold_tolerate_min - return too_large or too_small - - def adapt_if_necessary(self, tx): - tx_median = np.median(np.abs(tx)) - - if self._have_to_adapt(tx_median): - # Calculate new values - new_txgain, delta_db = self._calc_new_tx_gain(tx_median) - digital_gain, digital_gain_factor = \ - self._calc_digital_gain(delta_db) - - # Set new values. - # Avoid temorary increase of output power with correct order - if digital_gain_factor < 1: - self.adapt.set_digital_gain(digital_gain) - time.sleep(0.5) - txgain = self._set_tx_gain(new_txgain) - time.sleep(1) - else: - txgain = self._set_tx_gain(new_txgain) - time.sleep(1) - self.adapt.set_digital_gain(digital_gain) - time.sleep(0.5) - - logging.info( - "digital_gain = {}, txgain_new = {}, " \ - "delta_db = {}, tx_median {}, " \ - "digital_gain_factor = {}". - format(digital_gain, txgain, delta_db, - tx_median, digital_gain_factor)) - - return True - return False - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/__init__.py b/dpd/src/__init__.py deleted file mode 100644 index e69de29..0000000 --- a/dpd/src/__init__.py +++ /dev/null diff --git a/dpd/src/phase_align.py b/dpd/src/phase_align.py deleted file mode 100644 index 8654333..0000000 --- a/dpd/src/phase_align.py +++ /dev/null @@ -1,98 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, phase-align a signal against a reference. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file -import datetime -import os -import logging -import numpy as np -import matplotlib.pyplot as plt - - -def phase_align(sig, ref_sig, plot=False): - """Do phase alignment for sig relative to the reference signal - ref_sig. - - Returns the aligned signal""" - - angle_diff = (np.angle(sig) - np.angle(ref_sig)) % (2. * np.pi) - - real_diffs = np.cos(angle_diff) - imag_diffs = np.sin(angle_diff) - - if plot and self.c.plot_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = self.c.plot_location + "/" + dt + "_phase_align.png" - - plt.subplot(511) - plt.hist(angle_diff, bins=60, label="Angle Diff") - plt.xlabel("Angle") - plt.ylabel("Count") - plt.legend(loc=4) - - plt.subplot(512) - plt.hist(real_diffs, bins=60, label="Real Diff") - plt.xlabel("Real Part") - plt.ylabel("Count") - plt.legend(loc=4) - - plt.subplot(513) - plt.hist(imag_diffs, bins=60, label="Imaginary Diff") - plt.xlabel("Imaginary Part") - plt.ylabel("Count") - plt.legend(loc=4) - - plt.subplot(514) - plt.plot(np.angle(ref_sig[:128]), label="ref_sig") - plt.plot(np.angle(sig[:128]), label="sig") - plt.xlabel("Angle") - plt.ylabel("Sample") - plt.legend(loc=4) - - real_diff = np.median(real_diffs) - imag_diff = np.median(imag_diffs) - - angle = np.angle(real_diff + 1j * imag_diff) - - logging.debug( - "Compensating phase by {} rad, {} degree. real median {}, imag median {}".format( - angle, angle*180./np.pi, real_diff, imag_diff - )) - sig = sig * np.exp(1j * -angle) - - if logging.getLogger().getEffectiveLevel() == logging.DEBUG and plot: - plt.subplot(515) - plt.plot(np.angle(ref_sig[:128]), label="ref_sig") - plt.plot(np.angle(sig[:128]), label="sig") - plt.xlabel("Angle") - plt.ylabel("Sample") - plt.legend(loc=4) - plt.tight_layout() - plt.savefig(fig_path) - plt.close() - - return sig - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. diff --git a/dpd/src/subsample_align.py b/dpd/src/subsample_align.py deleted file mode 100755 index 20ae56b..0000000 --- a/dpd/src/subsample_align.py +++ /dev/null @@ -1,111 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, utility to do subsample alignment. -# -# http://www.opendigitalradio.org -# Licence: The MIT License, see notice at the end of this file -import datetime -import logging -import os -import numpy as np -from scipy import optimize -import matplotlib.pyplot as plt - -def gen_omega(length): - if (length % 2) == 1: - raise ValueError("Needs an even length array.") - - halflength = int(length / 2) - factor = 2.0 * np.pi / length - - omega = np.zeros(length, dtype=np.float) - for i in range(halflength): - omega[i] = factor * i - - for i in range(halflength, length): - omega[i] = factor * (i - length) - - return omega - - -def subsample_align(sig, ref_sig, plot_location=None): - """Do subsample alignment for sig relative to the reference signal - ref_sig. The delay between the two must be less than sample - - Returns the aligned signal""" - - n = len(sig) - if (n % 2) == 1: - raise ValueError("Needs an even length signal.") - halflen = int(n / 2) - - fft_sig = np.fft.fft(sig) - - omega = gen_omega(n) - - def correlate_for_delay(tau): - # A subsample offset between two signals corresponds, in the frequency - # domain, to a linearly increasing phase shift, whose slope - # corresponds to the delay. - # - # Here, we build this phase shift in rotate_vec, and multiply it with - # our signal. - - rotate_vec = np.exp(1j * tau * omega) - # zero-frequency is rotate_vec[0], so rotate_vec[N/2] is the - # bin corresponding to the [-1, 1, -1, 1, ...] time signal, which - # is both the maximum positive and negative frequency. - # I don't remember why we handle it differently. - rotate_vec[halflen] = np.cos(np.pi * tau) - - corr_sig = np.fft.ifft(rotate_vec * fft_sig) - - return -np.abs(np.sum(np.conj(corr_sig) * ref_sig)) - - optim_result = optimize.minimize_scalar(correlate_for_delay, bounds=(-1, 1), method='bounded', - options={'disp': True}) - - if optim_result.success: - best_tau = optim_result.x - - if plot_location is not None: - corr = np.vectorize(correlate_for_delay) - ixs = np.linspace(-1, 1, 100) - taus = corr(ixs) - - dt = datetime.datetime.now().isoformat() - tau_path = (plot_location + "/" + dt + "_tau.png") - plt.plot(ixs, taus) - plt.title("Subsample correlation, minimum is best: {}".format(best_tau)) - plt.savefig(tau_path) - plt.close() - - # Prepare rotate_vec = fft_sig with rotated phase - rotate_vec = np.exp(1j * best_tau * omega) - rotate_vec[halflen] = np.cos(np.pi * best_tau) - return np.fft.ifft(rotate_vec * fft_sig).astype(np.complex64) - else: - # print("Could not optimize: " + optim_result.message) - return np.zeros(0, dtype=np.complex64) - -# The MIT License (MIT) -# -# Copyright (c) 2017 Andreas Steger -# -# Permission is hereby granted, free of charge, to any person obtaining a copy -# of this software and associated documentation files (the "Software"), to deal -# in the Software without restriction, including without limitation the rights -# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -# copies of the Software, and to permit persons to whom the Software is -# furnished to do so, subject to the following conditions: -# -# The above copyright notice and this permission notice shall be included in all -# copies or substantial portions of the Software. -# -# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -# SOFTWARE. |