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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
   卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。
* b$ @7 b: h0 x
% W3 s- f  r% l( m2 F9 D
%%%%%%%%%%%%%%%%%%* b) H: ^# P5 u, s7 E- F: G- g
%功能说明:Kalman滤波用在一维温度数据测量系统中2 J) r- Q7 m, r+ X4 l
function main
1 C- n" e7 [8 p: V8 R; b, Q%%%%%%
; l' T: @7 K/ T0 H0 d  ?N = 120;                                    %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
' x5 _, q' S& ~4 ~0 m& ^CON = 25;                                %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动
5 `8 a% y0 B& z9 [$ s# U( R%ones(a,b) 产生a行b列值为1的矩阵
: W8 f! W3 H8 \. v/ d; o! ~; p%zeros(a,b)产生a行b列值为0的矩阵( J3 l) |8 C$ M
Xexpect = CON*ones(1,N);     %期望的温度是恒定的25度,但实际温度不可能会这样的) k  H$ s' k5 i1 r( d3 m" Z% l* D

+ V$ Q  _: N+ L# x# I/ _X = zeros(1,N);                         %房间各时刻真实温度值; H/ L: m4 W2 I( d! o7 W
Xkf = zeros(1,N);                      %Kalman滤波处理的状态,也叫估计值
; s3 o* r- h2 |Z = zeros(1,N);                         %温度计测量值% d  @. H/ I" A6 P2 n
P = zeros(1,N);/ ]6 J& C$ c$ E! o
%X(1)为数组的第一个元素
% _  s& |7 H7 Q) \- Q6 fX(1) = 25.1;     %假如初始值房间温度为25.1度: Y& R7 L9 u" F- L- `$ k: W
P(1) = 0.01;     %初始值的协方差   (测量值 - 真实值)^2
' O, Y3 Z) p" b2 _" T, a" g
$ Q# ^+ r5 W% T%产生0-1的随机数  模拟系统的随机数据/ `- c( R/ T, M* j: N/ y# G
Z = 24+rand(1,N)*10 - 5;      
6 l' J. t6 e* ?) k' T3 V) a% r5 @% [& W, o4 s; p8 T
Z(1) = 24.9;     %温度计测的第一个数据
/ C/ z, k. e. _. z, e0 `Xkf(1) = Z(1);  %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声
4 k& k6 @4 O0 q9 v8 ]9 D4 }/ w& u) r8 F$ k( H2 j) b3 P
Q = 0;        %扰动误差方差(不存在扰动误差只有测量误差)* L; t7 A9 r: m2 b% L) C- d
R = 0.25;        %测量误差方差5 b: L. x, R/ w! ~$ F
%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号, [+ j: N8 `0 }3 T1 P
%W为过程噪声& ]4 K; h8 z7 P& l3 F. U
%V为测量噪声
/ {* J/ e0 c  X2 F. H' y5 u# ^8 vW = sqrt(Q)*randn(1,N); %方差决定噪声的大小9 q" z/ ]7 z0 y5 C7 q, t
V = sqrt(R)*randn(1,N);%方差决定噪声的大小8 u6 s, ^& u, D& g. I
%系统矩阵
( f! z! w% q4 N, Z; E6 Y%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解! p& M1 w% O) C1 L% p
F = 1;     
+ l! D1 T8 f0 J0 C+ gG = 1;; i3 I6 u* i1 P$ j6 J
H = 1;1 s4 K9 Z4 w9 {
%eye产生m×n的单位矩阵 数值应该为12 S% k: l6 Z5 B/ O6 C0 S0 j% o
I = eye(1);  %本系统状态为一维8 _( d& c2 ]& J" R* W4 U
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' H% e6 k( ]4 \5 K; ]
%模拟房间温度和测量过程,并滤波
# `$ `* P/ i/ k1 c6 K
3 U! ?, t5 C. B; [& Z" @for k = 2:N- N; v, L  d7 M1 n8 k/ J
    %第一步:随时间推移,房间真实温度波动变化1 s& |# j. h# j7 k  A' X
    %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程
/ C, `/ e1 \: r" o    X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声) Z" T1 f- s6 S( R2 d& k7 O
    %第二步:随时间推移,获取实时数据* y& l% j- q9 K$ \
    %温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,8 Y! n, ?& {0 H1 T
    %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响
3 [( n; ]) m8 Y    %尽可能的逼近X(k),这也是Kalman滤波目的所在8 c/ r, }! [: u0 P! y5 w' n/ q) j
  z! u/ y3 F" L4 _4 i  ]6 Y8 @+ @
    %修改本次函数6 X% D7 P# N( v
    %Z(k) = H*X(k)+V(k);                %测量值为实际值叠加上测量噪声1 w! V; P6 {/ V6 n$ Q7 ]7 K
  `7 n/ N5 W$ e# u  l7 v
    %第三步:Kalman滤波& |7 w2 E/ g  i2 H7 h
    %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,
' V3 ?: n% v0 H# J  J    %读者可以对照(3.36)到式(3.40),理解滤波过程5 d! y* P+ Z) X
    X_pre = F*Xkf(k - 1);                        %状态预测          X_pre为上一次卡尔曼滤波值2 s# v! O  r8 M# L. |) r8 e% O
    P_pre = F*P(k - 1)*F + Q;                %协方差预测   
4 v& Q6 a+ g% m& u6 s    %inv()为求一个方阵的逆矩阵。+ K: E% {$ |) S) d& M0 I( y& C
    Kg = P_pre*inv(H * P_pre*H' + R);  %计算卡尔曼增益
* ]9 [4 ?( S8 _- g2 z$ f" A    e = Z(k) - H*X_pre;                           %新息                 本次测量值与上次预测值之差
5 O. |. o4 g' g# L    Xkf(k) = X_pre + Kg*e;                     %状态更新         本次预测值
. n; R& O0 `3 F- N8 W: E# x    P(k) = (I - Kg*H)*P_pre;                    %协方差更新1 _  P2 Y3 H# L% d! s5 Z
end
  G: g! E/ F* ]! X6 O4 ~2 S/ S: q8 `%计算误差
2 P9 U" \* L6 C( X" W- lErr_Messure = zeros(1,N);   %测量值与真实值之间的偏差
( R+ m& }) u, c8 TErr_Kalman = zeros(1,N);     %Kalman估计与真实值之间的偏差# f/ `1 i! I; Q  @+ C2 c' l2 ]+ h
for k = 1:N
3 j1 a( [0 f0 ~" @) L    Err_Messure(k) = abs(Z(k) - X(k));    %abs()为求解绝对值, R2 }, R: }* o7 L# j# S' H$ m
    Err_Kalman(k) = abs(Xkf(k) - X(k));
$ D  S* E  `. ^) R3 q: nend0 Q( M! p2 t% p$ w. W, H8 N% i
t = 1:N;
4 b' x( o3 M9 \; a9 C- Gfigure   %画图显示8 A- b: d' p  s3 Q2 o/ P
%依次输出理论值,叠加过程噪声(受波动影响)的真实值
7 v( y$ t: B2 @8 I& c5 H9 O%温度计测量值,Kalman估计值% O8 H" y9 F8 t8 P& ?3 o; F
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');9 i- q  e! w+ Y
legend('期望值','真实值','观测值','Kalman滤波值');
$ e- v/ h3 _6 c6 [/ G% ^xlabel('采样时间/s');
" M8 G' D+ C6 g+ P+ C5 M3 q2 q/ eylabel('温度值/C');
4 f: X, e, r4 K+ S& s%误差分析图, I0 Q+ K  P* P3 h# ?& O% I! D
figure  %画图显示4 n- _. R0 t* f, g
plot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');$ q; N& `9 I$ Y3 U9 z, K! T
legend('测量偏差','Kalman滤波偏差');  D4 v2 k0 r- ?
xlabel('采样时间/s');* U7 P6 c* q# u
ylabel('温度偏差值/C');
7 {) D# N" f4 f( a%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0 E: ~. y4 ]" U& X8 `2 k; i
& h# S' r: \: ?. B" P1 ]5 `5 ~! ~& q, p# X5 I6 d
0 m# c  s" o; j, u
1 m3 Y4 e/ ~# w: q) i( [4 \* j

4 S7 v) _6 b& N# W
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-8-4 14:39 , Processed in 0.109375 second(s), 23 queries , Gzip On.

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

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

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