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