3 # Copyright 2004,2005 Free Software Foundation, Inc.
5 # This file is part of GNU Radio
7 # GNU Radio is free software; you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation; either version 2, or (at your option)
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.
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.
26 # Pulsar receiver application
28 # Performs both harmonic folding analysis
29 # and epoch folding analysis
32 from gnuradio import gr, gru, blks, audio
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
49 class app_flow_graph(stdgui.gui_flow_graph):
50 def __init__(self, frame, panel, vbox, argv):
51 stdgui.gui_flow_graph.__init__(self)
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")
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()
93 self.show_debug_info = True
95 self.reflevel = options.reflevel
96 self.divbase = options.divbase
97 self.division = options.division
98 self.audiodev = options.audio_dev
100 # Low-pass cutoff for post-detector filter
101 # Set to 100Hz usually, since lots of pulsars fit in this
103 self.lowpass = options.lowpass
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
109 # What (dB) threshold to use in determining spectral candidates
110 self.threshold = options.threshold
112 # Filename prefix for recording file
113 self.prefix = options.prefix
115 # Dispersion Measure (DM)
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
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.
130 self.frequency = options.freq
131 if options.observing <= 0:
132 self.observing_freq = options.freq
134 self.observing_freq = options.observing
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))
141 # Recording file, in case we ever need to record baseband data
143 self.recording = gr.file_sink(gr.sizeof_char, "/dev/null")
144 self.recording_state = False
146 self.pulse_recording = gr.file_sink(gr.sizeof_short, "/dev/null")
147 self.pulse_recording_state = False
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()
156 # Need these two for converting 12-bit baseband signals to 8-bit
158 self.tofloat = gr.complex_to_float()
159 self.tochar = gr.float_to_char()
161 # Need this for recording pulses (post-detector)
162 self.toshort = gr.float_to_short()
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
170 # We set a sentimental value here
171 self.pulse_freq = options.pulsefreq
173 # Folder runs at this raw sample rate
174 self.folder_input_rate = 20000
176 # Each pulse in the epoch folder is sampled at 128 times the nominal
182 # Try to find candidate parameters for rational resampler
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):
194 # We didn't find anything, bail!
195 if (len(candidates) < 1):
196 print "Couldn't converge on resampler parameters"
200 # Now try to find candidate with the least sampling error
204 diff = self.pulse_freq * i
205 diff = diff - int(diff)
211 input_rate = self.folder_input_rate
212 output_rate = int(self.pulse_freq * save_i)
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
218 # Save optimized folding parameters, used later
219 self.folding = save_i
220 self.interp = int(interp)
221 self.decim = int(decim)
223 # So that we can view 4 pulses in the pulse viewer window
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)
230 # Compute raw input rate
231 input_rate = self.u.adc_freq() / self.u.decim_rate()
233 # BW==input_rate for complex data
237 # Set baseband filter bandwidth if DBS_RX:
239 if self.cardtype == usrp_dbid.DBS_RX:
243 self.subdev.set_bw(lbw)
246 # We use this as a crude volume control for the audio output
248 self.volume = gr.multiply_const_ff(10**(-1))
252 # Create location data for ephem package
254 self.locality = ephem.Observer()
255 self.locality.long = str(options.longitude)
256 self.locality.lat = str(options.latitude)
259 # What is the post-detector LPF cutoff for the FFT?
261 PULSAR_MAX_FREQ=int(options.lowpass)
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,
268 input_rate/FIRST_FACTOR,
269 input_rate/(FIRST_FACTOR*20),
270 gr.firdes.WIN_HAMMING)
272 # Second filter runs at the output rate of the first filter,
273 # And low-pass filters down to PULSAR_MAX_FREQ*10
275 second_input_rate = int(input_rate/(FIRST_FACTOR/2))
276 second_filter = gr.firdes.band_pass(1.0, second_input_rate,
280 gr.firdes.WIN_HAMMING)
282 # Third filter runs at PULSAR_MAX_FREQ*20
283 # and filters down to PULSAR_MAX_FREQ
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)
293 # Create the appropriate FFT scope
295 self.scope = ra_fftsink.ra_fft_sink_f (self, panel,
296 fft_size=int(options.fft_size), sample_rate=PULSAR_MAX_FREQ*2,
297 title="Post-detector spectrum",
298 cfunc=self.pulsarfunc, xydfunc=self.xydfunc, fft_rate=200)
301 # Tell scope we're looking from DC to PULSAR_MAX_FREQ
303 self.scope.set_baseband_freq (0.0)
307 # Setup stripchart for showing pulse profiles
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,
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)
320 # De-dispersion filter setup
322 # Do this here, just before creating the filter
323 # that will use the taps.
325 ntaps = self.compute_disp_ntaps(self.dm,self.bw,self.observing_freq)
327 # Taps for the de-dispersion filter
328 self.disp_taps = Numeric.zeros(ntaps,Numeric.Complex64)
330 # Compute the de-dispersion filter now
331 self.compute_dispfilter(self.dm,self.doppler,
332 self.bw,self.observing_freq)
335 # Call constructors for receive chains
339 # Now create the FFT filter using the computed taps
340 self.dispfilt = gr.fft_filter_ccc(1, self.disp_taps)
345 self.audio = audio.sink(second_input_rate, self.audiodev)
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.
353 self.first = gr.fir_filter_fff (FIRST_FACTOR/2, first_filter)
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)
359 # Split complex USRP stream into a pair of floats
360 self.splitter = gr.complex_to_float (1);
362 # I squarer (detector)
363 self.multI = gr.multiply_ff();
365 # Q squarer (detector)
366 self.multQ = gr.multiply_ff();
368 # Adding squared I and Q to produce instantaneous signal power
369 self.adder = gr.add_ff();
371 self.enable_comb_filter = False
372 # Epoch folder comb filter
373 if self.enable_comb_filter == True:
374 bogtaps = Numeric.zeros(512, Numeric.Float64)
375 self.folder_comb = gr.fft_filter_ccc(1,bogtaps)
378 self.folder_rr = blks.rational_resampler_fff(self, self.interp, self.decim)
380 # Epoch folder bandpass
381 bogtaps = Numeric.zeros(1, Numeric.Float64)
382 self.folder_bandpass = gr.fir_filter_fff (1, bogtaps)
384 # Epoch folder F2C/C2F
385 self.folder_f2c = gr.float_to_complex()
386 self.folder_c2f = gr.complex_to_float()
389 self.folder_s2p = gr.serial_to_parallel (gr.sizeof_float,
390 self.folding*FOLD_MULT)
392 # Epoch folder IIR Filter (produces average pulse profiles)
393 self.folder_iir = gr.single_pole_iir_filter_ff(1.0/options.favg,
394 self.folding*FOLD_MULT)
397 # Set all the epoch-folder goop up
399 self.set_folding_params()
402 # Start connecting configured modules in the receive chain
405 # Connect raw USRP to de-dispersion filter, complex->float splitter
406 self.connect(self.u, self.dispfilt, self.splitter)
408 # Connect splitter outputs to multipliers
410 self.connect((self.splitter, 0), (self.multI,0))
411 self.connect((self.splitter, 0), (self.multI,1))
414 self.connect((self.splitter, 1), (self.multQ,0))
415 self.connect((self.splitter, 1), (self.multQ,1))
417 # Then sum the squares
418 self.connect(self.multI, (self.adder,0))
419 self.connect(self.multQ, (self.adder,1))
421 # Connect detector/adder output to FIR LPF
422 # in two stages, followed by the FFT scope
423 self.connect(self.adder, self.first,
424 self.second, self.third, self.scope)
426 # Connect audio output
427 self.connect(self.first, self.volume)
428 self.connect(self.volume, (self.audio, 0))
429 self.connect(self.volume, (self.audio, 1))
431 # Connect epoch folder
432 if self.enable_comb_filter == True:
433 self.connect (self.first, self.folder_bandpass, self.folder_rr,
435 self.folder_comb, self.folder_c2f,
436 self.folder_s2p, self.folder_iir,
440 self.connect (self.first, self.folder_bandpass, self.folder_rr,
441 self.folder_s2p, self.folder_iir, self.chart)
443 # Connect baseband recording file (initially /dev/null)
444 self.connect(self.u, self.tofloat, self.tochar, self.recording)
446 # Connect pulse recording file (initially /dev/null)
447 self.connect(self.first, self.toshort, self.pulse_recording)
450 # Build the GUI elements
452 self._build_gui(vbox)
454 # Make GUI agree with command-line
455 self.myform['average'].set_value(int(options.avg))
456 self.myform['foldavg'].set_value(int(options.favg))
459 # Make spectral averager agree with command line
460 if options.avg != 1.0:
461 self.scope.set_avg_alpha(float(1.0/options.avg))
462 self.scope.set_average(True)
467 if options.gain is None:
468 # if no gain was specified, use the mid-point in dB
469 g = self.subdev.gain_range()
470 options.gain = float(g[0]+g[1])/2
472 if options.freq is None:
473 # if no freq was specified, use the mid-point
474 r = self.subdev.freq_range()
475 options.freq = float(r[0]+r[1])/2
477 self.set_gain(options.gain)
478 self.set_volume(-10.0)
480 if not(self.set_freq(options.freq)):
481 self._set_status_msg("Failed to set initial frequency")
483 self.myform['decim'].set_value(self.u.decim_rate())
484 self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
485 self.myform['dbname'].set_value(self.subdev.name())
486 self.myform['DM'].set_value(self.dm)
487 self.myform['Doppler'].set_value(self.doppler)
490 # Start the timer that shows current LMST on the GUI
492 self.lmst_timer.Start(1000)
495 def _set_status_msg(self, msg):
496 self.frame.GetStatusBar().SetStatusText(msg, 0)
498 def _build_gui(self, vbox):
500 def _form_set_freq(kv):
501 return self.set_freq(kv['freq'])
503 def _form_set_dm(kv):
504 return self.set_dm(kv['DM'])
506 def _form_set_doppler(kv):
507 return self.set_doppler(kv['Doppler'])
509 # Position the FFT or Waterfall
510 vbox.Add(self.scope.win, 5, wx.EXPAND)
511 vbox.Add(self.chart.win, 5, wx.EXPAND)
513 # add control area at the bottom
514 self.myform = myform = form.form()
515 hbox = wx.BoxSizer(wx.HORIZONTAL)
516 hbox.Add((7,0), 0, wx.EXPAND)
517 vbox1 = wx.BoxSizer(wx.VERTICAL)
518 myform['freq'] = form.float_field(
519 parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
520 callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
522 vbox1.Add((3,0), 0, 0)
524 # To show current Local Mean Sidereal Time
525 myform['lmst_high'] = form.static_text_field(
526 parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
527 vbox1.Add((3,0), 0, 0)
529 # To show current spectral cursor data
530 myform['spec_data'] = form.static_text_field(
531 parent=self.panel, sizer=vbox1, label="Pulse Freq", weight=1)
532 vbox1.Add((3,0), 0, 0)
534 # To show best pulses found in FFT output
535 myform['best_pulse'] = form.static_text_field(
536 parent=self.panel, sizer=vbox1, label="Best freq", weight=1)
537 vbox1.Add((3,0), 0, 0)
539 vboxBogus = wx.BoxSizer(wx.VERTICAL)
540 vboxBogus.Add ((2,0), 0, wx.EXPAND)
541 vbox2 = wx.BoxSizer(wx.VERTICAL)
542 g = self.subdev.gain_range()
543 myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
545 min=int(g[0]), max=int(g[1]),
546 callback=self.set_gain)
548 vbox2.Add((6,0), 0, 0)
549 myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
550 label="Spectral Averaging", weight=1, min=1, max=200, callback=self.set_averaging)
551 vbox2.Add((6,0), 0, 0)
552 myform['foldavg'] = form.slider_field(parent=self.panel, sizer=vbox2,
553 label="Folder Averaging", weight=1, min=1, max=20, callback=self.set_folder_averaging)
554 vbox2.Add((6,0), 0, 0)
555 myform['volume'] = form.quantized_slider_field(parent=self.panel, sizer=vbox2,
556 label="Audio Volume", weight=1, range=(-20, 0, 0.5), callback=self.set_volume)
557 vbox2.Add((6,0), 0, 0)
558 myform['DM'] = form.float_field(
559 parent=self.panel, sizer=vbox2, label="DM", weight=1,
560 callback=myform.check_input_and_call(_form_set_dm))
561 vbox2.Add((6,0), 0, 0)
562 myform['Doppler'] = form.float_field(
563 parent=self.panel, sizer=vbox2, label="Doppler", weight=1,
564 callback=myform.check_input_and_call(_form_set_doppler))
565 vbox2.Add((6,0), 0, 0)
568 # Baseband recording control
569 buttonbox = wx.BoxSizer(wx.HORIZONTAL)
570 self.record_control = form.button_with_callback(self.panel,
571 label="Recording baseband: Off ",
572 callback=self.toggle_recording)
573 self.record_pulse_control = form.button_with_callback(self.panel,
574 label="Recording pulses: Off ",
575 callback=self.toggle_pulse_recording)
577 buttonbox.Add(self.record_control, 0, wx.CENTER)
578 buttonbox.Add(self.record_pulse_control, 0, wx.CENTER)
579 vbox.Add(buttonbox, 0, wx.CENTER)
580 hbox.Add(vbox1, 0, 0)
581 hbox.Add(vboxBogus, 0, 0)
582 hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
583 vbox.Add(hbox, 0, wx.EXPAND)
585 self._build_subpanel(vbox)
587 self.lmst_timer = wx.PyTimer(self.lmst_timeout)
591 def _build_subpanel(self, vbox_arg):
592 # build a secondary information panel (sometimes hidden)
594 # FIXME figure out how to have this be a subpanel that is always
595 # created, but has its visibility controlled by foo.Show(True/False)
597 if not(self.show_debug_info):
604 #panel = wx.Panel(self.panel, -1)
605 #vbox = wx.BoxSizer(wx.VERTICAL)
607 hbox = wx.BoxSizer(wx.HORIZONTAL)
609 myform['decim'] = form.static_float_field(
610 parent=panel, sizer=hbox, label="Decim")
613 myform['fs@usb'] = form.static_float_field(
614 parent=panel, sizer=hbox, label="Fs@USB")
617 myform['dbname'] = form.static_text_field(
618 parent=panel, sizer=hbox)
621 myform['baseband'] = form.static_float_field(
622 parent=panel, sizer=hbox, label="Analog BB")
625 myform['ddc'] = form.static_float_field(
626 parent=panel, sizer=hbox, label="DDC")
629 vbox.Add(hbox, 0, wx.EXPAND)
633 def set_freq(self, target_freq):
635 Set the center frequency we're interested in.
637 @param target_freq: frequency in Hz
640 Tuning is a two step process. First we ask the front-end to
641 tune as close to the desired frequency as it can. Then we use
642 the result of that operation and our target_frequency to
643 determine the value for the digital down converter.
645 r = usrp.tune(self.u, 0, self.subdev, target_freq)
648 self.myform['freq'].set_value(target_freq) # update displayed value
649 self.myform['baseband'].set_value(r.baseband_freq)
650 self.myform['ddc'].set_value(r.dxc_freq)
651 # Adjust self.frequency, and self.observing_freq
652 # We pick up the difference between the current self.frequency
653 # and the just-programmed one, and use this to adjust
654 # self.observing_freq. We have to do it this way to
655 # make the dedispersion filtering work out properly.
656 delta = target_freq - self.frequency
657 self.frequency = target_freq
658 self.observing_freq += delta
660 # Now that we're adjusted, compute a new dispfilter, and
661 # set the taps for the FFT filter.
662 ntaps = self.compute_disp_ntaps(self.dm, self.bw, self.observing_freq)
663 self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
664 self.compute_dispfilter(self.dm,self.doppler,self.bw,
666 self.dispfilt.set_taps(self.disp_taps)
672 # Callback for gain-setting slider
673 def set_gain(self, gain):
674 self.myform['gain'].set_value(gain) # update displayed value
675 self.subdev.set_gain(gain)
678 def set_volume(self, vol):
679 self.myform['volume'].set_value(vol)
680 self.volume.set_k((10**(vol/10))/8192)
682 # Callback for spectral-averaging slider
683 def set_averaging(self, avval):
684 self.myform['average'].set_value(avval)
685 self.scope.set_avg_alpha(1.0/(avval))
686 self.scope.set_average(True)
688 def set_folder_averaging(self, avval):
689 self.myform['foldavg'].set_value(avval)
690 self.folder_iir.set_taps(1.0/avval)
692 # Timer callback to update LMST display
693 def lmst_timeout(self):
694 self.locality.date = ephem.now()
695 sidtime = self.locality.sidereal_time()
696 self.myform['lmst_high'].set_value(str(ephem.hours(sidtime)))
699 # Turn recording on/off
700 # Called-back by "Recording" button
702 def toggle_recording(self):
703 # Pick up current LMST
704 self.locality.date = ephem.now()
705 sidtime = self.locality.sidereal_time()
707 # Pick up localtime, for generating filenames
708 foo = time.localtime()
710 # Generate filenames for both data and header file
711 filename = "%04d%02d%02d%02d%02d.pdat" % (foo.tm_year, foo.tm_mon,
712 foo.tm_mday, foo.tm_hour, foo.tm_min)
713 hdrfilename = "%04d%02d%02d%02d%02d.phdr" % (foo.tm_year, foo.tm_mon,
714 foo.tm_mday, foo.tm_hour, foo.tm_min)
716 # Current recording? Flip state
717 if (self.recording_state == True):
718 self.recording_state = False
719 self.record_control.SetLabel("Recording baseband: Off ")
720 self.recording.close()
723 self.recording_state = True
724 self.record_control.SetLabel("Recording baseband to: "+filename)
726 # Cause gr_file_sink object to accept new filename
727 # note use of self.prefix--filename prefix from
728 # command line (defaults to ./)
730 self.recording.open (self.prefix+filename)
733 # We open the header file as a regular file, write header data,
735 hdrf = open(self.prefix+hdrfilename, "w")
736 hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
737 hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
738 hdrf.write("DM: "+str(self.dm)+"\n")
739 hdrf.write("doppler: "+str(self.doppler)+"\n")
741 hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
742 hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
743 hdrf.write("sample type: complex_char\n")
744 hdrf.write("sample size: "+str(gr.sizeof_char*2)+"\n")
747 # Turn recording on/off
748 # Called-back by "Recording" button
750 def toggle_pulse_recording(self):
751 # Pick up current LMST
752 self.locality.date = ephem.now()
753 sidtime = self.locality.sidereal_time()
755 # Pick up localtime, for generating filenames
756 foo = time.localtime()
758 # Generate filenames for both data and header file
759 filename = "%04d%02d%02d%02d%02d.padat" % (foo.tm_year, foo.tm_mon,
760 foo.tm_mday, foo.tm_hour, foo.tm_min)
761 hdrfilename = "%04d%02d%02d%02d%02d.pahdr" % (foo.tm_year, foo.tm_mon,
762 foo.tm_mday, foo.tm_hour, foo.tm_min)
764 # Current recording? Flip state
765 if (self.pulse_recording_state == True):
766 self.pulse_recording_state = False
767 self.record_pulse_control.SetLabel("Recording pulses: Off ")
768 self.pulse_recording.close()
771 self.pulse_recording_state = True
772 self.record_pulse_control.SetLabel("Recording pulses to: "+filename)
774 # Cause gr_file_sink object to accept new filename
775 # note use of self.prefix--filename prefix from
776 # command line (defaults to ./)
778 self.pulse_recording.open (self.prefix+filename)
781 # We open the header file as a regular file, write header data,
783 hdrf = open(self.prefix+hdrfilename, "w")
784 hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
785 hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
786 hdrf.write("DM: "+str(self.dm)+"\n")
787 hdrf.write("doppler: "+str(self.doppler)+"\n")
788 hdrf.write("pulse rate: "+str(self.pulse_freq)+"\n")
789 hdrf.write("pulse sps: "+str(self.pulse_freq*self.folding)+"\n")
790 hdrf.write("file sps: "+str(self.folder_input_rate)+"\n")
792 hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
793 hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
794 hdrf.write("sample type: short\n")
795 hdrf.write("sample size: 1\n")
798 # We get called at startup, and whenever the GUI "Set Folding params"
801 def set_folding_params(self):
802 if (self.pulse_freq <= 0):
805 # Compute required sample rate
806 self.sample_rate = int(self.pulse_freq*self.folding)
808 # And the implied decimation rate
809 required_decimation = int(self.folder_input_rate / self.sample_rate)
811 # We also compute a new FFT comb filter, based on the expected
812 # spectral profile of our pulse parameters
814 # FFT-based comb filter
816 N_COMB_TAPS=int(self.sample_rate*4)
817 if N_COMB_TAPS > 2000:
819 self.folder_comb_taps = Numeric.zeros(N_COMB_TAPS,Numeric.Complex64)
820 fincr = (self.sample_rate)/float(N_COMB_TAPS)
821 for i in range(0,len(self.folder_comb_taps)):
822 self.folder_comb_taps[i] = complex(0.0, 0.0)
825 harmonics = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
826 for i in range(0,len(self.folder_comb_taps)/2):
827 for j in range(0,len(harmonics)):
828 if abs(freq - harmonics[j]*self.pulse_freq) <= fincr:
829 self.folder_comb_taps[i] = complex(4.0, 0.0)
830 if harmonics[j] == 1.0:
831 self.folder_comb_taps[i] = complex(8.0, 0.0)
834 if self.enable_comb_filter == True:
835 # Set the just-computed FFT comb filter taps
836 self.folder_comb.set_taps(self.folder_comb_taps)
838 # And compute a new decimated bandpass filter, to go in front
839 # of the comb. Primary function is to decimate and filter down
840 # to an exact-ish multiple of the target pulse rate
842 self.folding_taps = gr.firdes_band_pass (1.0, self.folder_input_rate,
843 0.10, self.sample_rate/2, 10,
844 gr.firdes.WIN_HAMMING)
846 # Set the computed taps for the bandpass/decimate filter
847 self.folder_bandpass.set_taps (self.folding_taps)
849 # Record a spectral "hit" of a possible pulsar spectral profile
851 def record_hit(self,hits, hcavg, hcmax):
852 # Pick up current LMST
853 self.locality.date = ephem.now()
854 sidtime = self.locality.sidereal_time()
856 # Pick up localtime, for generating filenames
857 foo = time.localtime()
859 # Generate filenames for both data and header file
860 hitfilename = "%04d%02d%02d%02d.phit" % (foo.tm_year, foo.tm_mon,
861 foo.tm_mday, foo.tm_hour)
863 hitf = open(self.prefix+hitfilename, "a")
864 hitf.write("receiver center frequency: "+str(self.frequency)+"\n")
865 hitf.write("observing frequency: "+str(self.observing_freq)+"\n")
866 hitf.write("DM: "+str(self.dm)+"\n")
867 hitf.write("doppler: "+str(self.doppler)+"\n")
869 hitf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
870 hitf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
871 hitf.write("spectral peaks: "+str(hits)+"\n")
872 hitf.write("HCM: "+str(hcavg)+" "+str(hcmax)+"\n")
875 # This is a callback used by ra_fftsink.py (passed on creation of
877 # Whenever the user moves the cursor within the FFT display, this
878 # shows the coordinate data
880 def xydfunc(self,xyv):
881 s = "%.6fHz\n%.3fdB" % (xyv[0], xyv[1])
882 if self.lowpass >= 500:
883 s = "%.6fHz\n%.3fdB" % (xyv[0]*1000, xyv[1])
885 self.myform['spec_data'].set_value(s)
887 # This is another callback used by ra_fftsink.py (passed on creation
888 # of ra_fftsink). We pass this as our "calibrator" function, but
889 # we create interesting side-effects in the GUI.
891 # This function finds peaks in the FFT output data, and reports
892 # on them through the "Best" text object in the GUI
893 # It also computes the Harmonic Compliance Measure (HCM), and displays
896 def pulsarfunc(self,d,l):
898 incr = float(self.lowpass)/float(l)
905 # First, we need to find the average signal level
908 if (i * incr) > self.lowest_freq and (i*incr) < (self.lowpass-2):
911 # Set average signal level
916 # Then we find candidates that are greater than the user-supplied
919 # We try to cluster "hits" whose whole-number frequency is the
920 # same, and compute an average "hit" frequency.
928 # If frequency within bounds, and the (dB-avg) value is above our
930 if freq > self.lowest_freq and freq < self.lowpass-2 and (d[i] - avg) > self.threshold:
931 # If we're finding a new whole-number frequency
932 if lastint != int(freq):
933 # Record "center" of this hit, if this is a new hit
935 s2 += "%5.3fHz " % (freqavg/intcnt)
936 hits.append(freqavg/intcnt)
948 s2 += "%5.3fHz " % (freqavg/intcnt)
949 hits.append(freqavg/intcnt)
952 # Compute the HCM, by dividing each of the "hits" by each of the
953 # other hits, and comparing the difference between a "perfect"
954 # harmonic, and the observed frequency ratio.
961 for i in range(1,len(hits)):
962 meas = hits[i]/hits[0] - int(hits[i]/hits[0])
963 if abs((hits[i]-hits[i-1])-hits[0]) < 0.1:
964 avg_dist += hits[i]-hits[i-1]
966 if meas > 0.98 and meas < 1.0:
969 if meas >= max_measure:
979 s3="\nHCM: Avg %5.3fHz(%d) Max %5.3fHz Dist %5.3fHz(%d)" % (measure,mcnt,max_measure, avg_dist, acnt)
980 if max_measure < 0.5 and len(hits) >= 2:
981 self.record_hit(hits, measure, max_measure)
982 self.avg_dist = avg_dist
985 s4="\nAvg dB: %4.2f" % avg
986 self.myform['best_pulse'].set_value("("+s2+")"+s3+s4)
988 # Since we are nominally a calibrator function for ra_fftsink, we
989 # simply return what they sent us, untouched. A "real" calibrator
990 # function could monkey with the data before returning it to the
991 # FFT display function.
995 # Callback for the "DM" gui object
997 # We call compute_dispfilter() as appropriate to compute a new filter,
998 # and then set that new filter into self.dispfilt.
1000 def set_dm(self,dm):
1003 ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
1004 self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
1005 self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
1006 self.dispfilt.set_taps(self.disp_taps)
1007 self.myform['DM'].set_value(dm)
1011 # Callback for the "Doppler" gui object
1013 # We call compute_dispfilter() as appropriate to compute a new filter,
1014 # and then set that new filter into self.dispfilt.
1016 def set_doppler(self,doppler):
1017 self.doppler = doppler
1019 ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
1020 self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
1021 self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
1022 self.dispfilt.set_taps(self.disp_taps)
1023 self.myform['Doppler'].set_value(doppler)
1027 # Compute a de-dispersion filter
1028 # From Hankins, et al, 1975
1030 # This code translated from dedisp_filter.c from Swinburne
1031 # pulsar software repository
1033 def compute_dispfilter(self,dm,doppler,bw,centerfreq):
1034 npts = len(self.disp_taps)
1035 tmp = Numeric.zeros(npts, Numeric.Complex64)
1036 M_PI = 3.14159265358
1040 # Because astronomers are a crazy bunch, the "standard" calcultion
1041 # is in Mhz, rather than Hz
1043 centerfreq = centerfreq / 1.0e6
1046 isign = int(bw / abs (bw))
1048 # Center frequency may be doppler shifted
1049 cfreq = centerfreq / doppler
1051 # As well as the bandwidth..
1052 bandwidth = bw / doppler
1054 # Bandwidth divided among bins
1055 binwidth = bandwidth / npts
1057 # Delay is an "extra" parameter, in usecs, and largely
1058 # untested in the Swinburne code.
1061 # This determines the coefficient of the frequency response curve
1062 # Linear in DM, but quadratic in center frequency
1063 coeff = isign * 2.0*M_PI * DM / (cfreq*cfreq)
1067 for i in range(0,int(npts/2)):
1068 freq = (n + 0.5) * binwidth
1069 phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1070 tmp[i] = complex(math.cos(phi), math.sin(phi))
1076 for i in range(int(npts/2),npts):
1077 freq = (n + 0.5) * binwidth
1078 phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
1079 tmp[i] = complex(math.cos(phi), math.sin(phi))
1082 self.disp_taps = FFT.inverse_fft(tmp)
1083 return(self.disp_taps)
1086 # Compute minimum number of taps required in de-dispersion FFT filter
1088 def compute_disp_ntaps(self,dm,bw,freq):
1090 # Dt calculations are in Mhz, rather than Hz
1091 # crazy astronomers....
1095 f_lower = mfreq-(mbw/2)
1096 f_upper = mfreq+(mbw/2)
1098 # Compute smear time
1099 Dt = dm/2.41e-4 * (1.0/(f_lower*f_lower)-1.0/(f_upper*f_upper))
1101 # ntaps is now bandwidth*smeartime
1102 # Should be bandwidth*smeartime*2, but the Gnu Radio FFT filter
1103 # already expands it by a factor of 2
1110 app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY PULSAR RECEIVER: $Revision: 3626 $", nstatus=1)
1113 if __name__ == '__main__':