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

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

[复制链接]

该用户从未签到

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

EDA365欢迎您登录!

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

x
   卡尔曼滤波实现简单,滤波效果好 ,下面分享一个基于卡尔曼滤波的matlab算法,数据全部为一维状态,本人弥补的详细备注,供爱好者研究学习。9 L$ \6 B, H  c' C- j

& U- I) {' Z, T/ |7 t9 a%%%%%%%%%%%%%%%%%%% ^2 f. A5 h( U- x2 G
%功能说明:Kalman滤波用在一维温度数据测量系统中( I* N5 J2 a3 m' F, V' E) g
function main2 k$ b, Y) h) i# D( M) }
%%%%%%8 y! G  F1 p. r5 M* B0 b
N = 120;                                    %一共采样的点数,时间单位是分钟,可理解为实验进行了60分钟的测量
7 o* b5 h2 P- \$ X  v/ K  |& iCON = 25;                                %室内温度的理论值,在这个理论值的基础上受过程噪声会有波动* m& n% n2 U0 D) s9 D3 O% r$ ?
%ones(a,b) 产生a行b列值为1的矩阵8 b2 M; G- e1 t- }
%zeros(a,b)产生a行b列值为0的矩阵
9 _: ~$ T7 B! D. `6 AXexpect = CON*ones(1,N);     %期望的温度是恒定的25度,但实际温度不可能会这样的' {2 n$ ?2 a; r( C5 T- n

4 T7 X. ?+ Q! h/ |! N5 WX = zeros(1,N);                         %房间各时刻真实温度值
% z; V5 M1 O: v- P, l( Q' w9 f. nXkf = zeros(1,N);                      %Kalman滤波处理的状态,也叫估计值
  D' x  F- x2 Y8 I0 n1 E" oZ = zeros(1,N);                         %温度计测量值
- G% [2 H; @$ N. h; V) vP = zeros(1,N);( e- C# c4 T) V$ C
%X(1)为数组的第一个元素
$ y2 y1 Y& X+ V7 S$ g3 s" x2 dX(1) = 25.1;     %假如初始值房间温度为25.1度1 K$ z: i' z! d+ p  d) ?
P(1) = 0.01;     %初始值的协方差   (测量值 - 真实值)^2
7 v  A" G3 J, R" n
  J- l* W# s2 X2 k/ n%产生0-1的随机数  模拟系统的随机数据
. `3 T) u& @1 Y: {' DZ = 24+rand(1,N)*10 - 5;       # O( S! b7 O1 k

$ K/ f; P5 M1 F% g3 N8 w8 @Z(1) = 24.9;     %温度计测的第一个数据
' F8 t$ O. _5 v7 J7 {* sXkf(1) = Z(1);  %初始测量值为24.9度,可以作为滤波器的初始估计状态噪声& |# U4 ?; ^" ?  W: Q& Z* _, |
0 t3 `* L  u* t, w* b$ M+ F: w
Q = 0;        %扰动误差方差(不存在扰动误差只有测量误差)8 O5 o1 J4 c) Y' ?
R = 0.25;        %测量误差方差' T1 e  j& {5 X6 Y; F8 r' D6 `' g
%sqrt(Q)*randn(1,N)为产生方差为0.01的随机信号
% ~/ n+ a! ]5 l3 ~, ]%W为过程噪声# p- m  J, F" u) u# q. f2 x
%V为测量噪声
6 K; P& F; X( X+ H8 a9 M+ c4 ]W = sqrt(Q)*randn(1,N); %方差决定噪声的大小$ k+ C: B. X% l, q2 E7 [  V" n
V = sqrt(R)*randn(1,N);%方差决定噪声的大小
* g2 }! a7 \& T. c9 @# B%系统矩阵
; z' h9 o: k6 O6 k$ L  [%解释:因为该系统的数据都是一维的,故各变换矩阵都是1,原因自己找书理解6 s. z6 g7 X/ W9 r
F = 1;     / m1 z, Z5 e6 N# a0 i; T9 f$ C
G = 1;- a6 b/ _- s" k" F
H = 1;
$ j- w1 Y% e1 o. n6 h1 f  M%eye产生m×n的单位矩阵 数值应该为1
+ M: x; e# ^; s1 ?. [I = eye(1);  %本系统状态为一维
) d; s% v: E0 |7 M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0 [( k* k2 J& L" Y7 C% |! ~%模拟房间温度和测量过程,并滤波
0 X2 z6 c( e1 ^; x# C0 O* Y9 _9 e( t1 T+ j/ \& Z% B- F
for k = 2:N/ c3 K" ~4 Y# x3 b. r
    %第一步:随时间推移,房间真实温度波动变化
% ?/ X8 s, _% F8 Y% c; @# O    %k时刻房间的真实温度,对于温度计来说,这个真实值是不知道的但是它的存在又是客观事实,读者要深刻领悟这个计算机模拟过程9 ^, m$ ]) D2 I1 B) d. Y3 b4 ]+ w
    X(k) = F*X(k - 1)+G*W(k - 1);%实际值为理想值叠加上扰动噪声- K. V! b: c5 D; e
    %第二步:随时间推移,获取实时数据4 l6 C- w3 z/ w# m+ \
    %温度计对K时刻房间温度的测量,Kalman滤波是站在温度计角度进行的,/ e8 W/ s  ]3 J& z6 y2 T  m
    %它不知道此刻的真实状态X(k),只能利用本次测量值Z(k)和上一次估计值Xkf(k)来做处理,其目标是最大限度地降低测量噪声R的影响
! l4 P% y! c2 h    %尽可能的逼近X(k),这也是Kalman滤波目的所在  a5 K: Z* l1 x& m( p5 k
8 K  i" B! t+ j* C5 |
    %修改本次函数
; d3 N+ m% h7 C8 c3 z6 e    %Z(k) = H*X(k)+V(k);                %测量值为实际值叠加上测量噪声$ E* I# J( e( {+ ~% `, F& O+ B
3 w8 u6 c. ?" i: {! i
    %第三步:Kalman滤波6 d1 _6 h8 {# T8 x# J) x
    %有了k时刻的观测Z(k)和k-1时刻的状态,那么就可以进行滤波了,
, g2 O; C0 r2 \! y2 H7 O    %读者可以对照(3.36)到式(3.40),理解滤波过程. C: E3 F* B; R* s% o" J! R5 K
    X_pre = F*Xkf(k - 1);                        %状态预测          X_pre为上一次卡尔曼滤波值6 q0 b) H8 z3 d4 Z
    P_pre = F*P(k - 1)*F + Q;                %协方差预测   
: k) s) L5 E' U- @    %inv()为求一个方阵的逆矩阵。' k( }. Q2 V/ n7 y
    Kg = P_pre*inv(H * P_pre*H' + R);  %计算卡尔曼增益
, u4 T) U( z) @    e = Z(k) - H*X_pre;                           %新息                 本次测量值与上次预测值之差
+ ]. X8 R* M+ h  b( r8 q) w0 x- R    Xkf(k) = X_pre + Kg*e;                     %状态更新         本次预测值* G' a0 E& J1 B
    P(k) = (I - Kg*H)*P_pre;                    %协方差更新
