EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
最近因为要用到Kalman滤波器做东西,所以一直在学习这个东西,鉴于之前的仿真都是用matlab做的,所以呢,kalman滤波器的仿真程序也是用matlab编的。痛苦了几天,几天这个函数终于搞定了,具体的分析如下。function [zx,zy]=xyKalmanFliter(A,H,Rw,Rv,Rw_c,Rv_c,x0,p0,y)
8 x& ]) v# e1 W5 Q6 m7 j1 J$ P% T%----------------输入参数--------------------------------------------------) L$ N: A k& r& t; M: X
% A -- 系统矩阵. i1 N' R# h2 ^
% H -- 观察矩阵. l* l8 Y' k3 w9 L) F4 U1 k B
% Rw -- 扰动向量
+ ?6 Q8 J8 v. C7 [% Rv -- 测量噪声
8 X: y# e& x/ `0 N* V; Y0 r( i% Rw_c -- 扰动向量的协方差
+ ^# y6 ^# j& S+ r% Rv_c -- 测量噪声的协方差# q A& K! r9 h3 @, ]
% x0 -- 节点初始位置向量(x,y)'! f* {% ?1 c, P4 o2 b5 Z7 ^9 A
% p0 -- 初试协方差阵' ?! }" P* B% L% a+ {: Q |6 m
% y -- 采样周期
# O8 t' f1 T$ p1 q# X, n, l: {' P) d%-------------------------------------------------------------------------- & _8 X; f4 [, o; g( P
%--------------------------------------------------------------------------# o* n/ }& u" M7 a
% X(k) = A*X(k) + Rw(k) 噪声Q- S) t' |7 K) O3 C1 R
% Z(k) = H*X(k) + Rv(k) 噪声R
& A# `/ H/ V+ V7 o; A % x(k|k-1) = A*x(k-1|k-1) + B*U(k)
: S0 A# G! Z2 ]% l% P(k|k-1) = A*P(k-1|k-1)*A' + Q: L9 k7 ~& V4 D9 V5 L* U+ ~
% x(k|k) = x(k|k-1) + kg(k)*(z(k)-H*s(k|k-1))
+ A7 X8 a% H* w$ H% kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R)
. o% [9 Y! a; m' a% P(k|k) = (I-kg(k)*H)*P(k|k-1)
* C* W3 K* D$ \. `2 n%-------------------------------------------------------------------------- len = length( y ); %获取采样点数 r = size(A,1); %获取系统矩阵A的行数* D* @, R* w" _1 ^3 _2 B3 j
I= eye( r(1) ); %生成单位矩阵' B& \- }# E% w7 \/ I: A8 v
P1 = zeros( r(1),r(1) ); %初始化协方差阵 %初始化节点位置,协方差阵
/ R+ C3 l2 w3 q; ^7 l. X& k( w X = x0;6 l; ~" \0 p; G y9 L: ^. s* \
P = p0; len1 = size( H*X ,1);
: _- T$ d' [+ Z" v, r! D/ K. |9 V for s=1:1:len+ j# c* q& c/ |
z1 = A*X+ Rw; % X(k) = A*X(k) + Rw(k) 协方差 Rw_c {0 c" S; T9 X, O. R: t
zx(1:r(1),s) = z1(1:r(1)) ; . y3 b2 _ c( w
z2 = H*X + Rv; % Z(k) = H*X(k) + Rv(k) 协方差 Rv_c
0 L" z5 ]' h2 p" g2 c zy(1:len1,s) = z2(1:len1 );
! a4 U5 s7 t8 x & U- B* O n7 ~ z0 s! u! @3 F
P1 = A*P*A' + Rw_c; % P(k|k-1) = A*P(k-1|k-1)*A' + Q
& {( S% }8 Z: e; p7 T# [( J K = P1*H'*inv( H*P1*H' + Rv_c ); % kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R): R K4 U6 @4 S9 d* R
X = A*X + K*( zy(1:len1,s) - H*A*X ); % x(k|k) = x(k|k-1) + kg(k)*(z(k)-H*s(k|k-1))
8 s3 d8 A0 j* |8 g. b P = ( I - K*H ) * P1; % P(k|k) = (I-kg(k)*H)*P(k|k-1)9 e0 U! z, N, G2 @
end 8 F0 s S" o2 R0 q2 p2 q
return 原文地址:转:Kalman滤波的Matlab仿真程序解读作者:浩瀚
6 I( u0 g8 `- W" R) |% ` |