Home > Guides, TKJ Electronics > A practical approach to Kalman filter and how to implement it

## A practical approach to Kalman filter and how to implement it

I have for a long time been interrested in Kalman filers and how they work, I also used a Kalman filter for my Balancing robot, but I never explained how it actually was implemented. Actually I had never taken the time to sit down with a pen and a piece of paper and try to do the math by myself, so I actually did not know how it was implemented.
It turned out to be a good thing, as I actually discovered a mistake in the original code, but I will get back to that later.

I actually wrote about the Kalman filter as my master assignment in high school back in December 2011. But I only used the Kalman filter to calculate the true voltage of a DC signal modulated by known Gaussian white noise. My assignment can be found in the following zip file: http://www.tkjelectronics.dk/uploads/Kalman_SRP.zip. It is in danish, but you can properly use google translate to translate some of it. If you got any specific questions regarding the assignment, then ask in the comments below.

Okay, but back to the subject. As I sad I had never taken the time to sit down and do the math regarding the Kalman filter based on an accelerometer and a gyroscope. It was not as hard as I expected, but I must confess that I still have not studied the deeper theory behind, on why it actually works. But for me, and most people out there, I am more interrested in implementing the filter, than in the deeper theory behind and why the equations works.

Before we begin you must have some basic knowledge about matrices like multiplication of matrices and transposing of matrices. If not then please take a look at the following websites:

For those of you who do not know what a Kalman filter is, it is an algorithm which uses a series of measurements observed over time, in this context an accelerometer and a gyroscope. These measurements will contain noise that will contribute to the error of the measurement. The Kalman filter will then try to estimate the state of the system, based on the current and previous states, that tend to be more precise that than the measurements alone.

In this context the problem is that the accelerometer is in general very noise when it is used to measure the gravitational acceleration since the robot is moving back and forth. The problem with the gyro is that it drifts over time – just like a spinning wheel-gyro will start to fall down when it is losing speed.
In short you can say that you can only trust the gyroscope on a short term while you can only trust the accelerometer on a long term.

There is actually a very easy way to deal with this by using a complimentary filter, which basicly just consist of a digital low-pass filter on the accelerometer and digital high-pass filter on the gyroscope readings. But it is not as accurate as the Kalman filter, but other people have succesfully build balancing robots using a fine-tuned complimentary filter.

More information about gyroscopes, accelerometer and complimentary filters can be found in this pdf. A comparison between a complimentary filter and a Kalman filter can be found in the following blog post.

The Kalman filter operates by producing a statistically optimal estimate of the system state based upon the measurement(s). To do this it will need to know the noise of the input to the filter called the measurement noise, but also the noise of the system itself called the process noise. To do this the noise has to be Gaussian distributed and have a mean of zero, luckily for us most random noise have this characteristic.

The system state $\boldsymbol{x}_k$
The next of this article might seem very confusing for some, but I promise you if you grab a pen and a piece of paper and try to follow along it is not that hard if you are reasonable at math.

If you, like me, do not have a calculator or computer program that can work with matrices, then I recommend the free online calculator Wolfram Alpha. I used it for all the calculations in this article.

I will use the same notation as the wikipedia article, but I will like to note that when the matrixes are constants and does not depend on the current time you do not have to write the k after them. So for instance $F_k$ can be simplified to $F$.

Also I would like to write a small explanation of the other aspects of the notations.
First I will make a note about whats called the previous state:

$\boldsymbol{\hat{x}}_{k-1 | k-1}$

Which is the previous estimated state based on the previous state and the estimates of the states before it.

The next is the a priori state:

$\boldsymbol{\hat{x}}_{k | k-1}$

A priori means the estimate of the state matrix at the current time k based on the previous state of the system and the estimates of the states before it.

The last one is called a posteriori state:

$\boldsymbol{\hat{x}}_{k | k}$

Is the estimated of the state at time k given observations up to and including at time k.

The problem is that the system state itself is hidden and can only be observed through observation $z_k$. This is also called a Hidden Markov model.

This means that the state will be based upon the state at time k and all the previous states. That also means that you can not trust the estimate of the state before the Kalman filter has stabilized – take a look at the graph at the front page of my assignment.

The hat over the $\hat{x}$ means that is the estimate of the state. Unlike just a single $x$ which means the true state – the one we are trying to estimate.
So the notation for the state at time k is:

$\boldsymbol{x}_k$

The state of the system at time k if given by:

$\boldsymbol{x}_k = \boldsymbol{F}x_{k-1} + \boldsymbol{B}u_k + w_k$

Where $x_k$ is the state matrix which is given by:

$\boldsymbol{x}_k = \begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_k$

As you can see the output of the filter will be the angle $\theta$ but also the bias $\dot{\theta}_b$ based upon the measurements from the accelerometer and gyroscope. The bias is the amount the gyro has drifted. This means that one can get the true rate by subtracting the bias from the gyro measurement.