! u7 s9 f) V7 p$ W) Y8 e# Gend5 V' i( `0 k3 g9 F& d% |
%计算误差
5 e! r# @6 L# W& Z9 `Err_Messure = zeros(1,N);   %测量值与真实值之间的偏差8 D9 t& ~( l+ F1 S- B' p
Err_Kalman = zeros(1,N);     %Kalman估计与真实值之间的偏差
2 T# n+ A7 y) @9 ?% U1 r- ]for k = 1:N 4 D% N" e: e& ^# S
    Err_Messure(k) = abs(Z(k) - X(k));    %abs()为求解绝对值1 ]3 p. [* r) B* d/ o! H
    Err_Kalman(k) = abs(Xkf(k) - X(k));
( e  k! ], q1 d" _: oend& `1 z) f  l2 r( w2 V  Q7 ?/ T" ?
t = 1:N;
& S$ u4 e4 ]- R5 {8 u0 efigure   %画图显示  b6 }, ]. o- C' |5 \9 I5 E
%依次输出理论值,叠加过程噪声(受波动影响)的真实值. O, p+ W  P. P
%温度计测量值,Kalman估计值
" Q6 b/ x- Q0 @. L2 T, a2 mplot(t,Xexpect,'-b',t,X,'-r',t,Z,'-ko',t,Xkf,'-g*');6 i' P. g: ?( V+ _! Z
legend('期望值','真实值','观测值','Kalman滤波值');) z1 v7 K* q3 B
xlabel('采样时间/s');
$ ~2 w+ z5 i8 W* o0 m& yylabel('温度值/C');
8 L& _" N( e+ i7 G%误差分析图+ K. |! F9 @0 o& t( Z4 o' Y
figure  %画图显示
" E$ {% f( T6 M; r# l% Mplot(t,Err_Messure,'-b.',t,Err_Kalman,'-k*');
5 @5 j$ A' v9 n2 |8 }6 J, B5 H& ylegend('测量偏差','Kalman滤波偏差');
; [6 n' t  G$ Uxlabel('采样时间/s');2 s- N$ E3 J+ T2 T/ y' r; g& \; N. n
ylabel('温度偏差值/C');1 G" H: Z/ }! z. k
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
: H1 }* k6 }" a1 D

* ?! \4 B5 d# Z, }, B! V( E2 s  z+ O& j3 n' [# Y: u
. H( a7 y4 K+ ]8 j; `4 }# f, M
- m- ?7 R- g( G6 R, Y- O
9 V9 o2 U' x. |
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

EDA365公众号

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

GMT+8, 2025-11-23 17:57 , Processed in 0.140625 second(s), 23 queries , Gzip On.

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

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

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