kalman_update.m
上传用户:mozhenmi
上传日期:2008-02-18
资源大小:13k
文件大小:1k
源码类别:

其他小程序

开发平台:

Matlab

  1. function [xnew, Vnew, loglik, VVnew] = kalman_update(F, H, Q, R, y, x, V, initial)
  2. % KALMAN_UPDATE Do a one step update of the Kalman filter
  3. % [xnew, Vnew, loglik] = kalman_update(F, H, Q, R, y, x, V, initial)
  4. %
  5. % Given
  6. %  x(:) =   E[ X | Y(1:t-1) ] and
  7. %  V(:,:) = Var[ X(t-1) | Y(1:t-1) ],
  8. % compute 
  9. %  xnew(:) =   E[ X | Y(1:t-1) ] and
  10. %  Vnew(:,:) = Var[ X(t) | Y(1:t) ],
  11. %  VVnew(:,:) = Cov[ X(t), X(t-1) | Y(1:t) ],
  12. % using
  13. %  y(:)   - the observation at time t
  14. %  A(:,:) - the system matrix
  15. %  C(:,:) - the observation matrix
  16. %  Q(:,:) - the system covariance
  17. %  R(:,:) - the observation covariance
  18. %
  19. % If initial=true, x and V are taken as the initial conditions (and F and Q are ignored).
  20. % If there is no observation vector, set K = zeros(ss).
  21. if nargin < 8, initial = 0; end
  22. if initial
  23.   xpred = x;
  24.   Vpred = V;
  25. else
  26.   xpred = F*x;
  27.   Vpred = F*V*F' + Q;
  28. end
  29. e = y - H*xpred; % error (innovation)
  30. n = length(e);
  31. ss = length(F);
  32. S = H*Vpred*H' + R;
  33. Sinv = inv(S);
  34. ss = length(V);
  35. loglik = gaussian_prob(e, zeros(1,length(e)), S, 1);
  36. K = Vpred*H'*Sinv; % Kalman gain matrix
  37. xnew = xpred + K*e;
  38. Vnew = (eye(ss) - K*H)*Vpred;
  39. VVnew = (eye(ss) - K*H)*F*V;