|
|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。9 L$ \6 B, H c' C- j
& U- I) {' Z, T/ |7 t9 a%%%%%%%%%%%%%%%%%%% ^2 f. A5 h( U- x2 G
%功能说明:Kalman滤波用在一维温度数据测量系统中( I* N5 J2 a3 m' F, V' E) g
function main2 k$ b, Y) h) i# D( M) }
%%%%%%8 y! G F1 p. r5 M* B0 b
N = 120; %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
7 o* b5 h2 P- \$ X v/ K |& iCON = 25; %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动* m& n% n2 U0 D) s9 D3 O% r$ ?
%ones(a,b) 产生a行b列值为1的矩阵8 b2 M; G- e1 t- }
%zeros(a,b)产生a行b列值为0的矩阵
9 _: ~$ T7 B! D. `6 AXexpect = CON*ones(1,N); %期望的温度是恒定的25度,但实际温度不可能会这样的' {2 n$ ?2 a; r( C5 T- n
4 T7 X. ?+ Q! h/ |! N5 WX = zeros(1,N); %房间各时刻真实温度值
% z; V5 M1 O: v- P, l( Q' w9 f. nXkf = zeros(1,N); %Kalman滤波处理的状态,也叫估计值
D' x F- x2 Y8 I0 n1 E" oZ = zeros(1,N); %温度计测量值
- G% [2 H; @$ N. h; V) vP = zeros(1,N);( e- C# c4 T) V$ C
%X(1)为数组的第一个元素
$ y2 y1 Y& X+ V7 S$ g3 s" x2 dX(1) = 25.1; %假如初始值房间温度为25.1度1 K$ z: i' z! d+ p d) ?
P(1) = 0.01; %初始值的协方差 (测量值 - 真实值)^2
7 v A" G3 J, R" n
J- l* W# s2 X2 k/ n%产生0-1的随机数 模拟系统的随机数据
. `3 T) u& @1 Y: {' DZ = 24+rand(1,N)*10 - 5; # O( S! b7 O1 k
$ K/ f; P5 M1 F% g3 N8 w8 @Z(1) = 24.9; %温度计测的第一个数据
' F8 t$ O. _5 v7 J7 {* sXkf(1) = Z(1); %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声& |# U4 ?; ^" ? W: Q& Z* _, |
0 t3 `* L u* t, w* b$ M+ F: w
Q = 0; %扰动误差方差(不存在扰动误差只有测量误差)8 O5 o1 J4 c) Y' ?
R = 0.25; %测量误差方差' T1 e j& {5 X6 Y; F8 r' D6 `' g
%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号
% ~/ n+ a! ]5 l3 ~, ]%W为过程噪声# p- m J, F" u) u# q. f2 x
%V为测量噪声
6 K; P& F; X( X+ H8 a9 M+ c4 ]W = sqrt(Q)*randn(1,N); %方差决定噪声的大小$ k+ C: B. X% l, q2 E7 [ V" n
V = sqrt(R)*randn(1,N);%方差决定噪声的大小
* g2 }! a7 \& T. c9 @# B%系统矩阵
; z' h9 o: k6 O6 k$ L [%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解6 s. z6 g7 X/ W9 r
F = 1; / m1 z, Z5 e6 N# a0 i; T9 f$ C
G = 1;- a6 b/ _- s" k" F
H = 1;
$ j- w1 Y% e1 o. n6 h1 f M%eye产生m×n的单位矩阵 数值应该为1
+ M: x; e# ^; s1 ?. [I = eye(1); %本系统状态为一维
) d; s% v: E0 |7 M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0 [( k* k2 J& L" Y7 C% |! ~%模拟房间温度和测量过程,并滤波
0 X2 z6 c( e1 ^; x# C0 O* Y9 _9 e( t1 T+ j/ \& Z% B- F
for k = 2:N/ c3 K" ~4 Y# x3 b. r
%第一步:随时间推移,房间真实温度波动变化
% ?/ X8 s, _% F8 Y% c; @# O %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程9 ^, m$ ]) D2 I1 B) d. Y3 b4 ]+ w
X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声- K. V! b: c5 D; e
%第二步:随时间推移,获取实时数据4 l6 C- w3 z/ w# m+ \
%温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,/ e8 W/ s ]3 J& z6 y2 T m
%它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响
! l4 P% y! c2 h %尽可能的逼近X(k),这也是Kalman滤波目的所在 a5 K: Z* l1 x& m( p5 k
8 K i" B! t+ j* C5 |
%修改本次函数
; d3 N+ m% h7 C8 c3 z6 e %Z(k) = H*X(k)+V(k); %测量值为实际值叠加上测量噪声$ E* I# J( e( {+ ~% `, F& O+ B
3 w8 u6 c. ?" i: {! i
%第三步:Kalman滤波6 d1 _6 h8 {# T8 x# J) x
%有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,
, g2 O; C0 r2 \! y2 H7 O %读者可以对照(3.36)到式(3.40),理解滤波过程. C: E3 F* B; R* s% o" J! R5 K
X_pre = F*Xkf(k - 1); %状态预测 X_pre为上一次卡尔曼滤波值6 q0 b) H8 z3 d4 Z
P_pre = F*P(k - 1)*F + Q; %协方差预测
: k) s) L5 E' U- @ %inv()为求一个方阵的逆矩阵。' k( }. Q2 V/ n7 y
Kg = P_pre*inv(H * P_pre*H' + R); %计算卡尔曼增益
, u4 T) U( z) @ e = Z(k) - H*X_pre; %新息 本次测量值与上次预测值之差
+ ]. X8 R* M+ h b( r8 q) w0 x- R Xkf(k) = X_pre + Kg*e; %状态更新 本次预测值* G' a0 E& J1 B
P(k) = (I - Kg*H)*P_pre; %协方差更新
! u7 s9 f) V7 p$ W) Y8 e# Gend5 V' i( `0 k3 g9 F& d% |
%计算误差
5 e! r# @6 L# W& Z9 `Err_Messure = zeros(1,N); %测量值与真实值之间的偏差8 D9 t& ~( l+ F1 S- B' p
Err_Kalman = zeros(1,N); %Kalman估计与真实值之间的偏差
2 T# n+ A7 y) @9 ?% U1 r- ]for k = 1:N 4 D% N" e: e& ^# S
Err_Messure(k) = abs(Z(k) - X(k)); %abs()为求解绝对值1 ]3 p. [* r) B* d/ o! H
Err_Kalman(k) = abs(Xkf(k) - X(k));
( e k! ], q1 d" _: oend& `1 z) f l2 r( w2 V Q7 ?/ T" ?
t = 1:N;
& S$ u4 e4 ]- R5 {8 u0 efigure %画图显示 b6 }, ]. o- C' |5 \9 I5 E
%依次输出理论值,叠加过程噪声(受波动影响)的真实值. O, p+ W P. P
%温度计测量值,Kalman估计值
" Q6 b/ x- Q0 @. L2 T, a2 mplot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');6 i' P. g: ?( V+ _! Z
legend('期望值','真实值','观测值','Kalman滤波值');) z1 v7 K* q3 B
xlabel('采样时间/s');
$ ~2 w+ z5 i8 W* o0 m& yylabel('温度值/C');
8 L& _" N( e+ i7 G%误差分析图+ K. |! F9 @0 o& t( Z4 o' Y
figure %画图显示
" E$ {% f( T6 M; r# l% Mplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
5 @5 j$ A' v9 n2 |8 }6 J, B5 H& ylegend('测量偏差','Kalman滤波偏差');
; [6 n' t G$ Uxlabel('采样时间/s');2 s- N$ E3 J+ T2 T/ y' r; g& \; N. n
ylabel('温度偏差值/C');1 G" H: Z/ }! z. k
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%: H1 }* k6 }" a1 D
* ?! \4 B5 d# Z, }, B! V( E2 s z+ O& j3 n' [# Y: u
. H( a7 y4 K+ ]8 j; `4 }# f, M
- m- ?7 R- g( G6 R, Y- O
9 V9 o2 U' x. |
|
|