The next is the $F$ matrix, which is the state transition model which is applied to the prevouis state $x_{k-1}$.

In this case $F$ is defined as:

$\boldsymbol{F} = \begin{bmatrix} 1 & -\Delta t \\ 0 & 1 \end{bmatrix}$

I know that the $-\Delta t$ might seem confusing, but it will make sense later (take a look at my comment).

The next is the control input $u_k$, in this case it is the gyroscope measurement in degrees per second (°/s) at time k, this is also called the rate $\dot{\theta}$. We will actually rewrite the state equation as:

$\boldsymbol{x}_k = \boldsymbol{F}x_{k-1} + \boldsymbol{B}{\dot{\theta}_k} + w_k$

The next thing is the $B$ matrix. Which is called the control-input model, which is defined as:

$\boldsymbol{B} = \begin{bmatrix} \Delta t \\ 0 \end{bmatrix}$

This makes perfectly sense as you will get the angle $\theta$ when you multiply the rate $\dot{\theta}$ by the delta time $\Delta t$ and since we can not calculate the bias directly based on the rate we will set the bottom of the matrix to 0.

$w_k$ is process noise which is Gaussian distributed with a zero mean and with covariance $Q$ to the time k:

$\boldsymbol{w}_k \sim N\left ( 0, \boldsymbol{Q}_k \right )$

$Q_k$ is the process noise covariance matrix and in this case the covariance matrix of the state estimate of the accelerometer and bias. In this case we will consider the estimate of the bias and the accelerometer to be independent, so it’s actually just equal to the variance of the estimate of the accelerometer and bias.
The final matrix is defined as so:

$\boldsymbol{Q}_k = \begin{bmatrix} Q_\theta & 0 \\ 0 & Q_{\dot{\theta}_b} \end{bmatrix}\Delta t$

As you can see the $Q_k$ covariance matrix depends on the current time k, so the accelerometer variance $Q_\theta$ and the variance of the bias $Q_{\dot{\theta}_b}$ is multiplied by the delta time $\Delta t$.
This makes sense as the process noise will be larger as longer time it is since the last update of the state. For instance the gyro could have drifted.
You will have to know these constants for the Kalman filter to work.
Note if you set a larger value, the more noise in the estimation of the state. So for instance if the estimated angle starts to drift you have to increase the value of $Q_{\dot{\theta}_b}$. Otherwise if the estimate tends to be slow you are trusting the estimate of the angle too much and should try to decrease the value of $Q_\theta$ to make it more responsive.

The measurement $\boldsymbol{z}_k$
Now we will take a look at the observation or measurement $z_k$ of the true state $x_k$. The observation $z_k$ is given by:

$\boldsymbol{z}_k = \boldsymbol{H}x_{k} + v_k$

As you can see the measurement $z_k$ is given by the current state $x_k$ multiplied by the $H$ matrix plus the measurement noise $v_k$.

$H$ is called the observation model and is used to map the true state space into the observed space. The true state can not be observed. Since the measurement is just the measurement from the accelerometer, $H$ is given by:

$\boldsymbol{H} = \begin{bmatrix} 1 & 0 \end{bmatrix}$

The noise of the measurement have to be Gaussian distributed as well with a zero mean and $R$ as the covariance:

$\boldsymbol{v}_k \sim N\left ( 0, \boldsymbol{R} \right )$

But as $R$ is not a matrix the measurement noise is just equal to the variance of the measurement, since the covariance of the same variable is equal to the variance. See this page for more information.
Now we can define $R$ as so:

$\boldsymbol{R} = E \begin{bmatrix} v_k & {v_k}^T \end{bmatrix} = var(v_k)$

We will assume that the measurement noise is the same and does not depend on the time k:

$var(v_k) = var(v)$

Note that if you set the measurement noise variance $var(v)$ too high the filter will respond really slowly as it is trusting new measurements less, but if it is too small the value might overshoot and be noisy since we trust the accelerometer measurements too much.

So to round up you have to find the the process noise variances $Q_\theta$ and $Q_{\dot{\theta}_b}$ and the measurement variance of the measurement noise $var(v)$. There are multiple ways to find them, but it is out of the aspect of this article.

The Kalman filter equations
Okay now to the equations we will use to estimate the true state of the system at time k $\hat{x}_k$. Some clever guys came up with equations found below to estimate the state of the system.
The equations can be written more compact, but I prefer to have them stretched out, so it is easier to implement and understand the different steps.

Predict
In the first two equations we will try to predict the current state and the error covariance matrix at time k. First the filter will try to estimate the current state based on all the previous states and the gyro measurement:

$\boldsymbol{\hat{x}}_{k | k-1} = \boldsymbol{F}\hat{x}_{k-1 | k-1} + \boldsymbol{B}{\dot{\theta}_k}$

