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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
   卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。) {* C. c# R+ |$ Q" l9 Y& d2 v5 y
8 {% }% ~/ |2 O- ~4 w# [* J% I
%%%%%%%%%%%%%%%%%%' l2 c' T  k) C6 m) F8 S
%功能说明:Kalman滤波用在一维温度数据测量系统中! Z, k+ H' T" I8 v- Q- j: {! X$ I
function main
, d4 B( O# s) A4 Z/ w%%%%%%
  W) [( q( d! ]1 ?) rN = 120;                                    %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
3 Q8 E4 X8 T- b$ DCON = 25;                                %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动
& p$ `+ Q3 `( R2 F- g' V%ones(a,b) 产生a行b列值为1的矩阵
' ^" |2 A4 I; X- Y0 M  A- s0 U%zeros(a,b)产生a行b列值为0的矩阵5 F7 ~# E  M2 J0 i
Xexpect = CON*ones(1,N);     %期望的温度是恒定的25度,但实际温度不可能会这样的
1 Q' U* E+ q% h4 H1 L+ P
; `. E, u# y. z/ U% T1 n- gX = zeros(1,N);                         %房间各时刻真实温度值* K1 g8 J7 ^  ?+ b
Xkf = zeros(1,N);                      %Kalman滤波处理的状态,也叫估计值3 m# d, ~% J) b$ W; R5 z/ r3 N+ h) u
Z = zeros(1,N);                         %温度计测量值& h* }1 s+ `& `
P = zeros(1,N);
$ c3 v$ T) V5 Z  c' v$ `) M2 b5 w%X(1)为数组的第一个元素
! z8 @6 a- Y$ h8 ?  zX(1) = 25.1;     %假如初始值房间温度为25.1度
, _! P9 h0 T1 Z, n$ eP(1) = 0.01;     %初始值的协方差   (测量值 - 真实值)^2
9 s9 U/ T+ Y7 ^/ `
4 U7 G; C; U7 I1 K%产生0-1的随机数  模拟系统的随机数据) z1 ~+ k+ R0 Y% n: [
Z = 24+rand(1,N)*10 - 5;      
8 T& u: ?+ Q. U# W9 ~6 O
0 ]' c% V; v/ K% sZ(1) = 24.9;     %温度计测的第一个数据) C! ]0 i1 Y# b& n* w7 \7 v/ Z
Xkf(1) = Z(1);  %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声
7 H+ L, {# V: V
; E, b/ j. X5 \3 k9 M. M& r- i; UQ = 0;        %扰动误差方差(不存在扰动误差只有测量误差)8 X/ A$ Z# l# h( n( N, s
R = 0.25;        %测量误差方差
$ l) @+ c- @* {9 o" |%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号
, m. t) F* e0 M* r& E%W为过程噪声
( ?- p  U# c, [/ v, z%V为测量噪声
0 p4 ?: R" A- N7 @& vW = sqrt(Q)*randn(1,N); %方差决定噪声的大小
4 P1 p% z4 ~0 E3 B" `9 B4 VV = sqrt(R)*randn(1,N);%方差决定噪声的大小
5 v" P, B! A* A5 d- O1 {%系统矩阵
  v9 ?; B2 `$ b, \) e1 G%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解: j, L6 }0 [5 f$ J( f
F = 1;     
- w" l! \% P' C1 C2 wG = 1;  a  ^' Y* {' T/ h0 p2 W
H = 1;: d$ q& U3 x% R) S& g8 q: @0 K, @
%eye产生m×n的单位矩阵 数值应该为1
+ K! l: Y! w  b& wI = eye(1);  %本系统状态为一维
" \: U! B5 E+ @. t%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%. R8 k! E* \, N3 X/ G
%模拟房间温度和测量过程,并滤波
" c, h- B5 u2 t: _4 F
2 M1 p* {3 `, k' N9 `for k = 2:N
$ N  E: b4 P2 w$ P* Z    %第一步:随时间推移,房间真实温度波动变化/ X0 I. {% I& g7 J3 q( V
    %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程" F4 i- v) @0 h0 ]( V5 Z4 J
    X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声
: O9 D; {1 y3 [- m: N  k    %第二步:随时间推移,获取实时数据: k* n3 ]% g* B8 J8 E) u# I
    %温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,7 f1 n8 N$ T" T) B
    %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响5 \% c# W8 J  C. }- G
    %尽可能的逼近X(k),这也是Kalman滤波目的所在, s3 E! f& k* d0 B) x9 F

6 I. w( R: C- V3 R( t- |    %修改本次函数8 P2 g9 s* H  ?0 `0 j
    %Z(k) = H*X(k)+V(k);                %测量值为实际值叠加上测量噪声
& ?2 T7 d5 M! q: Y$ L; I7 [8 E# a' c2 \7 f+ m! S% T
    %第三步:Kalman滤波
* {- D- G9 @% \    %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,
) e0 Z* n! }+ ?& ~+ g    %读者可以对照(3.36)到式(3.40),理解滤波过程# E" _; ]! X( e4 o6 I! }- ?
    X_pre = F*Xkf(k - 1);                        %状态预测          X_pre为上一次卡尔曼滤波值9 I9 Q9 _1 z# Q1 ~7 r. T6 x
    P_pre = F*P(k - 1)*F + Q;                %协方差预测    % v2 n# o- l- b: k& ^
    %inv()为求一个方阵的逆矩阵。
) _2 }# v) K! A) D    Kg = P_pre*inv(H * P_pre*H' + R);  %计算卡尔曼增益: `( l7 O; D6 s" L8 P! x
    e = Z(k) - H*X_pre;                           %新息                 本次测量值与上次预测值之差$ e0 u- }$ D/ k
    Xkf(k) = X_pre + Kg*e;                     %状态更新         本次预测值
4 P7 N6 ?6 l+ A9 Q7 j. e/ i    P(k) = (I - Kg*H)*P_pre;                    %协方差更新/ m4 u0 Q7 `& I; U9 p% z( s
end
* ^' @: c( F+ h7 ^% I1 |, ?%计算误差
1 r$ C; d* W" o* yErr_Messure = zeros(1,N);   %测量值与真实值之间的偏差
% b8 m  l( G* R# l3 hErr_Kalman = zeros(1,N);     %Kalman估计与真实值之间的偏差
$ b* [) Z4 x7 [+ T, k2 o0 j+ dfor k = 1:N
! g7 |: G: U2 E    Err_Messure(k) = abs(Z(k) - X(k));    %abs()为求解绝对值
* K/ y9 m5 B! A    Err_Kalman(k) = abs(Xkf(k) - X(k));
+ s! J' Q2 n: L* _' e0 Send# D) \( ?4 C* r: o
t = 1:N;
1 q- U. F. J2 t7 W. tfigure   %画图显示
2 a) ?4 a& [6 z9 y9 ?+ g. F9 ]; N%依次输出理论值,叠加过程噪声(受波动影响)的真实值1 {* {+ g7 @* \8 y9 R4 O
%温度计测量值,Kalman估计值  ]  d& n+ w6 N# r
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');' ?: ^, t1 ~3 k! H, \, M
legend('期望值','真实值','观测值','Kalman滤波值');
: B& Y7 W! w+ l4 o$ E5 Gxlabel('采样时间/s');% f( c3 l6 z4 \7 `: o& F9 T1 s
ylabel('温度值/C');* f, a2 O, u0 r: @  p. S. J
%误差分析图
8 r! Q1 h, @5 U  V5 Z  afigure  %画图显示% J( L, g% H) Y9 |; a
plot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');" B! u) ^6 e3 [% j8 Q6 O) v
legend('测量偏差','Kalman滤波偏差');
( _, o5 X  Z. B) s4 jxlabel('采样时间/s');
9 Q3 V0 y- f2 w2 fylabel('温度偏差值/C');# G6 B" u2 d; {0 L; ?  z! I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
& }. y4 b2 U7 A) D5 T" O

& W. t* v6 a! k& s3 g1 u# P: U" c/ k  ~" `4 a5 N

% C  p4 f. Z( M. O6 |3 |
* D0 k6 p: u2 |& h) H& p
; X, k( G" U: f; l7 k
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-10-31 10:58 , Processed in 0.171875 second(s), 23 queries , Gzip On.

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

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

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