#!/usr/bin/env python # # Copyright 2005,2007 Free Software Foundation, Inc. # # This file is part of GNU Radio # # GNU Radio 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, or (at your option) # any later version. # # GNU Radio 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 GNU Radio; see the file COPYING. If not, write to # the Free Software Foundation, Inc., 51 Franklin Street, # Boston, MA 02110-1301, USA. # # usrp_spectrum_sense mod from Santix, discuss-gnuradio-l archive # modified: oswald berthold, 20100201-06 # XXX: # gps # speak # log full fft coeffs # log total mean from gnuradio import gr, gru, eng_notation, optfir, window from gnuradio import audio from gnuradio import usrp from gnuradio.eng_option import eng_option from optparse import OptionParser from usrpm import usrp_dbid import sys import math import struct import Gnuplot, Gnuplot.funcutils # Added to view the results import numpy as np import time import gps def dBm2mW(x): return 10**(x/10)/1000 class tune(gr.feval_dd): """ This class allows C++ code to callback into python. """ def __init__(self, tb): gr.feval_dd.__init__(self) self.tb = tb def eval(self, ignore): """ This method is called from gr.bin_statistics_f when it wants to change the center frequency. This method tunes the front end to the new center frequency, and returns the new frequency as its result. """ try: # We use this try block so that if something goes wrong from here # down, at least we'll have a prayer of knowing what went wrong. # Without this, you get a very mysterious: # # terminate called after throwing an instance of 'Swig::DirectorMethodException' # Aborted # # message on stderr. Not exactly helpful ;) new_freq = self.tb.set_next_freq() return new_freq except Exception, e: print "tune: Exception: ", e class parse_msg(object): def __init__(self, msg): self.center_freq = msg.arg1() self.vlen = int(msg.arg2()) assert(msg.length() == self.vlen * gr.sizeof_float) # FIXME consider using Numarray or NumPy vector t = msg.to_string() self.raw_data = t self.data = struct.unpack('%df' % (self.vlen,), t) class my_top_block(gr.top_block): def __init__(self): gr.top_block.__init__(self) usage = "usage: %prog [options] min_freq max_freq" # Example: ./widespectrum.py 2.23G 2.93G # that is the maximun range of the USRP Flex2400 device. parser = OptionParser(option_class=eng_option, usage=usage) parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=(0,0), help="select USRP Rx side A or B (default=A)") parser.add_option("-g", "--gain", type="eng_float", default=None, help="set gain in dB (default is midpoint)") parser.add_option("", "--tune-delay", type="eng_float", default=1e-3, metavar="SECS", help="time to delay (in seconds) after changing frequency [default=%default]") parser.add_option("", "--dwell-delay", type="eng_float", default=10e-3, metavar="SECS", help="time to dwell (in seconds) at a given frequncy [default=%default]") parser.add_option("-F", "--fft-size", type="int", default=256, help="specify number of FFT bins [default=%default]") parser.add_option("-d", "--decim", type="intx", default=64, help="set decimation to DECIM [default=%default]") parser.add_option("", "--real-time", action="store_true", default=False, help="Attempt to enable real-time scheduling") parser.add_option("-B", "--fusb-block-size", type="int", default=0, help="specify fast usb block size [default=%default]") parser.add_option("-N", "--fusb-nblocks", type="int", default=0, help="specify number of fast usb blocks [default=%default]") (options, args) = parser.parse_args() if len(args) != 2: parser.print_help() sys.exit(1) self.min_freq = eng_notation.str_to_num(args[0]) self.max_freq = eng_notation.str_to_num(args[1]) if self.min_freq > self.max_freq: self.min_freq, self.max_freq = self.max_freq, self.min_freq # swap them # # FIXME We set MANUALLY the physical limits of the device. In this case the USRP Flex2400 limits. # if self.min_freq < 2222000000: # print ("The minimum frequency of this device is 2.222GHz") # self.min_freq = 2222000000 # if self.max_freq < 2222000000: # print ("The minimum frequency of this device is 2.222GHz") # self.max_freq = 2222000000 # if self.min_freq > 2937000000: # print ("The maximun frequency of this device is 2.937GHz") # self.min_freq = 2937000000 # if self.max_freq > 2937000000: # print ("The maximun frequency of this device is 2.937GHz") # self.max_freq = 2937000000 # if self.min_freq == self.max_freq: # print ("Do not use this program for a single frecuency analysis please") # exit() self.fft_size = options.fft_size if not options.real_time: realtime = False else: # Attempt to enable realtime scheduling r = gr.enable_realtime_scheduling() if r == gr.RT_OK: realtime = True else: realtime = False print "Note: failed to enable realtime scheduling" # If the user hasn't set the fusb_* parameters on the command line, # pick some values that will reduce latency. if 1: if options.fusb_block_size == 0 and options.fusb_nblocks == 0: if realtime: # be more aggressive options.fusb_block_size = gr.prefs().get_long('fusb', 'rt_block_size', 1024) options.fusb_nblocks = gr.prefs().get_long('fusb', 'rt_nblocks', 16) else: options.fusb_block_size = gr.prefs().get_long('fusb', 'block_size', 4096) options.fusb_nblocks = gr.prefs().get_long('fusb', 'nblocks', 16) #print "fusb_block_size =", options.fusb_block_size #print "fusb_nblocks =", options.fusb_nblocks # build graph self.u = usrp.source_c(fusb_block_size=options.fusb_block_size, fusb_nblocks=options.fusb_nblocks) adc_rate = self.u.adc_rate() # 64 MS/s usrp_decim = options.decim self.u.set_decim_rate(usrp_decim) usrp_rate = adc_rate / usrp_decim self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec)) self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec) print "Using RX d'board %s" % (self.subdev.side_and_name(),) s2v = gr.stream_to_vector(gr.sizeof_gr_complex, self.fft_size) mywindow = window.blackmanharris(self.fft_size) fft = gr.fft_vcc(self.fft_size, True, mywindow) power = 0 for tap in mywindow: power += tap*tap c2mag = gr.complex_to_mag_squared(self.fft_size) # FIXME the log10 primitive is dog slow log = gr.nlog10_ff(10, self.fft_size, -20*math.log10(self.fft_size)-10*math.log10(power/self.fft_size)) # Set the freq_step to 75% of the actual data throughput. # This allows us to discard the bins on both ends of the spectrum. self.freq_step = 0.75 * usrp_rate self.min_center_freq = self.min_freq + self.freq_step/2 nsteps = math.ceil((self.max_freq - self.min_freq) / self.freq_step) self.max_center_freq = self.min_center_freq + (nsteps * self.freq_step) self.next_freq = self.min_center_freq # We define the minimum, maximum and frequency step in a global statement to use them later. global min_center_freq, max_center_freq, freq_step min_center_freq = self.min_center_freq max_center_freq = self.max_center_freq freq_step = self.freq_step tune_delay = max(0, int(round(options.tune_delay * usrp_rate / self.fft_size))) # in fft_frames dwell_delay = max(1, int(round(options.dwell_delay * usrp_rate / self.fft_size))) # in fft_frames self.msgq = gr.msg_queue(16) self._tune_callback = tune(self) # hang on to this to keep it from being GC'd stats = gr.bin_statistics_f(self.fft_size, self.msgq, self._tune_callback, tune_delay, dwell_delay) # FIXME leave out the log10 until we speed it up self.connect(self.u, s2v, fft, c2mag, log, stats) # self.connect(self.u, s2v, fft, c2mag, stats) if options.gain is None: # if no gain was specified, use the mid-point in dB g = self.subdev.gain_range() options.gain = float(g[0]+g[1])/2 self.set_gain(options.gain) print "gain =", options.gain def set_next_freq(self): target_freq = self.next_freq self.next_freq = self.next_freq + self.freq_step if self.next_freq >= self.max_center_freq: self.next_freq = self.min_center_freq if not self.set_freq(target_freq): print "Failed to set frequency to", target_freq return target_freq def set_freq(self, target_freq): """ Set the center frequency we're interested in. @param target_freq: frequency in Hz @rypte: bool Tuning is a two step process. First we ask the front-end to tune as close to the desired frequency as it can. Then we use the result of that operation and our target_frequency to determine the value for the digital down converter. """ return self.u.tune(0, self.subdev, target_freq) def set_gain(self, gain): self.subdev.set_gain(gain) def mean(data): # Returns the arithmetic mean of a numeric list return sum(data) / len(data) def main_loop(tb): # We give basic information about the Spectrum Analysis print "The start frequency is %s Hz" % min_center_freq print "The final frequency is %s Hz" % max_center_freq print "The frequency step is %s Hz" % freq_step #g = Gnuplot.Gnuplot(debug=1) # dont need this ftm ses = gps.gps() # get gps handle fr_count = 0 # full round count count = 0 # total slice count n = 0 # slice index modulo num_steps num_steps = (tb.max_center_freq - tb.min_center_freq)/tb.freq_step # num_steps: how many steps to go over the specified band num_chunks = 100 # dont use that full_spectrum = np.zeros(num_steps * tb.fft_size) # coefficient vector for the complete band print full_spectrum print "freq_step:", tb.freq_step print "num_steps:", num_steps now = time.localtime() timestr = time.strftime("%Y%m%d-%H%M%S", now) fname = "usrp_spectrum_navigation-" + timestr + "-" + str(num_steps) + "-" + str(num_chunks) + ".dat" f = open("data/" + fname, "w") hname = "usrp_spectrum_values-" + timestr + ".dat" h = open("data/" + hname, "w") num_rounds = 20 # navi: number of full rounds to use for mean mean_band = [] # navi: mean for the subband mean_total = [] # navi: mean over total band ses.send('admosy') # send gps data request ses.poll() # fetch gps answer while 1: # Get the next message sent from the C++ code (blocking call). # It contains the center frequency and the mag squared of the fft m = parse_msg(tb.msgq.delete_head()) # Print center freq so we know that something is happening... #print (m.center_freq) # FIXME do something useful with the data... # Mechanism to save in a file (power.dat) 2 columns, one for the frequencies and the other for the mean of the FFT_SIZE points of m.data # if m.center_freq == min_center_freq: # If we get the minimum frequency, it'll reset the power.dat file # power=open("power.dat", "w") # It will overwrite the power.dat file # power=open("power.dat", "a") # Each loop, it adds a dataline (append) p = str(m.center_freq) # with a frequency and the mean of the 256 FFT samples (Power in dB) media = mean(m.data) # mean_band.append(media) todo = p + " " + str(media) #print todo #power.write(todo) # # update gps ses.send('admosy') ses.poll() # If it gets to the final frecuency, do the 2nd order # averaging if m.center_freq == (max_center_freq-freq_step): a = np.array(mean_band) #print a #b = np.array(a) #print a, type(a) b = a.mean() mean_total.append(b) logm = "round #:" + str(fr_count) + \ ", fmin: " + str(tb.min_center_freq) + \ ", fmax: " + str(tb.max_center_freq) + \ ", dBm: " + str(b) + \ ", mW: " + str(dBm2mW(b)) #print logm # logm = str(tb.min_center_freq) + "\t" + \ # str(tb.max_center_freq) + "\t" + \ # str(b) + "\t" + \ # str(dBm2mW(b)) + "\t" + \ # str(ses.fix.latitude) + "\t" + \ # str(ses.fix.longitude) + "\t" + \ # str(ses.fix.time) + "\t" + \ # str(ses.fix.altitude) # h.write(logm + "\n") #c = np.mean(a) mean_band = [] fr_count += 1 # full-round count if fr_count % num_rounds == 0: d = mean(mean_total) print ses.fix.latitude, ses.fix.longitude, str(d) logm = str(tb.min_center_freq) + "\t" + \ str(tb.max_center_freq) + "\t" + \ str(d) + "\t" + \ str(dBm2mW(d)) + "\t" + \ str(ses.fix.latitude) + "\t" + \ str(ses.fix.longitude) + "\t" + \ str(ses.fix.time) + "\t" + \ str(ses.fix.altitude) h.write(logm + "\n") mean_total = [] # gps stuff # print 'latitude ' , ses.fix.latitude # print 'longitude ' , ses.fix.longitude # print 'time utc ' , ses.utc, ses.fix.time # print 'altitude ' , ses.fix.altitude #sys.exit(0) # # p=str(m.center_freq) # It'll write the last frecuency with its Power in the power.dat file # media=str(mean(m.data)) # # todo= p + " " + media + '\n' # # power.write(todo) # # g.load("plot.p") # Load the plot with the data obtained from URSP # power=open("power.dat", "a") # Without this line, the file will start with the last frecuency #g.hardcopy('spectrum.ps', enhanced=1, color=1) # It does a plot copy to the hard disk (I think there's not enough time to do it) # full log n = count % num_steps full_spectrum[n*tb.fft_size:(n+1)*tb.fft_size] = m.data count += 1 # stepped through all bands once if count % num_steps == 0 and count > 5: # do our stuff c2 = 0 for v in full_spectrum: f.write(str((c2 * tb.freq_step/tb.fft_size) + tb.min_center_freq) + "\t" \ + str(v) + "\t" + \ str(ses.fix.latitude) + "\t" + \ str(ses.fix.longitude) + "\t" + \ str(ses.fix.time) + "\t" + \ str(ses.fix.altitude) + "\t" + \ "\n") c2 += 1 f.write("\n\n") # for gnuplot # m.data in 'w' mode: only write, if it exist a file with the same name, it'll be overwrite. # 'a' to append # 'r+' for read and write # m.data are the mag_squared of the fft output (they are in the # standard order. I.e., bin 0 == DC.) # You'll probably want to do the equivalent of "fftshift" on them # m.raw_data is a string that contains the binary floats. # You could write this as binary to a file. if __name__ == '__main__': tb = my_top_block() try: tb.start() # start executing flow graph in another thread... main_loop(tb) except KeyboardInterrupt: pass