3 # Copyright 2004,2005,2007 Free Software Foundation, Inc.
5 # This file is part of GNU Radio
7 # GNU Radio is free software; you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation; either version 3, or (at your option)
12 # GNU Radio is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
17 # You should have received a copy of the GNU General Public License
18 # along with GNU Radio; see the file COPYING. If not, write to
19 # the Free Software Foundation, Inc., 51 Franklin Street,
20 # Boston, MA 02110-1301, USA.
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 stdgui2, ra_fftsink, ra_stripchartsink, ra_waterfallsink, form, slider
29 from optparse import OptionParser
37 class app_flow_graph(stdgui2.std_top_block):
38 def __init__(self, frame, panel, vbox, argv):
39 stdgui2.std_top_block.__init__(self, frame, panel, vbox, argv)
44 parser = OptionParser(option_class=eng_option)
45 parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=(0, 0),
46 help="select USRP Rx side A or B (default=A)")
47 parser.add_option("-d", "--decim", type="int", default=16,
48 help="set fgpa decimation rate to DECIM [default=%default]")
49 parser.add_option("-f", "--freq", type="eng_float", default=None,
50 help="set frequency to FREQ", metavar="FREQ")
51 parser.add_option("-a", "--avg", type="eng_float", default=1.0,
52 help="set spectral averaging alpha")
53 parser.add_option("-i", "--integ", type="eng_float", default=1.0,
54 help="set integration time")
55 parser.add_option("-g", "--gain", type="eng_float", default=None,
56 help="set gain in dB (default is midpoint)")
57 parser.add_option("-l", "--reflevel", type="eng_float", default=30.0,
58 help="Set Total power reference level")
59 parser.add_option("-y", "--division", type="eng_float", default=0.5,
60 help="Set Total power Y division size")
61 parser.add_option("-e", "--longitude", type="eng_float", default=-76.02,help="Set Observer Longitude")
62 parser.add_option("-c", "--latitude", type="eng_float", default=44.85,help="Set Observer Latitude")
63 parser.add_option("-o", "--observing", type="eng_float", default=0.0,
64 help="Set observing frequency")
65 parser.add_option("-x", "--ylabel", default="dB", help="Y axis label")
66 parser.add_option("-z", "--divbase", type="eng_float", default=0.025, help="Y Division increment base")
67 parser.add_option("-v", "--stripsize", type="eng_float", default=2400, help="Size of stripchart, in 2Hz samples")
68 parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT")
69 parser.add_option("-N", "--decln", type="eng_float", default=999.99, help="Observing declination")
70 parser.add_option("-X", "--prefix", default="./")
71 parser.add_option("-M", "--fft_rate", type="eng_float", default=8.0, help="FFT Rate")
72 parser.add_option("-A", "--calib_coeff", type="eng_float", default=1.0, help="Calibration coefficient")
73 parser.add_option("-B", "--calib_offset", type="eng_float", default=0.0, help="Calibration coefficient")
74 parser.add_option("-W", "--waterfall", action="store_true", default=False, help="Use Waterfall FFT display")
75 parser.add_option("-S", "--setimode", action="store_true", default=False, help="Enable SETI processing of spectral data")
76 parser.add_option("-K", "--setik", type="eng_float", default=1.5, help="K value for SETI analysis")
77 parser.add_option("-T", "--setibandwidth", type="eng_float", default=12500, help="Instantaneous SETI observing bandwidth--must be divisor of 250Khz")
78 parser.add_option("-Q", "--seti_range", type="eng_float", default=1.0e6, help="Total scan width, in Hz for SETI scans")
79 parser.add_option("-Z", "--dual_mode", action="store_true",
80 default=False, help="Dual-polarization mode")
81 parser.add_option("-I", "--interferometer", action="store_true", default=False, help="Interferometer mode")
82 parser.add_option("-D", "--switch_mode", action="store_true", default=False, help="Dicke Switching mode")
83 parser.add_option("-P", "--reference_divisor", type="eng_float", default=1.0, help="Reference Divisor")
84 parser.add_option("-U", "--ref_fifo", default="@@@@")
85 (options, args) = parser.parse_args()
87 self.setimode = options.setimode
88 self.dual_mode = options.dual_mode
89 self.interferometer = options.interferometer
90 self.normal_mode = False
91 self.switch_mode = options.switch_mode
93 self.reference_divisor = options.reference_divisor
94 self.ref_fifo = options.ref_fifo
96 if (self.ref_fifo != "@@@@"):
97 self.ref_fifo_file = open (self.ref_fifo, "w")
100 for modes in (self.dual_mode, self.interferometer):
102 modecount = modecount + 1
105 print "must select only 1 of --dual_mode, or --interferometer"
108 self.chartneeded = True
110 if (self.setimode == True):
111 self.chartneeded = False
113 if (self.setimode == True and self.interferometer == True):
114 print "can't pick both --setimode and --interferometer"
117 if (self.setimode == True and self.switch_mode == True):
118 print "can't pick both --setimode and --switch_mode"
121 if (self.interferometer == True and self.switch_mode == True):
122 print "can't pick both --interferometer and --switch_mode"
126 self.normal_mode = True
128 self.show_debug_info = True
130 # Pick up waterfall option
131 self.waterfall = options.waterfall
134 self.setimode = options.setimode
136 self.setik = options.setik
137 self.seti_fft_bandwidth = int(options.setibandwidth)
140 binwidth = self.seti_fft_bandwidth / options.fft_size
142 # Use binwidth, and knowledge of likely chirp rates to set reasonable
143 # values for SETI analysis code. We assume that SETI signals will
144 # chirp at somewhere between 0.10Hz/sec and 0.25Hz/sec.
146 # upper_limit is the "worst case"--that is, the case for which we have
147 # to wait the longest to actually see any drift, due to the quantizing
149 upper_limit = binwidth / 0.10
150 self.setitimer = int(upper_limit * 2.00)
153 # Calculate the CHIRP values based on Hz/sec
154 self.CHIRP_LOWER = 0.10 * self.setitimer
155 self.CHIRP_UPPER = 0.25 * self.setitimer
157 # Reset hit counters to 0
159 self.s1hitcounter = 0
160 self.s2hitcounter = 0
162 # We scan through 2Mhz of bandwidth around the chosen center freq
163 self.seti_freq_range = options.seti_range
164 # Calculate lower edge
165 self.setifreq_lower = options.freq - (self.seti_freq_range/2)
166 self.setifreq_current = options.freq
167 # Calculate upper edge
168 self.setifreq_upper = options.freq + (self.seti_freq_range/2)
170 # Maximum "hits" in a line
173 # Number of lines for analysis
176 # We change center frequencies based on nhitlines and setitimer
177 self.setifreq_timer = self.setitimer * (self.nhitlines * 5)
179 # Create actual timer
180 self.seti_then = time.time()
182 # The hits recording array
183 self.hits_array = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64)
184 self.hit_intensities = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64)
185 # Calibration coefficient and offset
186 self.calib_coeff = options.calib_coeff
187 self.calib_offset = options.calib_offset
188 if self.calib_offset < -750:
189 self.calib_offset = -750
190 if self.calib_offset > 750:
191 self.calib_offset = 750
193 if self.calib_coeff < 1:
195 if self.calib_coeff > 100:
196 self.calib_coeff = 100
198 self.integ = options.integ
199 self.avg_alpha = options.avg
200 self.gain = options.gain
201 self.decln = options.decln
203 # Set initial values for datalogging timed-output
204 self.continuum_then = time.time()
205 self.spectral_then = time.time()
210 self.subdev = [(0, 0), (0,0)]
213 # If SETI mode, we always run at maximum USRP decimation
218 if (self.dual_mode == False and self.interferometer == False):
219 self.u = usrp.source_c(decim_rate=options.decim)
220 self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec))
221 # determine the daughterboard subdevice we're using
222 self.subdev[0] = usrp.selected_subdev(self.u, options.rx_subdev_spec)
223 self.subdev[1] = self.subdev[0]
224 self.cardtype = self.subdev[0].dbid()
226 self.u=usrp.source_c(decim_rate=options.decim, nchan=2)
227 self.subdev[0] = usrp.selected_subdev(self.u, (0, 0))
228 self.subdev[1] = usrp.selected_subdev(self.u, (1, 0))
229 self.cardtype = self.subdev[0].dbid()
230 self.u.set_mux(0x32103210)
231 c1 = self.subdev[0].name()
232 c2 = self.subdev[1].name()
234 print "Must have identical cardtypes for --dual_mode or --interferometer"
241 format = self.u.make_format(width, shift)
242 r = self.u.set_format(format)
244 # Set initial declination
245 self.decln = options.decln
247 input_rate = self.u.adc_freq() / self.u.decim_rate()
250 # Set prefix for data files
252 self.prefix = options.prefix
255 # The lower this number, the fewer sample frames are dropped
256 # in computing the FFT. A sampled approach is taken to
257 # computing the FFT of the incoming data, which reduces
258 # sensitivity. Increasing sensitivity inreases CPU loading.
260 self.fft_rate = options.fft_rate
262 self.fft_size = int(options.fft_size)
264 # This buffer is used to remember the most-recent FFT display
265 # values. Used later by self.write_spectral_data() to write
266 # spectral data to datalogging files, and by the SETI analysis
269 self.fft_outbuf = Numeric.zeros(self.fft_size, Numeric.Float64)
272 # If SETI mode, only look at seti_fft_bandwidth
276 self.fft_input_rate = self.seti_fft_bandwidth
279 # Build a decimating bandpass filter
281 self.fft_input_taps = gr.firdes.complex_band_pass (1.0,
283 -(int(self.fft_input_rate/2)), int(self.fft_input_rate/2), 200,
284 gr.firdes.WIN_HAMMING, 0)
287 # Compute required decimation factor
289 decimation = int(input_rate/self.fft_input_rate)
290 self.fft_bandpass = gr.fir_filter_ccc (decimation,
293 self.fft_input_rate = input_rate
296 if self.waterfall == False:
297 self.scope = ra_fftsink.ra_fft_sink_c (panel,
298 fft_size=int(self.fft_size), sample_rate=self.fft_input_rate,
299 fft_rate=int(self.fft_rate), title="Spectral",
300 ofunc=self.fft_outfunc, xydfunc=self.xydfunc)
302 self.scope = ra_waterfallsink.waterfall_sink_c (panel,
303 fft_size=int(self.fft_size), sample_rate=self.fft_input_rate,
304 fft_rate=int(self.fft_rate), title="Spectral", ofunc=self.fft_outfunc, size=(1100, 600), xydfunc=self.xydfunc, ref_level=0, span=10)
306 # Set up ephemeris data
307 self.locality = ephem.Observer()
308 self.locality.long = str(options.longitude)
309 self.locality.lat = str(options.latitude)
311 # We make notes about Sunset/Sunrise in Continuum log files
312 self.sun = ephem.Sun()
315 # Set up stripchart display
317 if (self.dual_mode != False):
318 tit = "H+V Continuum"
319 if (self.interferometer != False):
320 tit = "East x West Correlation"
321 self.stripsize = int(options.stripsize)
322 if self.chartneeded == True:
323 self.chart = ra_stripchartsink.stripchart_sink_f (panel,
324 stripsize=self.stripsize,
326 xlabel="LMST Offset (Seconds)",
327 scaling=1.0, ylabel=options.ylabel,
328 divbase=options.divbase)
330 # Set center frequency
331 self.centerfreq = options.freq
333 # Set observing frequency (might be different from actual programmed
335 if options.observing == 0.0:
336 self.observing = options.freq
338 self.observing = options.observing
340 # Remember our input bandwidth
345 # The strip chart is fed at a constant 1Hz rate
349 # Call constructors for receive chains
352 if (self.dual_mode == True):
353 self.setup_dual (self.setimode)
355 if (self.interferometer == True):
356 self.setup_interferometer(self.setimode)
358 if (self.normal_mode == True):
359 self.setup_normal(self.setimode)
361 if (self.setimode == True):
364 self._build_gui(vbox)
366 # Make GUI agree with command-line
367 self.integ = options.integ
368 if self.setimode == False:
369 self.myform['integration'].set_value(int(options.integ))
370 self.myform['offset'].set_value(self.calib_offset)
371 self.myform['dcgain'].set_value(self.calib_coeff)
372 self.myform['average'].set_value(int(options.avg))
375 if self.setimode == False:
376 # Make integrator agree with command line
377 self.set_integration(int(options.integ))
379 self.avg_alpha = options.avg
381 # Make spectral averager agree with command line
382 if options.avg != 1.0:
383 self.scope.set_avg_alpha(float(1.0/options.avg))
384 self.scope.set_average(True)
386 if self.setimode == False:
388 self.chart.set_y_per_div(options.division)
389 # Set reference(MAX) level
390 self.chart.set_ref_level(options.reflevel)
394 if options.gain is None:
395 # if no gain was specified, use the mid-point in dB
396 g = self.subdev[0].gain_range()
397 options.gain = float(g[0]+g[1])/2
399 if options.freq is None:
400 # if no freq was specified, use the mid-point
401 r = self.subdev[0].freq_range()
402 options.freq = float(r[0]+r[1])/2
404 # Set the initial gain control
405 self.set_gain(options.gain)
407 if not(self.set_freq(options.freq)):
408 self._set_status_msg("Failed to set initial frequency")
411 self.set_decln (self.decln)
414 # RF hardware information
415 self.myform['decim'].set_value(self.u.decim_rate())
416 self.myform['USB BW'].set_value(self.u.adc_freq() / self.u.decim_rate())
417 if (self.dual_mode == True):
418 self.myform['dbname'].set_value(self.subdev[0].name()+'/'+self.subdev[1].name())
420 self.myform['dbname'].set_value(self.subdev[0].name())
422 # Set analog baseband filtering, if DBS_RX
423 if self.cardtype in (usrp_dbid.DBS_RX, usrp_dbid.DBS_RX_REV_2_1):
424 lbw = (self.u.adc_freq() / self.u.decim_rate()) / 2
427 self.subdev[0].set_bw(lbw)
428 self.subdev[1].set_bw(lbw)
430 # Start the timer for the LMST display and datalogging
431 self.lmst_timer.Start(1000)
432 if (self.switch_mode == True):
433 self.other_timer.Start(330)
436 def _set_status_msg(self, msg):
437 self.frame.GetStatusBar().SetStatusText(msg, 0)
439 def _build_gui(self, vbox):
441 def _form_set_freq(kv):
442 # Adjust current SETI frequency, and limits
443 self.setifreq_lower = kv['freq'] - (self.seti_freq_range/2)
444 self.setifreq_current = kv['freq']
445 self.setifreq_upper = kv['freq'] + (self.seti_freq_range/2)
447 # Reset SETI analysis timer
448 self.seti_then = time.time()
449 # Zero-out hits array when changing frequency
450 self.hits_array[:,:] = 0.0
451 self.hit_intensities[:,:] = -60.0
453 return self.set_freq(kv['freq'])
455 def _form_set_decln(kv):
456 return self.set_decln(kv['decln'])
458 # Position the FFT display
459 vbox.Add(self.scope.win, 15, wx.EXPAND)
461 if self.setimode == False:
462 # Position the Total-power stripchart
463 vbox.Add(self.chart.win, 15, wx.EXPAND)
465 # add control area at the bottom
466 self.myform = myform = form.form()
467 hbox = wx.BoxSizer(wx.HORIZONTAL)
468 hbox.Add((7,0), 0, wx.EXPAND)
469 vbox1 = wx.BoxSizer(wx.VERTICAL)
470 myform['freq'] = form.float_field(
471 parent=self.panel, sizer=vbox1, label="Center freq", weight=1,
472 callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg))
474 vbox1.Add((4,0), 0, 0)
476 myform['lmst_high'] = form.static_text_field(
477 parent=self.panel, sizer=vbox1, label="Current LMST", weight=1)
478 vbox1.Add((4,0), 0, 0)
480 if self.setimode == False:
481 myform['spec_data'] = form.static_text_field(
482 parent=self.panel, sizer=vbox1, label="Spectral Cursor", weight=1)
483 vbox1.Add((4,0), 0, 0)
485 vbox2 = wx.BoxSizer(wx.VERTICAL)
486 if self.setimode == False:
487 vbox3 = wx.BoxSizer(wx.VERTICAL)
488 g = self.subdev[0].gain_range()
489 myform['gain'] = form.slider_field(parent=self.panel, sizer=vbox2, label="RF Gain",
491 min=int(g[0]), max=int(g[1]),
492 callback=self.set_gain)
494 vbox2.Add((4,0), 0, 0)
495 if self.setimode == True:
499 myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2,
500 label="Spectral Averaging (FFT frames)", weight=1, min=1, max=max_savg, callback=self.set_averaging)
502 # Set up scan control button when in SETI mode
503 if (self.setimode == True):
504 # SETI scanning control
505 buttonbox = wx.BoxSizer(wx.HORIZONTAL)
506 self.scan_control = form.button_with_callback(self.panel,
508 callback=self.toggle_scanning)
510 buttonbox.Add(self.scan_control, 0, wx.CENTER)
511 vbox2.Add(buttonbox, 0, wx.CENTER)
513 vbox2.Add((4,0), 0, 0)
515 if self.setimode == False:
516 myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2,
517 label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration)
519 vbox2.Add((4,0), 0, 0)
521 myform['decln'] = form.float_field(
522 parent=self.panel, sizer=vbox2, label="Current Declination", weight=1,
523 callback=myform.check_input_and_call(_form_set_decln))
524 vbox2.Add((4,0), 0, 0)
526 if self.setimode == False:
527 myform['offset'] = form.slider_field(parent=self.panel, sizer=vbox3,
528 label="Post-Detector Offset", weight=1, min=-750, max=750,
529 callback=self.set_pd_offset)
530 vbox3.Add((2,0), 0, 0)
531 myform['dcgain'] = form.slider_field(parent=self.panel, sizer=vbox3,
532 label="Post-Detector Gain", weight=1, min=1, max=100,
533 callback=self.set_pd_gain)
534 vbox3.Add((2,0), 0, 0)
535 hbox.Add(vbox1, 0, 0)
536 hbox.Add(vbox2, wx.ALIGN_RIGHT, 0)
538 if self.setimode == False:
539 hbox.Add(vbox3, wx.ALIGN_RIGHT, 0)
541 vbox.Add(hbox, 0, wx.EXPAND)
543 self._build_subpanel(vbox)
545 self.lmst_timer = wx.PyTimer(self.lmst_timeout)
546 self.other_timer = wx.PyTimer(self.other_timeout)
549 def _build_subpanel(self, vbox_arg):
550 # build a secondary information panel (sometimes hidden)
552 # FIXME figure out how to have this be a subpanel that is always
553 # created, but has its visibility controlled by foo.Show(True/False)
555 if not(self.show_debug_info):
562 #panel = wx.Panel(self.panel, -1)
563 #vbox = wx.BoxSizer(wx.VERTICAL)
565 hbox = wx.BoxSizer(wx.HORIZONTAL)
567 myform['decim'] = form.static_float_field(
568 parent=panel, sizer=hbox, label="Decim")
571 myform['USB BW'] = form.static_float_field(
572 parent=panel, sizer=hbox, label="USB BW")
575 myform['dbname'] = form.static_text_field(
576 parent=panel, sizer=hbox)
579 myform['baseband'] = form.static_float_field(
580 parent=panel, sizer=hbox, label="Analog BB")
583 myform['ddc'] = form.static_float_field(
584 parent=panel, sizer=hbox, label="DDC")
587 vbox.Add(hbox, 0, wx.EXPAND)
591 def set_freq(self, target_freq):
593 Set the center frequency we're interested in.
595 @param target_freq: frequency in Hz
598 Tuning is a two step process. First we ask the front-end to
599 tune as close to the desired frequency as it can. Then we use
600 the result of that operation and our target_frequency to
601 determine the value for the digital down converter.
604 # Everything except BASIC_RX should support usrp.tune()
606 if not (self.cardtype == usrp_dbid.BASIC_RX):
607 r = usrp.tune(self.u, self.subdev[0]._which, self.subdev[0], target_freq)
608 r = usrp.tune(self.u, self.subdev[1]._which, self.subdev[1], target_freq)
610 r = self.u.set_rx_freq(0, target_freq)
611 f = self.u.rx_freq(0)
612 if abs(f-target_freq) > 2.0e3:
615 self.myform['freq'].set_value(target_freq) # update displayed value
617 # Make sure calibrator knows our target freq
620 # Remember centerfreq---used for doppler calcs
621 delta = self.centerfreq - target_freq
622 self.centerfreq = target_freq
623 self.observing -= delta
624 self.scope.set_baseband_freq (self.observing)
626 self.myform['baseband'].set_value(r.baseband_freq)
627 self.myform['ddc'].set_value(r.dxc_freq)
633 def set_decln(self, dec):
635 self.myform['decln'].set_value(dec) # update displayed value
637 def set_gain(self, gain):
638 self.myform['gain'].set_value(gain) # update displayed value
639 self.subdev[0].set_gain(gain)
640 self.subdev[1].set_gain(gain)
643 def set_averaging(self, avval):
644 self.myform['average'].set_value(avval)
645 self.scope.set_avg_alpha(1.0/(avval))
646 self.scope.set_average(True)
647 self.avg_alpha = avval
649 def set_integration(self, integval):
650 if self.setimode == False:
651 self.integrator.set_taps(1.0/((integval)*(self.bw/2)))
652 self.myform['integration'].set_value(integval)
653 self.integ = integval
657 # Used to update LMST display, as well as current
660 # We also write external data-logging files here
662 def lmst_timeout(self):
663 self.locality.date = ephem.now()
664 if self.setimode == False:
665 x = self.probe.level()
666 sidtime = self.locality.sidereal_time()
668 s = str(ephem.hours(sidtime)) + " " + self.sunstate
669 # Continuum detector value
670 if self.setimode == False:
672 s = s + "\nDet: " + str(sx)
674 sx = "%2d" % self.hitcounter
675 s1 = "%2d" % self.s1hitcounter
676 s2 = "%2d" % self.s2hitcounter
677 sa = "%4.2f" % self.avgdelta
678 sy = "%3.1f-%3.1f" % (self.CHIRP_LOWER, self.CHIRP_UPPER)
679 s = s + "\nHits: " + str(sx) + "\nS1:" + str(s1) + " S2:" + str(s2)
680 s = s + "\nAv D: " + str(sa) + "\nCh lim: " + str(sy)
682 self.myform['lmst_high'].set_value(s)
685 # Write data out to recording files
687 if self.setimode == False:
688 self.write_continuum_data(x,sidtime)
689 self.write_spectral_data(self.fft_outbuf,sidtime)
692 self.seti_analysis(self.fft_outbuf,sidtime)
694 if ((self.scanning == True) and ((now - self.seti_then) > self.setifreq_timer)):
696 self.setifreq_current = self.setifreq_current + self.fft_input_rate
697 if (self.setifreq_current > self.setifreq_upper):
698 self.setifreq_current = self.setifreq_lower
699 self.set_freq(self.setifreq_current)
700 # Make sure we zero-out the hits array when changing
702 self.hits_array[:,:] = 0.0
703 self.hit_intensities[:,:] = 0.0
705 def other_timeout(self):
706 if (self.switch_state == 0):
707 self.switch_state = 1
709 elif (self.switch_state == 1):
710 self.switch_state = 0
712 if (self.switch_state == 0):
714 self.cmute.set_n(int(1.0e9))
716 elif (self.switch_state == 1):
717 self.mute.set_n(int(1.0e9))
720 if (self.ref_fifo != "@@@@"):
721 self.ref_fifo_file.write(str(self.switch_state)+"\n")
722 self.ref_fifo_file.flush()
724 self.avg_reference_value = self.cprobe.level()
727 # Set reference value
729 self.reference_level.set_k(-1.0 * (self.avg_reference_value/self.reference_divisor))
731 def fft_outfunc(self,data,l):
734 def write_continuum_data(self,data,sidtime):
736 # Create localtime structure for producing filename
737 foo = time.localtime()
739 filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
740 foo.tm_mon, foo.tm_mday, foo.tm_hour)
742 # Open the data file, appending
743 continuum_file = open (filenamestr+".tpdat","a")
757 # If time to write full header info (saves storage this way)
759 if (now - self.continuum_then > 20):
760 self.sun.compute(self.locality)
764 if ((self.sun.rise_time < enow) and (enow < self.sun.set_time)):
767 self.continuum_then = now
769 continuum_file.write(str(ephem.hours(sidtime))+" "+flt+" Dn="+str(inter)+",")
770 continuum_file.write("Ti="+str(integ)+",Fc="+str(fc)+",Bw="+str(bw))
771 continuum_file.write(",Ga="+str(ga)+",Sun="+str(sun_insky)+"\n")
773 continuum_file.write(str(ephem.hours(sidtime))+" "+flt+"\n")
775 continuum_file.close()
778 def write_spectral_data(self,data,sidtime):
783 # If time to write out spectral data
784 # We don't write this out every time, in order to
785 # save disk space. Since the spectral data are
786 # typically heavily averaged, writing this data
787 # "once in a while" is OK.
789 if (now - self.spectral_then >= delta):
790 self.spectral_then = now
792 # Get localtime structure to make filename from
793 foo = time.localtime()
796 filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
797 foo.tm_mon, foo.tm_mday, foo.tm_hour)
800 spectral_file = open (filenamestr+".sdat","a")
802 # Setup data fields to be written
812 spectral_file.write("data:"+str(ephem.hours(sidtime))+" Dn="+str(inter)+",Fc="+str(fc)+",Bw="+str(bw)+",Av="+str(av))
813 spectral_file.write (" [ ")
815 spectral_file.write(" "+str(r))
817 spectral_file.write(" ]\n")
818 spectral_file.close()
823 def seti_analysis(self,fftbuf,sidtime):
828 if self.seticounter < self.setitimer:
829 self.seticounter = self.seticounter + 1
834 # Run through FFT output buffer, computing standard deviation (Sigma)
836 # First compute average
838 avg = avg + fftbuf[i]
842 # Then compute standard deviation (Sigma)
845 sigma = sigma + (d*d)
847 sigma = Numeric.sqrt(sigma/l)
850 # Snarfle through the FFT output buffer again, looking for
851 # outlying data points
853 start_f = self.observing - (self.fft_input_rate/2)
856 f_incr = self.fft_input_rate / l
860 for i in range(l/2,l):
862 # If current FFT buffer has an item that exceeds the specified
865 if ((fftbuf[i] - avg) > (self.setik * sigma)):
866 hits.append(current_f)
867 hit_intensities.append(fftbuf[i])
868 current_f = current_f + f_incr
871 for i in range(0,l/2):
873 # If current FFT buffer has an item that exceeds the specified
876 if ((fftbuf[i] - avg) > (self.setik * sigma)):
877 hits.append(current_f)
878 hit_intensities.append(fftbuf[i])
879 current_f = current_f + f_incr
887 # OK, so we have some hits in the FFT buffer
888 # They'll have a rather substantial gauntlet to run before
889 # being declared a real "hit"
893 self.s1hitcounter = self.s1hitcounter + len(hits)
895 # Weed out buffers with an excessive number of
896 # signals above Sigma
897 if (len(hits) > self.nhits):
901 # Weed out FFT buffers with apparent multiple narrowband signals
902 # separated significantly in frequency. This means that a
903 # single signal spanning multiple bins is OK, but a buffer that
904 # has multiple, apparently-separate, signals isn't OK.
908 for i in range(1,len(hits)):
909 if ((hits[i] - last) > (f_incr*3.0)):
914 self.s2hitcounter = self.s2hitcounter + ns2
917 # Run through all available hit buffers, computing difference between
918 # frequencies found there, if they're all within the chirp limits
922 f_ds = Numeric.zeros(self.nhitlines, Numeric.Float64)
925 for i in range(0,min(len(hits),len(self.hits_array[:,0]))):
926 f_ds[0] = abs(self.hits_array[i,0] - hits[i])
927 for j in range(1,len(f_ds)):
928 f_ds[j] = abs(self.hits_array[i,j] - self.hits_array[i,0])
929 avg_delta = avg_delta + f_ds[j]
932 if (self.seti_isahit (f_ds)):
934 self.hitcounter = self.hitcounter + 1
937 if (avg_delta/k < (self.seti_fft_bandwidth/2)):
938 self.avgdelta = avg_delta / k
940 # Save 'n shuffle hits
941 # Old hit buffers percolate through the hit buffers
942 # (there are self.nhitlines of these buffers)
943 # and then drop off the end
944 # A consequence is that while the nhitlines buffers are filling,
945 # you can get some absurd values for self.avgdelta, because some
946 # of the buffers are full of zeros
947 for i in range(self.nhitlines,1):
948 self.hits_array[:,i] = self.hits_array[:,i-1]
949 self.hit_intensities[:,i] = self.hit_intensities[:,i-1]
951 for i in range(0,len(hits)):
952 self.hits_array[i,0] = hits[i]
953 self.hit_intensities[i,0] = hit_intensities[i]
955 # Finally, write the hits/intensities buffer
957 self.write_hits(sidtime)
961 def seti_isahit(self,fdiffs):
964 for i in range(0,len(fdiffs)):
965 if (fdiffs[i] >= self.CHIRP_LOWER and fdiffs[i] <= self.CHIRP_UPPER):
966 truecount = truecount + 1
968 if truecount == len(fdiffs):
973 def write_hits(self,sidtime):
974 # Create localtime structure for producing filename
975 foo = time.localtime()
977 filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year,
978 foo.tm_mon, foo.tm_mday, foo.tm_hour)
980 # Open the data file, appending
981 hits_file = open (filenamestr+".seti","a")
983 # Write sidtime first
984 hits_file.write(str(ephem.hours(sidtime))+" "+str(self.decln)+" ")
987 # Then write the hits/hit intensities buffers with enough
988 # "syntax" to allow parsing by external (not yet written!)
991 for i in range(0,self.nhitlines):
993 for j in range(0,self.nhits):
994 hits_file.write(str(self.hits_array[j,i])+":")
995 hits_file.write(str(self.hit_intensities[j,i])+",")
996 hits_file.write("\n")
1000 def xydfunc(self,xyv):
1001 if self.setimode == True:
1003 magn = int(Numeric.log10(self.observing))
1004 if (magn == 6 or magn == 7 or magn == 8):
1006 dfreq = xyv[0] * pow(10.0,magn)
1007 ratio = self.observing / dfreq
1016 s = "%.6f%s\n%.3fdB" % (xyv[0], xhz, xyv[1])
1017 s2 = "\n%.3fkm/s" % vs
1018 self.myform['spec_data'].set_value(s+s2)
1020 def xydfunc_waterfall(self,pos):
1021 lower = self.observing - (self.seti_fft_bandwidth / 2)
1022 upper = self.observing + (self.seti_fft_bandwidth / 2)
1023 binwidth = self.seti_fft_bandwidth / 1024
1024 s = "%.6fMHz" % ((lower + (pos.x*binwidth)) / 1.0e6)
1025 self.myform['spec_data'].set_value(s)
1027 def toggle_cal(self):
1028 if (self.calstate == True):
1029 self.calstate = False
1030 self.u.write_io(0,0,(1<<15))
1031 self.calibrator.SetLabel("Calibration Source: Off")
1033 self.calstate = True
1034 self.u.write_io(0,(1<<15),(1<<15))
1035 self.calibrator.SetLabel("Calibration Source: On")
1037 def toggle_annotation(self):
1038 if (self.annotate_state == True):
1039 self.annotate_state = False
1040 self.annotation.SetLabel("Annotation: Off")
1042 self.annotate_state = True
1043 self.annotation.SetLabel("Annotation: On")
1045 # Turn scanning on/off
1046 # Called-back by "Recording" button
1048 def toggle_scanning(self):
1049 # Current scanning? Flip state
1050 if (self.scanning == True):
1051 self.scanning = False
1052 self.scan_control.SetLabel("Scan: Off")
1055 self.scanning = True
1056 self.scan_control.SetLabel("Scan: On ")
1058 def set_pd_offset(self,offs):
1059 self.myform['offset'].set_value(offs)
1060 self.calib_offset=offs
1061 x = self.calib_coeff / 100.0
1062 self.cal_offs.set_k(offs*(x*8000))
1064 def set_pd_gain(self,gain):
1065 self.myform['dcgain'].set_value(gain)
1066 self.cal_mult.set_k(gain*0.01)
1067 self.calib_coeff = gain
1069 self.cal_offs.set_k(self.calib_offset*(x*8000))
1071 def compute_notch_taps(self,notchlist):
1073 tmptaps = Numeric.zeros(NOTCH_TAPS,Numeric.Complex64)
1074 binwidth = self.bw / NOTCH_TAPS
1076 for i in range(0,NOTCH_TAPS):
1077 tmptaps[i] = complex(1.0,0.0)
1080 diff = i - self.observing
1084 idx = diff / binwidth
1086 if (idx < 0 or idx > (NOTCH_TAPS/2)):
1088 tmptaps[idx] = complex(0.0, 0.0)
1091 idx = -diff / binwidth
1092 idx = (NOTCH_TAPS/2) - idx
1093 idx = int(idx+(NOTCH_TAPS/2))
1094 if (idx < 0 or idx > (NOTCH_TAPS)):
1096 tmptaps[idx] = complex(0.0, 0.0)
1098 self.notch_taps = numpy.fft.ifft(tmptaps)
1101 # Setup common pieces of radiometer mode
1103 def setup_radiometer_common(self):
1104 # The IIR integration filter for post-detection
1105 self.integrator = gr.single_pole_iir_filter_ff(1.0)
1106 self.integrator.set_taps (1.0/self.bw)
1109 self.probe = gr.probe_signal_f()
1112 # Continuum calibration stuff
1114 x = self.calib_coeff/100.0
1115 self.cal_mult = gr.multiply_const_ff(self.calib_coeff/100.0)
1116 self.cal_offs = gr.add_const_ff(self.calib_offset*(x*8000))
1119 # Mega decimator after IIR filter
1121 if (self.switch_mode == False):
1122 self.keepn = gr.keep_one_in_n(gr.sizeof_float, self.bw)
1124 self.keepn = gr.keep_one_in_n(gr.sizeof_float, int(self.bw/2))
1127 # For the Dicke-switching scheme
1129 self.switch = gr.multiply_const_ff(1.0)
1132 if (self.switch_mode == True):
1133 self.vector = gr.vector_sink_f()
1134 self.swkeep = gr.keep_one_in_n(gr.sizeof_float, int(self.bw/3))
1135 self.mute = gr.keep_one_in_n(gr.sizeof_float, 1)
1136 self.cmute = gr.keep_one_in_n(gr.sizeof_float, int(1.0e9))
1137 self.cintegrator = gr.single_pole_iir_filter_ff(1.0/(self.bw/2))
1138 self.cprobe = gr.probe_signal_f()
1140 self.mute = gr.multiply_const_ff(1.0)
1143 self.avg_reference_value = 0.0
1144 self.reference_level = gr.add_const_ff(0.0)
1147 # Setup ordinary single-channel radiometer mode
1149 def setup_normal(self, setimode):
1154 if setimode == False:
1155 self.detector = gr.complex_to_mag_squared()
1156 self.setup_radiometer_common()
1157 self.connect(self.shead, self.scope)
1159 self.connect(self.head, self.detector, self.mute, self.reference_level,
1160 self.integrator, self.keepn, self.cal_mult, self.cal_offs, self.chart)
1162 self.connect(self.cal_offs, self.probe)
1165 # Add a side-chain detector chain, with a different integrator, for sampling
1166 # The reference channel data
1167 # This is used to derive the offset value for self.reference_level, used above
1169 if (self.switch_mode == True):
1170 self.connect(self.detector, self.cmute, self.cintegrator, self.swkeep, self.cprobe)
1175 # Setup dual-channel (two antenna, usual orthogonal polarity probes in the same waveguide)
1177 def setup_dual(self, setimode):
1179 self.di = gr.deinterleave(gr.sizeof_gr_complex)
1180 self.addchans = gr.add_cc ()
1181 self.detector = gr.add_ff ()
1182 self.h_power = gr.complex_to_mag_squared()
1183 self.v_power = gr.complex_to_mag_squared()
1184 self.connect (self.u, self.di)
1187 # For spectral, adding the two channels works, assuming no gross
1188 # phase or amplitude error
1189 self.connect ((self.di, 0), (self.addchans, 0))
1190 self.connect ((self.di, 1), (self.addchans, 1))
1193 # Connect heads of spectral and total-power chains
1196 self.shead = self.addchans
1198 if (setimode == False):
1200 self.setup_radiometer_common()
1203 # For dual-polarization mode, we compute the sum of the
1204 # powers on each channel, after they've been detected
1206 self.detector = gr.add_ff()
1209 # In dual-polarization mode, we compute things a little differently
1210 # In effect, we have two radiometer chains, terminating in an adder
1212 self.connect((self.di, 0), self.h_power)
1213 self.connect((self.di, 1), self.v_power)
1214 self.connect(self.h_power, (self.detector, 0))
1215 self.connect(self.v_power, (self.detector, 1))
1216 self.connect(self.detector, self.mute, self.reference_level,
1217 self.integrator, self.keepn, self.cal_mult, self.cal_offs, self.chart)
1218 self.connect(self.cal_offs, self.probe)
1219 self.connect(self.shead, self.scope)
1222 # Add a side-chain detector chain, with a different integrator, for sampling
1223 # The reference channel data
1224 # This is used to derive the offset value for self.reference_level, used above
1226 if (self.switch_mode == True):
1227 self.connect(self.detector, self.cmute, self.cintegrator, self.swkeep, self.cprobe)
1231 # Setup correlating interferometer mode
1233 def setup_interferometer(self, setimode):
1234 self.setup_radiometer_common()
1236 self.di = gr.deinterleave(gr.sizeof_gr_complex)
1237 self.connect (self.u, self.di)
1238 self.corr = gr.multiply_cc()
1239 self.c2f = gr.complex_to_float()
1241 self.shead = (self.di, 0)
1243 # Channel 0 to multiply port 0
1244 # Channel 1 to multiply port 1
1245 self.connect((self.di, 0), (self.corr, 0))
1246 self.connect((self.di, 1), (self.corr, 1))
1249 # Multiplier (correlator) to complex-to-float, followed by integrator, etc
1251 self.connect(self.corr, self.c2f, self.switch, self.integrator, self.keepn, self.cal_mult, self.cal_offs, self.chart)
1254 # FFT scope gets only 1 channel
1255 # FIX THIS, by cross-correlating the *outputs* of two different FFTs, then display
1258 self.connect(self.shead, self.scope)
1261 # Output of correlator/integrator chain to probe
1263 self.connect(self.cal_offs, self.probe)
1270 def setup_seti(self):
1271 self.connect (self.shead, self.fft_bandpass, self.scope)
1277 app = stdgui2.stdapp(app_flow_graph, "RADIO ASTRONOMY SPECTRAL/CONTINUUM RECEIVER: $Revision$", nstatus=1)
1280 if __name__ == '__main__':