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

   1 #!/usr/bin/env python
   2 #
   3 # Copyright 2004,2005 Free Software Foundation, Inc.
   4 #
   5 # This file is part of GNU Radio
   6 #
   7 # GNU Radio is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 2, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23 from gnuradio import gr, gru
  24 from gnuradio import usrp
  25 import usrp_dbid
  26 from gnuradio import eng_notation
  27 from gnuradio.eng_option import eng_option
  28 from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, ra_waterfallsink, form, slider
  29 from optparse import OptionParser
  30 import wx
  31 import sys
  32 import Numeric
  33 import time
  34 import FFT
  35 import ephem
  36
  37 class continuum_calibration(gr.feval_dd):
  38     def eval(self, x):
  39         str = globals()["calibration_codelet"]
  40         exec(str)
  41         return(x)
  42
  43 class app_flow_graph(stdgui.gui_flow_graph):
  44     def __init__(self, frame, panel, vbox, argv):
  45         stdgui.gui_flow_graph.__init__(self)
  46
  47         self.frame = frame
  48         self.panel = panel
  49
  50         parser = OptionParser(option_class=eng_option)
  51         parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=(0, 0),
  52                           help="select USRP Rx side A or B (default=A)")
  53         parser.add_option("-d", "--decim", type="int", default=16,
  54                           help="set fgpa decimation rate to DECIM [default=%default]")
  55         parser.add_option("-f", "--freq", type="eng_float", default=None,
  56                           help="set frequency to FREQ", metavar="FREQ")
  57         parser.add_option("-a", "--avg", type="eng_float", default=1.0,
  58                 help="set spectral averaging alpha")
  59         parser.add_option("-i", "--integ", type="eng_float", default=1.0,
  60                 help="set integration time")
  61         parser.add_option("-g", "--gain", type="eng_float", default=None,
  62                           help="set gain in dB (default is midpoint)")
  63         parser.add_option("-l", "--reflevel", type="eng_float", default=30.0,
  64                           help="Set Total power reference level")
  65         parser.add_option("-y", "--division", type="eng_float", default=0.5,
  66                           help="Set Total power Y division size")
  67         parser.add_option("-e", "--longitude", type="eng_float", default=-76.02,                          help="Set Observer Longitude")
  68         parser.add_option("-c", "--latitude", type="eng_float", default=44.85,                          help="Set Observer Latitude")
  69         parser.add_option("-o", "--observing", type="eng_float", default=0.0,
  70                         help="Set observing frequency")
  71         parser.add_option("-x", "--ylabel", default="dB", help="Y axis label")
  72         parser.add_option("-z", "--divbase", type="eng_float", default=0.025, help="Y Division increment base")
  73         parser.add_option("-v", "--stripsize", type="eng_float", default=2400, help="Size of stripchart, in 2Hz samples")
  74         parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT")
  75
  76         parser.add_option("-N", "--decln", type="eng_float", default=999.99, help="Observing declination")
  77         parser.add_option("-X", "--prefix", default="./")
  78         parser.add_option("-M", "--fft_rate", type="eng_float", default=8.0, help="FFT Rate")
  79         parser.add_option("-A", "--calib_coeff", type="eng_float", default=1.0, help="Calibration coefficient")
  80         parser.add_option("-B", "--calib_offset", type="eng_float", default=0.0, help="Calibration coefficient")
  81         parser.add_option("-Q", "--calib_eqn", default="x = x * 1.0", help="Calibration equation")
  82         parser.add_option("-W", "--waterfall", action="store_true", default=False, help="Use Waterfall FFT display")
  83         parser.add_option("-S", "--setimode", action="store_true", default=False, help="Enable SETI processing of spectral data")
  84         parser.add_option("-T", "--setitimer", type="eng_float", default=15.0, help="Timer for computing doppler chirp for SETI analysis")
  85         parser.add_option("-K", "--setik", type="eng_float", default=1.5, help="K value for SETI analysis")
  86         (options, args) = parser.parse_args()
  87         if len(args) != 0:
  88             parser.print_help()
  89             sys.exit(1)
  90
  91         self.show_debug_info = True
  92         self.waterfall = options.waterfall
  93         self.setimode = options.setimode
  94         self.seticounter = 0
  95         self.setitimer = int(options.setitimer)
  96         self.setik = options.setik
  97         self.hitcounter = 0
  98         self.CHIRP_LOWER = 15
  99         self.CHIRP_UPPER = 50
 100
 101         # Calibration coefficient and offset
 102         self.calib_coeff = options.calib_coeff
 103         self.calib_offset = options.calib_offset
 104
 105         self.calib_eqn = options.calib_eqn
 106         globals()["calibration_codelet"] = self.calib_eqn
 107
 108         self.integ = options.integ
 109         self.avg_alpha = options.avg
 110         self.gain = options.gain
 111         self.decln = options.decln
 112
 113         # Set initial values for datalogging timed-output
 114         self.continuum_then = time.time()
 115         self.spectral_then = time.time()
 116
 117         # build the graph
 118
 119         self.u = usrp.source_c(decim_rate=options.decim)
 120         self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
 121         self.cardtype = self.u.daughterboard_id(0)
 122         # Set initial declination
 123         self.decln = options.decln
 124
 125         # determine the daughterboard subdevice we're using
 126         self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
 127
 128         input_rate = self.u.adc_freq() / self.u.decim_rate()
 129
 130         #
 131         # Set prefix for data files
 132         #
 133         self.prefix = options.prefix
 134
 135         #
 136         # The lower this number, the fewer sample frames are dropped
 137         #  in computing the FFT.  A sampled approach is taken to
 138         #  computing the FFT of the incoming data, which reduces
 139         #  sensitivity.  Increasing sensitivity inreases CPU loading.
 140         #
 141         self.fft_rate = options.fft_rate
 142
 143         self.fft_size = options.fft_size
 144
 145         # This buffer is used to remember the most-recent FFT display
 146         #   values.  Used later by self.write_spectral_data() to write
 147         #   spectral data to datalogging files.
 148         self.fft_outbuf = Numeric.zeros(options.fft_size, Numeric.Float64)
 149         self.old_hits = Numeric.zeros(10, Numeric.Float64)
 150         self.older_hits = Numeric.zeros(10, Numeric.Float64)
 151
 152         # Set up FFT display
 153         if self.waterfall == False:
 154            self.scope = ra_fftsink.ra_fft_sink_c (self, panel,
 155                fft_size=int(self.fft_size), sample_rate=input_rate,
 156                fft_rate=int(self.fft_rate), title="Spectral",
 157                ofunc=self.fft_outfunc, xydfunc=self.xydfunc)
 158         else:
 159             self.scope = ra_waterfallsink.ra_waterfallsink_c (self, panel,
 160                 fft_size=int(self.fft_size), sample_rate=input_rate,
 161                 fft_rate=int(self.fft_rate), title="Spectral", ofunc=self.fft_outfunc, xydfunc=self.xydfunc)
 162
 163         # Set up ephemeris data
 164         self.locality = ephem.Observer()
 165         self.locality.long = str(options.longitude)
 166         self.locality.lat = str(options.latitude)
 167
 168         # Set up stripchart display
 169         self.stripsize = int(options.stripsize)
 170         self.chart = ra_stripchartsink.stripchart_sink_f (self, panel,
 171             stripsize=self.stripsize,
 172             title="Continuum",
 173             xlabel="LMST Offset (Seconds)",
 174             scaling=1.0, ylabel=options.ylabel,
 175             divbase=options.divbase)
 176
 177         # Set center frequency
 178         self.centerfreq = options.freq
 179
 180         # Set observing frequency (might be different from actual programmed
 181         #    RF frequency)
 182         if options.observing == 0.0:
 183             self.observing = options.freq
 184         else:
 185             self.observing = options.observing
 186
 187         self.bw = input_rate
 188
 189         # We setup the first two integrators to produce a fixed integration
 190         # Down to 1Hz, with output at 1 samples/sec
 191         N = input_rate/5000
 192
 193         # Second stage runs on decimated output of first
 194         M = (input_rate/N)
 195
 196         # Create taps for first integrator
 197         t = range(0,N-1)
 198         tapsN = []
 199         for i in t:
 200              tapsN.append(1.0/N)
 201
 202         # Create taps for second integrator
 203         t = range(0,M-1)
 204         tapsM = []
 205         for i in t:
 206             tapsM.append(1.0/M)
 207
 208         #
 209         # The 3rd integrator is variable, and user selectable at runtime
 210         # This integrator doesn't decimate, but is used to set the
 211         #  final integration time based on the constant 1Hz input samples
 212         # The strip chart is fed at a constant 1Hz rate as a result
 213         #
 214
 215         #
 216         # Call constructors for receive chains
 217         #
 218
 219         # The three integrators--two FIR filters, and an IIR final filter
 220         self.integrator1 = gr.fir_filter_fff (N, tapsN)
 221         self.integrator2 = gr.fir_filter_fff (M, tapsM)
 222         self.integrator3 = gr.single_pole_iir_filter_ff(1.0)
 223
 224         # Split complex USRP stream into a pair of floats
 225         self.splitter = gr.complex_to_float (1);
 226
 227         # I squarer (detector)
 228         self.multI = gr.multiply_ff();
 229
 230         # Q squarer (detector)
 231         self.multQ = gr.multiply_ff();
 232
 233         # Adding squared I and Q to produce instantaneous signal power
 234         self.adder = gr.add_ff();
 235
 236         # Signal probe
 237         self.probe = gr.probe_signal_f();
 238
 239         #
 240         # Continuum calibration stuff
 241         #
 242         self.cal_mult = gr.multiply_const_ff(self.calib_coeff);
 243         self.cal_offs = gr.add_const_ff(self.calib_offset);
 244
 245         #self.cal_eqn = continuum_calibration();
 246
 247         #
 248         # Start connecting configured modules in the receive chain
 249         #
 250         self.connect(self.u, self.scope)
 251         self.connect(self.u, self.splitter)
 252
 253         # Connect splitter outputs to multipliers
 254         # First do I^2
 255         self.connect((self.splitter, 0), (self.multI,0))
 256         self.connect((self.splitter, 0), (self.multI,1))
 257
 258         # Then do Q^2
 259         self.connect((self.splitter, 1), (self.multQ,0))
 260         self.connect((self.splitter, 1), (self.multQ,1))
 261
 262         # Then sum the squares
 263         self.connect(self.multI, (self.adder,0))
 264         self.connect(self.multQ, (self.adder,1))
 265
 266         # Connect adder output to two-stages of FIR integrator
 267         #   followed by a single stage IIR integrator, and
 268         #   the calibrator
 269         self.connect(self.adder, self.integrator1,
 270            self.integrator2, self.integrator3, self.cal_mult,
 271            self.cal_offs, self.chart)
 272
 273         # Connect calibrator to probe
 274         # SPECIAL NOTE:  I'm setting the ground work here
 275         #   for completely changing the way local_calibrator
 276         #   works, including removing some horrible kludges for
 277         #   recording data.
 278         # But for now, self.probe() will be used to display the
 279         #  current instantaneous integrated detector value
 280         self.connect(self.cal_offs, self.probe)
 281
 282         self._build_gui(vbox)
 283
 284         # Make GUI agree with command-line
 285         self.integ = options.integ
 286         self.myform['integration'].set_value(int(options.integ))
 287         self.myform['average'].set_value(int(options.avg))
 288
 289         # Make integrator agree with command line
 290         self.set_integration(int(options.integ))
 291
 292         self.avg_alpha = options.avg
 293
 294         # Make spectral averager agree with command line
 295         if options.avg != 1.0:
 296             self.scope.set_avg_alpha(float(1.0/options.avg))
 297             self.scope.set_average(True)
 298
 299         # Set division size
 300         self.chart.set_y_per_div(options.division)
 301
 302         # Set reference(MAX) level
 303         self.chart.set_ref_level(options.reflevel)
 304
 305         # set initial values
 306
 307         if options.gain is None:
 308             # if no gain was specified, use the mid-point in dB
 309             g = self.subdev.gain_range()
 310             options.gain = float(g[0]+g[1])/2
 311
 312         if options.freq is None:
 313             # if no freq was specified, use the mid-point
 314             r = self.subdev.freq_range()
 315             options.freq = float(r[0]+r[1])/2
 316
 317         # Set the initial gain control
 318         self.set_gain(options.gain)
 319
 320         if not(self.set_freq(options.freq)):
 321             self._set_status_msg("Failed to set initial frequency")
 322
 323         # Set declination
 324         self.set_decln (self.decln)
 325
 326
 327         # RF hardware information
 328         self.myform['decim'].set_value(self.u.decim_rate())
 329         self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
 330         self.myform['dbname'].set_value(self.subdev.name())
 331
 332         # Set analog baseband filtering, if DBS_RX
 333         if self.cardtype == usrp_dbid.DBS_RX:
 334             lbw = (self.u.adc_freq() / self.u.decim_rate()) / 2
 335             if lbw < 1.0e6:
 336                 lbw = 1.0e6
 337             self.subdev.set_bw(lbw)
 338
 339         # Start the timer for the LMST display and datalogging
 340         self.lmst_timer.Start(1000)
 341
 342
 343     def _set_status_msg(self, msg):
 344         self.frame.GetStatusBar().SetStatusText(msg, 0)
 345
 346     def _build_gui(self, vbox):
 347
 348         def _form_set_freq(kv):
 349             return self.set_freq(kv['freq'])
 350
 351         def _form_set_decln(kv):
 352             return self.set_decln(kv['decln'])
 353
 354         # Position the FFT display
 355         vbox.Add(self.scope.win, 15, wx.EXPAND)
 356
 357         # Position the Total-power stripchart
 358         vbox.Add(self.chart.win, 15, wx.EXPAND)
 359
 360         # add control area at the bottom
 361         self.myform = myform = form.form()
 362         hbox = wx.BoxSizer(wx.HORIZONTAL)
 363         hbox.Add((7,0), 0, wx.EXPAND)
 364         vbox1 = wx.BoxSizer(wx.VERTICAL)
 365         myform['freq'] = form.float_field(
 366             parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
 367             callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
 368
 369         vbox1.Add((4,0), 0, 0)
 370
 371         myform['lmst_high'] = form.static_text_field(
 372             parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
 373         vbox1.Add((4,0), 0, 0)
 374
 375         myform['spec_data'] = form.static_text_field(
 376             parent=self.panel, sizer=vbox1, label="Spectral Cursor", weight=1)
 377         vbox1.Add((4,0), 0, 0)
 378
 379         vbox2 = wx.BoxSizer(wx.VERTICAL)
 380         g = self.subdev.gain_range()
 381         myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
 382                                            weight=1,
 383                                            min=int(g[0]), max=int(g[1]),
 384                                            callback=self.set_gain)
 385
 386         vbox2.Add((4,0), 0, 0)
 387         myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
 388                     label="Spectral Averaging (FFT frames)", weight=1, min=1, max=2000, callback=self.set_averaging)
 389
 390         vbox2.Add((4,0), 0, 0)
 391
 392         myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2,
 393                label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration)
 394
 395         vbox2.Add((4,0), 0, 0)
 396         myform['decln'] = form.float_field(
 397             parent=self.panel, sizer=vbox2, label="Current Declination", weight=1,
 398             callback=myform.check_input_and_call(_form_set_decln))
 399         vbox2.Add((4,0), 0, 0)
 400
 401         buttonbox = wx.BoxSizer(wx.HORIZONTAL)
 402         vbox.Add(buttonbox, 0, wx.CENTER)
 403         hbox.Add(vbox1, 0, 0)
 404         hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
 405         vbox.Add(hbox, 0, wx.EXPAND)
 406
 407         self._build_subpanel(vbox)
 408
 409         self.lmst_timer = wx.PyTimer(self.lmst_timeout)
 410         self.lmst_timeout()
 411
 412
 413     def _build_subpanel(self, vbox_arg):
 414         # build a secondary information panel (sometimes hidden)
 415
 416         # FIXME figure out how to have this be a subpanel that is always
 417         # created, but has its visibility controlled by foo.Show(True/False)
 418
 419         if not(self.show_debug_info):
 420             return
 421
 422         panel = self.panel
 423         vbox = vbox_arg
 424         myform = self.myform
 425
 426         #panel = wx.Panel(self.panel, -1)
 427         #vbox = wx.BoxSizer(wx.VERTICAL)
 428
 429         hbox = wx.BoxSizer(wx.HORIZONTAL)
 430         hbox.Add((5,0), 0)
 431         myform['decim'] = form.static_float_field(
 432             parent=panel, sizer=hbox, label="Decim")
 433
 434         hbox.Add((5,0), 1)
 435         myform['fs@usb'] = form.static_float_field(
 436             parent=panel, sizer=hbox, label="Fs@USB")
 437
 438         hbox.Add((5,0), 1)
 439         myform['dbname'] = form.static_text_field(
 440             parent=panel, sizer=hbox)
 441
 442         hbox.Add((5,0), 1)
 443         myform['baseband'] = form.static_float_field(
 444             parent=panel, sizer=hbox, label="Analog BB")
 445
 446         hbox.Add((5,0), 1)
 447         myform['ddc'] = form.static_float_field(
 448             parent=panel, sizer=hbox, label="DDC")
 449
 450         hbox.Add((5,0), 0)
 451         vbox.Add(hbox, 0, wx.EXPAND)
 452
 453
 454
 455     def set_freq(self, target_freq):
 456         """
 457         Set the center frequency we're interested in.
 458
 459         @param target_freq: frequency in Hz
 460         @rypte: bool
 461
 462         Tuning is a two step process.  First we ask the front-end to
 463         tune as close to the desired frequency as it can.  Then we use
 464         the result of that operation and our target_frequency to
 465         determine the value for the digital down converter.
 466         """
 467         #
 468         # Everything except BASIC_RX should support usrp.tune()
 469         #
 470         if not (self.cardtype == usrp_dbid.BASIC_RX):
 471             r = usrp.tune(self.u, 0, self.subdev, target_freq)
 472         else:
 473             r = self.u.set_rx_freq(0, target_freq)
 474             f = self.u.rx_freq(0)
 475             if abs(f-target_freq) > 2.0e3:
 476                 r = 0
 477         if r:
 478             self.myform['freq'].set_value(target_freq)     # update displayed value
 479             #
 480             # Make sure calibrator knows our target freq
 481             #
 482
 483             # Remember centerfreq---used for doppler calcs
 484             delta = self.centerfreq - target_freq
 485             self.centerfreq = target_freq
 486             self.observing -= delta
 487             self.scope.set_baseband_freq (self.observing)
 488
 489             self.myform['baseband'].set_value(r.baseband_freq)
 490             self.myform['ddc'].set_value(r.dxc_freq)
 491
 492             return True
 493
 494         return False
 495
 496     def set_decln(self, dec):
 497         self.decln = dec
 498         self.myform['decln'].set_value(dec)     # update displayed value
 499
 500     def set_gain(self, gain):
 501         self.myform['gain'].set_value(gain)     # update displayed value
 502         self.subdev.set_gain(gain)
 503         self.gain = gain
 504
 505     def set_averaging(self, avval):
 506         self.myform['average'].set_value(avval)
 507         self.scope.set_avg_alpha(1.0/(avval))
 508         self.scope.set_average(True)
 509         self.avg_alpha = avval
 510
 511     def set_integration(self, integval):
 512         self.integrator3.set_taps(1.0/integval)
 513         self.myform['integration'].set_value(integval)
 514         self.integ = integval
 515
 516     #
 517     # Timeout function
 518     # Used to update LMST display, as well as current
 519     #  continuum value
 520     #
 521     # We also write external data-logging files here
 522     #
 523     def lmst_timeout(self):
 524          self.locality.date = ephem.now()
 525          x = self.probe.level()
 526          sidtime = self.locality.sidereal_time()
 527          # LMST
 528          s = str(ephem.hours(sidtime))
 529          # Continuum detector value
 530          sx = "%7.4f" % x
 531          s = s + "\nDet: " + str(sx)
 532          sx = "%2d" % self.hitcounter
 533          sy = "%2d" % self.CHIRP_LOWER
 534          s = s + "\nH: " + str(sx) + " Cl: " + str(sy)
 535          self.myform['lmst_high'].set_value(s)
 536
 537          #
 538          # Write data out to recording files
 539          #
 540          self.write_continuum_data(x,sidtime)
 541          self.write_spectral_data(self.fft_outbuf,sidtime)
 542
 543          if self.setimode == True:
 544              self.seti_analysis(self.fft_outbuf,sidtime)
 545
 546     def fft_outfunc(self,data,l):
 547         self.fft_outbuf=data
 548
 549     def write_continuum_data(self,data,sidtime):
 550
 551         # Create localtime structure for producing filename
 552         foo = time.localtime()
 553         pfx = self.prefix
 554         filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 555            foo.tm_mon, foo.tm_mday, foo.tm_hour)
 556
 557         # Open the data file, appending
 558         continuum_file = open (filenamestr+".tpdat","a")
 559
 560         flt = "%6.3f" % data
 561         inter = self.decln
 562         integ = self.integ
 563         fc = self.observing
 564         fc = fc / 1000000
 565         bw = self.bw
 566         bw = bw / 1000000
 567         ga = self.gain
 568
 569         now = time.time()
 570
 571         #
 572         # If time to write full header info (saves storage this way)
 573         #
 574         if (now - self.continuum_then > 20):
 575             self.continuum_then = now
 576
 577             continuum_file.write(str(ephem.hours(sidtime))+" "+flt+" Dn="+str(inter)+",")
 578             continuum_file.write("Ti="+str(integ)+",Fc="+str(fc)+",Bw="+str(bw))
 579             continuum_file.write(",Ga="+str(ga)+"\n")
 580         else:
 581             continuum_file.write(str(ephem.hours(sidtime))+" "+flt+"\n")
 582
 583         continuum_file.close()
 584         return(data)
 585
 586     def write_spectral_data(self,data,sidtime):
 587
 588         now = time.time()
 589         delta = 10
 590
 591         # If time to write out spectral data
 592         # We don't write this out every time, in order to
 593         #   save disk space.  Since the spectral data are
 594         #   typically heavily averaged, writing this data
 595         #   "once in a while" is OK.
 596         #
 597         if (now - self.spectral_then >= delta):
 598             self.spectral_then = now
 599
 600             # Get localtime structure to make filename from
 601             foo = time.localtime()
 602
 603             pfx = self.prefix
 604             filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 605                foo.tm_mon, foo.tm_mday, foo.tm_hour)
 606
 607             # Open the file
 608             spectral_file = open (filenamestr+".sdat","a")
 609
 610             # Setup data fields to be written
 611             r = data
 612             inter = self.decln
 613             fc = self.observing
 614             fc = fc / 1000000
 615             bw = self.bw
 616             bw = bw / 1000000
 617             av = self.avg_alpha
 618
 619             # Write those fields
 620             spectral_file.write("data:"+str(ephem.hours(sidtime))+" Dn="+str(inter)+",Fc="+str(fc)+",Bw="+str(bw)+",Av="+str(av))
 621             spectral_file.write(" "+str(r)+"\n")
 622             spectral_file.close()
 623             return(data)
 624
 625         return(data)
 626
 627     def seti_analysis(self,fftbuf,sidtime):
 628         l = len(fftbuf)
 629         x = 0
 630         hits = []
 631         if self.seticounter < self.setitimer:
 632             self.seticounter = self.seticounter + 1
 633             return
 634         else:
 635             self.seticounter = 0
 636
 637         # Run through FFT output buffer, computing standard deviation (Sigma)
 638         avg = 0
 639         # First compute average
 640         for i in range(0,l):
 641             avg = avg + fftbuf[i]
 642         avg = avg / l
 643
 644         sigma = 0.0
 645         # Then compute standard deviation (Sigma)
 646         for i in range(0,l):
 647             d = fftbuf[i] - avg
 648             sigma = sigma + (d*d)
 649
 650         sigma = Numeric.sqrt(sigma/l)
 651
 652         #
 653         # Snarfle through the FFT output buffer again, looking for
 654         #    outlying data points
 655
 656         start_f = self.observing - (self.bw/2)
 657         current_f = start_f
 658         f_incr = self.bw / l
 659         l = len(fftbuf)
 660         hit = -1
 661
 662         # -nyquist to DC
 663         for i in range(l/2,l):
 664             #
 665             # If current FFT buffer has an item that exceeds the specified
 666             #  sigma
 667             #
 668             if ((fftbuf[i] - avg) > (self.setik * sigma)):
 669                 hits.append(current_f)
 670             current_f = current_f + f_incr
 671
 672         # DC to nyquist
 673         for i in range(0,l/2):
 674             #
 675             # If current FFT buffer has an item that exceeds the specified
 676             #  sigma
 677             #
 678             if ((fftbuf[i] - avg) > (self.setik * sigma)):
 679                 hits.append(current_f)
 680             current_f = current_f + f_incr
 681
 682         if (len(hits) <= 0):
 683             return
 684
 685         if (len(hits) > len(self.old_hits)*2):
 686             return
 687
 688         #
 689         #
 690         # Calculate chirp limits from first principles
 691         #
 692         #
 693         earth_diam = 12500
 694         earth_circ = earth_diam * 3.14159
 695         surface_speed = earth_circ / 24
 696         surface_speed = surface_speed / 3600
 697
 698         c1 = (surface_speed/2) / 299792.0
 699         c2 = (surface_speed*5) / 299792.0
 700
 701         self.CHIRP_LOWER = (c1 * self.observing) * self.setitimer
 702         self.CHIRP_UPPER = (c2 * self.observing) * self.setitimer
 703
 704         #
 705         # Run through all three hit buffers, computing difference between
 706         #   frequencies found there, if they're all within the chirp limits
 707         #   declare a good hit
 708         #
 709         good_hit = 0
 710         for i in range(0,min(len(hits),len(self.old_hits))):
 711             f_diff = abs(self.old_hits[i] - hits[i])
 712             f_diff2 = abs(self.older_hits[i] - self.old_hits[i])
 713             # If frequency difference is within range, we have a hit
 714             if (f_diff >= self.CHIRP_LOWER and f_diff <= self.CHIRP_UPPER):
 715                 if (f_diff2 >= self.CHIRP_LOWER and f_diff <= self.CHIRP_UPPER):
 716                     good_hit = 1
 717                     self.hitcounter = self.hitcounter + 1
 718                     break
 719
 720         if (good_hit != 0):
 721             self.write_hits(hits,sidtime)
 722
 723         # Save old hits
 724         for i in range(0,len(self.older_hits)):
 725             self.older_hits[i] = self.old_hits[i]
 726             self.old_hits[i] = 0
 727         for i in range(0,min(len(hits),len(self.old_hits))):
 728             self.old_hits[i] = hits[i]
 729
 730         return
 731
 732     def write_hits(self,hits,sidtime):
 733         # Create localtime structure for producing filename
 734         foo = time.localtime()
 735         pfx = self.prefix
 736         filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 737            foo.tm_mon, foo.tm_mday, foo.tm_hour)
 738
 739         # Open the data file, appending
 740         hits_file = open (filenamestr+".seti","a")
 741         hits_file.write(str(ephem.hours(sidtime))+" "+str(hits)+"\n")
 742         hits_file.close()
 743         return
 744
 745     def xydfunc(self,xyv):
 746         magn = int(Numeric.log10(self.observing))
 747         if (magn == 6 or magn == 7 or magn == 8):
 748             magn = 6
 749         dfreq = xyv[0] * pow(10.0,magn)
 750         ratio = self.observing / dfreq
 751         vs = 1.0 - ratio
 752         vs *= 299792.0
 753         if magn >= 9:
 754            xhz = "Ghz"
 755         elif magn >= 6:
 756            xhz = "Mhz"
 757         elif magn <= 5:
 758            xhz =  "Khz"
 759         s = "%.6f%s\n%.3fdB" % (xyv[0], xhz, xyv[1])
 760         s2 = "\n%.3fkm/s" % vs
 761         self.myform['spec_data'].set_value(s+s2)
 762
 763     def toggle_cal(self):
 764         if (self.calstate == True):
 765           self.calstate = False
 766           self.u.write_io(0,0,(1<<15))
 767           self.calibrator.SetLabel("Calibration Source: Off")
 768         else:
 769           self.calstate = True
 770           self.u.write_io(0,(1<<15),(1<<15))
 771           self.calibrator.SetLabel("Calibration Source: On")
 772
 773     def toggle_annotation(self):
 774         if (self.annotate_state == True):
 775           self.annotate_state = False
 776           self.annotation.SetLabel("Annotation: Off")
 777         else:
 778           self.annotate_state = True
 779           self.annotation.SetLabel("Annotation: On")
 780
 781
 782 def main ():
 783     app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY SPECTRAL/CONTINUUM RECEIVER: $Revision$", nstatus=1)
 784     app.MainLoop()
 785
 786 if __name__ == '__main__':
 787     main ()