From 253be52c23528544d54a59b649a60193fffb2848 Mon Sep 17 00:00:00 2001 From: andreas128 Date: Thu, 28 Sep 2017 18:59:35 +0200 Subject: Cleanup --- dpd/apply_adapt_dumps.py | 8 +-- dpd/main.py | 20 ++++-- dpd/src/Adapt.py | 2 +- dpd/src/Agc.py | 165 ------------------------------------------- dpd/src/Const.py | 81 +++++++++++---------- dpd/src/Dab_Util.py | 44 ++++++------ dpd/src/ExtractStatistic.py | 78 ++++++++++---------- dpd/src/Measure_Shoulders.py | 4 ++ dpd/src/Model_AM.py | 48 +++++++------ dpd/src/RX_Agc.py | 165 +++++++++++++++++++++++++++++++++++++++++++ dpd/src/phase_align.py | 2 - dpd/src/test_dab_Util.py | 2 +- 12 files changed, 314 insertions(+), 305 deletions(-) delete mode 100644 dpd/src/Agc.py create mode 100644 dpd/src/RX_Agc.py diff --git a/dpd/apply_adapt_dumps.py b/dpd/apply_adapt_dumps.py index 4dac657..4ecdc08 100755 --- a/dpd/apply_adapt_dumps.py +++ b/dpd/apply_adapt_dumps.py @@ -42,7 +42,7 @@ import src.Measure as Measure import src.Model as Model import src.ExtractStatistic as ExtractStatistic import src.Adapt as Adapt -import src.Agc as Agc +import src.RX_Agc as Agc import src.TX_Agc as TX_Agc import argparse @@ -113,7 +113,7 @@ MS = src.Measure_Shoulders.Measure_Shoulders(c) meas = Measure.Measure(samplerate, port, num_req) extStat = ExtractStatistic.ExtractStatistic(c) adapt = Adapt.Adapt(port_rc, coef_path) -dpddata = adapt.get_predistorter() +dpddata = adapt.get_predistorter if cli_args.lut: model = Model.Lut(c) @@ -128,7 +128,7 @@ tx_gain = adapt.get_txgain() rx_gain = adapt.get_rxgain() digital_gain = adapt.get_digital_gain() -dpddata = adapt.get_predistorter() +dpddata = adapt.get_predistorter if dpddata[0] == "poly": coefs_am = dpddata[1] coefs_pm = dpddata[2] @@ -162,7 +162,7 @@ print(paths) for i, path in enumerate(paths): print(i, path) adapt.load(path) - dpddata_after = adapt.get_predistorter() + dpddata_after = adapt.get_predistorter coefs_am, coefs_pm = model.reset_coefs() adapt.set_predistorter(("poly", coefs_am, coefs_pm)) diff --git a/dpd/main.py b/dpd/main.py index 176f80f..4e9647a 100755 --- a/dpd/main.py +++ b/dpd/main.py @@ -5,6 +5,7 @@ # # http://www.opendigitalradio.org # Licence: The MIT License, see notice at the end of this file +# noinspection PyBroadException """This Python script is the main file for ODR-DabMod's DPD Computation Engine. This engine calculates and updates the parameter of the digital @@ -15,8 +16,7 @@ import os import argparse import matplotlib -matplotlib.use('GTKAgg') - +matplotlib.use('Agg') parser = argparse.ArgumentParser( description="DPD Computation Engine for ODR-DabMod") @@ -83,7 +83,8 @@ import logging dt = datetime.datetime.now().isoformat() logging_path = '/tmp/dpd_{}'.format(dt).replace('.', '_').replace(':', '-') -if name: logging_path += '_' + name +if name: + logging_path += '_' + name os.makedirs(logging_path) logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s', datefmt='%Y-%m-%d %H:%M:%S', @@ -109,7 +110,7 @@ import src.Heuristics as Heuristics from src.Measure import Measure from src.ExtractStatistic import ExtractStatistic from src.Adapt import Adapt -from src.Agc import Agc +from src.RX_Agc import Agc from src.TX_Agc import TX_Agc from src.Symbol_align import Symbol_align from src.Const import Const @@ -198,6 +199,11 @@ while i < num_iter: # Model elif state == 'model': # Calculate new model parameters and delete old measurements + if any([x is None for x in [tx, rx, phase_diff]]): + logging.error("No data to calculate model") + state = 'measure' + continue + lr = Heuristics.get_learning_rate(i) model.train(tx, rx, phase_diff, lr=lr) dpddata = model.get_dpd_data() @@ -224,8 +230,8 @@ while i < num_iter: rx_gain = adapt.get_rxgain() digital_gain = adapt.get_digital_gain() tx_median = np.median(np.abs(txframe_aligned)) - rx_shoulder_tuple = MS.average_shoulders(rxframe_aligned) if c.MS_enable else None - tx_shoulder_tuple = MS.average_shoulders(txframe_aligned) if c.MS_enable else None + rx_shoulder_tuple = MS.average_shoulders(rxframe_aligned) + tx_shoulder_tuple = MS.average_shoulders(txframe_aligned) # Generic logging logging.info(list((name, eval(name)) for name in @@ -253,7 +259,7 @@ while i < num_iter: i += 1 state = 'measure' - except Exception as e: + except: logging.error('Iteration {} failed.'.format(i)) logging.error(traceback.format_exc()) diff --git a/dpd/src/Adapt.py b/dpd/src/Adapt.py index cf8b6f5..1ab7710 100644 --- a/dpd/src/Adapt.py +++ b/dpd/src/Adapt.py @@ -219,7 +219,7 @@ class Adapt: "txgain": self.get_txgain(), "rxgain": self.get_rxgain(), "digital_gain": self.get_digital_gain(), - "predistorter": self.get_predistorter + "predistorter": self.get_predistorter() } with open(path, "w") as f: pickle.dump(d, f) diff --git a/dpd/src/Agc.py b/dpd/src/Agc.py deleted file mode 100644 index 670fbbb..0000000 --- a/dpd/src/Agc.py +++ /dev/null @@ -1,165 +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 -logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) - -import numpy as np -import matplotlib -matplotlib.use('agg') -import matplotlib.pyplot as plt - -import src.Adapt as Adapt -import src.Measure as Measure - -class Agc: - """The goal of the automatic gain control is to set the - RX gain to a value at which all received amplitudes can - be detected. This means that the maximum possible amplitude - should be quantized at the highest possible digital value. - - A problem we have to face, is that the estimation of the - maximum amplitude by applying the max() function is very - unstable. This is due to the maximum’s rareness. Therefore - we estimate a far more robust value, such as the median, - and then approximate the maximum amplitude from it. - - Given this, we tune the RX gain in such a way, that the - received signal fulfills our desired property, of having - all amplitudes properly quantized.""" - - def __init__(self, measure, adapt, 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.""" - self.adapt.set_rxgain(self.min_rxgain) - time.sleep(1) - - dt = datetime.datetime.now().isoformat() - fig_path = logging_path + "/" + dt + "_agc.svg" - fig, axs = plt.subplots(2, 2, figsize=(3*6,1*6)) - axs = axs.ravel() - - for j in range(5): - txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median =\ - self.measure.get_samples() - - rxframe_aligned_abs = np.abs(rxframe_aligned) - - x = np.arange(100, rxframe_aligned_abs.shape[0], dtype=int) - rx_max_until = [] - rx_median_until = [] - rx_mean_until = [] - for i in x: - rx_max_until.append(np.max(rxframe_aligned_abs[:i])) - rx_median_until.append(np.median(rxframe_aligned_abs[:i])) - rx_mean_until.append(np.mean(rxframe_aligned_abs[:i])) - - axs[0].plot(x, - rx_max_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.8,0.7)), - linestyle="-", linewidth=0.25) - axs[0].set_xlabel("Samples") - axs[0].set_ylabel("Amplitude") - axs[0].set_title("Estimation for Maximum RX Amplitude") - axs[0].legend() - - axs[1].plot(x, - rx_median_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), - linestyle="-", linewidth=0.25) - axs[1].set_xlabel("Samples") - axs[1].set_ylabel("Amplitude") - axs[1].set_title("Estimation for Median RX Amplitude") - axs[1].legend() - ylim_1 = axs[1].get_ylim() - - axs[2].plot(x, - rx_mean_until, - label="Run {}".format(j+1), - color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), - linestyle="-", linewidth=0.25) - axs[2].set_xlabel("Samples") - axs[2].set_ylabel("Amplitude") - axs[2].set_title("Estimation for Mean RX Amplitude") - ylim_2 = axs[2].get_ylim() - axs[2].legend() - - axs[1].set_ylim(min(ylim_1[0], ylim_2[0]), - max(ylim_1[1], ylim_2[1])) - - fig.tight_layout() - fig.savefig(fig_path) - - axs[3].hist(rxframe_aligned_abs, bins=60) - axs[3].set_xlabel("Amplitude") - axs[3].set_ylabel("Frequency") - axs[3].set_title("Histogram of Amplitudes") - axs[3].legend() - - fig.tight_layout() - fig.savefig(fig_path) - 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/Const.py b/dpd/src/Const.py index ebdfeb2..390dfbb 100644 --- a/dpd/src/Const.py +++ b/dpd/src/Const.py @@ -1,16 +1,21 @@ -# DAB Frame constants -# Sources: -# - etsi_EN_300_401_v010401p p145 -# - Measured with USRP B200 +# -*- coding: utf-8 -*- +# +# DPD Calculation Engine, constants +# +# 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 Const: def __init__(self, sample_rate, target_median, plot): + assert sample_rate == 8192000 # By now only constants for 8192000 self.sample_rate = sample_rate - self.tx_gain_max = 89 - + # DAB frame # Time domain self.T_F = sample_rate / 2048000 * 196608 # Transmission frame duration self.T_NULL = sample_rate / 2048000 * 2656 # Null symbol duration @@ -20,18 +25,10 @@ class Const: # Frequency Domain # example: np.delete(fft[3328:4865], 768) - self.FFT_delete = 768 self.FFT_delta = 1536 # Number of carrier frequencies - if sample_rate == 2048000: - self.FFT_start = 256 - self.FFT_end = 1793 - elif sample_rate == 8192000: - self.FFT_start = 3328 - self.FFT_end = 4865 - else: - raise RuntimeError("Sample Rate '{}' not supported".format( - sample_rate - )) + self.FFT_delete = 768 + self.FFT_start = 3328 + self.FFT_end = 4865 # Calculate sample offset from phase rotation # time per sample = 1 / sample_rate @@ -43,42 +40,44 @@ class Const: self.ES_plot = plot self.ES_start = 0.0 self.ES_end = 1.0 - self.ES_n_bins = 64 - self.ES_n_per_bin = 128 - - # Constants for TX_Agc - self.TAGC_max_txgain = 89 - self.TAGC_tx_median_target = 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 - - self.RAGC_min_rxgain = 25 - self.RAGC_rx_median_target = self.TAGC_tx_median_target - - # Constants for Model - self.MDL_plot = True or plot # Override default - - # Constants for MER - self.MER_plot = plot - - # Constants for Model_PM - self.MPM_tx_min = 0.1 + 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 = plot - assert sample_rate==8192000 - meas_offset = 976 # Offset from center frequency to measure shoulder [kHz] - meas_width = 100 # Size of frequency delta to measure shoulder [kHz] + + 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_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 = plot + + # Constants for Model + self.MDL_plot = True or plot # Override default + + # 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 specific + self.TAGC_tx_median_target = 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 specific + self.RAGC_rx_median_target = self.TAGC_tx_median_target diff --git a/dpd/src/Dab_Util.py b/dpd/src/Dab_Util.py index 37be5db..27a31ef 100644 --- a/dpd/src/Dab_Util.py +++ b/dpd/src/Dab_Util.py @@ -8,30 +8,40 @@ import datetime import os import logging + logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) import numpy as np -import scipy import matplotlib + matplotlib.use('agg') import matplotlib.pyplot as plt import src.subsample_align as sa import src.phase_align as pa from scipy import signal + +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, sample_rate, plot=False): """ :param sample_rate: sample rate [sample/sec] to use for calculations """ self.sample_rate = sample_rate - self.dab_bandwidth = 1536000 #Bandwidth of a dab signal - self.frame_ms = 96 #Duration of a Dab frame + self.dab_bandwidth = 1536000 # Bandwidth of a dab signal + self.frame_ms = 96 # Duration of a Dab frame - self.plot=plot + self.plot = plot def lag(self, sig_orig, sig_rec): """ @@ -56,10 +66,10 @@ class Dab_Util: 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) + sig_rec_up = signal.resample(sig_rec, sig_rec.shape[0] * n_up) else: sig_orig_up = sig_orig - sig_rec_up = sig_rec + sig_rec_up = sig_rec l = self.lag(sig_orig_up, sig_rec_up) l_orig = float(l) / n_up return l_orig @@ -69,7 +79,7 @@ class Dab_Util: 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]) + assert (sig_tx.shape[0] == sig_rx.shape[0]) if sig_tx.shape[0] % 2 == 1: sig_tx = sig_tx[:-1] @@ -124,11 +134,11 @@ class Dab_Util: 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)) + len(sig_tx), + sig_tx.dtype, + len(sig_rx), + sig_rx.dtype, + off_meas)) if off_meas > 0: sig_tx = sig_tx[:-off] @@ -211,16 +221,10 @@ class Dab_Util: 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)) + logging.debug( + "Sig1_cut: %d %s, Sig2_cut: %d %s, off: %d" % (len(sig_tx), sig_tx.dtype, len(sig_rx), sig_rx.dtype, off)) return sig_tx, sig_rx - def fromfile(self, filename, offset=0, length=None): - if length is None: - return np.memmap(filename, dtype=np.complex64, mode='r', offset=64/8*offset) - else: - return np.memmap(filename, dtype=np.complex64, mode='r', offset=64/8*offset, shape=length) - - # The MIT License (MIT) # # Copyright (c) 2017 Andreas Steger diff --git a/dpd/src/ExtractStatistic.py b/dpd/src/ExtractStatistic.py index 9df85bc..8ea849b 100644 --- a/dpd/src/ExtractStatistic.py +++ b/dpd/src/ExtractStatistic.py @@ -1,13 +1,12 @@ # -*- coding: utf-8 -*- # # DPD Calculation Engine, -# Extract statistic from data to use in Model +# 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 pickle import matplotlib.pyplot as plt import datetime @@ -30,6 +29,14 @@ def _check_input_extract(tx_dpd, 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""" @@ -37,31 +44,27 @@ class ExtractStatistic: 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.tx_values = self._tx_value_per_bin() - - self.rx_values = [] - for i in range(c.ES_n_bins): - self.rx_values.append(None) - self.plot = c.ES_plot - def _plot_and_log(self): + def _plot_and_log(self, tx_values, rx_values, phase_diffs_values, phase_diffs_values_lists): if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: - phase_diffs_values_lists = self._phase_diff_list_per_bin() - phase_diffs_values = self._phase_diff_value_per_bin(phase_diffs_values_lists) dt = datetime.datetime.now().isoformat() fig_path = logging_path + "/" + dt + "_ExtractStatistic.png" @@ -72,13 +75,13 @@ class ExtractStatistic: i_sub += 1 ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(self.tx_values, self.rx_values, + ax.plot(tx_values, rx_values, label="Estimated Values", color="red") - for i, tx_value in enumerate(self.tx_values): - rx_values = self.rx_values_lists[i] - ax.scatter(np.ones(len(rx_values)) * tx_value, - np.abs(rx_values), + 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") @@ -90,10 +93,10 @@ class ExtractStatistic: i_sub += 1 ax = plt.subplot(sub_rows, sub_cols, i_sub) - ax.plot(self.tx_values, np.rad2deg(phase_diffs_values), + ax.plot(tx_values, np.rad2deg(phase_diffs_values), label="Estimated Values", color="red") - for i, tx_value in enumerate(self.tx_values): + 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), @@ -101,14 +104,14 @@ class ExtractStatistic: color="black") ax.set_xlabel("TX Amplitude") ax.set_ylabel("Phase Difference") - ax.set_ylim(-60,60) + ax.set_ylim(-60, 60) ax.set_xlim(0, 1.1) ax.legend(loc=4) num = [] - for i, tx_value in enumerate(self.tx_values): - rx_values = self.rx_values_lists[i] - num.append(len(rx_values)) + 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) @@ -120,9 +123,6 @@ class ExtractStatistic: fig.savefig(fig_path) plt.close(fig) - pickle.dump(self.rx_values_lists, open("/tmp/rx_values", "wb")) - pickle.dump(self.tx_values, open("/tmp/tx_values", "wb")) - def _rx_value_per_bin(self): rx_values = [] for values in self.rx_values_lists: @@ -145,14 +145,9 @@ class ExtractStatistic: phase_values_lists.append(phase_diffs) return phase_values_lists - def _phase_diff_value_per_bin(self, 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 - 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 @@ -165,23 +160,22 @@ class ExtractStatistic: self.tx_values_lists[i] += \ list(tx_dpd[mask][:n_add]) - self.rx_values = self._rx_value_per_bin() - self.tx_values = self._tx_value_per_bin() - - self._plot_and_log() + 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) - # TODO cleanup phase_diffs_values_lists = self._phase_diff_list_per_bin() - phase_diffs_values = self._phase_diff_value_per_bin(phase_diffs_values_lists) + 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) - return np.array(self.tx_values, dtype=np.float32)[:idx_end], \ - np.array(self.rx_values, dtype=np.float32)[:idx_end], \ - np.array(phase_diffs_values, dtype=np.float32)[:idx_end], \ - n_per_bin + 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) # diff --git a/dpd/src/Measure_Shoulders.py b/dpd/src/Measure_Shoulders.py index 710e800..7d48a2b 100644 --- a/dpd/src/Measure_Shoulders.py +++ b/dpd/src/Measure_Shoulders.py @@ -101,6 +101,10 @@ class Measure_Shoulders: plt.close(fig) def average_shoulders(self, signal, n_avg=None): + if not self.c.MS_enable: + logging.debug("Shoulder Measurement disabled") + return None + assert signal.shape[0] > 4 * self.c.MS_FFT_size if n_avg is None: n_avg = self.c.MS_averaging_size diff --git a/dpd/src/Model_AM.py b/dpd/src/Model_AM.py index 85f6495..cdc3de1 100644 --- a/dpd/src/Model_AM.py +++ b/dpd/src/Model_AM.py @@ -15,14 +15,30 @@ 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_float32 = lambda x: (isinstance(x, np.ndarray) and - x.dtype == np.float32 and - x.flags.contiguous) - assert is_float32(tx_dpd), \ - "tx_dpd is not float32 but {}".format(tx_dpd[0].dtype) - assert is_float32(rx_received), \ - "rx_received is not float32 but {}".format(tx_dpd[0].dtype) + 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)[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: @@ -40,8 +56,8 @@ class Model_AM: def _plot(self, tx_dpd, rx_received, coefs_am, coefs_am_new): if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot: - tx_range, rx_est = self.calc_line(coefs_am, 0, 0.6) - tx_range_new, rx_est_new = self.calc_line(coefs_am_new, 0, 0.6) + 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 = logging_path + "/" + dt + "_Model_AM.svg" @@ -66,7 +82,6 @@ class Model_AM: ax.set_title("Model_AM") ax.set_xlabel("TX Amplitude") ax.set_ylabel("RX Amplitude") - xlim = ax.get_xlim() ax.set_xlim(-0.5, 1.5) ax.legend(loc=4) @@ -74,21 +89,10 @@ class Model_AM: fig.savefig(fig_path) plt.close(fig) - def poly(self, sig): - return np.array([sig ** i for i in range(1, 6)]).T - - def fit_poly(self, tx_abs, rx_abs): - return np.linalg.lstsq(self.poly(rx_abs), tx_abs)[0] - - def calc_line(self, coefs, min_amp, max_amp): - rx_range = np.linspace(min_amp, max_amp) - tx_est = np.sum(self.poly(rx_range) * coefs, axis=1) - return tx_est, rx_range - def get_next_coefs(self, tx_dpd, rx_received, coefs_am): check_input_get_next_coefs(tx_dpd, rx_received) - coefs_am_new = self.fit_poly(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) diff --git a/dpd/src/RX_Agc.py b/dpd/src/RX_Agc.py new file mode 100644 index 0000000..670fbbb --- /dev/null +++ b/dpd/src/RX_Agc.py @@ -0,0 +1,165 @@ +# -*- coding: utf-8 -*- +# +# Automatic Gain Control +# +# http://www.opendigitalradio.org +# Licence: The MIT License, see notice at the end of this file + +import datetime +import os +import logging +import time +logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) + +import numpy as np +import matplotlib +matplotlib.use('agg') +import matplotlib.pyplot as plt + +import src.Adapt as Adapt +import src.Measure as Measure + +class Agc: + """The goal of the automatic gain control is to set the + RX gain to a value at which all received amplitudes can + be detected. This means that the maximum possible amplitude + should be quantized at the highest possible digital value. + + A problem we have to face, is that the estimation of the + maximum amplitude by applying the max() function is very + unstable. This is due to the maximum’s rareness. Therefore + we estimate a far more robust value, such as the median, + and then approximate the maximum amplitude from it. + + Given this, we tune the RX gain in such a way, that the + received signal fulfills our desired property, of having + all amplitudes properly quantized.""" + + def __init__(self, measure, adapt, 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.""" + self.adapt.set_rxgain(self.min_rxgain) + time.sleep(1) + + dt = datetime.datetime.now().isoformat() + fig_path = logging_path + "/" + dt + "_agc.svg" + fig, axs = plt.subplots(2, 2, figsize=(3*6,1*6)) + axs = axs.ravel() + + for j in range(5): + txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median =\ + self.measure.get_samples() + + rxframe_aligned_abs = np.abs(rxframe_aligned) + + x = np.arange(100, rxframe_aligned_abs.shape[0], dtype=int) + rx_max_until = [] + rx_median_until = [] + rx_mean_until = [] + for i in x: + rx_max_until.append(np.max(rxframe_aligned_abs[:i])) + rx_median_until.append(np.median(rxframe_aligned_abs[:i])) + rx_mean_until.append(np.mean(rxframe_aligned_abs[:i])) + + axs[0].plot(x, + rx_max_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.8,0.7)), + linestyle="-", linewidth=0.25) + axs[0].set_xlabel("Samples") + axs[0].set_ylabel("Amplitude") + axs[0].set_title("Estimation for Maximum RX Amplitude") + axs[0].legend() + + axs[1].plot(x, + rx_median_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), + linestyle="-", linewidth=0.25) + axs[1].set_xlabel("Samples") + axs[1].set_ylabel("Amplitude") + axs[1].set_title("Estimation for Median RX Amplitude") + axs[1].legend() + ylim_1 = axs[1].get_ylim() + + axs[2].plot(x, + rx_mean_until, + label="Run {}".format(j+1), + color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)), + linestyle="-", linewidth=0.25) + axs[2].set_xlabel("Samples") + axs[2].set_ylabel("Amplitude") + axs[2].set_title("Estimation for Mean RX Amplitude") + ylim_2 = axs[2].get_ylim() + axs[2].legend() + + axs[1].set_ylim(min(ylim_1[0], ylim_2[0]), + max(ylim_1[1], ylim_2[1])) + + fig.tight_layout() + fig.savefig(fig_path) + + axs[3].hist(rxframe_aligned_abs, bins=60) + axs[3].set_xlabel("Amplitude") + axs[3].set_ylabel("Frequency") + axs[3].set_title("Histogram of Amplitudes") + axs[3].legend() + + fig.tight_layout() + fig.savefig(fig_path) + 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/phase_align.py b/dpd/src/phase_align.py index 7317d70..0d662ec 100644 --- a/dpd/src/phase_align.py +++ b/dpd/src/phase_align.py @@ -10,8 +10,6 @@ import logging logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename) import numpy as np -from scipy import signal, optimize -import sys import matplotlib.pyplot as plt diff --git a/dpd/src/test_dab_Util.py b/dpd/src/test_dab_Util.py index ec15586..86ee2d8 100644 --- a/dpd/src/test_dab_Util.py +++ b/dpd/src/test_dab_Util.py @@ -16,7 +16,7 @@ class TestDab_Util(TestCase): def test_subsample_align(self, sample_orig=r'../test_data/orig_rough_aligned.dat', sample_rec =r'../test_data/recored_rough_aligned.dat', length = 10240, max_size = 1000000): - du = DU.Dab_Util(8196000) + du = DU res1 = [] res2 = [] for i in range(10): -- cgit v1.2.3