As you may have noticed in some of my previous posts, I have dedicated some time to learn the GPS basics, which I would like to share with you today. Those with a matlab installed on their machine will be of course able to download all the source codes and run them (hopefully). And the best part: I have made a small GUI in matlab for signal processing :) I have even spent some time improving its appearance, robustness and speed, although some methods will still require your patience and in some cases a good processor as well :)

And because GPU programming is my hobby, I have reprogrammed the slowest method to run on the GPU ! :) Thus readers equipped with both 64-bit Matlab and new nVidia drivers for their graphics processor of GeForce 400+ will definitely feel like in heaven. I will start with some intro to the doppler effect,then I will add a few notes about GPS errors and then switch to my favorite part: Signal processing.

### Doppler Effect

Those familiar with doppler can skip this section. For others, believe me in saying that the satellites are moving across the sky. To be correct, GPS Satellites are located at the MEO (Medium Earth Orbit), which is located approximately 20 350km above the earth surface. Each satellite of the constellation orbits the Earth at a speed of about 3.9km/h. The Speed can be easily calculated, knowing that the satellite’s orbit time is 11h and 58 mins and taking into account the average Earth radius of 6378km. For simplicity, assume that the shape o the earth is a simple sphere. As the satellite is moving across the sky, it transmits a ranging signal on a frequency “f” (Which lies for GPS in L1-Band) and is about 1575,42 MHz. Because the satellite is moving with respect to the receiver, a so called Doppler Effect occurs. As a result the receiver receives a different frequency than the satellite transmitted. The same occurs when an ambulance driver passes nearby an observer,resulting in higher frequency as he drives towards and lower frequency as he moves away.

It is however important to note that the shift in frequency is dependent on the radial speed of the satellite, rather than on the tangential speed (3.9km/h). In other words it depends on the the speed between the receiver and the transmitter. Lets assume that the earth is not rotating and the receiver is static. I have created a small demo showing the doppler offset value in time when the satellite is visible on the sky:

As we can see, the distance to the satellite is highest when the satellite is reaching the horizon. The radial speed is highest at this point and is approximately 930m/s resulting in a doppler frequency offset of +-5KHz. If we take into account additional factor such as satellite orbit planes, rotation of the earth, position of the receiver( Imagine the difference between standing on the north pole and near the equator) and the phase noise of the local oscilator in the receiver, we can generally expect a higher Doppler offsets.

### Multipath and other Errors

There is a whole lot of problems, that can occur on the way from the satellite to the receiver. The receiver’s position is influenced by all kinds of errors, caused for example by imperfect GPS propagation models, or clock stability onboard the satellites (Although they are extremely precise). The position is generally dependent on the number of satellites visible and on their constellation, eg. their placement on the sky or constellation. To calculate a GPS position, at least 4 satellites are needed (x,y,z) coordinates and the last 4th to correct the time offset in the clock base of the receiver. One of the common problems in cities is a bad visibility of satellites (They are hidden behind buildings). This overall increases the error mainly due to a bad constellation settings (The best position can be obtained when a receiver is in the center and satellites are located perpendicular to each other).

With increasing number of the satellites,the receiver has more information and can thus calculate the position with better accuracy. A minor error also comes from the imperfect (Kepler-based) model of the trajectory of the satellites, which are extrapolated in time and updated each 2 Hrs from the Ground Control Center. Now lets consider a few propagation errors ;) One of the biggest problems was (Not mentioning disabled SA GPS mode – Selective Availability) a propagation error in the ionosphere. For simplicity this error is caused by the bending of the GPS signal due to change in relative permitivity in the Ionosphere. This error is however Frequency dependent and can be eliminated by receiving signals on at least 2 different frequencies (Which is common nowdays – L1 – 1500MHz and L2 – 1200 MHz).

Lets take a closer look also at the error caused by multipath propagation. This error is caused by a shift in the autocorrellation peak of the PRN sequence (Ranging code) due to the addition of the multipath signal to the direct signal at the receiver’s antenna. As a result, the correlation peak is somewhat shifted and the decoding algorithm (Lets assume the Early-Late Correlator) is unable to locate the peak exactly. By choosing a smaller distance between the Early and Late branches, we can estimate the peak’s position more accurately, however this choice results in a bad synchronization (The E and L are close). On the other hand, if we choose a higher distance of the E-L branches, the synchronization is more robust, although the eventual error rises:

As we can see, the eventual error can go as high as 80m. Due to the dilem between choosing the right EL distance, most of the market receivers are equipped with several correlators rather than with only 2 :) The Negative error is simulated by subtracting the multipath signal from the direct signal.

I believe I have not mentioned the most important thing about GPS, which is that the receiver’s position is calculated from measuring the distances to the visible satellites. The distance or virtual distance can be estimated from the correlation peak of the GPS ranging code. There are 37 available PRN sequences of length 1023 samples. The satellite broadcasts a signal modulated by this sequence so that each 1023 samples are repeated each 1ms, in other words 1.023Mbit/s . Each satellite is distinguishable by its own specific ranging code (Thus the GPS system is CDMA based unlike for example GLONASS, which is FDMA based).

In order to measure the distance, we need to find the correlation peak both in time and frequency (Due to the Doppler effect). The initial acquisition algorithm has therefore to scan the incoming signal for all PRN codes in 2D, which is time consuming, so most available market receivers are usually equipped with a build-in signal processor or some kind of accelerator.

Although there are other Acquisition methods, I will discuss only these:

- Maximum Likelihood (ML)
- FFT-Based Correlation
- Double Block Zero Padding (DBZP)

