|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
, Q( Y. @- I0 }7 n! \0 t& y
b9 {' r: [& B4 ?- ~/ `%%%%%%%%%%%%%%%%%%
9 a5 b. M. l: v/ b%功能说明:Kalman滤波用在一维温度数据测量系统中
# a: i' n: l [8 n- B% z3 tfunction main
5 V) C- ~# Q( ~$ L; _0 X J%%%%%%
* Z2 O" u9 _7 V# C' W, JN = 120; %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
6 R2 f t; l7 i+ [# jCON = 25; %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动) b" g& u" e z6 ]3 `" X% f; t
%ones(a,b) 产生a行b列值为1的矩阵
% P8 B0 E1 b/ c7 L% n%zeros(a,b)产生a行b列值为0的矩阵
1 c9 A4 O. m8 i1 y4 e# `" F. pXexpect = CON*ones(1,N); %期望的温度是恒定的25度,但实际温度不可能会这样的
' m2 C4 {% {. P( S4 G# H4 d' N F1 U3 ~ B- D5 @
X = zeros(1,N); %房间各时刻真实温度值+ W4 M5 X& B7 f8 y* c1 j
Xkf = zeros(1,N); %Kalman滤波处理的状态,也叫估计值5 x8 @5 `0 m ~4 d# R w3 Y
Z = zeros(1,N); %温度计测量值& l6 I, G- g/ `' N( T0 ]+ j: w$ g
P = zeros(1,N);# H* v. \5 `! N8 @7 ~4 B
%X(1)为数组的第一个元素, @7 ^9 B# S' ~& m2 \, u& J, O0 }
X(1) = 25.1; %假如初始值房间温度为25.1度5 ]! K6 T& U+ o6 w$ s- T
P(1) = 0.01; %初始值的协方差 (测量值 - 真实值)^28 }3 g* O4 a2 d3 u* l
* Y( G6 p8 H: ~
%产生0-1的随机数 模拟系统的随机数据% X3 R5 ~* E4 k+ o5 e4 C. Y
Z = 24+rand(1,N)*10 - 5;
* y' ~3 a9 Z& L4 }% s' O/ ?0 A
/ e4 h% z5 V" K. O7 V! E& E( ~- WZ(1) = 24.9; %温度计测的第一个数据# k7 Y* ]9 n9 L# f
Xkf(1) = Z(1); %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声
' u% \2 D& ^8 r5 N- L
; i6 d1 [6 A4 E( _! Y& sQ = 0; %扰动误差方差(不存在扰动误差只有测量误差) p# c( }0 X8 U5 X: [
R = 0.25; %测量误差方差
% c0 _" E* |+ I) b%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号
+ @- M3 L8 E( @9 l%W为过程噪声
1 i# k: ^ L0 y%V为测量噪声
) c5 s3 N' L7 }6 N( r; H# P% k9 NW = sqrt(Q)*randn(1,N); %方差决定噪声的大小
+ t! [/ B0 c4 e `V = sqrt(R)*randn(1,N);%方差决定噪声的大小( P5 a" ?: l, x0 q8 \/ C( ^
%系统矩阵
7 p5 @1 f- ^, g7 u( Z0 S%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解* F+ H: S$ r% \' A. D! H
F = 1;
0 `; z6 F$ s( Q B4 t, P. u* qG = 1;9 B! o+ s0 q6 R7 T) F/ [' O
H = 1;
* s$ ]' k# P) M% s/ `2 F%eye产生m×n的单位矩阵 数值应该为1" \+ v" d1 E& P
I = eye(1); %本系统状态为一维: W& K% m! ^8 o2 }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' {8 f- `7 `, K# z6 Y8 u
%模拟房间温度和测量过程,并滤波
: }# Z1 C5 A; U. Z. F
# E' C+ Q9 @) ~6 v8 q. p4 E0 Vfor k = 2:N0 G& c4 |0 l3 x4 l5 w# v
%第一步:随时间推移,房间真实温度波动变化
& d% s/ ?4 S; e %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程5 ?6 g1 l5 H9 E" E
X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声6 j# J# v. n1 d, _
%第二步:随时间推移,获取实时数据 G7 _' T6 k! N5 X! x9 S0 ]1 A- u, X
%温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,; X8 [! F, C3 y+ c) g
%它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响& U5 X/ x8 O- m; r6 ]
%尽可能的逼近X(k),这也是Kalman滤波目的所在$ T- F; `8 S/ |: l! l j5 Z# d
- D8 s6 z1 E$ Z8 r0 |2 h %修改本次函数% u7 a- T7 F' N: O, C" }7 B
%Z(k) = H*X(k)+V(k); %测量值为实际值叠加上测量噪声
/ P- F t( q @( c `1 b7 K. h2 m; f+ C0 q2 H
%第三步:Kalman滤波
! a$ L6 ^6 f2 H' } %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,, }1 e! }- |' l3 I l
%读者可以对照(3.36)到式(3.40),理解滤波过程
V1 b+ K' S. o9 Y X_pre = F*Xkf(k - 1); %状态预测 X_pre为上一次卡尔曼滤波值) {9 K V; x8 h3 q5 ^' W
P_pre = F*P(k - 1)*F + Q; %协方差预测
+ r* @ k- X# O! U %inv()为求一个方阵的逆矩阵。
- p6 o+ }' O: c, e6 B Kg = P_pre*inv(H * P_pre*H' + R); %计算卡尔曼增益/ ~* Z1 u1 e! y/ s
e = Z(k) - H*X_pre; %新息 本次测量值与上次预测值之差
4 I# m3 u& P, q3 i$ q8 K( z3 ~; E Xkf(k) = X_pre + Kg*e; %状态更新 本次预测值# ~0 }' ]9 _; f) d
P(k) = (I - Kg*H)*P_pre; %协方差更新2 t+ q1 k9 s( w4 A+ _1 n1 [3 G
end
; T2 \3 a$ o3 [, A9 t; `%计算误差
! Y9 w* g0 s/ a* X; D5 wErr_Messure = zeros(1,N); %测量值与真实值之间的偏差
8 m6 ^7 W4 X5 xErr_Kalman = zeros(1,N); %Kalman估计与真实值之间的偏差+ @! {; j! }& U7 M3 O4 _
for k = 1:N 0 C* M. N/ B' V+ ]% g' _, ?0 M' T0 F
Err_Messure(k) = abs(Z(k) - X(k)); %abs()为求解绝对值5 X" p6 m* w2 I" Q6 |1 U9 v5 {
Err_Kalman(k) = abs(Xkf(k) - X(k));
& P6 K6 s8 u6 `% |% |7 `end
; Q) p, \8 d, t9 A8 R! Y& E' W9 {t = 1:N;: L4 O1 ?, R- ]+ A1 G9 z
figure %画图显示2 R) k" v2 A7 a% ^ V0 B2 G
%依次输出理论值,叠加过程噪声(受波动影响)的真实值, g4 a- T/ U5 ?5 n6 k5 m
%温度计测量值,Kalman估计值
$ T; x" M: ]* G! A% aplot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');
; J* ]$ B$ H0 {& B+ klegend('期望值','真实值','观测值','Kalman滤波值');
0 [+ a* U8 U% m* o; B2 Uxlabel('采样时间/s');
! G$ {8 h& f! S Xylabel('温度值/C');) i* I; R8 b0 V! | p0 p" d& z) w
%误差分析图
C K" {% ]/ _& X0 p% Afigure %画图显示
( R$ x1 k$ ~& Vplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');( l% G: Y* l# ?( ]* B/ L+ W
legend('测量偏差','Kalman滤波偏差');
! ]. g( p. ?; S* y% hxlabel('采样时间/s');
0 P& t( N0 G* s4 ^; l- Qylabel('温度偏差值/C');. `* _: p8 x& O/ Y6 d# N
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0 T b" C0 W8 Q# i5 m+ F$ h% X
. L+ w; U6 L: |& w
4 w& u. i' _+ `, C9 S7 a, h3 f+ E; w6 z5 Y- P0 q/ F: s$ g
1 _6 {: L3 [5 P, a% [
3 E: z' P9 x" ^7 ^
|
|