156

I'm working with GPS data, getting values every second and displaying current position on a map. The problem is that sometimes (specially when accuracy is low) the values vary a lot, making the current position to "jump" between distant points in the map.

I was wondering about some easy enough method to avoid this. As a first idea, I thought about discarding values with accuracy beyond certain threshold, but I guess there are some other better ways to do. What's the usual way programs perform this?

P̲̳x͓L̳
  • 3,615
  • 3
  • 29
  • 37
Al.
  • 1,811
  • 4
  • 15
  • 8
  • I feel the bad effects of the "GPS noise" when trying to calculate associated (derivative) values like speed and slope, which are very discontinuous specially for high sample rate tracklogs (since time has integer [one second] resolution). – heltonbiker Sep 06 '12 at 13:45
  • 4
    (also, if you are navigating through main roads, you can use "snap to roads" algorithm provided you have a good [correct, precise] roadmap dataset. Just a thought) – heltonbiker Sep 06 '12 at 13:51
  • I am facing this issue for best accuracy also. – ViruMax Dec 08 '14 at 15:42

12 Answers12

99

Here's a simple Kalman filter that could be used for exactly this situation. It came from some work I did on Android devices.

General Kalman filter theory is all about estimates for vectors, with the accuracy of the estimates represented by covariance matrices. However, for estimating location on Android devices the general theory reduces to a very simple case. Android location providers give the location as a latitude and longitude, together with an accuracy which is specified as a single number measured in metres. This means that instead of a covariance matrix, the accuracy in the Kalman filter can be measured by a single number, even though the location in the Kalman filter is a measured by two numbers. Also the fact that the latitude, longitude and metres are effectively all different units can be ignored, because if you put scaling factors into the Kalman filter to convert them all into the same units, then those scaling factors end up cancelling out when converting the results back into the original units.

The code could be improved, because it assumes that the best estimate of current location is the last known location, and if someone is moving it should be possible to use Android's sensors to produce a better estimate. The code has a single free parameter Q, expressed in metres per second, which describes how quickly the accuracy decays in the absence of any new location estimates. A higher Q parameter means that the accuracy decays faster. Kalman filters generally work better when the accuracy decays a bit quicker than one might expect, so for walking around with an Android phone I find that Q=3 metres per second works fine, even though I generally walk slower than that. But if travelling in a fast car a much larger number should obviously be used.

public class KalmanLatLong {
    private final float MinAccuracy = 1;

    private float Q_metres_per_second;    
    private long TimeStamp_milliseconds;
    private double lat;
    private double lng;
    private float variance; // P matrix.  Negative means object uninitialised.  NB: units irrelevant, as long as same units used throughout

    public KalmanLatLong(float Q_metres_per_second) { this.Q_metres_per_second = Q_metres_per_second; variance = -1; }

    public long get_TimeStamp() { return TimeStamp_milliseconds; }
    public double get_lat() { return lat; }
    public double get_lng() { return lng; }
    public float get_accuracy() { return (float)Math.sqrt(variance); }

    public void SetState(double lat, double lng, float accuracy, long TimeStamp_milliseconds) {
        this.lat=lat; this.lng=lng; variance = accuracy * accuracy; this.TimeStamp_milliseconds=TimeStamp_milliseconds;
    }

