[this post is a complete rewrite of the initial article , lost in my blog crash]
I. Introduction : What is a VOR, my motivation
The VHF Omni Range (VOR) system was designed in 1937 and deployed by 1946 to provide a solution for airplanes to estimate their position and to follow their route by indicating bearing angle towards ground radio transmitters. At that time the system had to be completely based on analogue hardware and the design engineers created a sophisticated modulation scheme for easy processing at the receiver. There are now around 3000 VOR transmitters over the world, and this number is now decreasing, progressively replaced by GNSS systems like GPS or the coming European Galileo. Operating in the 108 to 118 MHz band, these transmitters have been broadcasting their dual-modulation signals towards the aircrafts for more than 50 years and are still widely used by all types of aircrafts (see [1]). For several works on passive radar [see refs 2,3 and 4] I had to study in depth the signal transmitted by this system. One can find different details on the signal structure, my motivation was to check how accurate they were. In this post I will describe the technique I used to process the signals acquired using a DVB-T SDR dongle and MATLAB.
II. Signal model – How it works
A typical VOR transmitter is shown in figure 1 (image from Wikipedia)
Figure 1 – VOR ground transmitter
A central antenna transmits the REF signal, with a 30 Hz frequency-modulated signal. Meanwhile, other antennas (initially mechanically rotating and now using phase shifting techniques) transmit on a second carrier an amplitude modulated signal so that, from a given angle, the aircraft receives the two signals with a phase difference proportional to bearing to the station (figure 2).
Figure 2
Figure 3 shows the internal structure of a VOR receiver, as they were built in the “analogue era” before using digital signal processing. An envelope detector extracts the AM modulation at Intermediate Frequency (IF). The resulting signal is split in two branches:
Figure 3
- A first branch, containing the REF signal (showing the North direction) is band-pass filtered and the FM sub carrier extracted,
- A second branch, containing the VAR signal (indicating the current angle to the VOR) is low-pass filtered.
The two resulting signals phases are then compared and the resulting difference is amplified to drive the instrument needle. We can find in the literature the following equation (see figure 4) to describe the time domain signals transmitted by the VOR ground station:
Figure 4 – VOR equation
To check this model is accurate… lets find a signal, analyse it and compare with this equation !
III Field validation – collecting RF data and analysis
I leave a few miles away the Rambouillet VOR transmitter (callsign is RBT), active on 114.7 MHz. I decided to go next to it with a RTL-SDR dongle and collect some samples to see what’s inside the signal. Figure 5 shows my car location respective to the transmitter.
Figure 5 – VOR transmitter position
I used HDSDR and the small antenna shipped with the dongle. Figure 6 below shows the waterfall spectrum I got :
Figure 6 – HDSDR VOR Waterfall
The VOR signals were acquired using the 1M samples per second setting, I manually set the gain to have a good SNR but no saturation (I was very close to the VOR, usually the transmitters use around 200 Watts).
YOU CAN DOWNLOAD THE RAW HDSDR file (1MHz bandwidth) HERE : http://www.f4gkr.org/wp-content/uploads/public/VOR/HDSDR_20140205_185101Z_114684kHz_RF.wav
What can we say from the waterfall (from center to side bands) ?
- We have one narrow band signal in the middle, consistent with the VAR signal definition;
- We have one audio carrier, amplitude modulated (symmetric). Use the waterfall to decode the signal ;-) “BT” is the end of the VOR callsign very slow keying
- One other ‘multi trace’ symmetric carrier approx 9600 Hz away from central peak, should be the REF signal.
Let’s see how we can decode and process this signal with Matlab (or Octave).
First, we need to load the audio (RF) samples from the wav file into the working space :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
close all; clear all; %% load the wav file recorded with HDSDR and RTL-SDR dongle center = 114.684e6 ; % our center frequency [raw_samples, Fs] = wavread('HDSDR_20140205_185101Z_114684kHz_RF.wav'); % we want complex samples, this wavread generates Left and Right data baseband = raw_samples(:,1)+sqrt(-1)*raw_samples(:,2); clear raw_samples; % to limit memory usage... Te = 1/Fs; %% give a look to the spectrum L = length( baseband); freqs = (-Fs/2:Fs/L:Fs/2-Fs/L)+center; figure; plot( freqs, fftshift( 20*log10( 1/L*abs(fft(baseband))))); title 'complete spectrum'; xlabel 'freq'; ylabel 'magnitude'; grid on; |
this loads (be patient…) the complete file, processes the Left and Right to complex samples and plots the spectrum of the complete file. You should get something like this:
VOR spectrum
We collected much more than required, and the file is huge… but lets zoom around 114.7:
VOR spectrum – zoom
Okay, consistent with HDSDR display. Next step is to remove extra band and just keep the useful bandwidth. If the math showed before is correct, we should only keep from -9600+REF_WIDTH to +9600+REF_WIDTH. That’s around 30 KHz wide.
We will use a 15 KHz low pass filter, it will remove all but what we want to keep. Let’s go and see what we get (we start with a simple 128 taps filter)
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
%% design a filter to remove anything else but our -15 to +15 KHz Band = 15e3 ; taps = 128 ; lpf=fir1( taps,Band/Fs); % apply the filter to the data filtered_baseband = filter( lpf, 1, baseband); % plot result figure; plot( freqs, fftshift( 20*log10( 1/L*abs(fft(filtered_baseband))))); title 'filtered spectrum'; xlabel 'freq'; ylabel 'magnitude'; grid on; |
at first glance that looks nice… but zoom in the center, you should have that :
After low-pass filtering
It looks like we are a bit off the center of the spectrum, so we are also removing wanted signal… we have to “shift left” the spectrum, so we need to mix it (multiply) with a complex signal of frequency -offset. I estimated visually on spectrum the offset to be 9950 Hz. Our software mixer is like follows:
|
1 2 3 4 5 6 7 8 9 10 11 |
%% recenter our signal offset = -9950; t = (1:L)*Te; % generate our mixing signal LO = exp( sqrt(-1)*2*pi*offset* t.' ); % mix ! baseband_shifted = baseband .* LO ; clear baseband; % give fresh air to matlab |
apply our filtering like before, plot and look what it is now:
VOR spectrum after mixing – zoom
That’s much better, but it seems :
- The removal of extra signals is not that good,
- The shape of the noise is not that flat …
this is because our initial filter cutting at 15 KHz with 128 taps is not sharp enough. As we have time and we are rich … we change the settings of the filter to :
|
1 2 3 4 5 |
%% design a filter to remove anything else but our -16 to +16 KHz Band = 16e3 ; taps = 1024 ; lpf=fir1( taps,Band/Fs); |
(feel free to try other values of taps, you’ll see how long it takes ;-)). This leads to:
Zoom after extending the FIR low-pass filter number of taps
As our signal is now band-limited, we have in fact too much samples… we can “downsample” (lower sampling rate) and trash most of them. In fact we have an initial sampling rate of 1028571 samples per second, but we only need 30000 samples per second (our useful bandwidth). this means we can keep only 30000/1028571 samples. To make it simple, we’ll keep 1/32 of the band.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
%% now decimate to remove unwanted samples decim = filtered_baseband(1:32:end); % update internal vars accordingly Fsd = Fs/32 ; L = length(decim); freqs = (-Fsd/2:Fsd/L:Fsd/2-Fsd/L); % clear unused variables clear filtered_baseband; clear t; clear LO; |
OK ! we have our base band samples, and now Matlab breathes again…
Here is the final Matlab code for the VOR signal extraction :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 |
%%=========================================================================== % This Software is released under the "Simplified BSD License" % Copyright 2014 Sylvain AZARIAN - F4GKR. All rights reserved. % Redistribution and use in source and binary forms, with or without % modification, are permitted provided that the following conditions are met: % % 1. Redistributions of source code must retain the above copyright notice, % this list of conditions and the following disclaimer. % % 2. Redistributions in binary form must reproduce the above copyright notice, % this list of conditions and the following disclaimer in the % documentation and/or other materials % provided with the distribution. % % THIS SOFTWARE IS PROVIDED BY Sylvain AZARIAN ``AS IS'' AND ANY EXPRESS % OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES % OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. % IN NO EVENT SHALL Sylvain AZARIAN OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, % INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES % (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR % SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) % HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, % STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) % ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF % ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. % % The views and conclusions contained in the software and documentation % are those of the authors and should not be interpreted as representing % official policies, either expressed or implied, of Sylvain AZARIAN. %========================================================================================== % code to analyse VOR signals collected by F4GKR at Rambouillet VOR % VOR is transmitting on 114.7 MHz with ID RBT % recorded using RTL-SDR dongle, HDSDR %% close all; clear all; %% load the wav file recorded with HDSDR and RTL-SDR dongle center = 114.684e6 ; % our center frequency [raw_samples, Fs] = wavread('HDSDR_20140205_185101Z_114684kHz_RF.wav'); % we want complex samples, this wavread generates Left and Right data baseband = raw_samples(:,1)+sqrt(-1)*raw_samples(:,2); clear raw_samples; % matlab needs memory ... remove unused Te = 1/Fs; %% give a look to the spectrum L = length( baseband); freqs = (-Fs/2:Fs/L:Fs/2-Fs/L) + center ; figure; plot( freqs, fftshift( 20*log10( 1/L*abs(fft(baseband))))); title 'complete spectrum'; xlabel 'freq'; ylabel 'magnitude'; grid on; %% recenter our signal offset = -9950; t = (1:L)*Te; LO = exp( sqrt(-1)*2*pi*offset* t.' ); baseband_shifted = baseband .* LO ; clear baseband; %% design a filter to remove anything else but our -16 to +16 KHz Band = 16e3 ; taps = 1024 ; lpf=fir1( taps,Band/Fs); % apply the filter to the data filtered_baseband = filter( lpf, 1, baseband_shifted); clear baseband_shifted ; % give fresh air to matlab, remove unused data % plot result figure; plot( freqs, fftshift( 20*log10( 1/L*abs(fft(filtered_baseband))))); title 'filtered+centered spectrum'; xlabel 'freq'; ylabel 'magnitude'; grid on; %% now decimate to remove unwanted samples decim = filtered_baseband(1:32:end); % update internal vars accordingly Fsd = Fs/32 ; L = length(decim); freqs = (-Fsd/2:Fsd/L:Fsd/2-Fsd/L); % clear unused variables clear filtered_baseband; clear t; clear LO; % plot result figure; plot( freqs, fftshift( 20*log10( 1/L*abs(fft(decim))))); title 'decimated spectrum'; xlabel 'freq'; ylabel 'magnitude'; grid on; |
Stay tuned for tutorial continuation !
References
[1] Nordian – “VOR and Doppler VOR” – http://www.nordian.net/
[2] Automatic real-time collection of RCS of airplanes in a real bistatic low-frequency configuration using a software defined passive radar based on illuminators of opportunity, J. Pisane, S. Azarian
[3] Experimental RCS acquisition system: Using software defined radio to build a classification dataset, Sylvain Azarian, Jonathan Pisane, Radar (Radar), 2013 International Conference on
[4] Automatic Target Recognition for Passive Radar; J. Pisane, S. Azarian, Aerospace and Electronic Systems, IEEE Transactions on, 2014/1