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