aboutsummaryrefslogtreecommitdiffstats
path: root/python/gui/dpd
diff options
context:
space:
mode:
Diffstat (limited to 'python/gui/dpd')
-rw-r--r--python/gui/dpd/Align.py166
-rw-r--r--python/gui/dpd/Capture.py253
-rw-r--r--python/gui/dpd/__init__.py93
3 files changed, 512 insertions, 0 deletions
diff --git a/python/gui/dpd/Align.py b/python/gui/dpd/Align.py
new file mode 100644
index 0000000..1634ec8
--- /dev/null
+++ b/python/gui/dpd/Align.py
@@ -0,0 +1,166 @@
+# -*- 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/python/gui/dpd/Capture.py b/python/gui/dpd/Capture.py
new file mode 100644
index 0000000..7d95f90
--- /dev/null
+++ b/python/gui/dpd/Capture.py
@@ -0,0 +1,253 @@
+# -*- 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/python/gui/dpd/__init__.py b/python/gui/dpd/__init__.py
new file mode 100644
index 0000000..9009436
--- /dev/null
+++ b/python/gui/dpd/__init__.py
@@ -0,0 +1,93 @@
+# -*- 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))
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-21 14:53:56 +0000