找回密码
 注册
关于网站域名变更的通知
查看: 1760|回复: 3
打印 上一主题 下一主题

利用matlab实现卡尔曼滤波算法

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2019-8-13 09:49 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

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' c
3 \' K8 i3 v1 U- J8 v
1 {, K1 v1 ~2 B  _- q5 T

9 B1 s+ c, }  T2 q! O0 ^5 v
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-7-20 08:51 , Processed in 0.125000 second(s), 23 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表