git.gag.com Git - debian/gnuradio/blob - gr-radio-astronomy/src/python/usrp_psr_receiver.py

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2004,2005 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 2, 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
  24 #
  25 #
  26 # Pulsar receiver application
  27 #
  28 # Performs both harmonic folding analysis
  29 #  and epoch folding analysis
  30 #
  31 #
  32 from gnuradio import gr, gru, blks, audio
  33 import usrp_dbid
  34 from gnuradio import usrp, optfir
  35 from gnuradio import eng_notation
  36 from gnuradio.eng_option import eng_option
  37 from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, form, slider
  38 from optparse import OptionParser
  39 import wx
  40 import sys
  41 import Numeric
  42 import FFT
  43 import ephem
  44 import time
  45 import os
  46 import math
  47
  48
  49 class app_flow_graph(stdgui.gui_flow_graph):
  50     def __init__(self, frame, panel, vbox, argv):
  51         stdgui.gui_flow_graph.__init__(self)
  52
  53         self.frame = frame
  54         self.panel = panel
  55
  56         parser = OptionParser(option_class=eng_option)
  57         parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=(0, 0),
  58                           help="select USRP Rx side A or B (default=A)")
  59         parser.add_option("-d", "--decim", type="int", default=16,
  60                           help="set fgpa decimation rate to DECIM [default=%default]")
  61         parser.add_option("-f", "--freq", type="eng_float", default=None,
  62                           help="set frequency to FREQ", metavar="FREQ")
  63         parser.add_option("-Q", "--observing", type="eng_float", default=0.0,
  64                           help="set observing frequency to FREQ")
  65         parser.add_option("-a", "--avg", type="eng_float", default=1.0,
  66                 help="set spectral averaging alpha")
  67         parser.add_option("-V", "--favg", type="eng_float", default=2.0,
  68                 help="set folder averaging alpha")
  69         parser.add_option("-g", "--gain", type="eng_float", default=None,
  70                           help="set gain in dB (default is midpoint)")
  71         parser.add_option("-l", "--reflevel", type="eng_float", default=30.0,
  72                           help="Set pulse display reference level")
  73         parser.add_option("-L", "--lowest", type="eng_float", default=1.5,
  74                           help="Lowest valid frequency bin")
  75         parser.add_option("-e", "--longitude", type="eng_float", default=-76.02,                          help="Set Observer Longitude")
  76         parser.add_option("-c", "--latitude", type="eng_float", default=44.85,                          help="Set Observer Latitude")
  77         parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT")
  78
  79         parser.add_option ("-t", "--threshold", type="eng_float", default=2.5, help="pulsar threshold")
  80         parser.add_option("-p", "--lowpass", type="eng_float", default=100, help="Pulse spectra cutoff freq")
  81         parser.add_option("-P", "--prefix", default="./", help="File prefix")
  82         parser.add_option("-u", "--pulsefreq", type="eng_float", default=0.748, help="Observation pulse rate")
  83         parser.add_option("-D", "--dm", type="eng_float", default=1.0e-5, help="Dispersion Measure")
  84         parser.add_option("-O", "--doppler", type="eng_float", default=1.0, help="Doppler ratio")
  85         parser.add_option("-B", "--divbase", type="eng_float", default=20, help="Y/Div menu base")
  86         parser.add_option("-I", "--division", type="eng_float", default=100, help="Y/Div")
  87         parser.add_option("-A", "--audio_source", default="plughw:0,0", help="Audio input device spec")
  88         (options, args) = parser.parse_args()
  89         if len(args) != 0:
  90             parser.print_help()
  91             sys.exit(1)
  92
  93         self.show_debug_info = True
  94
  95         self.reflevel = options.reflevel
  96         self.divbase = options.divbase
  97         self.division = options.division
  98         self.audiodev = options.audio_dev
  99
 100         # Low-pass cutoff for post-detector filter
 101         # Set to 100Hz usually, since lots of pulsars fit in this
 102         #   range
 103         self.lowpass = options.lowpass
 104
 105         # What is lowest valid frequency bin in post-detector FFT?
 106         # There's some pollution very close to DC
 107         self.lowest_freq = options.lowest
 108
 109         # What (dB) threshold to use in determining spectral candidates
 110         self.threshold = options.threshold
 111
 112         # Filename prefix for recording file
 113         self.prefix = options.prefix
 114
 115         # Dispersion Measure (DM)
 116         self.dm = options.dm
 117
 118         # Doppler shift, as a ratio
 119         #  1.0 == no doppler shift
 120         #  1.005 == a little negative shift
 121         #  0.995 == a little positive shift
 122         self.doppler = options.doppler
 123
 124         #
 125         # Input frequency and observing frequency--not necessarily the
 126         #   same thing, if we're looking at the IF of some downconverter
 127         #   that's ahead of the USRP and daughtercard.  This distinction
 128         #   is important in computing the correct de-dispersion filter.
 129         #
 130         self.frequency = options.freq
 131         if options.observing <= 0:
 132             self.observing_freq = options.freq
 133         else:
 134             self.observing_freq = options.observing
 135
 136         # build the graph
 137         self.u = usrp.source_c(decim_rate=options.decim)
 138         self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
 139
 140         #
 141         # Recording file, in case we ever need to record baseband data
 142         #
 143         self.recording = gr.file_sink(gr.sizeof_char, "/dev/null")
 144         self.recording_state = False
 145
 146         self.pulse_recording = gr.file_sink(gr.sizeof_short, "/dev/null")
 147         self.pulse_recording_state = False
 148
 149         #
 150         # We come up with recording turned off, but the user may
 151         #  request recording later on
 152         self.recording.close()
 153         self.pulse_recording.close()
 154
 155         #
 156         # Need these two for converting 12-bit baseband signals to 8-bit
 157         #
 158         self.tofloat = gr.complex_to_float()
 159         self.tochar = gr.float_to_char()
 160
 161         # Need this for recording pulses (post-detector)
 162         self.toshort = gr.float_to_short()
 163
 164
 165         #
 166         # The spectral measurer sets this when it has a valid
 167         #   average spectral peak-to-peak distance
 168         # We can then use this to program the parameters for the epoch folder
 169         #
 170         # We set a sentimental value here
 171         self.pulse_freq = options.pulsefreq
 172
 173         # Folder runs at this raw sample rate
 174         self.folder_input_rate = 20000
 175
 176         # Each pulse in the epoch folder is sampled at 128 times the nominal
 177         #  pulse rate
 178         self.folding = 128
 179
 180
 181         #
 182         # Try to find candidate parameters for rational resampler
 183         #
 184         save_i = 0
 185         candidates = []
 186         for i in range(20,300):
 187             input_rate = self.folder_input_rate
 188             output_rate = int(self.pulse_freq * i)
 189             interp = gru.lcm(input_rate, output_rate) / input_rate
 190             decim = gru.lcm(input_rate, output_rate) / output_rate
 191             if (interp < 500 and decim < 250000):
 192                  candidates.append(i)
 193
 194         # We didn't find anything, bail!
 195         if (len(candidates) < 1):
 196             print "Couldn't converge on resampler parameters"
 197             sys.exit(1)
 198
 199         #
 200         # Now try to find candidate with the least sampling error
 201         #
 202         mindiff = 999.999
 203         for i in candidates:
 204             diff = self.pulse_freq * i
 205             diff = diff - int(diff)
 206             if (diff < mindiff):
 207                 mindiff = diff
 208                 save_i = i
 209
 210         # Recompute rates
 211         input_rate = self.folder_input_rate
 212         output_rate = int(self.pulse_freq * save_i)
 213
 214         # Compute new interp and decim, based on best candidate
 215         interp = gru.lcm(input_rate, output_rate) / input_rate
 216         decim = gru.lcm(input_rate, output_rate) / output_rate
 217
 218         # Save optimized folding parameters, used later
 219         self.folding = save_i
 220         self.interp = int(interp)
 221         self.decim = int(decim)
 222
 223         # So that we can view 4 pulses in the pulse viewer window
 224         FOLD_MULT=1
 225
 226         # determine the daughterboard subdevice we're using
 227         self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
 228         self.cardtype = self.u.daughterboard_id(0)
 229
 230         # Compute raw input rate
 231         input_rate = self.u.adc_freq() / self.u.decim_rate()
 232
 233         # BW==input_rate for complex data
 234         self.bw = input_rate
 235
 236         #
 237         # Set baseband filter bandwidth if DBS_RX:
 238         #
 239         if self.cardtype == usrp_dbid.DBS_RX:
 240             lbw = input_rate / 2
 241             if lbw < 1.0e6:
 242                 lbw = 1.0e6
 243             self.subdev.set_bw(lbw)
 244
 245         #
 246         # We use this as a crude volume control for the audio output
 247         #
 248         self.volume = gr.multiply_const_ff(10**(-1))
 249
 250
 251         #
 252         # Create location data for ephem package
 253         #
 254         self.locality = ephem.Observer()
 255         self.locality.long = str(options.longitude)
 256         self.locality.lat = str(options.latitude)
 257
 258         #
 259         # What is the post-detector LPF cutoff for the FFT?
 260         #
 261         PULSAR_MAX_FREQ=int(options.lowpass)
 262
 263         # First low-pass filters down to input_rate/FIRST_FACTOR
 264         #   and decimates appropriately
 265         FIRST_FACTOR=int(input_rate/(self.folder_input_rate/2))
 266         first_filter = gr.firdes.low_pass (1.0,
 267                                           input_rate,
 268                                           input_rate/FIRST_FACTOR,
 269                                           input_rate/(FIRST_FACTOR*20),
 270                                           gr.firdes.WIN_HAMMING)
 271
 272         # Second filter runs at the output rate of the first filter,
 273         #  And low-pass filters down to PULSAR_MAX_FREQ*10
 274         #
 275         second_input_rate =  int(input_rate/(FIRST_FACTOR/2))
 276         second_filter = gr.firdes.band_pass(1.0, second_input_rate,
 277                                           0.10,
 278                                           PULSAR_MAX_FREQ*10,
 279                                           PULSAR_MAX_FREQ*1.5,
 280                                           gr.firdes.WIN_HAMMING)
 281
 282         # Third filter runs at PULSAR_MAX_FREQ*20
 283         #   and filters down to PULSAR_MAX_FREQ
 284         #
 285         third_input_rate = PULSAR_MAX_FREQ*20
 286         third_filter = gr.firdes_band_pass(1.0, third_input_rate,
 287                                            0.10, PULSAR_MAX_FREQ,
 288                                            PULSAR_MAX_FREQ/10.0,
 289                                            gr.firdes.WIN_HAMMING)
 290
 291
 292         #
 293         # Create the appropriate FFT scope
 294         #
 295         self.scope = ra_fftsink.ra_fft_sink_f (self, panel,
 296            fft_size=int(options.fft_size), sample_rate=PULSAR_MAX_FREQ*2,
 297            title="Post-detector spectrum",
 298            cfunc=self.pulsarfunc, xydfunc=self.xydfunc, fft_rate=200)
 299
 300         #
 301         # Tell scope we're looking from DC to PULSAR_MAX_FREQ
 302         #
 303         self.scope.set_baseband_freq (0.0)
 304
 305
 306         #
 307         # Setup stripchart for showing pulse profiles
 308         #
 309         hz = "%5.3fHz " % self.pulse_freq
 310         per = "(%5.3f sec)" % (1.0/self.pulse_freq)
 311         sr = "%d sps" % (int(self.pulse_freq*self.folding))
 312         self.chart = ra_stripchartsink.stripchart_sink_f (self, panel,
 313                sample_rate=1,
 314                stripsize=self.folding*FOLD_MULT, parallel=True, title="Pulse Profiles: "+hz+per,
 315                xlabel="Seconds @ "+sr, ylabel="Level", autoscale=True,
 316                divbase=self.divbase, scaling=1.0/(self.folding*self.pulse_freq))
 317         self.chart.set_ref_level(self.reflevel)
 318         self.chart.set_y_per_div(self.division)
 319
 320         # De-dispersion filter setup
 321         #
 322         # Do this here, just before creating the filter
 323         #  that will use the taps.
 324         #
 325         ntaps = self.compute_disp_ntaps(self.dm,self.bw,self.observing_freq)
 326
 327         # Taps for the de-dispersion filter
 328         self.disp_taps = Numeric.zeros(ntaps,Numeric.Complex64)
 329
 330         # Compute the de-dispersion filter now
 331         self.compute_dispfilter(self.dm,self.doppler,
 332             self.bw,self.observing_freq)
 333
 334         #
 335         # Call constructors for receive chains
 336         #
 337
 338         #
 339         # Now create the FFT filter using the computed taps
 340         self.dispfilt = gr.fft_filter_ccc(1, self.disp_taps)
 341
 342         #
 343         # Audio sink
 344         #
 345         self.audio = audio.sink(second_input_rate, self.audiodev)
 346
 347         #
 348         # The three post-detector filters
 349         # Done this way to allow an audio path (up to 10Khz)
 350         # ...and also because going from xMhz down to ~100Hz
 351         # In a single filter doesn't seem to work.
 352         #
 353         self.first = gr.fir_filter_fff (FIRST_FACTOR/2, first_filter)
 354
 355         p = second_input_rate / (PULSAR_MAX_FREQ*20)
 356         self.second = gr.fir_filter_fff (int(p), second_filter)
 357         self.third = gr.fir_filter_fff (10, third_filter)
 358
 359         # Split complex USRP stream into a pair of floats
 360         self.splitter = gr.complex_to_float (1);
 361
 362         # I squarer (detector)
 363         self.multI = gr.multiply_ff();
 364
 365         # Q squarer (detector)
 366         self.multQ = gr.multiply_ff();
 367
 368         # Adding squared I and Q to produce instantaneous signal power
 369         self.adder = gr.add_ff();
 370
 371         self.enable_comb_filter = False
 372         # Epoch folder comb filter
 373         if self.enable_comb_filter == True:
 374             bogtaps = Numeric.zeros(512, Numeric.Float64)
 375             self.folder_comb = gr.fft_filter_ccc(1,bogtaps)
 376
 377         # Rational resampler
 378         self.folder_rr = blks.rational_resampler_fff(self, self.interp, self.decim)
 379
 380         # Epoch folder bandpass
 381         bogtaps = Numeric.zeros(1, Numeric.Float64)
 382         self.folder_bandpass = gr.fir_filter_fff (1, bogtaps)
 383
 384         # Epoch folder F2C/C2F
 385         self.folder_f2c = gr.float_to_complex()
 386         self.folder_c2f = gr.complex_to_float()
 387
 388         # Epoch folder S2P
 389         self.folder_s2p = gr.serial_to_parallel (gr.sizeof_float,
 390              self.folding*FOLD_MULT)
 391
 392         # Epoch folder IIR Filter (produces average pulse profiles)
 393         self.folder_iir = gr.single_pole_iir_filter_ff(1.0/options.favg,
 394              self.folding*FOLD_MULT)
 395
 396         #
 397         # Set all the epoch-folder goop up
 398         #
 399         self.set_folding_params()
 400
 401         #
 402         # Start connecting configured modules in the receive chain
 403         #
 404
 405         # Connect raw USRP to de-dispersion filter, complex->float splitter
 406         self.connect(self.u, self.dispfilt, self.splitter)
 407
 408         # Connect splitter outputs to multipliers
 409         # First do I^2
 410         self.connect((self.splitter, 0), (self.multI,0))
 411         self.connect((self.splitter, 0), (self.multI,1))
 412
 413         # Then do Q^2
 414         self.connect((self.splitter, 1), (self.multQ,0))
 415         self.connect((self.splitter, 1), (self.multQ,1))
 416
 417         # Then sum the squares
 418         self.connect(self.multI, (self.adder,0))
 419         self.connect(self.multQ, (self.adder,1))
 420
 421         # Connect detector/adder output to FIR LPF
 422         #  in two stages, followed by the FFT scope
 423         self.connect(self.adder, self.first,
 424             self.second, self.third, self.scope)
 425
 426         # Connect audio output
 427         self.connect(self.first, self.volume)
 428         self.connect(self.volume, (self.audio, 0))
 429         self.connect(self.volume, (self.audio, 1))
 430
 431         # Connect epoch folder
 432         if self.enable_comb_filter == True:
 433             self.connect (self.first, self.folder_bandpass, self.folder_rr,
 434                 self.folder_f2c,
 435                 self.folder_comb, self.folder_c2f,
 436                 self.folder_s2p, self.folder_iir,
 437                 self.chart)
 438
 439         else:
 440             self.connect (self.first, self.folder_bandpass, self.folder_rr,
 441                 self.folder_s2p, self.folder_iir, self.chart)
 442
 443         # Connect baseband recording file (initially /dev/null)
 444         self.connect(self.u, self.tofloat, self.tochar, self.recording)
 445
 446         # Connect pulse recording file (initially /dev/null)
 447         self.connect(self.first, self.toshort, self.pulse_recording)
 448
 449         #
 450         # Build the GUI elements
 451         #
 452         self._build_gui(vbox)
 453
 454         # Make GUI agree with command-line
 455         self.myform['average'].set_value(int(options.avg))
 456         self.myform['foldavg'].set_value(int(options.favg))
 457
 458
 459         # Make spectral averager agree with command line
 460         if options.avg != 1.0:
 461             self.scope.set_avg_alpha(float(1.0/options.avg))
 462             self.scope.set_average(True)
 463
 464
 465         # set initial values
 466
 467         if options.gain is None:
 468             # if no gain was specified, use the mid-point in dB
 469             g = self.subdev.gain_range()
 470             options.gain = float(g[0]+g[1])/2
 471
 472         if options.freq is None:
 473             # if no freq was specified, use the mid-point
 474             r = self.subdev.freq_range()
 475             options.freq = float(r[0]+r[1])/2
 476
 477         self.set_gain(options.gain)
 478         self.set_volume(-10.0)
 479
 480         if not(self.set_freq(options.freq)):
 481             self._set_status_msg("Failed to set initial frequency")
 482
 483         self.myform['decim'].set_value(self.u.decim_rate())
 484         self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
 485         self.myform['dbname'].set_value(self.subdev.name())
 486         self.myform['DM'].set_value(self.dm)
 487         self.myform['Doppler'].set_value(self.doppler)
 488
 489         #
 490         # Start the timer that shows current LMST on the GUI
 491         #
 492         self.lmst_timer.Start(1000)
 493
 494
 495     def _set_status_msg(self, msg):
 496         self.frame.GetStatusBar().SetStatusText(msg, 0)
 497
 498     def _build_gui(self, vbox):
 499
 500         def _form_set_freq(kv):
 501             return self.set_freq(kv['freq'])
 502
 503         def _form_set_dm(kv):
 504             return self.set_dm(kv['DM'])
 505
 506         def _form_set_doppler(kv):
 507             return self.set_doppler(kv['Doppler'])
 508
 509         # Position the FFT or Waterfall
 510         vbox.Add(self.scope.win, 5, wx.EXPAND)
 511         vbox.Add(self.chart.win, 5, wx.EXPAND)
 512
 513         # add control area at the bottom
 514         self.myform = myform = form.form()
 515         hbox = wx.BoxSizer(wx.HORIZONTAL)
 516         hbox.Add((7,0), 0, wx.EXPAND)
 517         vbox1 = wx.BoxSizer(wx.VERTICAL)
 518         myform['freq'] = form.float_field(
 519             parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
 520             callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
 521
 522         vbox1.Add((3,0), 0, 0)
 523
 524         # To show current Local Mean Sidereal Time
 525         myform['lmst_high'] = form.static_text_field(
 526             parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
 527         vbox1.Add((3,0), 0, 0)
 528
 529         # To show current spectral cursor data
 530         myform['spec_data'] = form.static_text_field(
 531             parent=self.panel, sizer=vbox1, label="Pulse Freq", weight=1)
 532         vbox1.Add((3,0), 0, 0)
 533
 534         # To show best pulses found in FFT output
 535         myform['best_pulse'] = form.static_text_field(
 536             parent=self.panel, sizer=vbox1, label="Best freq", weight=1)
 537         vbox1.Add((3,0), 0, 0)
 538
 539         vboxBogus = wx.BoxSizer(wx.VERTICAL)
 540         vboxBogus.Add ((2,0), 0, wx.EXPAND)
 541         vbox2 = wx.BoxSizer(wx.VERTICAL)
 542         g = self.subdev.gain_range()
 543         myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
 544                                            weight=1,
 545                                            min=int(g[0]), max=int(g[1]),
 546                                            callback=self.set_gain)
 547
 548         vbox2.Add((6,0), 0, 0)
 549         myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
 550                     label="Spectral Averaging", weight=1, min=1, max=200, callback=self.set_averaging)
 551         vbox2.Add((6,0), 0, 0)
 552         myform['foldavg'] = form.slider_field(parent=self.panel, sizer=vbox2,
 553                     label="Folder Averaging", weight=1, min=1, max=20, callback=self.set_folder_averaging)
 554         vbox2.Add((6,0), 0, 0)
 555         myform['volume'] = form.quantized_slider_field(parent=self.panel, sizer=vbox2,
 556                     label="Audio Volume", weight=1, range=(-20, 0, 0.5), callback=self.set_volume)
 557         vbox2.Add((6,0), 0, 0)
 558         myform['DM'] = form.float_field(
 559             parent=self.panel, sizer=vbox2, label="DM", weight=1,
 560             callback=myform.check_input_and_call(_form_set_dm))
 561         vbox2.Add((6,0), 0, 0)
 562         myform['Doppler'] = form.float_field(
 563             parent=self.panel, sizer=vbox2, label="Doppler", weight=1,
 564             callback=myform.check_input_and_call(_form_set_doppler))
 565         vbox2.Add((6,0), 0, 0)
 566
 567
 568         # Baseband recording control
 569         buttonbox = wx.BoxSizer(wx.HORIZONTAL)
 570         self.record_control = form.button_with_callback(self.panel,
 571               label="Recording baseband: Off                           ",
 572               callback=self.toggle_recording)
 573         self.record_pulse_control = form.button_with_callback(self.panel,
 574               label="Recording pulses: Off                              ",
 575               callback=self.toggle_pulse_recording)
 576
 577         buttonbox.Add(self.record_control, 0, wx.CENTER)
 578         buttonbox.Add(self.record_pulse_control, 0, wx.CENTER)
 579         vbox.Add(buttonbox, 0, wx.CENTER)
 580         hbox.Add(vbox1, 0, 0)
 581         hbox.Add(vboxBogus, 0, 0)
 582         hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
 583         vbox.Add(hbox, 0, wx.EXPAND)
 584
 585         self._build_subpanel(vbox)
 586
 587         self.lmst_timer = wx.PyTimer(self.lmst_timeout)
 588         self.lmst_timeout()
 589
 590
 591     def _build_subpanel(self, vbox_arg):
 592         # build a secondary information panel (sometimes hidden)
 593
 594         # FIXME figure out how to have this be a subpanel that is always
 595         # created, but has its visibility controlled by foo.Show(True/False)
 596
 597         if not(self.show_debug_info):
 598             return
 599
 600         panel = self.panel
 601         vbox = vbox_arg
 602         myform = self.myform
 603
 604         #panel = wx.Panel(self.panel, -1)
 605         #vbox = wx.BoxSizer(wx.VERTICAL)
 606
 607         hbox = wx.BoxSizer(wx.HORIZONTAL)
 608         hbox.Add((5,0), 0)
 609         myform['decim'] = form.static_float_field(
 610             parent=panel, sizer=hbox, label="Decim")
 611
 612         hbox.Add((5,0), 1)
 613         myform['fs@usb'] = form.static_float_field(
 614             parent=panel, sizer=hbox, label="Fs@USB")
 615
 616         hbox.Add((5,0), 1)
 617         myform['dbname'] = form.static_text_field(
 618             parent=panel, sizer=hbox)
 619
 620         hbox.Add((5,0), 1)
 621         myform['baseband'] = form.static_float_field(
 622             parent=panel, sizer=hbox, label="Analog BB")
 623
 624         hbox.Add((5,0), 1)
 625         myform['ddc'] = form.static_float_field(
 626             parent=panel, sizer=hbox, label="DDC")
 627
 628         hbox.Add((5,0), 0)
 629         vbox.Add(hbox, 0, wx.EXPAND)
 630
 631
 632
 633     def set_freq(self, target_freq):
 634         """
 635         Set the center frequency we're interested in.
 636
 637         @param target_freq: frequency in Hz
 638         @rypte: bool
 639
 640         Tuning is a two step process.  First we ask the front-end to
 641         tune as close to the desired frequency as it can.  Then we use
 642         the result of that operation and our target_frequency to
 643         determine the value for the digital down converter.
 644         """
 645         r = usrp.tune(self.u, 0, self.subdev, target_freq)
 646
 647         if r:
 648             self.myform['freq'].set_value(target_freq)     # update displayed value
 649             self.myform['baseband'].set_value(r.baseband_freq)
 650             self.myform['ddc'].set_value(r.dxc_freq)
 651             # Adjust self.frequency, and self.observing_freq
 652             # We pick up the difference between the current self.frequency
 653             #   and the just-programmed one, and use this to adjust
 654             #   self.observing_freq.  We have to do it this way to
 655             #   make the dedispersion filtering work out properly.
 656             delta = target_freq - self.frequency
 657             self.frequency = target_freq
 658             self.observing_freq += delta
 659
 660             # Now that we're adjusted, compute a new dispfilter, and
 661             #   set the taps for the FFT filter.
 662             ntaps = self.compute_disp_ntaps(self.dm, self.bw, self.observing_freq)
 663             self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
 664             self.compute_dispfilter(self.dm,self.doppler,self.bw,
 665                 self.observing_freq)
 666             self.dispfilt.set_taps(self.disp_taps)
 667
 668             return True
 669
 670         return False
 671
 672     # Callback for gain-setting slider
 673     def set_gain(self, gain):
 674         self.myform['gain'].set_value(gain)     # update displayed value
 675         self.subdev.set_gain(gain)
 676
 677
 678     def set_volume(self, vol):
 679         self.myform['volume'].set_value(vol)
 680         self.volume.set_k((10**(vol/10))/8192)
 681
 682     # Callback for spectral-averaging slider
 683     def set_averaging(self, avval):
 684         self.myform['average'].set_value(avval)
 685         self.scope.set_avg_alpha(1.0/(avval))
 686         self.scope.set_average(True)
 687
 688     def set_folder_averaging(self, avval):
 689         self.myform['foldavg'].set_value(avval)
 690         self.folder_iir.set_taps(1.0/avval)
 691
 692     # Timer callback to update LMST display
 693     def lmst_timeout(self):
 694          self.locality.date = ephem.now()
 695          sidtime = self.locality.sidereal_time()
 696          self.myform['lmst_high'].set_value(str(ephem.hours(sidtime)))
 697
 698     #
 699     # Turn recording on/off
 700     # Called-back by "Recording" button
 701     #
 702     def toggle_recording(self):
 703         # Pick up current LMST
 704         self.locality.date = ephem.now()
 705         sidtime = self.locality.sidereal_time()
 706
 707         # Pick up localtime, for generating filenames
 708         foo = time.localtime()
 709
 710         # Generate filenames for both data and header file
 711         filename = "%04d%02d%02d%02d%02d.pdat" % (foo.tm_year, foo.tm_mon,
 712            foo.tm_mday, foo.tm_hour, foo.tm_min)
 713         hdrfilename = "%04d%02d%02d%02d%02d.phdr" % (foo.tm_year, foo.tm_mon,
 714            foo.tm_mday, foo.tm_hour, foo.tm_min)
 715
 716         # Current recording?  Flip state
 717         if (self.recording_state == True):
 718           self.recording_state = False
 719           self.record_control.SetLabel("Recording baseband: Off                           ")
 720           self.recording.close()
 721         # Not recording?
 722         else:
 723           self.recording_state = True
 724           self.record_control.SetLabel("Recording baseband to: "+filename)
 725
 726           # Cause gr_file_sink object to accept new filename
 727           #   note use of self.prefix--filename prefix from
 728           #   command line (defaults to ./)
 729           #
 730           self.recording.open (self.prefix+filename)
 731
 732           #
 733           # We open the header file as a regular file, write header data,
 734           #   then close
 735           hdrf = open(self.prefix+hdrfilename, "w")
 736           hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
 737           hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
 738           hdrf.write("DM: "+str(self.dm)+"\n")
 739           hdrf.write("doppler: "+str(self.doppler)+"\n")
 740
 741           hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 742           hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 743           hdrf.write("sample type: complex_char\n")
 744           hdrf.write("sample size: "+str(gr.sizeof_char*2)+"\n")
 745           hdrf.close()
 746     #
 747     # Turn recording on/off
 748     # Called-back by "Recording" button
 749     #
 750     def toggle_pulse_recording(self):
 751         # Pick up current LMST
 752         self.locality.date = ephem.now()
 753         sidtime = self.locality.sidereal_time()
 754
 755         # Pick up localtime, for generating filenames
 756         foo = time.localtime()
 757
 758         # Generate filenames for both data and header file
 759         filename = "%04d%02d%02d%02d%02d.padat" % (foo.tm_year, foo.tm_mon,
 760            foo.tm_mday, foo.tm_hour, foo.tm_min)
 761         hdrfilename = "%04d%02d%02d%02d%02d.pahdr" % (foo.tm_year, foo.tm_mon,
 762            foo.tm_mday, foo.tm_hour, foo.tm_min)
 763
 764         # Current recording?  Flip state
 765         if (self.pulse_recording_state == True):
 766           self.pulse_recording_state = False
 767           self.record_pulse_control.SetLabel("Recording pulses: Off                           ")
 768           self.pulse_recording.close()
 769         # Not recording?
 770         else:
 771           self.pulse_recording_state = True
 772           self.record_pulse_control.SetLabel("Recording pulses to: "+filename)
 773
 774           # Cause gr_file_sink object to accept new filename
 775           #   note use of self.prefix--filename prefix from
 776           #   command line (defaults to ./)
 777           #
 778           self.pulse_recording.open (self.prefix+filename)
 779
 780           #
 781           # We open the header file as a regular file, write header data,
 782           #   then close
 783           hdrf = open(self.prefix+hdrfilename, "w")
 784           hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
 785           hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
 786           hdrf.write("DM: "+str(self.dm)+"\n")
 787           hdrf.write("doppler: "+str(self.doppler)+"\n")
 788           hdrf.write("pulse rate: "+str(self.pulse_freq)+"\n")
 789           hdrf.write("pulse sps: "+str(self.pulse_freq*self.folding)+"\n")
 790           hdrf.write("file sps: "+str(self.folder_input_rate)+"\n")
 791
 792           hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 793           hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 794           hdrf.write("sample type: short\n")
 795           hdrf.write("sample size: 1\n")
 796           hdrf.close()
 797
 798     # We get called at startup, and whenever the GUI "Set Folding params"
 799     #   button is pressed
 800     #
 801     def set_folding_params(self):
 802         if (self.pulse_freq <= 0):
 803             return
 804
 805         # Compute required sample rate
 806         self.sample_rate = int(self.pulse_freq*self.folding)
 807
 808         # And the implied decimation rate
 809         required_decimation = int(self.folder_input_rate / self.sample_rate)
 810
 811         # We also compute a new FFT comb filter, based on the expected
 812         #  spectral profile of our pulse parameters
 813         #
 814         # FFT-based comb filter
 815         #
 816         N_COMB_TAPS=int(self.sample_rate*4)
 817         if N_COMB_TAPS > 2000:
 818             N_COMB_TAPS = 2000
 819         self.folder_comb_taps = Numeric.zeros(N_COMB_TAPS,Numeric.Complex64)
 820         fincr = (self.sample_rate)/float(N_COMB_TAPS)
 821         for i in range(0,len(self.folder_comb_taps)):
 822             self.folder_comb_taps[i] = complex(0.0, 0.0)
 823
 824         freq = 0.0
 825         harmonics = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
 826         for i in range(0,len(self.folder_comb_taps)/2):
 827             for j in range(0,len(harmonics)):
 828                  if abs(freq - harmonics[j]*self.pulse_freq) <= fincr:
 829                      self.folder_comb_taps[i] = complex(4.0, 0.0)
 830                      if harmonics[j] == 1.0:
 831                          self.folder_comb_taps[i] = complex(8.0, 0.0)
 832             freq += fincr
 833
 834         if self.enable_comb_filter == True:
 835             # Set the just-computed FFT comb filter taps
 836             self.folder_comb.set_taps(self.folder_comb_taps)
 837
 838         # And compute a new decimated bandpass filter, to go in front
 839         #   of the comb.  Primary function is to decimate and filter down
 840         #   to an exact-ish multiple of the target pulse rate
 841         #
 842         self.folding_taps = gr.firdes_band_pass (1.0, self.folder_input_rate,
 843             0.10, self.sample_rate/2, 10,
 844             gr.firdes.WIN_HAMMING)
 845
 846         # Set the computed taps for the bandpass/decimate filter
 847         self.folder_bandpass.set_taps (self.folding_taps)
 848     #
 849     # Record a spectral "hit" of a possible pulsar spectral profile
 850     #
 851     def record_hit(self,hits, hcavg, hcmax):
 852         # Pick up current LMST
 853         self.locality.date = ephem.now()
 854         sidtime = self.locality.sidereal_time()
 855
 856         # Pick up localtime, for generating filenames
 857         foo = time.localtime()
 858
 859         # Generate filenames for both data and header file
 860         hitfilename = "%04d%02d%02d%02d.phit" % (foo.tm_year, foo.tm_mon,
 861            foo.tm_mday, foo.tm_hour)
 862
 863         hitf = open(self.prefix+hitfilename, "a")
 864         hitf.write("receiver center frequency: "+str(self.frequency)+"\n")
 865         hitf.write("observing frequency: "+str(self.observing_freq)+"\n")
 866         hitf.write("DM: "+str(self.dm)+"\n")
 867         hitf.write("doppler: "+str(self.doppler)+"\n")
 868
 869         hitf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 870         hitf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 871         hitf.write("spectral peaks: "+str(hits)+"\n")
 872         hitf.write("HCM: "+str(hcavg)+" "+str(hcmax)+"\n")
 873         hitf.close()
 874
 875     # This is a callback used by ra_fftsink.py (passed on creation of
 876     #   ra_fftsink)
 877     # Whenever the user moves the cursor within the FFT display, this
 878     #   shows the coordinate data
 879     #
 880     def xydfunc(self,xyv):
 881         s = "%.6fHz\n%.3fdB" % (xyv[0], xyv[1])
 882         if self.lowpass >= 500:
 883             s = "%.6fHz\n%.3fdB" % (xyv[0]*1000, xyv[1])
 884
 885         self.myform['spec_data'].set_value(s)
 886
 887     # This is another callback used by ra_fftsink.py (passed on creation
 888     #   of ra_fftsink).  We pass this as our "calibrator" function, but
 889     #   we create interesting side-effects in the GUI.
 890     #
 891     # This function finds peaks in the FFT output data, and reports
 892     #  on them through the "Best" text object in the GUI
 893     #  It also computes the Harmonic Compliance Measure (HCM), and displays
 894     #  that also.
 895     #
 896     def pulsarfunc(self,d,l):
 897        x = range(0,l)
 898        incr = float(self.lowpass)/float(l)
 899        incr = incr * 2.0
 900        bestdb = -50.0
 901        bestfreq = 0.0
 902        avg = 0
 903        dcnt = 0
 904        #
 905        # First, we need to find the average signal level
 906        #
 907        for i in x:
 908            if (i * incr) > self.lowest_freq and (i*incr) < (self.lowpass-2):
 909                avg += d[i]
 910                dcnt += 1
 911        # Set average signal level
 912        avg /= dcnt
 913        s2=" "
 914        findcnt = 0
 915        #
 916        # Then we find candidates that are greater than the user-supplied
 917        #   threshold.
 918        #
 919        # We try to cluster "hits" whose whole-number frequency is the
 920        #   same, and compute an average "hit" frequency.
 921        #
 922        lastint = 0
 923        hits=[]
 924        intcnt = 0
 925        freqavg = 0
 926        for i in x:
 927            freq = i*incr
 928            # If frequency within bounds, and the (dB-avg) value is above our
 929            #   threshold
 930            if freq > self.lowest_freq and freq < self.lowpass-2 and (d[i] - avg) > self.threshold:
 931                # If we're finding a new whole-number frequency
 932                if lastint != int(freq):
 933                    # Record "center" of this hit, if this is a new hit
 934                    if lastint != 0:
 935                        s2 += "%5.3fHz " % (freqavg/intcnt)
 936                        hits.append(freqavg/intcnt)
 937                        findcnt += 1
 938                    lastint = int(freq)
 939                    intcnt = 1
 940                    freqavg = freq
 941                else:
 942                    intcnt += 1
 943                    freqavg += freq
 944            if (findcnt >= 14):
 945                break
 946
 947        if intcnt > 1:
 948            s2 += "%5.3fHz " % (freqavg/intcnt)
 949            hits.append(freqavg/intcnt)
 950
 951        #
 952        # Compute the HCM, by dividing each of the "hits" by each of the
 953        #   other hits, and comparing the difference between a "perfect"
 954        #   harmonic, and the observed frequency ratio.
 955        #
 956        measure = 0
 957        max_measure=0
 958        mcnt = 0
 959        avg_dist = 0
 960        acnt = 0
 961        for i in range(1,len(hits)):
 962            meas = hits[i]/hits[0] - int(hits[i]/hits[0])
 963            if abs((hits[i]-hits[i-1])-hits[0]) < 0.1:
 964                avg_dist += hits[i]-hits[i-1]
 965                acnt += 1
 966            if meas > 0.98 and meas < 1.0:
 967                meas = 1.0 - meas
 968            meas *= hits[0]
 969            if meas >= max_measure:
 970                max_measure = meas
 971            measure += meas
 972            mcnt += 1
 973        if mcnt > 0:
 974            measure /= mcnt
 975            if acnt > 0:
 976                avg_dist /= acnt
 977        if len(hits) > 1:
 978            measure /= mcnt
 979            s3="\nHCM: Avg %5.3fHz(%d) Max %5.3fHz Dist %5.3fHz(%d)" % (measure,mcnt,max_measure, avg_dist, acnt)
 980            if max_measure < 0.5 and len(hits) >= 2:
 981                self.record_hit(hits, measure, max_measure)
 982                self.avg_dist = avg_dist
 983        else:
 984            s3="\nHCM: --"
 985        s4="\nAvg dB: %4.2f" % avg
 986        self.myform['best_pulse'].set_value("("+s2+")"+s3+s4)
 987
 988        # Since we are nominally a calibrator function for ra_fftsink, we
 989        #  simply return what they sent us, untouched.  A "real" calibrator
 990        #  function could monkey with the data before returning it to the
 991        #  FFT display function.
 992        return(d)
 993
 994     #
 995     # Callback for the "DM" gui object
 996     #
 997     # We call compute_dispfilter() as appropriate to compute a new filter,
 998     #   and then set that new filter into self.dispfilt.
 999     #
