Summary of project to integrate GPS1A
board from gpscreations
with the software receiver tools in the GPSTk. Initial Project Description
A sample GPS1A data file (gnssGood.bin), a nav file (rin08_207.08n), and sample commands and output (gnssGoodInfo.txt) are available for testing and quickly getting started with the swrx software. Contact Jon Little for their location.
The board uses the SiGe SE4110L
GPS RF front end IC to collect data from the L1 band. The GPS1A configuration uses a 16.3676 MHz sampling rate, down-converts the signal to a 4.1304 MHz intermediate frequency (1575.42 - 96*clock), and outputs 2-bit real samples on sign and magnitude pins. A cypress IC interfaces with a USB connection. The low-level software (ogusb-lite.exe) provided with the board takes the samples and dumps them to a file as 8 bit values...inefficient but easy to work with. The USB interface can be downloaded from Sourceforge
(gn3s-linux). OSGPS also has a software receiver available for download there, and a simulation data file. OSGPS Homepage
Because the hardware for this project only deals with the L1 band, there is no reference to PRN code types other than C/A. However, the algorithms described can be modified for use with P code and the new codes.
The gpsSim simulator is one of the most useful tools in the GPSTk, because it allows us tight control over the attributes of a simulated signal. It can output on both the L1 and L2 bands with user defined codes, Doppler, noise and dummy nav bits. The sim can output at a lower IF by generating local oscillator samples and multiplying their conjugate by the simulated GPS signal. The resulting IF is the L1 freq minus the local oscillator frequency. The first part of this project included parameterizing the IF and sample rate of the simulator to be able to duplicate the signal generated by the GPS1A board.
Things to be careful about when using the simulator:
- Be sure to match output quantization to input type of next program in pipe. (Default quantization is usually float).
- The sim outputs 2 bands by outputting a sample from one band then the other. Since the GPS1A board is L1 band only, the acquisition code needs to know to throw away the second band when taking input from simulator. When working with the simulator, set the acquisition and tracking programs to -b 2 (default), when working with hilbert.cpp or other data from the GPS1A set them to -b 1.
- If using several periods of code, make sure to set the -t parameter (in ms) to generate enough samples.
- The nav signals on the sim are currently dummy signals that switch from 0 to 1 or 1 to 0 every (??) periods of C/A code.
- Code offsets for programs in the swrx folder are generally in microseconds, not chips.
We decided to work with both inphase and quadrature samples early in the project. A Hilbert transform pulls the quadrature component from the real GPS1A data. Using I and Q samples helps with code reuse and algorithm commonality in the toolkit. Also, since our acquisition method is FFT based, this can improve processing time by halving the sampling rate (N points of real data generate N/2 points of useful I and Q data, and in our FFT based acquisition we will already be dealing with complex values after the first FFT). The transform works by discarding half of the frequency components of the real input:
Our hilbert program reads real data from a SiGe output file (-i parameter) and outputs I and Q samples using the IQ2Stream class, a complex two bit stream that can be used as input to the other swrx tools.
is a free C library for performing FFT's. It is distributed under the terms of the GNU General Public License. The Hilbert transform and GPS Acquisition programs both link to FFTW. Because of this, I have not included them in the jamfile. The commands to compile and link them using g++ can be found at the top of each .cpp file.
Before passing data to the tracking program, we need to know the Doppler and code shift of the signal. The traditional approach is to correlate the input signal with a locally generated code/carrier signal for Doppler bins spanning at least +/-10KHz. The resulting [# shifts * # of samples to be multiplied * # of frequency bins] operations are time consuming in software.
We therefore used an FFT based acquisition (Parallel Code Phase Search) described by Tsui, to eliminate the code shifting process. Convolution in the time domain is the result of multiplication in the frequency domain, and correlation is similar to convolution, with a sign change. Therefore, the correlation can be performed with only [(#multiplies +FFT cost) * # freq bins] operations.
- Generate local code/carrier replicas for each frequency bin and convert them to frequency domain using an FFT.
- Convert input samples to frequency domain using an FFT and take their complex conjugate.
- Multiply the results, then convert back to time domain using an inverse FFT.
- The highest value's location will be the start of the C/A code and the correct Doppler frequency bin.
- Visual Example
Note: another acquisition method (Parallel Frequency Space Search) eliminates the need to search frequency bins. The input signal is be multiplied by C/A code replicas at all possible shift values. If the C/A code shift matches the input code shift, the spreading code will be removed and the result will be the carrier wave, which has a strong frequency component at the desired frequency. We have not implemented this method, because the Parallel Code Phase Search is usually faster for GPS. Visual of this type of acquisition
Example usage with live sky data:
hilbert -i data.bin | acquire -c 30 -p 3 -q 2 -b 1 -x 4.13 -r 8.184 -w 20000 -f 100
= acquire PRN 30, use three periods of C/A code, 2-bit quantization, 4.13 MHz IF, 16.368 MHz input sampling rate (divided in two by hilbert), +/-10 KHz Doppler search and 100 Hertz spacing of frequency bins.
The user can enter -c 0 parameter to acquire all 32 PRNs at once.
PRN: 30 - Doppler: -8600 Offset: 416.789 Height: 46.4584
- Tracker Input: -c c:1:30:416.789:-8600
Acquire can also output data from the correlation curve of the desired frequency bin for graphing purposes (more periods can be added to get a sharper peak / better SNR):
Sample correlation curve
As you can see, the results for one period are good enough to verify if the satellite is visible or not, but I would recommend using at least 3 periods to ensure correct frequency is passed to the tracker. In the example above, the actual Doppler shift is about -2280, and using 5 periods (same bin width) gives us -2300 instead of the -2000 obtained from one period.
Example of acquisition of a signal generated with the simulator:
gpsSim -q 2 -c c:1:1:50:5250:c -x 4.13 -r 16.368 -t 4 | acquire -c 1 -p 5 -q 2 -b 2 -x 4.13 -r 16.368 -w 12000 -f 50
PRN: 1 - Doppler: 5250 Offset: 50.0367 Height: 54.9682
- Tracker Input: -c c:1:1:50.0367:5250
The simulator is set for a 50 microsecond code shift and 5250 Hz Doppler. The simulator outputs 2 bands (acquire just ignores the 2nd) and we are now looking at 5 periods of PRN 1 with new frequency search width.
Our side-lobes from correlation of GPS1A data do not compare to theoretical C/A code correlation as well as we would have hoped, but there are similarities. At some point, I would like to get data using a better antenna to reduce noise and multipath. However, we were able to verify that our algorithm gives the same results as a time based correlation by acquiring a signal without noise from the gpsSim.
Comparison of Acquisition to Theoretical C/A Code Autocorrelation
Matlab File for Theoretical Plot
Guidelines for Choosing Frequency Bin Parameters
When choosing frequency bin size, we must consider the data length of the input signal (how many periods we choose to correlate), or else we may not see correlation. A good rule of thumb is to make the bin width:
If the local code/carrier signal is different than the input signal (frequency separation) by more than one cycle over a given data length, there will be no correlation between signals. There is partial correlation if the signals are within a cycle of each other. A 1 KHz change in signal will change the frequency separation by 1 cycle in 1 ms. So, if we use a 1KHz bin and one period of code, the greatest frequency separation will be 0.5 cycles (if the input is exactly in the middle of a bin, it will differ from the replicas by 500 Hz). We arbitrarily decide to make 0.5 cycles a cutoff point, so to keep the frequency separation at less than that value, our frequency bins should always be equal to or less than 1KHz/periods.
Bin Width Selection Examples
Frequency search width:
The maximum Doppler shift caused by satellite velocity is about +/- 5 KHz. Also, the error in local oscillator or the velocity of the receiver can cause the same effect to the signal. A good Doppler search width is +/- 10 KHz ('-w 20000' input parameter). If one already has a general idea of the Doppler shift (maybe a 1 period acquisition has already been done), narrowing the search width will decrease computation time.
The tools in the swrx folder can potentially be combined into a real time GPS receiver, although that is not their stated purpose. The FFT based acquisition method described above is a block method, in that the samples are not correlated directly on being generated, but gathered into a group and then correlated. We achieve O(N log N) performance, where N is the number of samples. A time based method is actually linear, but the curve is so steep that it is much slower for all but very small N. With the input at 16.368 MHz, a period of C/A code can be put through the hilbert transform and its shift and doppler acquired in less than a second (which is fast enough for real-time operation). The Doppler is not always precise enough to pass on to the tracking program, so sometimes more C/A periods are needed, but the speed can be improved because the frequency search width can be greatly reduced if centered around the Doppler from 1 period. The performance bottleneck is tracking, which currently takes about 5 minutes for a 1 minute data sample.
The tracking program uses an early minus late code tracker and a Costas carrier tracker to follow the code shift and frequency change and to de-modulate nav data. The tracker needs the Doppler within tens of Hertz (works fine as wide as 50) and the code phase within a microsecond or two. In verbose mode, the tracker will continue to output code and carrier info as well as the current nav bit.
The EML code tracker (Delay Locked Loop) correlates input samples with samples from three shifted replicas of the code/carrier (timed to be before, at, and after the expected code shift). The DLL error is calculated by:
and the code phase and carrier frequency offsets are then shifted by user defined multiples of this error (dllAlpha & dllBeta parameters).
The Costas tracker (Phase Locked Loop) uses the following formula:
and the carrier phase and carrier frequency offsets are then shifted by user defined multiples of this error (pllAlpha & pllBeta parameters).
hilbert -i data.bin | tracker -c c:1:4:170:-2300 -q 2 -b 1 -x 4.13 -r 8.184 -v
= code parameters (code:band:prn:offset:doppler), quantization, bands, IF, sample rate, verbose
This will output current Code Phase and Carrier Frequency offset, current nav bit, and other tracking loop information.
The tracker also demodulates nav data from the signal and checks for subframe parity using the EngNav class (called from NavFramer.cpp in the swrx file). Sub-frames occur every 6 seconds, and we have successfully demodulated all 9 complete sub-frames from visible satellites in a 60 second data file created with the GPS1A board. The tracker currently outputs subframe information every time a subframe is found, even if verbose mode is not on.
The tracker is now parallel at a very high level, using pthreads. (trackerMT.cpp, the old tracker is left as well to show the original algorithm). The user can enter multiple -c commands with the code parameters for each PRN to track. The input data is read into a buffer (right now holding one period of C/A code at a time) and a new thread is started to track that block for each PRN. Once all the threads are finished, the buffer is filled again with new input data.
Presentation of tasks completed over the summer
* End-of-Summer Presentation
Now included in position.cpp and RX.cpp
Based on: Tsui, James, Fundamentals of Global Positioning System Receivers, 2000 Section 9.9.
- Find the data points for the starts of a certain subframe for visible satellites.
- Guess a transmission time and generate estimated pseudoranges (all of which will have a "clock error" because of estimate)
- Correct this clock error using ephemeris data from nav bits (method is well documented)
- Calculate pseudoranges based on corrected transmission time.
We have added the position.cpp file, which contains the algorithm for finding position from tracker output (currently the user must enter data point values directly into the code). We are currently finding position accurate to within 100 meters or so and plan to continue improving/testing the algorithm. (position.cpp is included in the Jamfile).
Position.cpp Example usage:
position -e s011273a.08n -z 149478 -w 1499
Where -e is a RINEX ephemeris file, -z is the z-count obtained from the tracking program and -w is the week number.
Right now, the user must enter the data points from the tracking program directly into the code, like this:
// subframe 4 data points from gnss.bin: (z=360204)
Where the index is the PRN-1 and the dataPoint is the dp from the subframe output by the tracking program.
In RX.cpp, the tracker automatically populates the dataPoints arrays with subframe information from the tracker and finds a position solution every subframe (every 6 seconds).
RX.cpp Example usage:
hilbert -i /home/mdavis/docs/CarData/CarData.bin | ./RX -b 1 -q 2 -x 4.13 -r 8.184 -c c:1:30:416.789:-8800 -c c:1:2:410.191:-8300 -c c:1:10:836.144:-9000 -c c:1:15:174.609:-4000 -c c:1:18:255.254:-3500 -c c:1:24:183.284:-5700 -c c:1:26:907.747:-4300 -c c:1:29:355.327:-5400 -p 1 -e rin269.08n -w 1498
The -p parameter enables position solution, and if it is used, a rinex nav file and gps week number must be provided. The -c parameters were obtained using acquire.cpp:
hilbert -i data.bin | acquire -c 0 -p 5 -q 2 -b 1 -x 4.13 -r 8.184 -w 20000 -f 100
Position (ECEF): -756754.223563 -5465457.252978 3189012.200619
Time: 09/25/2008 04:03:12
# good SV's: 8
RMSResidual: 55.798621 meters
Position (ECEF): -756775.361255 -5465566.545293 3189208.578611
Time: 09/25/2008 04:03:18
# good SV's: 8
RMSResidual: 51.598548 meters
List of Information Needed to Get Position Solution:
1. .bin file from GPS1A board (or other software receiver hardware, provided the correct parameters are used)
2. GPS Week when data was taken.
3. RINEX nav file containing the time the data was taken (eventually we will generate our own nav data from subframes, but not yet implemented).
A sample GPS1A data file (gnssGood.bin), a nav file (rin08_207.08n), and sample commands and output (gnssGoodInfo.txt) are available for testing and getting started with the swrx software. Contact Jon Little for their location.
There seems to be a bug in NavFramer?
.cpp that sometimes gives us an inverted HOW and subframe #. NavFramer?
.cpp can be edited to output all the bits of a subframe, and I've seen that the preamble is the correct polarity, the parity checks, but somehow the HOW and SF# are still inverted. This has happened on 2 .bin files during tracking. In RX.cpp this will cause "text 0:Invalid week/seconds-of-week:" to be thrown.
Sometimes 1 or 3 periods is not enough to obtain an exact doppler shift (5 is usually sufficient but takes a long time). If satellites are not being tracked, use the -v parameter to check the tracker:
hilbert -i /home/mdavis/docs/CarData/CarData.bin | trackerMT -b 1 -q 2 -x 4.13 -r 8.184 -c c:1:30:416.789:-8600 -v
#h time dllErr codePO codeFO pllErr carrPO carrFO nav cp iad ely pmt lat pmtI pmtQ prn
#u ms % us Hz cyc cyc Hz - -- cnt % % % cnt cnt
# Taking input from (1 samples/epoch)
# Rx gain level: 1
1.0 -1.59 416.743 -0.08 -0.776 -8.8 -8794.41 0 21 8183 8.17 9.75 6.58 -8689 -7355 30
1.6 1.74 416.793 0.07 -0.774 -13.9 -8784.75 0 21 4775 5.94 10.32 7.68 -5337 -4574 30
2.6 -1.92 416.737 -0.03 -0.722 -22.6 -8777.78 0 21 8183 6.49 8.49 4.57 -6382 -7636 30
3. -0.49 416.723 -0.05 -0.855 -31.4 -8774.16 0 21 8183 7.24 9.70 6.75 -10199 -4989 30
The cp column should read 21 to indicate a lock and the nav bit should be changing at most every 10 lines. If not, acquire over more periods, or adjust the doppler up or down by a hundred Hz or so to find the correct doppler shift.
RX.cpp does not always exit gracefully when coming to the end of a file, throws a gpstk error code.