EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
最近因为要用到Kalman滤波器做东西,所以一直在学习这个东西,鉴于之前的仿真都是用matlab做的,所以呢,kalman滤波器的仿真程序也是用matlab编的。痛苦了几天,几天这个函数终于搞定了,具体的分析如下。function [zx,zy]=xyKalmanFliter(A,H,Rw,Rv,Rw_c,Rv_c,x0,p0,y)
: z) }5 |1 C4 H# U8 o/ w%----------------输入参数--------------------------------------------------
+ _ y" Q/ [! r% A -- 系统矩阵
/ l$ z! X# j& R h0 O& P; ~) B% ~% H -- 观察矩阵
: c+ l' h: U! R4 @6 L% Rw -- 扰动向量
, I+ C3 a9 E5 ~; ]: ?' A6 S- ]% Rv -- 测量噪声4 V" C# w5 ^0 v9 W0 C6 a! @
% Rw_c -- 扰动向量的协方差% ]# b( \' f* H# P" p, {: t
% Rv_c -- 测量噪声的协方差" C/ j! R5 r# M5 C! c
% x0 -- 节点初始位置向量(x,y)'- i4 W; A/ M; r3 t' b; \+ B
% p0 -- 初试协方差阵7 c3 D6 r( s' j9 M8 p
% y -- 采样周期/ P6 v1 n5 a* U% V8 N
%-------------------------------------------------------------------------- \/ G2 B3 Y! D; c% X
%--------------------------------------------------------------------------$ E; R/ ^4 o" `( X2 t
% X(k) = A*X(k) + Rw(k) 噪声Q/ k% j, g g: |) E" d
% Z(k) = H*X(k) + Rv(k) 噪声R/ T, s% ^$ ^2 e
% x(k|k-1) = A*x(k-1|k-1) + B*U(k)! G) O- b5 f% r6 V( W; S8 A
% P(k|k-1) = A*P(k-1|k-1)*A' + Q8 A5 G1 f U C0 ?) Q- h! C
% x(k|k) = x(k|k-1) + kg(k)*(z(k)-H*s(k|k-1))
0 N) g( Z# O& \! @% F2 G' B% kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R)
$ Z8 w: y& x/ Q) ?5 ]% P(k|k) = (I-kg(k)*H)*P(k|k-1)
1 r$ p" O" ~; S6 D2 Q%-------------------------------------------------------------------------- len = length( y ); %获取采样点数 r = size(A,1); %获取系统矩阵A的行数2 X7 M1 P* i% m
I= eye( r(1) ); %生成单位矩阵# `7 N2 ~ e2 }2 _& H
P1 = zeros( r(1),r(1) ); %初始化协方差阵 %初始化节点位置,协方差阵
) p: T% |: v- O: a X = x0;2 T& d, ^6 z& d7 c' N
P = p0; len1 = size( H*X ,1);( o, k( E7 E6 C5 J( c! Z3 p! x3 n. b
for s=1:1:len9 G4 t& V4 s6 t) k s8 B
z1 = A*X+ Rw; % X(k) = A*X(k) + Rw(k) 协方差 Rw_c
5 B! }- Y$ m& o zx(1:r(1),s) = z1(1:r(1)) ;
1 {% c& q3 u- O- p( p" H8 @ z2 = H*X + Rv; % Z(k) = H*X(k) + Rv(k) 协方差 Rv_c0 M. q$ Q: ]3 }' r! j. D5 S$ s
zy(1:len1,s) = z2(1:len1 );& L6 n8 J1 h! L, i& B: l
4 D. ^. b4 L0 V) y7 F P1 = A*P*A' + Rw_c; % P(k|k-1) = A*P(k-1|k-1)*A' + Q v8 h; \$ N, E( K! k7 q
K = P1*H'*inv( H*P1*H' + Rv_c ); % kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R)
; f: r0 [; U$ _ B+ L 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 a) }$ P1 Y. w( a
P = ( I - K*H ) * P1; % P(k|k) = (I-kg(k)*H)*P(k|k-1)
( ?* R, b' {- n; ?5 R end
: D/ D3 G A: h/ B0 F0 Qreturn
原文地址:转:Kalman滤波的Matlab仿真程序解读作者:浩瀚 ) B' ^2 @5 c* t
|