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