That is also why it is called a control input, since we use it as an extra input to estimate the state at the current time k called the a priori state $\hat{x}_{k | k-1}$ as described in the beginning of the article.

The next thing is that we will try to estimate the a priori error covariance matrix $P_{k | k-1}$ based on the previous error covariance matrix $P_{k-1 | k-1}$, which is defined as:

$\boldsymbol{P}_{k | k-1} = \boldsymbol{F}\boldsymbol{P}_{k-1 | k-1}\boldsymbol{F}^T + \boldsymbol{Q}_k$

This matrix is used to estimate how much we trust the current values of the estimated state. The smaller the more we trust the current estimated state. The principle of the equation above is actually pretty easy to understand, as it is pretty obvious that the error covariance will increase since we last updated the estimate of the state, therefore we multiplied the error covariance matrix by the state transition model $F$ and the transpose of that $F^T$ and add the current process noise $Q_k$ at time k.

The error covariance matrix $P$ in our case is a 2×2 matrix:

$\boldsymbol{P} = \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}$

Update
The fist thing we will do is to compute the difference between the measurement $z_k$ and the a priori state $x_{k | k-1}$, this is also called the innovation:

$\boldsymbol{\tilde{y}}_k = \boldsymbol{z}_k - \boldsymbol{H}\hat{x}_{k | k-1}$

The observation model $H$ is used to map the a priori state $x_{k | k-1}$ into the observed space which is the measurement from the accelerometer, therefore the innovation is not a matrix

$\boldsymbol{\tilde{y}}_k = \begin{bmatrix} \boldsymbol{\tilde{y}} \end{bmatrix}_k$

The next thing we will do is calculate whats called the innovation covariance:

$\boldsymbol{S}_k = \boldsymbol{H} \boldsymbol{P}_{k | k-1} \boldsymbol{H}^T + \boldsymbol{R}$

What it does is that it tries to predict how much we should trust the measurement based on the a priori error covariance matrix $P_{k | k-1}$ and the measurement covariance matrix $R$. The observation model $H$ is used to map the a priori error covariance matrix $P_{k | k-1}$ into observed space.

The bigger the value of the measurement noise the larger the value of $S$, this means that we do not trust the incoming measurement that much.
In this case $S$ is not a matrix and is just written as:

$\boldsymbol{S}_k = \begin{bmatrix} \boldsymbol{S} \end{bmatrix}_k$

The next step is to calculate the Kalman gain. The Kalman gain is used to to indicate how much we trust the innovation and is defined as:

$\boldsymbol{K}_k = \boldsymbol{P}_{k | k-1} \boldsymbol{H}^T \boldsymbol{S}_k^{-1}$

You can see that if we do not trust the innovation that much the innovation covariance $S$ will be high and if we trust the estimate of the state then the error covariance matrix $P$ will be small the Kalman gain will therefore be small and oppesite if we trust the innovation but does not trust the estimation of the current state.

If you take a deeper look you can see that the transpose of the observation model $H$ is used to map the state of the error covariance matrix $P$ into observed space. We then compare the error covariance matrix by multiplying with the inverse of the innovation covariance $S$.

This make sense as we will use the observation model $H$ to extract data from the state error covariance and compare that with the current estimate of the innovation covariance.

Note that if you do not know the state at startup you can set the error covariance matrix like so:

$\boldsymbol{P} = \begin{bmatrix} L & 0 \\ 0 & L \end{bmatrix}$

Where $L$ represent a large number.

For my balancing robot I know the starting angle and I find the bias of the gyro at startup by calibrating, so I assume that the state will be known at startup, so I initialize the error covariance matrix like so:

$\boldsymbol{P} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}$

In this case the Kalman gain is a 2×1 matrix:

$\boldsymbol{K} = \begin{bmatrix} K_0 \\ K_1 \end{bmatrix}$

Now we can update the a posteriori estimate of the current state:

$\boldsymbol{\hat{x}}_{k | k} = \boldsymbol{\hat{x}}_{k | k-1} + \boldsymbol{K}_k \; \boldsymbol{\tilde{y}}_k$

This is done by adding the a priori state $\hat{x}_{k | k-1}$ with the Kalman gain multiplied by the innovation $\tilde{y}_k$.
Remember that the innovation $\tilde{y}_k$ is the difference between the measurement $z_k$ and the estimated priori state $\hat{x}_{k | k-1}$, so the innovation can both be positive and negative.

A little simplified the equation can be understood as we simply correct the estimate of the a priori state $\hat{x}_{k | k-1}$, that was calculated using the previous state and the gyro measurement, with the measurement – in this case the accelerometer.

The last thing we will do is update the a posteriori error covariance matrix:

