diff options
Diffstat (limited to 'dpd/show_spectrum.py')
| -rwxr-xr-x | dpd/show_spectrum.py | 132 |
1 files changed, 114 insertions, 18 deletions
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") |