diff --git a/CygnusRFI.py b/CygnusRFI.py index 9074f8c..ac4b42a 100644 --- a/CygnusRFI.py +++ b/CygnusRFI.py @@ -1,99 +1,99 @@ -#!/usr/bin/env python -import os -import sys -import argparse -from time import sleep - -# Define argparse arguments -parser = argparse.ArgumentParser() - -parser.add_argument("-b", "--bandwidth", type=str, required=True) -parser.add_argument("-c", "--channels", type=str, default=1024) -parser.add_argument("-t", "--t_int", type=str, default=0.1) -parser.add_argument("-d", "--duration", type=str, required=True) -parser.add_argument("-f", "--fmin", type=float, required=True) -parser.add_argument("-F", "--fmax", type=float, required=True) - -args = parser.parse_args() - -bandwidth=args.bandwidth -channels=args.channels -t_int=args.t_int -duration=args.duration -fmin=args.fmin -fmax=args.fmax - -#Define GRC observation parameters -bandwidth = str(bandwidth) -channels = str(channels) -t_int = str(t_int) -duration = str(duration) -fmin = float(fmin) -fmax = float(fmax) - -nbins = str(int(float(t_int) * float(bandwidth)/float(channels))) - -#Delete pre-existing observation.dat file -for file_name in os.listdir('.'): - if file_name.endswith('.dat'): - try: - os.remove(file_name) - except OSError: - pass -try: - os.remove('rfi_plot.png') -except OSError: - pass - -print('\n\033[1;33;48m+=================================================================+') -print('\033[1;33;48m| \033[1;36;48mCygnusRFI\033[1;33;48m:\033[1;32;48m An open-source Radio Frequency Interference analyzer\033[1;33;48m |') -print('\033[1;33;48m+=================================================================+') -sleep(0.5) -print('\n\033[1;33;48m\033[4;33;48mRFI Measurement Parameters:\033[0;32;48m') -sleep(0.15) -print(('\033[1;32;48mFrequency range to scan: \033[1;36;48m'+str(float(fmin)/1000000)+'-'+str(float(fmax)/1000000)+' MHz')) -sleep(0.15) -print(('\033[1;32;48mBandwidth per spectrum: \033[1;36;48m'+str(float(bandwidth)/1000000)+' MHz')) -sleep(0.15) -print(('\033[1;32;48mIntegration time per spectrum: \033[1;36;48m'+duration+' sec')) -sleep(0.15) -print(('\033[1;32;48mNumber of channels per spectrum (FFT Size should be a power of 2): \033[1;36;48m'+str(channels))) -sleep(0.15) -print(('\033[1;32;48mIntegration time per FFT sample: \033[1;36;48m'+t_int+' sec')) -sleep(0.5) -print(("\n\033[1;32;48mEstimated completion time: \033[1;36;48m"+str(float(duration)*float(fmax-fmin)/float(bandwidth))+" sec")) - -sleep(0.5) -#proceed = eval(input("\n\033[1;36;48mProceed to measurement? [Y/n]: \033[1;33;48m")) -proceed = 'Y' - -if proceed.lower() != 'n' and proceed.lower() != 'no': - print('\n\033[1;33;48m+=================================================================+') - print('\033[1;33;48m| [+] \033[1;32;48m Starting measurement...\033[1;33;48m |') - print('\033[1;33;48m+=================================================================+\n') - - q=0 - for freq in range(int(fmin), int(fmax), int(float(bandwidth))): - print(("\033[1;33;48m\n---------------------------------------------------------------------------\n \033[1;33;48m[*] \033[1;32;48mCurrently monitoring f_center = "+str(0.000001*freq)+" +/- "+str(float(float(bandwidth)*0.000001)/2)+" MHz (iteration: "+str(q)+")...\n\033[1;33;48m---------------------------------------------------------------------------")) - - #Define observation frequency - f_center = str(freq) - - #Execute top_block.py with parameters - print('\033[0m') - sys.argv = ['top_block.py', '--c-freq='+f_center, '--samp-rate='+bandwidth, '--nchan='+channels, '--nbin='+nbins, '--obs-time='+duration] - exec(compile(open('top_block.py', "rb").read(), 'top_block.py', 'exec')) - os.rename('observation.dat', str(q)+'.dat') - - q = q+1 - print('\n\033[1;33;48m+=================================================================+') - print('\033[1;33;48m| \033[1;32;48mMeasurement finished!\033[1;33;48m |') - print('\033[1;33;48m+=================================================================+\n') - - f_center = str(fmin) - sys.argv = ['rfi_plotter.py', 'freq='+f_center, 'samp_rate='+bandwidth, 'nchan='+channels, 'nbin='+nbins, 'n='+str(q), 'dur='+duration, 'fminimum='+str(fmin), 'fmaximum='+str(fmax)] - exec(compile(open('rfi_plotter.py', "rb").read(), 'rfi_plotter.py', 'exec')) - - print('\033[1;32;48mYour data has been saved as \033[1;36;48mrfi_plot.png\033[1;32;48m.') -else: - print('\n\033[1;31;48m[!] Exiting...\n') +#!/usr/bin/env python +import os +import sys +import argparse +from time import sleep + +# Define argparse arguments +parser = argparse.ArgumentParser() + +parser.add_argument("-b", "--bandwidth", type=str, required=True) +parser.add_argument("-c", "--channels", type=str, default=1024) +parser.add_argument("-t", "--t_int", type=str, default=0.1) +parser.add_argument("-d", "--duration", type=str, required=True) +parser.add_argument("-f", "--fmin", type=float, required=True) +parser.add_argument("-F", "--fmax", type=float, required=True) + +args = parser.parse_args() + +bandwidth=args.bandwidth +channels=args.channels +t_int=args.t_int +duration=args.duration +fmin=args.fmin +fmax=args.fmax + +#Define GRC observation parameters +bandwidth = str(bandwidth) +channels = str(channels) +t_int = str(t_int) +duration = str(duration) +fmin = float(fmin) +fmax = float(fmax) + +nbins = str(int(float(t_int) * float(bandwidth)/float(channels))) + +#Delete pre-existing observation.dat file +for file_name in os.listdir('.'): + if file_name.endswith('.dat'): + try: + os.remove(file_name) + except OSError: + pass +try: + os.remove('rfi_plot.png') +except OSError: + pass + +print('\n\033[1;33;48m+=================================================================+') +print('\033[1;33;48m| \033[1;36;48mCygnusRFI\033[1;33;48m:\033[1;32;48m An open-source Radio Frequency Interference analyzer\033[1;33;48m |') +print('\033[1;33;48m+=================================================================+') +sleep(0.5) +print('\n\033[1;33;48m\033[4;33;48mRFI Measurement Parameters:\033[0;32;48m') +sleep(0.15) +print(('\033[1;32;48mFrequency range to scan: \033[1;36;48m'+str(float(fmin)/1000000)+'-'+str(float(fmax)/1000000)+' MHz')) +sleep(0.15) +print(('\033[1;32;48mBandwidth per spectrum: \033[1;36;48m'+str(float(bandwidth)/1000000)+' MHz')) +sleep(0.15) +print(('\033[1;32;48mIntegration time per spectrum: \033[1;36;48m'+duration+' sec')) +sleep(0.15) +print(('\033[1;32;48mNumber of channels per spectrum (FFT Size should be a power of 2): \033[1;36;48m'+str(channels))) +sleep(0.15) +print(('\033[1;32;48mIntegration time per FFT sample: \033[1;36;48m'+t_int+' sec')) +sleep(0.5) +print(("\n\033[1;32;48mEstimated completion time: \033[1;36;48m"+str(float(duration)*float(fmax-fmin)/float(bandwidth))+" sec")) + +sleep(0.5) +#proceed = eval(input("\n\033[1;36;48mProceed to measurement? [Y/n]: \033[1;33;48m")) +proceed = 'Y' + +if proceed.lower() != 'n' and proceed.lower() != 'no': + print('\n\033[1;33;48m+=================================================================+') + print('\033[1;33;48m| [+] \033[1;32;48m Starting measurement...\033[1;33;48m |') + print('\033[1;33;48m+=================================================================+\n') + + q=0 + for freq in range(int(fmin), int(fmax), int(float(bandwidth))): + print(("\033[1;33;48m\n---------------------------------------------------------------------------\n \033[1;33;48m[*] \033[1;32;48mCurrently monitoring f_center = "+str(0.000001*freq)+" +/- "+str(float(float(bandwidth)*0.000001)/2)+" MHz (iteration: "+str(q)+")...\n\033[1;33;48m---------------------------------------------------------------------------")) + + #Define observation frequency + f_center = str(freq) + + #Execute top_block.py with parameters + print('\033[0m') + sys.argv = ['top_block.py', '--c-freq='+f_center, '--samp-rate='+bandwidth, '--nchan='+channels, '--nbin='+nbins, '--obs-time='+duration] + exec(compile(open('top_block.py', "rb").read(), 'top_block.py', 'exec')) + os.rename('observation.dat', str(q)+'.dat') + + q = q+1 + print('\n\033[1;33;48m+=================================================================+') + print('\033[1;33;48m| \033[1;32;48mMeasurement finished!\033[1;33;48m |') + print('\033[1;33;48m+=================================================================+\n') + + f_center = str(fmin) + sys.argv = ['rfi_plotter.py', 'freq='+f_center, 'samp_rate='+bandwidth, 'nchan='+channels, 'nbin='+nbins, 'n='+str(q), 'dur='+duration, 'fminimum='+str(fmin), 'fmaximum='+str(fmax)] + exec(compile(open('rfi_plotter.py', "rb").read(), 'rfi_plotter.py', 'exec')) + + print('\033[1;32;48mYour data has been saved as \033[1;36;48mrfi_plot.png\033[1;32;48m.') +else: + print('\n\033[1;31;48m[!] Exiting...\n') diff --git a/rfi_plotter.py b/rfi_plotter.py index cb9de01..03fa0d0 100644 --- a/rfi_plotter.py +++ b/rfi_plotter.py @@ -1,91 +1,91 @@ -#!/usr/bin/env python -try: - import matplotlib - matplotlib.use('Agg') - - import numpy as np - import matplotlib.pyplot as plt - plt.rcParams.update({'font.size': 32}) - import argparse - from matplotlib.gridspec import GridSpec - - parser = argparse.ArgumentParser() - parser.add_argument('freq') - parser.add_argument('samp_rate') - parser.add_argument('nchan') - parser.add_argument('nbin') - parser.add_argument('n') - parser.add_argument('dur') - parser.add_argument('fminimum') - parser.add_argument('fmaximum') - args = parser.parse_args() - - def decibel(x): - return 10.0*np.log10(x) - #return x - - if __name__ == "__main__": - #Observation parameters - exec(args.freq) - exec(args.samp_rate) - exec(args.nchan) - exec(args.nbin) - exec(args.n) - exec(args.dur) - exec(args.fminimum) - exec(args.fmaximum) - - frequency=freq - - #Load data - z=1000*(np.fromfile("0.dat", dtype="float32").reshape(-1, nchan)/nbin) - allarrays=[] - - #N = amount of .dat files - for i in range(int(n)): - final = 1000*(np.fromfile(str(i)+".dat", dtype="float32").reshape(-1, nchan)/nbin) - allarrays.append(final) - z_final = np.concatenate(allarrays,axis=1) - - #Number of sub-integrations - nsub = z_final.shape[0] - - #Compute average spectrum - zmean = np.mean(z_final, axis=0) - - #Compute time axis - tint = float(nbin*nchan)/samp_rate - t = tint*np.arange(nsub) - - #Compute frequency axis (convert to MHz) - allfreq=[] - for i in range(int(n)): - freq = np.linspace((frequency+samp_rate*i)-0.5*samp_rate, (frequency+samp_rate*i)+0.5*samp_rate, nchan, endpoint=False)*1e-6 - allfreq.append(freq) - freq = np.concatenate(allfreq) - - #Initialize plot - fig = plt.figure(figsize=(5*n,20.25)) - gs = GridSpec(1,1) - - #Plot average spectrum - ax1 = fig.add_subplot(gs[0,0]) - ax1.plot(freq, decibel(zmean), '#3182bd') - ax1.fill_between(freq, decibel(zmean), color='#deebf7') - ax1.set_xlim(np.min(freq), np.max(freq)) #np.min(freq), np.max(freq) - ax1.set_ylim(np.amin(decibel(zmean))-0.5, np.amin(decibel(zmean))+15) - ax1.ticklabel_format(useOffset=False) - ax1.set_xlabel("Frequency (MHz)", labelpad=25) - ax1.set_ylabel("Relative Power (dB)", labelpad=25) - try: - ax1.set_title("\nAveraged RFI Spectrum", pad=25) - except: - ax1.set_title("\nAveraged RFI Spectrum") - ax1.annotate('Frequency range scanned: '+str(float(fminimum)/1000000)+'-'+str(float(fmaximum)/1000000)+' MHz ($\\Delta\\nu$ = '+str(float((fmaximum)-float(fminimum))/1000000)+' MHz)\nBandwidth per spectrum: '+str(float(samp_rate)/1000000)+' MHz\nIntegration time per spectrum: '+str(dur)+" sec\nNumber of channels per spectrum (FFT size): "+str(nchan), xy=(17, 1179), xycoords='axes points', size=32, ha='left', va='top', color='brown') - ax1.grid() - - plt.tight_layout() - plt.savefig("rfi_plot.png") -except Exception as e: - print(e) - pass +#!/usr/bin/env python +try: + import matplotlib + matplotlib.use('Agg') + + import numpy as np + import matplotlib.pyplot as plt + plt.rcParams.update({'font.size': 32}) + import argparse + from matplotlib.gridspec import GridSpec + + parser = argparse.ArgumentParser() + parser.add_argument('freq') + parser.add_argument('samp_rate') + parser.add_argument('nchan') + parser.add_argument('nbin') + parser.add_argument('n') + parser.add_argument('dur') + parser.add_argument('fminimum') + parser.add_argument('fmaximum') + args = parser.parse_args() + + def decibel(x): + return 10.0*np.log10(x) + #return x + + if __name__ == "__main__": + #Observation parameters + exec(args.freq) + exec(args.samp_rate) + exec(args.nchan) + exec(args.nbin) + exec(args.n) + exec(args.dur) + exec(args.fminimum) + exec(args.fmaximum) + + frequency=freq + + #Load data + z=1000*(np.fromfile("0.dat", dtype="float32").reshape(-1, nchan)/nbin) + allarrays=[] + + #N = amount of .dat files + for i in range(int(n)): + final = 1000*(np.fromfile(str(i)+".dat", dtype="float32").reshape(-1, nchan)/nbin) + allarrays.append(final) + z_final = np.concatenate(allarrays,axis=1) + + #Number of sub-integrations + nsub = z_final.shape[0] + + #Compute average spectrum + zmean = np.mean(z_final, axis=0) + + #Compute time axis + tint = float(nbin*nchan)/samp_rate + t = tint*np.arange(nsub) + + #Compute frequency axis (convert to MHz) + allfreq=[] + for i in range(int(n)): + freq = np.linspace((frequency+samp_rate*i)-0.5*samp_rate, (frequency+samp_rate*i)+0.5*samp_rate, nchan, endpoint=False)*1e-6 + allfreq.append(freq) + freq = np.concatenate(allfreq) + + #Initialize plot + fig = plt.figure(figsize=(5*n,20.25)) + gs = GridSpec(1,1) + + #Plot average spectrum + ax1 = fig.add_subplot(gs[0,0]) + ax1.plot(freq, decibel(zmean), '#3182bd') + ax1.fill_between(freq, decibel(zmean), color='#deebf7') + ax1.set_xlim(np.min(freq), np.max(freq)) #np.min(freq), np.max(freq) + ax1.set_ylim(np.amin(decibel(zmean))-0.5, np.amin(decibel(zmean))+15) + ax1.ticklabel_format(useOffset=False) + ax1.set_xlabel("Frequency (MHz)", labelpad=25) + ax1.set_ylabel("Relative Power (dB)", labelpad=25) + try: + ax1.set_title("\nAveraged RFI Spectrum", pad=25) + except: + ax1.set_title("\nAveraged RFI Spectrum") + ax1.annotate('Frequency range scanned: '+str(float(fminimum)/1000000)+'-'+str(float(fmaximum)/1000000)+' MHz ($\\Delta\\nu$ = '+str(float((fmaximum)-float(fminimum))/1000000)+' MHz)\nBandwidth per spectrum: '+str(float(samp_rate)/1000000)+' MHz\nIntegration time per spectrum: '+str(dur)+" sec\nNumber of channels per spectrum (FFT size): "+str(nchan), xy=(17, 1179), xycoords='axes points', size=32, ha='left', va='top', color='brown') + ax1.grid() + + plt.tight_layout() + plt.savefig("rfi_plot.png") +except Exception as e: + print(e) + pass