Hi everybody!
I have modified usrp_spectrum_sense.py to plot the results with gnuplot.
There are two files: widespectrum.py and plot.p
I would like everybody to test it and report me the errors and how can I
improve it.
I’ve used USRPv1 + Flex2400.
Thanks in advance!
Here it goes…
WIDESPECTRUM.PY:
#!/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.
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
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,
-20math.log10(self.fft_size)-10math.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)
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=str(mean(m.data)) #
todo= p + " " + media + ‘\n’ #
power.write(todo) #
if m.center_freq == (max_center_freq-freq_step): # If it gets the
final frecuency
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)
# 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
PLOT.P*
set autoscale
unset logscale
unset label
set xtic auto
set ytic auto
set title “Wideband Spectrum Analyzer”
set xlabel “Frecuency”
set ylabel “Power (dB)”
set grid
plot “power.dat” using 1:2 title ‘Mean power’ with linespoints