|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
; ? R3 a% c; O9 d9 e. S! h3 J; i2 p$ s+ A: j+ g
%%%%%%%%%%%%%%%%%%
# v" a6 ?' d1 s/ |2 Q& q( f, w%功能说明:Kalman滤波用在一维温度数据测量系统中! e$ x. y' d5 {- Z" X! |) L
function main
! U3 [' H& _+ `- T( h%%%%%%
9 S- n6 t, b3 o$ j# s+ Y" hN = 120; %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量" W3 `* ^2 \. Z# H. F& L
CON = 25; %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动
, [* Z' c' t: f# J" u: d: z%ones(a,b) 产生a行b列值为1的矩阵. m% t0 [9 ^& Y. ~2 x- P
%zeros(a,b)产生a行b列值为0的矩阵
. a$ a2 M3 ^1 \) }% }( `' UXexpect = CON*ones(1,N); %期望的温度是恒定的25度,但实际温度不可能会这样的& @ m m+ G9 s: o2 g
* {) e3 ^9 J6 {5 j8 w; YX = zeros(1,N); %房间各时刻真实温度值
, T, c3 _ F" L5 l& W1 mXkf = zeros(1,N); %Kalman滤波处理的状态,也叫估计值
+ H- V- Q6 d. k4 N- r$ E, c. f; rZ = zeros(1,N); %温度计测量值$ ?% r# Q+ A0 }& b% a
P = zeros(1,N);
$ a3 H9 k6 W. F: @! h%X(1)为数组的第一个元素
8 K* y6 B% L. m( eX(1) = 25.1; %假如初始值房间温度为25.1度3 [* L# W5 c9 r }
P(1) = 0.01; %初始值的协方差 (测量值 - 真实值)^2* C. I. t" K j. t, Q2 }! D( ]4 W( Q
7 `! E `& B1 G$ E
%产生0-1的随机数 模拟系统的随机数据* z: K$ r ?' G D; L2 z: W
Z = 24+rand(1,N)*10 - 5; " G4 h/ y/ K q
3 @' J3 n ^1 S+ ]Z(1) = 24.9; %温度计测的第一个数据
+ z1 s1 C* ^- w: y2 o7 H4 v, m) ^% L% ZXkf(1) = Z(1); %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声6 q2 u6 b$ h. `( a1 B5 j, f2 m4 X9 k
5 L, H' I4 C3 A4 Q3 t1 }# zQ = 0; %扰动误差方差(不存在扰动误差只有测量误差)
- s8 L8 l6 h4 D6 DR = 0.25; %测量误差方差$ n8 h* H4 ? J% O
%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号. Y- N7 g* ^$ S3 a! f0 |
%W为过程噪声( S3 m- }+ T& b% @, Q2 D: ?1 w
%V为测量噪声% y5 }. z# d/ ~! T4 \ h$ H
W = sqrt(Q)*randn(1,N); %方差决定噪声的大小8 P% p% W Q4 k `9 B; _% X2 ?
V = sqrt(R)*randn(1,N);%方差决定噪声的大小
0 ^$ t/ t7 ^4 a" k9 a Q6 {. V%系统矩阵
! p% q0 [6 G7 F% w% A* G%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解, T" u8 k6 v6 a# O" N
F = 1; ! ]+ @ V5 s6 p$ A6 l" G! F
G = 1;* u- y. P( B* H8 k0 I2 G( q0 s3 a. q
H = 1;: i, k0 c% Q, o) u4 i* r
%eye产生m×n的单位矩阵 数值应该为1
/ L9 L% n5 V7 J. m. ]3 OI = eye(1); %本系统状态为一维
1 g. z+ o2 d) Z; Y1 g, }* y%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%6 f/ ]- s$ F7 L) {4 {/ o4 l( K
%模拟房间温度和测量过程,并滤波
" e( M, Z, Q6 o% ]: l8 j! v
% n7 p: j& @! T, B" kfor k = 2:N7 I; o: a& ^: W
%第一步:随时间推移,房间真实温度波动变化0 ?- P1 J1 r" } ~: c: {- }
%k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程+ p" m4 V4 s; V0 _( ?5 T+ D
X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声
" u; }, ]# j9 ?- r" H9 l9 W+ X %第二步:随时间推移,获取实时数据; j( f8 `/ I6 h& U( Z
%温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,' \$ @' A9 [- p
%它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响; U" a9 _* q) H( X. a- P
%尽可能的逼近X(k),这也是Kalman滤波目的所在) @" r" o" z2 Y4 ~- [+ U
2 X* P/ v9 c8 b/ h/ _ %修改本次函数
, H0 z- ~4 ] m7 g %Z(k) = H*X(k)+V(k); %测量值为实际值叠加上测量噪声
2 o/ [, i" y* g5 L1 z! A: i
0 a& |: _1 [! M %第三步:Kalman滤波
7 _1 F9 g! Y# I( f %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,# I# l4 e6 p. t2 T: j+ T9 G
%读者可以对照(3.36)到式(3.40),理解滤波过程% h7 M1 T3 h0 l; I
X_pre = F*Xkf(k - 1); %状态预测 X_pre为上一次卡尔曼滤波值# Z% K4 y1 Y/ B0 y( Q5 r
P_pre = F*P(k - 1)*F + Q; %协方差预测
, F. V$ b( X7 k6 p0 ?* R; l2 _ %inv()为求一个方阵的逆矩阵。! l, Y8 ^: k4 r
Kg = P_pre*inv(H * P_pre*H' + R); %计算卡尔曼增益3 Q) R! }' U1 U2 H& g3 v
e = Z(k) - H*X_pre; %新息 本次测量值与上次预测值之差: e4 r7 f2 V0 e/ F# |, U
Xkf(k) = X_pre + Kg*e; %状态更新 本次预测值
; c, W) K. k2 x4 o. b, r* C P(k) = (I - Kg*H)*P_pre; %协方差更新8 k$ M' T4 N3 D; z% v0 P' _% z) _, ?
end) @& I. Q5 Q6 O4 i- s
%计算误差
# k+ Y$ L- ~+ @9 y4 s' ]0 YErr_Messure = zeros(1,N); %测量值与真实值之间的偏差
& a) F+ }4 p* @. EErr_Kalman = zeros(1,N); %Kalman估计与真实值之间的偏差0 C/ w g2 O/ k: R
for k = 1:N
r% u- |: p: e' T& v4 I9 { Err_Messure(k) = abs(Z(k) - X(k)); %abs()为求解绝对值
0 ?$ s1 u+ w& `4 P Err_Kalman(k) = abs(Xkf(k) - X(k));
! W: s, }3 s* G; n9 i+ g8 qend
( \6 c; }, g! o8 G6 at = 1:N;
; d! ^ R0 s6 X* lfigure %画图显示
2 |% b. v9 f4 N%依次输出理论值,叠加过程噪声(受波动影响)的真实值
1 f* p/ Y( M5 s, a8 u%温度计测量值,Kalman估计值6 M( ]: P2 A' W; d _( t5 B* e* a
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');% [( h5 o7 K$ Z
legend('期望值','真实值','观测值','Kalman滤波值');
% _3 A$ _( g; C" dxlabel('采样时间/s');
4 b% _: x/ v5 l* K) Q$ Bylabel('温度值/C');
8 a# e1 X9 f; g- Y8 ]% T%误差分析图
* m! J& F' |! V6 T+ ]figure %画图显示6 M3 A2 o% U, C, M3 d4 {
plot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
" H; n) `. q$ v0 ulegend('测量偏差','Kalman滤波偏差');# k. H. w# @. }. U- g" A* [
xlabel('采样时间/s');/ s6 G3 d! m4 n2 ]9 X
ylabel('温度偏差值/C');
9 i+ k/ x3 s$ o8 ~) m- \# O) T5 G( t%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/ m6 S9 z- q& C: w
/ E& t6 @6 M. v- g& ^3 ?2 Y2 s: x
, K3 _: ? p! Y& O' c3 \' K8 i3 v1 U- J8 v
1 {, K1 v1 ~2 B _- q5 T
9 B1 s+ c, } T2 q! O0 ^5 v |
|