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 3, or (at your option)
  10 # any later version.
  11 #
  12 # GNU Radio is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with GNU Radio; see the file COPYING.  If not, write to
  19 # the Free Software Foundation, Inc., 51 Franklin Street,
  20 # Boston, MA 02110-1301, USA.
  21 #
  22
  23 from gnuradio import gr, gru
  24 from gnuradio import usrp
  25 from usrpm 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, waterfallsink
  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("-W", "--waterfall", action="store_true", default=False, help="Use Waterfall FFT display")
  82         parser.add_option("-S", "--setimode", action="store_true", default=False, help="Enable SETI processing of spectral data")
  83         parser.add_option("-K", "--setik", type="eng_float", default=1.5, help="K value for SETI analysis")
  84         parser.add_option("-T", "--setibandwidth", type="eng_float", default=12500, help="Instantaneous SETI observing bandwidth--must be divisor of 250Khz")
  85         parser.add_option("-n", "--notches", action="store_true", default=False,
  86             help="Notches appear after all other arguments")
  87         parser.add_option("-Q", "--seti_range", type="eng_float", default=1.0e6, help="Total scan width, in Hz for SETI scans")
  88         (options, args) = parser.parse_args()
  89
  90         self.notches = Numeric.zeros(64,Numeric.Float64)
  91         if len(args) != 0 and options.notches == False:
  92             parser.print_help()
  93             sys.exit(1)
  94
  95         if len(args) == 0 and options.notches != False:
  96             parser.print_help()
  97             sys.exit()
  98
  99         self.use_notches = options.notches
 100
 101         # Get notch locations
 102         j = 0
 103         for i in args:
 104             self.notches[j] = float(i)
 105             j = j+1
 106
 107         self.notch_count = j
 108
 109         self.show_debug_info = True
 110
 111         # Pick up waterfall option
 112         self.waterfall = options.waterfall
 113
 114         # SETI mode stuff
 115         self.setimode = options.setimode
 116         self.seticounter = 0
 117         self.setik = options.setik
 118         self.seti_fft_bandwidth = int(options.setibandwidth)
 119
 120         # Calculate binwidth
 121         binwidth = self.seti_fft_bandwidth / options.fft_size
 122
 123         # Use binwidth, and knowledge of likely chirp rates to set reasonable
 124         #  values for SETI analysis code.   We assume that SETI signals will
 125         #  chirp at somewhere between 0.10Hz/sec and 0.25Hz/sec.
 126         #
 127         # upper_limit is the "worst case"--that is, the case for which we have
 128         #  to wait the longest to actually see any drift, due to the quantizing
 129         #  on FFT bins.
 130         upper_limit = binwidth / 0.10
 131         self.setitimer = int(upper_limit * 2.00)
 132         self.scanning = True
 133
 134         # Calculate the CHIRP values based on Hz/sec
 135         self.CHIRP_LOWER = 0.10 * self.setitimer
 136         self.CHIRP_UPPER = 0.25 * self.setitimer
 137
 138         # Reset hit counters to 0
 139         self.hitcounter = 0
 140         self.s1hitcounter = 0
 141         self.s2hitcounter = 0
 142         self.avgdelta = 0
 143         # We scan through 2Mhz of bandwidth around the chosen center freq
 144         self.seti_freq_range = options.seti_range
 145         # Calculate lower edge
 146         self.setifreq_lower = options.freq - (self.seti_freq_range/2)
 147         self.setifreq_current = options.freq
 148         # Calculate upper edge
 149         self.setifreq_upper = options.freq + (self.seti_freq_range/2)
 150
 151         # Maximum "hits" in a line
 152         self.nhits = 20
 153
 154         # Number of lines for analysis
 155         self.nhitlines = 4
 156
 157         # We change center frequencies based on nhitlines and setitimer
 158         self.setifreq_timer = self.setitimer * (self.nhitlines * 5)
 159         print "setifreq_timer", self.setifreq_timer
 160
 161         # Create actual timer
 162         self.seti_then = time.time()
 163
 164         # The hits recording array
 165         self.hits_array = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64)
 166         self.hit_intensities = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64)
 167         # Calibration coefficient and offset
 168         self.calib_coeff = options.calib_coeff
 169         self.calib_offset = options.calib_offset
 170         if self.calib_offset < -750:
 171             self.calib_offset = -750
 172         if self.calib_offset > 750:
 173             self.calib_offset = 750
 174
 175         if self.calib_coeff < 1:
 176             self.calib_coeff = 1
 177         if self.calib_coeff > 100:
 178             self.calib_coeff = 100
 179
 180         self.integ = options.integ
 181         self.avg_alpha = options.avg
 182         self.gain = options.gain
 183         self.decln = options.decln
 184
 185         # Set initial values for datalogging timed-output
 186         self.continuum_then = time.time()
 187         self.spectral_then = time.time()
 188
 189         # build the graph
 190
 191         #
 192         # If SETI mode, we always run at maximum USRP decimation
 193         #
 194         if (self.setimode):
 195             options.decim = 256
 196
 197         self.u = usrp.source_c(decim_rate=options.decim)
 198         self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
 199         # Set initial declination
 200         self.decln = options.decln
 201
 202         # determine the daughterboard subdevice we're using
 203         self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec)
 204         self.cardtype = self.subdev.dbid()
 205
 206         input_rate = self.u.adc_freq() / self.u.decim_rate()
 207
 208         #
 209         # Set prefix for data files
 210         #
 211         self.prefix = options.prefix
 212
 213         #
 214         # The lower this number, the fewer sample frames are dropped
 215         #  in computing the FFT.  A sampled approach is taken to
 216         #  computing the FFT of the incoming data, which reduces
 217         #  sensitivity.  Increasing sensitivity inreases CPU loading.
 218         #
 219         self.fft_rate = options.fft_rate
 220
 221         self.fft_size = int(options.fft_size)
 222
 223         # This buffer is used to remember the most-recent FFT display
 224         #   values.  Used later by self.write_spectral_data() to write
 225         #   spectral data to datalogging files, and by the SETI analysis
 226         #   function.
 227         #
 228         self.fft_outbuf = Numeric.zeros(self.fft_size, Numeric.Float64)
 229
 230         #
 231         # If SETI mode, only look at seti_fft_bandwidth
 232         #   at a time.
 233         #
 234         if (self.setimode):
 235             self.fft_input_rate = self.seti_fft_bandwidth
 236
 237             #
 238             # Build a decimating bandpass filter
 239             #
 240             self.fft_input_taps = gr.firdes.complex_band_pass (1.0,
 241                input_rate,
 242                -(int(self.fft_input_rate/2)), int(self.fft_input_rate/2), 200,
 243                gr.firdes.WIN_HAMMING, 0)
 244
 245             #
 246             # Compute required decimation factor
 247             #
 248             decimation = int(input_rate/self.fft_input_rate)
 249             self.fft_bandpass = gr.fir_filter_ccc (decimation,
 250                 self.fft_input_taps)
 251         else:
 252             self.fft_input_rate = input_rate
 253
 254         # Set up FFT display
 255         if self.waterfall == False:
 256            self.scope = ra_fftsink.ra_fft_sink_c (self, panel,
 257                fft_size=int(self.fft_size), sample_rate=self.fft_input_rate,
 258                fft_rate=int(self.fft_rate), title="Spectral",
 259                ofunc=self.fft_outfunc, xydfunc=self.xydfunc)
 260         else:
 261             self.scope = ra_waterfallsink.waterfall_sink_c (self, panel,
 262                 fft_size=int(self.fft_size), sample_rate=self.fft_input_rate,
 263                 fft_rate=int(self.fft_rate), title="Spectral", ofunc=self.fft_outfunc, size=(1100, 600), xydfunc=self.xydfunc, ref_level=0, span=10)
 264
 265         # Set up ephemeris data
 266         self.locality = ephem.Observer()
 267         self.locality.long = str(options.longitude)
 268         self.locality.lat = str(options.latitude)
 269         # We make notes about Sunset/Sunrise in Continuum log files
 270         self.sun = ephem.Sun()
 271         self.sunstate = "??"
 272
 273         # Set up stripchart display
 274         self.stripsize = int(options.stripsize)
 275         if self.setimode == False:
 276             self.chart = ra_stripchartsink.stripchart_sink_f (self, panel,
 277                 stripsize=self.stripsize,
 278                 title="Continuum",
 279                 xlabel="LMST Offset (Seconds)",
 280                 scaling=1.0, ylabel=options.ylabel,
 281                 divbase=options.divbase)
 282
 283         # Set center frequency
 284         self.centerfreq = options.freq
 285
 286         # Set observing frequency (might be different from actual programmed
 287         #    RF frequency)
 288         if options.observing == 0.0:
 289             self.observing = options.freq
 290         else:
 291             self.observing = options.observing
 292
 293         self.bw = input_rate
 294
 295         # We setup the first two integrators to produce a fixed integration
 296         # Down to 1Hz, with output at 1 samples/sec
 297         N = input_rate/5000
 298
 299         # Second stage runs on decimated output of first
 300         M = (input_rate/N)
 301
 302         # Create taps for first integrator
 303         t = range(0,N-1)
 304         tapsN = []
 305         for i in t:
 306              tapsN.append(1.0/N)
 307
 308         # Create taps for second integrator
 309         t = range(0,M-1)
 310         tapsM = []
 311         for i in t:
 312             tapsM.append(1.0/M)
 313
 314         #
 315         # The 3rd integrator is variable, and user selectable at runtime
 316         # This integrator doesn't decimate, but is used to set the
 317         #  final integration time based on the constant 1Hz input samples
 318         # The strip chart is fed at a constant 1Hz rate as a result
 319         #
 320
 321         #
 322         # Call constructors for receive chains
 323         #
 324
 325         if self.setimode == False:
 326             # The three integrators--two FIR filters, and an IIR final filter
 327             self.integrator1 = gr.fir_filter_fff (N, tapsN)
 328             self.integrator2 = gr.fir_filter_fff (M, tapsM)
 329             self.integrator3 = gr.single_pole_iir_filter_ff(1.0)
 330
 331             # The detector
 332             self.detector = gr.complex_to_mag_squared()
 333
 334
 335             # Signal probe
 336             self.probe = gr.probe_signal_f();
 337
 338             #
 339             # Continuum calibration stuff
 340             #
 341             x = self.calib_coeff/100.0
 342             self.cal_mult = gr.multiply_const_ff(self.calib_coeff/100.0);
 343             self.cal_offs = gr.add_const_ff(self.calib_offset*(x*8000));
 344
 345             if self.use_notches == True:
 346                 self.compute_notch_taps(self.notches)
 347                 self.notch_filt = gr.fft_filter_ccc(1, self.notch_taps)
 348
 349         #
 350         # Start connecting configured modules in the receive chain
 351         #
 352
 353         # The scope--handle SETI mode
 354         if (self.setimode == False):
 355             if (self.use_notches == True):
 356                 self.connect(self.u, self.notch_filt, self.scope)
 357             else:
 358                 self.connect(self.u, self.scope)
 359         else:
 360             if (self.use_notches == True):
 361                 self.connect(self.u, self.notch_filt,
 362                     self.fft_bandpass, self.scope)
 363             else:
 364                 self.connect(self.u, self.fft_bandpass, self.scope)
 365
 366         if self.setimode == False:
 367             if (self.use_notches == True):
 368                 self.connect(self.notch_filt, self.detector,
 369                     self.integrator1, self.integrator2,
 370                     self.integrator3, self.cal_mult, self.cal_offs, self.chart)
 371             else:
 372                 self.connect(self.u, self.detector,
 373                     self.integrator1, self.integrator2,
 374                     self.integrator3, self.cal_mult, self.cal_offs, self.chart)
 375
 376             #  current instantaneous integrated detector value
 377             self.connect(self.cal_offs, self.probe)
 378
 379         self._build_gui(vbox)
 380
 381         # Make GUI agree with command-line
 382         self.integ = options.integ
 383         if self.setimode == False:
 384             self.myform['integration'].set_value(int(options.integ))
 385             self.myform['offset'].set_value(self.calib_offset)
 386             self.myform['dcgain'].set_value(self.calib_coeff)
 387         self.myform['average'].set_value(int(options.avg))
 388
 389
 390         if self.setimode == False:
 391             # Make integrator agree with command line
 392             self.set_integration(int(options.integ))
 393
 394         self.avg_alpha = options.avg
 395
 396         # Make spectral averager agree with command line
 397         if options.avg != 1.0:
 398             self.scope.set_avg_alpha(float(1.0/options.avg))
 399             self.scope.set_average(True)
 400
 401         if self.setimode == False:
 402             # Set division size
 403             self.chart.set_y_per_div(options.division)
 404             # Set reference(MAX) level
 405             self.chart.set_ref_level(options.reflevel)
 406
 407         # set initial values
 408
 409         if options.gain is None:
 410             # if no gain was specified, use the mid-point in dB
 411             g = self.subdev.gain_range()
 412             options.gain = float(g[0]+g[1])/2
 413
 414         if options.freq is None:
 415             # if no freq was specified, use the mid-point
 416             r = self.subdev.freq_range()
 417             options.freq = float(r[0]+r[1])/2
 418
 419         # Set the initial gain control
 420         self.set_gain(options.gain)
 421
 422         if not(self.set_freq(options.freq)):
 423             self._set_status_msg("Failed to set initial frequency")
 424
 425         # Set declination
 426         self.set_decln (self.decln)
 427
 428
 429         # RF hardware information
 430         self.myform['decim'].set_value(self.u.decim_rate())
 431         self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate())
 432         self.myform['dbname'].set_value(self.subdev.name())
 433
 434         # Set analog baseband filtering, if DBS_RX
 435         if self.cardtype in (usrp_dbid.DBS_RX, usrp_dbid.DBS_RX_REV_2_1):
 436             lbw = (self.u.adc_freq() / self.u.decim_rate()) / 2
 437             if lbw < 1.0e6:
 438                 lbw = 1.0e6
 439             self.subdev.set_bw(lbw)
 440
 441         # Start the timer for the LMST display and datalogging
 442         self.lmst_timer.Start(1000)
 443
 444
 445     def _set_status_msg(self, msg):
 446         self.frame.GetStatusBar().SetStatusText(msg, 0)
 447
 448     def _build_gui(self, vbox):
 449
 450         def _form_set_freq(kv):
 451             # Adjust current SETI frequency, and limits
 452             self.setifreq_lower = kv['freq'] - (self.seti_freq_range/2)
 453             self.setifreq_current = kv['freq']
 454             self.setifreq_upper = kv['freq'] + (self.seti_freq_range/2)
 455
 456             # Reset SETI analysis timer
 457             self.seti_then = time.time()
 458             # Zero-out hits array when changing frequency
 459             self.hits_array[:,:] = 0.0
 460             self.hit_intensities[:,:] = -60.0
 461
 462             return self.set_freq(kv['freq'])
 463
 464         def _form_set_decln(kv):
 465             return self.set_decln(kv['decln'])
 466
 467         # Position the FFT display
 468         vbox.Add(self.scope.win, 15, wx.EXPAND)
 469
 470         if self.setimode == False:
 471             # Position the Total-power stripchart
 472             vbox.Add(self.chart.win, 15, wx.EXPAND)
 473
 474         # add control area at the bottom
 475         self.myform = myform = form.form()
 476         hbox = wx.BoxSizer(wx.HORIZONTAL)
 477         hbox.Add((7,0), 0, wx.EXPAND)
 478         vbox1 = wx.BoxSizer(wx.VERTICAL)
 479         myform['freq'] = form.float_field(
 480             parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
 481             callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
 482
 483         vbox1.Add((4,0), 0, 0)
 484
 485         myform['lmst_high'] = form.static_text_field(
 486             parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
 487         vbox1.Add((4,0), 0, 0)
 488
 489         if self.setimode == False:
 490             myform['spec_data'] = form.static_text_field(
 491                 parent=self.panel, sizer=vbox1, label="Spectral Cursor", weight=1)
 492             vbox1.Add((4,0), 0, 0)
 493
 494         vbox2 = wx.BoxSizer(wx.VERTICAL)
 495         if self.setimode == False:
 496             vbox3 = wx.BoxSizer(wx.VERTICAL)
 497         g = self.subdev.gain_range()
 498         myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
 499                                            weight=1,
 500                                            min=int(g[0]), max=int(g[1]),
 501                                            callback=self.set_gain)
 502
 503         vbox2.Add((4,0), 0, 0)
 504         if self.setimode == True:
 505             max_savg = 100
 506         else:
 507             max_savg = 3000
 508         myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
 509                     label="Spectral Averaging (FFT frames)", weight=1, min=1, max=max_savg, callback=self.set_averaging)
 510
 511         # Set up scan control button when in SETI mode
 512         if (self.setimode == True):
 513                 # SETI scanning control
 514                 buttonbox = wx.BoxSizer(wx.HORIZONTAL)
 515                 self.scan_control = form.button_with_callback(self.panel,
 516                       label="Scan: On ",
 517                       callback=self.toggle_scanning)
 518
 519                 buttonbox.Add(self.scan_control, 0, wx.CENTER)
 520                 vbox2.Add(buttonbox, 0, wx.CENTER)
 521
 522         vbox2.Add((4,0), 0, 0)
 523
 524         if self.setimode == False:
 525             myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2,
 526                    label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration)
 527
 528             vbox2.Add((4,0), 0, 0)
 529
 530         myform['decln'] = form.float_field(
 531             parent=self.panel, sizer=vbox2, label="Current Declination", weight=1,
 532             callback=myform.check_input_and_call(_form_set_decln))
 533         vbox2.Add((4,0), 0, 0)
 534
 535         if self.setimode == False:
 536             myform['offset'] = form.slider_field(parent=self.panel, sizer=vbox3,
 537                 label="Post-Detector Offset", weight=1, min=-750, max=750,
 538                 callback=self.set_pd_offset)
 539             vbox3.Add((2,0), 0, 0)
 540             myform['dcgain'] = form.slider_field(parent=self.panel, sizer=vbox3,
 541                 label="Post-Detector Gain", weight=1, min=1, max=100,
 542                 callback=self.set_pd_gain)
 543             vbox3.Add((2,0), 0, 0)
 544         hbox.Add(vbox1, 0, 0)
 545         hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
 546
 547         if self.setimode == False:
 548             hbox.Add(vbox3, wx.ALIGN_RIGHT, 0)
 549
 550         vbox.Add(hbox, 0, wx.EXPAND)
 551
 552         self._build_subpanel(vbox)
 553
 554         self.lmst_timer = wx.PyTimer(self.lmst_timeout)
 555         #self.lmst_timeout()
 556
 557
 558     def _build_subpanel(self, vbox_arg):
 559         # build a secondary information panel (sometimes hidden)
 560
 561         # FIXME figure out how to have this be a subpanel that is always
 562         # created, but has its visibility controlled by foo.Show(True/False)
 563
 564         if not(self.show_debug_info):
 565             return
 566
 567         panel = self.panel
 568         vbox = vbox_arg
 569         myform = self.myform
 570
 571         #panel = wx.Panel(self.panel, -1)
 572         #vbox = wx.BoxSizer(wx.VERTICAL)
 573
 574         hbox = wx.BoxSizer(wx.HORIZONTAL)
 575         hbox.Add((5,0), 0)
 576         myform['decim'] = form.static_float_field(
 577             parent=panel, sizer=hbox, label="Decim")
 578
 579         hbox.Add((5,0), 1)
 580         myform['fs@usb'] = form.static_float_field(
 581             parent=panel, sizer=hbox, label="Fs@USB")
 582
 583         hbox.Add((5,0), 1)
 584         myform['dbname'] = form.static_text_field(
 585             parent=panel, sizer=hbox)
 586
 587         hbox.Add((5,0), 1)
 588         myform['baseband'] = form.static_float_field(
 589             parent=panel, sizer=hbox, label="Analog BB")
 590
 591         hbox.Add((5,0), 1)
 592         myform['ddc'] = form.static_float_field(
 593             parent=panel, sizer=hbox, label="DDC")
 594
 595         hbox.Add((5,0), 0)
 596         vbox.Add(hbox, 0, wx.EXPAND)
 597
 598
 599
 600     def set_freq(self, target_freq):
 601         """
 602         Set the center frequency we're interested in.
 603
 604         @param target_freq: frequency in Hz
 605         @rypte: bool
 606
 607         Tuning is a two step process.  First we ask the front-end to
 608         tune as close to the desired frequency as it can.  Then we use
 609         the result of that operation and our target_frequency to
 610         determine the value for the digital down converter.
 611         """
 612         #
 613         # Everything except BASIC_RX should support usrp.tune()
 614         #
 615         if not (self.cardtype == usrp_dbid.BASIC_RX):
 616             r = usrp.tune(self.u, 0, self.subdev, target_freq)
 617         else:
 618             r = self.u.set_rx_freq(0, target_freq)
 619             f = self.u.rx_freq(0)
 620             if abs(f-target_freq) > 2.0e3:
 621                 r = 0
 622         if r:
 623             self.myform['freq'].set_value(target_freq)     # update displayed value
 624             #
 625             # Make sure calibrator knows our target freq
 626             #
 627
 628             # Remember centerfreq---used for doppler calcs
 629             delta = self.centerfreq - target_freq
 630             self.centerfreq = target_freq
 631             self.observing -= delta
 632             self.scope.set_baseband_freq (self.observing)
 633
 634             self.myform['baseband'].set_value(r.baseband_freq)
 635             self.myform['ddc'].set_value(r.dxc_freq)
 636
 637             if self.use_notches == True:
 638                 self.compute_notch_taps(self.notches)
 639                 self.notch_filt.set_taps(self.notch_taps)
 640
 641             return True
 642
 643         return False
 644
 645     def set_decln(self, dec):
 646         self.decln = dec
 647         self.myform['decln'].set_value(dec)     # update displayed value
 648
 649     def set_gain(self, gain):
 650         self.myform['gain'].set_value(gain)     # update displayed value
 651         self.subdev.set_gain(gain)
 652         self.gain = gain
 653
 654     def set_averaging(self, avval):
 655         self.myform['average'].set_value(avval)
 656         self.scope.set_avg_alpha(1.0/(avval))
 657         self.scope.set_average(True)
 658         self.avg_alpha = avval
 659
 660     def set_integration(self, integval):
 661         if self.setimode == False:
 662             self.integrator3.set_taps(1.0/integval)
 663         self.myform['integration'].set_value(integval)
 664         self.integ = integval
 665
 666     #
 667     # Timeout function
 668     # Used to update LMST display, as well as current
 669     #  continuum value
 670     #
 671     # We also write external data-logging files here
 672     #
 673     def lmst_timeout(self):
 674          self.locality.date = ephem.now()
 675          if self.setimode == False:
 676              x = self.probe.level()
 677          sidtime = self.locality.sidereal_time()
 678          # LMST
 679          s = str(ephem.hours(sidtime)) + " " + self.sunstate
 680          # Continuum detector value
 681          if self.setimode == False:
 682              sx = "%7.4f" % x
 683              s = s + "\nDet: " + str(sx)
 684          else:
 685              sx = "%2d" % self.hitcounter
 686              s1 = "%2d" % self.s1hitcounter
 687              s2 = "%2d" % self.s2hitcounter
 688              sa = "%4.2f" % self.avgdelta
 689              sy = "%3.1f-%3.1f" % (self.CHIRP_LOWER, self.CHIRP_UPPER)
 690              s = s + "\nHits: " + str(sx) + "\nS1:" + str(s1) + " S2:" + str(s2)
 691              s = s + "\nAv D: " + str(sa) + "\nCh lim: " + str(sy)
 692
 693          self.myform['lmst_high'].set_value(s)
 694
 695          #
 696          # Write data out to recording files
 697          #
 698          if self.setimode == False:
 699              self.write_continuum_data(x,sidtime)
 700              self.write_spectral_data(self.fft_outbuf,sidtime)
 701
 702          else:
 703              self.seti_analysis(self.fft_outbuf,sidtime)
 704              now = time.time()
 705              if ((self.scanning == True) and ((now - self.seti_then) > self.setifreq_timer)):
 706                  self.seti_then = now
 707                  self.setifreq_current = self.setifreq_current + self.fft_input_rate
 708                  if (self.setifreq_current > self.setifreq_upper):
 709                      self.setifreq_current = self.setifreq_lower
 710                  self.set_freq(self.setifreq_current)
 711                  # Make sure we zero-out the hits array when changing
 712                  #   frequency.
 713                  self.hits_array[:,:] = 0.0
 714                  self.hit_intensities[:,:] = 0.0
 715
 716     def fft_outfunc(self,data,l):
 717         self.fft_outbuf=data
 718
 719     def write_continuum_data(self,data,sidtime):
 720
 721         # Create localtime structure for producing filename
 722         foo = time.localtime()
 723         pfx = self.prefix
 724         filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 725            foo.tm_mon, foo.tm_mday, foo.tm_hour)
 726
 727         # Open the data file, appending
 728         continuum_file = open (filenamestr+".tpdat","a")
 729
 730         flt = "%6.3f" % data
 731         inter = self.decln
 732         integ = self.integ
 733         fc = self.observing
 734         fc = fc / 1000000
 735         bw = self.bw
 736         bw = bw / 1000000
 737         ga = self.gain
 738
 739         now = time.time()
 740
 741         #
 742         # If time to write full header info (saves storage this way)
 743         #
 744         if (now - self.continuum_then > 20):
 745             self.sun.compute(self.locality)
 746             enow = ephem.now()
 747             sun_insky = "Down"
 748             self.sunstate = "Dn"
 749             if ((self.sun.rise_time < enow) and (enow < self.sun.set_time)):
 750                sun_insky = "Up"
 751                self.sunstate = "Up"
 752             self.continuum_then = now
 753
 754             continuum_file.write(str(ephem.hours(sidtime))+" "+flt+" Dn="+str(inter)+",")
 755             continuum_file.write("Ti="+str(integ)+",Fc="+str(fc)+",Bw="+str(bw))
 756             continuum_file.write(",Ga="+str(ga)+",Sun="+str(sun_insky)+"\n")
 757         else:
 758             continuum_file.write(str(ephem.hours(sidtime))+" "+flt+"\n")
 759
 760         continuum_file.close()
 761         return(data)
 762
 763     def write_spectral_data(self,data,sidtime):
 764
 765         now = time.time()
 766         delta = 10
 767
 768         # If time to write out spectral data
 769         # We don't write this out every time, in order to
 770         #   save disk space.  Since the spectral data are
 771         #   typically heavily averaged, writing this data
 772         #   "once in a while" is OK.
 773         #
 774         if (now - self.spectral_then >= delta):
 775             self.spectral_then = now
 776
 777             # Get localtime structure to make filename from
 778             foo = time.localtime()
 779
 780             pfx = self.prefix
 781             filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 782                foo.tm_mon, foo.tm_mday, foo.tm_hour)
 783
 784             # Open the file
 785             spectral_file = open (filenamestr+".sdat","a")
 786
 787             # Setup data fields to be written
 788             r = data
 789             inter = self.decln
 790             fc = self.observing
 791             fc = fc / 1000000
 792             bw = self.bw
 793             bw = bw / 1000000
 794             av = self.avg_alpha
 795
 796             # Write those fields
 797             spectral_file.write("data:"+str(ephem.hours(sidtime))+" Dn="+str(inter)+",Fc="+str(fc)+",Bw="+str(bw)+",Av="+str(av))
 798             spectral_file.write (" [ ")
 799             for r in data:
 800                 spectral_file.write(" "+str(r))
 801
 802             spectral_file.write(" ]\n")
 803             spectral_file.close()
 804             return(data)
 805
 806         return(data)
 807
 808     def seti_analysis(self,fftbuf,sidtime):
 809         l = len(fftbuf)
 810         x = 0
 811         hits = []
 812         hit_intensities = []
 813         if self.seticounter < self.setitimer:
 814             self.seticounter = self.seticounter + 1
 815             return
 816         else:
 817             self.seticounter = 0
 818
 819         # Run through FFT output buffer, computing standard deviation (Sigma)
 820         avg = 0
 821         # First compute average
 822         for i in range(0,l):
 823             avg = avg + fftbuf[i]
 824         avg = avg / l
 825
 826         sigma = 0.0
 827         # Then compute standard deviation (Sigma)
 828         for i in range(0,l):
 829             d = fftbuf[i] - avg
 830             sigma = sigma + (d*d)
 831
 832         sigma = Numeric.sqrt(sigma/l)
 833
 834         #
 835         # Snarfle through the FFT output buffer again, looking for
 836         #    outlying data points
 837
 838         start_f = self.observing - (self.fft_input_rate/2)
 839         current_f = start_f
 840         l = len(fftbuf)
 841         f_incr = self.fft_input_rate / l
 842         hit = -1
 843
 844         # -nyquist to DC
 845         for i in range(l/2,l):
 846             #
 847             # If current FFT buffer has an item that exceeds the specified
 848             #  sigma
 849             #
 850             if ((fftbuf[i] - avg) > (self.setik * sigma)):
 851                 hits.append(current_f)
 852                 hit_intensities.append(fftbuf[i])
 853             current_f = current_f + f_incr
 854
 855         # DC to nyquist
 856         for i in range(0,l/2):
 857             #
 858             # If current FFT buffer has an item that exceeds the specified
 859             #  sigma
 860             #
 861             if ((fftbuf[i] - avg) > (self.setik * sigma)):
 862                 hits.append(current_f)
 863                 hit_intensities.append(fftbuf[i])
 864             current_f = current_f + f_incr
 865
 866         # No hits
 867         if (len(hits) <= 0):
 868             return
 869
 870
 871         #
 872         # OK, so we have some hits in the FFT buffer
 873         #   They'll have a rather substantial gauntlet to run before
 874         #   being declared a real "hit"
 875         #
 876
 877         # Update stats
 878         self.s1hitcounter = self.s1hitcounter + len(hits)
 879
 880         # Weed out buffers with an excessive number of
 881         #   signals above Sigma
 882         if (len(hits) > self.nhits):
 883             return
 884
 885
 886         # Weed out FFT buffers with apparent multiple narrowband signals
 887         #   separated significantly in frequency.  This means that a
 888         #   single signal spanning multiple bins is OK, but a buffer that
 889         #   has multiple, apparently-separate, signals isn't OK.
 890         #
 891         last = hits[0]
 892         ns2 = 1
 893         for i in range(1,len(hits)):
 894             if ((hits[i] - last) > (f_incr*3.0)):
 895                 return
 896             last = hits[i]
 897             ns2 = ns2 + 1
 898
 899         self.s2hitcounter = self.s2hitcounter + ns2
 900
 901         #
 902         # Run through all available hit buffers, computing difference between
 903         #   frequencies found there, if they're all within the chirp limits
 904         #   declare a good hit
 905         #
 906         good_hit = False
 907         f_ds = Numeric.zeros(self.nhitlines, Numeric.Float64)
 908         avg_delta = 0
 909         k = 0
 910         for i in range(0,min(len(hits),len(self.hits_array[:,0]))):
 911             f_ds[0] = abs(self.hits_array[i,0] - hits[i])
 912             for j in range(1,len(f_ds)):
 913                f_ds[j] = abs(self.hits_array[i,j] - self.hits_array[i,0])
 914                avg_delta = avg_delta + f_ds[j]
 915                k = k + 1
 916
 917             if (self.seti_isahit (f_ds)):
 918                 good_hit = True
 919                 self.hitcounter = self.hitcounter + 1
 920                 break
 921
 922         if (avg_delta/k < (self.seti_fft_bandwidth/2)):
 923             self.avgdelta = avg_delta / k
 924
 925         # Save 'n shuffle hits
 926         #  Old hit buffers percolate through the hit buffers
 927         #  (there are self.nhitlines of these buffers)
 928         #  and then drop off the end
 929         #  A consequence is that while the nhitlines buffers are filling,
 930         #  you can get some absurd values for self.avgdelta, because some
 931         #  of the buffers are full of zeros
 932         for i in range(self.nhitlines,1):
 933             self.hits_array[:,i] = self.hits_array[:,i-1]
 934             self.hit_intensities[:,i] = self.hit_intensities[:,i-1]
 935
 936         for i in range(0,len(hits)):
 937             self.hits_array[i,0] = hits[i]
 938             self.hit_intensities[i,0] = hit_intensities[i]
 939
 940         # Finally, write the hits/intensities buffer
 941         if (good_hit):
 942             self.write_hits(sidtime)
 943
 944         return
 945
 946     def seti_isahit(self,fdiffs):
 947         truecount = 0
 948
 949         for i in range(0,len(fdiffs)):
 950             if (fdiffs[i] >= self.CHIRP_LOWER and fdiffs[i] <= self.CHIRP_UPPER):
 951                 truecount = truecount + 1
 952
 953         if truecount == len(fdiffs):
 954             return (True)
 955         else:
 956             return (False)
 957
 958     def write_hits(self,sidtime):
 959         # Create localtime structure for producing filename
 960         foo = time.localtime()
 961         pfx = self.prefix
 962         filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
 963            foo.tm_mon, foo.tm_mday, foo.tm_hour)
 964
 965         # Open the data file, appending
 966         hits_file = open (filenamestr+".seti","a")
 967
 968         # Write sidtime first
 969         hits_file.write(str(ephem.hours(sidtime))+" "+str(self.decln)+" ")
 970
 971         #
 972         # Then write the hits/hit intensities buffers with enough
 973         #   "syntax" to allow parsing by external (not yet written!)
 974         #   "stuff".
 975         #
 976         for i in range(0,self.nhitlines):
 977             hits_file.write(" ")
 978             for j in range(0,self.nhits):
 979                 hits_file.write(str(self.hits_array[j,i])+":")
 980                 hits_file.write(str(self.hit_intensities[j,i])+",")
 981         hits_file.write("\n")
 982         hits_file.close()
 983         return
 984
 985     def xydfunc(self,xyv):
 986         if self.setimode == True:
 987             return
 988         magn = int(Numeric.log10(self.observing))
 989         if (magn == 6 or magn == 7 or magn == 8):
 990             magn = 6
 991         dfreq = xyv[0] * pow(10.0,magn)
 992         ratio = self.observing / dfreq
 993         vs = 1.0 - ratio
 994         vs *= 299792.0
 995         if magn >= 9:
 996            xhz = "Ghz"
 997         elif magn >= 6:
 998            xhz = "Mhz"
 999         elif magn <= 5:
