Houston, we have a trunk.

git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@3122 221aa14e-8319-0410-a670-987f0aec2ac5

1 files changed, 1101 insertions, 0 deletions

diff --git a/gr-radio-astronomy/src/python/usrp_psr_receiver.py b/gr-radio-astronomy/src/python/usrp_psr_receiver.py
new file mode 100755
index 000000000..8a5593706
--- /dev/null
+++ b/gr-radio-astronomy/src/python/usrp_psr_receiver.py
@@ -0,0 +1,1101 @@
+#!/usr/bin/env python
+#
+# Copyright 2004,2005 Free Software Foundation, Inc.
+# 
+# This file is part of GNU Radio
+# 
+# GNU Radio is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2, or (at your option)
+# any later version.
+# 
+# GNU Radio is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# 
+# You should have received a copy of the GNU General Public License
+# along with GNU Radio; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+# Boston, MA 02111-1307, USA.
+# 
+
+
+#
+#
+# Pulsar receiver application
+#
+# Performs both harmonic folding analysis
+#  and epoch folding analysis
+#
+#
+from gnuradio import gr, gru, blks, audio
+from gnuradio import usrp, optfir
+from gnuradio import eng_notation
+from gnuradio.eng_option import eng_option
+from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, form, slider
+from optparse import OptionParser
+import wx
+import sys
+import Numeric
+import FFT
+import ephem
+import time
+import os
+import math
+
+
+class app_flow_graph(stdgui.gui_flow_graph):
+    def __init__(self, frame, panel, vbox, argv):
+        stdgui.gui_flow_graph.__init__(self)
+
+        self.frame = frame
+        self.panel = panel
+        
+        parser = OptionParser(option_class=eng_option)
+        parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=(0, 0),
+                          help="select USRP Rx side A or B (default=A)")
+        parser.add_option("-d", "--decim", type="int", default=16,
+                          help="set fgpa decimation rate to DECIM [default=%default]")
+        parser.add_option("-f", "--freq", type="eng_float", default=None,
+                          help="set frequency to FREQ", metavar="FREQ")
+        parser.add_option("-Q", "--observing", type="eng_float", default=0.0,
+                          help="set observing frequency to FREQ")
+        parser.add_option("-a", "--avg", type="eng_float", default=1.0,
+		help="set spectral averaging alpha")
+        parser.add_option("-V", "--favg", type="eng_float", default=2.0,
+                help="set folder averaging alpha")
+        parser.add_option("-g", "--gain", type="eng_float", default=None,
+                          help="set gain in dB (default is midpoint)")
+        parser.add_option("-l", "--reflevel", type="eng_float", default=30.0,
+                          help="Set pulse display reference level")
+        parser.add_option("-L", "--lowest", type="eng_float", default=1.5,
+                          help="Lowest valid frequency bin")
+        parser.add_option("-e", "--longitude", type="eng_float", default=-76.02,                          help="Set Observer Longitude")
+        parser.add_option("-c", "--latitude", type="eng_float", default=44.85,                          help="Set Observer Latitude")
+        parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT")
+
+        parser.add_option ("-t", "--threshold", type="eng_float", default=2.5, help="pulsar threshold")
+        parser.add_option("-p", "--lowpass", type="eng_float", default=100, help="Pulse spectra cutoff freq")
+        parser.add_option("-P", "--prefix", default="./", help="File prefix")
+        parser.add_option("-u", "--pulsefreq", type="eng_float", default=0.748, help="Observation pulse rate")
+        parser.add_option("-D", "--dm", type="eng_float", default=1.0e-5, help="Dispersion Measure")
+        parser.add_option("-O", "--doppler", type="eng_float", default=1.0, help="Doppler ratio")
+        parser.add_option("-B", "--divbase", type="eng_float", default=20, help="Y/Div menu base")
+        parser.add_option("-I", "--division", type="eng_float", default=100, help="Y/Div")
+        (options, args) = parser.parse_args()
+        if len(args) != 0:
+            parser.print_help()
+            sys.exit(1)
+
+        self.show_debug_info = True
+
+        self.reflevel = options.reflevel
+        self.divbase = options.divbase
+        self.division = options.division
+
+        # Low-pass cutoff for post-detector filter
+        # Set to 100Hz usually, since lots of pulsars fit in this
+        #   range
+        self.lowpass = options.lowpass
+
+        # What is lowest valid frequency bin in post-detector FFT?
+        # There's some pollution very close to DC
+        self.lowest_freq = options.lowest
+
+        # What (dB) threshold to use in determining spectral candidates
+        self.threshold = options.threshold
+
+        # Filename prefix for recording file
+        self.prefix = options.prefix
+
+        # Dispersion Measure (DM)
+        self.dm = options.dm
+
+        # Doppler shift, as a ratio
+        #  1.0 == no doppler shift
+        #  1.005 == a little negative shift
+        #  0.995 == a little positive shift
+        self.doppler = options.doppler
+
+        #
+        # Input frequency and observing frequency--not necessarily the
+        #   same thing, if we're looking at the IF of some downconverter
+        #   that's ahead of the USRP and daughtercard.  This distinction
+        #   is important in computing the correct de-dispersion filter.
+        #
+        self.frequency = options.freq
+        if options.observing <= 0:
+            self.observing_freq = options.freq
+        else:
+            self.observing_freq = options.observing
+        
+        # build the graph
+        self.u = usrp.source_c(decim_rate=options.decim)
+        self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
+
+        #
+        # Recording file, in case we ever need to record baseband data
+        #
+        self.recording = gr.file_sink(gr.sizeof_char, "/dev/null")
+        self.recording_state = False
+
+        self.pulse_recording = gr.file_sink(gr.sizeof_short, "/dev/null")
+        self.pulse_recording_state = False
+
+        #
+        # We come up with recording turned off, but the user may
+        #  request recording later on
+        self.recording.close()
+        self.pulse_recording.close()
+
+        #
+        # Need these two for converting 12-bit baseband signals to 8-bit
+        #
+        self.tofloat = gr.complex_to_float()
+        self.tochar = gr.float_to_char()
+
+        # Need this for recording pulses (post-detector)
+        self.toshort = gr.float_to_short()
+
+
+        #
+        # The spectral measurer sets this when it has a valid
+        #   average spectral peak-to-peak distance
+        # We can then use this to program the parameters for the epoch folder
+        #
+        # We set a sentimental value here
+        self.pulse_freq = options.pulsefreq
+
+        # Folder runs at this raw sample rate
+        self.folder_input_rate = 20000
+
+        # Each pulse in the epoch folder is sampled at 128 times the nominal
+        #  pulse rate
+        self.folding = 128
+
+ 
+        #
+        # Try to find candidate parameters for rational resampler
+        #
+        save_i = 0
+        candidates = []
+        for i in range(20,300):
+            input_rate = self.folder_input_rate
+            output_rate = int(self.pulse_freq * i)
+            interp = gru.lcm(input_rate, output_rate) / input_rate
+            decim = gru.lcm(input_rate, output_rate) / output_rate
+            if (interp < 500 and decim < 250000):
+                 candidates.append(i)
+
+        # We didn't find anything, bail!
+        if (len(candidates) < 1):
+            print "Couldn't converge on resampler parameters"
+            sys.exit(1)
+
+        #
+        # Now try to find candidate with the least sampling error
+        #
+        mindiff = 999.999
+        for i in candidates:
+            diff = self.pulse_freq * i
+            diff = diff - int(diff)
+            if (diff < mindiff):
+                mindiff = diff
+                save_i = i
+
+        # Recompute rates
+        input_rate = self.folder_input_rate
+        output_rate = int(self.pulse_freq * save_i)
+
+        # Compute new interp and decim, based on best candidate
+        interp = gru.lcm(input_rate, output_rate) / input_rate
+        decim = gru.lcm(input_rate, output_rate) / output_rate
+
+        # Save optimized folding parameters, used later
+        self.folding = save_i
+        self.interp = int(interp)
+        self.decim = int(decim)
+
+        # So that we can view 4 pulses in the pulse viewer window
+        FOLD_MULT=4
+
+        # determine the daughterboard subdevice we're using
+        self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
+
+        # Compute raw input rate
+        input_rate = self.u.adc_freq() / self.u.decim_rate()
+
+        # BW==input_rate for complex data
+        self.bw = input_rate
+
+        #
+        # We use this as a crude volume control for the audio output
+        #
+        self.volume = gr.multiply_const_ff(10**(-1))
+        
+
+        #
+        # Create location data for ephem package
+        #
+        self.locality = ephem.Observer()
+        self.locality.long = str(options.longitude)
+        self.locality.lat = str(options.latitude)
+
+        #
+        # What is the post-detector LPF cutoff for the FFT?
+        #
+        PULSAR_MAX_FREQ=int(options.lowpass)
+
+        # First low-pass filters down to input_rate/FIRST_FACTOR
+        #   and decimates appropriately
+        FIRST_FACTOR=int(input_rate/(self.folder_input_rate/2))
+        first_filter = gr.firdes.low_pass (1.0,
+                                          input_rate,
+                                          input_rate/FIRST_FACTOR,
+                                          input_rate/(FIRST_FACTOR*20),         
+                                          gr.firdes.WIN_HAMMING)
+
+        # Second filter runs at the output rate of the first filter,
+        #  And low-pass filters down to PULSAR_MAX_FREQ*10
+        #
+        second_input_rate =  int(input_rate/(FIRST_FACTOR/2))
+        second_filter = gr.firdes.band_pass(1.0, second_input_rate,
+                                          0.10,
+                                          PULSAR_MAX_FREQ*10,
+                                          PULSAR_MAX_FREQ*1.5,
+                                          gr.firdes.WIN_HAMMING)
+
+        # Third filter runs at PULSAR_MAX_FREQ*20
+        #   and filters down to PULSAR_MAX_FREQ
+        #
+        third_input_rate = PULSAR_MAX_FREQ*20
+        third_filter = gr.firdes_band_pass(1.0, third_input_rate,
+                                           0.10, PULSAR_MAX_FREQ,
+                                           PULSAR_MAX_FREQ/10.0,
+                                           gr.firdes.WIN_HAMMING)
+
+
+        #
+        # Create the appropriate FFT scope
+        #
+        self.scope = ra_fftsink.ra_fft_sink_f (self, panel, 
+           fft_size=int(options.fft_size), sample_rate=PULSAR_MAX_FREQ*2,
+           title="Post-detector spectrum",  
+           cfunc=self.pulsarfunc, xydfunc=self.xydfunc, fft_rate=200)
+
+        #
+        # Tell scope we're looking from DC to PULSAR_MAX_FREQ
+        #
+        self.scope.set_baseband_freq (0.0)
+
+
+        #
+        # Setup stripchart for showing pulse profiles
+        #
+        hz = "%5.3fHz " % self.pulse_freq
+        per = "(%5.3f sec)" % (1.0/self.pulse_freq)
+        sr = "%d sps" % (int(self.pulse_freq*self.folding))
+        self.chart = ra_stripchartsink.stripchart_sink_f (self, panel,
+               sample_rate=1,
+               stripsize=self.folding*FOLD_MULT, parallel=True, title="Pulse Profiles: "+hz+per, 
+               xlabel="Seconds @ "+sr, ylabel="Level", autoscale=True,
+               divbase=self.divbase, scaling=1.0/(self.folding*self.pulse_freq))
+        self.chart.set_ref_level(self.reflevel)
+        self.chart.set_y_per_div(self.division)
+
+        # De-dispersion filter setup
+        #
+        # Do this here, just before creating the filter
+        #  that will use the taps.
+        #
+        ntaps = self.compute_disp_ntaps(self.dm,self.bw,self.observing_freq)
+
+        # Taps for the de-dispersion filter
+        self.disp_taps = Numeric.zeros(ntaps,Numeric.Complex64)
+
+        # Compute the de-dispersion filter now
+        self.compute_dispfilter(self.dm,self.doppler,
+            self.bw,self.observing_freq)
+
+        #
+        # Call constructors for receive chains
+        #
+
+        #
+        # Now create the FFT filter using the computed taps
+        self.dispfilt = gr.fft_filter_ccc(1, self.disp_taps)
+
+        #
+        # Audio sink
+        #
+        self.audio = audio.sink(second_input_rate, "plughw:0,0")
+
+        #
+        # The three post-detector filters
+        # Done this way to allow an audio path (up to 10Khz)
+        # ...and also because going from xMhz down to ~100Hz
+        # In a single filter doesn't seem to work.
+        #
+        self.first = gr.fir_filter_fff (FIRST_FACTOR/2, first_filter)
+
+        p = second_input_rate / (PULSAR_MAX_FREQ*20)
+        self.second = gr.fir_filter_fff (int(p), second_filter)
+        self.third = gr.fir_filter_fff (10, third_filter)
+
+        # Split complex USRP stream into a pair of floats
+        self.splitter = gr.complex_to_float (1);
+
+        # I squarer (detector)
+        self.multI = gr.multiply_ff();
+
+        # Q squarer (detector)
+        self.multQ = gr.multiply_ff();
+
+        # Adding squared I and Q to produce instantaneous signal power
+        self.adder = gr.add_ff();
+
+        self.enable_comb_filter = False
+        # Epoch folder comb filter
+        if self.enable_comb_filter == True:
+            bogtaps = Numeric.zeros(512, Numeric.Float64)
+            self.folder_comb = gr.fft_filter_ccc(1,bogtaps)
+
+        # Rational resampler
+        self.folder_rr = blks.rational_resampler_fff(self, self.interp, self.decim)
+
+        # Epoch folder bandpass
+        bogtaps = Numeric.zeros(1, Numeric.Float64)
+        self.folder_bandpass = gr.fir_filter_fff (1, bogtaps)
+
+        # Epoch folder F2C/C2F
+        self.folder_f2c = gr.float_to_complex()
+        self.folder_c2f = gr.complex_to_float()
+
+        # Epoch folder S2P
+        self.folder_s2p = gr.serial_to_parallel (gr.sizeof_float, 
+             self.folding*FOLD_MULT)
+
+        # Epoch folder IIR Filter (produces average pulse profiles)
+        self.folder_iir = gr.single_pole_iir_filter_ff(1.0/options.favg,
+             self.folding*FOLD_MULT)
+
+        #
+        # Set all the epoch-folder goop up
+        #
+        self.set_folding_params()
+
+        # 
+        # Start connecting configured modules in the receive chain
+        #
+
+        # Connect raw USRP to de-dispersion filter, complex->float splitter
+        self.connect(self.u, self.dispfilt, self.splitter)
+
+        # Connect splitter outputs to multipliers
+        # First do I^2
+        self.connect((self.splitter, 0), (self.multI,0))
+        self.connect((self.splitter, 0), (self.multI,1))
+
+        # Then do Q^2
+        self.connect((self.splitter, 1), (self.multQ,0))
+        self.connect((self.splitter, 1), (self.multQ,1))
+
+        # Then sum the squares
+        self.connect(self.multI, (self.adder,0))
+        self.connect(self.multQ, (self.adder,1))
+
+        # Connect detector/adder output to FIR LPF
+        #  in two stages, followed by the FFT scope
+        self.connect(self.adder, self.first,
+            self.second, self.third, self.scope)
+
+        # Connect audio output
+        self.connect(self.first, self.volume)
+        self.connect(self.volume, (self.audio, 0))
+        self.connect(self.volume, (self.audio, 1))
+
+        # Connect epoch folder
+        if self.enable_comb_filter == True:
+            self.connect (self.first, self.folder_bandpass, self.folder_rr,
+                self.folder_f2c,
+                self.folder_comb, self.folder_c2f,
+                self.folder_s2p, self.folder_iir,
+                self.chart)
+
+        else:
+            self.connect (self.first, self.folder_bandpass, self.folder_rr,
+                self.folder_s2p, self.folder_iir, self.chart)
+
+        # Connect baseband recording file (initially /dev/null)
+        self.connect(self.u, self.tofloat, self.tochar, self.recording)
+
+        # Connect pulse recording file (initially /dev/null)
+        self.connect(self.first, self.toshort, self.pulse_recording)
+
+        #
+        # Build the GUI elements
+        #
+        self._build_gui(vbox)
+
+        # Make GUI agree with command-line
+        self.myform['average'].set_value(int(options.avg))
+        self.myform['foldavg'].set_value(int(options.favg))
+
+
+        # Make spectral averager agree with command line
+        if options.avg != 1.0:
+            self.scope.set_avg_alpha(float(1.0/options.avg))
+            self.scope.set_average(True)
+
+
+        # set initial values
+
+        if options.gain is None:
+            # if no gain was specified, use the mid-point in dB
+            g = self.subdev.gain_range()
+            options.gain = float(g[0]+g[1])/2
+
+        if options.freq is None:
+            # if no freq was specified, use the mid-point
+            r = self.subdev.freq_range()
+            options.freq = float(r[0]+r[1])/2
+
+        self.set_gain(options.gain)
+        self.set_volume(-10.0)
+
+        if not(self.set_freq(options.freq)):
+            self._set_status_msg("Failed to set initial frequency")
+
+        self.myform['decim'].set_value(self.u.decim_rate())
+        self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
+        self.myform['dbname'].set_value(self.subdev.name())
+        self.myform['DM'].set_value(self.dm)
+        self.myform['Doppler'].set_value(self.doppler)
+
+        #
+        # Start the timer that shows current LMST on the GUI
+        #
+        self.lmst_timer.Start(1000)
+
+
+    def _set_status_msg(self, msg):
+        self.frame.GetStatusBar().SetStatusText(msg, 0)
+
+    def _build_gui(self, vbox):
+
+        def _form_set_freq(kv):
+            return self.set_freq(kv['freq'])
+
+        def _form_set_dm(kv):
+            return self.set_dm(kv['DM'])
+
+        def _form_set_doppler(kv):
+            return self.set_doppler(kv['Doppler'])
+
+        # Position the FFT or Waterfall
+        vbox.Add(self.scope.win, 5, wx.EXPAND)
+        vbox.Add(self.chart.win, 5, wx.EXPAND)
+
+        # add control area at the bottom
+        self.myform = myform = form.form()
+        hbox = wx.BoxSizer(wx.HORIZONTAL)
+        hbox.Add((7,0), 0, wx.EXPAND)
+        vbox1 = wx.BoxSizer(wx.VERTICAL)
+        myform['freq'] = form.float_field(
+            parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
+            callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
+
+        vbox1.Add((3,0), 0, 0)
+
+        # To show current Local Mean Sidereal Time
+        myform['lmst_high'] = form.static_text_field(
+            parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
+        vbox1.Add((3,0), 0, 0)
+
+        # To show current spectral cursor data
+        myform['spec_data'] = form.static_text_field(
+            parent=self.panel, sizer=vbox1, label="Pulse Freq", weight=1)
+        vbox1.Add((3,0), 0, 0)
+
+        # To show best pulses found in FFT output
+        myform['best_pulse'] = form.static_text_field(
+            parent=self.panel, sizer=vbox1, label="Best freq", weight=1)
+        vbox1.Add((3,0), 0, 0)
+
+        vboxBogus = wx.BoxSizer(wx.VERTICAL)
+        vboxBogus.Add ((2,0), 0, wx.EXPAND)
+        vbox2 = wx.BoxSizer(wx.VERTICAL)
+        g = self.subdev.gain_range()
+        myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
+                                           weight=1,
+                                           min=int(g[0]), max=int(g[1]),
+                                           callback=self.set_gain)
+
+        vbox2.Add((6,0), 0, 0)
+        myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2, 
+                    label="Spectral Averaging", weight=1, min=1, max=200, callback=self.set_averaging)
+        vbox2.Add((6,0), 0, 0)
+        myform['foldavg'] = form.slider_field(parent=self.panel, sizer=vbox2,
+                    label="Folder Averaging", weight=1, min=1, max=20, callback=self.set_folder_averaging)
+        vbox2.Add((6,0), 0, 0)
+        myform['volume'] = form.quantized_slider_field(parent=self.panel, sizer=vbox2,
+                    label="Audio Volume", weight=1, range=(-20, 0, 0.5), callback=self.set_volume)
+        vbox2.Add((6,0), 0, 0)
+        myform['DM'] = form.float_field(
+            parent=self.panel, sizer=vbox2, label="DM", weight=1,
+            callback=myform.check_input_and_call(_form_set_dm))
+        vbox2.Add((6,0), 0, 0)
+        myform['Doppler'] = form.float_field(
+            parent=self.panel, sizer=vbox2, label="Doppler", weight=1,
+            callback=myform.check_input_and_call(_form_set_doppler))
+        vbox2.Add((6,0), 0, 0)
+
+
+        # Baseband recording control
+        buttonbox = wx.BoxSizer(wx.HORIZONTAL)
+        self.record_control = form.button_with_callback(self.panel,
+              label="Recording baseband: Off                           ",
+              callback=self.toggle_recording)
+        self.record_pulse_control = form.button_with_callback(self.panel,
+              label="Recording pulses: Off                              ",
+              callback=self.toggle_pulse_recording)
+
+        buttonbox.Add(self.record_control, 0, wx.CENTER)
+        buttonbox.Add(self.record_pulse_control, 0, wx.CENTER)
+        vbox.Add(buttonbox, 0, wx.CENTER)
+        hbox.Add(vbox1, 0, 0)
+        hbox.Add(vboxBogus, 0, 0)
+	hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
+        vbox.Add(hbox, 0, wx.EXPAND)
+
+        self._build_subpanel(vbox)
+
+        self.lmst_timer = wx.PyTimer(self.lmst_timeout)
+        self.lmst_timeout()
+
+
+    def _build_subpanel(self, vbox_arg):
+        # build a secondary information panel (sometimes hidden)
+
+        # FIXME figure out how to have this be a subpanel that is always
+        # created, but has its visibility controlled by foo.Show(True/False)
+        
+        if not(self.show_debug_info):
+            return
+
+        panel = self.panel
+        vbox = vbox_arg
+        myform = self.myform
+
+        #panel = wx.Panel(self.panel, -1)
+        #vbox = wx.BoxSizer(wx.VERTICAL)
+
+        hbox = wx.BoxSizer(wx.HORIZONTAL)
+        hbox.Add((5,0), 0)
+        myform['decim'] = form.static_float_field(
+            parent=panel, sizer=hbox, label="Decim")
+
+        hbox.Add((5,0), 1)
+        myform['fs@usb'] = form.static_float_field(
+            parent=panel, sizer=hbox, label="Fs@USB")
+
+        hbox.Add((5,0), 1)
+        myform['dbname'] = form.static_text_field(
+            parent=panel, sizer=hbox)
+
+        hbox.Add((5,0), 1)
+        myform['baseband'] = form.static_float_field(
+            parent=panel, sizer=hbox, label="Analog BB")
+
+        hbox.Add((5,0), 1)
+        myform['ddc'] = form.static_float_field(
+            parent=panel, sizer=hbox, label="DDC")
+
+        hbox.Add((5,0), 0)
+        vbox.Add(hbox, 0, wx.EXPAND)
+
+        
+        
+    def set_freq(self, target_freq):
+        """
+        Set the center frequency we're interested in.
+
+        @param target_freq: frequency in Hz
+        @rypte: bool
+
+        Tuning is a two step process.  First we ask the front-end to
+        tune as close to the desired frequency as it can.  Then we use
+        the result of that operation and our target_frequency to
+        determine the value for the digital down converter.
+        """
+        r = usrp.tune(self.u, 0, self.subdev, target_freq)
+
+        if r:
+            self.myform['freq'].set_value(target_freq)     # update displayed value
+            self.myform['baseband'].set_value(r.baseband_freq)
+            self.myform['ddc'].set_value(r.dxc_freq)
+            # Adjust self.frequency, and self.observing_freq
+            # We pick up the difference between the current self.frequency
+            #   and the just-programmed one, and use this to adjust
+            #   self.observing_freq.  We have to do it this way to
+            #   make the dedispersion filtering work out properly.
+            delta = target_freq - self.frequency
+            self.frequency = target_freq
+            self.observing_freq += delta
+
+            # Now that we're adjusted, compute a new dispfilter, and
+            #   set the taps for the FFT filter.
+            ntaps = self.compute_disp_ntaps(self.dm, self.bw, self.observing_freq)
+            self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
+            self.compute_dispfilter(self.dm,self.doppler,self.bw,
+                self.observing_freq)
+            self.dispfilt.set_taps(self.disp_taps)
+
+            return True
+
+        return False
+
+    # Callback for gain-setting slider
+    def set_gain(self, gain):
+        self.myform['gain'].set_value(gain)     # update displayed value
+        self.subdev.set_gain(gain)
+
+
+    def set_volume(self, vol):
+        self.myform['volume'].set_value(vol)
+        self.volume.set_k((10**(vol/10))/8192)
+
+    # Callback for spectral-averaging slider
+    def set_averaging(self, avval):
+        self.myform['average'].set_value(avval)
+        self.scope.set_avg_alpha(1.0/(avval))
+        self.scope.set_average(True)
+
+    def set_folder_averaging(self, avval):
+        self.myform['foldavg'].set_value(avval)
+        self.folder_iir.set_taps(1.0/avval)
+
+    # Timer callback to update LMST display
+    def lmst_timeout(self):
+         self.locality.date = ephem.now()
+         sidtime = self.locality.sidereal_time()
+         self.myform['lmst_high'].set_value(str(ephem.hours(sidtime)))
+
+    #
+    # Turn recording on/off
+    # Called-back by "Recording" button
+    #
+    def toggle_recording(self):
+        # Pick up current LMST
+        self.locality.date = ephem.now()
+        sidtime = self.locality.sidereal_time()
+
+        # Pick up localtime, for generating filenames
+        foo = time.localtime()
+
+        # Generate filenames for both data and header file
+        filename = "%04d%02d%02d%02d%02d.pdat" % (foo.tm_year, foo.tm_mon,
+           foo.tm_mday, foo.tm_hour, foo.tm_min)
+        hdrfilename = "%04d%02d%02d%02d%02d.phdr" % (foo.tm_year, foo.tm_mon,
+           foo.tm_mday, foo.tm_hour, foo.tm_min)
+
+        # Current recording?  Flip state
+        if (self.recording_state == True):
+          self.recording_state = False
+          self.record_control.SetLabel("Recording baseband: Off                           ")
+          self.recording.close()
+        # Not recording?
+        else:
+          self.recording_state = True
+          self.record_control.SetLabel("Recording baseband to: "+filename)
+
+          # Cause gr_file_sink object to accept new filename
+          #   note use of self.prefix--filename prefix from
+          #   command line (defaults to ./)
+          #
+          self.recording.open (self.prefix+filename)
+
+          #
+          # We open the header file as a regular file, write header data,
+          #   then close
+          hdrf = open(self.prefix+hdrfilename, "w")
+          hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
+          hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
+          hdrf.write("DM: "+str(self.dm)+"\n")
+          hdrf.write("doppler: "+str(self.doppler)+"\n")
+
+          hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
+          hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
+          hdrf.write("sample type: complex_char\n")
+          hdrf.write("sample size: "+str(gr.sizeof_char*2)+"\n")
+          hdrf.close()
+    #
+    # Turn recording on/off
+    # Called-back by "Recording" button
+    #
+    def toggle_pulse_recording(self):
+        # Pick up current LMST
+        self.locality.date = ephem.now()
+        sidtime = self.locality.sidereal_time()
+
+        # Pick up localtime, for generating filenames
+        foo = time.localtime()
+
+        # Generate filenames for both data and header file
+        filename = "%04d%02d%02d%02d%02d.padat" % (foo.tm_year, foo.tm_mon,
+           foo.tm_mday, foo.tm_hour, foo.tm_min)
+        hdrfilename = "%04d%02d%02d%02d%02d.pahdr" % (foo.tm_year, foo.tm_mon,
+           foo.tm_mday, foo.tm_hour, foo.tm_min)
+
+        # Current recording?  Flip state
+        if (self.pulse_recording_state == True):
+          self.pulse_recording_state = False
+          self.record_pulse_control.SetLabel("Recording pulses: Off                           ")
+          self.pulse_recording.close()
+        # Not recording?
+        else:
+          self.pulse_recording_state = True
+          self.record_pulse_control.SetLabel("Recording pulses to: "+filename)
+
+          # Cause gr_file_sink object to accept new filename
+          #   note use of self.prefix--filename prefix from
+          #   command line (defaults to ./)
+          #
+          self.pulse_recording.open (self.prefix+filename)
+
+          #
+          # We open the header file as a regular file, write header data,
+          #   then close
+          hdrf = open(self.prefix+hdrfilename, "w")
+          hdrf.write("receiver center frequency: "+str(self.frequency)+"\n")
+          hdrf.write("observing frequency: "+str(self.observing_freq)+"\n")
+          hdrf.write("DM: "+str(self.dm)+"\n")
+          hdrf.write("doppler: "+str(self.doppler)+"\n")
+          hdrf.write("pulse rate: "+str(self.pulse_freq)+"\n")
+          hdrf.write("pulse sps: "+str(self.pulse_freq*self.folding)+"\n")
+          hdrf.write("file sps: "+str(self.folder_input_rate)+"\n")
+
+          hdrf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
+          hdrf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
+          hdrf.write("sample type: short\n")
+          hdrf.write("sample size: 1\n")
+          hdrf.close()
+
+    # We get called at startup, and whenever the GUI "Set Folding params"
+    #   button is pressed
+    #
+    def set_folding_params(self):
+        if (self.pulse_freq <= 0):
+            return
+
+        # Compute required sample rate
+        self.sample_rate = int(self.pulse_freq*self.folding)
+
+        # And the implied decimation rate
+        required_decimation = int(self.folder_input_rate / self.sample_rate)
+
+        # We also compute a new FFT comb filter, based on the expected
+        #  spectral profile of our pulse parameters
+        #
+        # FFT-based comb filter
+        #
+        N_COMB_TAPS=int(self.sample_rate*4)
+        if N_COMB_TAPS > 2000:
+            N_COMB_TAPS = 2000
+        self.folder_comb_taps = Numeric.zeros(N_COMB_TAPS,Numeric.Complex64)
+        fincr = (self.sample_rate)/float(N_COMB_TAPS)
+        for i in range(0,len(self.folder_comb_taps)):
+            self.folder_comb_taps[i] = complex(0.0, 0.0)
+
+        freq = 0.0
+        harmonics = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
+        for i in range(0,len(self.folder_comb_taps)/2):
+            for j in range(0,len(harmonics)):
+                 if abs(freq - harmonics[j]*self.pulse_freq) <= fincr:
+                     self.folder_comb_taps[i] = complex(4.0, 0.0)
+                     if harmonics[j] == 1.0:
+                         self.folder_comb_taps[i] = complex(8.0, 0.0)
+            freq += fincr
+
+        if self.enable_comb_filter == True:
+            # Set the just-computed FFT comb filter taps
+            self.folder_comb.set_taps(self.folder_comb_taps)
+
+        # And compute a new decimated bandpass filter, to go in front
+        #   of the comb.  Primary function is to decimate and filter down
+        #   to an exact-ish multiple of the target pulse rate
+        #
+        self.folding_taps = gr.firdes_band_pass (1.0, self.folder_input_rate,
+            0.10, self.sample_rate/2, 10, 
+            gr.firdes.WIN_HAMMING)
+
+        # Set the computed taps for the bandpass/decimate filter
+        self.folder_bandpass.set_taps (self.folding_taps)
+    #
+    # Record a spectral "hit" of a possible pulsar spectral profile
+    #
+    def record_hit(self,hits, hcavg, hcmax):
+        # Pick up current LMST
+        self.locality.date = ephem.now()
+        sidtime = self.locality.sidereal_time()
+
+        # Pick up localtime, for generating filenames
+        foo = time.localtime()
+
+        # Generate filenames for both data and header file
+        hitfilename = "%04d%02d%02d%02d.phit" % (foo.tm_year, foo.tm_mon,
+           foo.tm_mday, foo.tm_hour)
+
+        hitf = open(self.prefix+hitfilename, "a")
+        hitf.write("receiver center frequency: "+str(self.frequency)+"\n")
+        hitf.write("observing frequency: "+str(self.observing_freq)+"\n")
+        hitf.write("DM: "+str(self.dm)+"\n")
+        hitf.write("doppler: "+str(self.doppler)+"\n")
+
+        hitf.write("sidereal: "+str(ephem.hours(sidtime))+"\n")
+        hitf.write("bandwidth: "+str(self.u.adc_freq() / self.u.decim_rate())+"\n")
+        hitf.write("spectral peaks: "+str(hits)+"\n")
+        hitf.write("HCM: "+str(hcavg)+" "+str(hcmax)+"\n")
+        hitf.close()
+
+    # This is a callback used by ra_fftsink.py (passed on creation of
+    #   ra_fftsink)
+    # Whenever the user moves the cursor within the FFT display, this
+    #   shows the coordinate data
+    #
+    def xydfunc(self,xyv):
+        s = "%.6fHz\n%.3fdB" % (xyv[0], xyv[1])
+        if self.lowpass >= 500:
+            s = "%.6fHz\n%.3fdB" % (xyv[0]*1000, xyv[1])
+        
+        self.myform['spec_data'].set_value(s)
+
+    # This is another callback used by ra_fftsink.py (passed on creation
+    #   of ra_fftsink).  We pass this as our "calibrator" function, but
+    #   we create interesting side-effects in the GUI.
+    #
+    # This function finds peaks in the FFT output data, and reports
+    #  on them through the "Best" text object in the GUI
+    #  It also computes the Harmonic Compliance Measure (HCM), and displays
+    #  that also.
+    #
+    def pulsarfunc(self,d,l):
+       x = range(0,l)
+       incr = float(self.lowpass)/float(l)
+       incr = incr * 2.0
+       bestdb = -50.0
+       bestfreq = 0.0
+       avg = 0
+       dcnt = 0
+       #
+       # First, we need to find the average signal level
+       #
+       for i in x:
+           if (i * incr) > self.lowest_freq and (i*incr) < (self.lowpass-2):
+               avg += d[i]
+               dcnt += 1
+       # Set average signal level
+       avg /= dcnt
+       s2=" "
+       findcnt = 0
+       #
+       # Then we find candidates that are greater than the user-supplied
+       #   threshold.
+       #
+       # We try to cluster "hits" whose whole-number frequency is the
+       #   same, and compute an average "hit" frequency.
+       #
+       lastint = 0
+       hits=[]
+       intcnt = 0
+       freqavg = 0
+       for i in x:
+           freq = i*incr
+           # If frequency within bounds, and the (dB-avg) value is above our
+           #   threshold
+           if freq > self.lowest_freq and freq < self.lowpass-2 and (d[i] - avg) > self.threshold:
+               # If we're finding a new whole-number frequency
+               if lastint != int(freq):
+                   # Record "center" of this hit, if this is a new hit
+                   if lastint != 0:
+                       s2 += "%5.3fHz " % (freqavg/intcnt)
+                       hits.append(freqavg/intcnt)
+                       findcnt += 1
+                   lastint = int(freq)
+                   intcnt = 1
+                   freqavg = freq
+               else:
+                   intcnt += 1
+                   freqavg += freq
+           if (findcnt >= 14):
+               break
+
+       if intcnt > 1:
+           s2 += "%5.3fHz " % (freqavg/intcnt)
+           hits.append(freqavg/intcnt)
+
+       #
+       # Compute the HCM, by dividing each of the "hits" by each of the
+       #   other hits, and comparing the difference between a "perfect"
+       #   harmonic, and the observed frequency ratio.
+       #
+       measure = 0
+       max_measure=0
+       mcnt = 0
+       avg_dist = 0
+       acnt = 0
+       for i in range(1,len(hits)):
+           meas = hits[i]/hits[0] - int(hits[i]/hits[0])
+           if abs((hits[i]-hits[i-1])-hits[0]) < 0.1:
+               avg_dist += hits[i]-hits[i-1]
+               acnt += 1
+           if meas > 0.98 and meas < 1.0:
+               meas = 1.0 - meas
+           meas *= hits[0]
+           if meas >= max_measure:
+               max_measure = meas
+           measure += meas
+           mcnt += 1
+       if mcnt > 0:
+           measure /= mcnt
+           if acnt > 0:
+               avg_dist /= acnt
+       if len(hits) > 1:
+           measure /= mcnt
+           s3="\nHCM: Avg %5.3fHz(%d) Max %5.3fHz Dist %5.3fHz(%d)" % (measure,mcnt,max_measure, avg_dist, acnt)
+           if max_measure < 0.5 and len(hits) >= 2:
+               self.record_hit(hits, measure, max_measure)
+               self.avg_dist = avg_dist
+       else:
+           s3="\nHCM: --"
+       s4="\nAvg dB: %4.2f" % avg
+       self.myform['best_pulse'].set_value("("+s2+")"+s3+s4)
+
+       # Since we are nominally a calibrator function for ra_fftsink, we
+       #  simply return what they sent us, untouched.  A "real" calibrator
+       #  function could monkey with the data before returning it to the
+       #  FFT display function.
+       return(d)
+
+    #
+    # Callback for the "DM" gui object
+    #
+    # We call compute_dispfilter() as appropriate to compute a new filter,
+    #   and then set that new filter into self.dispfilt.
+    #
+    def set_dm(self,dm):
+       self.dm = dm
+
+       ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
+       self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
+       self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
+       self.dispfilt.set_taps(self.disp_taps)
+       self.myform['DM'].set_value(dm)
+       return(dm)
+
+    #
+    # Callback for the "Doppler" gui object
+    #
+    # We call compute_dispfilter() as appropriate to compute a new filter,
+    #   and then set that new filter into self.dispfilt.
+    #
+    def set_doppler(self,doppler):
+       self.doppler = doppler
+
+       ntaps = self.compute_disp_ntaps (self.dm, self.bw, self.observing_freq)
+       self.disp_taps = Numeric.zeros(ntaps, Numeric.Complex64)
+       self.compute_dispfilter(self.dm,self.doppler,self.bw,self.observing_freq)
+       self.dispfilt.set_taps(self.disp_taps)
+       self.myform['Doppler'].set_value(doppler)
+       return(doppler)
+
+    #
+    # Compute a de-dispersion filter
+    #  From Hankins, et al, 1975
+    #
+    # This code translated from dedisp_filter.c from Swinburne
+    #   pulsar software repository
+    #
+    def compute_dispfilter(self,dm,doppler,bw,centerfreq):
+        npts = len(self.disp_taps)
+        tmp = Numeric.zeros(npts, Numeric.Complex64)
+        M_PI = 3.14159265358
+        DM = dm/2.41e-10
+
+        #
+        # Because astronomers are a crazy bunch, the "standard" calcultion
+        #   is in Mhz, rather than Hz
+        #
+        centerfreq = centerfreq / 1.0e6
+        bw = bw / 1.0e6
+        
+        isign = int(bw / abs (bw))
+        
+        # Center frequency may be doppler shifted
+        cfreq     = centerfreq / doppler
+
+        # As well as the bandwidth..
+        bandwidth = bw / doppler
+
+        # Bandwidth divided among bins
+        binwidth  = bandwidth / npts
+
+        # Delay is an "extra" parameter, in usecs, and largely
+        #  untested in the Swinburne code.
+        delay = 0.0
+        
+        # This determines the coefficient of the frequency response curve
+        # Linear in DM, but quadratic in center frequency
+        coeff = isign * 2.0*M_PI * DM / (cfreq*cfreq)
+        
+        # DC to nyquist/2
+        n = 0
+        for i in range(0,int(npts/2)):
+            freq = (n + 0.5) * binwidth
+            phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
+            tmp[i] = complex(math.cos(phi), math.sin(phi))
+            n += 1
+
+        # -nyquist/2 to DC
+        n = int(npts/2)
+        n *= -1
+        for i in range(int(npts/2),npts):
+            freq = (n + 0.5) * binwidth
+            phi = coeff*freq*freq/(cfreq+freq) + (2.0*M_PI*freq*delay)
+            tmp[i] = complex(math.cos(phi), math.sin(phi))
+            n += 1
+        
+        self.disp_taps = FFT.inverse_fft(tmp)
+        return(self.disp_taps)
+
+    #
+    # Compute minimum number of taps required in de-dispersion FFT filter
+    #
+    def compute_disp_ntaps(self,dm,bw,freq):
+        #
+        # Dt calculations are in Mhz, rather than Hz
+        #    crazy astronomers....
+        mbw = bw/1.0e6
+        mfreq = freq/1.0e6
+
+        f_lower = mfreq-(mbw/2)
+        f_upper = mfreq+(mbw/2)
+
+        # Compute smear time
+        Dt = dm/2.41e-4 * (1.0/(f_lower*f_lower)-1.0/(f_upper*f_upper))
+
+        # ntaps is now bandwidth*smeartime
+        # Should be bandwidth*smeartime*2, but the Gnu Radio FFT filter
+        #   already expands it by a factor of 2
+        ntaps = bw*Dt
+        if ntaps < 64:
+            ntaps = 64
+        return(int(ntaps))
+
+def main ():
+    app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY PULSAR RECEIVER: $Revision$", nstatus=1)
+    app.MainLoop()
+
+if __name__ == '__main__':
+    main ()