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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
   卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
" V1 w9 C! e- ?8 V6 n6 g
, C( x. ]# g- C' d4 w2 J0 n
%%%%%%%%%%%%%%%%%%% W2 O9 t0 h$ z  V
%功能说明:Kalman滤波用在一维温度数据测量系统中
+ q  `- \6 ?9 `/ afunction main: O$ l6 f+ w7 s& C% A) z
%%%%%%
. K: |' K+ Z' y' }5 {$ `% K- a3 }N = 120;                                    %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量9 |. ?4 i' S5 J% E7 O
CON = 25;                                %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动2 n" l: y( `5 I) }& @3 c2 m
%ones(a,b) 产生a行b列值为1的矩阵
  U. S, ^: ^, v$ Q. a%zeros(a,b)产生a行b列值为0的矩阵) D4 D% y+ }$ }
Xexpect = CON*ones(1,N);     %期望的温度是恒定的25度,但实际温度不可能会这样的
$ Q* K/ ^7 k2 H+ R8 d" L0 z, Z* F8 Y( U4 e% P5 y
X = zeros(1,N);                         %房间各时刻真实温度值+ i: k( ]# z% K# ^$ o% w5 v
Xkf = zeros(1,N);                      %Kalman滤波处理的状态,也叫估计值
$ n, W+ \$ t' M5 q0 l( O4 Z# pZ = zeros(1,N);                         %温度计测量值5 N  C/ t" E& J0 @% _3 U
P = zeros(1,N);# w* q" }" h2 n) G
%X(1)为数组的第一个元素
, `1 I, u& d+ sX(1) = 25.1;     %假如初始值房间温度为25.1度
1 @. P' p: Q/ V/ H  RP(1) = 0.01;     %初始值的协方差   (测量值 - 真实值)^2) Y* x/ C# F  V2 ~, J

: w; ?0 M' n( V& D& }%产生0-1的随机数  模拟系统的随机数据4 f7 h- q7 O* ]2 D
Z = 24+rand(1,N)*10 - 5;      
  c9 V& {; A0 ]: R! s, I% `" [! c# H/ T/ _: i* E$ E, w1 M
