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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
   卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。2 \# u, v, \/ b0 s3 E
1 ?$ A+ X2 J/ z! y  e* @# S* I$ G
%%%%%%%%%%%%%%%%%%1 B: r6 h: C+ i
%功能说明:Kalman滤波用在一维温度数据测量系统中
, \6 L  w7 E' Z$ p+ d+ J% Ffunction main7 n) E% f4 K7 i! y9 J
%%%%%%
: F4 K7 V1 r" iN = 120;                                    %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
' n2 q+ g7 x' c; G3 Y- cCON = 25;                                %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动
8 \, \  k5 l$ a  w%ones(a,b) 产生a行b列值为1的矩阵
- V7 Q- w4 b' X. Z. L0 C%zeros(a,b)产生a行b列值为0的矩阵
, E) P4 M: N8 m; EXexpect = CON*ones(1,N);     %期望的温度是恒定的25度,但实际温度不可能会这样的. ]" z- x( l; s7 v& ]& j5 T( a9 W8 I
( i% W% ?* f2 S0 r
X = zeros(1,N);                         %房间各时刻真实温度值2 u! P& F# p+ I* A9 E) j
Xkf = zeros(1,N);                      %Kalman滤波处理的状态,也叫估计值* S. ^# w2 z" J
Z = zeros(1,N);                         %温度计测量值
; O  E$ q* m6 _4 g- J2 gP = zeros(1,N);
& Z5 O; g! h; F+ ~1 `%X(1)为数组的第一个元素3 L3 t- z# h# g9 c! R8 x! g! W
X(1) = 25.1;     %假如初始值房间温度为25.1度8 f4 u* ?/ ]$ P, u9 I
P(1) = 0.01;     %初始值的协方差   (测量值 - 真实值)^2' n+ h# p( w- a. X9 x6 Z
2 z" F  q0 [; y6 |3 w
%产生0-1的随机数  模拟系统的随机数据! r0 C2 P- @8 U/ f; {0 |9 u5 W  w
Z = 24+rand(1,N)*10 - 5;      
* T" Q' i% a2 J4 j% g0 |7 H% W: v# k$ x3 J0 A& o
Z(1) = 24.9;     %温度计测的第一个数据( ?8 i+ V; H+ g
Xkf(1) = Z(1);  %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声
  s0 O$ ^/ S, Y. k4 i- {$ B' H1 s
Q = 0;        %扰动误差方差(不存在扰动误差只有测量误差)
) g( V$ I: R* H1 ^5 b$ F3 ZR = 0.25;        %测量误差方差
5 S/ x3 P. I) V/ R& {%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号0 a" y. V2 P& x4 v; K5 j
%W为过程噪声3 Q- {+ c# Q$ M. A. R( C
%V为测量噪声/ y+ r3 H3 |* l) Z3 i  b
W = sqrt(Q)*randn(1,N); %方差决定噪声的大小
1 W( f( p; ~! u; LV = sqrt(R)*randn(1,N);%方差决定噪声的大小0 R% o# Y% K3 X4 n
%系统矩阵' F7 k, Q2 _2 t+ X, Y2 J
%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解
" }2 _6 t3 n5 HF = 1;     ) J8 [! N: k- u& G9 }
G = 1;
. M0 p4 a* v3 i) y- |2 N7 wH = 1;
4 u2 l, ]1 M% ^" N; q! Z%eye产生m×n的单位矩阵 数值应该为1) l# N+ s* j3 Z# l* J1 Y. y/ t. _4 B
I = eye(1);  %本系统状态为一维+ V' j! E. Q" n6 W- [
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%) {4 n4 n; r  Q* B0 U
%模拟房间温度和测量过程,并滤波3 N  J1 U) U1 |6 n8 _7 X
5 k4 }4 ]& }9 f6 c5 x& S" w4 W
for k = 2:N
& J3 e0 ?, j! F* z" V( O    %第一步:随时间推移,房间真实温度波动变化) I, D' \( y, X- j) J7 e
    %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程5 }9 i- _) L9 \
    X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声9 Q+ K4 D8 N& e# t
    %第二步:随时间推移,获取实时数据
5 D7 _2 ~% a9 U3 W" S9 n    %温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,
- {) }% y8 T- n    %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响# a$ W& h8 h; \. u$ I# C
    %尽可能的逼近X(k),这也是Kalman滤波目的所在
' z9 Q8 ~/ n/ E& [8 J( r
2 [$ \& ?3 c% m- A& W6 h, O    %修改本次函数# D6 u$ m7 P+ j
    %Z(k) = H*X(k)+V(k);                %测量值为实际值叠加上测量噪声8 I3 o& r, J0 Z! J) Z8 `+ A
" ]9 ~6 w+ A+ K8 Z' s, Q  W
    %第三步:Kalman滤波
2 E! r% I! o/ \% z  i    %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,: i6 |. Q5 l8 I2 H' C
    %读者可以对照(3.36)到式(3.40),理解滤波过程
5 d/ r5 k" M9 e* f. U& C9 d" K    X_pre = F*Xkf(k - 1);                        %状态预测          X_pre为上一次卡尔曼滤波值* E5 V8 V7 k& O% X; [4 x' D: U
    P_pre = F*P(k - 1)*F + Q;                %协方差预测    ' c# J$ _3 c# D3 x, V. @( g& m
    %inv()为求一个方阵的逆矩阵。
/ y7 t' x( X7 B9 a7 Y* w    Kg = P_pre*inv(H * P_pre*H' + R);  %计算卡尔曼增益" ^$ c0 H3 n/ i, i" L) v
    e = Z(k) - H*X_pre;                           %新息                 本次测量值与上次预测值之差5 T# J1 `+ ?  A' A- ^9 _
    Xkf(k) = X_pre + Kg*e;                     %状态更新         本次预测值' E  X  Z. Z: h, j- ?6 i5 x0 g
    P(k) = (I - Kg*H)*P_pre;                    %协方差更新
