Search code examples
gpskalman-filter

How to use Kalman filter with degrees 0-360


I have to filter heading provided by GPS. Heading is provided as degrees 0-360. When there is transition between 0-360 or 360-0 Kalman goes crazy and my filtered heading is completely wrong.

Is There any possibility to use Kalman with that kind of data?


Solution

  • I've dealt with similar problems by judicious use of the (C math library) function remainder. remainder(a,b) gives the remainder on dividing a by b, and is between -b/2 and b/2. NOT between 0 and b.

    For example in the simple case where you have just the heading as a state, the general kalman measurement update which I'll write as

    X^ = X~ + K*(Y-H*X~)
    

    (where X~ is the updated state, X^ the predicted state at time of measurement, Y the measurement, H the measurement model matrix, K the Kalman gain).

    becomes

    dY = remainder( Y-X~, 360)
    X^ = remainder( X~ + K*dY, 360)
    

    There are two remainders to avoid the filtered heading increasing for ever, but instead remaining in [-180,180]. Though it's not strictly necessary to do that I find it odd when my filter says heading is say 360 million degrees plus 15, instead of 15.

    If you are modelling the rate of change of heading as well, so that your state is say h,hdot (the heading and rate of change of heading) its worth doing a similar thing in the time update too. That is you should have when updating the time by dt say

    h~ = remainder( h + dt*hdot, 360)
    hdot~ = hdot
    

    Again, this is to keep h in [-180, 180].

    The updates for the state covariance, and how K is computed, are unchanged.

    If you are using some package that implements the kalman filter equations, you can implement the same thing by altering the observed value before passing it to the measurement update. If Y is the observed value, and X the state, change Y to

    Y^ = X + remainder( Y-X, 360)
    

    We will then have

    Y^-X = remainder( Y-X, 360) 
    

    For the reasons above, I'd also recommend 're-normalising' X^ by

    X^ = remainder( X^, 360)
    

    after the measurement update