    /// <summary>
    /// Kalman filter processing for lattitude and longitude
    /// </summary>
    /// <param name="lat_measurement_degrees">new measurement of lattidude</param>
    /// <param name="lng_measurement">new measurement of longitude</param>
    /// <param name="accuracy">measurement of 1 standard deviation error in metres</param>
    /// <param name="TimeStamp_milliseconds">time of measurement</param>
    /// <returns>new state</returns>
    public void Process(double lat_measurement, double lng_measurement, float accuracy, long TimeStamp_milliseconds) {
        if (accuracy < MinAccuracy) accuracy = MinAccuracy;
        if (variance < 0) {
            // if variance < 0, object is unitialised, so initialise with current values
            this.TimeStamp_milliseconds = TimeStamp_milliseconds;
            lat=lat_measurement; lng = lng_measurement; variance = accuracy*accuracy; 
        } else {
            // else apply Kalman filter methodology

            long TimeInc_milliseconds = TimeStamp_milliseconds - this.TimeStamp_milliseconds;
            if (TimeInc_milliseconds > 0) {
                // time has moved on, so the uncertainty in the current position increases
                variance += TimeInc_milliseconds * Q_metres_per_second * Q_metres_per_second / 1000;
                this.TimeStamp_milliseconds = TimeStamp_milliseconds;
                // TO DO: USE VELOCITY INFORMATION HERE TO GET A BETTER ESTIMATE OF CURRENT POSITION
            }

            // Kalman gain matrix K = Covarariance * Inverse(Covariance + MeasurementVariance)
            // NB: because K is dimensionless, it doesn't matter that variance has different units to lat and lng
            float K = variance / (variance + accuracy * accuracy);
            // apply K
            lat += K * (lat_measurement - lat);
            lng += K * (lng_measurement - lng);
            // new Covarariance  matrix is (IdentityMatrix - K) * Covarariance 
            variance = (1 - K) * variance;
        }
    }
}
user1003545
  • 2,078
  • 1
  • 15
  • 15
