Merge branch 'next_memless' into next

4 files changed, 401 insertions, 0 deletions

diff --git a/dpd/README.md b/dpd/README.md
new file mode 100644
index 0000000..ec7cec2
--- /dev/null
+++ b/dpd/README.md
@@ -0,0 +1,20 @@
+Digital Predistortion for ODR-DabMod
+====================================
+
+This folder contains work in progress for digital predistortion. It requires:
+
+- USRP B200.
+- Power amplifier.
+- A feedback connection from the power amplifier output, at an appropriate power level for the B200.
+  Usually this is done with a directional coupler.
+- ODR-DabMod with enabled dpd_port, and with a samplerate of 8192000 samples per second.
+- Synchronous=1 so that the USRP has the timestamping set properly, internal refclk and pps
+  are sufficient for this example.
+- A live mux source with TIST enabled.
+
+See dpd/dpd.ini for an example.
+
+TODO
+----
+
+Implement a PA model that updates the predistorter.
diff --git a/dpd/dpd.ini b/dpd/dpd.ini
new file mode 100644
index 0000000..5e809e5
--- /dev/null
+++ b/dpd/dpd.ini
@@ -0,0 +1,44 @@
+[remotecontrol]
+telnet=1
+telnetport=2121
+zmqctrl=1
+zmqctrlendpoint=tcp://127.0.0.1:9400
+
+[log]
+syslog=0
+filelog=0
+filename=/dev/stderr
+
+[input]
+transport=tcp
+source=localhost:9200
+
+[modulator]
+digital_gain=0.8
+rate=8192000
+
+[firfilter]
+enabled=0
+
+[output]
+output=uhd
+
+[uhdoutput]
+device=
+master_clock_rate=32768000
+type=b200
+txgain=75
+channel=13C
+refclk_source=internal
+pps_source=none
+behaviour_refclk_lock_lost=ignore
+max_gps_holdover_time=600
+dpd_port=50055
+rxgain=0
+
+[delaymanagement]
+; Use synchronous=1 so that the USRP time is set. This works
+; even in the absence of a reference clk and PPS
+synchronous=1
+mutenotimestamps=1
+offset=4.0
diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py
new file mode 100755
index 0000000..e92c1d0
--- /dev/null
+++ b/dpd/show_spectrum.py
@@ -0,0 +1,180 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This is an example tool that shows how to connect to ODR-DabMod's dpd TCP server
+# and get samples from there.
+#
+# Since the TX and RX samples are not perfectly aligned, the tool has to align them properly,
+# which is done in two steps: First on sample-level using a correlation, then with subsample
+# accuracy using a FFT approach.
+#
+# It requires SciPy and matplotlib.
+#
+# Copyright (C) 2017 Matthias P. Braendli
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import numpy as np
+import matplotlib.pyplot as pp
+from matplotlib.animation import FuncAnimation
+import argparse
+
+SIZEOF_SAMPLE = 8 # complex floats
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot the spectrum of ODR-DabMod's DPD feedback")
+    parser.add_argument('--samps', default='10240', help='Number of samples to request at once',
+            required=False)
+    parser.add_argument('--port', default='50055',
+            help='port to connect to ODR-DabMod DPD (default: 50055)',
+            required=False)
+
+    parser.add_argument('--animated', action='store_true', help='Enable real-time animation')
+
+    cli_args = parser.parse_args()
+
+    port = int(cli_args.port)
+    num_samps_to_request = int(cli_args.samps)
+
+    if cli_args.animated:
+        plot_spectrum_animated(port, num_samps_to_request)
+    else:
+        plot_spectrum_once(port, num_samps_to_request)
+
+def recv_exact(sock, num_bytes):
+    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 get_samples(port, num_samps_to_request):
+    """Connect to ODR-DabMod, retrieve TX and RX samples, load
+    into numpy arrays, and return a tuple
+    (tx_timestamp, tx_samples, rx_timestamp, rx_samples)
+    where the timestamps are doubles, and the samples are numpy
+    arrays of complex floats, both having the same size
+    """
+
+    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+    s.connect(('localhost', port))
+
+    print("Send version");
+    s.sendall(b"\x01")
+
+    print("Send request for {} samples".format(num_samps_to_request))
+    s.sendall(struct.pack("=I", num_samps_to_request))
+
+    print("Wait for TX metadata")
+    num_samps, tx_second, tx_pps = struct.unpack("=III", recv_exact(s, 12))
+    tx_ts = tx_second + tx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} TX samples".format(num_samps))
+        txframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        txframe = np.fromstring(txframe_bytes, dtype=np.complex64)
+    else:
+        txframe = np.array([], dtype=np.complex64)
+
+    print("Wait for RX metadata")
+    rx_second, rx_pps = struct.unpack("=II", recv_exact(s, 8))
+    rx_ts = rx_second + rx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} RX samples".format(num_samps))
+        rxframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64)
+    else:
+        rxframe = np.array([], dtype=np.complex64)
+
+    print("Disconnecting")
+    s.close()
+
+    return (tx_ts, txframe, rx_ts, rxframe)
+
+def get_spectrum(port, num_samps_to_request):
+    tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request)
+
+    # convert to complex doubles for more dynamic range
+    txframe = txframe.astype(np.complex128)
+    rxframe = rxframe.astype(np.complex128)
+
+    print("Received {} & {} frames at {} and {}".format(
+        len(txframe), len(rxframe), tx_ts, rx_ts))
+
+    print("Calculate TX and RX spectrum assuming 8192000 samples per second")
+    tx_spectrum = np.fft.fftshift(np.fft.fft(txframe, fft_size))
+    tx_power = 20*np.log10(np.abs(tx_spectrum))
+
+    rx_spectrum = np.fft.fftshift(np.fft.fft(rxframe, fft_size))
+    rx_power = 20*np.log10(np.abs(rx_spectrum))
+    return tx_power, rx_power
+
+
+sampling_rate = 8192000
+fft_size = 4096
+freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate))
+
+def plot_spectrum_once(port, num_samps_to_request):
+    tx_power, rx_power = get_spectrum(port, num_samps_to_request)
+    fig = pp.figure()
+
+    fig.suptitle("TX and RX spectrum")
+    ax1 = fig.add_subplot(211)
+    ax1.set_title("TX")
+    ax1.plot(freqs, tx_power, 'r')
+    ax2 = fig.add_subplot(212)
+    ax2.set_title("RX")
+    ax2.plot(freqs, rx_power, 'b')
+    pp.show()
+
+def plot_spectrum_animated(port, num_samps_to_request):
+    fig, axes = pp.subplots(2, sharex=True)
+    line1, = axes[0].plot(freqs, np.ones(len(freqs)), 'r', animated=True)
+    axes[0].set_title("TX")
+    line2, = axes[1].plot(freqs, np.ones(len(freqs)), 'b', animated=True)
+    axes[1].set_title("RX")
+    lines = [line1, line2]
+
+    axes[0].set_ylim(-30, 50)
+    axes[1].set_ylim(-60, 40)
+
+    def update(frame):
+        tx_power, rx_power = get_spectrum(port, num_samps_to_request)
+
+        lines[0].set_ydata(tx_power)
+        lines[1].set_ydata(rx_power)
+        return lines
+
+    ani = FuncAnimation(fig, update, blit=True)
+    pp.show()
+
+main()
+
+# The MIT License (MIT)
+#
+# 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/dpd/store_received.py b/dpd/store_received.py
new file mode 100755
index 0000000..902f607
--- /dev/null
+++ b/dpd/store_received.py
@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This is an example tool that shows how to connect to ODR-DabMod's dpd TCP server
+# and get samples from there.
+#
+# Since the TX and RX samples are not perfectly aligned, the tool has to align them properly,
+# which is done in two steps: First on sample-level using a correlation, then with subsample
+# accuracy using a FFT approach.
+#
+# It requires SciPy and matplotlib.
+#
+# Copyright (C) 2017 Matthias P. Braendli
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import numpy as np
+import matplotlib.pyplot as pp
+from matplotlib.animation import FuncAnimation
+import argparse
+import os
+import time
+import src.dab_util as du
+
+SIZEOF_SAMPLE = 8 # complex floats
+
+def main():
+    parser = argparse.ArgumentParser(description="Plot the spectrum of ODR-DabMod's DPD feedback")
+    parser.add_argument('--samps', default='10240', help='Number of samples to request at once',
+            required=False)
+    parser.add_argument('--port', default='50055',
+            help='port to connect to ODR-DabMod DPD (default: 50055)',
+            required=False)
+    parser.add_argument('--out_dir', default='/tmp/record', help='Output directory',
+            required=False)
+    parser.add_argument('--count', default='1', help='Number of recordings',
+            required=False)
+    parser.add_argument('--verbose', type=int, default=0, help='Level of verbosity',
+            required=False)
+
+    parser.add_argument('--animated', action='store_true', help='Enable real-time animation')
+
+
+    cli_args = parser.parse_args()
+
+    if not os.path.isdir(cli_args.out_dir):
+        os.mkdir(cli_args.out_dir)
+
+    port = int(cli_args.port)
+    num_samps_to_request = int(cli_args.samps)
+
+    for i in range(int(cli_args.count)):
+        if i>0:
+            time.sleep(0.1)
+
+        tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request)
+
+        txframe_aligned, rxframe_aligned = du.subsample_align(txframe, rxframe)
+
+        if cli_args.verbose >= 1:
+            n_up = 32
+            lag = du.lag_upsampling(txframe, rxframe, n_up)
+            lag_aligned = du.lag_upsampling(txframe_aligned, rxframe_aligned, n_up)
+            print("Lag from %d times oversampled signal:" % n_up)
+            print("Before alignment: %.2f" % lag)
+            print("After alignment: %.2f" % lag_aligned)
+            print("")
+
+        txframe_aligned.tofile("%s/%d_tx_record.iq" % (cli_args.out_dir, i))
+        rxframe_aligned.tofile("%s/%d_rx_record.iq" % (cli_args.out_dir, i))
+
+
+def recv_exact(sock, num_bytes):
+    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 get_samples(port, num_samps_to_request):
+    """Connect to ODR-DabMod, retrieve TX and RX samples, load
+    into numpy arrays, and return a tuple
+    (tx_timestamp, tx_samples, rx_timestamp, rx_samples)
+    where the timestamps are doubles, and the samples are numpy
+    arrays of complex floats, both having the same size
+    """
+
+    s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+    s.connect(('localhost', port))
+
+    print("Send version");
+    s.sendall(b"\x01")
+
+    print("Send request for {} samples".format(num_samps_to_request))
+    s.sendall(struct.pack("=I", num_samps_to_request))
+
+    print("Wait for TX metadata")
+    num_samps, tx_second, tx_pps = struct.unpack("=III", recv_exact(s, 12))
+    tx_ts = tx_second + tx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} TX samples".format(num_samps))
+        txframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        txframe = np.fromstring(txframe_bytes, dtype=np.complex64)
+    else:
+        txframe = np.array([], dtype=np.complex64)
+
+    print("Wait for RX metadata")
+    rx_second, rx_pps = struct.unpack("=II", recv_exact(s, 8))
+    rx_ts = rx_second + rx_pps / 16384000.0
+
+    if num_samps > 0:
+        print("Receiving {} RX samples".format(num_samps))
+        rxframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE)
+        rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64)
+    else:
+        rxframe = np.array([], dtype=np.complex64)
+
+    print("Disconnecting")
+    s.close()
+
+    return (tx_ts, txframe, rx_ts, rxframe)
+
+
+sampling_rate = 8192000
+fft_size = 4096
+freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate))
+
+main()
+
+# The MIT License (MIT)
+#
+# 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.