diff options
Diffstat (limited to 'gr-radio-astronomy/src/python')
| -rw-r--r-- | gr-radio-astronomy/src/python/Makefile.am | 56 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/local_calibrator.py | 149 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/qa_ra.py | 38 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/ra_fftsink.py | 520 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/ra_stripchartsink.py | 399 | ||||
| -rw-r--r-- | gr-radio-astronomy/src/python/run_tests.in | 47 | ||||
| -rw-r--r-- | gr-radio-astronomy/src/python/usrp_psr_receiver.help | 111 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/usrp_psr_receiver.py | 1101 | ||||
| -rw-r--r-- | gr-radio-astronomy/src/python/usrp_ra_receiver.help | 90 | ||||
| -rwxr-xr-x | gr-radio-astronomy/src/python/usrp_ra_receiver.py | 584 |
10 files changed, 3095 insertions, 0 deletions
diff --git a/gr-radio-astronomy/src/python/Makefile.am b/gr-radio-astronomy/src/python/Makefile.am new file mode 100644 index 000000000..62706fced --- /dev/null +++ b/gr-radio-astronomy/src/python/Makefile.am @@ -0,0 +1,56 @@ +# +# Copyright 2004,2006 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. +# + +include $(top_srcdir)/Makefile.common + +# Install this stuff so that it ends up as the gnuradio.ra module +# This usually ends up at: +# ${prefix}/lib/python${python_version}/site-packages/gnuradio + +ourpythondir = $(grpythondir) +ourlibdir = $(grpyexecdir) + +# wxgui stuff here +wxguipythondir = $(grpythondir)/wxgui +wxguilibdir = $(grpyexecdir)/wxgui + +EXTRA_DIST = run_tests.in + + +TESTS = \ + run_tests + + +noinst_PYTHON = \ + qa_ra.py + +ourpython_PYTHON = \ + local_calibrator.py + +wxguipython_PYTHON = \ + ra_stripchartsink.py \ + ra_fftsink.py + + +# and here for applications you want installed in prefix/bin +bin_SCRIPTS = \ + usrp_ra_receiver.py \ + usrp_psr_receiver.py diff --git a/gr-radio-astronomy/src/python/local_calibrator.py b/gr-radio-astronomy/src/python/local_calibrator.py new file mode 100755 index 000000000..6c5fb4f1f --- /dev/null +++ b/gr-radio-astronomy/src/python/local_calibrator.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python +# +# Copyright 2003,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. +# + +import Numeric +import math +import ephem +import time + +# +# Simple class for allowing local definition of a calibration function +# for raw samples coming from the RA detector chain. Each observatory +# is different, and rather than hacking up the main code in usrp_ra_receiver +# we define the appropriate function here. +# +# For example, one could calibrate the output in Janskys, rather than +# dB. +# +# + +def calib_default_total_power(data): + r = 10.0*math.log10(data) + return(r) + +def calib_numogate_ridge_observatory_total_power(data): + + me = ephem.Observer() + + # + # PyEphem wants lat/long as strings, rather than floats--took me quite + # a long time to figure that out. If they don't arrive as strings, + # the calculations for sidereal time are complete garbage + # + me.long = str(-76.043) + me.lat = str(44.967) + + me.date = ephem.now() + sidtime = me.sidereal_time() + + foo = time.localtime() + if not "calib_prefix" in globals(): + pfx = "./" + else: + pfx = globals()["calib_prefix"] + filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year, + foo.tm_mon, foo.tm_mday, foo.tm_hour) + + rainbow_file = open (filenamestr+".tpdat","a") + + r = (math.sqrt(data) / 2048) * 1000.0 + #r = calib_default_total_power(data) + inter = globals()["calib_decln"] + rainbow_file.write(str(ephem.hours(sidtime))+" "+str(r)+" "+str(inter)+"\n") + rainbow_file.close() + return(r) + +def calib_numogate_ridge_observatory_fft(data,l): + + me = ephem.Observer() + + # + # PyEphem wants lat/long as strings, rather than floats--took me quite + # a long time to figure that out. If they don't arrive as strings, + # the calculations for sidereal time are complete garbage + # + me.long = str(-76.043) + me.lat = str(44.967) + + me.date = ephem.now() + sidtime = me.sidereal_time() + + foo = time.localtime() + + if not "calib_prefix" in globals(): + pfx = "./" + else: + pfx = globals()["calib_prefix"] + filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year, + foo.tm_mon, foo.tm_mday, foo.tm_hour) + + now = time.time() + + if not "calib_then" in globals(): + globals()["calib_then"] = now + + delta = 5 + + if (now - globals()["calib_then"]) >= delta: + + globals()["calib_then"] = now + rainbow_file = open (filenamestr+".sdat","a") + + r = calib_default_fft(data,l) + inter = globals()["calib_decln"] + rainbow_file.write("data:"+str(ephem.hours(sidtime))+" "+str(inter)+" "+str(r)+"\n") + rainbow_file.close() + return(r) + + return(data) + +def calib_default_fft(db,l): + return(db) + +# +# We capture various parameters from the receive chain here, because +# they can affect the calibration equations. +# +# +def calib_set_gain(gain): + globals()["calib_gain_setting"] = gain + +def calib_set_integ(integ): + globals()["calib_integ_setting"] = integ + +def calib_set_bw(bw): + globals()["calib_bw_setting"] = bw + +def calib_set_freq(freq): + globals()["calib_freq_setting"] = freq + +def calib_set_avg_alpha(alpha): + globals()["calib_avg_alpha"] = alpha + +def calib_set_interesting(inter): + globals()["calib_is_interesting"] = inter + +def calib_set_decln(dec): + globals()["calib_decln"] = dec + +def calib_set_prefix(pfx): + globals()["calib_prefix"] = pfx diff --git a/gr-radio-astronomy/src/python/qa_ra.py b/gr-radio-astronomy/src/python/qa_ra.py new file mode 100755 index 000000000..4b24a1eb2 --- /dev/null +++ b/gr-radio-astronomy/src/python/qa_ra.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python +# +# Copyright 2004,2006 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. +# + +from gnuradio import gr, gr_unittest +import ra + +class qa_ra (gr_unittest.TestCase): + + def setUp (self): + self.fg = gr.flow_graph () + + def tearDown (self): + self.fg = None + + def test_000_(self): # ensure that we can load the module + pass + +if __name__ == '__main__': + gr_unittest.main () diff --git a/gr-radio-astronomy/src/python/ra_fftsink.py b/gr-radio-astronomy/src/python/ra_fftsink.py new file mode 100755 index 000000000..213f4a72a --- /dev/null +++ b/gr-radio-astronomy/src/python/ra_fftsink.py @@ -0,0 +1,520 @@ +#!/usr/bin/env python +# +# Copyright 2003,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. +# + +from gnuradio import gr, gru, window +from gnuradio.wxgui import stdgui +import wx +import gnuradio.wxgui.plot as plot +import Numeric +import threading +import math +import random + +default_ra_fftsink_size = (640,140) + + +def default_cfunc(db,l): + return(db) + + +class ra_fft_sink_base(object): + def __init__(self, input_is_real=False, baseband_freq=0, y_per_div=10, sc_y_per_div=0.5, ref_level=50, sc_ref_level=20, + sample_rate=1, fft_size=512, fft_rate=15, + average=False, avg_alpha=None, title='', peak_hold=False, cfunc=default_cfunc, xydfunc=None, interfunc=None): + + # initialize common attributes + self.baseband_freq = baseband_freq + self.y_divs = 8 + self.y_per_div=y_per_div + self.sc_y_per_div=sc_y_per_div + self.ref_level = ref_level + self.autoscale = False + self.sc_ref_level = sc_ref_level + self.sample_rate = sample_rate + self.fft_size = fft_size + self.fft_rate = fft_rate + self.binwidth = float(sample_rate/fft_size) + self.average = average + self.cfunc = cfunc + self.xydfunc = xydfunc + self.interfunc = interfunc + if avg_alpha is None: + self.avg_alpha = 2.0 / fft_rate + else: + self.avg_alpha = avg_alpha + self.title = title + self.peak_hold = peak_hold + self.input_is_real = input_is_real + self.msgq = gr.msg_queue(2) # queue that holds a maximum of 2 messages + + def set_y_per_div(self, y_per_div): + self.y_per_div = y_per_div + + + def set_ref_level(self, ref_level): + self.ref_level = ref_level + + def set_average(self, average): + self.average = average + if average: + self.avg.set_taps(self.avg_alpha) + self.set_peak_hold(False) + else: + self.avg.set_taps(1.0) + + def set_peak_hold(self, enable): + self.peak_hold = enable + if enable: + self.set_average(False) + self.win.set_peak_hold(enable) + + def set_autoscale(self, auto): + self.autoscale = auto + + def set_avg_alpha(self, avg_alpha): + self.avg_alpha = avg_alpha + + def set_baseband_freq(self, baseband_freq): + self.baseband_freq = baseband_freq + + +class ra_fft_sink_f(gr.hier_block, ra_fft_sink_base): + def __init__(self, fg, parent, baseband_freq=0, + y_per_div=10, sc_y_per_div=0.5, sc_ref_level=40, ref_level=50, sample_rate=1, fft_size=512, + fft_rate=15, average=False, avg_alpha=None, title='', + size=default_ra_fftsink_size, peak_hold=False, cfunc=default_cfunc, xydfunc=None, interfunc=None): + + ra_fft_sink_base.__init__(self, input_is_real=True, baseband_freq=baseband_freq, + y_per_div=y_per_div, sc_y_per_div=sc_y_per_div, + sc_ref_level=sc_ref_level, ref_level=ref_level, + sample_rate=sample_rate, fft_size=fft_size, + fft_rate=fft_rate, + average=average, avg_alpha=avg_alpha, title=title, + peak_hold=peak_hold, cfunc=cfunc, + xydfunc=xydfunc, interfunc=interfunc) + + self.binwidth = float(sample_rate/2.0)/float(fft_size) + s2p = gr.serial_to_parallel(gr.sizeof_float, fft_size) + one_in_n = gr.keep_one_in_n(gr.sizeof_float * fft_size, + max(1, int(sample_rate/fft_size/fft_rate))) + mywindow = window.blackmanharris(fft_size) + fft = gr.fft_vfc(fft_size, True, mywindow) + c2mag = gr.complex_to_mag(fft_size) + self.avg = gr.single_pole_iir_filter_ff(1.0, fft_size) + log = gr.nlog10_ff(20, fft_size, -20*math.log10(fft_size)) + sink = gr.message_sink(gr.sizeof_float * fft_size, self.msgq, True) + + fg.connect (s2p, one_in_n, fft, c2mag, self.avg, log, sink) + gr.hier_block.__init__(self, fg, s2p, sink) + + self.win = fft_window(self, parent, size=size) + self.set_average(self.average) + +class ra_fft_sink_c(gr.hier_block, ra_fft_sink_base): + def __init__(self, fg, parent, baseband_freq=0, + y_per_div=10, sc_y_per_div=0.5, sc_ref_level=40, + ref_level=50, sample_rate=1, fft_size=512, + fft_rate=15, average=False, avg_alpha=None, title='', + size=default_ra_fftsink_size, peak_hold=False, cfunc=default_cfunc, xydfunc=None, interfunc=None): + + ra_fft_sink_base.__init__(self, input_is_real=False, baseband_freq=baseband_freq, + y_per_div=y_per_div, sc_y_per_div=sc_y_per_div, + sc_ref_level=sc_ref_level, ref_level=ref_level, + sample_rate=sample_rate, fft_size=fft_size, + fft_rate=fft_rate, + average=average, avg_alpha=avg_alpha, + title=title, + peak_hold=peak_hold, cfunc=cfunc, + xydfunc=xydfunc, interfunc=interfunc) + + s2p = gr.serial_to_parallel(gr.sizeof_gr_complex, fft_size) + one_in_n = gr.keep_one_in_n(gr.sizeof_gr_complex * fft_size, + max(1, int(sample_rate/fft_size/fft_rate))) + mywindow = window.blackmanharris(fft_size) + fft = gr.fft_vcc(fft_size, True, mywindow) + c2mag = gr.complex_to_mag(fft_size) + self.avg = gr.single_pole_iir_filter_ff(1.0, fft_size) + log = gr.nlog10_ff(20, fft_size, -20*math.log10(fft_size)) + sink = gr.message_sink(gr.sizeof_float * fft_size, self.msgq, True) + + fg.connect(s2p, one_in_n, fft, c2mag, self.avg, log, sink) + gr.hier_block.__init__(self, fg, s2p, sink) + + self.win = fft_window(self, parent, size=size) + self.set_average(self.average) + + +# ------------------------------------------------------------------------ + +myDATA_EVENT = wx.NewEventType() +EVT_DATA_EVENT = wx.PyEventBinder (myDATA_EVENT, 0) + + +class DataEvent(wx.PyEvent): + def __init__(self, data): + wx.PyEvent.__init__(self) + self.SetEventType (myDATA_EVENT) + self.data = data + + def Clone (self): + self.__class__ (self.GetId()) + + +class input_watcher (threading.Thread): + def __init__ (self, msgq, fft_size, event_receiver, **kwds): + threading.Thread.__init__ (self, **kwds) + self.setDaemon (1) + self.msgq = msgq + self.fft_size = fft_size + self.event_receiver = event_receiver + self.keep_running = True + self.start () + + def run (self): + while (self.keep_running): + msg = self.msgq.delete_head() # blocking read of message queue + itemsize = int(msg.arg1()) + nitems = int(msg.arg2()) + + s = msg.to_string() # get the body of the msg as a string + + # There may be more than one FFT frame in the message. + # If so, we take only the last one + if nitems > 1: + start = itemsize * (nitems - 1) + s = s[start:start+itemsize] + + complex_data = Numeric.fromstring (s, Numeric.Float32) + de = DataEvent (complex_data) + wx.PostEvent (self.event_receiver, de) + del de + + +class fft_window (plot.PlotCanvas): + def __init__ (self, ra_fftsink, parent, id = -1, + pos = wx.DefaultPosition, size = wx.DefaultSize, + style = wx.DEFAULT_FRAME_STYLE, name = ""): + plot.PlotCanvas.__init__ (self, parent, id, pos, size, style, name) + + self.y_range = None + self.ra_fftsink = ra_fftsink + self.peak_hold = False + self.peak_vals = None + + self.SetEnableGrid (True) + # self.SetEnableZoom (True) + # self.SetBackgroundColour ('black') + + self.build_popup_menu() + + EVT_DATA_EVENT (self, self.set_data) + wx.EVT_CLOSE (self, self.on_close_window) + self.Bind(wx.EVT_RIGHT_UP, self.on_right_click) + self.Bind(wx.EVT_MOTION, self.on_motion) + self.Bind(wx.EVT_LEFT_UP, self.on_left_click) + + self.input_watcher = input_watcher(ra_fftsink.msgq, ra_fftsink.fft_size, self) + + + def on_close_window (self, event): + print "fft_window:on_close_window" + self.keep_running = False + + + def set_data (self, evt): + calc_min = 99e10 + calc_max = -99e10 + dB = evt.data + L = len (dB) + + calc_min = min(dB) + calc_max = max(dB) + + dB = self.ra_fftsink.cfunc(dB, L) + + if self.peak_hold: + if self.peak_vals is None: + self.peak_vals = dB + else: + self.peak_vals = Numeric.maximum(dB, self.peak_vals) + dB = self.peak_vals + + x = max(abs(self.ra_fftsink.sample_rate), abs(self.ra_fftsink.baseband_freq)) + if x >= 1e9: + sf = 1e-9 + units = "GHz" + elif x >= 1e6: + sf = 1e-6 + units = "MHz" + elif x >= 1e3: + sf = 1e-3 + units = "kHz" + else: + sf = 1.0 + units = "Hz" + + if self.ra_fftsink.input_is_real: # only plot 1/2 the points + x_vals = ((Numeric.arrayrange (L/2) + * (self.ra_fftsink.sample_rate * sf / L)) + + self.ra_fftsink.baseband_freq * sf) + points = Numeric.zeros((len(x_vals), 2), Numeric.Float64) + points[:,0] = x_vals + points[:,1] = dB[0:L/2] + else: + # the "negative freqs" are in the second half of the array + x_vals = ((Numeric.arrayrange (-L/2, L/2) + * (self.ra_fftsink.sample_rate * sf / L)) + + self.ra_fftsink.baseband_freq * sf) + points = Numeric.zeros((len(x_vals), 2), Numeric.Float64) + points[:,0] = x_vals + points[:,1] = Numeric.concatenate ((dB[L/2:], dB[0:L/2])) + + lines = plot.PolyLine (points, colour='BLUE') + graphics = plot.PlotGraphics ([lines], + title=self.ra_fftsink.title, + xLabel = units, yLabel = "dB") + + self.Draw (graphics, xAxis=None, yAxis=self.y_range) + d = calc_max - calc_min + d *= 0.1 + if self.ra_fftsink.autoscale == True: + self.y_range = self._axisInterval ('min', calc_min-d, calc_max+d) + else: + self.update_y_range () + + def set_peak_hold(self, enable): + self.peak_hold = enable + self.peak_vals = None + + def update_y_range (self): + ymax = self.ra_fftsink.ref_level + ymin = self.ra_fftsink.ref_level - self.ra_fftsink.y_per_div * self.ra_fftsink.y_divs + self.y_range = self._axisInterval ('min', ymin, ymax) + + def on_average(self, evt): + # print "on_average" + self.ra_fftsink.set_average(evt.IsChecked()) + + def on_peak_hold(self, evt): + # print "on_peak_hold" + self.ra_fftsink.set_peak_hold(evt.IsChecked()) + + def on_autoscale(self, evt): + self.ra_fftsink.set_autoscale(evt.IsChecked()) + + def on_incr_ref_level(self, evt): + # print "on_incr_ref_level" + self.ra_fftsink.set_ref_level(self.ra_fftsink.ref_level + + self.ra_fftsink.y_per_div) + + def on_decr_ref_level(self, evt): + # print "on_decr_ref_level" + self.ra_fftsink.set_ref_level(self.ra_fftsink.ref_level + - self.ra_fftsink.y_per_div) + + def on_incr_y_per_div(self, evt): + # print "on_incr_y_per_div" + self.ra_fftsink.set_y_per_div(next_up(self.ra_fftsink.y_per_div, (0.5,1,2,5,10))) + + def on_decr_y_per_div(self, evt): + # print "on_decr_y_per_div" + self.ra_fftsink.set_y_per_div(next_down(self.ra_fftsink.y_per_div, (0.5,1,2,5,10))) + + def on_y_per_div(self, evt): + # print "on_y_per_div" + Id = evt.GetId() + if Id == self.id_y_per_div_1: + self.ra_fftsink.set_y_per_div(0.5) + elif Id == self.id_y_per_div_2: + self.ra_fftsink.set_y_per_div(1.0) + elif Id == self.id_y_per_div_5: + self.ra_fftsink.set_y_per_div(2.0) + elif Id == self.id_y_per_div_10: + self.ra_fftsink.set_y_per_div(5.0) + elif Id == self.id_y_per_div_20: + self.ra_fftsink.set_y_per_div(10) + + + def on_right_click(self, event): + menu = self.popup_menu + for id, pred in self.checkmarks.items(): + item = menu.FindItemById(id) + item.Check(pred()) + self.PopupMenu(menu, event.GetPosition()) + + def on_motion(self, event): + if not self.ra_fftsink.xydfunc == None: + xy = self.GetXY(event) + self.ra_fftsink.xydfunc (xy) + + def on_left_click(self, event): + if not self.ra_fftsink.interfunc == None: + xy = self.GetXY(event) + self.ra_fftsink.interfunc (xy[0]) + + def build_popup_menu(self): + self.id_incr_ref_level = wx.NewId() + self.id_decr_ref_level = wx.NewId() + self.id_autoscale = wx.NewId() + self.id_incr_y_per_div = wx.NewId() + self.id_decr_y_per_div = wx.NewId() + self.id_y_per_div_1 = wx.NewId() + self.id_y_per_div_2 = wx.NewId() + self.id_y_per_div_5 = wx.NewId() + self.id_y_per_div_10 = wx.NewId() + self.id_y_per_div_20 = wx.NewId() + self.id_average = wx.NewId() + self.id_peak_hold = wx.NewId() + + self.Bind(wx.EVT_MENU, self.on_average, id=self.id_average) + self.Bind(wx.EVT_MENU, self.on_peak_hold, id=self.id_peak_hold) + self.Bind(wx.EVT_MENU, self.on_autoscale, id=self.id_autoscale) + self.Bind(wx.EVT_MENU, self.on_incr_ref_level, id=self.id_incr_ref_level) + self.Bind(wx.EVT_MENU, self.on_decr_ref_level, id=self.id_decr_ref_level) + self.Bind(wx.EVT_MENU, self.on_incr_y_per_div, id=self.id_incr_y_per_div) + self.Bind(wx.EVT_MENU, self.on_decr_y_per_div, id=self.id_decr_y_per_div) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_1) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_2) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_5) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_10) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_20) + + + # make a menu + menu = wx.Menu() + self.popup_menu = menu + menu.AppendCheckItem(self.id_average, "Average") + menu.AppendCheckItem(self.id_peak_hold, "Peak Hold") + menu.Append(self.id_incr_ref_level, "Incr Ref Level") + menu.Append(self.id_decr_ref_level, "Decr Ref Level") + # menu.Append(self.id_incr_y_per_div, "Incr dB/div") + # menu.Append(self.id_decr_y_per_div, "Decr dB/div") + menu.AppendSeparator() + # we'd use RadioItems for these, but they're not supported on Mac + menu.AppendCheckItem(self.id_autoscale, "Autoscale") + menu.AppendCheckItem(self.id_y_per_div_1, "0.5 dB/div") + menu.AppendCheckItem(self.id_y_per_div_2, "1.0 dB/div") + menu.AppendCheckItem(self.id_y_per_div_5, "2.0 dB/div") + menu.AppendCheckItem(self.id_y_per_div_10, "5.0 dB/div") + menu.AppendCheckItem(self.id_y_per_div_20, "10.0 dB/div") + + self.checkmarks = { + self.id_average : lambda : self.ra_fftsink.average, + self.id_peak_hold : lambda : self.ra_fftsink.peak_hold, + self.id_autoscale : lambda : self.ra_fftsink.autoscale, + self.id_y_per_div_1 : lambda : self.ra_fftsink.y_per_div == 0.5, + self.id_y_per_div_2 : lambda : self.ra_fftsink.y_per_div == 1.0, + self.id_y_per_div_5 : lambda : self.ra_fftsink.y_per_div == 2.0, + self.id_y_per_div_10 : lambda : self.ra_fftsink.y_per_div == 5.0, + self.id_y_per_div_20 : lambda : self.ra_fftsink.y_per_div == 10.0, + } + + +def next_up(v, seq): + """ + Return the first item in seq that is > v. + """ + for s in seq: + if s > v: + return s + return v + +def next_down(v, seq): + """ + Return the last item in seq that is < v. + """ + rseq = list(seq[:]) + rseq.reverse() + + for s in rseq: + if s < v: + return s + return v + + +# ---------------------------------------------------------------- +# Deprecated interfaces +# ---------------------------------------------------------------- + +# returns (block, win). +# block requires a single input stream of float +# win is a subclass of wxWindow + +def make_ra_fft_sink_f(fg, parent, title, fft_size, input_rate, ymin = 0, ymax=50): + + block = ra_fft_sink_f(fg, parent, title=title, fft_size=fft_size, sample_rate=input_rate, + y_per_div=(ymax - ymin)/8, ref_level=ymax) + return (block, block.win) + +# returns (block, win). +# block requires a single input stream of gr_complex +# win is a subclass of wxWindow + +def make_ra_fft_sink_c(fg, parent, title, fft_size, input_rate, ymin=0, ymax=50): + block = ra_fft_sink_c(fg, parent, title=title, fft_size=fft_size, sample_rate=input_rate, + y_per_div=(ymax - ymin)/8, ref_level=ymax) + return (block, block.win) + + +# ---------------------------------------------------------------- +# Standalone test app +# ---------------------------------------------------------------- + +class test_app_flow_graph (stdgui.gui_flow_graph): + def __init__(self, frame, panel, vbox, argv): + stdgui.gui_flow_graph.__init__ (self, frame, panel, vbox, argv) + + fft_size = 256 + + # build our flow graph + input_rate = 20.000e3 + + # Generate a complex sinusoid + src1 = gr.sig_source_c (input_rate, gr.GR_SIN_WAVE, 5.75e3, 1000) + #src1 = gr.sig_source_c (input_rate, gr.GR_CONST_WAVE, 5.75e3, 1000) + + # We add these throttle blocks so that this demo doesn't + # suck down all the CPU available. Normally you wouldn't use these. + thr1 = gr.throttle(gr.sizeof_gr_complex, input_rate) + + sink1 = ra_fft_sink_c (self, panel, title="Complex Data", fft_size=fft_size, + sample_rate=input_rate, baseband_freq=100e3, + ref_level=60, y_per_div=10) + vbox.Add (sink1.win, 1, wx.EXPAND) + self.connect (src1, thr1, sink1) + + src2 = gr.sig_source_f (input_rate, gr.GR_SIN_WAVE, 5.75e3, 1000) + #src2 = gr.sig_source_f (input_rate, gr.GR_CONST_WAVE, 5.75e3, 1000) + thr2 = gr.throttle(gr.sizeof_float, input_rate) + sink2 = ra_fft_sink_f (self, panel, title="Real Data", fft_size=fft_size*2, + sample_rate=input_rate, baseband_freq=100e3, + ref_level=60, y_per_div=10) + vbox.Add (sink2.win, 1, wx.EXPAND) + self.connect (src2, thr2, sink2) + +def main (): + app = stdgui.stdapp (test_app_flow_graph, + "FFT Sink Test App") + app.MainLoop () + +if __name__ == '__main__': + main () diff --git a/gr-radio-astronomy/src/python/ra_stripchartsink.py b/gr-radio-astronomy/src/python/ra_stripchartsink.py new file mode 100755 index 000000000..112a0d698 --- /dev/null +++ b/gr-radio-astronomy/src/python/ra_stripchartsink.py @@ -0,0 +1,399 @@ +#!/usr/bin/env python +# +# Copyright 2003,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. +# + +from gnuradio import gr, gru +from gnuradio.wxgui import stdgui +import wx +import gnuradio.wxgui.plot as plot +import Numeric +import threading +import math +import ephem +import time + +default_stripchartsink_size = (640,140) +global_yvalues = [] + +def default_cfunc(datum): + return(datum) + +class stripchart_sink_base(object): + def __init__(self, input_is_real=False, y_per_div=10, ref_level=50, + sample_rate=1, stripsize=4, + title='',xlabel="X", ylabel="Y", divbase=0.025, + cfunc=default_cfunc, parallel=False, scaling=1.0, autoscale=False): + + # initialize common attributes + self.y_divs = 8 + self.y_per_div=y_per_div + self.ref_level = ref_level + self.autoscale = autoscale + self.sample_rate = sample_rate + self.parallel = parallel + self.title = title + self.xlabel = xlabel + self.ylabel = ylabel + self.divbase = divbase + self.scaling = scaling + self.cfunc = cfunc + self.input_is_real = input_is_real + self.msgq = gr.msg_queue(2) # queue that holds a maximum of 2 messages + self.vector=Numeric.zeros(stripsize,Numeric.Float64) + self.wcnt = 0 + self.timecnt = 0 + self.stripsize=stripsize + + def set_y_per_div(self, y_per_div): + self.y_per_div = y_per_div + + def set_ref_level(self, ref_level): + self.ref_level = ref_level + + def set_autoscale(self, auto): + self.autoscale = auto + +class stripchart_sink_f(gr.hier_block, stripchart_sink_base): + def __init__(self, fg, parent, + y_per_div=10, ref_level=50, sample_rate=1, + title='', stripsize=4, + size=default_stripchartsink_size,xlabel="X", + ylabel="Y", divbase=0.025, cfunc=default_cfunc, + parallel=False, scaling=1.0, autoscale=False): + + stripchart_sink_base.__init__(self, input_is_real=True, + y_per_div=y_per_div, ref_level=ref_level, + sample_rate=sample_rate, + stripsize=stripsize, + xlabel=xlabel, ylabel=ylabel, + divbase=divbase, title=title, + cfunc=cfunc, parallel=parallel, + scaling=scaling, autoscale=autoscale) + + if (parallel == True): + one = gr.keep_one_in_n (gr.sizeof_float*stripsize, 1) + sink = gr.message_sink(gr.sizeof_float*stripsize, self.msgq, True) + else: + one = gr.keep_one_in_n (gr.sizeof_float, 1) + sink = gr.message_sink(gr.sizeof_float, self.msgq, True) + fg.connect (one, sink) + + gr.hier_block.__init__(self, fg, one, sink) + + self.win = stripchart_window(self, parent, size=size) + + + +# ------------------------------------------------------------------------ + +myDATA_EVENT = wx.NewEventType() +EVT_DATA_EVENT = wx.PyEventBinder (myDATA_EVENT, 0) + + +class DataEvent(wx.PyEvent): + def __init__(self, data): + wx.PyEvent.__init__(self) + self.SetEventType (myDATA_EVENT) + self.data = data + + def Clone (self): + self.__class__ (self.GetId()) + + +class input_watcher (threading.Thread): + def __init__ (self, msgq, evsize, event_receiver, **kwds): + threading.Thread.__init__ (self, **kwds) + self.setDaemon (1) + self.msgq = msgq + self.evsize = evsize + self.event_receiver = event_receiver + self.keep_running = True + self.start () + + def run (self): + while (self.keep_running): + msg = self.msgq.delete_head() # blocking read of message queue + itemsize = int(msg.arg1()) + nitems = int(msg.arg2()) + + s = msg.to_string() # get the body of the msg as a string + + # There may be more than one frame in the message. + # If so, we take only the last one + if nitems > 1: + start = itemsize * (nitems - 1) + s = s[start:start+itemsize] + + complex_data = Numeric.fromstring (s, Numeric.Float32) + de = DataEvent (complex_data) + wx.PostEvent (self.event_receiver, de) + del de + +class stripchart_window(plot.PlotCanvas): + def __init__ (self, stripchartsink, parent, id = -1, + pos = wx.DefaultPosition, size = wx.DefaultSize, + style = wx.DEFAULT_FRAME_STYLE, name = ""): + plot.PlotCanvas.__init__ (self, parent, id, pos, size, style, name) + + self.y_range = None + self.stripchartsink = stripchartsink + + self.SetEnableGrid (True) + # self.SetEnableZoom (True) + # self.SetBackgroundColour ('black') + + self.build_popup_menu() + + EVT_DATA_EVENT (self, self.set_data) + + wx.EVT_CLOSE (self, self.on_close_window) + self.Bind(wx.EVT_RIGHT_UP, self.on_right_click) + + self.input_watcher = input_watcher(stripchartsink.msgq, 1, self) + + + def on_close_window (self, event): + print "stripchart_window:on_close_window" + self.keep_running = False + + + def set_data (self, evt): + indata = evt.data + L = len (indata) + + calc_min = min(indata) + calc_max = max(indata) + d = calc_max - calc_min + d *= 0.1 + if self.stripchartsink.autoscale == True and self.stripchartsink.parallel == True: + self.y_range = self._axisInterval ('min', calc_min-d, calc_max+d) + + if (self.stripchartsink.parallel != True): + indata = self.stripchartsink.cfunc(indata) + + N = self.stripchartsink.stripsize + if self.stripchartsink.parallel != True: + for i in range(1,N): + pooey = N-i + self.stripchartsink.vector[pooey] = self.stripchartsink.vector[pooey-1] + + self.stripchartsink.vector[0] = indata + + else: + self.stripchartsink.vector = indata + + if self.stripchartsink.parallel == True: + avg = 0 + for i in range(0,self.stripchartsink.stripsize): + if self.stripchartsink.vector[i] > 0: + avg += self.stripchartsink.vector[i] + if self.stripchartsink.vector[i] < calc_min: + calc_min = self.stripchartsink.vector[i] + if self.stripchartsink.vector[i] > calc_max: + calc_max = self.stripchartsink.vector[i] + + avg /= self.stripchartsink.stripsize + markers = [] + placedmarkers = 0 + for i in range(0,self.stripchartsink.stripsize): + if (self.stripchartsink.vector[i] > 0 and + self.stripchartsink.vector[i] > (avg*5)): + markers.append((i*self.stripchartsink.scaling, + self.stripchartsink.vector[i])) + placedmarkers += 1 + + points = Numeric.zeros((N,2), Numeric.Float64) + for i in range(0,N): + if self.stripchartsink.scaling == 1.0: + points[i,0] = i + else: + points[i,0] = i * self.stripchartsink.scaling + points[i,1] = self.stripchartsink.vector[i] + + if self.stripchartsink.parallel == True and placedmarkers > 1: + for i in range(0,N): + self.stripchartsink.vector[i] = 0 + + marks = plot.PolyMarker(markers, colour='BLACK', marker='triangle_down') + + lines = plot.PolyLine (points, colour='RED') + + # Temporary--I'm find the markers distracting + placedmarkers = 0 + xlab = self.stripchartsink.xlabel + ylab = self.stripchartsink.ylabel + if (self.stripchartsink.parallel == False) or (placedmarkers <= 1): + graphics = plot.PlotGraphics ([lines], + title=self.stripchartsink.title, + xLabel = xlab, yLabel = ylab) + + else: + graphics = plot.PlotGraphics ([lines,marks], + title=self.stripchartsink.title, + xLabel = xlab, yLabel = ylab) + + self.Draw (graphics, xAxis=None, yAxis=self.y_range) + + if self.stripchartsink.autoscale == False or self.stripchartsink.parallel == False: + self.update_y_range () + + + def update_y_range (self): + ymax = self.stripchartsink.ref_level + ymin = self.stripchartsink.ref_level - self.stripchartsink.y_per_div * self.stripchartsink.y_divs + self.y_range = self._axisInterval ('min', ymin, ymax) + + def on_incr_ref_level(self, evt): + # print "on_incr_ref_level" + self.stripchartsink.set_ref_level(self.stripchartsink.ref_level + + self.stripchartsink.y_per_div) + + def on_decr_ref_level(self, evt): + # print "on_decr_ref_level" + self.stripchartsink.set_ref_level(self.stripchartsink.ref_level + - self.stripchartsink.y_per_div) + + def on_autoscale(self, evt): + self.stripchartsink.set_autoscale(evt.IsChecked()) + + def on_incr_y_per_div(self, evt): + divbase = self.stripchartsink.divbase + x1 = 1 * divbase + x2 = 2 * divbase + x4 = 4 * divbase + x10 = 10 * divbase + x20 = 20 * divbase + # print "on_incr_y_per_div" + self.stripchartsink.set_y_per_div(next_up(self.stripchartsink.y_per_div, (x1,x2,x4,x10,x20))) + + def on_decr_y_per_div(self, evt): + # print "on_decr_y_per_div" + divbase = self.stripchartsink.divbase + x1 = 1 * divbase + x2 = 2 * divbase + x4 = 4 * divbase + x10 = 10 * divbase + x20 = 20 * divbase + self.stripchartsink.set_y_per_div(next_down(self.stripchartsink.y_per_div, (x1,x2,x4,x10,x20))) + + def on_y_per_div(self, evt): + # print "on_y_per_div" + divbase=self.stripchartsink.divbase + Id = evt.GetId() + if Id == self.id_y_per_div_1: + self.stripchartsink.set_y_per_div(1*divbase) + elif Id == self.id_y_per_div_2: + self.stripchartsink.set_y_per_div(2*divbase) + elif Id == self.id_y_per_div_5: + self.stripchartsink.set_y_per_div(4*divbase) + elif Id == self.id_y_per_div_10: + self.stripchartsink.set_y_per_div(10*divbase) + elif Id == self.id_y_per_div_20: + self.stripchartsink.set_y_per_div(20*divbase) + + + def on_right_click(self, event): + menu = self.popup_menu + for id, pred in self.checkmarks.items(): + item = menu.FindItemById(id) + item.Check(pred()) + self.PopupMenu(menu, event.GetPosition()) + + + def build_popup_menu(self): + divbase=self.stripchartsink.divbase + self.id_incr_ref_level = wx.NewId() + self.id_decr_ref_level = wx.NewId() + self.id_autoscale = wx.NewId() + self.id_incr_y_per_div = wx.NewId() + self.id_decr_y_per_div = wx.NewId() + self.id_y_per_div_1 = wx.NewId() + self.id_y_per_div_2 = wx.NewId() + self.id_y_per_div_5 = wx.NewId() + self.id_y_per_div_10 = wx.NewId() + self.id_y_per_div_20 = wx.NewId() + + self.Bind(wx.EVT_MENU, self.on_incr_ref_level, id=self.id_incr_ref_level) + self.Bind(wx.EVT_MENU, self.on_decr_ref_level, id=self.id_decr_ref_level) + self.Bind(wx.EVT_MENU, self.on_autoscale, id=self.id_autoscale) + self.Bind(wx.EVT_MENU, self.on_incr_y_per_div, id=self.id_incr_y_per_div) + self.Bind(wx.EVT_MENU, self.on_decr_y_per_div, id=self.id_decr_y_per_div) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_1) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_2) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_5) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_10) + self.Bind(wx.EVT_MENU, self.on_y_per_div, id=self.id_y_per_div_20) + + + # make a menu + menu = wx.Menu() + self.popup_menu = menu + menu.Append(self.id_incr_ref_level, "Incr Ref Level") + menu.Append(self.id_decr_ref_level, "Decr Ref Level") + menu.AppendSeparator() + menu.AppendCheckItem(self.id_autoscale, "Auto Scale") + # we'd use RadioItems for these, but they're not supported on Mac + v = 1.0*divbase + s = "%.3f" % v + menu.AppendCheckItem(self.id_y_per_div_1, s) + v = 2.0*divbase + s = "%.3f" % v + menu.AppendCheckItem(self.id_y_per_div_2, s) + v = 4.0*divbase + s = "%.3f" % v + menu.AppendCheckItem(self.id_y_per_div_5, s) + v = 10*divbase + s = "%.3f" % v + menu.AppendCheckItem(self.id_y_per_div_10, s) + v = 20*divbase + s = "%.3f" % v + menu.AppendCheckItem(self.id_y_per_div_20, s) + + self.checkmarks = { + self.id_autoscale : lambda : self.stripchartsink.autoscale, + self.id_y_per_div_1 : lambda : self.stripchartsink.y_per_div == 1*divbase, + self.id_y_per_div_2 : lambda : self.stripchartsink.y_per_div == 2*divbase, + self.id_y_per_div_5 : lambda : self.stripchartsink.y_per_div == 4*divbase, + self.id_y_per_div_10 : lambda : self.stripchartsink.y_per_div == 10*divbase, + self.id_y_per_div_20 : lambda : self.stripchartsink.y_per_div == 20*divbase, + } + + +def next_up(v, seq): + """ + Return the first item in seq that is > v. + """ + for s in seq: + if s > v: + return s + return v + +def next_down(v, seq): + """ + Return the last item in seq that is < v. + """ + rseq = list(seq[:]) + rseq.reverse() + + for s in rseq: + if s < v: + return s + return v diff --git a/gr-radio-astronomy/src/python/run_tests.in b/gr-radio-astronomy/src/python/run_tests.in new file mode 100644 index 000000000..51568827c --- /dev/null +++ b/gr-radio-astronomy/src/python/run_tests.in @@ -0,0 +1,47 @@ +#!/bin/sh + +# All this strange PYTHONPATH manipulation is required to run our +# tests using our just built shared library and swig-generated python +# code prior to installation. + +# build tree == src tree unless you're doing a VPATH build. +# If you don't know what a VPATH build is, you're not doing one. Relax... + +prefix=@prefix@ +exec_prefix=@exec_prefix@ + +# Where to look in the build tree for our shared library +libbld=@abs_top_builddir@/gr-radio-astronomy/src/lib +# Where to look in the src tree for swig generated python code +libsrc=@abs_top_srcdir@/gr-radio-astronomy/src/lib +# Where to look in the src tree for hand written python code +py=@abs_top_srcdir@/gr-radio-astronomy/src/python + +# Where to look for GNU Radio python modules in current build tree +# FIXME this is wrong on a distcheck. We really need to ask gnuradio-core +# where it put its python files. +grpythonbld=@abs_top_builddir@/gnuradio-core/src/python/:@abs_top_builddir@/gnuradio-core/src/lib/swig/:@abs_top_builddir@/gnuradio-core/src/lib/swig/.libs + +PYTHONPATH="$grpythonbld:$libbld:$libbld/.libs:$libsrc:$py:$PYTHONPATH" +export PYTHONPATH + +# +# This is the simple part... +# Run everything that matches qa_*.py and return the final result. +# + +ok=yes +for file in @srcdir@/qa_*.py +do + if ! $file + then + ok=no + fi +done + +if [ $ok = yes ] +then + exit 0 +else + exit 1 +fi diff --git a/gr-radio-astronomy/src/python/usrp_psr_receiver.help b/gr-radio-astronomy/src/python/usrp_psr_receiver.help new file mode 100644 index 000000000..5801f3fbb --- /dev/null +++ b/gr-radio-astronomy/src/python/usrp_psr_receiver.help @@ -0,0 +1,111 @@ +This program is used to analyse pulsars of known parameters. It contains + both a post-detector spectral display, and a "pulse profile" display. + It has a built-in de-dispersion filter that will work up to DM=100 for + 21cm observing, and up to DM=5 for 327Mhz observing. + +The program takes the following options: + + --rx-subdev-spec which USRP Rx side? A or B + + --decim USRP decimation rate use either 64 or 128 + + --freq USRP daughtercard frequency + + --observing Actual observing frequency (default is to use the + setting for --freq) + + --avg Averaging setting for spectral display--higher numbers + equal more averaging. 25 to 40 is typical. + + --favg Pulse folding averaging. 2 to 5 is typical. + + --gain USRP daughtercard gain control + + --reflevel Reference level on pulse profile display + + --lowest Lowest spectral bin that is considered valid, in Hz + + --longitude Observer longitude: West is negative + + --latitude Observer latitude: South is negative + + --fft_size Size of FFT for post-detector spectrum: default is 1024 + + --threshold Threshold (dB) to be considered a spectral "peak" + This is relative to the average spectral level + + --lowpass Low pass frequency for post-detector spectral display + 20-100 is typical + + --prefix Filename prefix to use for recording files + Default is ./ + + --pulsefreq The frequency of the expected pulses + For sentimental reasons, this defaults to 0.748Hz + + --dm The DM + + --doppler The doppler shift, as a ratio + + --divbase The base of the Y/Div menu in pulsar display + + --division The initial Y/Div in pulsar display + +DM, Doppler, Gain, Frequency, and the averaging parameters can all be + changed using the GUI at runtime. + +If latitude and longitude are set correctly, and the system time is + correct, then the current LMST is displayed below the frequency + input, updated once per second. + +Moving the mouse in the post-detector spectrum display shows you that + point in the post-detector spectrum, both frequency and signal level. + +The post-detector spectrum is analysed, with results shown below + "Best freq". It shows the spectral peaks, and computes their relationship. + It shows the harmonic compliance among the peaks, as well as the average + peak-to-peak distance. + + +Here's a complete example for observing a pulsar with a frequency of + 1.35Hz, at 431.5Mhz, using an IF of 10.7Mhz, and a DM of 12.431, using + 1Mhz observing bandwidth: + +./usrp_psr_receiver.py --freq 10.7e6 --decim 64 --dm 12.431 --avg 35 \ + --pulsefreq 1.35 --fft_size 2048 --lowest 1.00 --gain 75 --threshold 11.5 \ + --observing 431.5e6 --reflevel 200 --division 100 --divbase 10 --favg 3 \ + --lowpass 20 --longitude -76.02 --latitude 44.95 + +Since the observed pulsar is at 1.35Hz, a lowpass cutoff for the + post-detector spectral display of 20Hz will be adequate. We + tell the spectral analyser to use a threshold of 11.5dB above + average when analysing spectral data, and set the epoch folder + averager (pulse profile display) to use an average from 3 samples. + Notice that our actual USRP/Daughtercard frequency is 10.7Mhz, while + our observing frequency is 431.5Mhz--this is important in order for + the DM de-dispersion calculations to be correct. We also set our + latitude and longitude, so that logfiles and the LMST display + will have the correct LMST in them. + +The entire complex baseband can be recorded, if the "Recording baseband" + button is pressed. Filenames are generated dynamically, and a header + file is produced giving observation parameters. The baseband data are + recorded as octet pairs: one for I and one for Q. Pressing the button again + turns off baseband recording. This baseband is "raw", so it will + not have been de-dispersed. The data rate will be whatever the + USRP was programmed to at the time (based on --decim). + + The files are: YYYYMMDDHHMM.pdat and YYYYMMDDHHMM.phdr + + The .phdr file contains ASCII header information describing the + contents of the .pdat file. + +Similarly the raw, pre-folded, band-limited post-detector "audio" data can be + recorded using the "Record Pulses" button. The data rate for these is + currently 20Khz, recorded as short integers. Just like baseband recording, + pressing the button again turns off pulse recording. + + The files are: YYYYMMDDHHMM.padat and YYMMDDHHMM.pahdr + + The .pahdr file is ascii text providing information about the contents + of the corresponding .padat file. 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 () diff --git a/gr-radio-astronomy/src/python/usrp_ra_receiver.help b/gr-radio-astronomy/src/python/usrp_ra_receiver.help new file mode 100644 index 000000000..45a21e297 --- /dev/null +++ b/gr-radio-astronomy/src/python/usrp_ra_receiver.help @@ -0,0 +1,90 @@ +This program is used to take spectra and total power measurements. + It records spectral and total-power data to external datalogging + files. + +The program takes the following options: + + --rx-subdev-spec which USRP Rx side? A or B + + --decim USRP decimation rate: 8, 16, 32, and 64 are good + (8Mhz, 4Mhz, 2Mhz, and 1Mhz bandwidth) + + --freq USRP daughtercard frequency + + --observing Actual observing frequency (default is to use the + setting for --freq) + + --avg Averaging setting for spectral display--higher numbers + equal more averaging. 25 to 40 is typical. + + --integ Total power integration time: seconds + + --gain USRP daughtercard gain control + + --reflevel Reference level on pulse profile display + + --longitude Observer longitude: West is negative + + --latitude Observer latitude: South is negative + + --fft_size Size of FFT for post-detector spectrum: default is 1024 + + --prefix Filename prefix to use for data logging files + Default is ./ + + --divbase The base of the Y/Div menu in pulsar display + + --division The initial Y/Div in pulsar display + + --ylabel Y axis label + + --cfunc The function name prefix for the spectral and + calibrator functions + + --waterfall Use waterfall, rather than regular spectral display + NOT TESTED IN THIS APPLICATION + + --stripsize Size of the total-power stripchart, in samples + +There are two windows--a spectral window, and the total-power window. + Moving the cursor around in the spectral window shows you the + corresponding frequency and doppler shift. Left clicking in this + window sets an interference marker, which sets a "zero" in the + interference filter. Use the "clear interference" button to clear this. + +The total power window is updated at a fixed 2Hz rate, and grows from + the left of the display. + +If latitude and longitude are set correctly, and system time is correct, + then the current LMST is displayed, updated once per second. + +Averaging parameters, gain, and frequency can all be set from the GUI using + the appropriate controls. You can also enter the current declination, which + will appear in the datalogging files. This is useful both for mapping, + and housekeeping purposes, particularly when you haven't looked at a datafile + for quite some time. + +There are two datalog files produced by this program: + + YYYYMMDDHH.tpdat Total power data + + The date/time portion of the filename is referred to local time, + rather than UTC or sidereal. + + First field is sidereal time when sample was taken + Second field is total power datum + Third field is declination in decimal degrees + + Samples are written once per second + + YYYYMMDDHH.sdat Spectral data + + The date/time portion of the filename is referred to local time, + rather than UTC or sidereal. + + First field is sidereal time when spectrum was taken + Second field is declination in decimal degrees + Third field is complex spectral data--in the same order that FFTW3 library + places bins: DC to bandwidth/2, followed by -bandwidth/2 to DC. + + Spectral snapshots are written once every 5 seconds diff --git a/gr-radio-astronomy/src/python/usrp_ra_receiver.py b/gr-radio-astronomy/src/python/usrp_ra_receiver.py new file mode 100755 index 000000000..37a1ebfff --- /dev/null +++ b/gr-radio-astronomy/src/python/usrp_ra_receiver.py @@ -0,0 +1,584 @@ +#!/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. +# + +from gnuradio import gr, gru +from gnuradio import usrp +import usrp_dbid +from gnuradio import eng_notation +from gnuradio.eng_option import eng_option +from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, waterfallsink, form, slider +from optparse import OptionParser +import wx +import sys +from Numeric import * +import FFT +import ephem +from gnuradio.local_calibrator import * + +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("-a", "--avg", type="eng_float", default=1.0, + help="set spectral averaging alpha") + parser.add_option("-i", "--integ", type="eng_float", default=1.0, + help="set integration time") + 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 Total power reference level") + parser.add_option("-y", "--division", type="eng_float", default=0.5, + help="Set Total power Y division size") + 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("-o", "--observing", type="eng_float", default=0.0, + help="Set observing frequency") + parser.add_option("-x", "--ylabel", default="dB", help="Y axis label") + parser.add_option("-C", "--cfunc", default="default", help="Calibration function name") + parser.add_option("-z", "--divbase", type="eng_float", default=0.025, help="Y Division increment base") + parser.add_option("-v", "--stripsize", type="eng_float", default=2400, help="Size of stripchart, in 2Hz samples") + parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT") + + parser.add_option("-N", "--decln", type="eng_float", default=999.99, help="Observing declination") + parser.add_option("-I", "--interfilt", action="store_true", default=False) + parser.add_option("-X", "--prefix", default="./") + (options, args) = parser.parse_args() + if len(args) != 0: + parser.print_help() + sys.exit(1) + + self.show_debug_info = True + + # 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)) + self.cardtype = self.u.daughterboard_id(0) + # Set initial declination + self.decln = options.decln + + # Turn off interference filter by default + self.use_interfilt = options.interfilt + + # determine the daughterboard subdevice we're using + self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec) + + input_rate = self.u.adc_freq() / self.u.decim_rate() + + tpstr="calib_"+options.cfunc+"_total_power" + sstr="calib_"+options.cfunc+"_fft" + self.tpcfunc=eval(tpstr) + self.scfunc=eval(sstr) + + # + # Set prefix for data files + # + self.prefix = options.prefix + calib_set_prefix(self.prefix) + + # Set up FFT display + self.scope = ra_fftsink.ra_fft_sink_c (self, panel, + fft_size=int(options.fft_size), sample_rate=input_rate, + fft_rate=8, title="Spectral", + cfunc=self.scfunc, xydfunc=self.xydfunc, interfunc=self.interference) + + # Set up ephemeris data + self.locality = ephem.Observer() + self.locality.long = str(options.longitude) + self.locality.lat = str(options.latitude) + + # Set up stripchart display + self.stripsize = int(options.stripsize) + self.chart = ra_stripchartsink.stripchart_sink_f (self, panel, + stripsize=self.stripsize, + title="Continuum", + xlabel="LMST Offset (Seconds)", + scaling=1.0, ylabel=options.ylabel, + divbase=options.divbase, cfunc=self.tpcfunc) + + # Set center frequency + self.centerfreq = options.freq + + # Set observing frequency (might be different from actual programmed + # RF frequency) + if options.observing == 0.0: + self.observing = options.freq + else: + self.observing = options.observing + + self.bw = input_rate + + # + # Produce a default interference map + # May not actually get used, unless --interfilt was specified + # + self.intmap = Numeric.zeros(256,Numeric.Complex64) + for i in range(0,len(self.intmap)): + self.intmap[i] = complex(1.0, 0.0) + + # We setup the first two integrators to produce a fixed integration + # Down to 1Hz, with output at 1 samples/sec + N = input_rate/5000 + + # Second stage runs on decimated output of first + M = (input_rate/N) + + # Create taps for first integrator + t = range(0,N-1) + tapsN = [] + for i in t: + tapsN.append(1.0/N) + + # Create taps for second integrator + t = range(0,M-1) + tapsM = [] + for i in t: + tapsM.append(1.0/M) + + # + # The 3rd integrator is variable, and user selectable at runtime + # This integrator doesn't decimate, but is used to set the + # final integration time based on the constant 1Hz input samples + # The strip chart is fed at a constant 1Hz rate as a result + # + + # + # Call constructors for receive chains + # + + # + # This is the interference-zapping filter + # + # The GUI is used to set/clear inteference zones in + # the filter. The non-interfering zones are set to + # 1.0. + # + if 0: + self.interfilt = gr.fft_filter_ccc(1,self.intmap) + tmp = FFT.inverse_fft(self.intmap) + self.interfilt.set_taps(tmp) + + # The three integrators--two FIR filters, and an IIR final filter + self.integrator1 = gr.fir_filter_fff (N, tapsN) + self.integrator2 = gr.fir_filter_fff (M, tapsM) + self.integrator3 = gr.single_pole_iir_filter_ff(1.0) + + # Split complex USRP stream into a pair of floats + self.splitter = gr.complex_to_float (1); + self.toshort = gr.float_to_short(); + + # 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(); + + # + # Start connecting configured modules in the receive chain + # + + # Connect interference-filtered USRP input to selected scope function + if self.use_interfilt == True: + self.connect(self.u, self.interfilt, self.scope) + + # Connect interference-filtered USRP to a complex->float splitter + self.connect(self.interfilt, self.splitter) + + else: + self.connect(self.u, self.scope) + self.connect(self.u, 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 adder output to three-stages of FIR integrator + self.connect(self.adder, self.integrator1, + self.integrator2, self.integrator3, self.chart) + + + self._build_gui(vbox) + + # Make GUI agree with command-line + self.myform['integration'].set_value(int(options.integ)) + self.myform['average'].set_value(int(options.avg)) + + # Make integrator agree with command line + self.set_integration(int(options.integ)) + + # Make spectral averager agree with command line + if options.avg != 1.0: + self.scope.set_avg_alpha(float(1.0/options.avg)) + calib_set_avg_alpha(float(options.avg)) + self.scope.set_average(True) + + + # Set division size + self.chart.set_y_per_div(options.division) + + # Set reference(MAX) level + self.chart.set_ref_level(options.reflevel) + + # 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 + + # Set the initial gain control + self.set_gain(options.gain) + + if not(self.set_freq(options.freq)): + self._set_status_msg("Failed to set initial frequency") + + self.set_decln (self.decln) + calib_set_bw(self.decln) + + 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()) + + # Make sure calibrator knows what our bandwidth is + calib_set_bw(self.u.adc_freq() / self.u.decim_rate()) + + # Tell calibrator our declination as well + calib_set_decln(self.decln) + + # Start the timer for the LMST display + 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_decln(kv): + return self.set_decln(kv['decln']) + + # Position the FFT display + vbox.Add(self.scope.win, 15, wx.EXPAND) + + # Position the Total-power stripchart + vbox.Add(self.chart.win, 15, 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((4,0), 0, 0) + + myform['lmst_high'] = form.static_text_field( + parent=self.panel, sizer=vbox1, label="Current LMST", weight=1) + vbox1.Add((4,0), 0, 0) + + myform['spec_data'] = form.static_text_field( + parent=self.panel, sizer=vbox1, label="Spectral Cursor", weight=1) + vbox1.Add((4,0), 0, 0) + + 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((4,0), 0, 0) + myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2, + label="Spectral Averaging (FFT frames)", weight=1, min=1, max=2000, callback=self.set_averaging) + + vbox2.Add((4,0), 0, 0) + + myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2, + label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration) + + vbox2.Add((4,0), 0, 0) + myform['decln'] = form.float_field( + parent=self.panel, sizer=vbox2, label="Current Declination", weight=1, + callback=myform.check_input_and_call(_form_set_decln)) + vbox2.Add((4,0), 0, 0) + + buttonbox = wx.BoxSizer(wx.HORIZONTAL) + if self.use_interfilt == True: + self.doit = form.button_with_callback(self.panel, + label="Clear Interference List", + callback=self.clear_interferers) + if self.use_interfilt == True: + buttonbox.Add(self.doit, 0, wx.CENTER) + vbox.Add(buttonbox, 0, wx.CENTER) + hbox.Add(vbox1, 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. + """ + # + # Everything except BASIC_RX should support usrp.tune() + # + if not (self.cardtype == usrp_dbid.BASIC_RX): + r = usrp.tune(self.u, 0, self.subdev, target_freq) + else: + r = self.u.set_rx_freq(0, target_freq) + f = self.u.rx_freq(0) + if abs(f-target_freq) > 2.0e3: + r = 0 + if r: + self.myform['freq'].set_value(target_freq) # update displayed value + # + # Make sure calibrator knows our target freq + # + + # Remember centerfreq---used for doppler calcs + delta = self.centerfreq - target_freq + self.centerfreq = target_freq + self.observing -= delta + self.scope.set_baseband_freq (self.observing) + calib_set_freq(self.observing) + + # Clear interference list + self.clear_interferers() + + self.myform['baseband'].set_value(r.baseband_freq) + self.myform['ddc'].set_value(r.dxc_freq) + + return True + + return False + + def set_decln(self, dec): + self.decln = dec + self.myform['decln'].set_value(dec) # update displayed value + calib_set_decln(dec) + + def set_gain(self, gain): + self.myform['gain'].set_value(gain) # update displayed value + self.subdev.set_gain(gain) + + # + # Make sure calibrator knows our gain setting + # + calib_set_gain(gain) + + def set_averaging(self, avval): + self.myform['average'].set_value(avval) + self.scope.set_avg_alpha(1.0/(avval)) + calib_set_avg_alpha(avval) + self.scope.set_average(True) + + def set_integration(self, integval): + self.integrator3.set_taps(1.0/integval) + self.myform['integration'].set_value(integval) + + # + # Make sure calibrator knows our integration time + # + calib_set_integ(integval) + + def lmst_timeout(self): + self.locality.date = ephem.now() + sidtime = self.locality.sidereal_time() + self.myform['lmst_high'].set_value(str(ephem.hours(sidtime))) + + def xydfunc(self,xyv): + magn = int(log10(self.observing)) + if (magn == 6 or magn == 7 or magn == 8): + magn = 6 + dfreq = xyv[0] * pow(10.0,magn) + ratio = self.observing / dfreq + vs = 1.0 - ratio + vs *= 299792.0 + if magn >= 9: + xhz = "Ghz" + elif magn >= 6: + xhz = "Mhz" + elif magn <= 5: + xhz = "Khz" + s = "%.6f%s\n%.3fdB" % (xyv[0], xhz, xyv[1]) + s2 = "\n%.3fkm/s" % vs + self.myform['spec_data'].set_value(s+s2) + + def interference(self,x): + if self.use_interfilt == False: + return + magn = int(log10(self.observing)) + dfreq = x * pow(10.0,magn) + delta = dfreq - self.observing + fincr = self.bw / len(self.intmap) + l = len(self.intmap) + if delta > 0: + offset = delta/fincr + else: + offset = (l) - int((abs(delta)/fincr)) + + offset = int(offset) + + if offset >= len(self.intmap) or offset < 0: + print "interference offset is invalid--", offset + return + + # + # Zero out the region around the selected interferer + # + self.intmap[offset-2] = complex (0.5, 0.0) + self.intmap[offset-1] = complex (0.25, 0.0) + self.intmap[offset] = complex (0.0, 0.0) + self.intmap[offset+1] = complex(0.25, 0.0) + self.intmap[offset+2] = complex(0.5, 0.0) + + # + # Set new taps + # + tmp = FFT.inverse_fft(self.intmap) + self.interfilt.set_taps(tmp) + + def clear_interf(self): + self.clear_interferers() + + def clear_interferers(self): + for i in range(0,len(self.intmap)): + self.intmap[i] = complex(1.0,0.0) + tmp = FFT.inverse_fft(self.intmap) + if self.use_interfilt == True: + self.interfilt.set_taps(tmp) + + + + def toggle_cal(self): + if (self.calstate == True): + self.calstate = False + self.u.write_io(0,0,(1<<15)) + self.calibrator.SetLabel("Calibration Source: Off") + else: + self.calstate = True + self.u.write_io(0,(1<<15),(1<<15)) + self.calibrator.SetLabel("Calibration Source: On") + + def toggle_annotation(self): + if (self.annotate_state == True): + self.annotate_state = False + self.annotation.SetLabel("Annotation: Off") + else: + self.annotate_state = True + self.annotation.SetLabel("Annotation: On") + calib_set_interesting(self.annotate_state) + + +def main (): + app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY SPECTRAL/CONTINUUM RECEIVER: $Revision$", nstatus=1) + app.MainLoop() + +if __name__ == '__main__': + main () |