aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore119
-rw-r--r--INSTALL20
-rw-r--r--README.md2
-rw-r--r--dpd/README.md48
-rw-r--r--dpd/dpd.ini25
-rwxr-xr-xdpd/iq_file_server.py120
-rwxr-xr-xdpd/main.py176
-rw-r--r--dpd/poly.coef11
-rwxr-xr-xdpd/show_spectrum.py132
-rw-r--r--dpd/src/Adapt.py190
-rw-r--r--dpd/src/Agc.py165
-rw-r--r--dpd/src/Dab_Util.py244
-rw-r--r--dpd/src/MER.py134
-rw-r--r--dpd/src/Measure.py145
-rw-r--r--dpd/src/Model.py355
-rw-r--r--dpd/src/Symbol_align.py191
-rw-r--r--dpd/src/TX_Agc.py109
-rw-r--r--dpd/src/Test_data.py136
-rw-r--r--dpd/src/__init__.py0
-rw-r--r--dpd/src/const.py38
-rw-r--r--dpd/src/phase_align.py102
-rwxr-xr-xdpd/src/subsample_align.py112
-rw-r--r--dpd/src/test_dab_Util.py92
-rw-r--r--dpd/src/test_measure.py62
-rw-r--r--src/ConfigParser.cpp5
-rw-r--r--src/ConfigParser.h2
-rw-r--r--src/DabModulator.cpp9
-rw-r--r--src/MemlessPoly.cpp244
-rw-r--r--src/MemlessPoly.h56
-rw-r--r--src/OutputUHDFeedback.cpp287
-rw-r--r--src/OutputUHDFeedback.h1
31 files changed, 3093 insertions, 239 deletions
diff --git a/.gitignore b/.gitignore
index 7fac4b1..d797fb1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,3 +19,122 @@ odr-dabmod
.*.swp
*~
+
+
+# Created by https://www.gitignore.io/api/vim,python
+
+### Python ###
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+# Usually these files are written by a python script from a template
+# before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# SageMath parsed files
+*.sage.py
+
+# dotenv
+.env
+
+# virtualenv
+.venv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+### Vim ###
+# swap
+[._]*.s[a-v][a-z]
+[._]*.sw[a-p]
+[._]s[a-v][a-z]
+[._]sw[a-p]
+# session
+Session.vim
+# temporary
+.netrwhist
+*~
+# auto-generated tag files
+tags
+
+# End of https://www.gitignore.io/api/vim,python
diff --git a/INSTALL b/INSTALL
index 5b67803..a11fd63 100644
--- a/INSTALL
+++ b/INSTALL
@@ -28,6 +28,10 @@ The configure script can be launch with a variety of options:
--disable-output-uhd Disable the binding to the UHD driver for USRPs
--enable-fast-math Compile using the -ffast-math option that gives a substantial
speedup at the cost of correctness.
+ --disable-native Do not compile ODR-DabMod with -march=native compiler option.
+ This is meant for distribution package maintainers who want to
+ use their own march option, and for people running into compilation
+ issues due to -march=native. (e.g. GCC bug 70132 on ARM systems)
Debugging options: You should not enable debug if you need good performance.
By default, debug is disabled.
@@ -38,6 +42,22 @@ For more information, call:
% ./configure --help
+Performance optimisation
+------------------------
+While the performance of modern systems is in most cases good enough to
+run ODR-DabMod, it is sometimes necessary to increase the compilation
+optimisation if all features are used or on slow systems.
+
+Tricks for best performance:
+
+ * Do not use --disable-native
+ * Use --enable-fast-math
+ * Add -O3 to compiler flags
+
+Applying all together:
+
+ % CFLAGS="-O3" CXXFLAGS="-O3" ./configure --enable-fast-math
+
Nearly as simple install procedure using repository:
====================================================
diff --git a/README.md b/README.md
index 2150395..3f72fe6 100644
--- a/README.md
+++ b/README.md
@@ -31,6 +31,8 @@ Short list of features:
- A Telnet and ZeroMQ remote-control that can be used to change
some parameters during runtime
- ZeroMQ PUB and REP output.
+- Ongoing work about digital predistortion for PA linearisation.
+ See dpd/README.md
The src/ directory contains the source code of ODR-DabMod.
diff --git a/dpd/README.md b/dpd/README.md
index ec7cec2..b5a6b81 100644
--- a/dpd/README.md
+++ b/dpd/README.md
@@ -1,20 +1,56 @@
-Digital Predistortion for ODR-DabMod
-====================================
+Digital Predistortion Calculation Engine for ODR-DabMod
+=======================================================
-This folder contains work in progress for digital predistortion. It requires:
+This folder contains work in progress for digital predistortion.
+
+Concept
+-------
+
+ODR-DabMod makes outgoing TX samples and feedback RX samples available for an external tool. This
+external tool can request a buffer of samples for analysis, can calculate coefficients for the
+polynomial predistorter in ODR-DabMod and load the new coefficients using the remote control.
+
+The *dpd/main.py* script is the entry point for the *DPD Calculation Engine* into which these
+features will be implemented. The tool uses modules from the *dpd/src/* folder:
+
+- Sample transfer and time alignment with subsample accuracy is done by *Measure.py*
+- Estimating the effects of the PA using some model and calculation of the updated
+ polynomial coefficients is done in *Model.py*
+- Finally, *Adapt.py* loads them into ODR-DabMod.
+
+These modules themselves use additional helper scripts in the *dpd/src/* folder.
+
+Requirements
+------------
- 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.
+ Usually this is done with a directional coupler and additional attenuators.
+- 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.
+The DPD server port can be tested with the *dpd/show_spectrum.py* helper tool, which can also display
+a constellation diagram.
+
+File format for coefficients
+----------------------------
+The coef file contains the polynomial coefficients used in the predistorter. The file format is
+very similar to the filtertaps file used in the FIR filter. It is a text-based format that can
+easily be inspected and edited in a text editor.
+
+The first line contains the number of coefficients as an integer. The second and third lines contain
+the real, respectively the imaginary parts of the first coefficient. Fourth and fifth lines give the
+second coefficient, and so on. The file therefore contains 2xN + 1 lines if it contains N
+coefficients.
+
TODO
----
-Implement a PA model that updates the predistorter.
+Implement a PA model.
+Implement cases for different oversampling for FFT bin choice.
+Fix loads of missing and buggy aspects of the implementation.
diff --git a/dpd/dpd.ini b/dpd/dpd.ini
index 5e809e5..1bc51de 100644
--- a/dpd/dpd.ini
+++ b/dpd/dpd.ini
@@ -6,35 +6,48 @@ zmqctrlendpoint=tcp://127.0.0.1:9400
[log]
syslog=0
-filelog=0
-filename=/dev/stderr
+filelog=1
+filename=/tmp/dabmod.log
[input]
transport=tcp
source=localhost:9200
[modulator]
-digital_gain=0.8
+gainmode=var
rate=8192000
[firfilter]
-enabled=0
+enabled=1
+
+[poly]
+enabled=1
+polycoeffile=dpd/poly.coef
+
+# How many threads to use for the predistorter.
+# If not set, detect automatically.
+#num_threads=2
[output]
+# to prepare a file for the dpd/iq_file_server.py script,
+# use output=file
output=uhd
+[fileoutput]
+filename=dpd.iq
+
[uhdoutput]
device=
master_clock_rate=32768000
type=b200
-txgain=75
+txgain=15
channel=13C
refclk_source=internal
pps_source=none
behaviour_refclk_lock_lost=ignore
max_gps_holdover_time=600
dpd_port=50055
-rxgain=0
+rxgain=15
[delaymanagement]
; Use synchronous=1 so that the USRP time is set. This works
diff --git a/dpd/iq_file_server.py b/dpd/iq_file_server.py
new file mode 100755
index 0000000..7a4e570
--- /dev/null
+++ b/dpd/iq_file_server.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# This example server simulates the ODR-DabMod's
+# DPD server, taking samples from an IQ file
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import sys
+import socket
+import struct
+import argparse
+import numpy as np
+from datetime import datetime
+
+SIZEOF_SAMPLE = 8 # complex floats
+# Constants for TM 1
+NbSymbols = 76
+NbCarriers = 1536
+Spacing = 2048
+NullSize = 2656
+SymSize = 2552
+FicSizeOut = 288
+FrameSize = NullSize + NbSymbols*SymSize
+
+def main():
+ parser = argparse.ArgumentParser(description="Simulate ODR-DabMod DPD server")
+ parser.add_argument('--port', default='50055',
+ help='port to listen on (default: 50055)',
+ required=False)
+ parser.add_argument('--file', help='I/Q File to read from (complex floats)',
+ required=True)
+ parser.add_argument('--samplerate', default='8192000', help='Sample rate',
+ required=False)
+
+ cli_args = parser.parse_args()
+
+ serve_file(cli_args)
+
+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 serve_file(options):
+ oversampling = int(int(options.samplerate) / 2048000)
+ consumesamples = 8*FrameSize * oversampling
+ iq_data = np.fromfile(options.file, count=consumesamples, dtype=np.complex64)
+
+ print("Loaded {} samples of IQ data".format(len(iq_data)))
+
+ s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+ s.bind(('localhost', int(options.port)))
+ s.listen()
+
+ try:
+ while True:
+ sock, addr_info = s.accept()
+ print("Got a connection from {}".format(addr_info))
+
+ ver = recv_exact(sock, 1)
+ (num_samps,) = struct.unpack("=I", recv_exact(sock, 4))
+ num_bytes = num_samps * SIZEOF_SAMPLE
+
+ if num_bytes > len(iq_data):
+ print("Truncating length to {} samples".format(len(iq_data)))
+ num_samps = len(iq_data)
+ num_bytes = num_samps * 4
+
+ tx_sec = datetime.now().timestamp()
+ tx_pps = int(16384000 * (tx_sec - int(tx_sec)))
+ tx_second = int(tx_sec)
+
+ # TX metadata and data
+ sock.sendall(struct.pack("=III", num_samps, tx_second, tx_pps))
+ sock.sendall(iq_data[-num_samps:].tobytes())
+
+ # RX metadata and data
+ rx_second = tx_second + 1
+ rx_pps = tx_pps
+ sock.sendall(struct.pack("=III", num_samps, rx_second, rx_pps))
+ sock.sendall(iq_data[-num_samps:].tobytes())
+
+ print("Sent {} samples".format(num_samps))
+
+ sock.close()
+ finally:
+ s.close()
+ raise
+
+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/main.py b/dpd/main.py
new file mode 100755
index 0000000..5e67c90
--- /dev/null
+++ b/dpd/main.py
@@ -0,0 +1,176 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine main file.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+"""This Python script is the main file for ODR-DabMod's DPD Computation Engine.
+This engine calculates and updates the parameter of the digital
+predistortion module of ODR-DabMod."""
+
+import datetime
+import os
+import time
+
+import matplotlib
+matplotlib.use('GTKAgg')
+
+import logging
+
+dt = datetime.datetime.now().isoformat()
+logging_path = "/tmp/dpd_{}".format(dt).replace(".", "_").replace(":", "-")
+os.makedirs(logging_path)
+logging.basicConfig(format='%(asctime)s - %(module)s - %(levelname)s - %(message)s',
+ datefmt='%Y-%m-%d %H:%M:%S',
+ filename='{}/dpd.log'.format(logging_path),
+ filemode='w',
+ level=logging.DEBUG)
+
+# also log up to INFO to console
+console = logging.StreamHandler()
+console.setLevel(logging.INFO)
+# set a format which is simpler for console use
+formatter = logging.Formatter('%(asctime)s - %(module)s - %(levelname)s - %(message)s')
+# tell the handler to use this format
+console.setFormatter(formatter)
+# add the handler to the root logger
+logging.getLogger('').addHandler(console)
+
+import numpy as np
+import traceback
+import src.Measure as Measure
+import src.Model as Model
+import src.Adapt as Adapt
+import src.Agc as Agc
+import src.TX_Agc as TX_Agc
+import src.Symbol_align
+import src.const
+import src.MER
+import argparse
+
+parser = argparse.ArgumentParser(
+ description="DPD Computation Engine for ODR-DabMod")
+parser.add_argument('--port', default='50055',
+ help='port of DPD server to connect to (default: 50055)',
+ required=False)
+parser.add_argument('--rc-port', default='9400',
+ help='port of ODR-DabMod ZMQ Remote Control to connect to (default: 9400)',
+ required=False)
+parser.add_argument('--samplerate', default='8192000',
+ help='Sample rate',
+ required=False)
+parser.add_argument('--coefs', default='poly.coef',
+ help='File with DPD coefficients, which will be read by ODR-DabMod',
+ required=False)
+parser.add_argument('--txgain', default=71,
+ help='TX Gain',
+ required=False,
+ type=int)
+parser.add_argument('--rxgain', default=30,
+ help='TX Gain',
+ required=False,
+ type=int)
+parser.add_argument('--digital_gain', default=1,
+ help='Digital Gain',
+ required=False,
+ type=float)
+parser.add_argument('--samps', default='81920',
+ help='Number of samples to request from ODR-DabMod',
+ required=False)
+parser.add_argument('-i', '--iterations', default='1',
+ help='Number of iterations to run',
+ required=False)
+parser.add_argument('-l', '--load-poly',
+ help='Load existing polynomial',
+ action="store_true")
+
+cli_args = parser.parse_args()
+
+port = int(cli_args.port)
+port_rc = int(cli_args.rc_port)
+coef_path = cli_args.coefs
+digital_gain = cli_args.digital_gain
+txgain = cli_args.txgain
+rxgain = cli_args.rxgain
+num_req = int(cli_args.samps)
+samplerate = int(cli_args.samplerate)
+num_iter = int(cli_args.iterations)
+
+SA = src.Symbol_align.Symbol_align(samplerate)
+MER = src.MER.MER(samplerate)
+c = src.const.const(samplerate)
+
+meas = Measure.Measure(samplerate, port, num_req)
+
+adapt = Adapt.Adapt(port_rc, coef_path)
+coefs_am, coefs_pm = adapt.get_coefs()
+if cli_args.load_poly:
+ model = Model.Model(c, SA, MER, coefs_am, coefs_pm, plot=True)
+else:
+ model = Model.Model(c, SA, MER, [1.0, 0, 0, 0, 0], [0, 0, 0, 0, 0], plot=True)
+adapt.set_coefs(model.coefs_am, model.coefs_pm)
+adapt.set_digital_gain(digital_gain)
+adapt.set_txgain(txgain)
+adapt.set_rxgain(rxgain)
+
+tx_gain = adapt.get_txgain()
+rx_gain = adapt.get_rxgain()
+digital_gain = adapt.get_digital_gain()
+dpd_coefs_am, dpd_coefs_pm = adapt.get_coefs()
+logging.info(
+ "TX gain {}, RX gain {}, dpd_coefs_am {},"
+ " dpd_coefs_pm {}, digital_gain {}".format(
+ tx_gain, rx_gain, dpd_coefs_am, dpd_coefs_pm, digital_gain
+ )
+)
+
+tx_agc = TX_Agc.TX_Agc(adapt)
+
+# Automatic Gain Control
+agc = Agc.Agc(meas, adapt)
+agc.run()
+
+for i in range(num_iter):
+ try:
+ txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median = meas.get_samples()
+ logging.debug("tx_ts {}, rx_ts {}".format(tx_ts, rx_ts))
+ assert tx_ts - rx_ts < 1e-5, "Time stamps do not match."
+
+ if tx_agc.adapt_if_necessary(txframe_aligned):
+ continue
+
+ coefs_am, coefs_pm = model.get_next_coefs(txframe_aligned, rxframe_aligned)
+ adapt.set_coefs(coefs_am, coefs_pm)
+
+ off = SA.calc_offset(txframe_aligned)
+ tx_mer = MER.calc_mer(txframe_aligned[off:off + c.T_U])
+ rx_mer = MER.calc_mer(rxframe_aligned[off:off + c.T_U], debug=True)
+ logging.info("MER with lag in it. {}: TX {}, RX {}".
+ format(i, tx_mer, rx_mer))
+ except Exception as e:
+ logging.warning("Iteration {} failed.".format(i))
+ logging.warning(traceback.format_exc())
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger, 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/poly.coef b/dpd/poly.coef
new file mode 100644
index 0000000..f65e369
--- /dev/null
+++ b/dpd/poly.coef
@@ -0,0 +1,11 @@
+5
+1
+0
+0
+0
+0
+0
+0
+0
+0
+0
diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py
index e92c1d0..f23dba2 100755
--- a/dpd/show_spectrum.py
+++ b/dpd/show_spectrum.py
@@ -1,16 +1,15 @@
#!/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.
+# 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.
+# 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
@@ -21,9 +20,18 @@ import numpy as np
import matplotlib.pyplot as pp
from matplotlib.animation import FuncAnimation
import argparse
+from scipy.misc import imsave
SIZEOF_SAMPLE = 8 # complex floats
+# Constants for TM 1
+NbSymbols = 76
+NbCarriers = 1536
+Spacing = 2048
+NullSize = 2656
+SymSize = 2552
+FicSizeOut = 288
+
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',
@@ -31,18 +39,19 @@ def main():
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')
+ parser.add_argument('--constellation', action='store_true', help='Draw constellaton plot')
+ parser.add_argument('--samplerate', default='8192000', help='Sample rate',
+ required=False)
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)
+ if cli_args.constellation:
+ plot_constellation_once(cli_args)
+ elif cli_args.animated:
+ plot_spectrum_animated(cli_args)
else:
- plot_spectrum_once(port, num_samps_to_request)
+ plot_spectrum_once(cli_args)
def recv_exact(sock, num_bytes):
bufs = []
@@ -82,6 +91,7 @@ def get_samples(port, num_samps_to_request):
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
@@ -98,7 +108,7 @@ def get_samples(port, num_samps_to_request):
return (tx_ts, txframe, rx_ts, rxframe)
-def get_spectrum(port, num_samps_to_request):
+def recv_rxtx(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
@@ -107,6 +117,10 @@ def get_spectrum(port, num_samps_to_request):
print("Received {} & {} frames at {} and {}".format(
len(txframe), len(rxframe), tx_ts, rx_ts))
+ return tx_ts, txframe, rx_ts, rxframe
+
+def get_spectrum(port, num_samps_to_request):
+ tx_ts, txframe, rx_ts, rxframe = recv_rxtx(port, num_samps_to_request)
print("Calculate TX and RX spectrum assuming 8192000 samples per second")
tx_spectrum = np.fft.fftshift(np.fft.fft(txframe, fft_size))
@@ -116,12 +130,90 @@ def get_spectrum(port, num_samps_to_request):
rx_power = 20*np.log10(np.abs(rx_spectrum))
return tx_power, rx_power
+def remove_guard_intervals(frame, options):
+ """Remove the cyclic prefix. The frame needs to be aligned to the
+ end of the transmission frame. Transmission Mode 1 is assumed"""
+ oversample = int(int(options.samplerate) / 2048000)
+
+ # From the end, take 2048 samples, then skip 504 samples
+ frame = frame[::-1]
+
+ stride_len = Spacing * oversample
+ stride_advance = SymSize * oversample
+
+ # Truncate the frame to an integer length of strides
+ newlen = len(frame) - (len(frame) % stride_advance)
+ print("Truncating frame from {} to {}".format(len(frame), newlen))
+ frame = frame[:newlen]
+
+ # Remove the cyclic prefix
+ frame = frame.reshape(-1, stride_advance)[:,:stride_len].reshape(-1)
+
+ # Reverse again
+ return frame[::-1]
+
+
+def plot_constellation_once(options):
+ port = int(options.port)
+ num_samps_to_request = int(options.samps)
+
+ tx_ts, txframe, rx_ts, rxframe = recv_rxtx(port, num_samps_to_request)
+
+ frame = remove_guard_intervals(txframe, options)
+
+ oversample = int(int(options.samplerate) / 2048000)
+
+ n = Spacing * oversample # is also number of samples per symbol
+ if len(frame) % n != 0:
+ raise ValueError("Frame length doesn't contain exact number of symbols")
+ num_syms = int(len(frame) / n)
+ print("frame {} has {} symbols".format(len(frame), num_syms))
+ spectrums = np.array([np.fft.fftshift(np.fft.fft(frame[n*i:n*(i+1)], n)) for i in range(num_syms)])
+
+ def normalise(x):
+ """Normalise a real-valued array x to the range [0,1]"""
+ y = x + np.min(x)
+ return x / np.max(x)
+
+ imsave("spectrums.png", np.concatenate([
+ normalise(np.abs(spectrums)),
+ normalise(np.angle(spectrums))]))
+
+ # Only take bins that are supposed to contain energy
+ # i.e. the middle 1536 bins, excluding the bin at n/2
+ assert(n % 2 == 0)
+ n_half = int(n/2)
+ spectrums = np.concatenate(
+ [spectrums[...,n_half-768:n_half],
+ spectrums[...,n_half + 1:n_half + 769]], axis=1)
+
+ sym_indices = (np.tile(np.arange(num_syms-1).reshape(num_syms-1,1), (1,NbCarriers)) +
+ np.tile(np.linspace(-0.4, 0.4, NbCarriers), (num_syms-1, 1) ) )
+ sym_indices = sym_indices.reshape(-1)
+ diff_angles = np.mod(np.diff(np.angle(spectrums, deg=1), axis=0), 360)
+ #sym_points = spectrums[:-1].reshape(-1)
+ # Set amplitude and phase of low power points to zero, avoid cluttering diagram
+ #sym_points[np.abs(sym_points) < np.mean(np.abs(sym_points)) * 0.1] = 0
+
+ print("ix {} spec {} da {}".format(
+ sym_indices.shape, spectrums.shape, diff_angles.shape))
+
+ fig = pp.figure()
+
+ fig.suptitle("Constellation")
+ ax1 = fig.add_subplot(111)
+ ax1.set_title("TX")
+ ax1.scatter(sym_indices, diff_angles.reshape(-1), alpha=0.1)
+
+ pp.show()
-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 plot_spectrum_once(options):
+ port = int(options.port)
+ num_samps_to_request = int(options.samps)
+ freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./int(options.samplerate)))
+
tx_power, rx_power = get_spectrum(port, num_samps_to_request)
fig = pp.figure()
@@ -134,7 +226,11 @@ def plot_spectrum_once(port, num_samps_to_request):
ax2.plot(freqs, rx_power, 'b')
pp.show()
-def plot_spectrum_animated(port, num_samps_to_request):
+def plot_spectrum_animated(options):
+ port = int(options.port)
+ num_samps_to_request = int(options.samps)
+ freqs = np.fft.fftshift(np.fft.fftfreq(fft_size, d=1./int(options.samplerate)))
+
fig, axes = pp.subplots(2, sharex=True)
line1, = axes[0].plot(freqs, np.ones(len(freqs)), 'r', animated=True)
axes[0].set_title("TX")
diff --git a/dpd/src/Adapt.py b/dpd/src/Adapt.py
new file mode 100644
index 0000000..b4042d6
--- /dev/null
+++ b/dpd/src/Adapt.py
@@ -0,0 +1,190 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine: updates ODR-DabMod's predistortion block.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+"""
+This module is used to change settings of ODR-DabMod using
+the ZMQ remote control socket.
+"""
+
+import zmq
+import logging
+import numpy as np
+
+class Adapt:
+ """Uses the ZMQ remote control to change parameters of the DabMod
+
+ Parameters
+ ----------
+ port : int
+ Port at which the ODR-DabMod is listening to connect the
+ ZMQ remote control.
+ """
+
+ def __init__(self, port, coef_path):
+ logging.debug("Instantiate Adapt object")
+ self.port = port
+ self.coef_path = coef_path
+ self.host = "localhost"
+ self._context = zmq.Context()
+
+ def _connect(self):
+ """Establish the connection to ODR-DabMod using
+ a ZMQ socket that is in request mode (Client).
+ Returns a socket"""
+ sock = self._context.socket(zmq.REQ)
+ sock.connect("tcp://%s:%d" % (self.host, self.port))
+
+ sock.send(b"ping")
+ data = [el.decode() for el in sock.recv_multipart()]
+
+ if data != ['ok']:
+ raise RuntimeError(
+ "Could not ping server at %s %d: %s" %
+ (self.host, self.port, data))
+
+ return sock
+
+ def send_receive(self, message):
+ """Send a message to ODR-DabMod. It always
+ returns the answer ODR-DabMod sends back.
+
+ An example message could be
+ "get uhd txgain" or "set uhd txgain 50"
+
+ Parameter
+ ---------
+ message : str
+ The message string that will be sent to the receiver.
+ """
+ sock = self._connect()
+ logging.debug("Send message: %s" % message)
+ msg_parts = message.split(" ")
+ for i, part in enumerate(msg_parts):
+ if i == len(msg_parts) - 1:
+ f = 0
+ else:
+ f = zmq.SNDMORE
+
+ sock.send(part.encode(), flags=f)
+
+ data = [el.decode() for el in sock.recv_multipart()]
+ logging.debug("Received message: %s" % message)
+ return data
+
+ def set_txgain(self, gain):
+ """Set a new txgain for the ODR-DabMod.
+
+ Parameters
+ ----------
+ gain : int
+ new TX gain, in the same format as ODR-DabMod's config file
+ """
+ # TODO this is specific to the B200
+ if gain < 0 or gain > 89:
+ raise ValueError("Gain has to be in [0,89]")
+ return self.send_receive("set uhd txgain %.4f" % float(gain))
+
+ def get_txgain(self):
+ """Get the txgain value in dB for the ODR-DabMod."""
+ # TODO handle failure
+ return float(self.send_receive("get uhd txgain")[0])
+
+ def set_rxgain(self, gain):
+ """Set a new rxgain for the ODR-DabMod.
+
+ Parameters
+ ----------
+ gain : int
+ new RX gain, in the same format as ODR-DabMod's config file
+ """
+ # TODO this is specific to the B200
+ if gain < 0 or gain > 89:
+ raise ValueError("Gain has to be in [0,89]")
+ return self.send_receive("set uhd rxgain %.4f" % float(gain))
+
+ def get_rxgain(self):
+ """Get the rxgain value in dB for the ODR-DabMod."""
+ # TODO handle failure
+ return float(self.send_receive("get uhd rxgain")[0])
+
+ def set_digital_gain(self, gain):
+ """Set a new rxgain for the ODR-DabMod.
+
+ Parameters
+ ----------
+ gain : int
+ new RX gain, in the same format as ODR-DabMod's config file
+ """
+ msg = "set gain digital %.5f" % gain
+ return self.send_receive(msg)
+
+ def get_digital_gain(self):
+ """Get the rxgain value in dB for the ODR-DabMod."""
+ # TODO handle failure
+ return float(self.send_receive("get gain digital")[0])
+
+ def _read_coef_file(self, path):
+ """Load the coefficients from the file in the format given in the README,
+ return ([AM coef], [PM coef])"""
+ coefs_am_out = []
+ coefs_pm_out = []
+ f = open(path, 'r')
+ lines = f.readlines()
+ n_coefs = int(lines[0])
+ coefs = [float(l) for l in lines[1:]]
+ i = 0
+ for c in coefs:
+ if i < n_coefs:
+ coefs_am_out.append(c)
+ elif i < 2*n_coefs:
+ coefs_pm_out.append(c)
+ else:
+ raise ValueError(
+ "Incorrect coef file format: too many coefficients in {}, should be {}, coefs are {}"
+ .format(path, n_coefs, coefs))
+ i += 1
+ f.close()
+ return (coefs_am_out, coefs_pm_out)
+
+ def get_coefs(self):
+ return self._read_coef_file(self.coef_path)
+
+ def _write_coef_file(self, coefs_am, coefs_pm, path):
+ assert(len(coefs_am) == len(coefs_pm))
+
+ f = open(path, 'w')
+ f.write("{}\n".format(len(coefs_am)))
+ for coef in coefs_am:
+ f.write("{}\n".format(coef))
+ for coef in coefs_pm:
+ f.write("{}\n".format(coef))
+ f.close()
+
+ def set_coefs(self, coefs_am, coefs_pm):
+ self._write_coef_file(coefs_am, coefs_pm, self.coef_path)
+ self.send_receive("set memlesspoly coeffile {}".format(self.coef_path))
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger, 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/src/Agc.py b/dpd/src/Agc.py
new file mode 100644
index 0000000..978b607
--- /dev/null
+++ b/dpd/src/Agc.py
@@ -0,0 +1,165 @@
+# -*- coding: utf-8 -*-
+#
+# Automatic Gain Control
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+import time
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
+import src.Adapt as Adapt
+import src.Measure as Measure
+
+class Agc:
+ """The goal of the automatic gain control is to set the
+ RX gain to a value at which all received amplitudes can
+ be detected. This means that the maximum possible amplitude
+ should be quantized at the highest possible digital value.
+
+ A problem we have to face, is that the estimation of the
+ maximum amplitude by applying the max() function is very
+ unstable. This is due to the maximum’s rareness. Therefore
+ we estimate a far more robust value, such as the median,
+ and then approximate the maximum amplitude from it.
+
+ Given this, we tune the RX gain in such a way, that the
+ received signal fulfills our desired property, of having
+ all amplitudes properly quantized."""
+
+ def __init__(self, measure, adapt, min_rxgain=25, peak_to_median=20):
+ assert isinstance(measure, Measure.Measure)
+ assert isinstance(adapt, Adapt.Adapt)
+ self.measure = measure
+ self.adapt = adapt
+ self.min_rxgain = min_rxgain
+ self.rxgain = self.min_rxgain
+ self.peak_to_median = peak_to_median
+
+ def run(self):
+ self.adapt.set_rxgain(self.rxgain)
+
+ for i in range(2):
+ # Measure
+ txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median= \
+ self.measure.get_samples()
+
+ # Estimate Maximum
+ rx_peak = self.peak_to_median * rx_median
+ correction_factor = 20*np.log10(1/rx_peak)
+ self.rxgain = self.rxgain + correction_factor
+
+ assert self.min_rxgain <= self.rxgain, ("Desired RX Gain is {} which is smaller than the minimum of {}".format(
+ self.rxgain, self.min_rxgain))
+
+ logging.info("RX Median {:1.4f}, estimated peak {:1.4f}, correction factor {:1.4f}, new RX gain {:1.4f}".format(
+ rx_median, rx_peak, correction_factor, self.rxgain
+ ))
+
+ self.adapt.set_rxgain(self.rxgain)
+ time.sleep(0.5)
+
+ def plot_estimates(self):
+ """Plots the estimate of for Max, Median, Mean for different
+ number of samples."""
+ self.adapt.set_rxgain(self.min_rxgain)
+ time.sleep(1)
+
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_agc.svg"
+ fig, axs = plt.subplots(2, 2, figsize=(3*6,1*6))
+ axs = axs.ravel()
+
+ for j in range(5):
+ txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median =\
+ self.measure.get_samples()
+
+ rxframe_aligned_abs = np.abs(rxframe_aligned)
+
+ x = np.arange(100, rxframe_aligned_abs.shape[0], dtype=int)
+ rx_max_until = []
+ rx_median_until = []
+ rx_mean_until = []
+ for i in x:
+ rx_max_until.append(np.max(rxframe_aligned_abs[:i]))
+ rx_median_until.append(np.median(rxframe_aligned_abs[:i]))
+ rx_mean_until.append(np.mean(rxframe_aligned_abs[:i]))
+
+ axs[0].plot(x,
+ rx_max_until,
+ label="Run {}".format(j+1),
+ color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.8,0.7)),
+ linestyle="-", linewidth=0.25)
+ axs[0].set_xlabel("Samples")
+ axs[0].set_ylabel("Amplitude")
+ axs[0].set_title("Estimation for Maximum RX Amplitude")
+ axs[0].legend()
+
+ axs[1].plot(x,
+ rx_median_until,
+ label="Run {}".format(j+1),
+ color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)),
+ linestyle="-", linewidth=0.25)
+ axs[1].set_xlabel("Samples")
+ axs[1].set_ylabel("Amplitude")
+ axs[1].set_title("Estimation for Median RX Amplitude")
+ axs[1].legend()
+ ylim_1 = axs[1].get_ylim()
+
+ axs[2].plot(x,
+ rx_mean_until,
+ label="Run {}".format(j+1),
+ color=matplotlib.colors.hsv_to_rgb((1./(j+1.),0.9,0.7)),
+ linestyle="-", linewidth=0.25)
+ axs[2].set_xlabel("Samples")
+ axs[2].set_ylabel("Amplitude")
+ axs[2].set_title("Estimation for Mean RX Amplitude")
+ ylim_2 = axs[2].get_ylim()
+ axs[2].legend()
+
+ axs[1].set_ylim(min(ylim_1[0], ylim_2[0]),
+ max(ylim_1[1], ylim_2[1]))
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+
+ axs[3].hist(rxframe_aligned_abs, bins=60)
+ axs[3].set_xlabel("Amplitude")
+ axs[3].set_ylabel("Frequency")
+ axs[3].set_title("Histogram of Amplitudes")
+ axs[3].legend()
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/Dab_Util.py b/dpd/src/Dab_Util.py
new file mode 100644
index 0000000..e3dbfe3
--- /dev/null
+++ b/dpd/src/Dab_Util.py
@@ -0,0 +1,244 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, utilities for working with DAB signals.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+import scipy
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+import src.subsample_align as sa
+import src.phase_align as pa
+from scipy import signal
+
+class Dab_Util:
+ """Collection of methods that can be applied to an array
+ complex IQ samples of a DAB signal
+ """
+ def __init__(self, sample_rate, plot=False):
+ """
+ :param sample_rate: sample rate [sample/sec] to use for calculations
+ """
+ self.sample_rate = sample_rate
+ self.dab_bandwidth = 1536000 #Bandwidth of a dab signal
+ self.frame_ms = 96 #Duration of a Dab frame
+
+ self.plot=plot
+
+ def lag(self, sig_orig, sig_rec):
+ """
+ Find lag between two signals
+ Args:
+ sig_orig: The signal that has been sent
+ sig_rec: The signal that has been recored
+ """
+ off = sig_rec.shape[0]
+ c = np.abs(signal.correlate(sig_orig, sig_rec))
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ corr_path = (logging_path + "/" + dt + "_tx_rx_corr.svg")
+ plt.plot(c, label="corr")
+ plt.legend()
+ plt.savefig(corr_path)
+ plt.clf()
+
+ return np.argmax(c) - off + 1
+
+ def lag_upsampling(self, sig_orig, sig_rec, n_up):
+ if n_up != 1:
+ sig_orig_up = signal.resample(sig_orig, sig_orig.shape[0] * n_up)
+ sig_rec_up = signal.resample(sig_rec, sig_rec.shape[0] * n_up)
+ else:
+ sig_orig_up = sig_orig
+ sig_rec_up = sig_rec
+ l = self.lag(sig_orig_up, sig_rec_up)
+ l_orig = float(l) / n_up
+ return l_orig
+
+ def subsample_align_upsampling(self, sig_tx, sig_rx, n_up=32):
+ """
+ Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment
+ Using upsampling
+ """
+ assert(sig_tx.shape[0] == sig_rx.shape[0])
+
+ if sig_tx.shape[0] % 2 == 1:
+ sig_tx = sig_tx[:-1]
+ sig_rx = sig_rx[:-1]
+
+ sig1_up = signal.resample(sig_tx, sig_tx.shape[0] * n_up)
+ sig2_up = signal.resample(sig_rx, sig_rx.shape[0] * n_up)
+
+ off_meas = self.lag_upsampling(sig2_up, sig1_up, n_up=1)
+ off = int(abs(off_meas))
+
+ if off_meas > 0:
+ sig1_up = sig1_up[:-off]
+ sig2_up = sig2_up[off:]
+ elif off_meas < 0:
+ sig1_up = sig1_up[off:]
+ sig2_up = sig2_up[:-off]
+
+ sig_tx = signal.resample(sig1_up, sig1_up.shape[0] / n_up).astype(np.complex64)
+ sig_rx = signal.resample(sig2_up, sig2_up.shape[0] / n_up).astype(np.complex64)
+ return sig_tx, sig_rx
+
+ def subsample_align(self, sig_tx, sig_rx):
+ """
+ Returns an aligned version of sig_tx and sig_rx by cropping and subsample alignment
+ """
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_sync_raw.svg"
+
+ fig, axs = plt.subplots(2)
+ axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame")
+ axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame")
+ axs[0].set_title("Raw Data")
+ axs[0].set_ylabel("Amplitude")
+ axs[0].set_xlabel("Samples")
+ axs[0].legend(loc=4)
+
+ axs[1].plot(np.real(sig_tx[:128]), label="TX Frame")
+ axs[1].plot(np.real(sig_rx[:128]), label="RX Frame")
+ axs[1].set_title("Raw Data")
+ axs[1].set_ylabel("Real Part")
+ axs[1].set_xlabel("Samples")
+ axs[1].legend(loc=4)
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+ off_meas = self.lag_upsampling(sig_rx, sig_tx, n_up=1)
+ off = int(abs(off_meas))
+
+ logging.debug("sig_tx_orig: {} {}, sig_rx_orig: {} {}, offset {}".format(
+ len(sig_tx),
+ sig_tx.dtype,
+ len(sig_rx),
+ sig_rx.dtype,
+ 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]
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_sync_sample_aligned.svg"
+
+ fig, axs = plt.subplots(2)
+ axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame")
+ axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame")
+ axs[0].set_title("Sample Aligned Data")
+ axs[0].set_ylabel("Amplitude")
+ axs[0].set_xlabel("Samples")
+ axs[0].legend(loc=4)
+
+ axs[1].plot(np.real(sig_tx[:128]), label="TX Frame")
+ axs[1].plot(np.real(sig_rx[:128]), label="RX Frame")
+ axs[1].set_ylabel("Real Part")
+ axs[1].set_xlabel("Samples")
+ axs[1].legend(loc=4)
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+ sig_rx = sa.subsample_align(sig_rx, sig_tx)
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_sync_subsample_aligned.svg"
+
+ fig, axs = plt.subplots(2)
+ axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame")
+ axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame")
+ axs[0].set_title("Subsample Aligned")
+ axs[0].set_ylabel("Amplitude")
+ axs[0].set_xlabel("Samples")
+ axs[0].legend(loc=4)
+
+ axs[1].plot(np.real(sig_tx[:128]), label="TX Frame")
+ axs[1].plot(np.real(sig_rx[:128]), label="RX Frame")
+ axs[1].set_ylabel("Real Part")
+ axs[1].set_xlabel("Samples")
+ axs[1].legend(loc=4)
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+ sig_rx = pa.phase_align(sig_rx, sig_tx)
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_sync_phase_aligned.svg"
+
+ fig, axs = plt.subplots(2)
+ axs[0].plot(np.abs(sig_tx[:128]), label="TX Frame")
+ axs[0].plot(np.abs(sig_rx[:128]), label="RX Frame")
+ axs[0].set_title("Phase Aligned")
+ axs[0].set_ylabel("Amplitude")
+ axs[0].set_xlabel("Samples")
+ axs[0].legend(loc=4)
+
+ axs[1].plot(np.real(sig_tx[:128]), label="TX Frame")
+ axs[1].plot(np.real(sig_rx[:128]), label="RX Frame")
+ axs[1].set_ylabel("Real Part")
+ axs[1].set_xlabel("Samples")
+ axs[1].legend(loc=4)
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+ logging.debug("Sig1_cut: %d %s, Sig2_cut: %d %s, off: %d" % (len(sig_tx), sig_tx.dtype, len(sig_rx), sig_rx.dtype, off))
+ return sig_tx, sig_rx
+
+ def fromfile(self, filename, offset=0, length=None):
+ if length is None:
+ return np.memmap(filename, dtype=np.complex64, mode='r', offset=64/8*offset)
+ else:
+ return np.memmap(filename, dtype=np.complex64, mode='r', offset=64/8*offset, shape=length)
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/MER.py b/dpd/src/MER.py
new file mode 100644
index 0000000..00fcc23
--- /dev/null
+++ b/dpd/src/MER.py
@@ -0,0 +1,134 @@
+# -*- coding: utf-8 -*-
+#
+# Modulation Error Rate
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+try:
+ logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+except:
+ logging_path = "/tmp/"
+
+import src.const
+import numpy as np
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
+class MER:
+ def __init__(self, sample_rate):
+ self.c = src.const.const(sample_rate)
+
+ def _calc_spectrum(self, tx):
+ fft = np.fft.fftshift(np.fft.fft(tx))
+ return np.delete(fft[self.c.FFT_start:self.c.FFT_end],
+ self.c.FFT_delete)
+
+ def _split_in_carrier(self, x, y):
+ if 0.5 < np.mean((np.abs(np.abs(x) - np.abs(y)) /
+ np.abs(np.abs(x) + np.abs(y)))):
+ # Constellation points are axis aligned on the Im/Re plane
+ x1 = x[(y < x) & (y > -x)]
+ y1 = y[(y < x) & (y > -x)]
+
+ x2 = x[(y > x) & (y > -x)]
+ y2 = y[(y > x) & (y > -x)]
+
+ x3 = x[(y > x) & (y < -x)]
+ y3 = y[(y > x) & (y < -x)]
+
+ x4 = x[(y < x) & (y < -x)]
+ y4 = y[(y < x) & (y < -x)]
+ else:
+ # Constellation points are on the diagonal or Im/Re plane
+ x1 = x[(+x > 0) & (+y > 0)]
+ y1 = y[(+x > 0) & (+y > 0)]
+
+ x2 = x[(-x > 0) & (+y > 0)]
+ y2 = y[(-x > 0) & (+y > 0)]
+
+ x3 = x[(-x > 0) & (-y > 0)]
+ y3 = y[(-x > 0) & (-y > 0)]
+
+ x4 = x[(+x > 0) & (-y > 0)]
+ y4 = y[(+x > 0) & (-y > 0)]
+ return (x1, y1), (x2, y2), (x3, y3), (x4, y4)
+
+ def _calc_mer_for_isolated_constellation_point(self, x, y):
+ """Calculate MER for one constellation point"""
+ x_mean = np.mean(x)
+ y_mean = np.mean(y)
+
+ U_RMS = np.sqrt(x_mean ** 2 + y_mean ** 2)
+ U_ERR = np.mean(np.sqrt((x - x_mean) ** 2 + (y - y_mean) ** 2))
+ MER = 20 * np.log10(U_ERR / U_RMS)
+
+ return x_mean, y_mean, U_RMS, U_ERR, MER
+
+ def calc_mer(self, tx, debug=False):
+ assert tx.shape[0] == self.c.T_U,\
+ "Wrong input length"
+
+ spectrum = self._calc_spectrum(tx)
+
+ if debug:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_MER.svg"
+
+ MERs = []
+ for i, (x, y) in enumerate(self._split_in_carrier(
+ np.real(spectrum),
+ np.imag(spectrum))):
+ x_mean, y_mean, U_RMS, U_ERR, MER =\
+ self._calc_mer_for_isolated_constellation_point(x, y)
+ MERs.append(MER)
+
+ tau = np.linspace(0, 2 * np.pi, num=100)
+ x_err = U_ERR * np.sin(tau) + x_mean
+ y_err = U_ERR * np.cos(tau) + y_mean
+
+ if debug:
+ ax = plt.subplot(221 + i)
+ ax.scatter(x, y, s=0.2, color='black')
+ ax.plot(x_mean, y_mean, 'p', color='red')
+ ax.plot(x_err, y_err, linewidth=2, color='blue')
+ ax.text(0.1, 0.1, "MER {:.1f}dB".format(MER), transform=ax.transAxes)
+ ax.set_xlabel("Real")
+ ax.set_ylabel("Imag")
+ ylim = ax.get_ylim()
+ ax.set_ylim(ylim[0] - (ylim[1] - ylim[0]) * 0.1, ylim[1])
+
+ if debug:
+ plt.tight_layout()
+ plt.savefig(fig_path)
+ plt.show()
+ plt.clf()
+
+ MER_res = 20 * np.log10(np.mean([10 ** (MER / 20) for MER in MERs]))
+ return MER_res
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/Measure.py b/dpd/src/Measure.py
new file mode 100644
index 0000000..2450b8a
--- /dev/null
+++ b/dpd/src/Measure.py
@@ -0,0 +1,145 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, Measure signal using ODR-DabMod's
+# DPD Server.
+#
+# 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 logging
+import src.Dab_Util as DU
+import os
+import datetime
+
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+class Measure:
+ """Collect Measurement from DabMod"""
+ def __init__(self, samplerate, port, num_samples_to_request):
+ logging.info("Instantiate Measure object")
+ self.samplerate = samplerate
+ self.sizeof_sample = 8 # complex floats
+ self.port = port
+ self.num_samples_to_request = num_samples_to_request
+
+ 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.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 %f, max %f, median %f" %
+ (np.min(np.abs(txframe)),
+ np.max(np.abs(txframe)),
+ np.median(np.abs(txframe))))
+
+ logging.debug("rxframe: min %f, max %f, median %f" %
+ (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
+ (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
+ """
+
+ txframe, tx_ts, rxframe, rx_ts = self.receive_tcp()
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG:
+ dt = datetime.datetime.now().isoformat()
+ txframe.tofile(logging_path + "/txframe_" + dt + ".iq")
+
+ # Normalize received signal with sent signal
+ rx_median = np.median(np.abs(rxframe))
+ rxframe = rxframe / rx_median * np.median(np.abs(txframe))
+
+ du = DU.Dab_Util(self.samplerate)
+ txframe_aligned, rxframe_aligned = du.subsample_align(txframe, rxframe)
+
+ logging.info(
+ "Measurement done, tx %d %s, rx %d %s, tx aligned %d %s, rx aligned %d %s"
+ % (len(txframe), txframe.dtype, len(rxframe), rxframe.dtype,
+ len(txframe_aligned), txframe_aligned.dtype, len(rxframe_aligned), rxframe_aligned.dtype) )
+
+ return txframe_aligned, tx_ts, rxframe_aligned, rx_ts, rx_median
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/Model.py b/dpd/src/Model.py
new file mode 100644
index 0000000..827027a
--- /dev/null
+++ b/dpd/src/Model.py
@@ -0,0 +1,355 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, model implementation.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+import matplotlib.pyplot as plt
+from sklearn import linear_model
+
+class Model:
+ """Calculates new coefficients using the measurement and the old
+ coefficients"""
+
+ def __init__(self,
+ c,
+ SA,
+ MER,
+ coefs_am,
+ coefs_pm,
+ learning_rate_am=1.,
+ learning_rate_pm=1.,
+ plot=False):
+ self.c = c
+ self.SA = SA
+ self.MER = MER
+
+ self.learning_rate_am = learning_rate_am
+ self.learning_rate_pm = learning_rate_pm
+
+ self.coefs_am = coefs_am
+ self.coefs_am_history = [coefs_am, ]
+ self.mses_am = []
+ self.errs_am = []
+
+ self.tx_mers = []
+ self.rx_mers = []
+
+ self.coefs_pm = coefs_pm
+ self.coefs_pm_history = [coefs_pm, ]
+ self.errs_pm = []
+
+ self.plot = plot
+
+ def sample_uniformly(self, tx_dpd, rx_received, n_bins=5):
+ """This function returns tx and rx samples in a way
+ that the tx amplitudes have an approximate uniform
+ distribution with respect to the tx_dpd amplitudes"""
+ mask = np.logical_and((np.abs(tx_dpd) > 0.01), (np.abs(rx_received) > 0.01))
+ tx_dpd = tx_dpd[mask]
+ rx_received = rx_received[mask]
+
+ txframe_aligned_abs = np.abs(tx_dpd)
+ ccdf_min = 0
+ ccdf_max = np.max(txframe_aligned_abs)
+ tx_hist, ccdf_edges = np.histogram(txframe_aligned_abs,
+ bins=n_bins,
+ range=(ccdf_min, ccdf_max))
+ n_choise = np.min(tx_hist)
+ tx_choice = np.zeros(n_choise * n_bins, dtype=np.complex64)
+ rx_choice = np.zeros(n_choise * n_bins, dtype=np.complex64)
+
+ for idx, bin in enumerate(tx_hist):
+ indices = np.where((txframe_aligned_abs >= ccdf_edges[idx]) &
+ (txframe_aligned_abs <= ccdf_edges[idx + 1]))[0]
+ indices_choise = np.random.choice(indices,
+ n_choise,
+ replace=False)
+ rx_choice[idx * n_choise:(idx + 1) * n_choise] = \
+ rx_received[indices_choise]
+ tx_choice[idx * n_choise:(idx + 1) * n_choise] = \
+ tx_dpd[indices_choise]
+
+ assert isinstance(rx_choice[0], np.complex64), \
+ "rx_choice is not complex64 but {}".format(rx_choice[0].dtype)
+ assert isinstance(tx_choice[0], np.complex64), \
+ "tx_choice is not complex64 but {}".format(tx_choice[0].dtype)
+
+ return tx_choice, rx_choice
+
+ def dpd_amplitude(self, sig, coefs=None):
+ if coefs is None:
+ coefs = self.coefs_am
+ assert isinstance(sig[0], np.complex64), "Sig is not complex64 but {}".format(sig[0].dtype)
+ sig_abs = np.abs(sig)
+ A_sig = np.vstack([np.ones(sig_abs.shape),
+ sig_abs ** 1,
+ sig_abs ** 2,
+ sig_abs ** 3,
+ sig_abs ** 4,
+ ]).T
+ sig_dpd = sig * np.sum(A_sig * coefs, axis=1)
+ return sig_dpd, A_sig
+
+ def dpd_phase(self, sig, coefs=None):
+ if coefs is None:
+ coefs = self.coefs_pm
+ assert isinstance(sig[0], np.complex64), "Sig is not complex64 but {}".format(sig[0].dtype)
+ sig_abs = np.abs(sig)
+ A_phase = np.vstack([np.ones(sig_abs.shape),
+ sig_abs ** 1,
+ sig_abs ** 2,
+ sig_abs ** 3,
+ sig_abs ** 4,
+ ]).T
+ phase_diff_est = np.sum(A_phase * coefs, axis=1)
+ return phase_diff_est, A_phase
+
+ def _next_am_coefficent(self, tx_choice, rx_choice):
+ """Calculate new coefficients for AM/AM correction"""
+ rx_dpd, rx_A = self.dpd_amplitude(rx_choice)
+ rx_dpd = rx_dpd * (
+ np.median(np.abs(tx_choice)) /
+ np.median(np.abs(rx_dpd)))
+ err = np.abs(rx_dpd) - np.abs(tx_choice)
+ mse = np.mean(np.abs((rx_dpd - tx_choice) ** 2))
+ self.mses_am.append(mse)
+ self.errs_am.append(np.mean(err**2))
+
+ reg = linear_model.Ridge(alpha=0.00001)
+ reg.fit(rx_A, err)
+ a_delta = reg.coef_
+ new_coefs_am = self.coefs_am - self.learning_rate_am * a_delta
+ new_coefs_am = new_coefs_am * (self.coefs_am[0] / new_coefs_am[0])
+ return new_coefs_am
+
+ def _next_pm_coefficent(self, tx_choice, rx_choice):
+ """Calculate new coefficients for AM/PM correction
+ Assuming deviations smaller than pi/2"""
+ phase_diff_choice = np.angle(
+ (rx_choice * tx_choice.conjugate()) /
+ (np.abs(rx_choice) * np.abs(tx_choice))
+ )
+ plt.hist(phase_diff_choice)
+ plt.savefig('/tmp/hist_' + str(np.random.randint(0,1000)) + '.svg')
+ plt.clf()
+ phase_diff_est, phase_A = self.dpd_phase(rx_choice)
+ err_phase = phase_diff_est - phase_diff_choice
+ self.errs_pm.append(np.mean(np.abs(err_phase ** 2)))
+
+ reg = linear_model.Ridge(alpha=0.00001)
+ reg.fit(phase_A, err_phase)
+ p_delta = reg.coef_
+ new_coefs_pm = self.coefs_pm - self.learning_rate_pm * p_delta
+
+ return new_coefs_pm, phase_diff_choice
+
+ def get_next_coefs(self, tx_dpd, rx_received):
+ # Check data type
+ assert tx_dpd[0].dtype == np.complex64, \
+ "tx_dpd is not complex64 but {}".format(tx_dpd[0].dtype)
+ assert rx_received[0].dtype == np.complex64, \
+ "rx_received is not complex64 but {}".format(rx_received[0].dtype)
+ # Check if signals have same normalization
+ normalization_error = np.abs(np.median(np.abs(tx_dpd)) -
+ np.median(np.abs(rx_received))) / (
+ np.median(np.abs(tx_dpd)) + np.median(np.abs(rx_received)))
+ assert normalization_error < 0.01, "Non normalized signals"
+
+ tx_choice, rx_choice = self.sample_uniformly(tx_dpd, rx_received)
+ new_coefs_am = self._next_am_coefficent(tx_choice, rx_choice)
+ new_coefs_pm, phase_diff_choice = self._next_pm_coefficent(tx_choice, rx_choice)
+
+ logging.debug('txframe: min {:.2f}, max {:.2f}, ' \
+ 'median {:.2f}; rxframe: min {:.2f}, max {:.2f}, ' \
+ 'median {:.2f}; new coefs_am {};' \
+ 'new_coefs_pm {}'.format(
+ np.min(np.abs(tx_dpd)),
+ np.max(np.abs(tx_dpd)),
+ np.median(np.abs(tx_dpd)),
+ np.min(np.abs(rx_choice)),
+ np.max(np.abs(rx_choice)),
+ np.median(np.abs(rx_choice)),
+ new_coefs_am,
+ new_coefs_pm))
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ off = self.SA.calc_offset(tx_dpd)
+ tx_mer = self.MER.calc_mer(tx_dpd[off:off + self.c.T_U])
+ rx_mer = self.MER.calc_mer(rx_received[off:off + self.c.T_U], debug=True)
+ self.tx_mers.append(tx_mer)
+ self.rx_mers.append(rx_mer)
+
+ if logging.getLogger().getEffectiveLevel() == logging.DEBUG and self.plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_Model.svg"
+
+ fig = plt.figure(figsize=(2 * 6, 2 * 6))
+
+ i_sub = 1
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ ax.plot(np.abs(tx_dpd[:128]),
+ label="TX sent",
+ linestyle=":")
+ ax.plot(np.abs(rx_received[:128]),
+ label="RX received",
+ color="red")
+ ax.set_title("Synchronized Signals of Iteration {}"
+ .format(len(self.coefs_am_history)))
+ ax.set_xlabel("Samples")
+ ax.set_ylabel("Amplitude")
+ ax.text(0, 0, "TX (max {:01.3f}, mean {:01.3f}, " \
+ "median {:01.3f})".format(
+ np.max(np.abs(tx_dpd)),
+ np.mean(np.abs(tx_dpd)),
+ np.median(np.abs(tx_dpd))
+ ), size=8)
+ ax.legend(loc=4)
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ ccdf_min, ccdf_max = 0, 1
+ tx_hist, ccdf_edges = np.histogram(np.abs(tx_dpd),
+ bins=60,
+ range=(ccdf_min, ccdf_max))
+ tx_hist_normalized = tx_hist.astype(float) / np.sum(tx_hist)
+ ccdf = 1.0 - np.cumsum(tx_hist_normalized)
+ ax.semilogy(ccdf_edges[:-1], ccdf, label="CCDF")
+ ax.semilogy(ccdf_edges[:-1],
+ tx_hist_normalized,
+ label="Histogram",
+ drawstyle='steps')
+ ax.legend(loc=4)
+ ax.set_ylim(1e-5, 2)
+ ax.set_title("Complementary Cumulative Distribution Function")
+ ax.set_xlabel("TX Amplitude")
+ ax.set_ylabel("Ratio of Samples larger than x")
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ ax.semilogy(np.array(self.mses_am) + 1e-10, label="log(MSE)")
+ ax.semilogy(np.array(self.errs_am) + 1e-10, label="log(ERR)")
+ ax.legend(loc=4)
+ ax.set_title("MSE History")
+ ax.set_xlabel("Iterations")
+ ax.set_ylabel("MSE")
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ ax.plot(self.tx_mers, label="TX MER")
+ ax.plot(self.rx_mers, label="RX MER")
+ ax.legend(loc=4)
+ ax.set_title("MER History")
+ ax.set_xlabel("Iterations")
+ ax.set_ylabel("MER")
+
+ ax = plt.subplot(4, 2, i_sub)
+ rx_range = np.linspace(0, 1, num=100, dtype=np.complex64)
+ rx_range_dpd = self.dpd_amplitude(rx_range)[0]
+ rx_range_dpd_new = self.dpd_amplitude(rx_range, new_coefs_am)[0]
+ i_sub += 1
+ ax.scatter(
+ np.abs(tx_choice),
+ np.abs(rx_choice),
+ s=0.1)
+ ax.plot(rx_range_dpd / self.coefs_am[0], rx_range, linewidth=0.25, label="current")
+ ax.plot(rx_range_dpd_new / self.coefs_am[0], rx_range, linewidth=0.25, label="next")
+ ax.set_ylim(0, 1)
+ ax.set_xlim(0, 1)
+ ax.legend()
+ ax.set_title("Amplifier Characteristic")
+ ax.set_xlabel("TX Amplitude")
+ ax.set_ylabel("RX Amplitude")
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ coefs_am_history = np.array(self.coefs_am_history)
+ for idx, coef_hist in enumerate(coefs_am_history.T):
+ ax.plot(coef_hist,
+ label="Coef {}".format(idx),
+ linewidth=0.5)
+ ax.legend(loc=4)
+ ax.set_title("AM/AM Coefficient History")
+ ax.set_xlabel("Iterations")
+ ax.set_ylabel("Coefficient Value")
+
+ phase_range = np.linspace(0, 1, num=100, dtype=np.complex64)
+ phase_range_dpd = self.dpd_phase(phase_range)[0]
+ phase_range_dpd_new = self.dpd_phase(phase_range,
+ coefs=new_coefs_pm)[0]
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ ax.scatter(
+ np.abs(tx_choice),
+ np.rad2deg(phase_diff_choice),
+ s=0.1)
+ ax.plot(
+ np.abs(phase_range),
+ np.rad2deg(phase_range_dpd),
+ linewidth=0.25,
+ label="current")
+ ax.plot(
+ np.abs(phase_range),
+ np.rad2deg(phase_range_dpd_new),
+ linewidth=0.25,
+ label="next")
+ ax.set_ylim(-60, 60)
+ ax.set_xlim(0, 1)
+ ax.legend()
+ ax.set_title("Amplifier Characteristic")
+ ax.set_xlabel("TX Amplitude")
+ ax.set_ylabel("Phase Difference")
+
+ ax = plt.subplot(4, 2, i_sub)
+ i_sub += 1
+ coefs_pm_history = np.array(self.coefs_pm_history)
+ for idx, coef_phase_hist in enumerate(coefs_pm_history.T):
+ ax.plot(coef_phase_hist,
+ label="Coef {}".format(idx),
+ linewidth=0.5)
+ ax.legend(loc=4)
+ ax.set_title("AM/PM Coefficient History")
+ ax.set_xlabel("Iterations")
+ ax.set_ylabel("Coefficient Value")
+
+ fig.tight_layout()
+ fig.savefig(fig_path)
+ fig.clf()
+
+ self.coefs_am = new_coefs_am
+ self.coefs_am_history.append(self.coefs_am)
+ self.coefs_pm = new_coefs_pm
+ self.coefs_pm_history.append(self.coefs_pm)
+ return self.coefs_am, self.coefs_pm
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/Symbol_align.py b/dpd/src/Symbol_align.py
new file mode 100644
index 0000000..05a9049
--- /dev/null
+++ b/dpd/src/Symbol_align.py
@@ -0,0 +1,191 @@
+# -*- coding: utf-8 -*-
+#
+# Modulation Error Rate
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+import time
+try:
+ logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+except:
+ logging_path = "/tmp/"
+
+import numpy as np
+import src.const
+import scipy
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
+class Symbol_align:
+ """
+ Find the phase offset to the start of the DAB symbols in an
+ unaligned dab signal.
+ """
+ def __init__(self, sample_rate):
+ self.c = src.const.const(sample_rate)
+ pass
+
+ def _calc_offset_to_first_symbol_without_prefix(self, tx, debug=False):
+ tx_orig = tx[0:-self.c.T_U]
+ tx_cut_prefix = tx[self.c.T_U:]
+
+ tx_product = np.abs(tx_orig - tx_cut_prefix)
+ tx_product_avg = np.correlate(
+ tx_product,
+ np.ones(self.c.T_C),
+ mode='valid')
+ tx_product_avg_min_filt = \
+ scipy.ndimage.filters.minimum_filter1d(
+ tx_product_avg,
+ int(1.5 * self.c.T_S)
+ )
+ peaks = np.ravel(np.where(tx_product_avg == tx_product_avg_min_filt))
+
+ offset = peaks[np.argmin([tx_product_avg[peak] for peak in peaks])]
+
+ if debug:
+ fig = plt.figure(figsize=(9, 9))
+
+ ax = fig.add_subplot(4, 1, 1)
+ ax.plot(tx_product)
+ ylim = ax.get_ylim()
+ for peak in peaks:
+ ax.plot((peak, peak), (ylim[0], ylim[1]))
+ if peak == offset:
+ ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
+ else:
+ ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
+ ax.set_xlabel("Sample")
+ ax.set_ylabel("Conj. Product")
+ ax.set_title("Difference with shifted self")
+
+ ax = fig.add_subplot(4, 1, 2)
+ ax.plot(tx_product_avg)
+ ylim = ax.get_ylim()
+ for peak in peaks:
+ ax.plot((peak, peak), (ylim[0], ylim[1]))
+ if peak == offset:
+ ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
+ else:
+ ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
+ ax.set_xlabel("Sample")
+ ax.set_ylabel("Conj. Product")
+ ax.set_title("Moving Average")
+
+ ax = fig.add_subplot(4, 1, 3)
+ ax.plot(tx_product_avg_min_filt)
+ ylim = ax.get_ylim()
+ for peak in peaks:
+ ax.plot((peak, peak), (ylim[0], ylim[1]))
+ if peak == offset:
+ ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
+ else:
+ ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
+ ax.set_xlabel("Sample")
+ ax.set_ylabel("Conj. Product")
+ ax.set_title("Min Filter")
+
+ ax = fig.add_subplot(4, 1, 4)
+ tx_product_crop = tx_product[peaks[0]-50:peaks[0]+50]
+ x = range(tx_product.shape[0])[peaks[0]-50:peaks[0]+50]
+ ax.plot(x, tx_product_crop)
+ ylim = ax.get_ylim()
+ ax.plot((peaks[0], peaks[0]), (ylim[0], ylim[1]))
+ ax.set_xlabel("Sample")
+ ax.set_ylabel("Conj. Product")
+ ax.set_title("Difference with shifted self")
+ fig.tight_layout()
+
+ # "offset" measures where the shifted signal matches the
+ # original signal. Therefore we have to subtract the size
+ # of the shift to find the offset of the symbol start.
+ return (offset + self.c.T_C) % self.c.T_S
+
+ def _remove_outliers(self, x, stds=5):
+ deviation_from_mean = np.abs(x - np.mean(x))
+ inlier_idxs = deviation_from_mean < stds * np.std(x)
+ x = x[inlier_idxs]
+ return x
+
+ def _calc_delta_angle(self, fft, debug=False):
+ # Introduce invariance against carrier
+ angles = np.angle(fft) % (np.pi / 2.)
+
+ # Calculate Angle difference and compensate jumps
+ deltas_angle = np.diff(angles)
+ deltas_angle[deltas_angle > np.pi/4.] =\
+ deltas_angle[deltas_angle > np.pi/4.] - np.pi/2.
+ deltas_angle[-deltas_angle > np.pi/4.] = \
+ deltas_angle[-deltas_angle > np.pi/4.] + np.pi/2.
+ deltas_angle = self._remove_outliers(deltas_angle)
+
+ delta_angle = np.mean(deltas_angle)
+ delta_angle_std = np.std(deltas_angle)
+ if debug:
+ plt.subplot(211)
+ plt.plot(angles, 'p')
+ plt.subplot(212)
+ plt.plot(deltas_angle, 'p')
+ return delta_angle
+
+ def _delta_angle_to_samples(self, angle):
+ return - angle / self.c.phase_offset_per_sample
+
+ def _calc_sample_offset(self, sig, debug=False):
+ assert sig.shape[0] == self.c.T_U,\
+ "Input length is not a Symbol without cyclic prefix"
+
+ fft = np.fft.fftshift(np.fft.fft(sig))
+ fft_crop = np.delete(fft[self.c.FFT_start:self.c.FFT_end], self.c.FFT_delete)
+ delta_angle = self._calc_delta_angle(fft_crop, debug=debug)
+ delta_sample = self._delta_angle_to_samples(delta_angle)
+ delta_sample_int = np.round(delta_sample).astype(int)
+ error = np.abs(delta_sample_int - delta_sample)
+ if error > 0.1:
+ raise RuntimeError("Could not calculate sample offset")
+ return delta_sample_int
+
+ def calc_offset(self, tx):
+ off_sym = self._calc_offset_to_first_symbol_without_prefix(
+ tx, debug=False)
+ off_sam = self._calc_sample_offset(
+ tx[off_sym:off_sym + self.c.T_U])
+ off = (off_sym + off_sam) % self.c.T_S
+
+ assert self._calc_sample_offset(tx[off:off + self.c.T_U]) == 0, \
+ "Failed to calculate offset"
+ return off
+
+ def crop_symbol_without_cyclic_prefix(self, tx):
+ off = self.calc_offset(tx)
+ return tx[
+ off:
+ off+self.c.T_U
+ ]
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/TX_Agc.py b/dpd/src/TX_Agc.py
new file mode 100644
index 0000000..9274612
--- /dev/null
+++ b/dpd/src/TX_Agc.py
@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+#
+# Automatic Gain Control
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+import time
+
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+import matplotlib
+
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
+import src.Adapt as Adapt
+
+#TODO fix for float tx_gain
+class TX_Agc:
+ def __init__(self,
+ adapt,
+ max_txgain=89,
+ tx_median_target=0.1,
+ tx_median_threshold_max=0.12,
+ tx_median_threshold_min=0.08):
+ """
+ In order to avoid digital clipping, this class increases the
+ TX gain and reduces the digital gain. Digital clipping happens
+ when the digital analog converter receives values greater than
+ it's maximal output. This class solves that problem by adapting
+ the TX gain in a way that the peaks of the TX signal are in a
+ specified range. The TX gain is adapted accordingly. The TX peaks
+ are approximated by estimating it based on the signal median.
+
+ :param adapt: Instance of Adapt Class to update
+ txgain and coefficients
+ :param max_txgain: limit for TX gain
+ :param tx_median_threshold_max: if the median of TX is larger
+ than this value, then the digital gain is reduced
+ :param tx_median_threshold_min: if the median of TX is smaller
+ than this value, then the digital gain is increased
+ :param tx_median_target: The digital gain is reduced in a way that
+ the median TX value is expected to be lower than this value.
+ """
+ assert isinstance(adapt, Adapt.Adapt)
+ self.adapt = adapt
+ self.max_txgain = max_txgain
+ self.txgain = self.max_txgain
+
+ assert tx_median_threshold_max > tx_median_target,\
+ "The tolerated tx_median has to be larger then the goal tx_median"
+ self.tx_median_threshold_tolerate_max = tx_median_threshold_max
+ self.tx_median_threshold_tolerate_min = tx_median_threshold_min
+ self.tx_median_target = tx_median_target
+
+ def adapt_if_necessary(self, tx):
+ tx_median = np.median(np.abs(tx))
+ if tx_median > self.tx_median_threshold_tolerate_max or\
+ tx_median < self.tx_median_threshold_tolerate_min:
+ delta_db = 20 * np.log10(self.tx_median_target / tx_median)
+ new_txgain = self.adapt.get_txgain() - delta_db
+ assert new_txgain < self.max_txgain,\
+ "TX_Agc failed. New TX gain of {} is too large.".format(
+ new_txgain
+ )
+ self.adapt.set_txgain(new_txgain)
+ txgain = self.adapt.get_txgain()
+
+ digital_gain_factor = 10 ** (delta_db / 20.)
+ digital_gain = self.adapt.get_digital_gain() * digital_gain_factor
+ self.adapt.set_digital_gain(digital_gain)
+
+ logging.info(
+ "digital_gain = {}, txgain_new = {}, "\
+ "delta_db = {}, tx_median {}, "\
+ "digital_gain_factor = {}".
+ format(digital_gain, txgain, delta_db,
+ tx_median, digital_gain_factor))
+
+ time.sleep(1)
+ return True
+ return False
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/Test_data.py b/dpd/src/Test_data.py
new file mode 100644
index 0000000..9dd0913
--- /dev/null
+++ b/dpd/src/Test_data.py
@@ -0,0 +1,136 @@
+# -*- coding: utf-8 -*-
+#
+# Modulation Error Rate
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+import datetime
+import os
+import logging
+import time
+try:
+ logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+except:
+ logging_path = "/tmp/"
+
+import src.const
+import src.Dab_Util
+import numpy as np
+import matplotlib
+matplotlib.use('agg')
+import matplotlib.pyplot as plt
+
+
+class Test_data:
+ def __init__(self, sample_rate, type):
+ """
+ Standardized access to complex64 test data files.
+
+ :param sample_rate:
+ :param type:
+
+ Testing:
+ TD = src.Test_data.Test_data(8192000, 'file')
+ tx_orig = TD.get_symbol(0,0)
+ fig = plt.figure(figsize=(9,6))
+ ax = fig.add_subplot(2,1,1)
+ ax.plot(tx_orig)
+ ax = fig.add_subplot(2,1,2)
+ plt.plot(np.angle(np.fft.fftshift(np.fft.fft(tx_orig))), 'p')
+ """
+
+ self.c = src.const.const(sample_rate)
+ self.du = src.Dab_Util.Dab_Util(sample_rate)
+
+ self.file_paths = {
+ (2048000, 'file'):
+ ("./test_data/odr-dabmod_to_file_2048_NoFir_noDPD.iq",
+ (
+ self.c.T_F + # Pipelineing
+ self.c.T_NULL + # NULL Symbol
+ self.c.T_S + # Synchronization Symbol
+ self.c.T_C # Cyclic Prefix
+ )),
+ (8192000, 'file'):
+ ("./test_data/odr-dabmod_to_file_8192_NoFir_noDPD.iq",
+ (
+ self.c.T_F + # Pipelining
+ self.c.T_U + # Pipelining Resampler TODO(?)
+ self.c.T_NULL + # NULL Symbol
+ self.c.T_S + # Synchronization Symbol
+ self.c.T_C # Cyclic Prefix
+ )),
+ (8192000, 'rec_noFir'):
+ ("./test_data/odr-dabmod_reconding_8192_NoFir_DPD_2104.iq",
+ ( 64 )),
+ (8192000, 'rec_fir'):
+ ("./test_data/odr-dabmod_reconding_8192_Fir_DPD_2104.iq",
+ ( 232 )),
+ }
+
+ config = (sample_rate, type)
+ if not config in self.file_paths.keys():
+ raise RuntimeError("Configuration not found, possible are:\n {}".
+ format(self.file_paths))
+
+ self.path, self.file_offset = self.file_paths[(sample_rate, type)]
+
+ def _load_from_file(self, offset, length):
+ print(offset, length, self.file_offset)
+ return self.du.fromfile(
+ self.path,
+ length=length,
+ offset=offset + self.file_offset)
+
+ def get_symbol_without_prefix(self,
+ frame_idx=0,
+ symbol_idx=0,
+ off=0):
+ return self._load_from_file(
+ frame_idx*self.c.T_F +
+ symbol_idx*self.c.T_S +
+ off,
+ self.c.T_U)
+
+ def get_symbol_with_prefix(self,
+ frame_idx=0,
+ symbol_idx=0,
+ n_symbols=1,
+ off=0):
+ offset = (
+ frame_idx*self.c.T_F +
+ symbol_idx*self.c.T_S -
+ self.c.T_C +
+ off)
+ return self._load_from_file( offset, self.c.T_S * n_symbols)
+
+ def get_file_length_in_symbols(self):
+ symbol_size = float(
+ 64/8 * # complex 64
+ self.c.T_S # Symbol Size
+ )
+ return os.path.getsize(self.path) / symbol_size
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/__init__.py b/dpd/src/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/dpd/src/__init__.py
diff --git a/dpd/src/const.py b/dpd/src/const.py
new file mode 100644
index 0000000..1011c6c
--- /dev/null
+++ b/dpd/src/const.py
@@ -0,0 +1,38 @@
+# DAB Frame constants
+# Sources:
+# - etsi_EN_300_401_v010401p p145
+# - Measured with USRP B200
+
+import numpy as np
+
+class const:
+ def __init__(self, sample_rate):
+ self.sample_rate = sample_rate
+
+ # Time domain
+ self.T_F = sample_rate / 2048000 * 196608 # Transmission frame duration
+ self.T_NULL = sample_rate / 2048000 * 2656 # Null symbol duration
+ self.T_S = sample_rate / 2048000 * 2552 # Duration of OFDM symbols of indices l = 1, 2, 3,... L;
+ self.T_U = sample_rate / 2048000 * 2048 # Inverse of carrier spacing
+ self.T_C = sample_rate / 2048000 * 504 # Duration of cyclic prefix
+
+ # Frequency Domain
+ # example: np.delete(fft[3328:4865], 768)
+ self.FFT_delete = 768
+ self.FFT_delta = 1536 # Number of carrier frequencies
+ if sample_rate == 2048000:
+ self.FFT_start = 256
+ self.FFT_end = 1793
+ elif sample_rate == 8192000:
+ self.FFT_start = 3328
+ self.FFT_end = 4865
+ else:
+ raise RuntimeError("Sample Rate '{}' not supported".format(
+ sample_rate
+ ))
+
+ # Calculate sample offset from phase rotation
+ # time per sample = 1 / sample_rate
+ # frequency per bin = 1kHz
+ # phase difference per sample offset = delta_t * 2 * pi * delta_freq
+ self.phase_offset_per_sample = 1. / sample_rate * 2 * np.pi * 1000
diff --git a/dpd/src/phase_align.py b/dpd/src/phase_align.py
new file mode 100644
index 0000000..7f82392
--- /dev/null
+++ b/dpd/src/phase_align.py
@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, phase-align a signal against a reference.
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+from scipy import signal, optimize
+import sys
+import matplotlib.pyplot as plt
+
+
+def phase_align(sig, ref_sig, plot=False):
+ """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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and plot:
+ dt = datetime.datetime.now().isoformat()
+ fig_path = logging_path + "/" + dt + "_phase_align.svg"
+
+ 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 logging.getLogger().getEffectiveLevel() == logging.DEBUG and plot:
+ 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.clf()
+
+ return sig
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/subsample_align.py b/dpd/src/subsample_align.py
new file mode 100755
index 0000000..6d1cd2a
--- /dev/null
+++ b/dpd/src/subsample_align.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, utility to do subsample alignment
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+import datetime
+import os
+import logging
+logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
+
+import numpy as np
+from scipy import signal, 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=False):
+ """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:
+ corr = np.vectorize(correlate_for_delay)
+ ixs = np.linspace(-1, 1, 100)
+ taus = corr(ixs)
+
+ dt = datetime.datetime.now().isoformat()
+ tau_path = (logging_path + "/" + dt + "_tau.svg")
+ plt.plot(ixs, taus)
+ plt.title("Subsample correlation, minimum is best: {}".format(best_tau))
+ plt.savefig(tau_path)
+ plt.clf()
+
+ # 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)
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/test_dab_Util.py b/dpd/src/test_dab_Util.py
new file mode 100644
index 0000000..ec15586
--- /dev/null
+++ b/dpd/src/test_dab_Util.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+#
+# Test code for DAB util
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+
+from unittest import TestCase
+
+import numpy as np
+import pandas as pd
+import src.Dab_Util as DU
+
+class TestDab_Util(TestCase):
+
+ def test_subsample_align(self, sample_orig=r'../test_data/orig_rough_aligned.dat',
+ sample_rec =r'../test_data/recored_rough_aligned.dat',
+ length = 10240, max_size = 1000000):
+ du = DU.Dab_Util(8196000)
+ res1 = []
+ res2 = []
+ for i in range(10):
+ start = np.random.randint(50, max_size)
+ r = np.random.randint(-50, 50)
+
+ s1 = du.fromfile(sample_orig, offset=start+r, length=length)
+ s2 = du.fromfile(sample_rec, offset=start, length=length)
+
+ res1.append(du.lag_upsampling(s2, s1, 32))
+
+ s1_aligned, s2_aligned = du.subsample_align(s1, s2)
+
+ res2.append(du.lag_upsampling(s2_aligned, s1_aligned, 32))
+
+ error_rate = np.mean(np.array(res2) != 0)
+ self.assertEqual(error_rate, 0.0, "The error rate for aligning was %.2f%%"
+ % error_rate * 100)
+
+#def test_using_aligned_pair(sample_orig=r'../data/orig_rough_aligned.dat', sample_rec =r'../data/recored_rough_aligned.dat', length = 10240, max_size = 1000000):
+# res = []
+# for i in tqdm(range(100)):
+# start = np.random.randint(50, max_size)
+# r = np.random.randint(-50, 50)
+#
+# s1 = du.fromfile(sample_orig, offset=start+r, length=length)
+# s2 = du.fromfile(sample_rec, offset=start, length=length)
+#
+# res.append({'offset':r,
+# '1':r - du.lag_upsampling(s2, s1, n_up=1),
+# '2':r - du.lag_upsampling(s2, s1, n_up=2),
+# '3':r - du.lag_upsampling(s2, s1, n_up=3),
+# '4':r - du.lag_upsampling(s2, s1, n_up=4),
+# '8':r - du.lag_upsampling(s2, s1, n_up=8),
+# '16':r - du.lag_upsampling(s2, s1, n_up=16),
+# '32':r - du.lag_upsampling(s2, s1, n_up=32),
+# })
+# df = pd.DataFrame(res)
+# df = df.reindex_axis(sorted(df.columns), axis=1)
+# print(df.describe())
+#
+#
+#print("Align using upsampling")
+#for n_up in [1, 2, 3, 4, 7, 8, 16]:
+# correct_ratio = test_phase_offset(lambda x,y: du.lag_upsampling(x,y,n_up), tol=1./n_up)
+# print("%.1f%% of the tested offsets were measured within tolerance %.4f for n_up = %d" % (correct_ratio * 100, 1./n_up, n_up))
+#test_using_aligned_pair()
+#
+#print("Phase alignment")
+#test_subsample_alignment()
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/test_measure.py b/dpd/src/test_measure.py
new file mode 100644
index 0000000..ee3cfdb
--- /dev/null
+++ b/dpd/src/test_measure.py
@@ -0,0 +1,62 @@
+# -*- coding: utf-8 -*-
+#
+# DPD Calculation Engine, test case for measure
+#
+# http://www.opendigitalradio.org
+# Licence: The MIT License, see notice at the end of this file
+from unittest import TestCase
+from Measure import Measure
+import socket
+
+
+class TestMeasure(TestCase):
+
+ def _open_socks(self):
+ sock_server = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+ sock_server.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
+ sock_server.bind(('localhost', 1234))
+ sock_server.listen(1)
+
+ sock_client = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+ sock_client.connect(('localhost', 1234))
+
+ conn_server, addr_server = sock_server.accept()
+ return conn_server, sock_client
+
+ def test__recv_exact(self):
+ m = Measure(1234, 1)
+ payload = b"test payload"
+
+ conn_server, sock_client = self._open_socks()
+ conn_server.send(payload)
+ rec = m._recv_exact(sock_client, len(payload))
+
+ self.assertEqual(rec, payload,
+ "Did not receive the same message as sended. (%s, %s)" %
+ (rec, payload))
+
+ def test_get_samples(self):
+ self.fail()
+
+
+# The MIT License (MIT)
+#
+# Copyright (c) 2017 Andreas Steger
+#
+# 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/src/ConfigParser.cpp b/src/ConfigParser.cpp
index 459811f..9ac1280 100644
--- a/src/ConfigParser.cpp
+++ b/src/ConfigParser.cpp
@@ -171,7 +171,10 @@ static void parse_configfile(
// Poly coefficients:
if (pt.get("poly.enabled", 0) == 1) {
mod_settings.polyCoefFilename =
- pt.get<std::string>("poly.polycoeffile", "default");
+ pt.get<std::string>("poly.polycoeffile", "dpd/poly.coef");
+
+ mod_settings.polyNumThreads =
+ pt.get<int>("poly.num_threads", 0);
}
// Output options
diff --git a/src/ConfigParser.h b/src/ConfigParser.h
index 22a4fc5..89f0fb7 100644
--- a/src/ConfigParser.h
+++ b/src/ConfigParser.h
@@ -75,7 +75,7 @@ struct mod_settings_t {
std::string filterTapsFilename = "";
std::string polyCoefFilename = "";
-
+ unsigned polyNumThreads = 0;
#if defined(HAVE_OUTPUT_UHD)
OutputUHDConfig outputuhd_conf;
diff --git a/src/DabModulator.cpp b/src/DabModulator.cpp
index 34d8e66..cc2642a 100644
--- a/src/DabModulator.cpp
+++ b/src/DabModulator.cpp
@@ -205,13 +205,8 @@ int DabModulator::process(Buffer* dataOut)
shared_ptr<MemlessPoly> cifPoly;
if (not m_settings.polyCoefFilename.empty()) {
- cifPoly = make_shared<MemlessPoly>(m_settings.polyCoefFilename);
- etiLog.level(debug) << m_settings.polyCoefFilename << "\n";
- etiLog.level(debug) << cifPoly->m_coefs[0] << " " <<
- cifPoly->m_coefs[1] << " "<< cifPoly->m_coefs[2] << " "<<
- cifPoly->m_coefs[3] << " "<< cifPoly->m_coefs[4] << " "<<
- cifPoly->m_coefs[5] << " "<< cifPoly->m_coefs[6] << " "<<
- cifPoly->m_coefs[7] << "\n";
+ cifPoly = make_shared<MemlessPoly>(m_settings.polyCoefFilename,
+ m_settings.polyNumThreads);
rcs.enrol(cifPoly.get());
}
diff --git a/src/MemlessPoly.cpp b/src/MemlessPoly.cpp
index 7e074eb..d7f9a96 100644
--- a/src/MemlessPoly.cpp
+++ b/src/MemlessPoly.cpp
@@ -29,6 +29,8 @@
along with ODR-DabMod. If not, see <http://www.gnu.org/licenses/>.
*/
+#pragma GCC optimize ("O3")
+
#include "MemlessPoly.h"
#include "PcDebug.h"
#include "Utils.h"
@@ -36,31 +38,60 @@
#include <stdio.h>
#include <stdexcept>
+#include <future>
#include <array>
#include <iostream>
#include <fstream>
#include <memory>
+#include <complex>
using namespace std;
+// Number of AM/AM coefs, identical to number of AM/PM coefs
+#define NUM_COEFS 5
-// By default the signal is unchanged
-static const std::array<float, 8> default_coefficients({
- 1, 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0, 0.0
- });
-
-
-MemlessPoly::MemlessPoly(const std::string& coefs_file) :
+MemlessPoly::MemlessPoly(const std::string& coefs_file, unsigned int num_threads) :
PipelinedModCodec(),
RemoteControllable("memlesspoly"),
- m_coefs_file(coefs_file)
+ m_coefs_am(),
+ m_coefs_pm(),
+ m_coefs_file(coefs_file),
+ m_coefs_mutex()
{
PDEBUG("MemlessPoly::MemlessPoly(%s) @ %p\n",
coefs_file.c_str(), this);
RC_ADD_PARAMETER(ncoefs, "(Read-only) number of coefficients.");
- RC_ADD_PARAMETER(coeffile, "Filename containing coefficients. When written to, the new file gets automatically loaded.");
+ RC_ADD_PARAMETER(coeffile, "Filename containing coefficients. "
+ "When set, the file gets loaded.");
+
+ if (num_threads == 0) {
+ const unsigned int hw_concurrency = std::thread::hardware_concurrency();
+ etiLog.level(info) << "Polynomial Predistorter will use " <<
+ hw_concurrency << " threads (auto detected)";
+
+ for (size_t i = 0; i < hw_concurrency; i++) {
+ m_workers.emplace_back();
+ }
+
+ for (auto& worker : m_workers) {
+ worker.thread = std::thread(
+ &MemlessPoly::worker_thread, &worker);
+ }
+ }
+ else {
+ etiLog.level(info) << "Polynomial Predistorter will use " <<
+ num_threads << " threads (set in config file)";
+
+ for (size_t i = 0; i < num_threads; i++) {
+ m_workers.emplace_back();
+ }
+
+ for (auto& worker : m_workers) {
+ worker.thread = std::thread(
+ &MemlessPoly::worker_thread, &worker);
+ }
+ }
load_coefficients(m_coefs_file);
@@ -69,76 +100,165 @@ MemlessPoly::MemlessPoly(const std::string& coefs_file) :
void MemlessPoly::load_coefficients(const std::string &coefFile)
{
- std::vector<float> coefs;
- if (coefFile == "default") {
- std::copy(default_coefficients.begin(), default_coefficients.end(),
- std::back_inserter(coefs));
+ std::vector<float> coefs_am;
+ std::vector<float> coefs_pm;
+ std::ifstream coef_fstream(coefFile.c_str());
+ if (!coef_fstream) {
+ throw std::runtime_error("MemlessPoly: Could not open file with coefs!");
}
- else {
- std::ifstream coef_fstream(coefFile.c_str());
- if(!coef_fstream) {
- fprintf(stderr, "MemlessPoly: file %s could not be opened !\n", coefFile.c_str());
- throw std::runtime_error("MemlessPoly: Could not open file with coefs! ");
- }
- int n_coefs;
- coef_fstream >> n_coefs;
+ int n_coefs;
+ coef_fstream >> n_coefs;
- if (n_coefs <= 0) {
- fprintf(stderr, "MemlessPoly: warning: coefs file has invalid format\n");
- throw std::runtime_error("MemlessPoly: coefs file has invalid format.");
- }
+ if (n_coefs <= 0) {
+ throw std::runtime_error("MemlessPoly: coefs file has invalid format.");
+ }
+ else if (n_coefs != NUM_COEFS) {
+ throw std::runtime_error("MemlessPoly: invalid number of coefs: " +
+ std::to_string(n_coefs) + " expected " + std::to_string(NUM_COEFS));
+ }
- if (n_coefs != 8) {
- throw std::runtime_error( "MemlessPoly: error: coefs file does not have 8 coefs\n");
- }
+ const int n_entries = 2 * n_coefs;
- fprintf(stderr, "MemlessPoly: Reading %d coefs...\n", n_coefs);
+ etiLog.log(debug, "MemlessPoly: Reading %d coefs...", n_entries);
- coefs.resize(n_coefs);
+ coefs_am.resize(n_coefs);
+ coefs_pm.resize(n_coefs);
- int n;
- for (n = 0; n < n_coefs; n++) {
- coef_fstream >> coefs[n];
- PDEBUG("MemlessPoly: coef: %f\n", coefs[n] );
- if (coef_fstream.eof()) {
- fprintf(stderr, "MemlessPoly: file %s should contains %d coefs, but EOF reached "\
- "after %d coefs !\n", coefFile.c_str(), n_coefs, n);
- throw std::runtime_error("MemlessPoly: coefs file invalid ! ");
- }
+ for (int n = 0; n < n_entries; n++) {
+ float a;
+ coef_fstream >> a;
+
+ if (n < n_coefs) {
+ coefs_am[n] = a;
+ }
+ else {
+ coefs_pm[n - n_coefs] = a;
+ }
+
+ if (coef_fstream.eof()) {
+ etiLog.log(error, "MemlessPoly: file %s should contains %d coefs, "
+ "but EOF reached after %d coefs !",
+ coefFile.c_str(), n_entries, n);
+ throw std::runtime_error("MemlessPoly: coefs file invalid !");
}
}
{
std::lock_guard<std::mutex> lock(m_coefs_mutex);
- m_coefs = coefs;
+ m_coefs_am = coefs_am;
+ m_coefs_pm = coefs_pm;
}
}
+/* The restrict keyword is C99, g++ and clang++ however support __restrict
+ * instead, and this allows the compiler to auto-vectorize the loop.
+ */
+static void apply_coeff(
+ const float *__restrict coefs_am, const float *__restrict coefs_pm,
+ const complexf *__restrict in, size_t start, size_t stop,
+ complexf *__restrict out)
+{
+ for (size_t i = start; i < stop; i+=1) {
+
+ float in_mag_sq = in[i].real() * in[i].real() + in[i].imag() * in[i].imag();
+
+ float amplitude_correction =
+ ( coefs_am[0] + in_mag_sq *
+ ( coefs_am[1] + in_mag_sq *
+ ( coefs_am[2] + in_mag_sq *
+ ( coefs_am[3] + in_mag_sq *
+ coefs_am[4]))));
+
+ float phase_correction = -1 *
+ ( coefs_pm[0] + in_mag_sq *
+ ( coefs_pm[1] + in_mag_sq *
+ ( coefs_pm[2] + in_mag_sq *
+ ( coefs_pm[3] + in_mag_sq *
+ coefs_pm[4]))));
+
+ float phase_correction_sq = phase_correction * phase_correction;
+
+ // Approximation for Cosinus 1 - 1/2 x^2 + 1/24 x^4 - 1/720 x^6
+ float re = (1.0f - phase_correction_sq *
+ ( -0.5f + phase_correction_sq *
+ ( 0.486666f + phase_correction_sq *
+ ( -0.00138888f))));
+
+ // Approximation for Sinus x + 1/6 x^3 + 1/120 x^5
+ float im = phase_correction *
+ (1.0f + phase_correction_sq *
+ (0.166666f + phase_correction_sq *
+ (0.00833333f)));
+
+ out[i] = in[i] * amplitude_correction * complex<float>(re, im);
+ }
+}
+
+void MemlessPoly::worker_thread(MemlessPoly::worker_t *workerdata)
+{
+ while (true) {
+ worker_t::input_data_t in_data;
+ workerdata->in_queue.wait_and_pop(in_data);
+
+ if (in_data.terminate) {
+ break;
+ }
+
+ apply_coeff(in_data.coefs_am, in_data.coefs_pm,
+ in_data.in, in_data.start, in_data.stop,
+ in_data.out);
+
+ workerdata->out_queue.push(1);
+ }
+}
int MemlessPoly::internal_process(Buffer* const dataIn, Buffer* dataOut)
{
- const float* in = reinterpret_cast<const float*>(dataIn->getData());
- float* out = reinterpret_cast<float*>(dataOut->getData());
- size_t sizeIn = dataIn->getLength() / sizeof(float);
-
- {
- std::lock_guard<std::mutex> lock(m_coefs_mutex);
- for (size_t i = 0; i < sizeIn; i += 1) {
- float mag = std::abs(in[i]);
- //out[i] = in[i];
- out[i] = in[i] * (
- m_coefs[0] +
- m_coefs[1] * mag +
- m_coefs[2] * mag*mag +
- m_coefs[3] * mag*mag*mag +
- m_coefs[4] * mag*mag*mag*mag +
- m_coefs[5] * mag*mag*mag*mag*mag +
- m_coefs[6] * mag*mag*mag*mag*mag*mag +
- m_coefs[7] * mag*mag*mag*mag*mag*mag*mag
- );
+ dataOut->setLength(dataIn->getLength());
+
+ const complexf* in = reinterpret_cast<const complexf*>(dataIn->getData());
+ complexf* out = reinterpret_cast<complexf*>(dataOut->getData());
+ size_t sizeOut = dataOut->getLength() / sizeof(complexf);
+
+ {
+ std::lock_guard<std::mutex> lock(m_coefs_mutex);
+ const size_t num_threads = m_workers.size();
+
+ if (num_threads > 0) {
+ const size_t step = sizeOut / num_threads;
+
+ size_t start = 0;
+ for (auto& worker : m_workers) {
+ worker_t::input_data_t dat;
+ dat.terminate = false;
+ dat.coefs_am = m_coefs_am.data();
+ dat.coefs_pm = m_coefs_pm.data();
+ dat.in = in;
+ dat.start = start;
+ dat.stop = start + step;
+ dat.out = out;
+
+ worker.in_queue.push(dat);
+
+ start += step;
+ }
+
+ // Do the last in this thread
+ apply_coeff(m_coefs_am.data(), m_coefs_pm.data(),
+ in, start, sizeOut, out);
+
+ // Wait for completion of the tasks
+ for (auto& worker : m_workers) {
+ int ret;
+ worker.out_queue.wait_and_pop(ret);
}
}
+ else {
+ apply_coeff(m_coefs_am.data(), m_coefs_pm.data(),
+ in, 0, sizeOut, out);
+ }
+ }
return dataOut->getLength();
}
@@ -172,9 +292,9 @@ const string MemlessPoly::get_parameter(const string& parameter) const
{
stringstream ss;
if (parameter == "ncoefs") {
- ss << m_coefs.size();
+ ss << m_coefs_am.size();
}
- else if (parameter == "coefFile") {
+ else if (parameter == "coeffile") {
ss << m_coefs_file;
}
else {
diff --git a/src/MemlessPoly.h b/src/MemlessPoly.h
index 210b4b4..57c0924 100644
--- a/src/MemlessPoly.h
+++ b/src/MemlessPoly.h
@@ -52,7 +52,7 @@ typedef std::complex<float> complexf;
class MemlessPoly : public PipelinedModCodec, public RemoteControllable
{
public:
- MemlessPoly(const std::string& coefs_file);
+ MemlessPoly(const std::string& coefs_file, unsigned int num_threads);
virtual const char* name() { return "MemlessPoly"; }
@@ -63,16 +63,58 @@ public:
virtual const std::string get_parameter(
const std::string& parameter) const;
-//TODO to protected
- std::vector<float> m_coefs;
-
-
-protected:
+private:
int internal_process(Buffer* const dataIn, Buffer* dataOut);
void load_coefficients(const std::string &coefFile);
+ struct worker_t {
+ struct input_data_t {
+ bool terminate = false;
+
+ const float *coefs_am = nullptr;
+ const float *coefs_pm = nullptr;
+ const complexf *in = nullptr;
+ size_t start = 0;
+ size_t stop = 0;
+ complexf *out = nullptr;
+ };
+
+ worker_t() {}
+ worker_t(const worker_t& other) = delete;
+ worker_t operator=(const worker_t& other) = delete;
+ worker_t operator=(worker_t&& other) = delete;
+
+ // The move constructor creates a new in_queue and out_queue,
+ // because ThreadsafeQueue is neither copy- nor move-constructible.
+ // Not an issue because creating the workers happens at startup, before
+ // the first work item.
+ worker_t(worker_t&& other) :
+ in_queue(),
+ out_queue(),
+ thread(std::move(other.thread)) {}
+
+ ~worker_t() {
+ if (thread.joinable()) {
+ input_data_t terminate_tag;
+ terminate_tag.terminate = true;
+ in_queue.push(terminate_tag);
+ thread.join();
+ }
+ }
+
+ ThreadsafeQueue<input_data_t> in_queue;
+ ThreadsafeQueue<int> out_queue;
+
+ std::thread thread;
+ };
+
+ std::vector<worker_t> m_workers;
+
+ static void worker_thread(worker_t *workerdata);
+
+ std::vector<float> m_coefs_am; // AM/AM coefficients
+ std::vector<float> m_coefs_pm; // AM/PM coefficients
std::string m_coefs_file;
-
mutable std::mutex m_coefs_mutex;
};
diff --git a/src/OutputUHDFeedback.cpp b/src/OutputUHDFeedback.cpp
index 2a99e6b..b370885 100644
--- a/src/OutputUHDFeedback.cpp
+++ b/src/OutputUHDFeedback.cpp
@@ -43,6 +43,7 @@ DESCRIPTION:
#include <sys/socket.h>
#include <errno.h>
#include <poll.h>
+#include <boost/date_time/posix_time/posix_time.hpp>
#include "OutputUHDFeedback.h"
#include "Utils.h"
@@ -218,152 +219,176 @@ static ssize_t sendall(int socket, const void *buffer, size_t buflen)
return buflen;
}
-void OutputUHDFeedback::ServeFeedbackThread()
+void OutputUHDFeedback::ServeFeedback()
{
- set_thread_name("uhdservefeedback");
+ if ((m_server_sock = socket(PF_INET, SOCK_STREAM, 0)) < 0) {
+ throw std::runtime_error("Can't create TCP socket");
+ }
+
+ struct sockaddr_in addr;
+ addr.sin_family = AF_INET;
+ addr.sin_port = htons(m_port);
+ addr.sin_addr.s_addr = htonl(INADDR_ANY);
+
+ const int reuse = 1;
+ if (setsockopt(m_server_sock, SOL_SOCKET, SO_REUSEADDR, &reuse, sizeof(reuse))
+ < 0) {
+ throw std::runtime_error("Can't reuse address for TCP socket");
+ }
+
+ if (bind(m_server_sock, (struct sockaddr*)&addr, sizeof(addr)) < 0) {
+ close(m_server_sock);
+ throw std::runtime_error("Can't bind TCP socket");
+ }
+
+ if (listen(m_server_sock, 1) < 0) {
+ close(m_server_sock);
+ throw std::runtime_error("Can't listen TCP socket");
+ }
+
+ etiLog.level(info) << "DPD Feedback server listening on port " << m_port;
- try {
- if ((m_server_sock = socket(PF_INET, SOCK_STREAM, 0)) < 0) {
- throw std::runtime_error("Can't create TCP socket");
+ while (m_running) {
+ struct sockaddr_in client;
+ int client_sock = accept_with_timeout(m_server_sock, 1000, &client);
+
+ if (client_sock == -1) {
+ close(m_server_sock);
+ throw runtime_error("Could not establish new connection");
+ }
+ else if (client_sock == -2) {
+ continue;
}
- struct sockaddr_in addr;
- addr.sin_family = AF_INET;
- addr.sin_port = htons(m_port);
- addr.sin_addr.s_addr = htonl(INADDR_ANY);
+ uint8_t request_version = 0;
+ ssize_t read = recv(client_sock, &request_version, 1, 0);
+ if (!read) break; // done reading
+ if (read < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client read request version failed: " << strerror(errno);
+ break;
+ }
- if (bind(m_server_sock, (struct sockaddr*)&addr, sizeof(addr)) < 0) {
- throw std::runtime_error("Can't bind TCP socket");
+ if (request_version != 1) {
+ etiLog.level(info) << "DPD Feedback Server wrong request version";
+ break;
}
- if (listen(m_server_sock, 1) < 0) {
- throw std::runtime_error("Can't listen TCP socket");
+ uint32_t num_samples = 0;
+ read = recv(client_sock, &num_samples, 4, 0);
+ if (!read) break; // done reading
+ if (read < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client read num samples failed";
+ break;
}
- etiLog.level(info) << "DPD Feedback server listening on port " << m_port;
-
- while (m_running) {
- struct sockaddr_in client;
- int client_sock = accept_with_timeout(m_server_sock, 1000, &client);
-
- if (client_sock == -1) {
- throw runtime_error("Could not establish new connection");
- }
- else if (client_sock == -2) {
- continue;
- }
-
- uint8_t request_version = 0;
- ssize_t read = recv(client_sock, &request_version, 1, 0);
- if (!read) break; // done reading
- if (read < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client read request version failed: " << strerror(errno);
- break;
- }
-
- if (request_version != 1) {
- etiLog.level(info) << "DPD Feedback Server wrong request version";
- break;
- }
-
- uint32_t num_samples = 0;
- read = recv(client_sock, &num_samples, 4, 0);
- if (!read) break; // done reading
- if (read < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client read num samples failed";
- break;
- }
-
- // We are ready to issue the request now
- {
- boost::mutex::scoped_lock lock(burstRequest.mutex);
- burstRequest.num_samples = num_samples;
- burstRequest.state = BurstRequestState::SaveTransmitFrame;
-
- lock.unlock();
- }
-
- // Wait for the result to be ready
+ // We are ready to issue the request now
+ {
boost::mutex::scoped_lock lock(burstRequest.mutex);
- while (burstRequest.state != BurstRequestState::Acquired) {
- if (not m_running) break;
- burstRequest.mutex_notification.wait(lock);
- }
+ burstRequest.num_samples = num_samples;
+ burstRequest.state = BurstRequestState::SaveTransmitFrame;
- burstRequest.state = BurstRequestState::None;
lock.unlock();
+ }
+
+ // Wait for the result to be ready
+ boost::mutex::scoped_lock lock(burstRequest.mutex);
+ while (burstRequest.state != BurstRequestState::Acquired) {
+ if (not m_running) break;
+ burstRequest.mutex_notification.wait(lock);
+ }
+
+ 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)));
- burstRequest.num_samples = std::min(burstRequest.num_samples,
- std::min(
- burstRequest.tx_samples.size() / sizeof(complexf),
- burstRequest.rx_samples.size() / sizeof(complexf)));
-
- uint32_t num_samples_32 = burstRequest.num_samples;
- if (sendall(client_sock, &num_samples_32, sizeof(num_samples_32)) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send num_samples failed";
- break;
- }
-
- if (sendall(client_sock,
- &burstRequest.tx_second,
- sizeof(burstRequest.tx_second)) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send tx_second failed";
- break;
- }
-
- if (sendall(client_sock,
- &burstRequest.tx_pps,
- sizeof(burstRequest.tx_pps)) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send tx_pps failed";
- break;
- }
-
- const size_t frame_bytes = burstRequest.num_samples * sizeof(complexf);
-
- assert(burstRequest.tx_samples.size() >= frame_bytes);
- if (sendall(client_sock,
- &burstRequest.tx_samples[0],
- frame_bytes) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send tx_frame failed";
- break;
- }
-
- if (sendall(client_sock,
- &burstRequest.rx_second,
- sizeof(burstRequest.rx_second)) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send rx_second failed";
- break;
- }
-
- if (sendall(client_sock,
- &burstRequest.rx_pps,
- sizeof(burstRequest.rx_pps)) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send rx_pps failed";
- break;
- }
-
- assert(burstRequest.rx_samples.size() >= frame_bytes);
- if (sendall(client_sock,
- &burstRequest.rx_samples[0],
- frame_bytes) < 0) {
- etiLog.level(info) <<
- "DPD Feedback Server Client send rx_frame failed";
- break;
- }
-
- close(client_sock);
+ uint32_t num_samples_32 = burstRequest.num_samples;
+ if (sendall(client_sock, &num_samples_32, sizeof(num_samples_32)) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send num_samples failed";
+ break;
}
+
+ if (sendall(client_sock,
+ &burstRequest.tx_second,
+ sizeof(burstRequest.tx_second)) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send tx_second failed";
+ break;
+ }
+
+ if (sendall(client_sock,
+ &burstRequest.tx_pps,
+ sizeof(burstRequest.tx_pps)) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send tx_pps failed";
+ break;
+ }
+
+ const size_t frame_bytes = burstRequest.num_samples * sizeof(complexf);
+
+ assert(burstRequest.tx_samples.size() >= frame_bytes);
+ if (sendall(client_sock,
+ &burstRequest.tx_samples[0],
+ frame_bytes) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send tx_frame failed";
+ break;
+ }
+
+ if (sendall(client_sock,
+ &burstRequest.rx_second,
+ sizeof(burstRequest.rx_second)) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send rx_second failed";
+ break;
+ }
+
+ if (sendall(client_sock,
+ &burstRequest.rx_pps,
+ sizeof(burstRequest.rx_pps)) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send rx_pps failed";
+ break;
+ }
+
+ assert(burstRequest.rx_samples.size() >= frame_bytes);
+ if (sendall(client_sock,
+ &burstRequest.rx_samples[0],
+ frame_bytes) < 0) {
+ etiLog.level(info) <<
+ "DPD Feedback Server Client send rx_frame failed";
+ break;
+ }
+
+ close(client_sock);
}
- catch (runtime_error &e) {
- etiLog.level(error) << "DPD Feedback Server fault: " << e.what();
+}
+
+void OutputUHDFeedback::ServeFeedbackThread()
+{
+ set_thread_name("uhdservefeedback");
+
+ while (m_running) {
+ try {
+ ServeFeedback();
+ }
+ catch (const runtime_error &e) {
+ etiLog.level(error) << "DPD Feedback Server runtime error: " << e.what();
+ }
+ catch (const std::exception &e) {
+ etiLog.level(error) << "DPD Feedback Server exception: " << e.what();
+ }
+ catch (...) {
+ etiLog.level(error) << "DPD Feedback Server unknown exception!";
+ }
+
+ boost::this_thread::sleep(boost::posix_time::seconds(5));
}
m_running = false;
diff --git a/src/OutputUHDFeedback.h b/src/OutputUHDFeedback.h
index 32668b6..c68f4c2 100644
--- a/src/OutputUHDFeedback.h
+++ b/src/OutputUHDFeedback.h
@@ -101,6 +101,7 @@ class OutputUHDFeedback {
// Thread that listens for requests over TCP to get TX and RX feedback
void ServeFeedbackThread(void);
+ void ServeFeedback(void);
boost::thread rx_burst_thread;
boost::thread burst_tcp_thread;
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 03:11:15 +0000