Stochastically
  • 7,616
  • 5
  • 30
  • 58
  • 1
    Shouldn't the variance calculation be: variance += TimeInc_milliseconds * TimeInc_milliseconds * Q_metres_per_second * Q_metres_per_second / 1000000 – Horacio Sep 03 '13 at 08:22
  • 5
    @Horacio, I know why you think that, but no! Mathematically, the uncertainty here is being modelled by a Wiener process (see http://en.wikipedia.org/wiki/Wiener_process ) and with a Wiener process the variance grows linearly with time. The variable `Q_metres_per_second` corresponds to the variable `sigma` in the section "Related processes" in that Wikipedia article. `Q_metres_per_second` is a standard deviation and it's measured in metres, so metres and not metres/seconds are its units. It corresponds to the standard deviation of the distribution after 1 second has elapsed. – Stochastically Sep 10 '13 at 06:56
  • Thanks for the explanation. Can you elaborate on how to account velocity information to get a better estimate? can I simply set Q_meters_per_second equal to the velocity of the user (in meters per 1 sec) multiplied by some factor? say 1.5? – Horacio Oct 03 '13 at 05:22
  • @Horacio what you suggest won't work at all in terms of incorporating velocity to get a better estimate. I was trying to answer that question myself, and ended up asking a couple of questions elsewhere, see http://dsp.stackexchange.com/questions/8860/more-on-kalman-filter-for-position-and-velocity . In fact Kalman filters are quite sophisticated, and at the moment I don't think there's a simple answer to your question. The code that I posted here is a bit of an anomaly in the sense that for this special case, everything works out very nicely and it looks quite simple! – Stochastically Oct 03 '13 at 07:14
  • 3
    I tried this approach and the code, but it ended up shortening the total distance too much. Made it too imprecise. – Andreas Rudolph Sep 12 '14 at 08:14
  • How to use your class? Should I process (with your code) already received coordinates in onLocationChanged() method? – user2999943 Mar 24 '15 at 15:01
  • 1
    @user2999943 yes, use the code to process the coordinates that you get from onLocationChanged(). – Stochastically Apr 01 '15 at 07:38
  • Why do you think that this will improve the positions? The GPS chip internally heavilly uses the kalman filter, the output is already kalman filtered. – AlexWien Jul 06 '15 at 16:36
  • @AlexWien the above code is a simple Kalman filter :-). Whether it *improves* anything or not is subjective, and will depend on the objectives. Also, whether GPS readings are already Kalman filtered (and to what extent) is presumably hardware dependent. Finally, I haven't looked into the detailed mathematical theory but I wouldn't be surprised if someone told me that a Kalman filter applied to a Kalman filter is just a Kalman filter which is less sensitive to the inputs. If so, even if GPS readings are already filtered, if one wants less noisy readings another filter may be useful. – Stochastically Jul 06 '15 at 22:17
  • @Stochastically, thanks for this. It works really well for my use case. I'm still unclear on why you say adding the noisy velocity measurements (from the GPS points) to the model won't work. Maybe you mean calculating it is "baked in"? (I've been following http://dsp.stackexchange.com/questions/8860/more-on-kalman-filter-for-position-and-velocity but don't have the rep. to comment there). I'm on iOS, but I'm sure it's quite similar to Android. They provide a speed and course, seemingly without notion of accuracy. Could those measurements be fused with the position measurements? Thanks again! – pickaxe May 01 '17 at 22:05
  • @Stochastically, one more question: using the Weiner process here means that the variance will just increase to infinitum for streaming data, K --> 1, and eventually the filtered location will just be the input. does this poses a problem for your use case? Or do you reset variance at some point? – pickaxe May 03 '17 at 20:57
  • @pickaxe I'm not sure you understand what a Kalman filter is meant to do. If you you have perfect accuracy of your measurements then K=1 and the position is set to exactly the measurement. For completely imperfect measurements the accuracy value is huge so K=0 and the position is unaffected. – Stochastically May 04 '17 at 11:54
  • @Stochastically, yes, I understand. I don't think I said anything to the contrary regarding the behavior of K. But... I overlooked the re-calc of the covariance when trying to understand how it behaves over time. What I learned so far is that this is a really nice approach for when GPS data is sparse and using a physical model of motion produces terrible results. Thanks again. – pickaxe May 04 '17 at 16:21
  • 1
    @Stochastically - Thank you for this code. This is excellent. Do you have usage code (an example of how to use this?) I'm just ramping up on Kalman filters. – Naikrovek May 19 '17 at 15:51
  • @Naikrovek I used this an an android App which needed to know the location of a phone, but I haven't worked with this since 2013 so I don't have any more code to add. Good luck! – Stochastically May 19 '17 at 15:55
  • @Naikrovek: something like this: `filter.Process(geo.getLatitude(), geo.getLongitude(), (float) point.getLocation().getAccuracy(), point.getTimestamp().getMillis());` ` ` `AddPoint(filter.get_lat(), filter.get_lng(), point.getTimestamp();` – Glebka May 28 '18 at 11:13
  • @GLeBaTi what if we don't have the accuracy info? – Koray Dec 13 '18 at 12:55
  • 2
    @Koray if you don't have accuracy info then you can't use a Kalman filter. It's completely fundamental to what the Kalman filter is trying to do. – Stochastically Dec 14 '18 at 13:52
  • @Stochastically I've used your code with a static accuracy value and modified a little (ie. calculate acceleration and eliminate the points that creates more than a threshold) and the results look good. I was thinking to use a Douglas Peucker or a spline smoothing, but my results with your code were just fine. Do you think its ok to use your code like that? – Koray Dec 17 '18 at 08:51
  • 1
    @Koray a Kalman Filter has a current position estimate with some variance, and then a measurement comes along for a different position, also with some variance (N.B. variance = accuracy parameter). The formulas in the filter are derived so that the old position is combined with the new position estimate so as to produce an estimate with the lowest possible variance. So if there's no real accuracy parameter, the variances are wrong and the results will end up being a bit arbitrary. It will work if your variances are constant(ish) and you use that constant, but unreliable otherwise. – Stochastically Dec 17 '18 at 10:48
  • The way I understand it, the code is dampening the changes in lat/lon by taking only a part (controlled by K) of the difference to the previous lat/lon pair. This is a form of exponential averaging. Does it matter that this happening in the lat/lon domain and not in the euclidian domain? – Christian Lindig Apr 13 '19 at 12:44
  • 1
    @ChristianLindig I suggest you look at the mathematical basis for the Kalman filter. The idea is that the existing best guess location has an error, and the new measurement has an error too. The choice of K is calculated so as to combine the existing best guess location with the new measurement so as to minimize the error of the new best guess location :-) . – Stochastically Apr 14 '19 at 16:47
  • 1. Does anyone has examples of what parameters Q, accuracy for a car gps track? – Baglay Vyacheslav Jul 31 '19 at 11:41
  • How does your implementation handle crossing the international date line? Seems like the discontinuity in longitude might not be properly handled. eg 179.99999W -> 179.99999E – Matt Oct 16 '20 at 13:58
  • Without doubt, @Matt, this implementation does not handle the international date line. – Stochastically Oct 18 '20 at 01:21
  • 1
    Thanks for the answer @Stochastically. Would it make sense to update q_meters_per_second as you go? For example: keep track of total distance and time and divide to get a new q as you go (or you could take the last x points and use those). Before I implement, I'd like your thoughts. The reason I ask is the vastly different speeds people might, for example, run. – Daniel K Oct 07 '22 at 17:32
  • 2
    Wow, @DanielKanaan, a comment on an answer of mine that's almost 10 years old! You could certainly try a dynamic q_meters_per_second, but if you're not careful it might make the results unstable. Also, what may work well in one situation may not be so good in another, so testing in all situations that you're interested in is important. And as I said 10 years ago, a larger q_meters_per_second than one expects generally works well . Good luck! – Stochastically Oct 09 '22 at 02:34
  • The answer says it doesn't matter about working in degrees lat and lon because `K` is dimensionless. But in my location one degree latitude is about 1.5 times the distance of one degree longitude. So would it be better to convert into local OSGB easting and northing in metres, before applying the filter? – Weather Vane Feb 14 '23 at 20:20
  • @WeatherVane Not necessary. As I say in my answer "the fact that the latitude, longitude and metres are effectively all different units can be ignored, because if you put scaling factors into the Kalman filter to convert them all into the same units, then those scaling factors end up cancelling out when converting the results back into the original units." – Stochastically Feb 16 '23 at 02:13