$\boldsymbol{P}_{k | k} = (\boldsymbol{I} - \boldsymbol{K}_k \boldsymbol{H}) \boldsymbol{P}_{k | k-1}$

Where $I$ is called the identity matrix and is defined as:

$\boldsymbol{I} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$

What the filter is doing is that it is basically self-correcting the error covariance matrix based on how much we corrected the estimate. This make sense as we corrected the state based the a priori error covariance matrix $P_{k | k-1}$, but also the innovation covariance $S_k$.

Implementing the filter
In this section I will use the equation from above to implement the filter into a simple c++ code that can be used for balancing robots, quadcopters and other applications where you need to compute the angle, bias or rate.

In case you want the code next to you, it can be found at github: https://github.com/TKJElectronics/KalmanFilter.

I will simply write the equations at the top of each step and then simplify them after that I will write how it is can be done i C and finally I will link to calculations done in Wolfram Alpha in the bottom of each step, as I used them to do the calculation.

Step 1:

$\boldsymbol{\hat{x}}_{k | k-1} = \boldsymbol{F}\hat{x}_{k-1 | k-1} + \boldsymbol{B}{\dot{\theta}_k}$

$\begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k | k-1} = \begin{bmatrix} 1 & -\Delta t \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k-1 | k-1} + \begin{bmatrix} \Delta t \\ 0 \end{bmatrix} \dot{\theta}_k$

$= \begin{bmatrix} \theta - \dot{\theta}_b \Delta t \\ \dot{\theta}_b \end{bmatrix}_{k-1 | k-1} + \begin{bmatrix} \Delta t \\ 0 \end{bmatrix} \dot{\theta}_k$

$= \begin{bmatrix} \theta - \dot{\theta}_b \Delta t + \dot{\theta} \Delta t \\ \dot{\theta}_b \end{bmatrix}$

$= \begin{bmatrix} \theta + \Delta t (\dot{\theta} - \dot{\theta}_b) \\ \dot{\theta}_b \end{bmatrix}$

As you can see the a priori estimate of the angle is $\hat{\theta}_{k | k-1}$ is equal to the estimate of the previous state $\hat{\theta}_{k-1 | k-1}$ plus the unbiased rate times the delta time $\Delta t$.
Since we can not directly measure the bias the estimate of the a priori bias is just equal to the previous one.

This can be written in C like so:

rate = newRate - bias;
angle += dt * rate;

Note that I calculate the unbiased rate, so it can be be used by the user as well.

Step 2:

$\boldsymbol{P}_{k | k-1} = \boldsymbol{F}\boldsymbol{P}_{k-1 | k-1}\boldsymbol{F}^T + \boldsymbol{Q}_k$

$\begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1} = \begin{bmatrix} 1 & -\Delta t \\ 0 & 1 \end{bmatrix} \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k-1 | k-1} \begin{bmatrix} 1 & 0 \\ -\Delta t & 1 \end{bmatrix} + \begin{bmatrix} Q_\theta & 0 \\ 0 & Q_{\dot{\theta}_b} \end{bmatrix}\Delta t$

$= \begin{bmatrix} P_{00} - \Delta t P_{10} & P_{01} - \Delta t P_{11} \\ P_{10} & P_{11} \end{bmatrix}_{k-1 | k-1} \begin{bmatrix} 1 & 0 \\ -\Delta t & 1 \end{bmatrix} + \begin{bmatrix} Q_\theta & 0 \\ 0 & Q_{\dot{\theta}_b} \end{bmatrix}\Delta t$

$= \begin{bmatrix} P_{00} -\Delta t P_{10} -\Delta t (P_{01} - \Delta t P_{11}) & P_{01} -\Delta t P_{11} \\ P_{10} -\Delta t P_{11} & P_{11} \end{bmatrix}_{k-1 | k-1} + \begin{bmatrix} Q_\theta & 0 \\ 0 & Q_{\dot{\theta}_b} \end{bmatrix}\Delta t$

$= \begin{bmatrix} P_{00} -\Delta t P_{10} -\Delta t (P_{01} - \Delta t P_{11}) + Q_\theta \Delta t & P_{01} -\Delta t P_{11} \\ P_{10} -\Delta t P_{11} & P_{11} + Q_{\dot{\theta}_b} \Delta t \end{bmatrix}$

$= \begin{bmatrix} P_{00} + \Delta t (\Delta t P_{11} - P_{01} - P_{10} + Q_\theta) & P_{01} -\Delta t P_{11} \\ P_{10} -\Delta t P_{11} & P_{11} + Q_{\dot{\theta}_b} \Delta t \end{bmatrix}$

The equations above can be written in C like so:

P[0][0] += dt * (dt*P[1][1] - P[0][1] - P[1][0] + Q_angle);
P[0][1] -= dt * P[1][1];
P[1][0] -= dt * P[1][1];
P[1][1] += Q_gyroBias * dt;

