Add TX constellation plot to show_spectrum

1 files changed, 95 insertions, 12 deletions

diff --git a/dpd/show_spectrum.py b/dpd/show_spectrum.py
index e92c1d0..0ed0ebf 100755
--- a/dpd/show_spectrum.py
+++ b/dpd/show_spectrum.py
@@ -21,9 +21,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 +40,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 = []
@@ -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,77 @@ 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)])
+    #imsave("spectrums.png", np.abs(spectrums))
+
+    # Only take bins that are supposed to contain energy
+    spectrums = np.concatenate([spectrums[...,256:1024], spectrums[...,1025:1793]], axis=1)
+
+    sym_indices = (np.tile(np.arange(num_syms-1).reshape(num_syms-1,1), (1,NbCarriers)) +
+                   np.tile(np.linspace(-0.25, 0.25, 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 +213,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")