Z(1) = 24.9;     %温度计测的第一个数据
$ u) y% N1 k7 Z: DXkf(1) = Z(1);  %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声( O: p" x; k- V3 }

( r4 u- V( b% O5 tQ = 0;        %扰动误差方差(不存在扰动误差只有测量误差)) H; c% T; B, x/ h! Q
R = 0.25;        %测量误差方差
, f! j+ s! l, M* K0 ~%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号
4 p* C( p" p- M%W为过程噪声
4 F6 A# B' A" k% i: t; R%V为测量噪声6 b2 I4 A5 t4 L0 G
W = sqrt(Q)*randn(1,N); %方差决定噪声的大小9 ?& C1 p4 J2 T: L- g
V = sqrt(R)*randn(1,N);%方差决定噪声的大小0 O7 y3 H5 h' L
%系统矩阵( @5 y: d: [$ h! j! D
%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解
5 S. o% M2 t( J/ XF = 1;     
0 F/ [) J1 w  ?  N  `+ X. hG = 1;
* b9 O% D. q! y, K# @* i  n% Z; n1 nH = 1;" I  t: D8 r% n9 W' x+ q7 u
%eye产生m×n的单位矩阵 数值应该为1
2 m" g4 O7 W3 T. Y2 gI = eye(1);  %本系统状态为一维+ j. |$ }4 Y. n( V5 L! P
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
' t2 m2 M, z' U- ~7 d+ f%模拟房间温度和测量过程,并滤波
& g6 g  o0 \' W- ~2 M- O
# |: t% \# v/ d/ Ffor k = 2:N
# @4 T. t4 g4 j0 R    %第一步:随时间推移,房间真实温度波动变化+ G, Q+ i( `/ m6 P, l5 ~
    %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程
( S3 X8 n) M" l- \) G+ G    X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声1 N  @9 F% L. Z- c
    %第二步:随时间推移,获取实时数据
6 h) j7 w4 {, g) i, U; i# ]    %温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,
3 ^. N, e/ Y; O7 j6 Z& H' w    %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响
9 r8 C$ a2 z8 H0 n( _    %尽可能的逼近X(k),这也是Kalman滤波目的所在8 h% x: f! |4 D

* m( R1 @$ Z/ t% G( Q    %修改本次函数
+ L) R/ S! \5 }1 X. x    %Z(k) = H*X(k)+V(k);                %测量值为实际值叠加上测量噪声* X: @4 @! S. @
) [( Z; U1 j/ K0 d6 B' K6 g
    %第三步:Kalman滤波
  i$ A9 U& h! m! _8 F7 E2 A, s    %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,
6 y8 c1 ~$ r8 M( o9 M2 s1 e    %读者可以对照(3.36)到式(3.40),理解滤波过程8 H9 j, v  u/ U2 D; F) V0 j0 ]
    X_pre = F*Xkf(k - 1);                        %状态预测          X_pre为上一次卡尔曼滤波值
/ x' [) z1 h2 v, r( Y, A- v% J    P_pre = F*P(k - 1)*F + Q;                %协方差预测   
# V+ b. T( o/ `/ w( S    %inv()为求一个方阵的逆矩阵。
7 A: C- X1 u( }8 P: J7 U( C1 }    Kg = P_pre*inv(H * P_pre*H' + R);  %计算卡尔曼增益. x' N+ u) T; C% J/ L8 t0 {( K7 `$ c
    e = Z(k) - H*X_pre;                           %新息                 本次测量值与上次预测值之差! U0 U$ W' c4 Y; r( i: q
    Xkf(k) = X_pre + Kg*e;                     %状态更新         本次预测值7 e' L0 @$ ^: z7 I$ ~# y& J; u
    P(k) = (I - Kg*H)*P_pre;                    %协方差更新
% I( d+ q1 ~4 }; R$ Nend
8 T1 i7 P! p/ W) v; h. C%计算误差
% l# J, N5 R; r$ f: w4 c( k0 I* Y  cErr_Messure = zeros(1,N);   %测量值与真实值之间的偏差
( _# d# G' t) a9 s, g) Q: G# \2 _Err_Kalman = zeros(1,N);     %Kalman估计与真实值之间的偏差
9 t! o+ J2 X* S+ |, ~8 s: ^for k = 1:N 7 V+ s6 z$ `: Y
    Err_Messure(k) = abs(Z(k) - X(k));    %abs()为求解绝对值( Z' t- A* |; \! B8 A; J+ a
    Err_Kalman(k) = abs(Xkf(k) - X(k));4 u% e' o( e: i; H; \; c* E
end
3 O: g0 o4 M$ L* q% o. \% M! p7 ~t = 1:N;
: a" {9 Z7 u- w( ]% o  }$ t2 D2 s/ Afigure   %画图显示
, D8 ~9 r' x+ Q; f) Y%依次输出理论值,叠加过程噪声(受波动影响)的真实值
" S" k" E1 \' c" C2 g- Q2 @1 ^4 c. U%温度计测量值,Kalman估计值- G: x9 @, s, {& n) ~( r8 o
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');0 q4 m0 C  o; ?# K. ~. V
legend('期望值','真实值','观测值','Kalman滤波值');/ I6 n1 _  ?3 L: E# K) Z4 K, y
xlabel('采样时间/s');  J* l( C# Z% M) ]
ylabel('温度值/C');
3 Q# b1 N! K( `! S, R( ]%误差分析图& U* Q  e6 b+ l; c9 P" t; F
figure  %画图显示
0 U" S8 V$ Y3 B/ oplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
0 C9 U8 t7 d9 w9 J" N5 J5 a# Slegend('测量偏差','Kalman滤波偏差');
1 o& r# p9 b9 G, `0 Bxlabel('采样时间/s');
6 d" \0 |& m: K2 kylabel('温度偏差值/C');
2 |: i6 g5 o" M4 |: s/ L  N# T%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
& K7 o0 o# g' T

" A% m* n  d, Z
. u- @# H  a$ k% F$ N; }% r

* l7 [6 l4 r+ ]# d/ Q( P. P- o1 l
  \) R2 ~# C" c; F3 }
( y( h. Z$ `( G8 i
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-11-23 19:53 , Processed in 0.156250 second(s), 24 queries , Gzip On.

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

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

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