EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
最近因为要用到Kalman滤波器做东西,所以一直在学习这个东西,鉴于之前的仿真都是用matlab做的,所以呢,kalman滤波器的仿真程序也是用matlab编的。痛苦了几天,几天这个函数终于搞定了,具体的分析如下。function [zx,zy]=xyKalmanFliter(A,H,Rw,Rv,Rw_c,Rv_c,x0,p0,y)
) w5 s9 {$ U) `+ z1 I%----------------输入参数--------------------------------------------------
4 B9 x! o1 v% g2 r# ~! f% A -- 系统矩阵
1 h( w4 W% e" l' l% H -- 观察矩阵
% S) b/ u( b1 S9 y: g0 U% Rw -- 扰动向量
& u% k8 T$ u5 l6 l% Rv -- 测量噪声( E! B- z& X K7 w' O
% Rw_c -- 扰动向量的协方差8 n* [! e. e% j" z
% Rv_c -- 测量噪声的协方差
0 K" \$ C8 m" `( `/ ?- p0 r% x0 -- 节点初始位置向量(x,y)'+ _* g5 u* l2 z8 F1 U7 n
% p0 -- 初试协方差阵5 Y# l ^9 }9 J% G- P1 {) c" |
% y -- 采样周期0 t" f& y& M1 Z# r! k. W
%--------------------------------------------------------------------------
+ ~+ j% A0 S! [, l0 |3 X* J%--------------------------------------------------------------------------. u5 w+ C, h. l" d/ ~
% X(k) = A*X(k) + Rw(k) 噪声Q
" s: [8 z* A. Q! \8 g% Z(k) = H*X(k) + Rv(k) 噪声R N; ~, v) a& Q* L& C: {/ s a5 A
% x(k|k-1) = A*x(k-1|k-1) + B*U(k)( @& H+ [9 M+ Q( x K
% P(k|k-1) = A*P(k-1|k-1)*A' + Q
4 Y6 N! Y" f( F- J3 W! w% x(k|k) = x(k|k-1) + kg(k)*(z(k)-H*s(k|k-1))
3 F' z7 Z8 r/ S* d% kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R)
5 P: {' b. Z1 P7 w; Q; g$ F2 @% P(k|k) = (I-kg(k)*H)*P(k|k-1)
8 w' i: n; y* }% E) J& Z6 f%-------------------------------------------------------------------------- len = length( y ); %获取采样点数 r = size(A,1); %获取系统矩阵A的行数& l6 m$ O# d G) D* K2 c8 E) z o+ T
I= eye( r(1) ); %生成单位矩阵* D4 u4 F: L+ b9 g
P1 = zeros( r(1),r(1) ); %初始化协方差阵 %初始化节点位置,协方差阵3 L7 ?5 ^- z+ F' S' [, s" W
X = x0;' x6 }7 ]6 ^/ r: k( b5 `# }* m/ a
P = p0; len1 = size( H*X ,1);$ T5 Q# n3 u4 G6 S# P3 b
for s=1:1:len/ ^5 Q# l& T* P2 Q
z1 = A*X+ Rw; % X(k) = A*X(k) + Rw(k) 协方差 Rw_c
- e1 @1 v& R( K" } zx(1:r(1),s) = z1(1:r(1)) ; 0 q O$ X; \% E: R. W
z2 = H*X + Rv; % Z(k) = H*X(k) + Rv(k) 协方差 Rv_c% i0 O7 R! ~% f1 L4 h
zy(1:len1,s) = z2(1:len1 );
! Z! U$ u5 D. o! W6 X0 @4 e ) ]5 o' N) g7 M5 p8 f% U: @
P1 = A*P*A' + Rw_c; % P(k|k-1) = A*P(k-1|k-1)*A' + Q
5 H% @3 R+ C* I' R8 `: N' ]0 q K = P1*H'*inv( H*P1*H' + Rv_c ); % kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R), l$ ~9 `, X% r _8 b! T4 Y
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))3 ~, |5 D. G8 W; I- u* x+ I
P = ( I - K*H ) * P1; % P(k|k) = (I-kg(k)*H)*P(k|k-1)2 s/ B0 |3 l+ V6 V
end & J; ^6 X4 [0 H- A, ^* |
return 原文地址:转:Kalman滤波的Matlab仿真程序解读作者:浩瀚
4 b0 o5 X+ i7 a7 `2 \6 d; \ |