import matplotlib #import scipy.stats # Allow running headless from the command line matplotlib.use("agg") import pandas as pd from pylab import * import numpy as np import os import argparse import plotly # # Copyright Trimble Inc # # Provided to Geotronics Poland to demonstrate accessing Trimble receiver FFT data # # def ecdf(data): # Compute ECDF x = np.sort(data) n = x.size y = np.arange(1, n+1) / n return(x,y) # FFT - spectragram FFT # freqData - frequency of each FFT point in MHz # band - ASCII string for the band def analyzeFFT(FFT, FFTtime, freqData, band, fileBase): # Here is an example of a very crude way to analyze the FFT data and look for jammers. # In reality, a much more sophisticated detector will be needed! if(band == 'L1'): # For GPS, Galileo, and BeiDou B1C, the most critical frequencies # are around 1575.42MHz. Review what is happening +/-2MHz idx = np.where((freqData >= 1573.42) & (freqData <= 1576.42))[0] # Drop all data into a single array powerData = FFT[:,idx].flatten() # Find the median of the data. If the jammer is really powerful, it might have # changed the noise floor. medianPower = np.median(powerData) print(medianPower) # Adjust by the median powerData -= medianPower # Find any data > 15dB above the noise floor # We won't normally get any data above this threshold powerIndex = np.where(powerData > 15.0)[0] print('Number of points > 15dB above the noise floor: %d Percent %.4f%%' % (len(powerIndex), len(powerIndex)/len(powerData)*100)) # Plot powerData as a CDF fig, axs = plt.subplots(1, 1) # Matplotlib >= 3.8 has this #axs.ecdf(powerData, label='CDF') (x,y) = ecdf(powerData) axs.plot(x, y, label='CDF', color='blue') axs.set_title(fileBase + ' Power L1 +/- 2MHz CDF') axs.set_xlabel('Power [dB]') axs.set_ylabel('CDF') axs.grid(True) axs.set_xlim(-20,20) axs.set_ylim(0,1) axs.tick_params(axis='both', which='major', labelsize=6) plt.tight_layout() fig.savefig(fileBase + '-Jam-CDF.png',dpi=600) plt.close() FFT[:,idx] -= medianPower FFTZoom = FFT[:,idx] i = np.where(FFTZoom < 10.0) FFTZoom[i] = np.nan FreqAxis, TimeAxis = np.meshgrid(freq_mhz[idx], FFTtime) pltImage = axs.pcolormesh(TimeAxis, FreqAxis, FFTZoom, vmin=-10, vmax=10, cmap='rainbow') XMin = np.min(TimeAxis) XMax = np.max(TimeAxis) axs.set_xlim((XMin, XMax)) YMin = np.min(FreqAxis) YMax = np.max(FreqAxis) axs.set_ylim((YMin, YMax)) axs.set_title(fileBase + ' Power L1 +/- 2MHz Outliers') axs.set_xlabel('GPS Time [Sec]') axs.set_ylabel('Frequency [MHz]') # Set the font size of the X and Y axis text axs.tick_params(axis='both', which='major', labelsize=6) # Now add a colorbar to the plot cbar = fig.colorbar(pltImage, ax=axs) # Add a label to the colorbar rotated by 90 degrees cbar.set_label('Power [dB]', rotation=90) # Set the font size of the colorbar text cbar.ax.tick_params(labelsize=6) plt.tight_layout() fig.savefig(fileBase + '-Jam-Outliers.png',dpi=600) plt.close() # Main entry point if __name__ == '__main__': # Parse arguments parser = argparse.ArgumentParser(description='FFT analysis tool') parser.add_argument("filename", help="space delimited FFT file") parser.add_argument('-f','--force', help='Forces generation of a spectrum "median" calibration file',action='store_true') args = parser.parse_args() if(args.force): force = args.force else: force = False filename = args.filename filename = 'Auger5/2024-04-15-16-Auger5-L1.log' #filename = 'Alloy156/2024-04-15-16-Alloy156-L1.log' tokens = filename.split('/') inputFile = tokens[1] station = tokens[0] band = tokens[1].split('-')[-1][:2] # Load the file using pandas - this library implements a fast load fid = open(filename,'rb') csv_args = {'skiprows':0,'header':None,'delim_whitespace':True} #csv_args['on_bad_lines']='skip' FFT = pd.read_csv(fid,**csv_args).apply(pd.to_numeric,errors='coerce').values fid.close() # Unless something changes this should be 2048 # The 3 is due to # Col 0 = Time in milliseconds # Col 1 = Weeknumber # Col 2 = Center frequency FFTLen = len(FFT[0,3:]) # The FFTs are 50MHz binSize = 50.0/FFTLen # Get the center frequency. It should be repeated on every row, so # just get the first entry thisFreq = FFT[0,2] # In MHz freq_mhz = (np.r_[0:2048] - 1024) * binSize + thisFreq if( (force == True) or os.path.exists(station + '-' + band + '.cal') == False): # Get the median of the FFT FFTCal = np.median(FFT[:,3:],axis=0) # Save the median data to a calibration file in a space delimited file np.savetxt(station + '-' + band + '.cal',FFTCal,fmt='%.4f') # Now plot the calibration data fig, axs = plt.subplots(1, 1) pltImage = axs.plot(freq_mhz, FFTCal) axs.set_title(station + ' ' + band + ' Calibration Data') axs.set_xlabel('Frequency [MHz]') axs.set_ylabel('Power [dB]') axs.grid(True) axs.tick_params(axis='both', which='major', labelsize=6) plt.tight_layout() fig.savefig(station+'-'+band+'-cal.png',dpi=600) plt.close() else: # Read the calibration file for this station / band FFTCal = np.loadtxt(station + '-' + band + '.cal') #FreqAxis, TimeAxis = np.meshgrid(freq_mhz, FFT[:, 0]/1000) #thisPlot = 0 for thisPlot in range(2): #if(True): fig, axs = plt.subplots(1, 1) FreqAxis, TimeAxis = np.meshgrid(freq_mhz, FFT[:, 0]/1000) if(thisPlot == 0): titleStr = filename plotName = inputFile + '.png' pltImage = axs.pcolormesh(TimeAxis, FreqAxis, FFT[:,3:], vmin=20, vmax=70, cmap='rainbow') else: titleStr = filename + ' - Calibration' plotName = inputFile + '-cal.png' # Remove the calibration - this is the "expected" noise floor of the FFT. It will # vary across the band due to the shape of the receiver and antenna filters, along # with GNSS transmissions. For example, there is an expected "bump" around the L1 # C/A center frequency (1575.42MHz). By removing the expected noise floor, we can # look for anomalies. The user does need to be careful to ensure a jammer is not # active during the calibration! FFT[:,3:] -= FFTCal pltImage = axs.pcolormesh(TimeAxis, FreqAxis, FFT[:,3:], vmin=-10, vmax=10, cmap='rainbow') analyzeFFT(FFT[:,3:], FFT[:, 0]/1000, freq_mhz, band, inputFile) # Get the time-axis start/stop values XMin = np.min(TimeAxis) XMax = np.max(TimeAxis) axs.set_xlim((XMin, XMax)) axs.set_title(titleStr) axs.set_xlabel('GPS Time [Sec]') axs.set_ylabel('Frequency [MHz]') # Set the font size of the X and Y axis text axs.tick_params(axis='both', which='major', labelsize=6) # Now add a colorbar to the plot cbar = fig.colorbar(pltImage, ax=axs) # Add a label to the colorbar rotated by 90 degrees cbar.set_label('Power [dB]', rotation=90) # Set the font size of the colorbar text cbar.ax.tick_params(labelsize=6) plt.tight_layout() #plt.show() fig.savefig(plotName,dpi=600) plt.close() if(thisPlot == 0): #idx = np.where((freq_mhz >= 1573.42) & (freq_mhz <= 1576.42))[0] idx = np.where((freq_mhz >= 1572.42) & (freq_mhz <= 1577.42))[0] FreqAxis, TimeAxis = np.meshgrid(freq_mhz[idx], FFT[:, 0]/1000) # Now use plotly to generte an interactive 3D mesh plot of the # data and save to an HTML file fig = plotly.graph_objs.Figure(data=[plotly.graph_objs.Surface(z=FFT[:,idx+3], x=TimeAxis, y=FreqAxis, colorscale='Rainbow')]) # Ensure that the x and y axis tick labels are full numbers without commas and not abbreivated or offset # set the z range from 20 to 70 fig.update_layout(scene=dict(xaxis=dict(tickmode='linear', tick0=0, dtick=1000, tickformat='d'), yaxis=dict(tickmode='linear', tick0=0, dtick=1, tickformat='d'), zaxis=dict(range=[20,70]))) # Set the x-axis label fig.update_layout(scene=dict(xaxis_title='Time [GPS Sec]')) # Set the y-axis label fig.update_layout(scene=dict(yaxis_title='Frequency [MHz]')) # Set the z-axis label fig.update_layout(scene=dict(zaxis_title='Power [dB]')) # Set the title fig.update_layout(title=inputFile + ' FFT') print('Save figure to html') # Save without the embedded Javascript, reduces the file size, but requires a network connection plotly.offline.plot(fig, filename=inputFile + '.html', include_plotlyjs='cdn',auto_open=False)