From 201d2cd2e0431a5ea79fb69561c27555f3a03dc1 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 10:44:33 +0200 Subject: Add dpd example script to fetch samples --- dpd/README.md | 10 ++++++ dpd/show_spectrum.py | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+) create mode 100644 dpd/README.md create mode 100755 dpd/show_spectrum.py (limited to 'dpd') diff --git a/dpd/README.md b/dpd/README.md new file mode 100644 index 0000000..828a483 --- /dev/null +++ b/dpd/README.md @@ -0,0 +1,10 @@ +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. diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py new file mode 100755 index 0000000..6c489e0 --- /dev/null +++ b/dpd/show_spectrum.py @@ -0,0 +1,99 @@ +#!/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 as mp +import scipy.signal + +SIZEOF_SAMPLE = 8 # complex floats + +if len(sys.argv) != 3: + print("Usage: show_spectrum.py ") + sys.exit(1) + +port = int(sys.argv[1]) +num_samps_to_request = int(sys.argv[2]) + + +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", s.recv(12)) + tx_ts = tx_second + tx_pps / 16384000.0 + + print("Receiving {} TX samples".format(num_samps)) + txframe_bytes = s.recv(num_samps * SIZEOF_SAMPLE) + txframe = np.fromstring(txframe_bytes, dtype=np.complex64) + + print("Wait for RX metadata") + rx_second, rx_pps = struct.unpack("=II", s.recv(8)) + rx_ts = rx_second + rx_pps / 16384000.0 + + print("Receiving {} RX samples".format(num_samps)) + rxframe_bytes = s.recv(num_samps * SIZEOF_SAMPLE) + rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64) + + print("Disconnecting") + s.close() + + return (tx_ts, txframe, rx_ts, rxframe) + + +tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) + +print("Received {} & {} frames at {} and {}".format( + len(txframe), len(rxframe), tx_ts, rx_ts)) + + +# 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. -- cgit v1.2.3 From a656cee6c9c230bb921c6bb6be0f0180460a96b4 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 11:23:15 +0200 Subject: DPD: handle incomplete frames --- dpd/show_spectrum.py | 31 +++++++++++++++++++++++-------- src/OutputUHD.cpp | 2 +- src/OutputUHDFeedback.cpp | 25 ++++++++++++++++++++++--- 3 files changed, 46 insertions(+), 12 deletions(-) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index 6c489e0..c1d5fe5 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -30,6 +30,15 @@ if len(sys.argv) != 3: port = int(sys.argv[1]) num_samps_to_request = int(sys.argv[2]) +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 @@ -49,20 +58,26 @@ def get_samples(port, 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", s.recv(12)) + num_samps, tx_second, tx_pps = struct.unpack("=III", recv_exact(s, 12)) tx_ts = tx_second + tx_pps / 16384000.0 - print("Receiving {} TX samples".format(num_samps)) - txframe_bytes = s.recv(num_samps * SIZEOF_SAMPLE) - txframe = np.fromstring(txframe_bytes, dtype=np.complex64) + 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", s.recv(8)) + rx_second, rx_pps = struct.unpack("=II", recv_exact(s, 8)) rx_ts = rx_second + rx_pps / 16384000.0 - print("Receiving {} RX samples".format(num_samps)) - rxframe_bytes = s.recv(num_samps * SIZEOF_SAMPLE) - rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64) + 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: + txframe = np.array([], dtype=np.complex64) print("Disconnecting") s.close() diff --git a/src/OutputUHD.cpp b/src/OutputUHD.cpp index 6edf7df..5e9e17c 100644 --- a/src/OutputUHD.cpp +++ b/src/OutputUHD.cpp @@ -50,7 +50,7 @@ using namespace std; // Maximum number of frames that can wait in uwd.frames -static const size_t FRAMES_MAX_SIZE = 2; +static const size_t FRAMES_MAX_SIZE = 8; typedef std::complex complexf; diff --git a/src/OutputUHDFeedback.cpp b/src/OutputUHDFeedback.cpp index 09b73ba..3ef5648 100644 --- a/src/OutputUHDFeedback.cpp +++ b/src/OutputUHDFeedback.cpp @@ -138,6 +138,8 @@ void OutputUHDFeedback::ReceiveBurstThread() if (not m_running) break; + etiLog.level(debug) << "Prepare RX stream command for " << burstRequest.num_samples; + uhd::stream_cmd_t cmd( uhd::stream_cmd_t::stream_mode_t::STREAM_MODE_NUM_SAMPS_AND_DONE); cmd.num_samps = burstRequest.num_samples; @@ -146,16 +148,28 @@ void OutputUHDFeedback::ReceiveBurstThread() double pps = burstRequest.rx_pps / 16384000.0; cmd.time_spec = uhd::time_spec_t(burstRequest.rx_second, pps); + const double usrp_time = m_usrp->get_time_now().get_real_secs(); + const double cmd_time = cmd.time_spec.get_real_secs(); + + etiLog.level(debug) << + "RX stream command ts=" << std::fixed << cmd_time << " Delta=" << cmd_time - usrp_time; + rxStream->issue_stream_cmd(cmd); uhd::rx_metadata_t md; burstRequest.rx_samples.resize(burstRequest.num_samples * sizeof(complexf)); - rxStream->recv(&burstRequest.rx_samples[0], burstRequest.num_samples, md); + size_t samples_read = rxStream->recv(&burstRequest.rx_samples[0], burstRequest.num_samples, md); + assert(samples_read <= burstRequest.num_samples); + burstRequest.rx_samples.resize(samples_read * sizeof(complexf)); // The recv might have happened at another time than requested burstRequest.rx_second = md.time_spec.get_full_secs(); burstRequest.rx_pps = md.time_spec.get_frac_secs() * 16384000.0; + etiLog.level(debug) << "Read " << samples_read << " RX feedback samples " + << "at time " << std::fixed << burstRequest.tx_second << "." << + burstRequest.tx_pps / 16384000.0; + burstRequest.state = BurstRequestState::Acquired; lock.unlock(); @@ -258,6 +272,11 @@ void OutputUHDFeedback::ServeFeedbackThread() burstRequest.state = BurstRequestState::None; lock.unlock(); + burstRequest.num_samples = std::min(burstRequest.num_samples, + std::min( + burstRequest.tx_samples.size() / sizeof(complexf), + burstRequest.rx_samples.size() / sizeof(complexf))); + if (send(client_sock, &burstRequest.num_samples, sizeof(burstRequest.num_samples), @@ -287,7 +306,7 @@ void OutputUHDFeedback::ServeFeedbackThread() const size_t frame_bytes = burstRequest.num_samples * sizeof(complexf); - assert(burstRequest.tx_samples.size() == frame_bytes); + assert(burstRequest.tx_samples.size() >= frame_bytes); if (send(client_sock, &burstRequest.tx_samples[0], frame_bytes, @@ -315,7 +334,7 @@ void OutputUHDFeedback::ServeFeedbackThread() break; } - assert(burstRequest.rx_samples.size() == frame_bytes); + assert(burstRequest.rx_samples.size() >= frame_bytes); if (send(client_sock, &burstRequest.rx_samples[0], frame_bytes, -- cgit v1.2.3 From eefe0ff989243b5af65fc6af0448fa4578fc713e Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 11:25:39 +0200 Subject: Update DPD readme and add example --- dpd/README.md | 8 ++++++++ dpd/dpd.ini | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+) create mode 100644 dpd/dpd.ini (limited to 'dpd') diff --git a/dpd/README.md b/dpd/README.md index 828a483..96c4fb0 100644 --- a/dpd/README.md +++ b/dpd/README.md @@ -8,3 +8,11 @@ This folder contains work in progress for digital predistortion. It requires: - 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. + +See dpd/dpd.ini for an example. + +TODO +---- + +Fix timestamps and test if frame data is valid. diff --git a/dpd/dpd.ini b/dpd/dpd.ini new file mode 100644 index 0000000..906827b --- /dev/null +++ b/dpd/dpd.ini @@ -0,0 +1,41 @@ +[remotecontrol] +telnet=1 +telnetport=2121 + +[log] +syslog=0 +filelog=0 +filename=/dev/stderr + +[input] +transport=tcp +source=localhost:9200 + +[modulator] +digital_gain=0.9 +rate=8192000 + +[firfilter] +enabled=0 + +[output] +output=uhd + +[uhdoutput] +device= +master_clock_rate=32768000 +type=b200 +txgain=50 +channel=13C +refclk_source=internal +pps_source=none +behaviour_refclk_lock_lost=ignore +max_gps_holdover_time=600 +dpd_port=50055 + +[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 -- cgit v1.2.3 From 403ce1709cd8204769f43a2e6cc68c0286d0fb25 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 14:49:32 +0200 Subject: Fix typo in show_spectrum.py --- dpd/show_spectrum.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index c1d5fe5..477cd99 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -77,7 +77,7 @@ def get_samples(port, num_samps_to_request): rxframe_bytes = recv_exact(s, num_samps * SIZEOF_SAMPLE) rxframe = np.fromstring(rxframe_bytes, dtype=np.complex64) else: - txframe = np.array([], dtype=np.complex64) + rxframe = np.array([], dtype=np.complex64) print("Disconnecting") s.close() -- cgit v1.2.3 From f69b4a1da493c73192aa0c86749bd6bcf396422d Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 15:01:08 +0200 Subject: Actually plot a spectrum in show_spectrum.py --- dpd/show_spectrum.py | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index 477cd99..fcc682f 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -18,7 +18,7 @@ import sys import socket import struct import numpy as np -import matplotlib as mp +import matplotlib.pyplot as pp import scipy.signal SIZEOF_SAMPLE = 8 # complex floats @@ -90,6 +90,27 @@ tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) 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") +fft_size = 4096 +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)) + +sampling_rate = 8192000 +freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate)) + +fig = pp.figure() + +fig.suptitle("TX and RX spectrum") +ax1 = fig.add_subplot(211) +ax1.set_title("TX") +ax1.plot(freqs, tx_power) +ax2 = fig.add_subplot(212) +ax2.set_title("RX") +ax2.plot(freqs, rx_power) +pp.show() # The MIT License (MIT) # -- cgit v1.2.3 From 1deec1418cd49e27bc2f7ddd6cd22ac6607642b3 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 15:21:22 +0200 Subject: DPD: Set RX gain and frequency, add gain to RX and config --- dpd/dpd.ini | 3 ++ src/ConfigParser.cpp | 1 + src/OutputUHD.cpp | 70 +++++++++++++++++++++++++++++++---------------- src/OutputUHD.h | 1 + src/OutputUHDFeedback.cpp | 10 +++---- 5 files changed, 56 insertions(+), 29 deletions(-) (limited to 'dpd') diff --git a/dpd/dpd.ini b/dpd/dpd.ini index 906827b..625df73 100644 --- a/dpd/dpd.ini +++ b/dpd/dpd.ini @@ -1,6 +1,8 @@ [remotecontrol] telnet=1 telnetport=2121 +zmqctrl=1 +zmqctrlendpoint=tcp://127.0.0.1:9400 [log] syslog=0 @@ -32,6 +34,7 @@ pps_source=none behaviour_refclk_lock_lost=ignore max_gps_holdover_time=600 dpd_port=50055 +rxgain=10 [delaymanagement] ; Use synchronous=1 so that the USRP time is set. This works diff --git a/src/ConfigParser.cpp b/src/ConfigParser.cpp index 8892642..459811f 100644 --- a/src/ConfigParser.cpp +++ b/src/ConfigParser.cpp @@ -218,6 +218,7 @@ static void parse_configfile( } outputuhd_conf.txgain = pt.get("uhdoutput.txgain", 0.0); + outputuhd_conf.rxgain = pt.get("uhdoutput.rxgain", 0.0); outputuhd_conf.frequency = pt.get("uhdoutput.frequency", 0); std::string chan = pt.get("uhdoutput.channel", ""); outputuhd_conf.dabMode = mod_settings.dabMode; diff --git a/src/OutputUHD.cpp b/src/OutputUHD.cpp index 5e9e17c..f764fb8 100644 --- a/src/OutputUHD.cpp +++ b/src/OutputUHD.cpp @@ -80,6 +80,36 @@ void uhd_msg_handler(uhd::msg::type_t type, const std::string &msg) } } +static void tune_usrp_to( + uhd::usrp::multi_usrp::sptr usrp, + double lo_offset, + double frequency) +{ + if (lo_offset != 0.0) { + etiLog.level(info) << std::fixed << std::setprecision(3) << + "OutputUHD:Setting freq to " << frequency << + " with LO offset " << lo_offset << "..."; + + const auto tr = uhd::tune_request_t(frequency, lo_offset); + uhd::tune_result_t result = usrp->set_tx_freq(tr); + + etiLog.level(debug) << "OutputUHD:" << + std::fixed << std::setprecision(0) << + " Target RF: " << result.target_rf_freq << + " Actual RF: " << result.actual_rf_freq << + " Target DSP: " << result.target_dsp_freq << + " Actual DSP: " << result.actual_dsp_freq; + } + else { + //set the centre frequency + etiLog.level(info) << std::fixed << std::setprecision(3) << + "OutputUHD:Setting freq to " << frequency << "..."; + usrp->set_tx_freq(frequency); + } + + usrp->set_rx_freq(frequency); +} + // Check function for GPS TIMELOCK sensor from the ODR LEA-M8F board GPSDO bool check_gps_timelock(uhd::usrp::multi_usrp::sptr usrp) { @@ -165,6 +195,7 @@ OutputUHD::OutputUHD( /* register the parameters that can be remote controlled */ RC_ADD_PARAMETER(txgain, "UHD analog daughterboard TX gain"); + RC_ADD_PARAMETER(rxgain, "UHD analog daughterboard RX gain for DPD feedback"); RC_ADD_PARAMETER(freq, "UHD transmission frequency"); RC_ADD_PARAMETER(muting, "Mute the output by stopping the transmitter"); RC_ADD_PARAMETER(staticdelay, "Set static delay (uS) between 0 and 96000"); @@ -223,31 +254,14 @@ OutputUHD::OutputUHD( throw std::runtime_error("Cannot set USRP sample rate. Aborted."); } - if (myConf.lo_offset != 0.0) { - etiLog.level(info) << std::fixed << std::setprecision(3) << - "OutputUHD:Setting freq to " << myConf.frequency << - " with LO offset " << myConf.lo_offset << "..."; - - const auto tr = uhd::tune_request_t(myConf.frequency, myConf.lo_offset); - uhd::tune_result_t result = myUsrp->set_tx_freq(tr); - - etiLog.level(debug) << "OutputUHD:" << - std::fixed << std::setprecision(0) << - " Target RF: " << result.target_rf_freq << - " Actual RF: " << result.actual_rf_freq << - " Target DSP: " << result.target_dsp_freq << - " Actual DSP: " << result.actual_dsp_freq; - } - else { - //set the centre frequency - etiLog.level(info) << std::fixed << std::setprecision(3) << - "OutputUHD:Setting freq to " << myConf.frequency << "..."; - myUsrp->set_tx_freq(myConf.frequency); - } + tune_usrp_to(myUsrp, myConf.lo_offset, myConf.frequency); myConf.frequency = myUsrp->get_tx_freq(); etiLog.level(info) << std::fixed << std::setprecision(3) << - "OutputUHD:Actual frequency: " << myConf.frequency; + "OutputUHD:Actual TX frequency: " << myConf.frequency; + + etiLog.level(info) << std::fixed << std::setprecision(3) << + "OutputUHD:Actual RX frequency: " << myUsrp->get_tx_freq(); myUsrp->set_tx_gain(myConf.txgain); MDEBUG("OutputUHD:Actual TX Gain: %f ...\n", myUsrp->get_tx_gain()); @@ -284,6 +298,9 @@ OutputUHD::OutputUHD( SetDelayBuffer(myConf.dabMode); + myUsrp->set_rx_gain(myConf.rxgain); + MDEBUG("OutputUHD:Actual RX Gain: %f ...\n", myUsrp->get_rx_gain()); + uhdFeedback.setup(myUsrp, myConf.dpdFeedbackServerPort, myConf.sampleRate); MDEBUG("OutputUHD:UHD ready.\n"); @@ -910,9 +927,13 @@ void OutputUHD::set_parameter(const string& parameter, const string& value) ss >> myConf.txgain; myUsrp->set_tx_gain(myConf.txgain); } + else if (parameter == "rxgain") { + ss >> myConf.rxgain; + myUsrp->set_rx_gain(myConf.rxgain); + } else if (parameter == "freq") { ss >> myConf.frequency; - myUsrp->set_tx_freq(myConf.frequency); + tune_usrp_to(myUsrp, myConf.lo_offset, myConf.frequency); myConf.frequency = myUsrp->get_tx_freq(); } else if (parameter == "muting") { @@ -951,6 +972,9 @@ const string OutputUHD::get_parameter(const string& parameter) const if (parameter == "txgain") { ss << myConf.txgain; } + else if (parameter == "rxgain") { + ss << myConf.rxgain; + } else if (parameter == "freq") { ss << myConf.frequency; } diff --git a/src/OutputUHD.h b/src/OutputUHD.h index 1246fc5..c966c7e 100644 --- a/src/OutputUHD.h +++ b/src/OutputUHD.h @@ -189,6 +189,7 @@ struct OutputUHDConfig { double frequency = 0.0; double lo_offset = 0.0; double txgain = 0.0; + double rxgain = 0.0; bool enableSync = false; bool muteNoTimestamps = false; unsigned dabMode = 0; diff --git a/src/OutputUHDFeedback.cpp b/src/OutputUHDFeedback.cpp index 9e3aab2..788b0a9 100644 --- a/src/OutputUHDFeedback.cpp +++ b/src/OutputUHDFeedback.cpp @@ -153,9 +153,6 @@ void OutputUHDFeedback::ReceiveBurstThread() const double usrp_time = m_usrp->get_time_now().get_real_secs(); const double cmd_time = cmd.time_spec.get_real_secs(); - etiLog.level(debug) << - "RX stream command ts=" << std::fixed << cmd_time << " Delta=" << cmd_time - usrp_time; - rxStream->issue_stream_cmd(cmd); uhd::rx_metadata_t md; @@ -173,9 +170,10 @@ void OutputUHDFeedback::ReceiveBurstThread() burstRequest.rx_second = md.time_spec.get_full_secs(); burstRequest.rx_pps = md.time_spec.get_frac_secs() * 16384000.0; - etiLog.level(debug) << "Read " << samples_read << " RX feedback samples " - << "at time " << std::fixed << burstRequest.tx_second << "." << - burstRequest.tx_pps / 16384000.0; + etiLog.level(debug) << "DPD: acquired " << samples_read << " RX feedback samples " << + "at time " << burstRequest.tx_second << " + " << + std::fixed << burstRequest.tx_pps / 16384000.0 << + " Delta=" << cmd_time - usrp_time; burstRequest.state = BurstRequestState::Acquired; -- cgit v1.2.3 From d57fe2c74f2c9e0a76f4b2c577942837dfac0866 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 17:20:29 +0200 Subject: DPD: Set RX sample rate --- dpd/README.md | 6 ++++-- dpd/dpd.ini | 4 ++-- src/OutputUHD.cpp | 3 +++ 3 files changed, 9 insertions(+), 4 deletions(-) (limited to 'dpd') diff --git a/dpd/README.md b/dpd/README.md index 96c4fb0..ec7cec2 100644 --- a/dpd/README.md +++ b/dpd/README.md @@ -8,11 +8,13 @@ This folder contains work in progress for digital predistortion. It requires: - 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. +- 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 ---- -Fix timestamps and test if frame data is valid. +Implement a PA model that updates the predistorter. diff --git a/dpd/dpd.ini b/dpd/dpd.ini index 625df73..910f251 100644 --- a/dpd/dpd.ini +++ b/dpd/dpd.ini @@ -14,7 +14,7 @@ transport=tcp source=localhost:9200 [modulator] -digital_gain=0.9 +digital_gain=0.8 rate=8192000 [firfilter] @@ -34,7 +34,7 @@ pps_source=none behaviour_refclk_lock_lost=ignore max_gps_holdover_time=600 dpd_port=50055 -rxgain=10 +rxgain=0 [delaymanagement] ; Use synchronous=1 so that the USRP time is set. This works diff --git a/src/OutputUHD.cpp b/src/OutputUHD.cpp index 1a137d3..c2f985b 100644 --- a/src/OutputUHD.cpp +++ b/src/OutputUHD.cpp @@ -298,6 +298,9 @@ OutputUHD::OutputUHD( SetDelayBuffer(myConf.dabMode); + myUsrp->set_rx_rate(myConf.sampleRate); + MDEBUG("OutputUHD:Actual RX Rate: %f sps...\n", myUsrp->get_rx_rate()); + myUsrp->set_rx_antenna("RX2"); MDEBUG("OutputUHD:Set RX Antenna: %s ...\n", myUsrp->get_rx_antenna().c_str()); -- cgit v1.2.3 From a34f81d00b9a3f247fa6bd2e31d2c95d4ce7c12f Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Fri, 12 May 2017 17:53:03 +0200 Subject: Add animation support to show_spectrum --- dpd/show_spectrum.py | 99 ++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 73 insertions(+), 26 deletions(-) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index fcc682f..b9e0180 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -19,16 +19,30 @@ import socket import struct import numpy as np import matplotlib.pyplot as pp -import scipy.signal +from matplotlib.animation import FuncAnimation +import argparse SIZEOF_SAMPLE = 8 # complex floats -if len(sys.argv) != 3: - print("Usage: show_spectrum.py ") - sys.exit(1) +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) -port = int(sys.argv[1]) -num_samps_to_request = int(sys.argv[2]) + 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 = [] @@ -84,33 +98,66 @@ def get_samples(port, num_samps_to_request): 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)) -tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) +def plot_spectrum_once(port, num_samps_to_request): + tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) -print("Received {} & {} frames at {} and {}".format( - len(txframe), len(rxframe), tx_ts, rx_ts)) + 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") -fft_size = 4096 -tx_spectrum = np.fft.fftshift(np.fft.fft(txframe, fft_size)) -tx_power = 20*np.log10(np.abs(tx_spectrum)) + 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)) + rx_spectrum = np.fft.fftshift(np.fft.fft(rxframe, fft_size)) + rx_power = 20*np.log10(np.abs(rx_spectrum)) -sampling_rate = 8192000 -freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./sampling_rate)) + 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_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) + + 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)) + + lines[0].set_ydata(tx_power) + lines[1].set_ydata(rx_power) + return lines -fig = pp.figure() + ani = FuncAnimation(fig, update, blit=True) + pp.show() -fig.suptitle("TX and RX spectrum") -ax1 = fig.add_subplot(211) -ax1.set_title("TX") -ax1.plot(freqs, tx_power) -ax2 = fig.add_subplot(212) -ax2.set_title("RX") -ax2.plot(freqs, rx_power) -pp.show() +main() # The MIT License (MIT) # -- cgit v1.2.3 From 7e536b0169cd606dac6c2f4241a520e1c68bfa4a Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Sat, 13 May 2017 11:06:57 +0200 Subject: DPD: Move power calculation into function --- dpd/show_spectrum.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index b9e0180..83ecd8d 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -98,11 +98,7 @@ def get_samples(port, num_samps_to_request): 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)) - -def plot_spectrum_once(port, num_samps_to_request): +def get_spectrum(port, num_samps_to_request): tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) print("Received {} & {} frames at {} and {}".format( @@ -114,7 +110,15 @@ def plot_spectrum_once(port, num_samps_to_request): 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") @@ -138,17 +142,7 @@ def plot_spectrum_animated(port, num_samps_to_request): axes[1].set_ylim(-60, 40) def update(frame): - tx_ts, txframe, rx_ts, rxframe = get_samples(port, num_samps_to_request) - - 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)) + tx_power, rx_power = get_spectrum(port, num_samps_to_request) lines[0].set_ydata(tx_power) lines[1].set_ydata(rx_power) -- cgit v1.2.3 From acec874ecd00b9a1f87b0d1405f6fa1e6503d2b7 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Sat, 13 May 2017 15:27:35 +0200 Subject: DPD: calculate spectrum in complex double --- dpd/show_spectrum.py | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'dpd') diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py index 83ecd8d..e92c1d0 100755 --- a/dpd/show_spectrum.py +++ b/dpd/show_spectrum.py @@ -101,6 +101,10 @@ def get_samples(port, num_samps_to_request): 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)) -- cgit v1.2.3 From bfe0a3fc254bc3117955e64532d8ff7a0de0b1fe Mon Sep 17 00:00:00 2001 From: andreas128 Date: Mon, 29 May 2017 22:25:41 +0100 Subject: Add store_received with alignment --- dpd/store_received.py | 157 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100755 dpd/store_received.py (limited to 'dpd') 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. -- cgit v1.2.3 From 418eb33d0948bb12b7b2ed2179d43ad66258aa72 Mon Sep 17 00:00:00 2001 From: andreas128 Date: Mon, 29 May 2017 22:51:19 +0100 Subject: Configure txgrain --- dpd/dpd.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'dpd') diff --git a/dpd/dpd.ini b/dpd/dpd.ini index 910f251..5e809e5 100644 --- a/dpd/dpd.ini +++ b/dpd/dpd.ini @@ -27,7 +27,7 @@ output=uhd device= master_clock_rate=32768000 type=b200 -txgain=50 +txgain=75 channel=13C refclk_source=internal pps_source=none -- cgit v1.2.3