git.gag.com Git - debian/gnuradio/blob - gr-msdd6000/src/python-examples/msdd_spectrum_waterfall.py

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2008 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   7 # GNU Radio is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 3, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23 from gnuradio import gr, gru, eng_notation, optfir, window
  24 from gnuradio import msdd
  25 from gnuradio.eng_option import eng_option
  26 from optparse import OptionParser
  27 import sys
  28 import math
  29 import struct
  30 from pylab import *
  31 from numpy import array
  32 import time
  33
  34 matplotlib.interactive(True)
  35 matplotlib.use('TkAgg')
  36
  37 class tune(gr.feval_dd):
  38     """
  39     This class allows C++ code to callback into python.
  40     """
  41     def __init__(self, tb):
  42         gr.feval_dd.__init__(self)
  43         self.tb = tb
  44
  45     def eval(self, ignore):
  46         """
  47         This method is called from gr.bin_statistics_f when it wants to change
  48         the center frequency.  This method tunes the front end to the new center
  49         frequency, and returns the new frequency as its result.
  50         """
  51         try:
  52             # We use this try block so that if something goes wrong from here
  53             # down, at least we'll have a prayer of knowing what went wrong.
  54             # Without this, you get a very mysterious:
  55             #
  56             #   terminate called after throwing an instance of 'Swig::DirectorMethodException'
  57             #   Aborted
  58             #
  59             # message on stderr.  Not exactly helpful ;)
  60
  61             new_freq = self.tb.set_next_freq()
  62             return new_freq
  63
  64         except Exception, e:
  65             print "tune: Exception: ", e
  66
  67
  68 class parse_msg(object):
  69     def __init__(self, sample_rate, percent):
  70         self.axis_font_size = 16
  71         self.label_font_size = 18
  72         self.title_font_size = 20
  73         self.text_size = 22
  74
  75         # Set up figures and subplots
  76         self.fig = figure(1, facecolor="w", figsize=(12,9))
  77         self.sp  = self.fig.add_subplot(1,1,1)
  78         self.pl  = self.sp.matshow(100*[range(100),])
  79
  80         params = {'xtick.labelsize': self.axis_font_size,
  81                   'ytick.labelsize': self.axis_font_size}
  82         rcParams.update(params)
  83
  84         # Throw up some title info
  85         self.sp.set_title(("FFT"), fontsize=self.title_font_size, fontweight="bold")
  86         self.sp.set_xlabel("Frequency (Hz)", fontsize=self.label_font_size, fontweight="bold")
  87         self.sp.set_ylabel("Sample index (should be time)", fontsize=self.label_font_size, fontweight="bold")
  88
  89         self.freqrange = list()
  90         self.data = list()
  91         self.data3 = list()
  92
  93         self.index = 0
  94         self.last_cfreq = 0
  95
  96         # So we know how to splice the data
  97         self.sample_rate = sample_rate
  98         self.percent = (1.0-percent)/2.0
  99
 100     def parse(self, msg):
 101         self.center_freq = msg.arg1()    # read the current center frequency
 102         self.vlen = int(msg.arg2())      # read the length of the data set received
 103
 104         # wait until we wrap around before plotting the entire collected band
 105         if(self.center_freq < self.last_cfreq):
 106             #print "Plotting spectrum\n"
 107
 108             # If we have 100 sets, start dropping the oldest
 109             if(len(self.data3) > 100):
 110                 self.data3.pop(0)
 111             self.data3.append(self.data)
 112
 113             # add the new data to the plot
 114             self.pl.set_data(self.data3)
 115             draw()
 116
 117             # reset lists to collect next round
 118             self.index = 0
 119             del self.freqrange
 120             self.freqrange = list()
 121             del self.data
 122             self.data = list()
 123             #raw_input()
 124
 125         self.last_cfreq = self.center_freq
 126
 127         startind = int(self.percent * self.vlen)
 128         endind = int((1.0 - self.percent) * self.vlen)
 129
 130         fstep = self.sample_rate / self.vlen
 131         f = [self.center_freq - self.sample_rate/2.0 + i*fstep for i in range(startind, endind)]
 132         self.freqrange += f
 133
 134         t = msg.to_string();
 135
 136         d = struct.unpack('%df' % (self.vlen,), t)
 137
 138         self.data += [di for di in d[startind:endind]]
 139
 140
 141 class my_top_block(gr.top_block):
 142     def __init__(self):
 143         gr.top_block.__init__(self)
 144
 145         # Build an options parser to bring in information from the user on usage
 146         usage = "usage: %prog [options] host min_freq max_freq"
 147         parser = OptionParser(option_class=eng_option, usage=usage)
 148         parser.add_option("-g", "--gain", type="eng_float", default=32,
 149                           help="set gain in dB (default is midpoint)")
 150         parser.add_option("", "--tune-delay", type="eng_float", default=5e-5, metavar="SECS",
 151                           help="time to delay (in seconds) after changing frequency [default=%default]")
 152         parser.add_option("", "--dwell-delay", type="eng_float", default=50e-5, metavar="SECS",
 153                           help="time to dwell (in seconds) at a given frequncy [default=%default]")
 154         parser.add_option("-F", "--fft-size", type="int", default=256,
 155                           help="specify number of FFT bins [default=%default]")
 156         parser.add_option("-d", "--decim", type="intx", default=16,
 157                           help="set decimation to DECIM [default=%default]")
 158         parser.add_option("", "--real-time", action="store_true", default=False,
 159                           help="Attempt to enable real-time scheduling")
 160
 161         (options, args) = parser.parse_args()
 162         if len(args) != 3:
 163             parser.print_help()
 164             sys.exit(1)
 165
 166         # get user-provided info on address of MSDD and frequency to sweep
 167         self.address  = args[0]
 168         self.min_freq = eng_notation.str_to_num(args[1])
 169         self.max_freq = eng_notation.str_to_num(args[2])
 170
 171         self.decim = options.decim
 172         self.gain  = options.gain
 173
 174         if self.min_freq > self.max_freq:
 175             self.min_freq, self.max_freq = self.max_freq, self.min_freq   # swap them
 176
 177         self.fft_size = options.fft_size
 178
 179         if not options.real_time:
 180             realtime = False
 181         else:
 182             # Attempt to enable realtime scheduling
 183             r = gr.enable_realtime_scheduling()
 184             if r == gr.RT_OK:
 185                 realtime = True
 186             else:
 187                 realtime = False
 188                 print "Note: failed to enable realtime scheduling"
 189
 190         # Sampling rate is hardcoded and cannot be read off device
 191         adc_rate = 102.4e6
 192         self.int_rate = adc_rate / self.decim
 193         print "Sampling rate: ", self.int_rate
 194
 195         # build graph
 196         self.port = 10001   # required port for UDP packets
 197
 198         #       which board,    op mode,        adx,    port
 199 #        self.src = msdd.source_c(0, 1, self.address, self.port)  # build source object
 200
 201         self.conv = gr.interleaved_short_to_complex();
 202
 203         self.src = msdd.source_simple(self.address,self.port);
 204         self.src.set_decim_rate(self.decim)                      # set decimation rate
 205 #        self.src.set_desired_packet_size(0, 1460)                # set packet size to collect
 206
 207         self.set_gain(self.gain)                                 # set receiver's attenuation
 208         self.set_freq(self.min_freq)                             # set receiver's rx frequency
 209
 210         # restructure into vector format for FFT input
 211         s2v = gr.stream_to_vector(gr.sizeof_gr_complex, self.fft_size)
 212
 213         # set up FFT processing block
 214         mywindow = window.blackmanharris(self.fft_size)
 215         fft = gr.fft_vcc(self.fft_size, True, mywindow, True)
 216         power = 0
 217         for tap in mywindow:
 218             power += tap*tap
 219
 220         # calculate magnitude squared of output of FFT
 221         c2mag = gr.complex_to_mag_squared(self.fft_size)
 222
 223         # FIXME the log10 primitive is dog slow
 224         log = gr.nlog10_ff(10, self.fft_size,
 225                            -20*math.log10(self.fft_size)-10*math.log10(power/self.fft_size))
 226
 227         # Set the freq_step to % of the actual data throughput.
 228         # This allows us to discard the bins on both ends of the spectrum.
 229         self.percent = 0.4
 230
 231         # Calculate the frequency steps to use in the collection over the whole bandwidth
 232         self.freq_step = self.percent * self.int_rate
 233         self.min_center_freq = self.min_freq + self.freq_step/2
 234         nsteps = math.ceil((self.max_freq - self.min_freq) / self.freq_step)
 235         self.max_center_freq = self.min_center_freq + (nsteps * self.freq_step)
 236
 237         self.next_freq = self.min_center_freq
 238
 239         # use these values to set receiver settling time between samples and sampling time
 240         # the default values provided seem to work well with the MSDD over 100 Mbps ethernet
 241         tune_delay  = max(0, int(round(options.tune_delay * self.int_rate / self.fft_size)))  # in fft_frames
 242         dwell_delay = max(1, int(round(options.dwell_delay * self.int_rate / self.fft_size))) # in fft_frames
 243
 244         # set up message callback routine to get data from bin_statistics_f block
 245         self.msgq = gr.msg_queue(16)
 246         self._tune_callback = tune(self)        # hang on to this to keep it from being GC'd
 247
 248         # FIXME this block doesn't like to work with negatives because of the "d_max[i]=0" on line
 249         # 151 of gr_bin_statistics_f.cc file. Set this to -10000 or something to get it to work.
 250         stats = gr.bin_statistics_f(self.fft_size, self.msgq,
 251                                     self._tune_callback, tune_delay, dwell_delay)
 252
 253         # FIXME there's a concern over the speed of the log calculation
 254         # We can probably calculate the log inside the stats block
 255         self.connect(self.src, self.conv, s2v, fft, c2mag, log, stats)
 256
 257
 258     def set_next_freq(self):
 259         ''' Find and set the next frequency of the reciver. After going past the maximum frequency,
 260         the frequency is wrapped around to the start again'''
 261         target_freq = self.next_freq
 262         self.next_freq = self.next_freq + self.freq_step
 263         if self.next_freq >= self.max_center_freq:
 264             self.next_freq = self.min_center_freq
 265
 266         if not self.set_freq(target_freq):
 267             print "Failed to set frequency to", target_freq
 268
 269         return target_freq
 270
 271
 272     def set_freq(self, target_freq):
 273         """
 274         Set the center frequency we're interested in.
 275
 276         @param target_freq: frequency in Hz
 277         @rypte: bool
 278
 279         """
 280         return self.src.set_rx_freq(0, target_freq)
 281
 282
 283     def set_gain(self, gain):
 284         self.src.set_pga(0, gain)
 285
 286
 287 def main_loop(tb):
 288     # Set up parser to get data from stats block and display them.
 289     msgparser = parse_msg(tb.int_rate, tb.percent)
 290
 291     while 1:
 292         # Get the next message sent from the C++ code (blocking call).
 293         # It contains the center frequency and the mag squared of the fft
 294         d = tb.msgq.delete_head();
 295         print d.to_string();
 296         msgparser.parse(d)
 297         #print msgparser.center_freq
 298
 299 if __name__ == '__main__':
 300     tb = my_top_block()
 301     try:
 302         tb.start()              # start executing flow graph in another thread...
 303         main_loop(tb)
 304
 305     except KeyboardInterrupt:
 306         pass