From 74e4486598b3ea40e13312c355814372930e5b62 Mon Sep 17 00:00:00 2001 From: mleech Date: Tue, 10 Oct 2006 02:40:19 +0000 Subject: Improved SETI mode--added 1Mhz swath scanning, and reduced SETI mode FFT bandwidth to 12.5Khz, allowing sub-1Hz resolution. git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@3758 221aa14e-8319-0410-a670-987f0aec2ac5 --- gr-radio-astronomy/src/python/ra_waterfallsink.py | 4 +- gr-radio-astronomy/src/python/usrp_ra_receiver.py | 160 +++++++++++++++------- 2 files changed, 112 insertions(+), 52 deletions(-) (limited to 'gr-radio-astronomy') diff --git a/gr-radio-astronomy/src/python/ra_waterfallsink.py b/gr-radio-astronomy/src/python/ra_waterfallsink.py index b6c63ac25..f2244773f 100755 --- a/gr-radio-astronomy/src/python/ra_waterfallsink.py +++ b/gr-radio-astronomy/src/python/ra_waterfallsink.py @@ -193,7 +193,7 @@ class ra_waterfall_window (wx.Panel): wx.Panel.__init__(self, parent, id, pos, size, style, name) self.fftsink = fftsink - self.bm = wx.EmptyBitmap(self.fftsink.fft_size, 300, -1) + self.bm = wx.EmptyBitmap(1024, 300, -1) self.scale_factor = self.fftsink.scaling @@ -295,7 +295,7 @@ class ra_waterfall_window (wx.Panel): dc1 = wx.MemoryDC() dc1.SelectObject(self.bm) - dc1.Blit(0,1,self.fftsink.fft_size,300,dc1,0,0,wx.COPY,False,-1,-1) + dc1.Blit(0,1,1024,300,dc1,0,0,wx.COPY,False,-1,-1) x = max(abs(self.fftsink.sample_rate), abs(self.fftsink.baseband_freq)) if x >= 1e9: diff --git a/gr-radio-astronomy/src/python/usrp_ra_receiver.py b/gr-radio-astronomy/src/python/usrp_ra_receiver.py index cdc5b6e19..a4ffccf90 100755 --- a/gr-radio-astronomy/src/python/usrp_ra_receiver.py +++ b/gr-radio-astronomy/src/python/usrp_ra_receiver.py @@ -81,7 +81,6 @@ class app_flow_graph(stdgui.gui_flow_graph): parser.add_option("-Q", "--calib_eqn", default="x = x * 1.0", help="Calibration equation") 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("-T", "--setitimer", type="eng_float", default=15.0, help="Timer for computing doppler chirp for SETI analysis") parser.add_option("-K", "--setik", type="eng_float", default=1.5, help="K value for SETI analysis") (options, args) = parser.parse_args() if len(args) != 0: @@ -89,14 +88,29 @@ class app_flow_graph(stdgui.gui_flow_graph): sys.exit(1) self.show_debug_info = True + self.waterfall = options.waterfall + + # SETI mode stuff self.setimode = options.setimode self.seticounter = 0 - self.setitimer = int(options.setitimer) self.setik = options.setik + self.seti_fft_bandwidth = 12500 + binwidth = self.seti_fft_bandwidth / options.fft_size + upper_limit = binwidth / 0.10 + self.setitimer = int(upper_limit * 2.00) + self.CHIRP_LOWER = 0.10 * self.setitimer + self.CHIRP_UPPER = 0.25 * self.setitimer self.hitcounter = 0 - self.CHIRP_LOWER = 15 - self.CHIRP_UPPER = 50 + # We scan through 1Mhz of bandwidth around the chosen center freq + self.seti_freq_range = 1.0e6 + self.setifreq_lower = options.freq - (self.seti_freq_range/2) + self.setifreq_current = options.freq + self.setifreq_upper = options.freq + (self.seti_freq_range/2) + self.setifreq_timer = self.setitimer * 10 + self.seti_then = time.time() + self.nhits = 10 + self.hits_array = Numeric.zeros((self.nhits,3), Numeric.Float64) # Calibration coefficient and offset self.calib_coeff = options.calib_coeff @@ -116,6 +130,12 @@ class app_flow_graph(stdgui.gui_flow_graph): # 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) @@ -146,18 +166,38 @@ class app_flow_graph(stdgui.gui_flow_graph): # values. Used later by self.write_spectral_data() to write # spectral data to datalogging files. self.fft_outbuf = Numeric.zeros(options.fft_size, Numeric.Float64) - self.old_hits = Numeric.zeros(10, Numeric.Float64) - self.older_hits = Numeric.zeros(10, Numeric.Float64) + # + # If SETI mode, only look at a input_rate/8 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 if self.waterfall == False: self.scope = ra_fftsink.ra_fft_sink_c (self, panel, - fft_size=int(self.fft_size), sample_rate=input_rate, + 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=input_rate, + 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) # Set up ephemeris data @@ -247,7 +287,16 @@ class app_flow_graph(stdgui.gui_flow_graph): # # Start connecting configured modules in the receive chain # - self.connect(self.u, self.scope) + + # The scope--handle SETI mode + if (self.setimode == False): + self.connect(self.u, self.scope) + else: + self.connect(self.u, self.fft_bandpass, self.scope) + + # + # The head of the continuum chain + # self.connect(self.u, self.splitter) # Connect splitter outputs to multipliers @@ -346,6 +395,11 @@ class app_flow_graph(stdgui.gui_flow_graph): def _build_gui(self, vbox): def _form_set_freq(kv): + 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) + self.seti_then = time.time() + self.hits_array[:,:] = 0.0 return self.set_freq(kv['freq']) def _form_set_decln(kv): @@ -407,7 +461,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): @@ -530,18 +584,27 @@ class app_flow_graph(stdgui.gui_flow_graph): sx = "%7.4f" % x s = s + "\nDet: " + str(sx) sx = "%2d" % self.hitcounter - sy = "%2d" % self.CHIRP_LOWER - s = s + "\nH: " + str(sx) + " Cl: " + str(sy) + sy = "%3.1f/%3.1f" % (self.CHIRP_LOWER, self.CHIRP_UPPER) + s = s + "\nH: " + str(sx) + " C: " + str(sy) self.myform['lmst_high'].set_value(s) # # Write data out to recording files # - self.write_continuum_data(x,sidtime) - self.write_spectral_data(self.fft_outbuf,sidtime) + if self.setimode == False: + self.write_continuum_data(x,sidtime) + self.write_spectral_data(self.fft_outbuf,sidtime) if self.setimode == True: self.seti_analysis(self.fft_outbuf,sidtime) + now = time.time() + if ((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) + self.hits_array[:,:] = 0.0 def fft_outfunc(self,data,l): self.fft_outbuf=data @@ -653,9 +716,9 @@ class app_flow_graph(stdgui.gui_flow_graph): # Snarfle through the FFT output buffer again, looking for # outlying data points - start_f = self.observing - (self.bw/2) + start_f = self.observing - (self.fft_input_rate/2) current_f = start_f - f_incr = self.bw / l + f_incr = self.fft_input_rate / l l = len(fftbuf) hit = -1 @@ -682,24 +745,9 @@ class app_flow_graph(stdgui.gui_flow_graph): if (len(hits) <= 0): return - if (len(hits) > len(self.old_hits)*2): + if (len(hits) > self.nhits): return - # - # - # Calculate chirp limits from first principles - # - # - earth_diam = 12500 - earth_circ = earth_diam * 3.14159 - surface_speed = earth_circ / 24 - surface_speed = surface_speed / 3600 - - c1 = (surface_speed/2) / 299792.0 - c2 = (surface_speed*5) / 299792.0 - - self.CHIRP_LOWER = (c1 * self.observing) * self.setitimer - self.CHIRP_UPPER = (c2 * self.observing) * self.setitimer # # Run through all three hit buffers, computing difference between @@ -707,28 +755,40 @@ class app_flow_graph(stdgui.gui_flow_graph): # declare a good hit # good_hit = 0 - for i in range(0,min(len(hits),len(self.old_hits))): - f_diff = abs(self.old_hits[i] - hits[i]) - f_diff2 = abs(self.older_hits[i] - self.old_hits[i]) - # If frequency difference is within range, we have a hit - if (f_diff >= self.CHIRP_LOWER and f_diff <= self.CHIRP_UPPER): - if (f_diff2 >= self.CHIRP_LOWER and f_diff <= self.CHIRP_UPPER): - good_hit = 1 - self.hitcounter = self.hitcounter + 1 - break - - if (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 + + if (good_hit): self.write_hits(hits,sidtime) - # Save old hits - for i in range(0,len(self.older_hits)): - self.older_hits[i] = self.old_hits[i] - self.old_hits[i] = 0 - for i in range(0,min(len(hits),len(self.old_hits))): - self.old_hits[i] = hits[i] + # Save 'n shuffle hits + self.hits_array[:,2] = self.hits_array[:,1] + self.hits_array[:,1] = self.hits_array[:,0] + + for i in range(0,len(hits)): + self.hits_array[i,0] = hits[i] 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,hits,sidtime): # Create localtime structure for producing filename foo = time.localtime() @@ -738,7 +798,7 @@ class app_flow_graph(stdgui.gui_flow_graph): # Open the data file, appending hits_file = open (filenamestr+".seti","a") - hits_file.write(str(ephem.hours(sidtime))+" "+str(hits)+"\n") + hits_file.write(str(ephem.hours(sidtime))+" "+str(self.decln)+" "+str(hits)+"\n") hits_file.close() return -- cgit