Note that this is the part of the code that there was an error in in the original code that I used.

Step 3:

$\boldsymbol{\tilde{y}}_k = \boldsymbol{z}_k - \boldsymbol{H}\hat{x}_{k | k-1}$

$= \boldsymbol{z}_k - \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k | k-1}$

$= \boldsymbol{z}_k - \theta_{k | k-1}$

The innovation can be calculated in C like so:

y = newAngle - angle;

Step 4:

$\boldsymbol{S}_k = \boldsymbol{H} \boldsymbol{P}_{k | k-1} \boldsymbol{H}^T + \boldsymbol{R}$

$= \begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1} \begin{bmatrix} 1 \\ 0 \end{bmatrix} + \boldsymbol{R}$

$= {P_{00}}_{k | k-1} + \boldsymbol{R}$

$= {P_{00}}_{k | k-1} + var(v)$

Again the C code is pretty simple:

S = P[0][0] + R_angle;

Step 5:

$\boldsymbol{K}_k = \boldsymbol{P}_{k | k-1} \boldsymbol{H}^T \boldsymbol{S}_k^{-1}$

$\begin{bmatrix} K_0 \\ K_1 \end{bmatrix}_k = \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \boldsymbol{S}_k^{-1}$

$= \begin{bmatrix} P_{00} \\ P_{10} \end{bmatrix}_{k | k-1} \boldsymbol{S}_k^{-1}$

$= \dfrac{\begin{bmatrix} P_{00} \\ P_{10}\end{bmatrix}_{k | k-1}}{\boldsymbol{S}_k}$

Note that in other cases $S$ can be a matrix and you can not just simply divide $P$ by $S$. Instead you have to calculate the inverse of the matrix. See the following page for more information on how to do so.

The C implementation looks like this:

K[0] = P[0][0] / S;
K[1] = P[1][0] / S;

Step 6:

$\boldsymbol{\hat{x}}_{k | k} = \boldsymbol{\hat{x}}_{k | k-1} + \boldsymbol{K}_k \; \boldsymbol{\tilde{y}}_k$

$\begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k | k} = \begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k | k-1} + \begin{bmatrix} K_0 \\ K_1 \end{bmatrix}_k \boldsymbol{\tilde{y}}_k$

$= \begin{bmatrix} \theta \\ \dot{\theta}_b \end{bmatrix}_{k | k-1} + \begin{bmatrix} K_0 \; \boldsymbol{\tilde{y}} \\ K_1 \; \boldsymbol{\tilde{y}} \end{bmatrix}_k$

Yet again the equation end up pretty short, and can be written as so in C:

angle += K[0] * y;
bias += K[1] * y;

Step 7:

$\boldsymbol{P}_{k | k} = (\boldsymbol{I} - \boldsymbol{K}_k \boldsymbol{H}) \boldsymbol{P}_{k | k-1}$

$\begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k} = \left(\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} - \begin{bmatrix} K_0 \\ K_1 \end{bmatrix}_k \begin{bmatrix} 1 & 0 \end{bmatrix}\right) \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1}$

$= \left(\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} - \begin{bmatrix} K_0 & 0 \\ K_1 & 0 \end{bmatrix}_k \right) \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1}$

$= \begin{bmatrix} P_{00} & P_{01} \\ P_{10} & P_{11} \end{bmatrix}_{k | k-1} - \begin{bmatrix} K_0 \; P_{00} & K_0 \; P_{01} \\ K_1 \; P_{00} & K_1 \; P_{01} \end{bmatrix}$

Remember that we decrease the error covariance matrix again, since the error of the estimate of the state has been decreased.
The C code looks like this:

P[0][0] -= K[0] * P[0][0];
P[0][1] -= K[0] * P[0][1];
P[1][0] -= K[1] * P[0][0];
P[1][1] -= K[1] * P[0][1];

Note that I have found that the following variances works perfectly for most IMUs:

const double Q_angle = 0.001;
const double Q_gyroBias = 0.003;
const double R_angle = 0.03;

Remember that it’s very important to set the target angle at startup if you need to use the output at startup. For more information, see the calibration routine for my balancing robot.

In case you missed it here is the library I wrote a library that can be used by any microcontroller that supports floating math. The source code can be found at github: https://github.com/TKJElectronics/KalmanFilter.

If you prefer a video explanation about the Kalman filter, I recommend the following video series: http://www.youtube.com/watch?v=FkCT_LV9Syk.

Note that you can not use the library if you need to represent something in a full 3D orientations, as euler angles suffer from what is called Gimbal lock you will need to use Quaternions to do that, but that is a whole nother story. For now take a look at the following page.

This is all for know, I hope that you will find i helpfull, if you do or have any questions fell free to post a comment below – it supports LaTeX syntax as well, if you need to write equations.
If you spot any errors please let me know as well.

