|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
% `3 D8 a/ M( e
- o0 w1 V' L- a! }. C1 z%%%%%%%%%%%%%%%%%%
/ G( t# R+ W& ^# v( ~2 @( n%功能说明:Kalman滤波用在一维温度数据测量系统中
& E4 P# f& Z. ?4 O! |% |function main7 `: S; i( W/ o7 }4 A( N# w: p
%%%%%%
( b' ]9 u5 [, B# g6 o* ]N = 120; %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
6 L# Q9 ]7 q6 C) iCON = 25; %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动! l7 Y8 e" K. n4 z$ N& H: e
%ones(a,b) 产生a行b列值为1的矩阵
1 p+ z* s6 g- f' r; f" l# ^3 m%zeros(a,b)产生a行b列值为0的矩阵7 r# B1 O! S: g2 L
Xexpect = CON*ones(1,N); %期望的温度是恒定的25度,但实际温度不可能会这样的
t3 ?- |" _' Z; ~$ s
) X, y, f( Z" O8 h" c3 L" iX = zeros(1,N); %房间各时刻真实温度值( E' { `, v( p3 M* w) `
Xkf = zeros(1,N); %Kalman滤波处理的状态,也叫估计值. W y7 d$ @" {- I
Z = zeros(1,N); %温度计测量值( ^7 e2 v) Y% R j5 Y
P = zeros(1,N);' P9 }! i0 h( E/ d& x) y! g
%X(1)为数组的第一个元素
q& A% T. r. ~6 TX(1) = 25.1; %假如初始值房间温度为25.1度1 K8 @6 r7 }2 @( j, @
P(1) = 0.01; %初始值的协方差 (测量值 - 真实值)^2+ }0 [% \: C! A4 o
; s. Y% N8 R9 Q/ J%产生0-1的随机数 模拟系统的随机数据
. H f( w8 @8 V2 kZ = 24+rand(1,N)*10 - 5; ; J4 j1 l3 f8 `! P7 G+ h2 m
1 U% c2 k3 }" O6 w+ `1 `
Z(1) = 24.9; %温度计测的第一个数据6 a0 ]. K" B- Y7 u5 J
Xkf(1) = Z(1); %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声: d( t' d5 `7 |: @0 A1 R
3 d$ q/ f- E' M
Q = 0; %扰动误差方差(不存在扰动误差只有测量误差)" |1 w; w& \6 _2 L5 ~; O
R = 0.25; %测量误差方差
8 M1 \0 B) R5 a3 Q$ |%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号# S9 Q9 E; o$ c D% f3 c
%W为过程噪声' \& ~' [& a. R R/ u
%V为测量噪声
/ Q: \2 l4 ]: cW = sqrt(Q)*randn(1,N); %方差决定噪声的大小
8 q( V% k3 W. B( L, ZV = sqrt(R)*randn(1,N);%方差决定噪声的大小
' K) j9 A" b4 ~! Z* S%系统矩阵
! G {. x1 w; T% J6 r1 `( |+ M%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解
8 z" s8 ]% d0 LF = 1; 5 e6 ]9 n; w* U3 t2 ]5 `& l
G = 1;
- J2 V7 A7 `$ ~1 sH = 1;
* Q" k; f1 g8 c%eye产生m×n的单位矩阵 数值应该为18 F* [1 n# W) g" i9 Y2 C! \
I = eye(1); %本系统状态为一维
4 y; V* ^/ a" s%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2 J- E: m# w- Q, Y. {/ r5 X0 ~
%模拟房间温度和测量过程,并滤波
3 T" ^# v, u# V# `, P0 r) d! V7 ^$ C9 K0 ]! W$ r
for k = 2:N6 O+ K. s1 U: Y; m& q7 J4 Y
%第一步:随时间推移,房间真实温度波动变化: B, `) B) r+ a: g7 y
%k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程
; _9 x0 K) K& ~ X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声
% ]# R/ T& Y6 d' u$ ]7 n$ z( H! z %第二步:随时间推移,获取实时数据( U- f* V# A6 ^& P. Z
%温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,
2 _9 n3 C" W0 b6 b6 F, a8 K7 B %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响 C, D7 X8 o9 t- z
%尽可能的逼近X(k),这也是Kalman滤波目的所在* x3 M( q& J/ l
. E3 q* I6 N& \" S2 a1 \* y$ I %修改本次函数2 m& [0 a4 U8 _+ E
%Z(k) = H*X(k)+V(k); %测量值为实际值叠加上测量噪声$ J' `/ d/ S& Y; _" q1 _9 d6 M1 L
! F, L& y! m8 p9 F5 [) Q %第三步:Kalman滤波2 S1 `2 S/ x/ ^9 T% m4 c* [
%有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,& R! h( e, [! q- ]) t- @
%读者可以对照(3.36)到式(3.40),理解滤波过程7 _, s* m% _2 H; w4 S' ] g0 }& j
X_pre = F*Xkf(k - 1); %状态预测 X_pre为上一次卡尔曼滤波值
' I% f- @) b9 y" M4 k, y$ x7 b P_pre = F*P(k - 1)*F + Q; %协方差预测
* k! c2 u/ U1 W& U& f$ j %inv()为求一个方阵的逆矩阵。; u% L/ t" Q+ Y
Kg = P_pre*inv(H * P_pre*H' + R); %计算卡尔曼增益
6 B0 S6 o+ o+ h' r* K; f e = Z(k) - H*X_pre; %新息 本次测量值与上次预测值之差& d/ h9 q( f9 F3 {+ U/ k, R
Xkf(k) = X_pre + Kg*e; %状态更新 本次预测值
8 t" D, R, ?( G4 w1 A2 @) Z/ x3 ] P(k) = (I - Kg*H)*P_pre; %协方差更新6 w; I9 V$ H, S! a9 y! P
end7 p j- \. w) Z8 Y. F% O6 `
%计算误差4 ^! a2 r) O/ J. x% o I& l& t
Err_Messure = zeros(1,N); %测量值与真实值之间的偏差
# N! G& G! T" w! r$ g+ [; D4 {) TErr_Kalman = zeros(1,N); %Kalman估计与真实值之间的偏差
5 u& K# M) t6 K$ q$ ~for k = 1:N 8 u, n! k6 K8 M7 }! }! p5 I
Err_Messure(k) = abs(Z(k) - X(k)); %abs()为求解绝对值
9 d8 x, v7 b! E6 ` Err_Kalman(k) = abs(Xkf(k) - X(k));
7 @2 V" L% C# e* E( k# }% Wend
! Q: t$ J$ A5 O# D/ F' ~9 n3 et = 1:N;
( n# A& ]" p5 u( Y U" Z9 G( mfigure %画图显示
0 K3 f) h) E9 ^2 V* R7 V& y' [%依次输出理论值,叠加过程噪声(受波动影响)的真实值
) b* \- d8 ?6 P U" G%温度计测量值,Kalman估计值; C" c* }8 O9 ]7 I
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');
# K- ?+ P& u$ X; X0 `# D% t7 Elegend('期望值','真实值','观测值','Kalman滤波值');6 M$ z1 p, @ {
xlabel('采样时间/s');
" i( S9 s0 c( b. Sylabel('温度值/C');7 |: {9 N4 ?& ]$ ?/ F
%误差分析图
4 U7 |' s$ _" r& N4 P$ hfigure %画图显示
& L& `) J- c) Y& Yplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');6 O4 T0 U' C0 |6 @, y6 [5 c9 o# d3 X* t
legend('测量偏差','Kalman滤波偏差');5 h! I' m3 }2 ?: j
xlabel('采样时间/s');
1 P" J* K! I3 G5 k+ Z6 `9 jylabel('温度偏差值/C');+ B3 i t& [" Y) @9 [& z4 C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%. Y/ ^6 k( i6 A8 V
3 r& s5 \# O5 i1 m4 A
. n" W: F6 z, ~0 E. O
( |* Y' l; o0 t. W# Y4 t: n: U9 N7 }( h+ F* _
4 r& ^3 M; W8 V& f' \ |
|