My VOR receiver is almost there, but not quite... My code is below. I've demodulated the first 30 Hz signal from AM, and the second from the FM subcarrier. Both appear pretty solid on the O-scope (except for a DC offset on the AM signal), and visually I can see a constant difference in phase on the scope.
However, my attempts to get relative phase aren't working. First I tried exactly what was listed below, and I got either 0 or pi as a phase difference. Then I realized that the demodulators were returning a real signal, and type converting it back to complex was leaving off the imaginary portion. So I used a Hilbert transform to re-create the imaginary portion, and then fed those to the phase difference part of the graph. Unfortunately, I'm not getting a constant return, it's basically 'spinning', the returned value is constantly changing in a cyclical manner. I tried another method to determine the relative phase, but it didn't do anything different. I figure it's either got to do with the DC offset, and I'm not sure the best way to get rid of this because my signal is so low frequency, or I'm using the Hilbert filter wrong. Below is the source code, it's not cleaned up yet, but it's functional. To use you'd need to find a VOR station and get near it. They really aren't designed for ground use, so you may need to get close. But like I said, the signals look clean on the scope, so I'm baffled. Thoughts? Geof ------------------------ > You can get the phase difference by multiplying one signal by the > complex conjugate of the other, then take the arg (arctan) of that. > There are blocks for each of those operations. > > sig0 = ... > sig1 = .. > > conj = gr.conjugate_cc() > mult = gr.multiply_cc() > arg = gr.complex_to_arg() > > dst = ... # stream of float radians > > connect(sig0, conj, mult) > connect(sig1, (mult, 1)) > connect(mult, arg) > connect(arg, dst) -------------------- #!/usr/bin/env python # # Copyright 2005,2006,2007,2008 Free Software Foundation, Inc. # # This file is part of GNU Radio # # GNU Radio is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 3, or (at your option) # any later version. # # GNU Radio is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # 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., 51 Franklin Street, # Boston, MA 02110-1301, USA. # # This module tunes to VOR receivers and determines the bearing from the radio # to the receiver # # See http://www.doe.carleton.ca/courses/ELEC4504/nav_r7.pdf for a quick # VOR theory discussion from gnuradio import gr, gru, eng_notation, optfir from gnuradio import audio from gnuradio import usrp from gnuradio import blks2 from gnuradio.eng_option import eng_option from gnuradio.wxgui import slider, powermate, numbersink2 from gnuradio.wxgui import stdgui2, fftsink2, form, scopesink2 from optparse import OptionParser from usrpm import usrp_dbid import sys import math import wx def pick_subdevice(u): """ The user didn't specify a subdevice on the command line. Try for one of these, in order: BASIC_RX,TV_RX, BASIC_RX, whatever is on side A. @return a subdev_spec """ return usrp.pick_subdev(u, (usrp_dbid.BASIC_RX, usrp_dbid.LF_RX, usrp_dbid.TV_RX, usrp_dbid.TV_RX_REV_2, usrp_dbid.TV_RX_REV_3)) class vor_rx_block (stdgui2.std_top_block): def __init__(self,frame,panel,vbox,argv): stdgui2.std_top_block.__init__ (self,frame,panel,vbox,argv) parser=OptionParser(option_class=eng_option) parser.add_option("-R", "--rx-subdev-spec", type="subdev", default=None, help="select USRP Rx side A or B (default=A)") parser.add_option("-f", "--freq", type="eng_float", default=1008.0e3, help="set frequency to FREQ", metavar="FREQ") parser.add_option("-I", "--use-if-freq", action="store_true", default=False, help="use intermediate freq (compensates DC problems in quadrature boards)" ) parser.add_option("-g", "--gain", type="eng_float", default=None, help="set gain in dB (default is maximum)") parser.add_option("-V", "--volume", type="eng_float", default=None, help="set volume (default is midpoint)") parser.add_option("-O", "--audio-output", type="string", default="", help="pcm device name. E.g., hw:0,0 or surround51 or /dev/dsp") (options, args) = parser.parse_args() if len(args) != 0: parser.print_help() sys.exit(1) self.frame = frame self.panel = panel self.use_IF=options.use_if_freq if self.use_IF: self.IF_freq=64000.0 else: self.IF_freq=0.0 self.vol = 0 self.state = "FREQ" self.freq = 0 # build graph #TODO: add an AGC after the channel filter and before the AM_demod self.u = usrp.source_c() # usrp is data source adc_rate = self.u.adc_rate() # 64 MS/s usrp_decim = 250 self.u.set_decim_rate(usrp_decim) usrp_rate = adc_rate / usrp_decim # 256 kS/s chanfilt_decim = 8 demod_rate = usrp_rate / chanfilt_decim # 32 kHz audio_decimation = 128 audio_rate = demod_rate / audio_decimation # 500 Hz # This is the range of data which holds the carrier AM 30 Hz signal fm_decimation = 1 self.fm_rate = demod_rate / fm_decimation fm_audio_decimation = 128 self.fm_audio_rate = self.fm_rate / fm_audio_decimation# 500 Hz if options.rx_subdev_spec is None: options.rx_subdev_spec = pick_subdevice(self.u) self.u.set_mux(usrp.determine_rx_mux_value(self.u, options.rx_subdev_spec)) self.subdev = usrp.selected_subdev(self.u, options.rx_subdev_spec) print "Using RX d'board %s" % (self.subdev.side_and_name(),) # First bring in all the data we need to use, which includes the 30 Hz signal (AM modulated), # the 1020 Hz morse code identifier # and the 9960 Hz subcarrier which has a 480Hz width, so we'll cut off at 11000 Hz chan_filt_coeffs = optfir.low_pass (1, # gain usrp_rate, # sampling rate 11e3, # passband cutoff 12e3, # stopband cutoff 1.0, # passband ripple 60) # stopband attenuation if self.use_IF: # Turn If to baseband and filter. self.chan_filt = gr.freq_xlating_fir_filter_ccf (chanfilt_decim, chan_filt_coeffs, self.IF_freq, usrp_rate) else: self.chan_filt = gr.fir_filter_ccf (chanfilt_decim, chan_filt_coeffs) self.am_demod = gr.complex_to_mag() # Now we split off the 30 Hz signal after it's been demodulated audio_filt_coeffs = gr.firdes.band_pass (.5, # gain demod_rate, # sampling rate 20, # low cutoff 50, # high cutoff 150) # stopband cutoff self.audio_filt=gr.fir_filter_fff(audio_decimation,audio_filt_coeffs) # convert it back to complex with a Hilbert filter am_hilbert_coeffs = gr.firdes.hilbert (27) am_hilbert_filt = gr.fir_filter_fff(1, am_hilbert_coeffs) # that detected the main AM 30 Hz signal, now we need to decode the second signal # It's on a 9960Hz subcarrier, FM modulated second_chan_coeffs = gr.firdes.low_pass (1.0, # gain demod_rate, # sampling rate 1000, # low pass cutoff freq 1000, # width of trans. band gr.firdes.WIN_HANN) # filter type self.fm_filt = gr.freq_xlating_fir_filter_fcc(fm_decimation, # decimation rate second_chan_coeffs, # taps -9960, # frequency translation amount demod_rate) # input sample rate # Once the signal is brought down to baseband, then we can demodulate it self.fm_rx = gr.quadrature_demod_cf(1) # convert it back to complex with a Hilbert filter fm_hilbert_coeffs = gr.firdes.hilbert (27) fm_hilbert_filt = gr.fir_filter_fff(1, fm_hilbert_coeffs) # and once demodulated, then we can filter out everything but the 30Hz signal fm_audio_filt_coeffs = optfir.low_pass (10, # gain demod_rate / fm_decimation, # sampling rate 50, # passband cutoff 150, # stopband cutoff 0.1, # passband ripple 60) # stopband attenuation self.fm_audio_filt=gr.fir_filter_fff(fm_audio_decimation,fm_audio_filt_coeff s) # Now there are 2 30 Hz signal, one out of fm_audio_filt, another out of audio_filt. # The directional signal is found by comparing their phase self.fm_conv = gr.float_to_complex() self.am_conv = gr.float_to_complex() conj = gr.conjugate_cc() self.mult = gr.multiply_cc() self.arg = gr.complex_to_arg() self.fm_arg = gr.complex_to_arg() #self.sub = gr.sub_ff() # For the moment, no sound is output. However, we really should also take the chan_filt output, # and run it through a different filter/AM demod to get the 1020 Hz ID signal, and send that to the sound card # The code has been chopped for the moment, but we'll keep part of it here for future use self.volume_control = gr.multiply_const_ff(self.vol) # sound card as the final sink audio_sink = audio.sink (int (audio_rate), options.audio_output, False) # ok_to_block # now wire it all together # obtain our 2 real signals self.connect (self.u, self.chan_filt, self.am_demod, self.audio_filt) #, (self.volume_control, 0), (audio_sink,0)) self.connect (self.am_demod, self.fm_filt, self.fm_rx, self.fm_audio_filt) #, (audio_sink,1)) # convert the real signals back to complex w/ a hilbert filter self.connect(self.audio_filt, (self.am_conv,0)) self.connect(self.audio_filt, am_hilbert_filt, (self.am_conv,1)) self.connect(self.fm_audio_filt, (self.fm_conv,0)) self.connect(self.fm_audio_filt, fm_hilbert_filt, (self.fm_conv,1)) # now run them through the phase comparator self.connect(self.am_conv, conj, self.mult) self.connect(self.fm_conv, (self.mult, 1)) self.connect(self.mult, self.arg) # alternative phase comparator (comment out the mult/conj definitions above to use) #self.connect(self.am_conv, self.arg, (self.sub,0)) #self.connect(self.fm_conv, self.fm_arg, (self.sub,1)) self._build_gui(vbox, usrp_rate, demod_rate, audio_rate) if options.gain is None: g = self.subdev.gain_range() if True: # if no gain was specified, use the maximum gain available # (usefull for Basic_RX which is relatively deaf and the most probable board to be used for AM) # TODO: check db type to decide on default gain. #options.gain = float(g[1]) options.gain = float(g[0]+g[1])/2 #Changed with to 50% because the TV Tuner board doesn't need the full gain and that's what I'm using. else: # if no gain was specified, use the mid-point in dB options.gain = float(g[0]+g[1])/2 if options.volume is None: g = self.volume_range() options.volume = float(g[0]*3+g[1])/4 if abs(options.freq) < 1e3: options.freq *= 1e3 # set initial values self.set_gain(options.gain) self.set_vol(options.volume) if not(self.set_freq(options.freq)): self._set_status_msg("Failed to set initial frequency") def _set_status_msg(self, msg, which=0): self.frame.GetStatusBar().SetStatusText(msg, which) def _build_gui(self, vbox, usrp_rate, demod_rate, audio_rate): def _form_set_freq(kv): return self.set_freq(kv['freq']) if 0: self.src_fft = fftsink2.fft_sink_c(self.panel, title="Data from USRP", fft_size=512, sample_rate=usrp_rate, ref_scale=32768.0, ref_level=0.0, y_divs=12) self.connect (self.u, self.src_fft) vbox.Add (self.src_fft.win, 4, wx.EXPAND) if 0: self.post_filt_fft = fftsink2.fft_sink_c(self.panel, title="Post Channel filter", fft_size=512, sample_rate=demod_rate) self.connect (self.chan_filt, self.post_filt_fft) vbox.Add (self.post_filt_fft.win, 4, wx.EXPAND) if 0: pre_demod_fft = fftsink2.fft_sink_f(self.panel, title="Alt Signal pre-FM Demod", fft_size=512, sample_rate=demod_rate, y_per_div=10, ref_level=40) self.connect (self.am_demod, pre_demod_fft) vbox.Add (pre_demod_fft.win, 4, wx.EXPAND) if 0: post_demod_fft = fftsink2.fft_sink_f(self.panel, title="Alt Signal", fft_size=512, sample_rate=self.fm_audio_rate, y_per_div=10, ref_level=40) self.connect (self.fm_audio_filt, post_demod_fft) vbox.Add (post_demod_fft.win, 4, wx.EXPAND) if 1: oscope = scopesink2.scope_sink_c(self.panel, sample_rate=audio_rate, frame_decim=1, v_scale=2, t_scale=.05, num_inputs=2) #self.connect (self.mult, (oscope,1)) self.connect (self.fm_conv, (oscope,1)) self.connect (self.am_conv, (oscope,0)) vbox.Add (oscope.win, 4, wx.EXPAND) if 0: audio_fft = fftsink2.fft_sink_f(self.panel, title="Carrier AM signal", fft_size=512, sample_rate=audio_rate, y_per_div=10, ref_level=40) self.connect (self.audio_filt, audio_fft) vbox.Add (audio_fft.win, 4, wx.EXPAND) if 1: num_sink = numbersink2.number_sink_f (self.panel, unit='Radians',label="Angle", avg_alpha=1.0e-5,average=False, sample_rate=audio_rate, factor=1.0,base_value=0, minval=-7, maxval=7, ref_level=0, decimal_places=5,number_rate=15) self.connect(self.arg, num_sink) vbox.Add (num_sink.win, 1, wx.EXPAND) # control area form at bottom self.myform = myform = form.form() hbox = wx.BoxSizer(wx.HORIZONTAL) hbox.Add((5,0), 0) myform['freq'] = form.float_field( parent=self.panel, sizer=hbox, label="Freq", weight=1, callback=myform.check_input_and_call(_form_set_freq, self._set_status_msg)) hbox.Add((5,0), 0) myform['freq_slider'] = \ form.quantized_slider_field(parent=self.panel, sizer=hbox, weight=3, range=(108.0e6, 118.0e6, 1.0e5), callback=self.set_freq) hbox.Add((5,0), 0) vbox.Add(hbox, 0, wx.EXPAND) hbox = wx.BoxSizer(wx.HORIZONTAL) hbox.Add((5,0), 0) myform['volume'] = \ form.quantized_slider_field(parent=self.panel, sizer=hbox, label="Volume", weight=3, range=self.volume_range(), callback=self.set_vol) hbox.Add((5,0), 1) myform['gain'] = \ form.quantized_slider_field(parent=self.panel, sizer=hbox, label="Gain", weight=3, range=self.subdev.gain_range(), callback=self.set_gain) hbox.Add((5,0), 0) vbox.Add(hbox, 0, wx.EXPAND) try: self.knob = powermate.powermate(self.frame) self.rot = 0 powermate.EVT_POWERMATE_ROTATE (self.frame, self.on_rotate) powermate.EVT_POWERMATE_BUTTON (self.frame, self.on_button) except: print "FYI: No Powermate or Contour Knob found" def on_rotate (self, event): self.rot += event.delta if (self.state == "FREQ"): if self.rot >= 3: self.set_freq(self.freq + .1e6) self.rot -= 3 elif self.rot <=-3: self.set_freq(self.freq - .1e6) self.rot += 3 else: step = self.volume_range()[2] if self.rot >= 3: self.set_vol(self.vol + step) self.rot -= 3 elif self.rot <=-3: self.set_vol(self.vol - step) self.rot += 3 def on_button (self, event): if event.value == 0: # button up return self.rot = 0 if self.state == "FREQ": self.state = "VOL" else: self.state = "FREQ" self.update_status_bar () def set_vol (self, vol): g = self.volume_range() self.vol = max(g[0], min(g[1], vol)) self.volume_control.set_k(10**(self.vol/10)) self.myform['volume'].set_value(self.vol) self.update_status_bar () def set_freq(self, target_freq): """ Set the center frequency we're interested in. @param target_freq: frequency in Hz @rypte: bool Tuning is a two step process. First we ask the front-end to tune as close to the desired frequency as it can. Then we use the result of that operation and our target_frequency to determine the value for the digital down converter. """ r = usrp.tune(self.u, 0, self.subdev, target_freq + self.IF_freq) #TODO: check if db is inverting the spectrum or not to decide if we should do + self.IF_freq or - self.IF_freq if r: self.freq = target_freq self.myform['freq'].set_value(target_freq) # update displayed value self.myform['freq_slider'].set_value(target_freq) # update displayed value self.update_status_bar() self._set_status_msg("OK", 0) return True self._set_status_msg("Failed", 0) return False def set_gain(self, gain): self.myform['gain'].set_value(gain) # update displayed value self.subdev.set_gain(gain) def update_status_bar (self): msg = "Volume:%r Setting:%s" % (self.vol, self.state) self._set_status_msg(msg, 1) try: self.src_fft.set_baseband_freq(self.freq) except: None def volume_range(self): return (-40.0, 0.0, 0.5) if __name__ == '__main__': app = stdgui2.stdapp (vor_rx_block, "USRP VOR RX") app.MainLoop () _______________________________________________ Discuss-gnuradio mailing list Discuss-gnuradio@gnu.org http://lists.gnu.org/mailman/listinfo/discuss-gnuradio