Categories: Tags:
1. July 30th, 2013 at 23:56 | #1

@Krishna
Thank you.
Try to implement a low pass filter before the Kalman filter and see if it helps.

2. July 31st, 2013 at 00:21 | #2

@Lauszus
Hi,

Thanks for the reply but yes I have tried a low pass filter it helps a little but it slows down the response very much.

3. July 31st, 2013 at 00:48 | #3

@Krishna
Okay. What application are you using it for?

4. July 31st, 2013 at 11:21 | #4

@Lauszus
I am using it in a car to obtain pitch and roll.

5. August 1st, 2013 at 13:05 | #5

@Krishna
Okay. You might have to use an magnetometer too if you can’t find a compromise between noise and latency.

6. August 1st, 2013 at 23:22 | #6

@Lauszus
Magnetometer for pitch? I always thought they were for yaw. Ok I will look into it. Thanks for your help.

7. August 2nd, 2013 at 13:44 | #7

@krishna
I’m pretty sure it’s possible with a more advanced algorithm like the extended Kalman filter or a DCM algorithm. You need a gyroscope too of course.

8. August 2nd, 2013 at 23:44 | #8

Yes I was also thinking about extended kalman filter, but I am still wondering how the magnetometer works to measure pitch. I know it can me used for yaw as it gives heading with respect to North.

9. August 3rd, 2013 at 21:57 | #9

@krishna
Sorry I don’t know what I was thinking about. The magnetometer is used for yaw.

10. August 4th, 2013 at 23:18 | #10

hey lauszus, man thank you very much, this article was like painkiller for me. i understand kalman much more now.
but i need a little help with the code.
when it comes to combining the bytes to doubles, i didnt understand the logic. i found an explanation in the code like this;
// Gyro format is MSB first
gyroXrate = buffer[0] << 8 | buffer[1];
………
// Accel is LSB first.
accelX = buffer[7] << 8 | buffer[6];
……..

how do i know which byte is first and what does it exactly mean?

11. August 7th, 2013 at 10:46 | #11

Ok no problem. So how can this problem be solved? I tried the low pass filter but I have to use something with less than 0dB gain at low frequency to make it work. This is not traditionally a LPF as it is supposed to allow everything unaltered at lower frequencies. Any ideas?

Thanks

12. August 16th, 2013 at 20:55 | #12

@yusuf
First of all, sorry for the long delay!
It’s not converted to doubles, but words (16-bit).
For an explanation of MSB and LSB see the following pages: http://en.wikipedia.org/wiki/Most_significant_bit and http://en.wikipedia.org/wiki/Least_significant_bit.

@krishna
Sorry for the long response time to you too!
I think you need some more advanced algorithm. Try to google DCM or extended Kalman filter.

13. August 16th, 2013 at 22:37 | #13

Thanks for the help. I will look into it.

14. August 22nd, 2013 at 15:05 | #14

Hellow, thank you for your article. But I am a little confused with your measurement model? I wonder how you related the accelerometer reading with angle information. It seems to me that Z(k) denotes accelerometer reading, X(k) is angle information at state(k), How did you related them use H=1 by Z(k) = H*X(K)? thanks a lot!

15. August 22nd, 2013 at 22:59 | #15

@YANG
That is because I convert the accelerometer reading into degrees before it’s used as a input to the filter. See the following lines in the example code: https://github.com/TKJElectronics/KalmanFilter/blob/master/examples/MPU6050/MPU6050.ino#L89-L90.

16. August 26th, 2013 at 12:32 | #16

Hi,
Thank you for a very good kalman tutorial.

I have questions from your example.
https://github.com/TKJElectronics/KalmanFilter/blob/master/examples/MPU6050/MPU6050.ino

1. Why you use accY for finding accXangle and accX for accYangle.(Line : 89-90)

why not use accX for accXangle and accY for accYangle ?

2. Again from the link above on line 99-100.
Why you use gyroXrate to find compAngleX and gyroYrate to find compAngleY?

If we look in datasheet Acc X-axis will rotate on gyro Y-axis so we should use gyroXrate for finding compAngleY and gyroYrate for compAngleX.

Thank you.

17. August 26th, 2013 at 13:18 | #17

@Beran
Hi Beran.
I’ll give you some more details about the case and why you have been confused.

Have a look at the picture above. The 6-axis sensor, the MPU6050 fx, contains a 3-axis accelerometer that measures the acceleration i 3 different directions, X, Y and Z. In our case the acceleration is the gravity from the earth.
The other 3-axis is the gyroscope that measures the rotation around the axis. So when we start tilting the robot, the gyro values will increase and display the speed that we are rotating the robot with.

