|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
# r1 z2 W. k" T( X1 I7 E0 _) @7 ?3 _* e4 a, j8 ~9 M7 W
%%%%%%%%%%%%%%%%%%
1 P* O& S' k% p1 L" a%功能说明:Kalman滤波用在一维温度数据测量系统中
2 {( Y. W/ @. b" jfunction main
( k& z# X" _- l& O5 \: Y* U+ \%%%%%%4 [0 l! o2 Q8 ^
N = 120; %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
5 a6 z& k. B0 W4 F' W, Q2 U9 Q0 U% yCON = 25; %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动. r5 y4 a% u1 f% A/ p8 Q
%ones(a,b) 产生a行b列值为1的矩阵# k3 \2 v; x% Z- N1 Q
%zeros(a,b)产生a行b列值为0的矩阵
; @! f* X# R% {* VXexpect = CON*ones(1,N); %期望的温度是恒定的25度,但实际温度不可能会这样的$ U! Z% C: Y- \: z
3 L( _; @, q/ x$ _: T0 x
X = zeros(1,N); %房间各时刻真实温度值
3 _3 l8 z# H$ X% t3 C0 S( m; g7 Z5 KXkf = zeros(1,N); %Kalman滤波处理的状态,也叫估计值1 _# u- u3 y. Z- u- b3 A0 q, Z
Z = zeros(1,N); %温度计测量值5 g2 s% Q4 _( e+ F, w
P = zeros(1,N);# |3 p. i ? i- k, \& y/ C
%X(1)为数组的第一个元素
6 Y7 S( l" e! M# Y2 c5 kX(1) = 25.1; %假如初始值房间温度为25.1度
2 W2 a4 C7 g. u& [4 G; G4 VP(1) = 0.01; %初始值的协方差 (测量值 - 真实值)^2
8 C# E) Z) O5 T$ q" L+ a, p. L5 S) j
%产生0-1的随机数 模拟系统的随机数据
3 r/ @3 W9 {- L. w( TZ = 24+rand(1,N)*10 - 5;
. k, H% Z9 o4 O" G
. x2 L4 S) k# Y# WZ(1) = 24.9; %温度计测的第一个数据
5 s7 Z5 E: @+ MXkf(1) = Z(1); %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声+ x0 A+ U R) K* o- i% A8 F! F
1 a3 [( \) a4 Z7 E7 e
Q = 0; %扰动误差方差(不存在扰动误差只有测量误差)# C2 I6 W# m8 v
R = 0.25; %测量误差方差; ]2 j7 m; `/ H
%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号+ T) D$ `! N) h% I/ W3 J& h, r
%W为过程噪声" q1 F& E' g) H; ^
%V为测量噪声
+ i4 E% f, ?2 p/ }$ L0 Q8 Z" CW = sqrt(Q)*randn(1,N); %方差决定噪声的大小
: Z& V1 N2 M; x4 rV = sqrt(R)*randn(1,N);%方差决定噪声的大小0 q6 P" m- d1 K5 o- Y
%系统矩阵+ g& l5 o3 d0 v" c% c1 U' |
%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解
3 ^, A. z& w8 sF = 1; 3 `, q) ~: L! |
G = 1;
; m7 o$ O# o* w6 C# n& S6 NH = 1;( J% l# w& y& R. m
%eye产生m×n的单位矩阵 数值应该为1
( G$ ~* v" m* n5 m; t+ DI = eye(1); %本系统状态为一维7 ]8 R% l. J% |; }6 v
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%) D: e% |; S Q! c
%模拟房间温度和测量过程,并滤波) l2 U/ o! r1 K8 e
# C' [- u4 g& S
for k = 2:N
6 o, G; k" O# j %第一步:随时间推移,房间真实温度波动变化
" \+ o+ k% a' e: X1 c %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程% q0 U+ H: _9 K/ T6 _' I' m9 B
X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声
0 T9 Q6 U6 r( ^, x" J; D %第二步:随时间推移,获取实时数据1 j4 Z; X. i. n1 h
%温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,% V8 b* s0 h* _0 Q( I4 a( r
%它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响; ~ J' q5 P. ?8 k5 n
%尽可能的逼近X(k),这也是Kalman滤波目的所在
" X% K/ y' ?! e \5 A ]* h' ?
1 s; G ~: g6 j; p! x %修改本次函数
& i+ |2 l, C- ` |# n6 Y7 x %Z(k) = H*X(k)+V(k); %测量值为实际值叠加上测量噪声
5 v# }" |" m. j/ r" O7 ?& D6 w$ c; K% d. G+ B
%第三步:Kalman滤波" y! N$ x. P; p
%有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,) k& G# m/ T( X) z+ \
%读者可以对照(3.36)到式(3.40),理解滤波过程
! y4 I& h n, C4 A- r2 s6 Z0 E' q X_pre = F*Xkf(k - 1); %状态预测 X_pre为上一次卡尔曼滤波值& T5 h, m( |. Z4 s& z3 G, Z2 [
P_pre = F*P(k - 1)*F + Q; %协方差预测 ) [# v* V7 E4 y. k2 B
%inv()为求一个方阵的逆矩阵。
, |( j; {- L5 n. r) C$ x5 e Kg = P_pre*inv(H * P_pre*H' + R); %计算卡尔曼增益) }' T7 O" `+ S* B/ F& G
e = Z(k) - H*X_pre; %新息 本次测量值与上次预测值之差
! ]) s! m- @/ d. q+ E: s1 l/ e Xkf(k) = X_pre + Kg*e; %状态更新 本次预测值
1 v8 w' j9 p) @* h9 D9 D P(k) = (I - Kg*H)*P_pre; %协方差更新" K, u7 \# `7 W
end( B1 j* n& V) a) @: d# L* ^
%计算误差
$ X: L2 v; g3 Z% ^( v0 r2 oErr_Messure = zeros(1,N); %测量值与真实值之间的偏差
]' z( Q9 Z- NErr_Kalman = zeros(1,N); %Kalman估计与真实值之间的偏差
" P. \$ h( o0 `% P s/ G, r( z9 Rfor k = 1:N ' T% K- `4 E* M) T) a
Err_Messure(k) = abs(Z(k) - X(k)); %abs()为求解绝对值
# f7 l& y$ ^$ |. _6 _ Err_Kalman(k) = abs(Xkf(k) - X(k));
( b$ H# _$ O; |$ Y# ~end% \, v9 E1 M7 R0 y% I _
t = 1:N;
3 a( `' J# [9 g8 A! K. qfigure %画图显示% E% W: o. M, u) o
%依次输出理论值,叠加过程噪声(受波动影响)的真实值
) e9 \" C3 |/ \( J: ~# L D& v%温度计测量值,Kalman估计值+ O3 m; Y! W- z0 x+ J
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');5 T) a5 H8 M9 U
legend('期望值','真实值','观测值','Kalman滤波值');9 x) m, }" k- R% ^
xlabel('采样时间/s');
. m4 V+ W5 ^, i' \. d( [3 D" zylabel('温度值/C');! O4 W, t Q& I1 J0 U, d
%误差分析图
+ M# q5 K. `; u# k# B: B$ K1 [figure %画图显示
- [1 N/ k( w6 g2 j7 Qplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
; W3 W! v: Q3 elegend('测量偏差','Kalman滤波偏差');
Q3 ?4 D" y3 j% U# _* C/ j7 p6 n `xlabel('采样时间/s');
3 V' ]0 ^' c& |7 c1 C# Q/ fylabel('温度偏差值/C');0 T) s5 B% w0 n& k: V: b
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%, q5 ?$ X Z7 l6 ]" o
2 q& B5 m0 i0 _, _) c7 M3 u0 G) ]0 _% l7 z4 g8 I% y/ |' I, @# u
- P& U5 ~# `% v& ~: |7 q# \: a
4 V4 I: m! x" P, N, B
7 h) J" w9 t7 f1 C8 N: M2 m. o |
|