X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=gr-radio-astronomy%2Fsrc%2Fpython%2Fusrp_ra_receiver.py;h=a16f9d3f8b8f7eee905ede06008a1432843f1e43;hb=d53437a03ff992f1f3f9651248ac8aad94414fac;hp=5546a59e0f3b1d42f8aff11d87c4c92ac4839606;hpb=88d87e948315bc22ce02e635cb4771a836b73616;p=debian%2Fgnuradio diff --git a/gr-radio-astronomy/src/python/usrp_ra_receiver.py b/gr-radio-astronomy/src/python/usrp_ra_receiver.py index 5546a59e..a16f9d3f 100755 --- a/gr-radio-astronomy/src/python/usrp_ra_receiver.py +++ b/gr-radio-astronomy/src/python/usrp_ra_receiver.py @@ -16,23 +16,29 @@ # # You should have received a copy of the GNU General Public License # along with GNU Radio; see the file COPYING. If not, write to -# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, -# Boston, MA 02111-1307, USA. +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. # from gnuradio import gr, gru from gnuradio import usrp -import usrp_dbid +from usrpm import usrp_dbid from gnuradio import eng_notation from gnuradio.eng_option import eng_option -from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, waterfallsink, form, slider +from gnuradio.wxgui import stdgui, ra_fftsink, ra_stripchartsink, ra_waterfallsink, form, slider from optparse import OptionParser import wx import sys -from Numeric import * +import Numeric +import time import FFT import ephem -from gnuradio.local_calibrator import * + +class continuum_calibration(gr.feval_dd): + def eval(self, x): + str = globals()["calibration_codelet"] + exec(str) + return(x) class app_flow_graph(stdgui.gui_flow_graph): def __init__(self, frame, panel, vbox, argv): @@ -63,67 +69,188 @@ class app_flow_graph(stdgui.gui_flow_graph): parser.add_option("-o", "--observing", type="eng_float", default=0.0, help="Set observing frequency") parser.add_option("-x", "--ylabel", default="dB", help="Y axis label") - parser.add_option("-C", "--cfunc", default="default", help="Calibration function name") parser.add_option("-z", "--divbase", type="eng_float", default=0.025, help="Y Division increment base") parser.add_option("-v", "--stripsize", type="eng_float", default=2400, help="Size of stripchart, in 2Hz samples") parser.add_option("-F", "--fft_size", type="eng_float", default=1024, help="Size of FFT") parser.add_option("-N", "--decln", type="eng_float", default=999.99, help="Observing declination") - parser.add_option("-I", "--interfilt", action="store_true", default=False) parser.add_option("-X", "--prefix", default="./") + parser.add_option("-M", "--fft_rate", type="eng_float", default=8.0, help="FFT Rate") + parser.add_option("-A", "--calib_coeff", type="eng_float", default=1.0, help="Calibration coefficient") + parser.add_option("-B", "--calib_offset", type="eng_float", default=0.0, help="Calibration coefficient") + parser.add_option("-W", "--waterfall", action="store_true", default=False, help="Use Waterfall FFT display") + parser.add_option("-S", "--setimode", action="store_true", default=False, help="Enable SETI processing of spectral data") + parser.add_option("-K", "--setik", type="eng_float", default=1.5, help="K value for SETI analysis") + parser.add_option("-T", "--setibandwidth", type="eng_float", default=12500, help="Instantaneous SETI observing bandwidth--must be divisor of 250Khz") (options, args) = parser.parse_args() if len(args) != 0: parser.print_help() sys.exit(1) self.show_debug_info = True - + + # Pick up waterfall option + self.waterfall = options.waterfall + + # SETI mode stuff + self.setimode = options.setimode + self.seticounter = 0 + self.setik = options.setik + # Because we force the input rate to be 250Khz, 12.5Khz is + # exactly 1/20th of this, which makes building decimators + # easier. + # This also allows larger FFTs to be used without totally-gobbling + # CPU. With an FFT size of 16384, for example, this bandwidth + # yields a binwidth of 0.762Hz, and plenty of CPU left over + # for other things, like the SETI analysis code. + # + self.seti_fft_bandwidth = int(options.setibandwidth) + + # Calculate binwidth + binwidth = self.seti_fft_bandwidth / options.fft_size + + # Use binwidth, and knowledge of likely chirp rates to set reasonable + # values for SETI analysis code. We assume that SETI signals will + # chirp at somewhere between 0.10Hz/sec and 0.25Hz/sec. + # + # upper_limit is the "worst case"--that is, the case for which we have + # wait the longest to actually see any drift, due to the quantizing + # on FFT bins. + upper_limit = binwidth / 0.10 + self.setitimer = int(upper_limit * 2.00) + self.scanning = True + + # Calculate the CHIRP values based on Hz/sec + self.CHIRP_LOWER = 0.10 * self.setitimer + self.CHIRP_UPPER = 0.25 * self.setitimer + + # Reset hit counter to 0 + self.hitcounter = 0 + # We scan through 1Mhz of bandwidth around the chosen center freq + self.seti_freq_range = 1.0e6 + # Calculate lower edge + self.setifreq_lower = options.freq - (self.seti_freq_range/2) + self.setifreq_current = options.freq + # Calculate upper edge + self.setifreq_upper = options.freq + (self.seti_freq_range/2) + + # We change center frequencies every 10 self.setitimer intervals + self.setifreq_timer = self.setitimer * 10 + + # Create actual timer + self.seti_then = time.time() + + # The hits recording array + self.nhits = 10 + self.nhitlines = 3 + self.hits_array = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64) + self.hit_intensities = Numeric.zeros((self.nhits,self.nhitlines), Numeric.Float64) + # Calibration coefficient and offset + self.calib_coeff = options.calib_coeff + self.calib_offset = options.calib_offset + + self.integ = options.integ + self.avg_alpha = options.avg + self.gain = options.gain + self.decln = options.decln + + # Set initial values for datalogging timed-output + self.continuum_then = time.time() + self.spectral_then = time.time() + # build the graph + # + # If SETI mode, we always run at maximum USRP decimation + # + if (self.setimode): + options.decim = 256 + self.u = usrp.source_c(decim_rate=options.decim) self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec)) self.cardtype = self.u.daughterboard_id(0) # Set initial declination self.decln = options.decln - # Turn off interference filter by default - self.use_interfilt = options.interfilt - # determine the daughterboard subdevice we're using self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec) input_rate = self.u.adc_freq() / self.u.decim_rate() - tpstr="calib_"+options.cfunc+"_total_power" - sstr="calib_"+options.cfunc+"_fft" - self.tpcfunc=eval(tpstr) - self.scfunc=eval(sstr) - # # Set prefix for data files # self.prefix = options.prefix - calib_set_prefix(self.prefix) + + # + # The lower this number, the fewer sample frames are dropped + # in computing the FFT. A sampled approach is taken to + # computing the FFT of the incoming data, which reduces + # sensitivity. Increasing sensitivity inreases CPU loading. + # + self.fft_rate = options.fft_rate + + self.fft_size = int(options.fft_size) + + # This buffer is used to remember the most-recent FFT display + # values. Used later by self.write_spectral_data() to write + # spectral data to datalogging files, and by the SETI analysis + # function. + # + self.fft_outbuf = Numeric.zeros(self.fft_size, Numeric.Float64) + + # + # If SETI mode, only look at seti_fft_bandwidth (currently 12.5Khz) + # at a time. + # + if (self.setimode): + self.fft_input_rate = self.seti_fft_bandwidth + + # + # Build a decimating bandpass filter + # + self.fft_input_taps = gr.firdes.complex_band_pass (1.0, + input_rate, + -(int(self.fft_input_rate/2)), int(self.fft_input_rate/2), 200, + gr.firdes.WIN_HAMMING, 0) + + # + # Compute required decimation factor + # + decimation = int(input_rate/self.fft_input_rate) + self.fft_bandpass = gr.fir_filter_ccc (decimation, + self.fft_input_taps) + else: + self.fft_input_rate = input_rate # Set up FFT display - self.scope = ra_fftsink.ra_fft_sink_c (self, panel, - fft_size=int(options.fft_size), sample_rate=input_rate, - fft_rate=8, title="Spectral", - cfunc=self.scfunc, xydfunc=self.xydfunc, interfunc=self.interference) + if self.waterfall == False: + self.scope = ra_fftsink.ra_fft_sink_c (self, panel, + fft_size=int(self.fft_size), sample_rate=self.fft_input_rate, + fft_rate=int(self.fft_rate), title="Spectral", + ofunc=self.fft_outfunc, xydfunc=self.xydfunc) + else: + self.scope = ra_waterfallsink.ra_waterfallsink_c (self, panel, + fft_size=int(self.fft_size), sample_rate=self.fft_input_rate, + fft_rate=int(self.fft_rate), title="Spectral", ofunc=self.fft_outfunc, xydfunc=self.xydfunc_waterfall) # Set up ephemeris data self.locality = ephem.Observer() self.locality.long = str(options.longitude) self.locality.lat = str(options.latitude) + # We make notes about Sunset/Sunrise in Continuum log files + self.sun = ephem.Sun() + self.sunstate = "??" # Set up stripchart display self.stripsize = int(options.stripsize) - self.chart = ra_stripchartsink.stripchart_sink_f (self, panel, - stripsize=self.stripsize, - title="Continuum", - xlabel="LMST Offset (Seconds)", - scaling=1.0, ylabel=options.ylabel, - divbase=options.divbase, cfunc=self.tpcfunc) + if self.setimode == False: + self.chart = ra_stripchartsink.stripchart_sink_f (self, panel, + stripsize=self.stripsize, + title="Continuum", + xlabel="LMST Offset (Seconds)", + scaling=1.0, ylabel=options.ylabel, + divbase=options.divbase) # Set center frequency self.centerfreq = options.freq @@ -137,14 +264,6 @@ class app_flow_graph(stdgui.gui_flow_graph): self.bw = input_rate - # - # Produce a default interference map - # May not actually get used, unless --interfilt was specified - # - self.intmap = Numeric.zeros(256,Numeric.Complex64) - for i in range(0,len(self.intmap)): - self.intmap[i] = complex(1.0, 0.0) - # We setup the first two integrators to produce a fixed integration # Down to 1Hz, with output at 1 samples/sec N = input_rate/5000 @@ -175,90 +294,102 @@ class app_flow_graph(stdgui.gui_flow_graph): # Call constructors for receive chains # - # - # This is the interference-zapping filter - # - # The GUI is used to set/clear inteference zones in - # the filter. The non-interfering zones are set to - # 1.0. - # - if 0: - self.interfilt = gr.fft_filter_ccc(1,self.intmap) - tmp = FFT.inverse_fft(self.intmap) - self.interfilt.set_taps(tmp) - - # The three integrators--two FIR filters, and an IIR final filter - self.integrator1 = gr.fir_filter_fff (N, tapsN) - self.integrator2 = gr.fir_filter_fff (M, tapsM) - self.integrator3 = gr.single_pole_iir_filter_ff(1.0) - - # Split complex USRP stream into a pair of floats - self.splitter = gr.complex_to_float (1); - self.toshort = gr.float_to_short(); - - # I squarer (detector) - self.multI = gr.multiply_ff(); - - # Q squarer (detector) - self.multQ = gr.multiply_ff(); - - # Adding squared I and Q to produce instantaneous signal power - self.adder = gr.add_ff(); + if self.setimode == False: + # The three integrators--two FIR filters, and an IIR final filter + self.integrator1 = gr.fir_filter_fff (N, tapsN) + self.integrator2 = gr.fir_filter_fff (M, tapsM) + self.integrator3 = gr.single_pole_iir_filter_ff(1.0) + + # Split complex USRP stream into a pair of floats + self.splitter = gr.complex_to_float (1); + + # I squarer (detector) + self.multI = gr.multiply_ff(); + + # Q squarer (detector) + self.multQ = gr.multiply_ff(); + + # Adding squared I and Q to produce instantaneous signal power + self.adder = gr.add_ff(); + + # Signal probe + self.probe = gr.probe_signal_f(); + + # + # Continuum calibration stuff + # + self.cal_mult = gr.multiply_const_ff(self.calib_coeff); + self.cal_offs = gr.add_const_ff(self.calib_offset); # # Start connecting configured modules in the receive chain # - # Connect interference-filtered USRP input to selected scope function - if self.use_interfilt == True: - self.connect(self.u, self.interfilt, self.scope) - - # Connect interference-filtered USRP to a complex->float splitter - self.connect(self.interfilt, self.splitter) - - else: + # The scope--handle SETI mode + if (self.setimode == False): self.connect(self.u, self.scope) - self.connect(self.u, self.splitter) - - # Connect splitter outputs to multipliers - # First do I^2 - self.connect((self.splitter, 0), (self.multI,0)) - self.connect((self.splitter, 0), (self.multI,1)) - - # Then do Q^2 - self.connect((self.splitter, 1), (self.multQ,0)) - self.connect((self.splitter, 1), (self.multQ,1)) - - # Then sum the squares - self.connect(self.multI, (self.adder,0)) - self.connect(self.multQ, (self.adder,1)) - - # Connect adder output to three-stages of FIR integrator - self.connect(self.adder, self.integrator1, - self.integrator2, self.integrator3, self.chart) + else: + self.connect(self.u, self.fft_bandpass, self.scope) + if self.setimode == False: + # + # The head of the continuum chain + # + self.connect(self.u, self.splitter) + + # Connect splitter outputs to multipliers + # First do I^2 + self.connect((self.splitter, 0), (self.multI,0)) + self.connect((self.splitter, 0), (self.multI,1)) + + # Then do Q^2 + self.connect((self.splitter, 1), (self.multQ,0)) + self.connect((self.splitter, 1), (self.multQ,1)) + + # Then sum the squares + self.connect(self.multI, (self.adder,0)) + self.connect(self.multQ, (self.adder,1)) + + # Connect adder output to two-stages of FIR integrator + # followed by a single stage IIR integrator, and + # the calibrator + self.connect(self.adder, self.integrator1, + self.integrator2, self.integrator3, self.cal_mult, + self.cal_offs, self.chart) + + # Connect calibrator to probe + # SPECIAL NOTE: I'm setting the ground work here + # for completely changing the way local_calibrator + # works, including removing some horrible kludges for + # recording data. + # But for now, self.probe() will be used to display the + # current instantaneous integrated detector value + self.connect(self.cal_offs, self.probe) self._build_gui(vbox) # Make GUI agree with command-line - self.myform['integration'].set_value(int(options.integ)) + self.integ = options.integ + if self.setimode == False: + self.myform['integration'].set_value(int(options.integ)) self.myform['average'].set_value(int(options.avg)) - # Make integrator agree with command line - self.set_integration(int(options.integ)) + if self.setimode == False: + # Make integrator agree with command line + self.set_integration(int(options.integ)) + + self.avg_alpha = options.avg # Make spectral averager agree with command line if options.avg != 1.0: self.scope.set_avg_alpha(float(1.0/options.avg)) - calib_set_avg_alpha(float(options.avg)) self.scope.set_average(True) - - # Set division size - self.chart.set_y_per_div(options.division) - - # Set reference(MAX) level - self.chart.set_ref_level(options.reflevel) + if self.setimode == False: + # Set division size + self.chart.set_y_per_div(options.division) + # Set reference(MAX) level + self.chart.set_ref_level(options.reflevel) # set initial values @@ -278,24 +409,23 @@ class app_flow_graph(stdgui.gui_flow_graph): if not(self.set_freq(options.freq)): self._set_status_msg("Failed to set initial frequency") + # Set declination self.set_decln (self.decln) - calib_set_decln (self.decln) + + # RF hardware information self.myform['decim'].set_value(self.u.decim_rate()) self.myform['fs@usb'].set_value(self.u.adc_freq() / self.u.decim_rate()) self.myform['dbname'].set_value(self.subdev.name()) - # Make sure calibrator knows what our bandwidth is - calib_set_bw(self.u.adc_freq() / self.u.decim_rate()) - # Set analog baseband filtering, if DBS_RX if self.cardtype == usrp_dbid.DBS_RX: - self.subdev.set_bw((self.u.adc_freq() / self.u.decim_rate())/2) - - # Tell calibrator our declination as well - calib_set_decln(self.decln) + lbw = (self.u.adc_freq() / self.u.decim_rate()) / 2 + if lbw < 1.0e6: + lbw = 1.0e6 + self.subdev.set_bw(lbw) - # Start the timer for the LMST display + # Start the timer for the LMST display and datalogging self.lmst_timer.Start(1000) @@ -305,6 +435,17 @@ class app_flow_graph(stdgui.gui_flow_graph): def _build_gui(self, vbox): def _form_set_freq(kv): + # Adjust current SETI frequency, and limits + self.setifreq_lower = kv['freq'] - (self.seti_freq_range/2) + self.setifreq_current = kv['freq'] + self.setifreq_upper = kv['freq'] + (self.seti_freq_range/2) + + # Reset SETI analysis timer + self.seti_then = time.time() + # Zero-out hits array when changing frequency + self.hits_array[:,:] = 0.0 + self.hit_intensities[:,:] = -60.0 + return self.set_freq(kv['freq']) def _form_set_decln(kv): @@ -313,8 +454,9 @@ class app_flow_graph(stdgui.gui_flow_graph): # Position the FFT display vbox.Add(self.scope.win, 15, wx.EXPAND) - # Position the Total-power stripchart - vbox.Add(self.chart.win, 15, wx.EXPAND) + if self.setimode == False: + # Position the Total-power stripchart + vbox.Add(self.chart.win, 15, wx.EXPAND) # add control area at the bottom self.myform = myform = form.form() @@ -344,26 +486,32 @@ class app_flow_graph(stdgui.gui_flow_graph): vbox2.Add((4,0), 0, 0) myform['average'] = form.slider_field(parent=self.panel, sizer=vbox2, - label="Spectral Averaging (FFT frames)", weight=1, min=1, max=2000, callback=self.set_averaging) - + label="Spectral Averaging (FFT frames)", weight=1, min=1, max=3000, callback=self.set_averaging) + + # Set up scan control button when in SETI mode + if (self.setimode == True): + # SETI scanning control + buttonbox = wx.BoxSizer(wx.HORIZONTAL) + self.scan_control = form.button_with_callback(self.panel, + label="Scan: On ", + callback=self.toggle_scanning) + + buttonbox.Add(self.scan_control, 0, wx.CENTER) + vbox2.Add(buttonbox, 0, wx.CENTER) vbox2.Add((4,0), 0, 0) - myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2, - label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration) + if self.setimode == False: + myform['integration'] = form.slider_field(parent=self.panel, sizer=vbox2, + label="Continuum Integration Time (sec)", weight=1, min=1, max=180, callback=self.set_integration) + + vbox2.Add((4,0), 0, 0) - vbox2.Add((4,0), 0, 0) myform['decln'] = form.float_field( parent=self.panel, sizer=vbox2, label="Current Declination", weight=1, callback=myform.check_input_and_call(_form_set_decln)) vbox2.Add((4,0), 0, 0) buttonbox = wx.BoxSizer(wx.HORIZONTAL) - if self.use_interfilt == True: - self.doit = form.button_with_callback(self.panel, - label="Clear Interference List", - callback=self.clear_interferers) - if self.use_interfilt == True: - buttonbox.Add(self.doit, 0, wx.CENTER) vbox.Add(buttonbox, 0, wx.CENTER) hbox.Add(vbox1, 0, 0) hbox.Add(vbox2, wx.ALIGN_RIGHT, 0) @@ -372,7 +520,7 @@ class app_flow_graph(stdgui.gui_flow_graph): self._build_subpanel(vbox) self.lmst_timer = wx.PyTimer(self.lmst_timeout) - self.lmst_timeout() + #self.lmst_timeout() def _build_subpanel(self, vbox_arg): @@ -450,10 +598,6 @@ class app_flow_graph(stdgui.gui_flow_graph): self.centerfreq = target_freq self.observing -= delta self.scope.set_baseband_freq (self.observing) - calib_set_freq(self.observing) - - # Clear interference list - self.clear_interferers() self.myform['baseband'].set_value(r.baseband_freq) self.myform['ddc'].set_value(r.dxc_freq) @@ -465,39 +609,313 @@ class app_flow_graph(stdgui.gui_flow_graph): def set_decln(self, dec): self.decln = dec self.myform['decln'].set_value(dec) # update displayed value - calib_set_decln(dec) def set_gain(self, gain): self.myform['gain'].set_value(gain) # update displayed value self.subdev.set_gain(gain) - - # - # Make sure calibrator knows our gain setting - # - calib_set_gain(gain) + self.gain = gain def set_averaging(self, avval): self.myform['average'].set_value(avval) self.scope.set_avg_alpha(1.0/(avval)) - calib_set_avg_alpha(avval) self.scope.set_average(True) + self.avg_alpha = avval def set_integration(self, integval): - self.integrator3.set_taps(1.0/integval) + if self.setimode == False: + self.integrator3.set_taps(1.0/integval) self.myform['integration'].set_value(integval) + self.integ = integval + + # + # Timeout function + # Used to update LMST display, as well as current + # continuum value + # + # We also write external data-logging files here + # + def lmst_timeout(self): + self.locality.date = ephem.now() + if self.setimode == False: + x = self.probe.level() + sidtime = self.locality.sidereal_time() + # LMST + s = str(ephem.hours(sidtime)) + " " + self.sunstate + # Continuum detector value + if self.setimode == False: + sx = "%7.4f" % x + s = s + "\nDet: " + str(sx) + else: + sx = "%2d" % self.hitcounter + sy = "%3.1f-%3.1f" % (self.CHIRP_LOWER, self.CHIRP_UPPER) + s = s + "\nHits: " + str(sx) + "\nCh lim: " + str(sy) + self.myform['lmst_high'].set_value(s) + + # + # Write data out to recording files + # + if self.setimode == False: + self.write_continuum_data(x,sidtime) + self.write_spectral_data(self.fft_outbuf,sidtime) + + else: + self.seti_analysis(self.fft_outbuf,sidtime) + now = time.time() + if ((self.scanning == True) and ((now - self.seti_then) > self.setifreq_timer)): + self.seti_then = now + self.setifreq_current = self.setifreq_current + self.fft_input_rate + if (self.setifreq_current > self.setifreq_upper): + self.setifreq_current = self.setifreq_lower + self.set_freq(self.setifreq_current) + # Make sure we zero-out the hits array when changing + # frequency. + self.hits_array[:,:] = 0.0 + self.hit_intensities[:,:] = 0.0 + + def fft_outfunc(self,data,l): + self.fft_outbuf=data + + def write_continuum_data(self,data,sidtime): + + # Create localtime structure for producing filename + foo = time.localtime() + pfx = self.prefix + filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year, + foo.tm_mon, foo.tm_mday, foo.tm_hour) + + # Open the data file, appending + continuum_file = open (filenamestr+".tpdat","a") + + flt = "%6.3f" % data + inter = self.decln + integ = self.integ + fc = self.observing + fc = fc / 1000000 + bw = self.bw + bw = bw / 1000000 + ga = self.gain + + now = time.time() + + # + # If time to write full header info (saves storage this way) + # + if (now - self.continuum_then > 20): + self.sun.compute(self.locality) + enow = ephem.now() + sun_insky = "Down" + self.sunstate = "Dn" + if ((self.sun.rise_time < enow) and (enow < self.sun.set_time)): + sun_insky = "Up" + self.sunstate = "Up" + self.continuum_then = now + + continuum_file.write(str(ephem.hours(sidtime))+" "+flt+" Dn="+str(inter)+",") + continuum_file.write("Ti="+str(integ)+",Fc="+str(fc)+",Bw="+str(bw)) + continuum_file.write(",Ga="+str(ga)+",Sun="+str(sun_insky)+"\n") + else: + continuum_file.write(str(ephem.hours(sidtime))+" "+flt+"\n") + + continuum_file.close() + return(data) + + def write_spectral_data(self,data,sidtime): + + now = time.time() + delta = 10 + + # If time to write out spectral data + # We don't write this out every time, in order to + # save disk space. Since the spectral data are + # typically heavily averaged, writing this data + # "once in a while" is OK. + # + if (now - self.spectral_then >= delta): + self.spectral_then = now + + # Get localtime structure to make filename from + foo = time.localtime() + + pfx = self.prefix + filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year, + foo.tm_mon, foo.tm_mday, foo.tm_hour) + + # Open the file + spectral_file = open (filenamestr+".sdat","a") + + # Setup data fields to be written + r = data + inter = self.decln + fc = self.observing + fc = fc / 1000000 + bw = self.bw + bw = bw / 1000000 + av = self.avg_alpha + + # Write those fields + spectral_file.write("data:"+str(ephem.hours(sidtime))+" Dn="+str(inter)+",Fc="+str(fc)+",Bw="+str(bw)+",Av="+str(av)) + spectral_file.write(" "+str(r)+"\n") + spectral_file.close() + return(data) + + return(data) + + def seti_analysis(self,fftbuf,sidtime): + l = len(fftbuf) + x = 0 + hits = [] + hit_intensities = [] + if self.seticounter < self.setitimer: + self.seticounter = self.seticounter + 1 + return + else: + self.seticounter = 0 + + # Run through FFT output buffer, computing standard deviation (Sigma) + avg = 0 + # First compute average + for i in range(0,l): + avg = avg + fftbuf[i] + avg = avg / l + + sigma = 0.0 + # Then compute standard deviation (Sigma) + for i in range(0,l): + d = fftbuf[i] - avg + sigma = sigma + (d*d) + + sigma = Numeric.sqrt(sigma/l) + + # + # Snarfle through the FFT output buffer again, looking for + # outlying data points + + start_f = self.observing - (self.fft_input_rate/2) + current_f = start_f + f_incr = self.fft_input_rate / l + l = len(fftbuf) + hit = -1 + + # -nyquist to DC + for i in range(l/2,l): + # + # If current FFT buffer has an item that exceeds the specified + # sigma + # + if ((fftbuf[i] - avg) > (self.setik * sigma)): + hits.append(current_f) + hit_intensities.append(fftbuf[i]) + current_f = current_f + f_incr + + # DC to nyquist + for i in range(0,l/2): + # + # If current FFT buffer has an item that exceeds the specified + # sigma + # + if ((fftbuf[i] - avg) > (self.setik * sigma)): + hits.append(current_f) + hit_intensities.append(fftbuf[i]) + current_f = current_f + f_incr + + # No hits + if (len(hits) <= 0): + return # - # Make sure calibrator knows our integration time + # OK, so we have some hits in the FFT buffer + # They'll have a rather substantial gauntlet to run before + # being declared a real "hit" # - calib_set_integ(integval) - def lmst_timeout(self): - self.locality.date = ephem.now() - sidtime = self.locality.sidereal_time() - self.myform['lmst_high'].set_value(str(ephem.hours(sidtime))) + # Weed out buffers with an excessive number of strong signals + if (len(hits) > self.nhits): + return + + # Weed out FFT buffers with apparent multiple narrowband signals + # separated significantly in frequency. This means that a + # single signal spanning multiple bins is OK, but a buffer that + # has multiple, apparently-separate, signals isn't OK. + # + last = hits[0] + for i in range(1,len(hits)): + if ((hits[i] - last) > (f_incr*2.0)): + return + last = hits[i] + + # + # Run through all three hit buffers, computing difference between + # frequencies found there, if they're all within the chirp limits + # declare a good hit + # + good_hit = 0 + good_hit = False + for i in range(0,min(len(hits),len(self.hits_array[:,0]))): + f_d1 = abs(self.hits_array[i,0] - hits[i]) + f_d2 = abs(self.hits_array[i,1] - self.hits_array[i,0]) + f_d3 = abs(self.hits_array[i,2] - self.hits_array[i,1]) + if (self.seti_isahit ([f_d1, f_d2, f_d3])): + good_hit = True + self.hitcounter = self.hitcounter + 1 + break + + + # Save 'n shuffle hits + for i in range(self.nhitlines,1): + self.hits_array[:,i] = self.hits_array[:,i-1] + self.hit_intensities[:,i] = self.hit_intensities[:,i-1] + + for i in range(0,len(hits)): + self.hits_array[i,0] = hits[i] + self.hit_intensities[i,0] = hit_intensities[i] + + # Finally, write the hits/intensities buffer + if (good_hit): + self.write_hits(sidtime) + + return + + def seti_isahit(self,fdiffs): + truecount = 0 + + for i in range(0,len(fdiffs)): + if (fdiffs[i] >= self.CHIRP_LOWER and fdiffs[i] <= self.CHIRP_UPPER): + truecount = truecount + 1 + + if truecount == len(fdiffs): + return (True) + else: + return (False) + + def write_hits(self,sidtime): + # Create localtime structure for producing filename + foo = time.localtime() + pfx = self.prefix + filenamestr = "%s/%04d%02d%02d%02d" % (pfx, foo.tm_year, + foo.tm_mon, foo.tm_mday, foo.tm_hour) + + # Open the data file, appending + hits_file = open (filenamestr+".seti","a") + + # Write sidtime first + hits_file.write(str(ephem.hours(sidtime))+" "+str(self.decln)+" ") + + # + # Then write the hits/hit intensities buffers with enough + # "syntax" to allow parsing by external (not yet written!) + # "stuff". + # + for i in range(0,self.nhitlines): + hits_file.write(" ") + for j in range(0,self.nhits): + hits_file.write(str(self.hits_array[j,i])+":") + hits_file.write(str(self.hit_intensities[j,i])+",") + hits_file.write("\n") + hits_file.close() + return def xydfunc(self,xyv): - magn = int(log10(self.observing)) + magn = int(Numeric.log10(self.observing)) if (magn == 6 or magn == 7 or magn == 8): magn = 6 dfreq = xyv[0] * pow(10.0,magn) @@ -514,51 +932,12 @@ class app_flow_graph(stdgui.gui_flow_graph): s2 = "\n%.3fkm/s" % vs self.myform['spec_data'].set_value(s+s2) - def interference(self,x): - if self.use_interfilt == False: - return - magn = int(log10(self.observing)) - dfreq = x * pow(10.0,magn) - delta = dfreq - self.observing - fincr = self.bw / len(self.intmap) - l = len(self.intmap) - if delta > 0: - offset = delta/fincr - else: - offset = (l) - int((abs(delta)/fincr)) - - offset = int(offset) - - if offset >= len(self.intmap) or offset < 0: - print "interference offset is invalid--", offset - return - - # - # Zero out the region around the selected interferer - # - self.intmap[offset-2] = complex (0.5, 0.0) - self.intmap[offset-1] = complex (0.25, 0.0) - self.intmap[offset] = complex (0.0, 0.0) - self.intmap[offset+1] = complex(0.25, 0.0) - self.intmap[offset+2] = complex(0.5, 0.0) - - # - # Set new taps - # - tmp = FFT.inverse_fft(self.intmap) - self.interfilt.set_taps(tmp) - - def clear_interf(self): - self.clear_interferers() - - def clear_interferers(self): - for i in range(0,len(self.intmap)): - self.intmap[i] = complex(1.0,0.0) - tmp = FFT.inverse_fft(self.intmap) - if self.use_interfilt == True: - self.interfilt.set_taps(tmp) - - + def xydfunc_waterfall(self,pos): + lower = self.observing - (self.seti_fft_bandwidth / 2) + upper = self.observing + (self.seti_fft_bandwidth / 2) + binwidth = self.seti_fft_bandwidth / 1024 + s = "%.6fMHz" % ((lower + (pos.x*binwidth)) / 1.0e6) + self.myform['spec_data'].set_value(s) def toggle_cal(self): if (self.calstate == True): @@ -577,8 +956,19 @@ class app_flow_graph(stdgui.gui_flow_graph): else: self.annotate_state = True self.annotation.SetLabel("Annotation: On") - calib_set_interesting(self.annotate_state) - + # + # Turn scanning on/off + # Called-back by "Recording" button + # + def toggle_scanning(self): + # Current scanning? Flip state + if (self.scanning == True): + self.scanning = False + self.scan_control.SetLabel("Scan: Off") + # Not scanning + else: + self.scanning = True + self.scan_control.SetLabel("Scan: On ") def main (): app = stdgui.stdapp(app_flow_graph, "RADIO ASTRONOMY SPECTRAL/CONTINUUM RECEIVER: $Revision$", nstatus=1)