#!/usr/bin/env python """ Plot jamming mitigation data. Either processes a single T04 file or else all of the T04 files in a specified directory using wildcard. This loads rec35:50 ### If you have very old data using the obsolete ### ### rec35:30 then you'll need to use an older ### ### version of this script -- the last version to ### ### support rec35:30 was committed on 05/21/2025. ### Example: ./plt_jamming_mitigation.py /home/ddelore/DataDir/RX6-2209 Optional arguments: -f, --fft_blanking Plot FFT bins blanked vs. time --rt_band RT band number [defaults to #0] --ant_num Antenna number [defaults to #0] --start_time Start time for plotting [defaults to -inf] --stop_time Stop time for plotting [defaults to +inf] --save_figs Save to .png [defaults to OFF] Remember your rt_band mapping: band rf_band rt_band L1 0 0 L2 1 1 L5 6 4 E6 7 5 B1 8 6 S1 10 11 """ import os import argparse import numpy as np import matplotlib.pyplot as plt from pylab import subplots, where, log10, meshgrid, r_, show from mutils import vd2cls, rt_band_to_label T0X_MAXSTATES = 12 T0X_N_FIR = 6 def rt_band_to_center_freq(rt_band): switch = \ { 0: 1590, # LG1 #0: 1583, # LGB1 #0: 1585, # Argon 1: 1240, # LG2 4: 1190, # LE5 5: 1280, # E6 6: 1550, # B1 11: 2503 # S-band } return switch.get(rt_band,"Invalid rt_band!!") def plt_jamming_mitigation(args): path_name = args.path_name # Grab all T04 files in a directory if not path_name.endswith('.T04'): path_name = os.path.join(path_name, '*.T04') # Load rec35:50 d = vd2cls(path_name, rec='-d35:50') print('\n Rec35:sub50 version = ' + str(int(d[0,d.k['VERSION']])) + '\n') version = int(d[0,d.k['VERSION']]) n_bands = int(d[0,d.k['MAX_BANDS']]) n_comm = int([x for x in d.k.keys() if x.startswith('COMM_ATTEN_')][-1].split('_')[-1])+1 # Get the band center frequency if version >= 4: for i in range(n_bands): rt_band_idx = d.k['NF_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] if np.isfinite(rt_band) and int(rt_band) == args.rt_band: idx = d.k['NF_CENTER_FREQ_%2.2d' %i] center_freq_mhz = int(d[0,idx])/1e6 break else: center_freq_mhz = int(rt_band_to_center_freq(args.rt_band)) # Set the start & stop times for plotting if args.start_time >= 0: start_x = args.start_time else: start_x = d[0].GPS_SEC if args.stop_time >= 0: stop_x = args.stop_time else: stop_x = d[-1].GPS_SEC ################################################################### # # RSSI & Noise Power # ################################################################### # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['NF_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['NF_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots(2,1,sharex=True) ax[0].set_title(path_name) ax[0].set_xlabel('GPS TOW [sec]') ax[0].set_ylabel('RSSI [dBm]') ax[1].set_xlabel('GPS TOW [sec]') ax[1].set_ylabel('Noise Pwr [dB]') # Channel RSSI idx = d.k['CHAN_RSSI_%2.2d' %i] ax[0].plot( d[:].GPS_SEC, d[:,idx], 'b', label=rt_band_to_label(rt_band) ) # Loop over common-path entries for j in range(n_comm): if version >= 2: rt_band_idx = d.k['RT_BAND_%2.2d' %j] ant_num_idx = d.k['ANT_NUM_%2.2d' %j] # Is this L-band -- chan RT band != 11 but common RT band != 0 (L-band) if args.rt_band >= 0 and rt_band != 11 and version < 2: continue elif args.rt_band >= 0 and rt_band != 11 and int(d[0,rt_band_idx]) != 0: continue # Is this S-band -- chan RT band == 11 but common RT band != 11 (S-band) if args.rt_band >= 0 and rt_band == 11 and version < 2: continue elif args.rt_band >= 0 and rt_band == 11 and int(d[0,rt_band_idx]) != 11: continue # Is this the requested antenna if ant_num != int(d[0,ant_num_idx]): continue # Common RSSI idx = d.k['COMM_RSSI_%2.2d' %j] ax[0].plot( d[:].GPS_SEC, d[:,idx], 'r', label='common', linewidth=3 ) ax[0].grid() ax[0].legend(loc='best') ax[0].set_xlim(start_x,stop_x) ax[0].set_ylim(-60,+5) # Noise power idx = d.k['NF_NOISE_PWR_dB_%2.2d' %i] ax[1].plot( d[:].GPS_SEC, d[:,idx], 'b', label=rt_band_to_label(rt_band) ) ax[1].grid() ax[1].legend(loc='best') ax[1].set_xlim(start_x,stop_x) ax[1].set_ylim(min(20,ax[1].get_ylim()[0]), max(40,ax[1].get_ylim()[1])) ################################################################### # # AGC & Attenuation # ################################################################### # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['NF_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['NF_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots(3,1,sharex=True) ax[0].set_title(path_name) ax[0].set_xlabel('GPS TOW [sec]') ax[0].set_ylabel('Control Mode') ax[1].set_xlabel('GPS TOW [sec]') ax[1].set_ylabel('AGC [PWM count]') ax[2].set_xlabel('GPS TOW [sec]') ax[2].set_ylabel('Attenuation [dB]') # AGC Control Mode idx = d.k['AGC_CTRL_MODE_%2.2d' %i] ax[0].plot( d[:].GPS_SEC, d[:,idx], 'r', label=rt_band_to_label(rt_band) + ' ctrl mode' ) ax[0].grid() ax[0].legend(loc='best') ax[0].set_xlim(start_x,stop_x) ax[0].set_ylim(-0.5,+4.5) # AGC Setting idx = d.k['AGC_PWM_%2.2d' %i] ax[1].plot( d[:].GPS_SEC, d[:,idx], 'b', label=rt_band_to_label(rt_band) + ' agc pwm' ) ax[1].grid() ax[1].legend(loc='best') ax[1].set_xlim(start_x,stop_x) ax[1].set_ylim(min(200,ax[1].get_ylim()[0]), max(450,ax[1].get_ylim()[1])) # Channel attenuation idx = d.k['CHAN_ATTEN_%2.2d' %i] ax[2].plot( d[:].GPS_SEC, d[:,idx], 'b', label=rt_band_to_label(rt_band) ) # Loop over common-path entries for j in range(n_comm): if version >= 2: rt_band_idx = d.k['RT_BAND_%2.2d' %j] ant_num_idx = d.k['ANT_NUM_%2.2d' %j] # Is this L-band -- chan RT band != 11 but common RT band != 0 (L-band) if args.rt_band >= 0 and rt_band != 11 and version < 2: continue elif args.rt_band >= 0 and rt_band != 11 and int(d[0,rt_band_idx]) != 0: continue # Is this S-band -- chan RT band == 11 but common RT band != 11 (S-band) if args.rt_band >= 0 and rt_band == 11 and version < 2: continue elif args.rt_band >= 0 and rt_band == 11 and int(d[0,rt_band_idx]) != 11: continue # Is this the requested antenna if ant_num != int(d[0,ant_num_idx]): continue # Common attenuation idx = d.k['COMM_ATTEN_%2.2d' %j] ax[2].plot( d[:].GPS_SEC, d[:,idx], 'r', label='common', linewidth=3 ) ax[2].grid() ax[2].legend(loc='best') ax[2].set_xlim(start_x,stop_x) ax[2].set_ylim(-5,35) ################################################################### # # RFI States # # Note: The order that fields are written to a T04 file # means that the second antenna (if present) in the # write order will have the RFI states all jumbled up # in indexing. Thus the plot colors for frequency and # J/N will look all mixed up, and it's too much bother # to try to sort plot colors by ID -- just know that # under the hood the RFI states are persistent. # ################################################################### # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['NF_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['NF_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots(2,1,sharex=True) ax[0].set_title(path_name + ' - ' + rt_band_to_label(rt_band)) ax[0].set_xlabel('GPS TOW [sec]') ax[0].set_ylabel('RFI States [MHz]') ax[1].set_xlabel('GPS TOW [sec]') ax[1].set_ylabel('RFI J/N [dB]') # RFI Frequencies for j in range(n_bands * T0X_MAXSTATES): rfi_band_idx = d.k['RFI_RT_BAND_%3.3d' %j] rfi_ant_idx = d.k['RFI_ANT_NUM_%3.3d' %j] idx = where( (d[:,rfi_band_idx]==rt_band) & (d[:,rfi_ant_idx]==ant_num) )[0] if len(idx) > 0: rfi_freq_idx = d.k['RFI_FREQ_%3.3d' %j] ax[0].plot( d[idx].GPS_SEC, d[idx,rfi_freq_idx]*1e-6, '.' ) ax[0].grid() ax[0].set_xlim(start_x,stop_x) ax[0].set_ylim(center_freq_mhz-25,center_freq_mhz+25) # RFI J/N [dB] for j in range(n_bands * T0X_MAXSTATES): rfi_band_idx = d.k['RFI_RT_BAND_%3.3d' %j] rfi_ant_idx = d.k['RFI_ANT_NUM_%3.3d' %j] idx = where( (d[:,rfi_band_idx]==rt_band) & (d[:,rfi_ant_idx]==ant_num) )[0] if len(idx) > 0: rfi_j2n_idx = d.k['RFI_J2N_FILT_%3.3d' %j] ax[1].plot( d[idx].GPS_SEC, d[idx,rfi_j2n_idx], '.' ) ax[1].grid() ax[1].set_xlim(start_x,stop_x) ax[1].set_ylim(-5,max(20,ax[1].get_ylim()[1])) ################################################################### # # Notch Filters # # Note: The order that fields are written to a T04 file # means that the second antenna (if present) in the # write order will have the notch filters all jumbled # in indexing. Thus the plot colors for frequency will # look all mixed up, and it's too much bother to try to # sort plot colors by ID -- just know that under the hood # the notch filters are persistent. # ################################################################### # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['NF_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['NF_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots() ax.set_title(path_name + ' - ' + rt_band_to_label(rt_band)) ax.set_xlabel('GPS TOW [sec]') ax.set_ylabel('Notch Filter Placement [MHz]') # Notch Filters for j in range(n_bands * T0X_N_FIR): fir_band_idx = d.k['HW_ASSIST_RT_BAND_%3.3d' %j] fir_ant_idx = d.k['HW_ASSIST_ANT_NUM_%3.3d' %j] idx = where( (d[:,fir_band_idx]==rt_band) & (d[:,fir_ant_idx]==ant_num) )[0] if len(idx) > 0: alloc_idx = d.k['HW_ASSIST_ALLOC_BY_HW_%3.3d' %j] cmd_idx = d.k['HW_ASSIST_CMD_BY_SW_%3.3d' %j] fir_freq_idx = d.k['HW_ASSIST_FREQ_%3.3d' %j] ax.plot( d[idx].GPS_SEC, (d[idx,fir_freq_idx]*1e-6)*d[idx,cmd_idx], '.' ) ax.plot( d[idx].GPS_SEC, (d[idx,fir_freq_idx]*1e-6)*d[idx,alloc_idx], 'rx' ) ax.grid() ax.set_xlim(start_x,stop_x) ax.set_ylim(center_freq_mhz-25,center_freq_mhz+25) if args.save_figs: plt.savefig(os.path.splitext(path_name)[0] + "_notch_filters.png",dpi=150) ################################################################### # # FFT Filter # ################################################################### # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['FFT_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['FFT_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots(3,1,sharex=True) ax[0].set_title(path_name + ' - ' + rt_band_to_label(rt_band)) ax[0].set_xlabel('GPS TOW [sec]') ax[0].set_ylabel('FFT Blanking Fraction') #ax[0].set_ylabel('FFT Epochs') ax[1].set_xlabel('GPS TOW [sec]') ax[1].set_ylabel('FFT Blanking Threshold') ax[2].set_xlabel('GPS TOW [sec]') ax[2].set_ylabel('Override is Active') # FFT Blanking Fractions a_idx = d.k['FFT_N_BLANKED_BY_HW_%2.2d' %i] b_idx = d.k['FFT_N_GT_THRESHOLD_SW_%2.2d' %i] ax[0].plot( d[:].GPS_SEC, d[:,a_idx]/128, 'b', label='blanked by HW' ) #ax[0].plot( d[:].GPS_SEC, d[:,b_idx]/128, 'r', label='SW predicted' ) ax[0].grid() ax[0].legend(loc='best') ax[0].set_xlim(start_x,stop_x) ax[0].set_ylim(-0.10,+1.10) # FFT Blanking Thresholds candidate_1_idx = d.k['FFT_CANDIDATE_1_%2.2d' %i] candidate_2_idx = d.k['FFT_CANDIDATE_2_%2.2d' %i] thresh_1 = 10*log10(d[:,candidate_1_idx]) thresh_2 = 10*log10(d[:,candidate_2_idx]) ax[1].plot( d[:].GPS_SEC, thresh_1, 'b', label='sigma-based' ) ax[1].plot( d[:].GPS_SEC, thresh_2, 'r', label='median-based' ) ax[1].grid() ax[1].legend(loc='best') ax[1].set_xlim(start_x,stop_x) # FFT Over-ride States d0_idx = d.k['FFT_OVERRIDE_IS_ACTIVE_0_%2.2d' %i] d1_idx = d.k['FFT_OVERRIDE_IS_ACTIVE_1_%2.2d' %i] d2_idx = d.k['FFT_OVERRIDE_IS_ACTIVE_2_%2.2d' %i] if version >= 6: d3_idx = d.k['FFT_OVERRIDE_IS_ACTIVE_3_%2.2d' %i] ax[2].plot( d[:].GPS_SEC, d[:,d0_idx], label='override #0 is active', linewidth=3 ) ax[2].plot( d[:].GPS_SEC, d[:,d1_idx], label='override #1 is active', linewidth=3 ) ax[2].plot( d[:].GPS_SEC, d[:,d2_idx], label='override #2 is active', linewidth=3 ) if version >= 6: ax[2].plot( d[:].GPS_SEC, d[:,d3_idx], label='override #3 is active', linewidth=3 ) ax[2].grid() ax[2].legend(loc='best') ax[2].set_xlim(start_x,stop_x) ax[2].set_ylim(-0.10,+1.10) ################################################################### # # HW Actual FFT Blanking # ################################################################### # Plot FFT bins blanked vs. time if args.fft_blanking: # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['FFT_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['FFT_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots() ax.set_title(path_name + ' - ' + rt_band_to_label(rt_band)) ax.set_xlabel('GPS TOW [sec]') ax.set_ylabel('HW FFT Filter Blanking [MHz]') freq_mhz = (r_[0:128]/128. - 0.5)*50 + center_freq_mhz X,Y = meshgrid(freq_mhz, d[:].GPS_SEC) z = np.full_like(X,0) for j in range(0,z.shape[0]): for b in range(0,8): bitfield_idx = d.k['HW_BINS_BLANKED_%2.2d_%2.2d' %(i,b)] bitfield = int(d[j,bitfield_idx]) for n in range(0,16): z[j,16*b+n] = 1 & (bitfield >> (15-n)) ax.pcolormesh(Y,X,z) ax.set_xlim(start_x,stop_x) ax.set_ylim(center_freq_mhz-25,center_freq_mhz+25) if args.save_figs: plt.savefig(os.path.splitext(path_name)[0] + "_fft_blanking.png",dpi=150) ################################################################### # # SW Predicted FFT Blanking # ################################################################### # Plot FFT bins blanked vs. time if args.fft_blanking: # Loop over all data entries for i in range(n_bands): # Get the RT band for this entry rt_band_idx = d.k['FFT_RT_BAND_%2.2d' %i] rt_band = d[0,rt_band_idx] # Get the antenna number for this entry ant_num_idx = d.k['FFT_ANT_NUM_%2.2d' %i] ant_num = d[0,ant_num_idx] # Is this the requested RT band if np.isnan(rt_band) or int(rt_band) != args.rt_band: continue # Is this the requested antenna if np.isnan(ant_num) or int(ant_num) != args.ant_num: continue # Cast to integer rt_band = int(rt_band) ant_num = int(ant_num) # Open a new figure window fig,ax = subplots() ax.set_title(path_name + ' - ' + rt_band_to_label(rt_band)) ax.set_xlabel('GPS TOW [sec]') ax.set_ylabel('SW FFT Filter Blanking [MHz]') freq_mhz = (r_[0:128]/128. - 0.5)*50 + center_freq_mhz X,Y = meshgrid(freq_mhz, d[:].GPS_SEC) z = np.full_like(X,0) for j in range(0,z.shape[0]): for b in range(0,8): bitfield_idx = d.k['SW_BINS_BLANKED_%2.2d_%2.2d' %(i,b)] bitfield = int(d[j,bitfield_idx]) for n in range(0,16): z[j,16*b+n] = 1 & (bitfield >> (15-n)) ax.pcolormesh(Y,X,z) ax.set_xlim(start_x,stop_x) ax.set_ylim(center_freq_mhz-25,center_freq_mhz+25) ################################################################### # # All done, refresh the figure windows # ################################################################### show() ####################################################################### # # Command line arguments # ####################################################################### if __name__ == '__main__': parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, description=__doc__) parser.add_argument('path_name', help='name of or path to your T04 file(s)') parser.add_argument('-f', '--fft_blanking', action='store_true', help='(optional) plot FFT bins blanked vs. time') parser.add_argument('--rt_band', type=int, default=0, help='RT band number [defaults to LG1]') parser.add_argument('--ant_num', type=int, default=0, help='antenna number [defaults to #0]') parser.add_argument('--start_time', type=int, default=-1, help='start time for plotting [defaults to -inf]') parser.add_argument('--stop_time', type=int, default=-1, help='stop time for plotting [defaults to +inf]') parser.add_argument('--save_figs', action='store_true', help='(optional) save to .png [defaults to OFF]') args = parser.parse_args() plt_jamming_mitigation(args)