So to actually answer both of your questions, as they deal with the same confusion, look at the picture and think what happens to the values of the different axis when you tilt it around the x-axis.
By tilting the device around the x-axis, the gyroscope x-value will change and indicate the rotational speed.
The accelerometer x-axis however will stay stable, as that is still pointing out to you, as you are only tilting it around this axis.
Instead both the z-axis and y-axis values of the accelerometer will change according to the impact of the gravity on both of them.

I hope this clearified both of your questions, and you now understand why we had to chose and use the values from the axis as we did.

Kind regards
Thomas Jespersen

18. August 27th, 2013 at 05:09 | #18

Hi thank you.

http://s10.postimg.org/7ekox96h5/Snap2.jpg

From MPU6050 picture in link above, if AccX (+X axis) change +10 degrees from horizontal which gyro changed?
I think Gyro-Y should be changed? so why we not use gyroYrate in line 99?.

Here is your code line 99 :
compAngleX = (0.93*(compAngleX+(gyroXrate*(double)(micros()-timer)/1000000)))+(0.07*accXangle);

Thank you again.

19. August 27th, 2013 at 05:18 | #19

You code should be used with this board?
https://www.sparkfun.com/products/retired/10010

not for MPU-6050

sorry if I misunderstood.

20. August 30th, 2013 at 09:31 | #20

This is a very good article and first of all ı want to thank you so much.
ım new in this topic so if my question is nonsensical please forgive me .
what is the physical meaning of P00 or p01 … ? in other words if , for example , p00 is 5 or 10 , what is the meaning of this numbers ?

thanks you.

21. August 30th, 2013 at 09:47 | #21

Şükrü :
This is a very good article and first of all ı want to thank you so much.
ım new in this topic so if my question is nonsensical please forgive me .
what is the physical meaning of P00 or p01 … ? in other words if , for example , p00 is 5 or 10 , what is the meaning of this numbers ? and again ; what is the meaning of the kalman gain ? in your kalman code , you write P00 -= k0*P00 so if ı can understand this equations’s physical meaning , ı can understand this topic
thanks you.

22. August 30th, 2013 at 17:28 | #22

@Beran
Please read Thomas’ explanation again. Remember that the accelerometer measures the acceleration in a certain direction (think of it as vectors), while the gyro measures the rotational speed (measured in deg/s) around an axis.

I originally posted the code for the LPR530AL, LY530ALH and ADXL335 combo, but then afterward I create a similar code for the MPU-6050. See: https://github.com/TKJElectronics/Example-Sketch-for-IMU-including-Kalman-filter/tree/master/IMU6DOF/MPU6050.

@Şükrü
Those values are part of the error covariance matrix. See the text after the headline “Predict”. It explains what it is.

23. September 3rd, 2013 at 03:52 | #23

hi, Lauszus, I am sorry to tell you that all pictures in your artical can’t be seen at this moment. is that caused by github parsers or some reason else?

24. September 3rd, 2013 at 08:57 | #24

@bow
I can see the pictures just fine. Are you still affected by the problem?
Btw they are not hosted on Github, but are written in Latex and then automatically generated by a WordPress plugin called “WP LaTeX”

25. September 5th, 2013 at 15:59 | #25

Hello everyone, and thanks to Lauszus for this great tutorial! I got a question:

I´m using the Microchip Motion Sense board, it has the MPU6050 and AK8975 magnetometer. Microchip also gives libraries to use this board and obtain quaternion data (in 3D). I think is possible with some operations convert this quaternion data into yaw, pitch and roll angles, to use them in a drone.

But for this dynamic systems is necessary for example a complementary or Kalman filter to obtain a good solution, and Lauszus post in a good explain of this. But i´m a bit confuse… MPU6050 can work in combination with the magnetometer and Microchip did it for this board, i ask them to get the angles and answer me to filter accel (complementary filter)… So at the end i think i need to read the gyro, accel and magnetometer and fusion in a kalman filter, maybe extended kalman???

I´ve readed many info and i think i´m a bit mixed! Also Lauszus post a link in another forum for the kalman implementation in C # but link is dead, would be great have this info for C#.

Thanks a lot!

26. September 8th, 2013 at 12:56 | #26

@Andrew
You will have to use the magnetometer and gyro to estimate the yaw. Then you simply create another instance of the library and use the magnetometer angle as the first argument instead of the accelerometer.
But if you need full 3D measurement, then it would be better to use a extended Kalman filter or a DCM algorithm.

27. September 9th, 2013 at 10:24 | #27

I can´t remember now, you answered to another person something about the extended Kalman. Do you have some info about this extended version to start? Sorry but i don´t understand what you mean about the accelerometer…

Yaw -> Magnetometer and gyro.
Pitch and roll?

Thanks again Lauszus.

28. September 9th, 2013 at 15:33 | #28

@Andrew
No I don’t, sorry.
Use the magnetometer and gyro for yaw and then gyro and accelerometer for pitch and roll.