1000     def set_dm(self,dm):
1001        self.dm = dm
1002
1003        ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
1004        self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
1005        self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
1006        self.dispfilt.set_taps(self.disp_taps)
1007        self.myform['DM'].set_value(dm)
1008        return(dm)
1009
1010     #
1011     # Callback for the "Doppler" gui object
1012     #
1013     # We call compute_dispfilter() as appropriate to compute a new filter,
1014     #   and then set that new filter into self.dispfilt.
1015     #
1016     def set_doppler(self,doppler):
1017        self.doppler = doppler
1018
1019        ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
1020        self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
1021        self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
1022        self.dispfilt.set_taps(self.disp_taps)
1023        self.myform['Doppler'].set_value(doppler)
1024        return(doppler)
1025
1026     #
1027     # Compute a de-dispersion filter
1028     #  From Hankins, et al, 1975
1029     #
1030     # This code translated from dedisp_filter.c from Swinburne
1031     #   pulsar software repository
1032     #
1033     def compute_dispfilter(self,dm,doppler,bw,centerfreq):
1034         npts = len(self.disp_taps)
1035         tmp = Numeric.zeros(npts, Numeric.Complex64)
1036         M_PI = 3.14159265358
1037         DM = dm/2.41e-10
1038
1039         #
1040         # Because astronomers are a crazy bunch, the "standard" calcultion
1041         #   is in Mhz, rather than Hz
1042         #
1043         centerfreq = centerfreq / 1.0e6
1044         bw = bw / 1.0e6
1045
1046         isign = int(bw / abs (bw))
1047
1048         # Center frequency may be doppler shifted
1049         cfreq     = centerfreq / doppler
1050
1051         # As well as the bandwidth..
1052         bandwidth = bw / doppler
1053
1054         # Bandwidth divided among bins
1055         binwidth  = bandwidth / npts
1056
1057         # Delay is an "extra" parameter, in usecs, and largely
1058         #  untested in the Swinburne code.
1059         delay = 0.0
1060
1061         # This determines the coefficient of the frequency response curve
1062         # Linear in DM, but quadratic in center frequency
1063         coeff = isign * 2.0*M_PI * DM / (cfreq*cfreq)
1064
1065         # DC to nyquist/2
1066         n = 0
1067         for i in range(0,int(npts/2)):
1068             freq = (n + 0.5) * binwidth
1069             phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1070             tmp[i] = complex(math.cos(phi), math.sin(phi))
1071             n += 1
1072
1073         # -nyquist/2 to DC
1074         n = int(npts/2)
1075         n *= -1
1076         for i in range(int(npts/2),npts):
1077             freq = (n + 0.5) * binwidth
1078             phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1079             tmp[i] = complex(math.cos(phi), math.sin(phi))
1080             n += 1
1081
1082         self.disp_taps = FFT.inverse_fft(tmp)
1083         return(self.disp_taps)
1084
1085     #
1086     # Compute minimum number of taps required in de-dispersion FFT filter
1087     #
1088     def compute_disp_ntaps(self,dm,bw,freq):
1089         #
1090         # Dt calculations are in Mhz, rather than Hz
1091         #    crazy astronomers....
1092         mbw = bw/1.0e6
1093         mfreq = freq/1.0e6
1094
1095         f_lower = mfreq-(mbw/2)
1096         f_upper = mfreq+(mbw/2)
1097
1098         # Compute smear time
1099         Dt = dm/2.41e-4 * (1.0/(f_lower*f_lower)-1.0/(f_upper*f_upper))
1100
1101         # ntaps is now bandwidth*smeartime
1102         # Should be bandwidth*smeartime*2, but the Gnu Radio FFT filter
1103         #   already expands it by a factor of 2
1104         ntaps = bw*Dt
1105         if ntaps < 64:
1106             ntaps = 64
1107         return(int(ntaps))
1108
1109 def main ():
1110     app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY PULSAR RECEIVER: $Revision: 3626 $", nstatus=1)
1111     app.MainLoop()
1112
1113 if __name__ == '__main__':
1114     main ()