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