29. September 11th, 2013 at 05:21 | #29

https://github.com/TKJElectronics/Example-Sketch-for-IMU-including-Kalman-filter/blob/master
/IMU6DOF/MPU6050/MPU6050.ino#L89

Is the accXangle in the link above is the angle of projection of the X on YZ plane to Z-axis?
Any images are appreciated.

Thank you.

30. September 11th, 2013 at 12:57 | #30
31. October 4th, 2013 at 23:57 | #31

Hi,
first of all, I want to thank you so much for your article. it’s really helpful.
So, that’s my question: In this example:
https://github.com/TKJElectronics/KalmanFilter/blob/master/examples/MPU6050/MPU6050.ino
You calculate accXangle and accYangle that are respectively rotation about x-axis and Y-axis. But what about accZangle? can we get it (for example) from the accelerometer output only as accZangle = (atan2(accX,accY)+PI)*RAD_TO_DEG?

thank you

32. October 5th, 2013 at 11:11 | #32

@Yvan
No it is not possible to estimate yaw using an accelerometer. Remember that each axis all measure the gravitational force in their respective axis -- see them as vectors.
When you rotate the IMU along the z-axis none of the axis will change, as the gravitational force doesn’t change in any of the axis.
You will need a gyro and a magnetometer to estimate yaw.

33. October 9th, 2013 at 03:54 | #33

Thanks a lot for this guide Lauszus! It is incredibly well explained and the Arduino library is terrific as well! I managed to get it working with Pololu minIMU v2 even though I have almost no idea about different gyro and compass (they refer to the accelerometer and magnetometer as a compass) set-up modes -- in terms of refresh rates, sensitivity and so on. I will try implementing yaw as well using the magnetometer and gyro as you’ve explained in the comments before, but because I need it for a 3D application I’ll have to try using extended Kalman. I also wanted to mention that there is sample Arduino code for Pololu minIMU that uses DCM and there is a vPython program that shows 3D real time animation of the sensor -- maybe you or someone reading this might be interested.

Bottom line -- AMAZING GUIDE -- you belong in MIT!!!

34. October 9th, 2013 at 11:59 | #34

@Peter
Glad you found my guide useful
I might look into the extended Kalman filter in the future and DCM, but I’m too busy right now to dig into it!
Haha yeah maybe. I’m actually thinking of studying abroad next year, so I’m thinking of apply for MIT and see if I can get in.

35. October 10th, 2013 at 01:51 | #35

however i have another little problem. I tried to implement that code https://github.com/TKJElectronics/KalmanFilter/blob/master/examples/MPU6050/MPU6050.ino
it works but i get i an overshoot that i can’t disregard. for example, if i try to get 180°, i would get about 190°, then go slowly to 180°. i don’t know how i can fix it
Thanks, best regards

36. October 18th, 2013 at 21:59 | #36

@Yvan
It sounds like you are trusting the measurement of the accelerometer too much.

See the following part of the guide:

“Note that if you set the measurement noise variance var(v) too high the filter will respond really slowly as it is trusting new measurements less, but if it is too small the value might overshoot and be noisy since we trust the accelerometer measurements too much.”

37. October 30th, 2013 at 10:20 | #37

Hi, Thanks for the code on GitHub, I use it for my quadcopter controller (with a MPU 6050). I have got one question though. On this file: https://github.com/TKJElectronics/Example-Sketch-for-IMU-including-Kalman-filter/blob/master/IMU6DOF/MPU6050/MPU6050.ino
You set Gyro Full Scale Range to ±250deg/s. And after you write:
double gyroXrate = (double)gyroX/131.0;
I don’t see where that “131.0″ comes from?
I mean it works but I’d like to understand why
Thanks, Romain

38. October 30th, 2013 at 15:46 | #38

@Romain
See the Register Map datasheet: http://invensense.com/mems/gyro/documents/RM-MPU-6000A-00v4.2.pdf page 14 you can see that by writing a 0 to the FS_SEL register the range is set to ±250deg/s.

After that go see the Product Specification: http://invensense.com/mems/gyro/documents/PS-MPU-6000A-00v3.4.pdf at page 12 -- there you will see if the FS_SEL register is set to 0, then sensivity is: 131 LSB/(º/s)

39. November 12th, 2013 at 07:29 | #39

Just a thanks for taking the time to explain this. I hope I can use some of the knowledge when I try to build a robot.

40. November 12th, 2013 at 14:08 | #40

@Graham Stott
You are welcome

41. December 6th, 2013 at 16:22 | #41

Thank Lauszus so much, it’s helpful for me. You’re very nice.

Comment pages
1 2 3 2868
1. September 10th, 2012 at 19:47 | #1
2. September 11th, 2012 at 06:08 | #2
3. September 12th, 2012 at 01:10 | #3