Join GitHub today
GitHub is home to over 28 million developers working together to host and review code, manage projects, and build software together.Sign up
Using an extended Kalman filter for GPS position solution and clock oscillator calibration #147
this is not an issue but a possible enhancement.
As you have written here (#146 (comment)) it would be interesting to improve the accuracy of the GPS position+time solution.
I had a look at the GPS position solutions, placing the GPS patch antenna just outside my window on a short wooden boom at a position where more than half of the sky is obstructed with buildings. With this setup, GPS position solutions are available for much less than 50% of the time.
For about 10 hours I saved the GPS measurements, i.e.,
Note that is is very preliminary as
The EKF is used starting from the time a single point position solution (SPP) is first available.
I find that sometimes the GPS satellite transmit times
This looks like a bug in the GPS code as the time offsets are all exactly 20ms (with relative precision of < 5e-6); one possible location to start investigating this bug would be this one https://github.com/jks-prv/Beagle_SDR_GPS/blob/master/gps/solve.cpp#L157
For linear algebra I use the eigen library; I have checked that it compiles on the beagle-bone and does generate NEON SIMD code.
In terms of lines of code it would add ~200-300 lines and probably remove some of the lines for outlier filtering of the oscillator frequency offset.
Let me know what you think about this.
Hi Christoph. Top-notch work as always. I've known the term Kalman filtering for a long time but have never looked into the details.
This looks like a bug in the GPS code as the time offsets are all exactly 20ms (with relative precision of < 5e-6)
I would like to know more about this. Do you mean the outliers are all multiples of 20ms? And are these the individual uncorrected PRN times? i.e. t_tx I looked into the contributions of the return value of GetClock() but didn't see an obvious problem. Even for values where my clock_correction() routine determines there is an outlier. The return value would be a multiple of 20ms if both "chips" and "ca_phase" were zero. But this doesn't seem to happen.
The difference between the outlier values and the expected values are multiples of 20ms
When an outlier is present, the position solution fails, so you will probably not see it in clock_correction(). You should be able to see these outliers in t_tx in epochs with >=4 satellites in which the position solution fails.
Below are all the outliers I found in the 10hours of data:
Except in one case (marked with **), dz is a multiple of 0.02 seconds. As 0.02 sec = 1/50Hz I suspect that these outliers come from the contribution of
outlier idx= 0 z=0.097945745 zp=0.077945701 dz=0.020000045
Okay, I understand now -- thanks. The outliers I was referring to occur in the clock_correction() routine when the newly computed ADC clock value is outside the ~50 ppm window. This happens because the t_rx of the last GPS solution has an unexpected value. And if you print the lat/lon for the last solution the latitude is usually fine but the longitude is typically way off (tens of degrees).
So I wasn't looking at t_tx for failed position solutions. Let me see what I can find..
Here is some more information about the extended Kalman filter:
Below is the filter output (in red) and the single point position solver output (blue) for the same data used for the plots in the 1st comment, now with adjusted process noise covariance:
As you can see GPS reception is marginal with my setup but using the Kalman filter almost 100% efficiency can be reached.
As I have mentioned before I use the Eigen library for linear algebra. One possibility is to have it as a git submodule in Beagle_SDR_GPS/pkgs, see https://github.com/hcab14/Beagle_SDR_GPS (work in progress). Let me know what you think.
Very sorry for the delayed reply. Lots of things going on, not all Kiwi related.
Unfortunately, my understanding of the GPS solution algorithm (MLAT?) is very limited. Can you clear up some terminology for me? Does SPP refer to the classic solution obtained with four or more sats? I understand how application of an EKF could help smooth out the variations in SPP solutions. What I don’t understand is how you can continue to generate valid position solutions based on subsequent pseudo-ranges from less than 4 sats (or even just one) in the intervening times between SPPs.
I think it’s fine to use the Eigen lib. My only concern is about referring to that package as a submodule within the Kiwi project instead of including a “snapshot” of it in the pkgs directory. Since all Kiwis get their updates from github if the Eigen guys make a faulty commit I wouldn’t want that to brick all Kiwis worldwide on the next Kiwi update. But I don’t understand the details of how git submodules work exactly.
Also, the code you added in your latest commit to your fork didn't seem to be called from gps/solve.cpp:Solve() anywhere. Was it just added in preparation for eventual integration? I guess so since you said "work in progress".
One way of looking at it is the following: suppose that the receiver position is known, then only the receiver time is unknown and for this one satellite is enough.
The EKF update step combines the knowledge about the filter state with the new measurements into a new state in two steps:
I have to read the git submodule documentation in more detail but seem to recall that the submodule points to a specific commit. It still might be more safe make a snapshot in case the Eigen repository gets deleted or changes its name.
Yes it is not yet integrated into the KiwiSDR software - I did not have time for this yet
Okay, so this is like the "position and hold" mode of GPS timing receivers where the user position becomes a constant.
I see that now. That would work fine and makes it easier to switch to a new commit if they fix an essential bug.
I am still exploring the issue, but I think I may have found the cause of the occasional large outlier you noted above.
Andrew’s code takes care to “snapshot” the correct sat ephemeris data when the reception times for all sats are being computed. The largest contributor to the reception time is of course the TOW (time-of-week, in seconds) which is part of the ephemerides. These two processes, the decoding of the subframe data for each sat which stores the most recent TOW and the atomic computing of the reception time of all sats for each computed solution, are implemented as two separate asynchronous tasks (actually the processing of subframe data is a separate task per sat due to the bit-level processing involved). Therefore there has to be careful synchronization between them. If there is a problem updating the TOW for a sat, such as due to a parity error in the decode, the reception time code must be aware of that fact and prevent the sat from participating in the current positioning solution. Andrew’s code handles this case properly.
What seems not to be handled is the case where there is an error in the received 8-bit sync word (8 bits out of 300 per subframe) that causes the entire subframe not to be detected and hence no parity check even attempted. In this case the reception time code will use the TOW stored in the ephemerides from the last correctly decoded subframe which of course would lead to a huge error in the computed time.
Example: Suppose subframe 2 is decoded correctly and its TOW stored in the ephemerides. And then an asynchronous solution is attempted while subframe 3 is in the process of being received and decoded (NB: the solution task is run every 5 seconds and a subframe takes 6 seconds to transmit, i.e. 300 bits @50bps). This is the normal case and the code works fine. The reception time code uses the subframe 2 TOW and adds to it the exact amount of time that the sat has advanced into transmitting subframe 3 (by looking at how many bits have been received, the current position of the code generator in the FPGA etc.) If the decoding of subframe 3 turns up a parity error at any point then a flag is set (the “probation” flag in the code) which prevents the reception time code from using this sat in any solution computation for a few subframe times. But suppose one of the bits of the subframe 3 sync word is received incorrectly. The subframe 3 detection is missed entirely and its TOW will never be stored. And the reception time code is never told to not use the sat. So when a solution is attempted during the subframe 4 time the TOW from subframe 2 will be used instead of the correct subframe 3 which was never stored. The frequency of this situation is reduced somewhat because a sync pattern match can occur in the middle of normal subframe data. And that would trigger the parity checking code which would almost certainly fail and thus set the probation flag.
I discovered this while looking at the code in detail because Galileo, unlike Navstar, doesn’t transmit a TOW per subframe (page in Galileo terminology). And so the concept of current TOW will have to be handled differently for Galileo. The reception time code will possibly have to synthesize the current TOW based on the last one received and a local measurement of how long it’s been.