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,2007 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   7 # GNU Radio is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 3, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23
  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, blks2, audio
  33 from usrpm 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 numpy.fft
  43 import ephem
  44 import time
  45 import os
  46 import math
  47
  48
  49 class app_flow_graph(stdgui2.std_top_block):
  50     def __init__(self, frame, panel, vbox, argv):
  51         stdgui2.std_top_block.__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_source
  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 (panel,
 296            fft_size=int(options.fft_size), sample_rate=PULSAR_MAX_FREQ*2,
 297            title="Post-detector spectrum",
 298            ofunc=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         print "input_rate ", second_input_rate, "audiodev ", self.audiodev
 346         self.audio = audio.sink(second_input_rate, self.audiodev)
 347
 348         #
 349         # The three post-detector filters
 350         # Done this way to allow an audio path (up to 10Khz)
 351         # ...and also because going from xMhz down to ~100Hz
 352         # In a single filter doesn't seem to work.
 353         #
 354         self.first = gr.fir_filter_fff (FIRST_FACTOR/2, first_filter)
 355
 356         p = second_input_rate / (PULSAR_MAX_FREQ*20)
 357         self.second = gr.fir_filter_fff (int(p), second_filter)
 358         self.third = gr.fir_filter_fff (10, third_filter)
 359
 360         # Detector
 361         self.detector = gr.complex_to_mag_squared()
 362
 363         self.enable_comb_filter = False
 364         # Epoch folder comb filter
 365         if self.enable_comb_filter == True:
 366             bogtaps = Numeric.zeros(512, Numeric.Float64)
 367             self.folder_comb = gr.fft_filter_ccc(1,bogtaps)
 368
 369         # Rational resampler
 370         self.folder_rr = blks2.rational_resampler_fff(self.interp, self.decim)
 371
 372         # Epoch folder bandpass
 373         bogtaps = Numeric.zeros(1, Numeric.Float64)
 374         self.folder_bandpass = gr.fir_filter_fff (1, bogtaps)
 375
 376         # Epoch folder F2C/C2F
 377         self.folder_f2c = gr.float_to_complex()
 378         self.folder_c2f = gr.complex_to_float()
 379
 380         # Epoch folder S2P
 381         self.folder_s2p = gr.serial_to_parallel (gr.sizeof_float,
 382              self.folding*FOLD_MULT)
 383
 384         # Epoch folder IIR Filter (produces average pulse profiles)
 385         self.folder_iir = gr.single_pole_iir_filter_ff(1.0/options.favg,
 386              self.folding*FOLD_MULT)
 387
 388         #
 389         # Set all the epoch-folder goop up
 390         #
 391         self.set_folding_params()
 392
 393         #
 394         # Start connecting configured modules in the receive chain
 395         #
 396
 397         # Connect raw USRP to de-dispersion filter, detector
 398         self.connect(self.u, self.dispfilt, self.detector)
 399
 400         # Connect detector output to FIR LPF
 401         #  in two stages, followed by the FFT scope
 402         self.connect(self.detector, self.first,
 403             self.second, self.third, self.scope)
 404
 405         # Connect audio output
 406         self.connect(self.first, self.volume)
 407         self.connect(self.volume, (self.audio, 0))
 408         self.connect(self.volume, (self.audio, 1))
 409
 410         # Connect epoch folder
 411         if self.enable_comb_filter == True:
 412             self.connect (self.first, self.folder_bandpass, self.folder_rr,
 413                 self.folder_f2c,
 414                 self.folder_comb, self.folder_c2f,
 415                 self.folder_s2p, self.folder_iir,
 416                 self.chart)
 417
 418         else:
 419             self.connect (self.first, self.folder_bandpass, self.folder_rr,
 420                 self.folder_s2p, self.folder_iir, self.chart)
 421
 422         # Connect baseband recording file (initially /dev/null)
 423         self.connect(self.u, self.tofloat, self.tochar, self.recording)
 424
 425         # Connect pulse recording file (initially /dev/null)
 426         self.connect(self.first, self.toshort, self.pulse_recording)
 427
 428         #
 429         # Build the GUI elements
 430         #
 431         self._build_gui(vbox)
 432
 433         # Make GUI agree with command-line
 434         self.myform['average'].set_value(int(options.avg))
 435         self.myform['foldavg'].set_value(int(options.favg))
 436
 437
 438         # Make spectral averager agree with command line
 439         if options.avg != 1.0:
 440             self.scope.set_avg_alpha(float(1.0/options.avg))
 441             self.scope.set_average(True)
 442
 443
 444         # set initial values
 445
 446         if options.gain is None:
 447             # if no gain was specified, use the mid-point in dB
 448             g = self.subdev.gain_range()
 449             options.gain = float(g[0]+g[1])/2
 450
 451         if options.freq is None:
 452             # if no freq was specified, use the mid-point
 453             r = self.subdev.freq_range()
 454             options.freq = float(r[0]+r[1])/2
 455
 456         self.set_gain(options.gain)
 457         self.set_volume(-10.0)
 458
 459         if not(self.set_freq(options.freq)):
 460             self._set_status_msg("Failed to set initial frequency")
 461
 462         self.myform['decim'].set_value(self.u.decim_rate())
 463         self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
 464         self.myform['dbname'].set_value(self.subdev.name())
 465         self.myform['DM'].set_value(self.dm)
 466         self.myform['Doppler'].set_value(self.doppler)
 467
 468         #
 469         # Start the timer that shows current LMST on the GUI
 470         #
 471         self.lmst_timer.Start(1000)
 472
 473
 474     def _set_status_msg(self, msg):
 475         self.frame.GetStatusBar().SetStatusText(msg, 0)
 476
 477     def _build_gui(self, vbox):
 478
 479         def _form_set_freq(kv):
 480             return self.set_freq(kv['freq'])
 481
 482         def _form_set_dm(kv):
 483             return self.set_dm(kv['DM'])
 484
 485         def _form_set_doppler(kv):
 486             return self.set_doppler(kv['Doppler'])
 487
 488         # Position the FFT or Waterfall
 489         vbox.Add(self.scope.win, 5, wx.EXPAND)
 490         vbox.Add(self.chart.win, 5, wx.EXPAND)
 491
 492         # add control area at the bottom
 493         self.myform = myform = form.form()
 494         hbox = wx.BoxSizer(wx.HORIZONTAL)
 495         hbox.Add((7,0), 0, wx.EXPAND)
 496         vbox1 = wx.BoxSizer(wx.VERTICAL)
 497         myform['freq'] = form.float_field(
 498             parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
 499             callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
 500
 501         vbox1.Add((3,0), 0, 0)
 502
 503         # To show current Local Mean Sidereal Time
 504         myform['lmst_high'] = form.static_text_field(
 505             parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
 506         vbox1.Add((3,0), 0, 0)
 507
 508         # To show current spectral cursor data
 509         myform['spec_data'] = form.static_text_field(
 510             parent=self.panel, sizer=vbox1, label="Pulse Freq", weight=1)
 511         vbox1.Add((3,0), 0, 0)
 512
 513         # To show best pulses found in FFT output
 514         myform['best_pulse'] = form.static_text_field(
 515             parent=self.panel, sizer=vbox1, label="Best freq", weight=1)
 516         vbox1.Add((3,0), 0, 0)
 517
 518         vboxBogus = wx.BoxSizer(wx.VERTICAL)
 519         vboxBogus.Add ((2,0), 0, wx.EXPAND)
 520         vbox2 = wx.BoxSizer(wx.VERTICAL)
 521         g = self.subdev.gain_range()
 522         myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
 523                                            weight=1,
 524                                            min=int(g[0]), max=int(g[1]),
 525                                            callback=self.set_gain)
 526
 527         vbox2.Add((6,0), 0, 0)
 528         myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
 529                     label="Spectral Averaging", weight=1, min=1, max=200, callback=self.set_averaging)
 530         vbox2.Add((6,0), 0, 0)
 531         myform['foldavg'] = form.slider_field(parent=self.panel, sizer=vbox2,
 532                     label="Folder Averaging", weight=1, min=1, max=20, callback=self.set_folder_averaging)
 533         vbox2.Add((6,0), 0, 0)
 534         myform['volume'] = form.quantized_slider_field(parent=self.panel, sizer=vbox2,
 535                     label="Audio Volume", weight=1, range=(-20, 0, 0.5), callback=self.set_volume)
 536         vbox2.Add((6,0), 0, 0)
 537         myform['DM'] = form.float_field(
 538             parent=self.panel, sizer=vbox2, label="DM", weight=1,
 539             callback=myform.check_input_and_call(_form_set_dm))
 540         vbox2.Add((6,0), 0, 0)
 541         myform['Doppler'] = form.float_field(
 542             parent=self.panel, sizer=vbox2, label="Doppler", weight=1,
 543             callback=myform.check_input_and_call(_form_set_doppler))
 544         vbox2.Add((6,0), 0, 0)
 545
 546
 547         # Baseband recording control
 548         buttonbox = wx.BoxSizer(wx.HORIZONTAL)
 549         self.record_control = form.button_with_callback(self.panel,
 550               label="Recording baseband: Off                           ",
 551               callback=self.toggle_recording)
 552         self.record_pulse_control = form.button_with_callback(self.panel,
 553               label="Recording pulses: Off                              ",
 554               callback=self.toggle_pulse_recording)
 555
 556         buttonbox.Add(self.record_control, 0, wx.CENTER)
 557         buttonbox.Add(self.record_pulse_control, 0, wx.CENTER)
 558         vbox.Add(buttonbox, 0, wx.CENTER)
 559         hbox.Add(vbox1, 0, 0)
 560         hbox.Add(vboxBogus, 0, 0)
 561         hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
 562         vbox.Add(hbox, 0, wx.EXPAND)
 563
 564         self._build_subpanel(vbox)
 565
 566         self.lmst_timer = wx.PyTimer(self.lmst_timeout)
 567         self.lmst_timeout()
 568
 569
 570     def _build_subpanel(self, vbox_arg):
 571         # build a secondary information panel (sometimes hidden)
 572
 573         # FIXME figure out how to have this be a subpanel that is always
 574         # created, but has its visibility controlled by foo.Show(True/False)
 575
 576         if not(self.show_debug_info):
 577             return
 578
 579         panel = self.panel
 580         vbox = vbox_arg
 581         myform = self.myform
 582
 583         #panel = wx.Panel(self.panel, -1)
 584         #vbox = wx.BoxSizer(wx.VERTICAL)
 585
 586         hbox = wx.BoxSizer(wx.HORIZONTAL)
 587         hbox.Add((5,0), 0)
 588         myform['decim'] = form.static_float_field(
 589             parent=panel, sizer=hbox, label="Decim")
 590
 591         hbox.Add((5,0), 1)
 592         myform['fs@usb'] = form.static_float_field(
 593             parent=panel, sizer=hbox, label="Fs@USB")
 594
 595         hbox.Add((5,0), 1)
 596         myform['dbname'] = form.static_text_field(
 597             parent=panel, sizer=hbox)
 598
 599         hbox.Add((5,0), 1)
 600         myform['baseband'] = form.static_float_field(
 601             parent=panel, sizer=hbox, label="Analog BB")
 602
 603         hbox.Add((5,0), 1)
 604         myform['ddc'] = form.static_float_field(
 605             parent=panel, sizer=hbox, label="DDC")
 606
 607         hbox.Add((5,0), 0)
 608         vbox.Add(hbox, 0, wx.EXPAND)
 609
 610
 611
 612     def set_freq(self, target_freq):
 613         """
 614         Set the center frequency we're interested in.
 615
 616         @param target_freq: frequency in Hz
 617         @rypte: bool
 618
 619         Tuning is a two step process.  First we ask the front-end to
 620         tune as close to the desired frequency as it can.  Then we use
 621         the result of that operation and our target_frequency to
 622         determine the value for the digital down converter.
 623         """
 624         r = usrp.tune(self.u, 0, self.subdev, target_freq)
 625
 626         if r:
 627             self.myform['freq'].set_value(target_freq)     # update displayed value
 628             self.myform['baseband'].set_value(r.baseband_freq)
 629             self.myform['ddc'].set_value(r.dxc_freq)
 630             # Adjust self.frequency, and self.observing_freq
 631             # We pick up the difference between the current self.frequency
 632             #   and the just-programmed one, and use this to adjust
 633             #   self.observing_freq.  We have to do it this way to
 634             #   make the dedispersion filtering work out properly.
 635             delta = target_freq - self.frequency
 636             self.frequency = target_freq
 637             self.observing_freq += delta
 638
 639             # Now that we're adjusted, compute a new dispfilter, and
 640             #   set the taps for the FFT filter.
 641             ntaps = self.compute_disp_ntaps(self.dm, self.bw, self.observing_freq)
 642             self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
 643             self.compute_dispfilter(self.dm,self.doppler,self.bw,
 644                 self.observing_freq)
 645             self.dispfilt.set_taps(self.disp_taps)
 646
 647             return True
 648
 649         return False
 650
 651     # Callback for gain-setting slider
 652     def set_gain(self, gain):
 653         self.myform['gain'].set_value(gain)     # update displayed value
 654         self.subdev.set_gain(gain)
 655
 656
 657     def set_volume(self, vol):
 658         self.myform['volume'].set_value(vol)
 659         self.volume.set_k((10**(vol/10))/8192)
 660
 661     # Callback for spectral-averaging slider
 662     def set_averaging(self, avval):
 663         self.myform['average'].set_value(avval)
 664         self.scope.set_avg_alpha(1.0/(avval))
 665         self.scope.set_average(True)
 666
 667     def set_folder_averaging(self, avval):
 668         self.myform['foldavg'].set_value(avval)
 669         self.folder_iir.set_taps(1.0/avval)
 670
 671     # Timer callback to update LMST display
 672     def lmst_timeout(self):
 673          self.locality.date = ephem.now()
 674          sidtime = self.locality.sidereal_time()
 675          self.myform['lmst_high'].set_value(str(ephem.hours(sidtime)))
 676
 677     #
 678     # Turn recording on/off
 679     # Called-back by "Recording" button
 680     #
 681     def toggle_recording(self):
 682         # Pick up current LMST
 683         self.locality.date = ephem.now()
 684         sidtime = self.locality.sidereal_time()
 685
 686         # Pick up localtime, for generating filenames
 687         foo = time.localtime()
 688
 689         # Generate filenames for both data and header file
 690         filename = "%04d%02d%02d%02d%02d.pdat" % (foo.tm_year, foo.tm_mon,
 691            foo.tm_mday, foo.tm_hour, foo.tm_min)
 692         hdrfilename = "%04d%02d%02d%02d%02d.phdr" % (foo.tm_year, foo.tm_mon,
 693            foo.tm_mday, foo.tm_hour, foo.tm_min)
 694
 695         # Current recording?  Flip state
 696         if (self.recording_state == True):
 697           self.recording_state = False
 698           self.record_control.SetLabel("Recording baseband: Off                           ")
 699           self.recording.close()
 700         # Not recording?
 701         else:
 702           self.recording_state = True
 703           self.record_control.SetLabel("Recording baseband to: "+filename)
 704
 705           # Cause gr_file_sink object to accept new filename
 706           #   note use of self.prefix--filename prefix from
 707           #   command line (defaults to ./)
 708           #
 709           self.recording.open (self.prefix+filename)
 710
 711           #
 712           # We open the header file as a regular file, write header data,
 713           #   then close
 714           hdrf = open(self.prefix+hdrfilename, "w")
 715           hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
 716           hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
 717           hdrf.write("DM: "+str(self.dm)+"\n")
 718           hdrf.write("doppler: "+str(self.doppler)+"\n")
 719
 720           hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 721           hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 722           hdrf.write("sample type: complex_char\n")
 723           hdrf.write("sample size: "+str(gr.sizeof_char*2)+"\n")
 724           hdrf.close()
 725     #
 726     # Turn recording on/off
 727     # Called-back by "Recording" button
 728     #
 729     def toggle_pulse_recording(self):
 730         # Pick up current LMST
 731         self.locality.date = ephem.now()
 732         sidtime = self.locality.sidereal_time()
 733
 734         # Pick up localtime, for generating filenames
 735         foo = time.localtime()
 736
 737         # Generate filenames for both data and header file
 738         filename = "%04d%02d%02d%02d%02d.padat" % (foo.tm_year, foo.tm_mon,
 739            foo.tm_mday, foo.tm_hour, foo.tm_min)
 740         hdrfilename = "%04d%02d%02d%02d%02d.pahdr" % (foo.tm_year, foo.tm_mon,
 741            foo.tm_mday, foo.tm_hour, foo.tm_min)
 742
 743         # Current recording?  Flip state
 744         if (self.pulse_recording_state == True):
 745           self.pulse_recording_state = False
 746           self.record_pulse_control.SetLabel("Recording pulses: Off                           ")
 747           self.pulse_recording.close()
 748         # Not recording?
 749         else:
 750           self.pulse_recording_state = True
 751           self.record_pulse_control.SetLabel("Recording pulses to: "+filename)
 752
 753           # Cause gr_file_sink object to accept new filename
 754           #   note use of self.prefix--filename prefix from
 755           #   command line (defaults to ./)
 756           #
 757           self.pulse_recording.open (self.prefix+filename)
 758
 759           #
 760           # We open the header file as a regular file, write header data,
 761           #   then close
 762           hdrf = open(self.prefix+hdrfilename, "w")
 763           hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
 764           hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
 765           hdrf.write("DM: "+str(self.dm)+"\n")
 766           hdrf.write("doppler: "+str(self.doppler)+"\n")
 767           hdrf.write("pulse rate: "+str(self.pulse_freq)+"\n")
 768           hdrf.write("pulse sps: "+str(self.pulse_freq*self.folding)+"\n")
 769           hdrf.write("file sps: "+str(self.folder_input_rate)+"\n")
 770
 771           hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 772           hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 773           hdrf.write("sample type: short\n")
 774           hdrf.write("sample size: 1\n")
 775           hdrf.close()
 776
 777     # We get called at startup, and whenever the GUI "Set Folding params"
 778     #   button is pressed
 779     #
 780     def set_folding_params(self):
 781         if (self.pulse_freq <= 0):
 782             return
 783
 784         # Compute required sample rate
 785         self.sample_rate = int(self.pulse_freq*self.folding)
 786
 787         # And the implied decimation rate
 788         required_decimation = int(self.folder_input_rate / self.sample_rate)
 789
 790         # We also compute a new FFT comb filter, based on the expected
 791         #  spectral profile of our pulse parameters
 792         #
 793         # FFT-based comb filter
 794         #
 795         N_COMB_TAPS=int(self.sample_rate*4)
 796         if N_COMB_TAPS > 2000:
 797             N_COMB_TAPS = 2000
 798         self.folder_comb_taps = Numeric.zeros(N_COMB_TAPS,Numeric.Complex64)
 799         fincr = (self.sample_rate)/float(N_COMB_TAPS)
 800         for i in range(0,len(self.folder_comb_taps)):
 801             self.folder_comb_taps[i] = complex(0.0, 0.0)
 802
 803         freq = 0.0
 804         harmonics = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
 805         for i in range(0,len(self.folder_comb_taps)/2):
 806             for j in range(0,len(harmonics)):
 807                  if abs(freq - harmonics[j]*self.pulse_freq) <= fincr:
 808                      self.folder_comb_taps[i] = complex(4.0, 0.0)
 809                      if harmonics[j] == 1.0:
 810                          self.folder_comb_taps[i] = complex(8.0, 0.0)
 811             freq += fincr
 812
 813         if self.enable_comb_filter == True:
 814             # Set the just-computed FFT comb filter taps
 815             self.folder_comb.set_taps(self.folder_comb_taps)
 816
 817         # And compute a new decimated bandpass filter, to go in front
 818         #   of the comb.  Primary function is to decimate and filter down
 819         #   to an exact-ish multiple of the target pulse rate
 820         #
 821         self.folding_taps = gr.firdes_band_pass (1.0, self.folder_input_rate,
 822             0.10, self.sample_rate/2, 10,
 823             gr.firdes.WIN_HAMMING)
 824
 825         # Set the computed taps for the bandpass/decimate filter
 826         self.folder_bandpass.set_taps (self.folding_taps)
 827     #
 828     # Record a spectral "hit" of a possible pulsar spectral profile
 829     #
 830     def record_hit(self,hits, hcavg, hcmax):
 831         # Pick up current LMST
 832         self.locality.date = ephem.now()
 833         sidtime = self.locality.sidereal_time()
 834
 835         # Pick up localtime, for generating filenames
 836         foo = time.localtime()
 837
 838         # Generate filenames for both data and header file
 839         hitfilename = "%04d%02d%02d%02d.phit" % (foo.tm_year, foo.tm_mon,
 840            foo.tm_mday, foo.tm_hour)
 841
 842         hitf = open(self.prefix+hitfilename, "a")
 843         hitf.write("receiver center frequency: "+str(self.frequency)+"\n")
 844         hitf.write("observing frequency: "+str(self.observing_freq)+"\n")
 845         hitf.write("DM: "+str(self.dm)+"\n")
 846         hitf.write("doppler: "+str(self.doppler)+"\n")
 847
 848         hitf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
 849         hitf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
 850         hitf.write("spectral peaks: "+str(hits)+"\n")
 851         hitf.write("HCM: "+str(hcavg)+" "+str(hcmax)+"\n")
 852         hitf.close()
 853
 854     # This is a callback used by ra_fftsink.py (passed on creation of
 855     #   ra_fftsink)
 856     # Whenever the user moves the cursor within the FFT display, this
 857     #   shows the coordinate data
 858     #
 859     def xydfunc(self,xyv):
 860         s = "%.6fHz\n%.3fdB" % (xyv[0], xyv[1])
 861         if self.lowpass >= 500:
 862             s = "%.6fHz\n%.3fdB" % (xyv[0]*1000, xyv[1])
 863
 864         self.myform['spec_data'].set_value(s)
 865
 866     # This is another callback used by ra_fftsink.py (passed on creation
 867     #   of ra_fftsink).  We pass this as our "calibrator" function, but
 868     #   we create interesting side-effects in the GUI.
 869     #
 870     # This function finds peaks in the FFT output data, and reports
 871     #  on them through the "Best" text object in the GUI
 872     #  It also computes the Harmonic Compliance Measure (HCM), and displays
 873     #  that also.
 874     #
 875     def pulsarfunc(self,d,l):
 876        x = range(0,l)
 877        incr = float(self.lowpass)/float(l)
 878        incr = incr * 2.0
 879        bestdb = -50.0
 880        bestfreq = 0.0
 881        avg = 0
 882        dcnt = 0
 883        #
 884        # First, we need to find the average signal level
 885        #
 886        for i in x:
 887            if (i * incr) > self.lowest_freq and (i*incr) < (self.lowpass-2):
 888                avg += d[i]
 889                dcnt += 1
 890        # Set average signal level
 891        avg /= dcnt
 892        s2=" "
 893        findcnt = 0
 894        #
 895        # Then we find candidates that are greater than the user-supplied
 896        #   threshold.
 897        #
 898        # We try to cluster "hits" whose whole-number frequency is the
 899        #   same, and compute an average "hit" frequency.
 900        #
 901        lastint = 0
 902        hits=[]
 903        intcnt = 0
 904        freqavg = 0
 905        for i in x:
 906            freq = i*incr
 907            # If frequency within bounds, and the (dB-avg) value is above our
 908            #   threshold
 909            if freq > self.lowest_freq and freq < self.lowpass-2 and (d[i] - avg) > self.threshold:
 910                # If we're finding a new whole-number frequency
 911                if lastint != int(freq):
 912                    # Record "center" of this hit, if this is a new hit
 913                    if lastint != 0:
 914                        s2 += "%5.3fHz " % (freqavg/intcnt)
 915                        hits.append(freqavg/intcnt)
 916                        findcnt += 1
 917                    lastint = int(freq)
 918                    intcnt = 1
 919                    freqavg = freq
 920                else:
 921                    intcnt += 1
 922                    freqavg += freq
 923            if (findcnt >= 14):
 924                break
 925
 926        if intcnt > 1:
 927            s2 += "%5.3fHz " % (freqavg/intcnt)
 928            hits.append(freqavg/intcnt)
 929
 930        #
 931        # Compute the HCM, by dividing each of the "hits" by each of the
 932        #   other hits, and comparing the difference between a "perfect"
 933        #   harmonic, and the observed frequency ratio.
 934        #
 935        measure = 0
 936        max_measure=0
 937        mcnt = 0
 938        avg_dist = 0
 939        acnt = 0
 940        for i in range(1,len(hits)):
 941            meas = hits[i]/hits[0] - int(hits[i]/hits[0])
 942            if abs((hits[i]-hits[i-1])-hits[0]) < 0.1:
 943                avg_dist += hits[i]-hits[i-1]
 944                acnt += 1
 945            if meas > 0.98 and meas < 1.0:
 946                meas = 1.0 - meas
 947            meas *= hits[0]
 948            if meas >= max_measure:
 949                max_measure = meas
 950            measure += meas
 951            mcnt += 1
 952        if mcnt > 0:
 953            measure /= mcnt
 954            if acnt > 0:
 955                avg_dist /= acnt
 956        if len(hits) > 1:
 957            measure /= mcnt
 958            s3="\nHCM: Avg %5.3fHz(%d) Max %5.3fHz Dist %5.3fHz(%d)" % (measure,mcnt,max_measure, avg_dist, acnt)
 959            if max_measure < 0.5 and len(hits) >= 2:
 960                self.record_hit(hits, measure, max_measure)
 961                self.avg_dist = avg_dist
 962        else:
 963            s3="\nHCM: --"
 964        s4="\nAvg dB: %4.2f" % avg
 965        self.myform['best_pulse'].set_value("("+s2+")"+s3+s4)
 966
 967        # Since we are nominally a calibrator function for ra_fftsink, we
 968        #  simply return what they sent us, untouched.  A "real" calibrator
 969        #  function could monkey with the data before returning it to the
 970        #  FFT display function.
 971        return(d)
 972
 973     #
 974     # Callback for the "DM" gui object
 975     #
 976     # We call compute_dispfilter() as appropriate to compute a new filter,
 977     #   and then set that new filter into self.dispfilt.
 978     #
 979     def set_dm(self,dm):
 980        self.dm = dm
 981
 982        ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
 983        self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
 984        self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
 985        self.dispfilt.set_taps(self.disp_taps)
 986        self.myform['DM'].set_value(dm)
 987        return(dm)
 988
 989     #
 990     # Callback for the "Doppler" gui object
 991     #
 992     # We call compute_dispfilter() as appropriate to compute a new filter,
 993     #   and then set that new filter into self.dispfilt.
 994     #
 995     def set_doppler(self,doppler):
 996        self.doppler = doppler
 997
 998        ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
 999        self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
1000        self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
1001        self.dispfilt.set_taps(self.disp_taps)
1002        self.myform['Doppler'].set_value(doppler)
1003        return(doppler)
1004
1005     #
1006     # Compute a de-dispersion filter
1007     #  From Hankins, et al, 1975
1008     #
1009     # This code translated from dedisp_filter.c from Swinburne
1010     #   pulsar software repository
1011     #
1012     def compute_dispfilter(self,dm,doppler,bw,centerfreq):
1013         npts = len(self.disp_taps)
1014         tmp = Numeric.zeros(npts, Numeric.Complex64)
1015         M_PI = 3.14159265358
1016         DM = dm/2.41e-10
1017
1018         #
1019         # Because astronomers are a crazy bunch, the "standard" calcultion
1020         #   is in Mhz, rather than Hz
1021         #
1022         centerfreq = centerfreq / 1.0e6
1023         bw = bw / 1.0e6
1024
1025         isign = int(bw / abs (bw))
1026
1027         # Center frequency may be doppler shifted
1028         cfreq     = centerfreq / doppler
1029
1030         # As well as the bandwidth..
1031         bandwidth = bw / doppler
1032
1033         # Bandwidth divided among bins
1034         binwidth  = bandwidth / npts
1035
1036         # Delay is an "extra" parameter, in usecs, and largely
1037         #  untested in the Swinburne code.
1038         delay = 0.0
1039
1040         # This determines the coefficient of the frequency response curve
1041         # Linear in DM, but quadratic in center frequency
1042         coeff = isign * 2.0*M_PI * DM / (cfreq*cfreq)
1043
1044         # DC to nyquist/2
1045         n = 0
1046         for i in range(0,int(npts/2)):
1047             freq = (n + 0.5) * binwidth
1048             phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1049             tmp[i] = complex(math.cos(phi), math.sin(phi))
1050             n += 1
1051
1052         # -nyquist/2 to DC
1053         n = int(npts/2)
1054         n *= -1
1055         for i in range(int(npts/2),npts):
1056             freq = (n + 0.5) * binwidth
1057             phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1058             tmp[i] = complex(math.cos(phi), math.sin(phi))
1059             n += 1
1060
1061         self.disp_taps = numpy.fft.ifft(tmp)
1062         return(self.disp_taps)
1063
1064     #
1065     # Compute minimum number of taps required in de-dispersion FFT filter
1066     #
1067     def compute_disp_ntaps(self,dm,bw,freq):
1068         #
1069         # Dt calculations are in Mhz, rather than Hz
1070         #    crazy astronomers....
1071         mbw = bw/1.0e6
1072         mfreq = freq/1.0e6
1073
1074         f_lower = mfreq-(mbw/2)
1075         f_upper = mfreq+(mbw/2)
1076
1077         # Compute smear time
1078         Dt = dm/2.41e-4 * (1.0/(f_lower*f_lower)-1.0/(f_upper*f_upper))
1079
1080         # ntaps is now bandwidth*smeartime
1081         # Should be bandwidth*smeartime*2, but the Gnu Radio FFT filter
1082         #   already expands it by a factor of 2
1083         ntaps = bw*Dt
1084         if ntaps < 64:
1085             ntaps = 64
1086         return(int(ntaps))
1087
1088 def main ():
1089     app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY PULSAR RECEIVER: $Revision$", nstatus=1)
1090     app.MainLoop()
1091
1092 if __name__ == '__main__':
1093     main ()