8 o1 [) }8 S1 t+ M, U! a8 jend
! G  z+ @& F! x0 v% T%计算误差
' A2 M5 @6 l& l* N8 `) ZErr_Messure = zeros(1,N);   %测量值与真实值之间的偏差
: h$ h1 c3 P$ S0 |Err_Kalman = zeros(1,N);     %Kalman估计与真实值之间的偏差
  w0 h/ |: @2 c' N8 h6 ^2 ?for k = 1:N 6 c2 f- C) h) a* ^
    Err_Messure(k) = abs(Z(k) - X(k));    %abs()为求解绝对值
: T/ H5 a6 a6 n. p& [    Err_Kalman(k) = abs(Xkf(k) - X(k));
- R4 |) \/ O, ~% X5 {end
. x- `" s5 J% z1 Ot = 1:N;1 y; ~5 F" b6 z1 ~8 [
figure   %画图显示) p/ Z9 l  B5 L7 J, g& {! w* Z$ Z
%依次输出理论值,叠加过程噪声(受波动影响)的真实值5 K6 d- V, r- s6 L4 `1 w
%温度计测量值,Kalman估计值( A9 B- W4 A2 Q
plot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');
" X) D/ x$ x0 `+ K# s" plegend('期望值','真实值','观测值','Kalman滤波值');
4 z6 t" Q: U3 l: a- Yxlabel('采样时间/s');
7 F. C; T! [" |9 V" A0 v2 S0 Sylabel('温度值/C');
; p* J8 v4 O; ]2 u  _%误差分析图8 {% J) Z7 F" B/ e# S) B& `
figure  %画图显示
6 z" ]& m! |6 S9 n+ Dplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
! {' o: F# k1 Plegend('测量偏差','Kalman滤波偏差');- }. e% k& S: K2 O
xlabel('采样时间/s');& M( N8 \% Y( q" o: ]' d
ylabel('温度偏差值/C');' z1 k$ W1 v, G* s! h1 j
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0 v7 [  J% i$ j3 G6 W* c  i8 j* {( w' U; w
5 X! O1 A; H4 `- c9 Z

2 H3 R/ I" p# F- t# k1 y
: W2 Z- c1 x) I- J
5 F8 m4 L! e; m, E' v- V" q
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

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

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

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

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