diff options
| -rw-r--r-- | analyze_aligned_rx_tx.py | 132 | ||||
| -rwxr-xr-x | grc/live_analyse_dab_poly.py | 51 | ||||
| -rw-r--r-- | src/dab_util.py | 28 |
3 files changed, 190 insertions, 21 deletions
diff --git a/analyze_aligned_rx_tx.py b/analyze_aligned_rx_tx.py new file mode 100644 index 0000000..2b9f70c --- /dev/null +++ b/analyze_aligned_rx_tx.py @@ -0,0 +1,132 @@ +""" +Generate analytic plots from aligned TX/RX recordings. +""" + +import argparse +import numpy as np +import matplotlib.pyplot as plt +import src.dab_util as du +import scipy.stats as st +from sklearn.neighbors import KernelDensity +import os + +SCATTER_SIZE = 0.0001 + +def kde2D(x, y, bandwidth, xbins=256j, ybins=256j, **kwargs): + xx, yy = np.mgrid[x.min():x.max():xbins, + y.min():y.max():ybins] + + xy_sample = np.vstack([yy.ravel(), xx.ravel()]).T + xy_train = np.vstack([y, x]).T + + kde_skl = KernelDensity(bandwidth=bandwidth, **kwargs) + kde_skl.fit(xy_train) + + z = np.exp(kde_skl.score_samples(xy_sample)) + return xx, yy, np.reshape(z, xx.shape) + +def plot_density(x, y, scatter=False, path=None): + x = np.abs(x) + y = np.abs(y) + + x = x - np.min(x) + x = x / np.max(x) + y = y - np.min(y) + y = y / np.max(y) + + max_val = max(np.max(x), np.max(y)) + min_val = max(np.min(x), np.min(y)) + + xx, yy, zz = kde2D(x, y, (max_val - min_val) * 0.01) + plt.pcolormesh(xx, yy, zz) + + plt.xlabel("Normalized Absolute TX Amplitude") + plt.ylabel("Normalized Absolute RX Amplitude") + + plt.savefig(path) + +def scatter(x_pos, y_pos, x_neg, y_neg, path=None): + x_pos, y_pos, x_neg, y_neg = np.abs(x_pos), np.abs(y_pos), np.abs(x_neg), np.abs(y_neg) + plt.scatter(x_pos, y_pos, s=SCATTER_SIZE, facecolor='blue', label="Positive TX/RX") + plt.scatter(x_neg, y_neg, s=SCATTER_SIZE, facecolor='red', label="Negative TX/RX") + plt.xlabel("Absolute TX Amplitude") + plt.ylabel("Absolute RX Amplitude") + plt.legend() + plt.savefig(path) + plt.clf() + +def scatter_phase(x_pos, y_pos, x_neg, y_neg, path=None): + x_pos_abs = np.abs(x_pos) + x_neg_abs = np.abs(x_neg) + + phase_diff_pos = np.angle(x_pos, deg=True) - np.angle(y_pos, deg=True) + phase_diff_neg = np.angle(x_neg, deg=True) - np.angle(y_neg, deg=True) + + phase_diff_pos = np.mod(phase_diff_pos, 180) + phase_diff_neg = np.mod(phase_diff_neg, 180) + + plt.scatter(x_pos_abs, phase_diff_pos, s=SCATTER_SIZE, facecolor='blue', label="Positive TX/RX") + plt.scatter(x_neg_abs, phase_diff_neg, s=SCATTER_SIZE, facecolor='red', label="Negative TX/TX") + + plt.ylabel("Phase difference") + plt.xlabel("Absolute Amplitude") + plt.legend() + + plt.savefig(path) + plt.clf() + +def plot_time(rx_rec, tx_rec, path=None, samples=256): + plt.plot(np.angle(rx_rec[:256]), c='blue', label="RX") + plt.plot(np.angle(tx_rec[:256]), c='red', label="TX") + plt.ylabel("Phase") + plt.xlabel("Sample") + plt.savefig(path) + plt.clf() + +def main(): + if not os.path.isdir(FLAGS.out_dir): + os.makedirs(FLAGS.out_dir) + + rx_rec = du.fromfile(filename=FLAGS.rx_path) + tx_rec = du.fromfile(filename=FLAGS.tx_path) + + sel_pos = (rx_rec > 0) & (tx_rec > 0) + rx_rec_pos = rx_rec[sel_pos] + tx_rec_pos = tx_rec[sel_pos] + + + sel_pos = (rx_rec < 0) & (tx_rec < 0) + rx_rec_neg = rx_rec[sel_pos] + tx_rec_neg = tx_rec[sel_pos] + + scatter(tx_rec_pos, rx_rec_pos, tx_rec_neg, rx_rec_neg, path=FLAGS.out_dir + '/am_am_scatter.pdf') + scatter_phase(tx_rec_pos, rx_rec_pos, tx_rec_neg, rx_rec_neg, path=FLAGS.out_dir + '/am_pm_pos_scatter.pdf') + + plot_time(rx_rec, tx_rec, path=FLAGS.out_dir + '/phase_over_time.pdf') + + plot_density(tx_rec_pos, rx_rec_pos, path=FLAGS.out_dir + '/am_am_pos.pdf') + plot_density(tx_rec_neg, rx_rec_neg, path=FLAGS.out_dir + '/am_am_neg.pdf') + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument( + '--rx_path', + type=str, + default='/tmp/record/2_rx_record.iq', + help="Path to complex64 rx recording" + ) + parser.add_argument( + '--tx_path', + type=str, + default='/tmp/record/2_tx_record.iq', + help="Path to complex64 tx recording" + ) + parser.add_argument( + '--out_dir', + type=str, + default='/tmp/analyze_aligned_rx_tx', + help="Output path" + ) + FLAGS, unparsed = parser.parse_known_args() + + main() diff --git a/grc/live_analyse_dab_poly.py b/grc/live_analyse_dab_poly.py index 126c43d..029b541 100755 --- a/grc/live_analyse_dab_poly.py +++ b/grc/live_analyse_dab_poly.py @@ -3,7 +3,7 @@ ################################################## # GNU Radio Python Flow Graph # Title: Live Analyse Dab Poly -# Generated: Tue May 23 20:21:44 2017 +# Generated: Wed Jun 7 21:22:24 2017 ################################################## if __name__ == '__main__': @@ -60,12 +60,13 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): ################################################## # Variables ################################################## + self.samp_rate = samp_rate = 8e6 self.txgain = txgain = 80 self.shift_freq = shift_freq = 1 - self.samp_rate_2 = samp_rate_2 = 2048000 - self.samp_rate_1 = samp_rate_1 = 8192000 self.rxgain = rxgain = 10 self.freq = freq = 222e6 + self.f2 = f2 = samp_rate / 3.875 + self.f1 = f1 = samp_rate / 4 self.ampl = ampl = 0.4 self.a_8 = a_8 = 0 self.a_7 = a_7 = 0 @@ -121,7 +122,7 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): channels=range(1), ), ) - self.uhd_usrp_source_0.set_samp_rate(samp_rate_2) + self.uhd_usrp_source_0.set_samp_rate(samp_rate) self.uhd_usrp_source_0.set_center_freq(freq, 0) self.uhd_usrp_source_0.set_gain(rxgain, 0) self.uhd_usrp_sink_0 = uhd.usrp_sink( @@ -131,7 +132,7 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): channels=range(1), ), ) - self.uhd_usrp_sink_0.set_samp_rate(samp_rate_1) + self.uhd_usrp_sink_0.set_samp_rate(samp_rate) self.uhd_usrp_sink_0.set_center_freq(freq, 0) self.uhd_usrp_sink_0.set_gain(txgain, 0) self.uhd_amsg_source_0 = uhd.amsg_source(device_addr="", msgq=uhd_amsg_source_0_msgq_out) @@ -142,7 +143,7 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): 16000, #size firdes.WIN_BLACKMAN_hARRIS, #wintype 0, #fc - samp_rate_2, #bw + samp_rate, #bw "", #name 1 #number of inputs ) @@ -205,6 +206,17 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): event.accept() + def get_samp_rate(self): + return self.samp_rate + + def set_samp_rate(self, samp_rate): + self.samp_rate = samp_rate + self.set_f1(self.samp_rate / 4) + self.set_f2(self.samp_rate / 3.875) + self.qtgui_freq_sink_x_0_0.set_frequency_range(0, self.samp_rate) + self.uhd_usrp_sink_0.set_samp_rate(self.samp_rate) + self.uhd_usrp_source_0.set_samp_rate(self.samp_rate) + def get_txgain(self): return self.txgain @@ -219,21 +231,6 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): def set_shift_freq(self, shift_freq): self.shift_freq = shift_freq - def get_samp_rate_2(self): - return self.samp_rate_2 - - def set_samp_rate_2(self, samp_rate_2): - self.samp_rate_2 = samp_rate_2 - self.qtgui_freq_sink_x_0_0.set_frequency_range(0, self.samp_rate_2) - self.uhd_usrp_source_0.set_samp_rate(self.samp_rate_2) - - def get_samp_rate_1(self): - return self.samp_rate_1 - - def set_samp_rate_1(self, samp_rate_1): - self.samp_rate_1 = samp_rate_1 - self.uhd_usrp_sink_0.set_samp_rate(self.samp_rate_1) - def get_rxgain(self): return self.rxgain @@ -250,6 +247,18 @@ class live_analyse_dab_poly(gr.top_block, Qt.QWidget): self.uhd_usrp_sink_0.set_center_freq(self.freq, 0) self.uhd_usrp_source_0.set_center_freq(self.freq, 0) + def get_f2(self): + return self.f2 + + def set_f2(self, f2): + self.f2 = f2 + + def get_f1(self): + return self.f1 + + def set_f1(self, f1): + self.f1 = f1 + def get_ampl(self): return self.ampl diff --git a/src/dab_util.py b/src/dab_util.py index 3187036..e33ea05 100644 --- a/src/dab_util.py +++ b/src/dab_util.py @@ -77,6 +77,34 @@ def lag_upsampling(sig_orig, sig_rec, n_up): l_orig = float(l) / n_up return l_orig +def subsample_align_upsampling(sig1, sig2, n_up=32): + """ + Returns an aligned version of sig1 and sig2 by cropping and subsample alignment + Using upsampling + """ + assert(sig1.shape[0] == sig2.shape[0]) + + if sig1.shape[0] % 2 == 1: + sig1 = sig1[:-1] + sig2 = sig2[:-1] + + sig1_up = signal.resample(sig1, sig1.shape[0] * n_up) + sig2_up = signal.resample(sig2, sig2.shape[0] * n_up) + + off_meas = 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] + + sig1 = signal.resample(sig1_up, sig1_up.shape[0] / n_up).astype(np.complex64) + sig2 = signal.resample(sig2_up, sig2_up.shape[0] / n_up).astype(np.complex64) + return sig1, sig2 + def subsample_align(sig1, sig2): """ Returns an aligned version of sig1 and sig2 by cropping and subsample alignment |