1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
|
# -*- coding: utf-8 -*-
#
# DPD Computation Engine, Capture TX signal and RX feedback using ODR-DabMod's
# DPD Server.
#
# Copyright (c) 2017 Andreas Steger
# Copyright (c) 2018 Matthias P. Braendli
#
# http://www.opendigitalradio.org
#
# This file is part of ODR-DabMod.
#
# ODR-DabMod is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# ODR-DabMod is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ODR-DabMod. If not, see <http://www.gnu.org/licenses/>.
import socket
import struct
import os.path
import logging
import numpy as np
from scipy import signal
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import io
from . import Align as sa
def correlation_coefficient(sig_tx, sig_rx):
return np.corrcoef(sig_tx, sig_rx)[0, 1]
def align_samples(sig_tx, sig_rx):
"""
Returns an aligned version of sig_tx and sig_rx by cropping, subsample alignment and
correct phase offset
"""
# Coarse sample-level alignment
c = np.abs(signal.correlate(sig_rx, sig_tx))
off_meas = np.argmax(c) - sig_tx.shape[0] + 1
off = int(abs(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]
# Fine subsample alignment and phase offset
sig_rx = sa.subsample_align(sig_rx, sig_tx)
sig_rx = sa.phase_align(sig_rx, sig_tx)
return sig_tx, sig_rx, abs(off_meas)
class Capture:
"""Capture samples from ODR-DabMod"""
def __init__(self, samplerate, port, num_samples_to_request, plot_dir):
self.samplerate = samplerate
self.sizeof_sample = 8 # complex floats
self.port = port
self.num_samples_to_request = num_samples_to_request
self.plot_dir = plot_dir
# Before we run the samples through the model, we want to accumulate
# them into bins depending on their amplitude, and keep only n_per_bin
# samples to avoid that the polynomial gets overfitted in the low-amplitude
# part, which is less interesting than the high-amplitude part, where
# non-linearities become apparent.
self.binning_n_bins = 64 # Number of bins between binning_start and binning_end
self.binning_n_per_bin = 128 # Number of measurements pre bin
self.rx_normalisation = 1.0
self.clear_accumulated()
def clear_accumulated(self):
self.binning_start = 0.0
self.binning_end = 1.0
# axis 0: bins
# axis 1: 0=tx, 1=rx
self.accumulated_bins = [np.zeros((0, 2), dtype=np.complex64) for i in range(self.binning_n_bins)]
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.settimeout(4)
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 {}, max {}, median {}'.format(
np.min(np.abs(txframe)),
np.max(np.abs(txframe)),
np.median(np.abs(txframe))))
logging.debug('rxframe: min {}, max {}, median {}'.format(
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 _plot_spectrum(self, signal, filename, title):
fig = plt.figure()
ax = plt.subplot(1, 1, 1)
fft = np.fft.fftshift(np.fft.fft(signal))
fft_db = 20 * np.log10(np.abs(fft))
ax.plot(fft_db)
ax.set_title(title)
fig.tight_layout()
fig.savefig(os.path.join(self.plot_dir, filename))
plt.close(fig)
def calibrate(self):
txframe, tx_ts, rxframe, rx_ts = self.receive_tcp()
# Normalize received signal with sent signal
tx_median = np.median(np.abs(txframe))
rx_median = np.median(np.abs(rxframe))
self.rx_normalisation = tx_median / rx_median
rxframe = rxframe * self.rx_normalisation
txframe_aligned, rxframe_aligned, coarse_offset = align_samples(txframe, rxframe)
self._plot_spectrum(rxframe[:8192], "rxframe.png", "RX Frame")
self._plot_spectrum(txframe[:8192], "txframe.png", "RX Frame")
return tx_ts, tx_median, rx_ts, rx_median, np.abs(coarse_offset), correlation_coefficient(txframe_aligned, rxframe_aligned)
def get_samples(self):
"""Connect to ODR-DabMod, retrieve TX and RX samples, load
into numpy arrays, and return a tuple
(txframe_aligned, tx_ts, tx_median, rxframe_aligned, rx_ts, rx_median)
"""
txframe, tx_ts, rxframe, rx_ts = self.receive_tcp()
# Normalize received signal with calibrated normalisation
rxframe = rxframe * self.rx_normalisation
txframe_aligned, rxframe_aligned, coarse_offset = align_samples(txframe, rxframe)
self._bin_and_accumulate(txframe_aligned, rxframe_aligned)
return txframe_aligned, tx_ts, tx_median, rxframe_aligned, rx_ts, rx_median
def bin_histogram(self):
return [b.shape[0] for b in self.accumulated_bins]
def pointcloud_png(self):
fig = plt.figure()
ax = plt.subplot(1, 1, 1)
for b in self.accumulated_bins:
if b:
ax.scatter(
np.abs(b[0]),
np.abs(b[1]),
s=0.1,
color="black")
ax.set_title("Captured and Binned Samples")
ax.set_xlabel("TX Amplitude")
ax.set_ylabel("RX Amplitude")
ax.set_ylim(0, 0.8)
ax.set_xlim(0, 1.1)
ax.legend(loc=4)
fig.tight_layout()
fig.savefig(os.path.join(self.plot_dir, "pointcloud.png"))
plt.close(fig)
def _bin_and_accumulate(self, txframe, rxframe):
"""Bin the samples and extend the accumulated samples"""
bin_edges = np.linspace(self.binning_start, self.binning_end, self.binning_n_bins)
minsize = self.num_samples_to_request
for i, (tx_start, tx_end) in enumerate(zip(bin_edges, bin_edges[1:])):
txframe_abs = np.abs(txframe)
indices = np.bitwise_and(tx_start < txframe_abs, txframe_abs <= tx_end)
txsamples = np.asmatrix(txframe[indices])
rxsamples = np.asmatrix(rxframe[indices])
binned_sample_pairs = np.concatenate((txsamples, rxsamples)).T
missing_in_bin = self.binning_n_per_bin - self.accumulated_bins[i].shape[0]
num_to_append = min(missing_in_bin, binned_sample_pairs.shape[0])
print("Handling bin {} {}-{}, {} available, {} missing".format(i, tx_start, tx_end, binned_sample_pairs.shape[0], missing_in_bin))
if num_to_append:
print("Appending {} to bin {} with shape {}".format(num_to_append, i, self.accumulated_bins[i].shape))
self.accumulated_bins[i] = np.concatenate((self.accumulated_bins[i], binned_sample_pairs[:num_to_append,...]))
print("{} now has shape {}".format(i, self.accumulated_bins[i].shape))
|