Among all of them, the least efficient is the ML algorithm, which integrates the signal in 2D space. The longer the integration time, the more accurate the results. However since GPS signal is modulated with a BPSK modulation of 50bits/sec (Broadcasting Navigation data such as the kepler parameters). The integration time is limited. This problem was solved by introducing a so called “pilot signals” without data modulation. The ML formula is quite easy:

In discrete form,we of course switch from integration to summation. Lets assume that the doppler is between -6000Hz and 6000,the sweeping interval is 100Hz and the signal time is 6600 samples. This results in integrating (12000/100 + 1)*6600 times to 6600. As you can guess, this is definitely not something you would ever want to calculate. My first implementation was so slow,that I could go cycling in the meanwhile. To be honest, it took about 35 minutes on my intel i5-750 processor (Although it wasnt optimized at all :D ). The ML is however extremely easy to implement on the nVidia CUDA architecture, so I reprogrammed the computation on the GPU, thus reducing the time from 35 minutes to just a few freaking ms! :D … But I had bad conscience about the CPU version, so that I have recently recoded the algorithm and removed 2 of 3 for-cycles from the code and implemented the computation via vector operations. Yet still, its its extremely slow in comparison with any other methods.

The FFT-Based Correlation method is very fast, in fact, its quite hard to set the settings so that you could walk away for a few minutes away from the computer. Instead of integrating, we calculate the autocorrelation function for each doppler offset specified. So that all integration steps for a given doppler offset are replaced only with a 2 (But can be only 1) FFT + 1 IFFT operations, which overall heavily decreases the computing demands. The correlation is done in matlab as simple as:

abs( ifft( fft(GPSSignal).*conj( fft(Replica) ) ) ).

I have also implemented another quite efficient method for the initial GPS acquisition (DBZP). This is actually the only method,that requires you to think a bit how to implement it :) But even then its not that hard, however the best explanation seems to be a graph or flowchart:

The input signal is divided into Nfd blocks of NTc size (Can be specified by user). The same applies for PRN ranging code. While the GPS Signal blocks are padded with another GPS Signal blocks to form a double-sized blocks, the PRN blocks are padded with zeros instead. A cyclic correlation is applied to each of these double-sized blocks and only the first half of the result is left and added to the first column of the final 2D correlation matrix. In the next step, we shift the PRN blocks by 1 (Blocks, not their data!), apply the same operations and fill the second column of the correlation matrix. And so on and on until we fill the whole 2D correlation matrix. In the end a 1D-FFT is applied to each of the columns of the 2D correlation matrix. Due to a relatively large amount of short FFTs, this method is not recommended for FPGA implementations. I found this algorithm quite interesting, although its a bit slower than the FFT-based correlation.

In the end,I would love to present my Self-Coded Matlab GUI (That was originally a part of my semestral work at the CTU), where you can play with all of these methods and signals. You can for example try to compute the crosscorrelation products (Use PRNx for generation and PRNy for decoding). Alternatively, you can play with integration time and SNR, which is modeled via gaussian noise. You can also try to sample in frequency domain with a low resolution and see what happens! :) Feel free to edit, share and post comments.

### Posittioning

Once we obtain the pseudoranges from at least 4 satellites, we are able to calculate the Receiver’s position (Also knowing the exact satellite positions). There are at least 3 Methods available for calculating the receiver’s position, which include: Kalman filtration, Direct search method and the **Least Squares Estimation**. We will talk more about the last option ,which unlike the direct search, can take into account more than just 4 satellites, which usually tends to increase the precision of the position. The Kalman filtration can be even more accurate, however its implementation is a bit harder to understand.

The LSE is an iterative method, which means, that we have to specify an initial position estimate, which in most scenarios is [0,0,0] – that is the middle of the Earth. Once we however know a better estimation,we can use that one to lower the number of iterations needed. Either way,once we got enough visible satellites,the method usually converges very fast anyway (Just a few iterations). During each Iteration,we create a so called **Direction Cosines Matrix A** and Calculate the **distance residues vector b** from the measured pseudoranges. Then we calculate the Position Correction, check against some predefined threshold and add the correction to our position. The exact equations are available here:

Where “**i**” denotes a satellite number, **di** is the satellite’s measured pseudorange and **Po** is the current receiver position (0,0,0 in the first iteration). If the correction is small enough, we can terminate the method and consider our position accurate, otherwise, we continue onto the next iteration. The next chart shows the corrections during each of the iterations:

We can see, that the position is almost accurate after just 4 iterations. This however depends on the constellation and number of visible satellites in general. Furthermore, we are ignoring all kind of errors, that will definitely influence the convergence. I have simulated the situation with receiver being located approximately at the Mount Mc. Kinley Mountain in Canada (Denali 6190 m.). And with a total of 8 satellites being visible (7 orbiting the equator and one directly above the north pole) – Please note that this is not a real-life scenario, as the situation completely ignores the real GPS orbits. I have also created a simple chart showing the receiver position along with the visible satellites and the Earth:

Feel free to download the matlab code and play with initial position and the number of satellites. Just one more note: My model of the Earth is quite simplified – its just a sphere, so a spherical coordinates can be used to set the Receiver’s position.

- Visual Studio 2013 MexFile – CUDA ML
- Doppler Script (Matlab)
- Delay Script (Matlab)
- Positioning (matlab)
- Whole GUI (Matlab)*

*May require you to update java, install MS Redistributables and or use newer version of Matlab (My own version is 64-Bit R2015a). The CUDA-ML is available only for 64-bit systems :( I am sorry for this.

Edit 7.1.2016: Minor bug fixes and Code update :)