git.gag.com Git - debian/gnuradio/blob - gr-msdd6000/src/python-examples/msdd_spectrum_sense.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, alpha=0.01):
  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         self.fig = figure(1, facecolor="w", figsize=(12,9))
  76         self.sp  = self.fig.add_subplot(1,1,1)
  77         self.pl  = self.sp.plot(range(100), 100*[1,])
  78
  79         params = {'backend': 'ps',
  80                   'xtick.labelsize': self.axis_font_size,
  81                   'ytick.labelsize': self.axis_font_size,
  82                   'text.usetex': False}
  83         rcParams.update(params)
  84
  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("Magnitude (dB)", fontsize=self.label_font_size, fontweight="bold")
  88         self.text_alpha = figtext(0.10, 0.94, ('Moving average alpha: %s' % alpha), weight="heavy", size=self.text_size)
  89
  90         self.cfreqs = list()
  91         self.freqrange = list()
  92         self.data = list() #array('f')
  93
  94         self.alpha = alpha
  95
  96         self.index = 0
  97         self.full = False
  98         self.last_cfreq = 0
  99
 100         self.sample_rate = sample_rate
 101         self.percent = (1.0-percent)/2.0
 102
 103     def parse(self, msg):
 104         self.center_freq = msg.arg1()
 105         self.vlen = int(msg.arg2())
 106         assert(msg.length() == self.vlen * gr.sizeof_float)
 107
 108
 109         if(self.center_freq < self.last_cfreq):
 110             print "Plotting spectrum\n"
 111             self.full = True
 112
 113             self.pl[0].set_data([self.freqrange, self.data])
 114             self.sp.set_ylim([min(self.data), max(self.data)])
 115             self.sp.set_xlim([min(self.freqrange), max(self.freqrange)])
 116             draw()
 117
 118             self.index = 0
 119             del self.freqrange
 120             self.freqrange = list()
 121             #raw_input()
 122
 123         self.last_cfreq = self.center_freq
 124
 125         startind = int(self.percent * self.vlen)
 126         endind = int((1.0 - self.percent) * self.vlen)
 127
 128         fstep = self.sample_rate / self.vlen
 129         f = [self.center_freq - self.sample_rate/2.0 + i*fstep for i in range(startind, endind)]
 130         self.freqrange += f
 131
 132         t = msg.to_string()
 133         d = struct.unpack('%df' % (self.vlen,), t)
 134
 135         if self.full:
 136             for i in range(startind, endind):
 137                 self.data[self.index] = (1.0-self.alpha)*self.data[self.index] + (self.alpha)*d[i]
 138                 self.index += 1
 139         else:
 140             self.data += [di for di in d[startind:endind]]
 141
 142
 143 class my_top_block(gr.top_block):
 144
 145     def __init__(self):
 146         gr.top_block.__init__(self)
 147
 148         usage = "usage: %prog [options] host min_freq max_freq"
 149         parser = OptionParser(option_class=eng_option, usage=usage)
 150         parser.add_option("-g", "--gain", type="eng_float", default=None,
 151                           help="set gain in dB (default is midpoint)")
 152         parser.add_option("", "--tune-delay", type="eng_float", default=5e-5, metavar="SECS",
 153                           help="time to delay (in seconds) after changing frequency [default=%default]")
 154         parser.add_option("", "--dwell-delay", type="eng_float", default=50e-5, metavar="SECS",
 155                           help="time to dwell (in seconds) at a given frequncy [default=%default]")
 156         parser.add_option("-F", "--fft-size", type="int", default=256,
 157                           help="specify number of FFT bins [default=%default]")
 158         parser.add_option("-d", "--decim", type="intx", default=16,
 159                           help="set decimation to DECIM [default=%default]")
 160         parser.add_option("", "--real-time", action="store_true", default=False,
 161                           help="Attempt to enable real-time scheduling")
 162
 163         (options, args) = parser.parse_args()
 164         if len(args) != 3:
 165             parser.print_help()
 166             sys.exit(1)
 167
 168         self.address  = args[0]
 169         self.min_freq = eng_notation.str_to_num(args[1])
 170         self.max_freq = eng_notation.str_to_num(args[2])
 171
 172         self.decim = options.decim
 173         self.gain  = options.gain
 174
 175         if self.min_freq > self.max_freq:
 176             self.min_freq, self.max_freq = self.max_freq, self.min_freq   # swap them
 177
 178         self.fft_size = options.fft_size
 179
 180         if not options.real_time:
 181             realtime = False
 182         else:
 183             # Attempt to enable realtime scheduling
 184             r = gr.enable_realtime_scheduling()
 185             if r == gr.RT_OK:
 186                 realtime = True
 187             else:
 188                 realtime = False
 189                 print "Note: failed to enable realtime scheduling"
 190
 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
 197         self.src = msdd.source_simple(self.address, self.port)
 198         self.src.set_decim_rate(self.decim)
 199
 200         self.set_gain(self.gain)
 201         self.set_freq(self.min_freq)
 202
 203         s2v = gr.stream_to_vector(gr.sizeof_gr_complex, self.fft_size)
 204
 205         mywindow = window.blackmanharris(self.fft_size)
 206         fft = gr.fft_vcc(self.fft_size, True, mywindow, True)
 207         power = 0
 208         for tap in mywindow:
 209             power += tap*tap
 210
 211         norm = gr.multiply_const_cc(1.0/self.fft_size)
 212         c2mag = gr.complex_to_mag_squared(self.fft_size)
 213
 214         # FIXME the log10 primitive is dog slow
 215         log = gr.nlog10_ff(10, self.fft_size,
 216                            -20*math.log10(self.fft_size)-10*math.log10(power/self.fft_size))
 217
 218         # Set the freq_step to % of the actual data throughput.
 219         # This allows us to discard the bins on both ends of the spectrum.
 220         self.percent = 0.4
 221
 222         self.freq_step = self.percent * self.int_rate
 223         self.min_center_freq = self.min_freq + self.freq_step/2
 224         nsteps = math.ceil((self.max_freq - self.min_freq) / self.freq_step)
 225         self.max_center_freq = self.min_center_freq + (nsteps * self.freq_step)
 226
 227         self.next_freq = self.min_center_freq
 228
 229         tune_delay  = max(0, int(round(options.tune_delay * self.int_rate / self.fft_size)))  # in fft_frames
 230         dwell_delay = max(1, int(round(options.dwell_delay * self.int_rate / self.fft_size))) # in fft_frames
 231
 232         self.msgq = gr.msg_queue(16)
 233         self._tune_callback = tune(self)        # hang on to this to keep it from being GC'd
 234         stats = gr.bin_statistics_f(self.fft_size, self.msgq,
 235                                     self._tune_callback, tune_delay, dwell_delay)
 236
 237         # FIXME leave out the log10 until we speed it up
 238         self.connect(self.src, s2v, fft, c2mag, log, stats)
 239
 240
 241     def set_next_freq(self):
 242         target_freq = self.next_freq
 243         self.next_freq = self.next_freq + self.freq_step
 244         if self.next_freq >= self.max_center_freq:
 245             self.next_freq = self.min_center_freq
 246
 247         if not self.set_freq(target_freq):
 248             print "Failed to set frequency to", target_freq
 249
 250         return target_freq
 251
 252
 253     def set_freq(self, target_freq):
 254         """
 255         Set the center frequency we're interested in.
 256
 257         @param target_freq: frequency in Hz
 258         @rypte: bool
 259
 260         """
 261         return self.src.set_rx_freq(0, target_freq)
 262
 263
 264     def set_gain(self, gain):
 265         self.src.set_pga(0, gain)
 266
 267
 268 def main_loop(tb):
 269     msgparser = parse_msg(tb.int_rate, tb.percent)
 270
 271     while 1:
 272
 273         # Get the next message sent from the C++ code (blocking call).
 274         # It contains the center frequency and the mag squared of the fft
 275         msgparser.parse(tb.msgq.delete_head())
 276
 277         # Print center freq so we know that something is happening...
 278         print msgparser.center_freq
 279
 280         # FIXME do something useful with the data...
 281
 282         # m.data are the mag_squared of the fft output (they are in the
 283         # standard order.  I.e., bin 0 == DC.)
 284         # You'll probably want to do the equivalent of "fftshift" on them
 285         # m.raw_data is a string that contains the binary floats.
 286         # You could write this as binary to a file.
 287
 288
 289 if __name__ == '__main__':
 290     tb = my_top_block()
 291     try:
 292         tb.start()              # start executing flow graph in another thread...
 293         main_loop(tb)
 294
 295     except KeyboardInterrupt:
 296         pass