1000            xhz =  "Khz"
1001         s = "%.6f%s\n%.3fdB" % (xyv[0], xhz, xyv[1])
1002         s2 = "\n%.3fkm/s" % vs
1003         self.myform['spec_data'].set_value(s+s2)
1004
1005     def xydfunc_waterfall(self,pos):
1006         lower = self.observing - (self.seti_fft_bandwidth / 2)
1007         upper = self.observing + (self.seti_fft_bandwidth / 2)
1008         binwidth = self.seti_fft_bandwidth / 1024
1009         s = "%.6fMHz" % ((lower + (pos.x*binwidth)) / 1.0e6)
1010         self.myform['spec_data'].set_value(s)
1011
1012     def toggle_cal(self):
1013         if (self.calstate == True):
1014           self.calstate = False
1015           self.u.write_io(0,0,(1<<15))
1016           self.calibrator.SetLabel("Calibration Source: Off")
1017         else:
1018           self.calstate = True
1019           self.u.write_io(0,(1<<15),(1<<15))
1020           self.calibrator.SetLabel("Calibration Source: On")
1021
1022     def toggle_annotation(self):
1023         if (self.annotate_state == True):
1024           self.annotate_state = False
1025           self.annotation.SetLabel("Annotation: Off")
1026         else:
1027           self.annotate_state = True
1028           self.annotation.SetLabel("Annotation: On")
1029     #
1030     # Turn scanning on/off
1031     # Called-back by "Recording" button
1032     #
1033     def toggle_scanning(self):
1034         # Current scanning?  Flip state
1035         if (self.scanning == True):
1036           self.scanning = False
1037           self.scan_control.SetLabel("Scan: Off")
1038         # Not scanning
1039         else:
1040           self.scanning = True
1041           self.scan_control.SetLabel("Scan: On ")
1042
1043     def set_pd_offset(self,offs):
1044          self.myform['offset'].set_value(offs)
1045          self.calib_offset=offs
1046          x = self.calib_coeff / 100.0
1047          self.cal_offs.set_k(offs*(x*8000))
1048
1049     def set_pd_gain(self,gain):
1050          self.myform['dcgain'].set_value(gain)
1051          self.cal_mult.set_k(gain*0.01)
1052          self.calib_coeff = gain
1053          self.cal_offs.set_k(self.calib_offset*(self.calib_coeff*8000))
1054
1055     def compute_notch_taps(self,notchlist):
1056          NOTCH_TAPS = 256
1057          tmptaps = Numeric.zeros(NOTCH_TAPS,Numeric.Complex64)
1058          binwidth = self.bw / NOTCH_TAPS
1059
1060          for i in range(0,NOTCH_TAPS):
1061              tmptaps[i] = complex(1.0,0.0)
1062
1063          for i in notchlist:
1064              diff = i - self.observing
1065              if i == 0:
1066                  break
1067              if (diff > 0):
1068                  idx = diff / binwidth
1069                  idx = int(idx)
1070                  if (idx < 0 or idx > (NOTCH_TAPS/2)):
1071                      break
1072                  tmptaps[idx] = complex(0.0, 0.0)
1073
1074              if (diff < 0):
1075                  idx = -diff / binwidth
1076                  idx = (NOTCH_TAPS/2) - idx
1077                  idx = int(idx+(NOTCH_TAPS/2))
1078                  if (idx < 0 or idx > (NOTCH_TAPS)):
1079                      break
1080                  tmptaps[idx] = complex(0.0, 0.0)
1081
1082          self.notch_taps = FFT.inverse_fft(tmptaps)
1083
1084 def main ():
1085     app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY SPECTRAL/CONTINUUM RECEIVER: $Revision$", nstatus=1)
1086     app.MainLoop()
1087
1088 if __name__ == '__main__':
1089     main ()