Add DPD page to web gui

9 files changed, 539 insertions, 2 deletions

diff --git a/gui/README.md b/gui/README.md
index 112f6d9..0094a8d 100644
--- a/gui/README.md
+++ b/gui/README.md
@@ -29,6 +29,7 @@ Todo
 * Integrate DPDCE
   * Show DPD settings and effect visually
   * Allow load/store of DPD settings
+  * Make ports configurable
 * Use Feedback Server interface and make spectrum and constellation plots
 * Get authentication to work
 * Read and write config file, and add forms to change ODR-DabMod configuration
diff --git a/gui/api/__init__.py b/gui/api/__init__.py
index 1bc95a6..959dedc 100755
--- a/gui/api/__init__.py
+++ b/gui/api/__init__.py
@@ -31,8 +31,9 @@ import io
 import datetime
 
 class API:
-    def __init__(self, mod_rc):
+    def __init__(self, mod_rc, dpd):
         self.mod_rc = mod_rc
+        self.dpd = dpd
 
     @cherrypy.expose
     def index(self):
@@ -60,3 +61,13 @@ class API:
             cherrypy.response.headers["Content-Type"] = "application/json"
             cherrypy.response.status = 400
             return json.dumps("POST only").encode()
+
+    @cherrypy.expose
+    def trigger_capture(self, **kwargs):
+        if cherrypy.request.method == 'POST':
+            cherrypy.response.headers["Content-Type"] = "text/plain"
+            return self.dpd.capture_samples()
+        else:
+            cherrypy.response.headers["Content-Type"] = "text/plain"
+            cherrypy.response.status = 400
+            return "POST only"
diff --git a/gui/dpd/Align.py b/gui/dpd/Align.py
new file mode 100644
index 0000000..1634ec8
--- /dev/null
+++ b/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/gui/dpd/Capture.py b/gui/dpd/Capture.py
new file mode 100644
index 0000000..d6d0307
--- /dev/null
+++ b/gui/dpd/Capture.py
@@ -0,0 +1,197 @@
+# -*- 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
+import logging
+import numpy as np
+from scipy import signal
+
+from . import Align as sa
+
+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
+
+class Capture:
+    """Capture samples from ODR-DabMod"""
+    def __init__(self, samplerate, port, num_samples_to_request):
+        self.samplerate = samplerate
+        self.sizeof_sample = 8 # complex floats
+        self.port = port
+        self.num_samples_to_request = num_samples_to_request
+
+        # 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_start = 0.0
+        self.binning_end = 1.0
+        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
+
+        # axis 1: 0=tx, 1=rx
+        self.accumulated_samples = np.zeros((0, 2), dtype=np.complex64)
+
+    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, tx_median, 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))
+        tx_median = np.median(np.abs(txframe))
+        rxframe = rxframe / rx_median * tx_median
+
+        txframe_aligned, rxframe_aligned = 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 num_accumulated(self):
+        return self.accumulated_samples.shape[0]
+
+    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)
+
+        binned_sample_pairs = {}
+
+        minsize = self.num_samples_to_request
+
+        for i, (tx_start, tx_end) in enumerate(zip(bin_edges, bin_edges[1:])):
+            indices = np.bitwise_and(tx_start < txframe, txframe <= tx_end)
+            binned_sample_pairs[i] = (txframe[indices], rxframe[indices])
+            len_bin = len(txframe[indices])
+
+            #TODO this doesn't work, the min is always 0
+            if len_bin < minsize:
+                minsize = len_bin
+
+        # axis 0: bins, axis 1: sample index, axis 2: tx(0) and rx(1)
+        samples = np.zeros((self.binning_n_bins, len_bin, 2), dtype=np.complex64)
+
+        for i in binned_sample_pairs:
+            tx, rx = binned_sample_pairs[i]
+            new_samples = np.array((tx[:minsize], rx[:minsize]), dtype=np.complex64)
+            self.accumulated_samples = np.concatenate((self.accumulated_samples, new_samples.T))
+
diff --git a/gui/dpd/__init__.py b/gui/dpd/__init__.py
new file mode 100644
index 0000000..c0c76cc
--- /dev/null
+++ b/gui/dpd/__init__.py
@@ -0,0 +1,60 @@
+# -*- 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
+
+class DPD:
+    def __init__(self, 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
+
+        port = 50055
+        samples_to_capture = 81920
+        self.capture = Capture.Capture(self.samplerate, port, samples_to_capture)
+
+    def capture_samples(self):
+        """Captures samples and store them in the accumulated samples,
+        returns a string with some info"""
+        txframe_aligned, tx_ts, tx_median, rxframe_aligned, rx_ts, rx_median = self.capture.get_samples()
+
+        num_acc = self.capture.num_accumulated()
+
+        return "Captured {} with median {}/{}, accumulated {} samples".format(
+                len(txframe_aligned), tx_median, rx_median, num_acc)
+
+        # 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))
diff --git a/gui/run.py b/gui/run.py
index c5a0568..94f6def 100755
--- a/gui/run.py
+++ b/gui/run.py
@@ -28,6 +28,7 @@ import argparse
 from jinja2 import Environment, FileSystemLoader
 from api import API
 import zmqrc
+import dpd
 
 env = Environment(loader=FileSystemLoader('templates'))
 
@@ -36,7 +37,8 @@ class Root:
         self.config_file = config_file
         self.conf = configuration.Configuration(self.config_file)
         self.mod_rc = zmqrc.ModRemoteControl("localhost")
-        self.api = API(self.mod_rc)
+        self.dpd = dpd.DPD()
+        self.api = API(self.mod_rc, self.dpd)
 
     @cherrypy.expose
     def index(self):
@@ -66,6 +68,12 @@ class Root:
         js = ["js/odr-modulator.js"]
         return tmpl.render(tab='modulator', js=js, is_login=False)
 
+    @cherrypy.expose
+    def predistortion(self):
+        tmpl = env.get_template("predistortion.html")
+        js = ["js/odr-predistortion.js"]
+        return tmpl.render(tab='predistortion', js=js, is_login=False)
+
 if __name__ == '__main__':
     parser = argparse.ArgumentParser(description='ODR-DabMod Web GUI')
     parser.add_argument('-c', '--config',
diff --git a/gui/static/js/odr-predistortion.js b/gui/static/js/odr-predistortion.js
new file mode 100644
index 0000000..d385f6c
--- /dev/null
+++ b/gui/static/js/odr-predistortion.js
@@ -0,0 +1,49 @@
+//   Copyright (C) 2018
+//   Matthias P. Braendli, matthias.braendli@mpb.li
+//
+//    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/>.
+
+$(function(){
+    $('#capturebutton').click(function() {
+        $.ajax({
+            type: "POST",
+            url: "/api/trigger_capture",
+            contentType: 'text/plain',
+            dataType: '',
+            data: '',
+
+            error: function(data) {
+                $.gritter.add({
+                    title: 'Capture trigger',
+                    text: "ERROR",
+                    image: '/fonts/warning.png',
+                    sticky: true,
+                });
+            },
+            success: function(data) {
+                $('#capturestatus').val(data);
+            },
+        });
+    });
+});
+
+
+// ToolTip init
+$(function(){
+    $('[data-toggle="tooltip"]').tooltip();
+});
diff --git a/gui/templates/body-nav.html b/gui/templates/body-nav.html
index 064bc5b..f403619 100644
--- a/gui/templates/body-nav.html
+++ b/gui/templates/body-nav.html
@@ -9,6 +9,7 @@
                                 <li{% if tab == 'home' %} class="active"{% endif %}><a href="/">Home</a></li>
                                 <li{% if tab == 'rcvalues' %} class="active"{% endif %}><a href="/rcvalues">RC values</a></li>
                                 <li{% if tab == 'modulator' %} class="active"{% endif %}><a href="/modulator">Modulator</a></li>
+                                <li{% if tab == 'predistortion' %} class="active"{% endif %}><a href="/predistortion">Predistortion</a></li>
                                 <li class="dropdown{% if tab == 'help' or tab == 'about' %} active"{% endif %}">
                                         <a class="dropdown-toggle" data-toggle="dropdown" href="#">Help
                                         <span class="caret"></span></a>
diff --git a/gui/templates/predistortion.html b/gui/templates/predistortion.html
new file mode 100644
index 0000000..6828f14
--- /dev/null
+++ b/gui/templates/predistortion.html
@@ -0,0 +1,44 @@
+<!--
+   Copyright (C) 2018
+   Matthias P. Braendli, matthias.braendli@mpb.li
+
+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/>.
+-->
+
+<!DOCTYPE html>
+<html lang="en">
+
+{% include 'head.html' %}
+
+<body>
+ {% include 'body-nav.html' %}
+
+ <div class="container-fluid">
+  <div class="panel-group">
+   <div class="panel panel-default">
+    <div class="panel-body">
+     <h3>Capture</h3>
+     <div class="form-group">
+      <label for="capturestatus">Capture status:</label>
+      <input type="text" class="form-control" id="capturestatus">
+      <button type="button" class="btn btn-sm btn-info" id="capturebutton">Capture</button>
+     </div>
+    </div>
+   </div>
+  </div>
+ </div>
+</body>
+</html>