81

What you are looking for is called a Kalman Filter. It's frequently used to smooth navigational data. It is not necessarily trivial, and there is a lot of tuning you can do, but it is a very standard approach and works well. There is a KFilter library available which is a C++ implementation.

My next fallback would be least squares fit. A Kalman filter will smooth the data taking velocities into account, whereas a least squares fit approach will just use positional information. Still, it is definitely simpler to implement and understand. It looks like the GNU Scientific Library may have an implementation of this.

Chris Arguin
  • 11,850
  • 4
  • 34
  • 50
  • 1
    Thanks Chris. Yes, I read about Kalman while doing some search, but it's certainly a bit beyond my math knowledge. Are you aware of any sample code easy to read (and understand!), or better yet, some implementation available? (C / C++ / Java) – Al. Jul 15 '09 at 23:12
  • 1
    @Al Unfortunately my only exposure with Kalman filters is through work, so I have some wonderfully elegant code I can't show you. – Chris Arguin Jul 15 '09 at 23:24
  • No problem :-) I tried looking but for some reason it seems this Kalman thing is black magic. Lots of theory pages but little to none code.. Thanks, will try the other methods. – Al. Jul 15 '09 at 23:28
  • 2
    http://kalman.sourceforge.net/index.php here is C++ implementation of Kalman filter. – Rostyslav Druzhchenko Aug 21 '14 at 09:04
  • @Rost Thanks, that's good to know. I added the link to the main body of the answer to make it more obvious. I'll have to play with that library now... – Chris Arguin Aug 23 '14 at 08:20
  • 1
    @ChrisArguin You are welcome. Let me know if the result is good please. – Rostyslav Druzhchenko Aug 23 '14 at 08:29
  • Is there any recommendation when some method is more suitable than the other? Like using in real time on device (with very few samples from the past) vs on server processing in background (many samples). – mauron85 Jun 15 '16 at 09:32
  • @Al. and anyone else who thinks K filters are black magic: accelerometers and gyros (inertial guidance systems) are jitter resistant but very prone to drift. GPS does not drift but is jittery. K filters apply statistical noise rejection math to a sliding window of data to integrate the GPS and inertial data into results that are neither drifty nor jittery. – Peter Wone Feb 20 '17 at 23:34
