diff options
Diffstat (limited to 'gui/dpd')
| -rw-r--r-- | gui/dpd/Align.py | 166 | ||||
| -rw-r--r-- | gui/dpd/Capture.py | 253 | ||||
| -rw-r--r-- | gui/dpd/__init__.py | 93 |
3 files changed, 0 insertions, 512 deletions
diff --git a/gui/dpd/Align.py b/gui/dpd/Align.py deleted file mode 100644 index 1634ec8..0000000 --- a/gui/dpd/Align.py +++ /dev/null @@ -1,166 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, utility to do subsample alignment. -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2018 Matthias P. Braendli -# -# http://www.opendigitalradio.org -# -# This file is part of ODR-DabMod. -# -# ODR-DabMod is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. -# -# ODR-DabMod is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODR-DabMod. If not, see <http://www.gnu.org/licenses/>. -import datetime -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) - -def phase_align(sig, ref_sig, plot_location=None): - """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_location is not None: - dt = datetime.datetime.now().isoformat() - fig_path = 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 plot_location is not None: - 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 diff --git a/gui/dpd/Capture.py b/gui/dpd/Capture.py deleted file mode 100644 index 7d95f90..0000000 --- a/gui/dpd/Capture.py +++ /dev/null @@ -1,253 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine, Capture TX signal and RX feedback using ODR-DabMod's -# DPD Server. -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2018 Matthias P. Braendli -# -# http://www.opendigitalradio.org -# -# This file is part of ODR-DabMod. -# -# ODR-DabMod is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. -# -# ODR-DabMod is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODR-DabMod. If not, see <http://www.gnu.org/licenses/>. - -import socket -import struct -import os.path -import logging -import numpy as np -from scipy import signal -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -import io - -from . import Align as sa - -def correlation_coefficient(sig_tx, sig_rx): - return np.corrcoef(sig_tx, sig_rx)[0, 1] - -def align_samples(sig_tx, sig_rx): - """ - Returns an aligned version of sig_tx and sig_rx by cropping, subsample alignment and - correct phase offset - """ - - # Coarse sample-level alignment - c = np.abs(signal.correlate(sig_rx, sig_tx)) - off_meas = np.argmax(c) - sig_tx.shape[0] + 1 - off = int(abs(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] - - # Fine subsample alignment and phase offset - sig_rx = sa.subsample_align(sig_rx, sig_tx) - sig_rx = sa.phase_align(sig_rx, sig_tx) - return sig_tx, sig_rx, abs(off_meas) - -class Capture: - """Capture samples from ODR-DabMod""" - def __init__(self, samplerate, port, num_samples_to_request, plot_dir): - self.samplerate = samplerate - self.sizeof_sample = 8 # complex floats - self.port = port - self.num_samples_to_request = num_samples_to_request - self.plot_dir = plot_dir - - # Before we run the samples through the model, we want to accumulate - # them into bins depending on their amplitude, and keep only n_per_bin - # samples to avoid that the polynomial gets overfitted in the low-amplitude - # part, which is less interesting than the high-amplitude part, where - # non-linearities become apparent. - self.binning_n_bins = 64 # Number of bins between binning_start and binning_end - self.binning_n_per_bin = 128 # Number of measurements pre bin - - self.rx_normalisation = 1.0 - - self.clear_accumulated() - - def clear_accumulated(self): - self.binning_start = 0.0 - self.binning_end = 1.0 - - # axis 0: bins - # axis 1: 0=tx, 1=rx - self.accumulated_bins = [np.zeros((0, 2), dtype=np.complex64) for i in range(self.binning_n_bins)] - - 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 _plot_spectrum(self, signal, filename, title): - fig = plt.figure() - ax = plt.subplot(1, 1, 1) - - fft = np.fft.fftshift(np.fft.fft(signal)) - fft_db = 20 * np.log10(np.abs(fft)) - - ax.plot(fft_db) - ax.set_title(title) - fig.tight_layout() - fig.savefig(os.path.join(self.plot_dir, filename)) - plt.close(fig) - - def calibrate(self): - txframe, tx_ts, rxframe, rx_ts = self.receive_tcp() - - # Normalize received signal with sent signal - tx_median = np.median(np.abs(txframe)) - rx_median = np.median(np.abs(rxframe)) - self.rx_normalisation = tx_median / rx_median - - rxframe = rxframe * self.rx_normalisation - txframe_aligned, rxframe_aligned, coarse_offset = align_samples(txframe, rxframe) - - self._plot_spectrum(rxframe[:8192], "rxframe.png", "RX Frame") - self._plot_spectrum(txframe[:8192], "txframe.png", "RX Frame") - - return tx_ts, tx_median, rx_ts, rx_median, np.abs(coarse_offset), correlation_coefficient(txframe_aligned, rxframe_aligned) - - 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, tx_median, rxframe_aligned, rx_ts, rx_median) - """ - - txframe, tx_ts, rxframe, rx_ts = self.receive_tcp() - - # Normalize received signal with calibrated normalisation - rxframe = rxframe * self.rx_normalisation - txframe_aligned, rxframe_aligned, coarse_offset = align_samples(txframe, rxframe) - self._bin_and_accumulate(txframe_aligned, rxframe_aligned) - return txframe_aligned, tx_ts, tx_median, rxframe_aligned, rx_ts, rx_median - - def bin_histogram(self): - return [b.shape[0] for b in self.accumulated_bins] - - def pointcloud_png(self): - fig = plt.figure() - ax = plt.subplot(1, 1, 1) - for b in self.accumulated_bins: - if b: - ax.scatter( - np.abs(b[0]), - np.abs(b[1]), - s=0.1, - color="black") - ax.set_title("Captured and Binned Samples") - 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) - fig.tight_layout() - fig.savefig(os.path.join(self.plot_dir, "pointcloud.png")) - plt.close(fig) - - def _bin_and_accumulate(self, txframe, rxframe): - """Bin the samples and extend the accumulated samples""" - - bin_edges = np.linspace(self.binning_start, self.binning_end, self.binning_n_bins) - - minsize = self.num_samples_to_request - - for i, (tx_start, tx_end) in enumerate(zip(bin_edges, bin_edges[1:])): - txframe_abs = np.abs(txframe) - indices = np.bitwise_and(tx_start < txframe_abs, txframe_abs <= tx_end) - txsamples = np.asmatrix(txframe[indices]) - rxsamples = np.asmatrix(rxframe[indices]) - binned_sample_pairs = np.concatenate((txsamples, rxsamples)).T - - missing_in_bin = self.binning_n_per_bin - self.accumulated_bins[i].shape[0] - num_to_append = min(missing_in_bin, binned_sample_pairs.shape[0]) - print("Handling bin {} {}-{}, {} available, {} missing".format(i, tx_start, tx_end, binned_sample_pairs.shape[0], missing_in_bin)) - if num_to_append: - print("Appending {} to bin {} with shape {}".format(num_to_append, i, self.accumulated_bins[i].shape)) - - self.accumulated_bins[i] = np.concatenate((self.accumulated_bins[i], binned_sample_pairs[:num_to_append,...])) - print("{} now has shape {}".format(i, self.accumulated_bins[i].shape)) - diff --git a/gui/dpd/__init__.py b/gui/dpd/__init__.py deleted file mode 100644 index 9009436..0000000 --- a/gui/dpd/__init__.py +++ /dev/null @@ -1,93 +0,0 @@ -# -*- coding: utf-8 -*- -# -# DPD Computation Engine module -# -# Copyright (c) 2017 Andreas Steger -# Copyright (c) 2018 Matthias P. Braendli -# -# http://www.opendigitalradio.org -# -# This file is part of ODR-DabMod. -# -# ODR-DabMod is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as -# published by the Free Software Foundation, either version 3 of the -# License, or (at your option) any later version. -# -# ODR-DabMod is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with ODR-DabMod. If not, see <http://www.gnu.org/licenses/>. - -from . import Capture -import numpy as np - -class DPD: - def __init__(self, plot_dir, samplerate=8192000): - self.samplerate = samplerate - - oversample = int(self.samplerate / 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 - - self.last_capture_info = {} - - port = 50055 - samples_to_capture = 81920 - self.capture = Capture.Capture(self.samplerate, port, samples_to_capture, plot_dir) - - def status(self): - r = {} - r['histogram'] = self.capture.bin_histogram() - r['capture'] = self.last_capture_info - return r - - def pointcloud_png(self): - return self.capture.pointcloud_png() - - def clear_accumulated(self): - return self.capture.clear_accumulated() - - def capture_calibration(self): - tx_ts, tx_median, rx_ts, rx_median, coarse_offset, correlation_coefficient = self.capture.calibrate() - result = {'status': "ok"} - result['tx_median'] = "{:.2f}dB".format(20*np.log10(tx_median)) - result['rx_median'] = "{:.2f}dB".format(20*np.log10(rx_median)) - result['tx_ts'] = tx_ts - result['rx_ts'] = rx_ts - result['coarse_offset'] = int(coarse_offset) - result['correlation'] = float(correlation_coefficient) - return result - - def capture_samples(self): - """Captures samples and store them in the accumulated samples, - returns a dict with some info""" - result = {} - try: - txframe_aligned, tx_ts, tx_median, rxframe_aligned, rx_ts, rx_median = self.capture.get_samples() - result['status'] = "ok" - result['length'] = len(txframe_aligned) - result['tx_median'] = float(tx_median) - result['rx_median'] = float(rx_median) - result['tx_ts'] = tx_ts - result['rx_ts'] = rx_ts - except ValueError as e: - result['status'] = "Capture failed: {}".format(e) - - self.last_capture_info = result - - # tx, rx, phase_diff, n_per_bin = extStat.extract(txframe_aligned, rxframe_aligned) - # off = SA.calc_offset(txframe_aligned) - # print("off {}".format(off)) - # tx_mer = MER.calc_mer(txframe_aligned[off:off + c.T_U], debug_name='TX') - # print("tx_mer {}".format(tx_mer)) - # rx_mer = MER.calc_mer(rxframe_aligned[off:off + c.T_U], debug_name='RX') - # print("rx_mer {}".format(rx_mer)) - # mse = np.mean(np.abs((txframe_aligned - rxframe_aligned) ** 2)) - # print("mse {}".format(mse)) |