1
$\begingroup$

I am relatively new to using Kalman Filtering. Currently I am trying to understand it and how to implement it in Matlab. I found a website with some nice examples that I would like to rewrite in Matlab using the unscentedKalmanFilter() function. The website is Kalman Filter examples and I am trying to rebuild the first example where depending on some measurments with noise the weight of a gold bar is estimated.

I've come up with some code but I doubt that my state transistion function $f$ is correct. From what I understand one needs the measurement function $y=h(x)$ calculating the state $x$ from the measurement $y$ and the state transition function $x_{n+1}=f(x_n,u)$ where $x_{n+1}$ is the prediction for the next state from the input $u$ and the current state $x_{n}$.

As the gold bar weight is directly measured the measurement function should simply be $h(x)=x+v_{k}$ with $v_k$ being the additive noise term. For the state transition function I don't have an idea. In the example they write $\hat x_{n}=\hat x_{n-1}+α_n(y_n−\hat x_{n-1})$ where $y_n$ is the last measurement and $α_n$ a factor. However I don't think the state transition function $f$ should include the measurement.

So maybe someone could enlighten me. Here is the Matlab code I have so far:

initialStateGuess = [1030];
n=length(initialStateGuess);

h=@(x)[x(1)];% Measurement function
f=@(x)[x(1)];% State transition function

ukf = unscentedKalmanFilter(...
    f,... % State transition function
    h,... % Measurement function
    initialStateGuess,...
    'HasAdditiveMeasurementNoise',true,...
    'HasAdditiveProcessNoise',true);

yMeas = [1030 989 1017 1009 1013 979 1008 1042 1012 1011];% Measurement Values
R = var(yMeas);%0.2; % Variance of the measurement noise v[k]
ukf.MeasurementNoise = R;
ukf.ProcessNoise = 60;

Nsteps = length(yMeas); % Number of time steps
xCorrectedUKF = zeros(Nsteps,n); % Corrected state estimates
PCorrected = zeros(Nsteps,n,n); % Corrected state estimation error covariances
e = zeros(Nsteps,1); % Residuals (or innovations)

for k=1:Nsteps
    % Let k denote the current time.
    %
    % Residuals (or innovations): Measured output - Predicted output
    e(k) = yMeas(k) - vdpMeasurementFcn(ukf.State); % ukf.State is x[k|k-1] at this point
    [xCorrectedUKF(k,:), PCorrected(k,:,:)] = correct(ukf,yMeas(k));
    predict(ukf);
end

plot(xCorrectedUKF')
hold on
plot(yMeas)
plot(1010*ones(1,Nsteps))
legend('Kalman FIlter state','Measurement Values','True Value')

10 steps of the Kalman Filter state implemented in my Matlab code

Edit:

As can be seen in the above picture the Matlab code I wrote does work. As more measurements come in the state gets closer to the true value. However the variance value for the Process Noise seems to be quite relevant to achieve a fast convergence. Even after playing around with it for a bit I don't manage to achieve that fast convergence as in the example:

10 steps to the Kalman Filter from the example source

So I wonder why my code convergeces so slowliy?

Edit2:

As suggested I have tried to set the process noise to 0 and the measurement noise to the variance of the measurement values. However this results in a very slow settling time. Might this have to da with me using an unscented Kalman filter? I thought that would be a good choise as I can use it also for other, nonlinear problems. ProcessNoise0 MeasurementNoise320

$\endgroup$
8
  • $\begingroup$ I didn't look at the example, but to me the equation containing $\alpha_n$ seems to be both the prediction and correction step combined in one equation. So the model is $x_n=x_{n-1}$ (maybe with some added noise). The $\alpha_n$ would then the correction gain obtained from the Kalman filter. $\endgroup$ – fibonatic Mar 4 at 21:22
  • $\begingroup$ They write $\hat x_{n}=\hat x_{n-1}+α_n(y_n − \hat x_{n-1})$. $\hat x$. The estimate of $x$. Not the actual $x$. Does it make sense now? $\endgroup$ – TimWescott Mar 5 at 1:55
  • $\begingroup$ @fibonatic: Thank you for your reply. Ok in that case I did model it corectly in Matlab? I assumed $h: x_n=y_n$ and $f: x_{n+1}=x_n$. Both equations have a noise term added to them. Then I wonder why my model convergece worse than the one in the example? Maybe I have to chose my Measurment and Process Noise Variance differently? $\endgroup$ – Matthias La Mar 5 at 10:54
  • $\begingroup$ @TimWescott: Thank you! Yes, I indeed forgot to add the $^$. However I still don't fully get it. Is my modelling in Matlab correct then with the two functions $f$ and $h$? $\endgroup$ – Matthias La Mar 5 at 10:56
  • $\begingroup$ Just looking at the first few lines of your code, it appears to be correct. But I wasn't getting to that, because you had errors in the math -- the one I mentioned, and you're calling out $h(y) = x$ when you should model $y$ as the dependent variable (and appear to in the code). I.e., $y = h(x)$. $\endgroup$ – TimWescott Mar 5 at 15:33
0
$\begingroup$

Hi Matthias La: I hate to be critical but the first example at the link you provided is actually quite poor because they end up using exponential smoothing (their update for $\hat{x}$ can re-written as an exponential smoothing update) for the update of the state, $x_t$, without deriving it. The way to derive the exponential smoothing update for the state (They also come up with a heuristic formula for the exponential smoothing factor, $\alpha_{n}$, of $\frac{1}{n}$ which is not optimal. It actually ends up being a function of the SNR) is to write the whole system in the following way:

$y_t = x_t + \epsilon_t $ # Measurement Equation

$x_t = x_{t-1} + \gamma_t$ # State Equation

where $\epsilon_t$ and $\gamma_t$ are independent normal error terms.

With the framework above, if you derive the updating equations for the kalman filter, you will obtain the exponential smoothing update for the state ( and $y_t = \hat{x_t}$ as update for the measurement). In the statistics-econometrics literature, the KF setup above is called a "random walk + noise" model. I'm not claiming that the link is completely useless but, as far as an intro to KF, that is not a good first example because they pull things out of the sky and don't explain where they come from. That can be confusing when going through the later examples.

I hope the setup above helps to show how to obtain the ES update. Oh, I recommend below for a dense look at the KF. Harvey has a terse style which can be difficult at first (it was for me) but, in the end ( after a few reads:) ), its hits you that he actually provides a very nice explanation. I've read various parts of it more than a few times.

https://www.amazon.com/Forecasting-Structural-Time-Harvey/dp/0521405734/ref=sr_1_1?dchild=1&keywords=Andrew+Harvey+Kalman+Filter&qid=1615252088&sr=8-1

EDIT: I didn't mention it above but I remembered later that you'll need a method for estimating the measurement variance and the state variance if you want to run the updating equations. Harvey uses prediction error decomposition but my impression is that dsp people have a different way of estimating them. Prediction error decomp is the only method I'm familiar with but you may know of a better one. If you do and can explain it to me, it's appreciated.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.