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

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