12

This might come a little late...

I wrote this KalmanLocationManager for Android, which wraps the two most common location providers, Network and GPS, kalman-filters the data, and delivers updates to a LocationListener (like the two 'real' providers).

I use it mostly to "interpolate" between readings - to receive updates (position predictions) every 100 millis for instance (instead of the maximum gps rate of one second), which gives me a better frame rate when animating my position.

Actually, it uses three kalman filters, on for each dimension: latitude, longitude and altitude. They're independent, anyway.

This makes the matrix math much easier: instead of using one 6x6 state transition matrix, I use 3 different 2x2 matrices. Actually in the code, I don't use matrices at all. Solved all equations and all values are primitives (double).

The source code is working, and there's a demo activity. Sorry for the lack of javadoc in some places, I'll catch up.

villoren
  • 319
  • 3
  • 5
  • 1
    I tried to use your lib code, I got some undesired results, I am not sure if I'm doing something wrong...(Below is image url, blue is filtered locations' path, orange are raw locations) https://app.box.com/s/w3uvaz007glp2utvgznmh8vlggvaiifk – umesh Apr 28 '15 at 13:43
  • The spikes you're seeing 'growing' from the mean (orange line) look like network provider updates. Can you try plotting both raw network and gps updates? Perhaps you would be better off with no network updates, depending on what you're trying to achieve. Btw, where are you getting those raw orange updates from? – villoren May 05 '15 at 16:19
  • 1
    orange points are from gps provider, and the blue are from Kalman, I plotted logs on map – umesh May 06 '15 at 09:30
  • Could you send me that data in some text format? Each location update has the Location.getProvider() field set. Just one file with all Location.toString(). – villoren Sep 22 '15 at 05:03
9

You should not calculate speed from position change per time. GPS may have inaccurate positions, but it has accurate speed (above 5km/h). So use the speed from GPS location stamp. And further you should not do that with course, although it works most of the times.

GPS positions, as delivered, are already Kalman filtered, you probably cannot improve, in postprocessing usually you have not the same information like the GPS chip.

You can smooth it, but this also introduces errors.

Just make sure that your remove the positions when the device stands still, this removes jumping positions, that some devices/Configurations do not remove.

Nat
  • 3,587
  • 20
  • 22
