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

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

[复制链接]

该用户从未签到

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

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 ^
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

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

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

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

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