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

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