AlexWien
  • 28,470
  • 6
  • 53
  • 83
  • 5
    Could you provide some references for this please? – Matt Upson Jul 09 '15 at 10:39
  • 1
    There is many info and much professional experience in that sentences, Which sentence exactly you want a reference for? for speed: search for doppler effect and GPS. internal Kalman? This is basic GPS knowlege, every paper or book describing how a GPS chip internaly works. smootig-errors: ever smoothing introduce erros. stand still? try it out. – AlexWien Jul 09 '15 at 12:56
  • 3
    The "jumping around" when standing still is not the only source of error. There are also signal reflections (e.g. from mountains) where the position jumps around. My GPS chips (e.g. Garmin Dakota 20, SonyEricsson Neo) haven't filtered this out... And what is really a joke is the elevation value of the GPS signals when not combined with barometric pressure. These values aren't filtered or I don't want to see the unfiltered values. – hgoebl Jul 28 '15 at 20:21
  • 1
    @AlexWien GPS computes distance from a point at a time to a tolerance giving you a sphere with thickness, a _shell_ centred around a satellite. You are somewhere in this shell volume. The intersection of three of these shell volumes gives you a position volume, the centroid of which is your computed position. If you have a set of reported positions and you know the sensor is at rest, computing the centroid effectively intersects a lot more shells, improving precision. Error in this case is _reduced_. – Peter Wone Feb 20 '17 at 23:51
  • 7
    "GPS positions, as delivered, are already Kalman filtered, you probably cannot improve". If you can point to a source that confirms this for modern smartphones (for example), that would be very useful. I cannot see evidence of it myself. Even simple Kalman filtering of a device's raw locations strongly suggests it is not true. The raw locations dance around erratically, while the filtered locations most often hold close to the real (known) location. – sobri May 29 '17 at 09:34
  • 1
    @sobri - you are correct. It is possible to compute the centroid of a dataset of n+1 positions from the centroid of the first n positions and the n+1 position, and this is essentially what your solution will do. More related info in my 5yo answer below. – Peter Wone Jun 14 '18 at 11:38
  • 1
    For reference, this is a [paper mentioning the doppler effect technique for very accurate speed](https://www.researchgate.net/publication/265991687_High_accuracy_speed_measurement_using_GPS_Global_Positioning_System) and acceleration measurements. I don't know, however, if this is implemented by GPS vendors. – jjmontes Sep 29 '18 at 14:29
  • @PeterWone your comment refers to how a GPS chip calculates its position. Speed is calculaetd differently, it is not only position change over time. – AlexWien Apr 12 '21 at 16:22
7

I usually use the accelerometers. A sudden change of position in a short period implies high acceleration. If this is not reflected in accelerometer telemetry it is almost certainly due to a change in the "best three" satellites used to compute position (to which I refer as GPS teleporting).

When an asset is at rest and hopping about due to GPS teleporting, if you progressively compute the centroid you are effectively intersecting a larger and larger set of shells, improving precision.

To do this when the asset is not at rest you must estimate its likely next position and orientation based on speed, heading and linear and rotational (if you have gyros) acceleration data. This is more or less what the famous K filter does. You can get the whole thing in hardware for about $150 on an AHRS containing everything but the GPS module, and with a jack to connect one. It has its own CPU and Kalman filtering on board; the results are stable and quite good. Inertial guidance is highly resistant to jitter but drifts with time. GPS is prone to jitter but does not drift with time, they were practically made to compensate each other.

Peter Wone
  • 17,965
  • 12
  • 82
  • 134
5

One method that uses less math/theory is to sample 2, 5, 7, or 10 data points at a time and determine those which are outliers. A less accurate measure of an outlier than a Kalman Filter is to to use the following algorithm to take all pair wise distances between points and throw out the one that is furthest from the the others. Typically those values are replaced with the value closest to the outlying value you are replacing

For example

Smoothing at five sample points A, B, C, D, E

ATOTAL = SUM of distances AB AC AD AE

BTOTAL = SUM of distances AB BC BD BE

CTOTAL = SUM of distances AC BC CD CE

DTOTAL = SUM of distances DA DB DC DE

ETOTAL = SUM of distances EA EB EC DE

If BTOTAL is largest you would replace point B with D if BD = min { AB, BC, BD, BE }

This smoothing determines outliers and can be augmented by using the midpoint of BD instead of point D to smooth the positional line. Your mileage may vary and more mathematically rigorous solutions exist.

ojblass
  • 21,146
  • 22
  • 83
  • 132
  • Thanks, I'll give it a shot too. Note that I want to smooth the current position, as it's the one being displayed and the one used to retrieve some data. I'm not interested in past points. My original idea was using weighted means, but I still have to see what's best. – Al. Jul 16 '09 at 01:29
  • 1
    Al, this appears to be a form of weighted means. You will need to use "past" points if you want to do any smoothing, because the system needs to have more than the current position in order to know where to smooth too. If your GPS is taking datapoints once per second and your user looks at the screen once per five seconds, you can use 5 datapoints without him noticing! A moving average would only be delayed by one dp also. – Karl Jul 31 '09 at 06:05
4

As for least squares fit, here are a couple other things to experiment with:

  1. Just because it's least squares fit doesn't mean that it has to be linear. You can least-squares-fit a quadratic curve to the data, then this would fit a scenario in which the user is accelerating. (Note that by least squares fit I mean using the coordinates as the dependent variable and time as the independent variable.)

  2. You could also try weighting the data points based on reported accuracy. When the accuracy is low weight those data points lower.

  3. Another thing you might want to try is rather than display a single point, if the accuracy is low display a circle or something indicating the range in which the user could be based on the reported accuracy. (This is what the iPhone's built-in Google Maps application does.)

Alex319
  • 3,818
  • 9
  • 34
  • 40
4

Going back to the Kalman Filters ... I found a C implementation for a Kalman filter for GPS data here: http://github.com/lacker/ikalman I haven't tried it out yet, but it seems promising.

KarateSnowMachine
  • 1,490
  • 2
  • 13
  • 17
3

You can also use a spline. Feed in the values you have and interpolate points between your known points. Linking this with a least-squares fit, moving average or kalman filter (as mentioned in other answers) gives you the ability to calculate the points inbetween your "known" points.

Being able to interpolate the values between your knowns gives you a nice smooth transition and a /reasonable/ approximation of what data would be present if you had a higher-fidelity. http://en.wikipedia.org/wiki/Spline_interpolation

Different splines have different characteristics. The one's I've seen most commonly used are Akima and Cubic splines.

Another algorithm to consider is the Ramer-Douglas-Peucker line simplification algorithm, it is quite commonly used in the simplification of GPS data. (http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm)

Aidos
  • 2,753
  • 3
  • 27
  • 31
2

Here's a Javascript implementation of @Stochastically's Java implementation for anyone needing it:

class GPSKalmanFilter {
  constructor (decay = 3) {
    this.decay = decay
    this.variance = -1
    this.minAccuracy = 1
  }

  process (lat, lng, accuracy, timestampInMs) {
    if (accuracy < this.minAccuracy) accuracy = this.minAccuracy

    if (this.variance < 0) {
      this.timestampInMs = timestampInMs
      this.lat = lat
      this.lng = lng
      this.variance = accuracy * accuracy
    } else {
      const timeIncMs = timestampInMs - this.timestampInMs

      if (timeIncMs > 0) {
        this.variance += (timeIncMs * this.decay * this.decay) / 1000
        this.timestampInMs = timestampInMs
      }

      const _k = this.variance / (this.variance + (accuracy * accuracy))
      this.lat += _k * (lat - this.lat)
      this.lng += _k * (lng - this.lng)

      this.variance = (1 - _k) * this.variance
    }

    return [this.lng, this.lat]
  }
}

Usage example:

   const kalmanFilter = new GPSKalmanFilter()
   const updatedCoords = []

    for (let index = 0; index < coords.length; index++) {
      const { lat, lng, accuracy, timestampInMs } = coords[index]
      updatedCoords[index] = kalmanFilter.process(lat, lng, accuracy, timestampInMs)
    }
jdixon04
  • 1,435
  • 16
  • 28
1

I have transformed the Java code from @Stochastically to Kotlin

class KalmanLatLong
{
    private val MinAccuracy: Float = 1f

    private var Q_metres_per_second: Float = 0f
    private var TimeStamp_milliseconds: Long = 0
    private var lat: Double = 0.toDouble()
    private var lng: Double = 0.toDouble()
    private var variance: Float =
        0.toFloat() // P matrix.  Negative means object uninitialised.  NB: units irrelevant, as long as same units used throughout

    fun KalmanLatLong(Q_metres_per_second: Float)
    {
        this.Q_metres_per_second = Q_metres_per_second
        variance = -1f
    }

    fun get_TimeStamp(): Long { return TimeStamp_milliseconds }
    fun get_lat(): Double { return lat }
    fun get_lng(): Double { return lng }
    fun get_accuracy(): Float { return Math.sqrt(variance.toDouble()).toFloat() }

    fun SetState(lat: Double, lng: Double, accuracy: Float, TimeStamp_milliseconds: Long)
    {
        this.lat = lat
        this.lng = lng
        variance = accuracy * accuracy
        this.TimeStamp_milliseconds = TimeStamp_milliseconds
    }

    /// <summary>
    /// Kalman filter processing for lattitude and longitude
    /// https://stackoverflow.com/questions/1134579/smooth-gps-data/15657798#15657798
    /// </summary>
    /// <param name="lat_measurement_degrees">new measurement of lattidude</param>
    /// <param name="lng_measurement">new measurement of longitude</param>
    /// <param name="accuracy">measurement of 1 standard deviation error in metres</param>
    /// <param name="TimeStamp_milliseconds">time of measurement</param>
    /// <returns>new state</returns>
    fun Process(lat_measurement: Double, lng_measurement: Double, accuracy: Float, TimeStamp_milliseconds: Long)
    {
        var accuracy = accuracy
        if (accuracy < MinAccuracy) accuracy = MinAccuracy

        if (variance < 0)
        {
            // if variance < 0, object is unitialised, so initialise with current values
            this.TimeStamp_milliseconds = TimeStamp_milliseconds
            lat = lat_measurement
            lng = lng_measurement
            variance = accuracy * accuracy
        }
        else
        {
            // else apply Kalman filter methodology

            val TimeInc_milliseconds = TimeStamp_milliseconds - this.TimeStamp_milliseconds

            if (TimeInc_milliseconds > 0)
            {
                // time has moved on, so the uncertainty in the current position increases
                variance += TimeInc_milliseconds.toFloat() * Q_metres_per_second * Q_metres_per_second / 1000
                this.TimeStamp_milliseconds = TimeStamp_milliseconds
                // TO DO: USE VELOCITY INFORMATION HERE TO GET A BETTER ESTIMATE OF CURRENT POSITION
            }

            // Kalman gain matrix K = Covarariance * Inverse(Covariance + MeasurementVariance)
            // NB: because K is dimensionless, it doesn't matter that variance has different units to lat and lng
            val K = variance / (variance + accuracy * accuracy)
            // apply K
            lat += K * (lat_measurement - lat)
            lng += K * (lng_measurement - lng)
            // new Covarariance  matrix is (IdentityMatrix - K) * Covarariance
            variance = (1 - K) * variance
        }
    }
}
PeterPazmandi
  • 533
  • 10
  • 13
0

Mapped to CoffeeScript if anyones interested. **edit -> sorry using backbone too, but you get the idea.

Modified slightly to accept a beacon with attribs

{latitude: item.lat,longitude: item.lng,date: new Date(item.effective_at),accuracy: item.gps_accuracy}

MIN_ACCURACY = 1

# mapped from http://stackoverflow.com/questions/1134579/smooth-gps-data

class v.Map.BeaconFilter

  constructor: ->
    _.extend(this, Backbone.Events)

  process: (decay,beacon) ->

    accuracy     = Math.max beacon.accuracy, MIN_ACCURACY

    unless @variance?
      # if variance nil, inititalise some values
      @variance     = accuracy * accuracy
      @timestamp_ms = beacon.date.getTime();
      @lat          = beacon.latitude
      @lng          = beacon.longitude

    else

      @timestamp_ms = beacon.date.getTime() - @timestamp_ms

      if @timestamp_ms > 0
        # time has moved on, so the uncertainty in the current position increases
        @variance += @timestamp_ms * decay * decay / 1000;
        @timestamp_ms = beacon.date.getTime();

      # Kalman gain matrix K = Covarariance * Inverse(Covariance + MeasurementVariance)
      # NB: because K is dimensionless, it doesn't matter that variance has different units to lat and lng
      _k  = @variance / (@variance + accuracy * accuracy)
      @lat = _k * (beacon.latitude  - @lat)
      @lng = _k * (beacon.longitude - @lng)

      @variance = (1 - _k) * @variance

    [@lat,@lng]
lucygenik
  • 1,698
  • 2
  • 16
  • 18
  • 1
    Tried to edit this, but there's a typo in the last lines where `@lat` and `@lng` are set. Should be `+=` rather than `=` – jdixon04 Dec 14 '19 at 00:56