文件一
% THIS PROGRAM IS FOR IMPLEMENTATION OF DISCRETE TIME PROCESS EXTENDED KALMAN FILTER
% FOR GAUSSIAN AND LINEAR STOCHASTIC DIFFERENCE EQUATION.
% By (R.C.R.C.R),SPLABS,MPL.
% (17 JULY 2005).
% Help by Aarthi Nadarajan is acknowledged.
% (drawback of EKF is when nonlinearity is high, we can extend the % approximation taking additional terms in Taylor's series).
clc; close all; clear all;
Xint_v = [1; 0; 0; 0; 0];
wk = [1 0 0 0 0];
vk = [1 0 0 0 0];
for ii = 1:1:length(Xint_v)
Ap(ii) = Xint_v(ii)*2;
W(ii) = 0;
H(ii) = ‐sin(Xint_v(ii));
V(ii) = 0;
Wk(ii) = 0;
end
Uk = randn(1,200);
Qu = cov(Uk);
Vk = randn(1,200);
Qv = cov(Vk);
C = [1 0 0 0 0];
n = 100;
[YY XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V);
for it = 1:1:length(XX)
MSE(it) = YY(it) ‐ XX(it);
end
tt = 1:1:length(XX);
figure(1); subplot(211); plot(XX); title('ORIGINAL SIGNAL'); subplot(212); plot(YY); title('ESTIMATED SIGNAL');
figure(2); plot(tt,XX,tt,YY); title('Combined plot');
legend('original','estimated');
figure(3); plot(MSE.^2); title('Mean square error');
子文件：：function [YY,XX] =
EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V);
Ap(2,:) = 0;
for ii = 1:1:length(Ap)‐1
Ap(ii+1,ii) = 1;
end
inx = 1;
UUk = [Uk(inx); 0; 0; 0; 0];
PPk = (Xint_v*Xint_v');
VVk = [Vk(inx); 0; 0; 0; 0];
Qv = V*V';
for ii = 1:1:length(Xint_v)
XKk(ii,1) =
Xint_v(ii)^2; % FIRST STEP end
PPk = Ap*PPk*Ap'; % SECOND STEP
Kk = PPk*C'*inv( (C*PPk*C') + (V*Qv*V') ); % THIRD STEP
for ii = 1:1:length(Xint_v)
XUPK(ii,1) = XKk(ii)^2 + UUk(ii); % UPPER EQUATIONS.
Zk(ii,1) = cos(XUPK(ii)) +
VVk(ii); % UPPER EQUATIONS.
end
for ii = 1:1:length(XKk)
XBARk(ii,1) = XKk(ii) + Kk(ii)*(Zk(ii) ‐
(cos(XKk(ii)))) ; % FOURTH STEP
end
II = eye(5,5);
Pk = ( II ‐ Kk*C)*PPk; % FIFTH STEP
%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
for ii = 1:1:n
UUk = [Uk(ii+1); 0; 0; 0; 0];
PPk = XBARk*XBARk';
VVk = [Vk(ii+1); 0; 0; 0; 0];
XKk = exp(‐XBARk); % FIRST STEP
PPkM = Ap*PPk*Ap'; % SECOND STEP
Kk = PPkM*C'*inv( (C*PPkM*C') + (V*Qv*V') ); % THIRD STEP
for nn = 1:1:length(XBARk)
XUPK(nn) = exp(‐XKk(nn)) + UUk(nn); % UPPER EQUATIONS.
Zk(nn) = cos(XUPK(nn)) + VVk(nn); % UPPER EQUATIONS.
end
for in = 1:1:length(XUPK)
XNEW(in) = XBARk(in) + Kk(in)*(Zk(in) ‐ cos(XBARk(in))); % FOURTH STEP
end
II = eye(5,5);
Pk = (II ‐ Kk*C)*PPkM; % FIFTH STEP
XBARk = XNEW;
OUTX(ii) = XBARk(1,1);
OUTY(ii) = Zk(1,1);
end
YY = OUTY;
XX = OUTX;
一个应用实例详解卡尔曼滤波及其算法实现 标签：算法filtermatlabalgorithm优化工作 2012-05-14 10:48 75511人阅读评论(25) 收藏举报分类： 数据结构及其算法（4） 为了可以更加容易的理解卡尔曼滤波器，这里会应用形象的描述方法来讲解，而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。 在介绍他的5条公式之前，先让我们来根据下面的例子一步一步的探索。 假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就是下一分钟的温度等于现在这一分钟的温度（假设我们用一分钟来做时间单位）。假设你对你的经验不是100%的相信，可能会有上下偏差几度。 我们把这些偏差看成是高斯白噪声（White Gaussian Noise），也就是这些偏差跟前后时间是没有关系的而且符合高斯分配（Gaussian Distribution）。另外，我们在房间里放一个温度计，但是这个温度计也不准确的，测量值会比实际值偏差。我们也把这些偏差看成是高斯白噪声。 好了，现在对于某一分钟我们有两个有关于该房间的温度值：你根据经验的预测值（系统的预测值）和温度计的值（测量值）。下面我们要用这两个值结合他们各自的噪声来估算出房间的实际温度值。 假如我们要估算k时刻的是实际温度值。首先你要根据k-1时刻的温度值，来预测k时刻的温度。因为你相信温度是恒定的，所以你会得到k时刻的温度预测值是跟k-1时刻一样的，假设是23度，同时该值的高斯噪声的偏差是5度（5是这样得到的：如果k-1时刻估算出的最优温度值的偏差是3，你对自己预测的不确定度是4度，他们平方相加再开方，就是5）。然后，你从温度计那里得到了k时刻的温度值，假设是25度，同时该值的偏差是4度。 由于我们用于估算k时刻的实际温度有两个温度值，分别是23 度和25度。究竟实际温度是多少呢？相信自己还是相信温度计呢？究竟相信谁多一点，我们可以用他们的covariance（协方差）来判断。因为Kg^2=5^2/(5^2+4^2)，所以Kg=0.78，我们可以估算出k时刻的实际温度值是：23+0.78*(25-23)=24.56度。 可以看出，因为温度计的covariance比较小（比较相信温度计），所以估算出的最优温度值偏向温度计的值。 现在我们已经得到k时刻的最优温度值了，下一步就是要进入k+1时刻，进行新的最优估算。到现在为止，好像还没看到什么自回归的东西出现。对了，在进入k+1时刻之前，我们还要算出k时刻那个最优值（24.56 度）的偏差。算法如下：((1-Kg)*5^2)^0.5=2.35。这里的5就是上面的k时刻你预测的那个23度
ardupilot(EKF)扩展卡尔曼滤波 一、初识卡尔曼滤波器 为了描述方便我从网上找了一张卡尔曼滤波器的5大公式的图片。篇幅所限，下图所示的是多维卡尔曼滤波器（因为EKF2是多维扩展卡尔曼滤波器，所以我们从多维说起），为了跟好的理解卡尔曼滤波器可以百度一下，从一维开始。 这5个公式之外还有一个观测模型，根据你实际的观测量来确定，它的主 要作用是根据实际情况来求观测矩阵H。 因为卡尔曼滤波器是线性滤波器，状态转移矩阵A和观测矩阵H是确定的。在维基百科上状态转移矩阵用F表示。在ardupilot EKF2算法中，状态转移矩阵也是用F表示的。下面是维基百科给出的线性卡尔曼滤波器的相关公式。
上述更新（后验）估计协方差的公式对任何增益K k都有效，有时称为约瑟夫形式。为了获得最佳卡尔曼增益，该公式进一步简化为P k|k=(I-K k H k)P k|k-1，它在哪种形式下应用最广泛。但是，必须记住它仅对最小化残差误差的最佳增益有效。 为了使用卡尔曼滤波器来估计仅给出一系列噪声观测过程的内部状态，必须根据卡尔曼滤波器的框架对过程进行建模，这意味着指定一下矩阵：
只要记住一点就行了，卡尔曼滤波器的作用就是输入一些包含噪声的数据，得到一些比较接近真是情况的数据。比如无人机所使用的陀螺仪和加速度计的 读值，他们的读值都是包含噪声的，比如明明真实的角速度是俯仰2°/s,陀螺 仪的读值却是2.5°/s。通过扩展卡尔曼之后的角速度值会变得更加接近2o/s 的真实值，有可能是2.1o/s。 二、扩展卡尔曼滤波器 因为卡尔曼滤波器针对的是线性系统，状态转移模型（说的白话一点就是知道上一时刻被估计量的值，通过状态转移模型的公式可以推算出当前时刻被 估计量的值）和观测模型。注：有的资料显示状态模型中有,有的没有，目前 我也不清楚是为什么，有可能和被估计的对象有关。但看多了你就会发现不管 网上给的公式有怎样的不同，但总体的流程是一样的，都是这5大步骤。我个 人觉得维基百科给的公式较为标准。 因为扩展卡尔曼滤波器（EKF，Extended Kalman filter）的使用场景为非线性系统。所以上面两公式改写为下面所示的样子，我个人的理解是，因为是 非线性系统，所以没有固定的状态转移矩阵和观测矩阵。到这儿为止卡尔曼滤 波器到扩展卡尔曼滤波器的过度就完成了（多说一句，因为传感器的数据采样 是有时间间隔的，算法的运行也是有间隔的，所以本文提到的KF和EKF都是离散型的）。下面是扩展卡尔曼滤波器的相关公式。
文件一 % THIS PROGRAM IS FOR IMPLEMENTATION OF DISCRETE TIME PROCESS EXTENDED KALMAN FILTER % FOR GAUSSIAN AND LINEAR STOCHASTIC DIFFERENCE EQUATION. % By (R.C.R.C.R),SPLABS,MPL. % (17 JULY 2005). % Help by Aarthi Nadarajan is acknowledged. % (drawback of EKF is when nonlinearity is high, we can extend the % approximation taking additional terms in Taylor's series). clc; close all; clear all; Xint_v = [1; 0; 0; 0; 0]; wk = [1 0 0 0 0]; vk = [1 0 0 0 0]; for ii = 1:1:length(Xint_v) Ap(ii) = Xint_v(ii)*2; W(ii) = 0; H(ii) = ‐sin(Xint_v(ii)); V(ii) = 0; Wk(ii) = 0; end Uk = randn(1,200); Qu = cov(Uk); Vk = randn(1,200); Qv = cov(Vk); C = [1 0 0 0 0]; n = 100; [YY XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); for it = 1:1:length(XX) MSE(it) = YY(it) ‐ XX(it); end tt = 1:1:length(XX); figure(1); subplot(211); plot(XX); title('ORIGINAL SIGNAL'); subplot(212); plot(YY); title('ESTIMATED SIGNAL'); figure(2); plot(tt,XX,tt,YY); title('Combined plot'); legend('original','estimated'); figure(3); plot(MSE.^2); title('Mean square error'); 子文件：：function [YY,XX] = EKLMNFTR1(Ap,Xint_v,Uk,Qu,Vk,Qv,C,n,Wk,W,V); Ap(2,:) = 0; for ii = 1:1:length(Ap)‐1 Ap(ii+1,ii) = 1;
卡尔曼滤波入门： 卡尔曼滤波是用来进行数据滤波用的，就是把含噪声的数据进行处理之后得出相对真值。卡尔曼滤波也可进行系统辨识。 卡尔曼滤波是一种基于统计学理论的算法，可以用来对含噪声数据进行在线处理，对噪声有特殊要求，也可以通过状态变量的增广形式实现系统辨识。 用上一个状态和当前状态的测量值来估计当前状态，这是因为上一个状态估计此时状态时会有误差，而测量的当前状态时也有一个测量误差，所以要根据这两个误差重新估计一个最接近真实状态的值。 信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，我们需要寻找一种所谓有最佳线性过滤特性的滤波器。这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。 维纳(Wiener)滤波与卡尔曼(Kalman)滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。 (1)过滤或滤波 - 从当前的和过去的观察值x(n)，x(n-1)，x(n-2)，…估计当前的信号值称为过滤或滤波; (2)预测或外推 - 从过去的观察值，估计当前的或将来的信号值称为预测或外推; (3)平滑或内插 - 从过去的观察值，估计过去的信号值称为平滑或内插; 因此，维纳过滤与卡尔曼过滤又常常被称为最佳线性过滤与预测或线性最优估计。这里所谓“最佳”与“最优”是以最小均方误差为准则的。 维纳过滤与卡尔曼过滤都是解决最佳线性过滤和预测问题，并且都是以均方误差最小为准则的。因此在平稳条件下，它们所得到的稳态结果是一致的。然而，它们解决的方法有很大区别。 维纳过滤是根据全部过去的和当前的观察数据来估计信号的当前值，它的解是以均方误差最小条件下所得到的系统的传递函数H(z)或单位样本响应h(n)的形式给出的，因此更常称这种系统为最佳线性过滤器或滤波器。 而卡尔曼过滤是用前一个估计值和最近一个观察数据(它不需要全部过去的观察数据)来估计信号的当前值，它是用状态方程和递推的方法进行估计的，它的解是以估计值(常常是状态变量值)形式给出的。因此更常称这种系统为线性最优估计器或滤波器。 维纳滤波器只适用于平稳随机过程，而卡尔曼滤波器却没有这个限制。维纳过滤中信号和噪声是用相关函数表示的，因此设计维纳滤波器要求已知信号和噪声的相关函数。 卡尔曼过滤中信号和噪声是状态方程和量测方程表示的，因此设计卡尔曼滤波器要求已知状态方程和量测方程(当然，相关函数与状态方程和量测方程之间会存在一定的关系。卡尔曼过滤方法看来似乎比维纳过滤方法优越，它用递推法计算，不需要知道全部过去的数据，从而运用计算机计算方便，而且它可用于平稳和不平稳的随机过程(信号)，非时变和时变的系统。 但从发展历史上来看维纳过滤的思想是40年代初提出来的，1949年正式以书的形式出版。卡尔曼过滤到60年代初才提出来，它是在维纳过滤的基础上发展起来的，虽然如上所述它比维纳过滤方法有不少优越的地方，但是最佳线性过滤问题是由维纳过滤首先解决的，维纳过滤的物理概念比较清楚，也可以认为卡尔曼滤波仅仅是对最佳线性过滤问题提出的一种新的算法。 卡尔曼滤波在数学上是一种统计估算方法，通过处理一系列带有误差的实际量测数据而得到的物理参数的最佳估算。例如在气象应用上，根据滤波的基本思想，利用前一时刻预报误差的反馈信息及时修正预报方程，以提高下一时刻预报精度。作温度预报一般只需要连续两个月的资料即可建立方程和递推关系。
clear all; %% 惯性-GPS组合导航模型参数初始化 we = 360/24/60/60*pi/180; %地球自转角速度，弧度/s psi = 10*pi/180; %psi角度/ 弧度 Tge = 0.12; Tgn = 0.10; Tgz = 0.10; %这三个参数的含义详见参考文献 sigma_ge=1; sigma_gn=1; sigma_gz=1; %% 连续空间系统状态方程 % X_dot(t) = A(t)*X(t) + B(t)*W(t) A=[0 we*sin(psi) -we*cos(psi) 1 0 0 1 0 0; -we*sin(psi) 0 0 0 1 0 0 1 0; we*cos(psi) 0 0 0 0 1 0 0 1; 0 0 0 -1/Tge 0 0 0 0 0; 0 0 0 0 -1/Tgn 0 0 0 0; 0 0 0 0 0 -1/Tgz 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0;]; %状态转移矩阵 B=[0 0 0 sigma_ge*sqrt(2/Tge) 0 0 0 0 0; 0 0 0 0 sigma_gn*sqrt(2/Tgn) 0 0 0 0; 0 0 0 0 0 sigma_gz*sqrt(2/Tgz) 0 0 0;]';%输入控制矩阵%% 转化为离散时间系统状态方程 % X(k+1) = F*X(k) + G*W(k) T = 0.1; [F,G]=c2d(A,B,T);
H=[1 0 0 0 0 0 0 0 0; 0 -sec(psi) 0 0 0 0 0 0 0;];%观测矩阵 %% 卡尔曼滤波器参数初始化 t=0:T:50-T; length=size(t,2); y=zeros(2,length); Q=0.5^2*eye(3); %系统噪声协方差 R=0.25^2*eye(2); %测量噪声协方差 y(1,:)=2*sin(pi*t*0.5); y(2,:)=2*cos(pi*t*0.5); Z=y+sqrt(R)*randn(2,length); %生成的含有噪声的假定观测值，2维X=zeros(9,length); %状态估计值，9维 X(:,1)=[0,0,0,0,0,0,0,0,0]'; %状态估计初始值设定 P=eye(9); %状态估计协方差 %% 卡尔曼滤波算法迭代过程 for n=2:length X(:,n)=F*X(:,n-1); P=F*P*F'+ G*Q*G'; Kg=P*H'/(H*P*H'+R); X(:,n)=X(:,n)+Kg*(Z(:,n)-H*X(:,n)); P=(eye(9,9)-Kg*H)*P; end %% 绘图代码 figure(1) plot(y(1,:)) hold on; plot(y(2,:)) hold off; title('理想的观测量'); figure(2)
卡尔曼滤波器及其简matlab仿真
卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems （线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。
卡尔曼滤波器及其简matlab仿真 一、卡尔曼滤波的起源 谈到信号的分析与处理，就离不开滤波两个字。通常，信号的频谱处于有限的频率范围内，而噪声的频谱则散布在很广的频率范围内，为了消除噪声，可以进行频域滤波。但在许多应用场合，需要直接进行时域滤波，从带噪声的信号中提取有用信号。虽然这样的过程其实也算是对信号的滤波，但其所依据的理论，即针对随机信号的估计理论，是自成体系的。人们对于随机信号干扰下的有用信号不能“确知”，只能“估计”。为了“估计”，要事先确定某种准则以评定估计的好坏程度。 1960年卡尔曼发表了用递归方法解决离散数据线性滤波问题的论文A New Approach to Linear Filtering and Prediction Problems（线性滤波与预测问题的新方法），在这篇文章里一种克服了维纳滤波缺点的新方法被提出来，这就是我们今天称之为卡尔曼滤波的方法。卡尔曼滤波应用广泛且功能强大，它可以估计信号的过去和当前状态甚至能估计将来的状态即使并不知道模型的确切性质。 其基本思想是以最小均方误差为最佳估计准则，采用信号与噪声的状态空间模型利用前一时刻的估计值和当前时刻的观测值来更新对状态变量的估计，求出当前时刻的估计值。算法根据建立的系统方程和观测方程对需要处理的信号做出满足最小均方误差的估计。 对于解决很大部分的问题，它是最优，效率最高甚至是最有用的。它的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻的储值的估计，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。 二、卡尔曼滤波的原理
扩展卡尔曼滤波器教程 在使用OpenPilot和Pixhawk飞控时，经常遇到扩展卡尔曼滤波(EKF)。从不同的网页和参考论文中搜索这个词，其中大部分都太深奥了。所以我决定创建自己学习教程。本教程从一些简单的例子和标准(线性)卡尔曼滤波器，通过对实际例子来理解卡尔曼滤波器。 Part 1: 一个简单的例子 想象一个飞机准备降落时，尽管我们可能会担心许多事情，像空速、燃料、等等，当然最明显是关注飞机的高度(海拔高度)。通过简单的近似，我们可以认为当前高度是之前的高度失去了一小部分。例如，当每次我们观察飞行高度时，认为飞机失去了2%的高度，那么它的当前高度是上一时刻高度的98%: altitude current_time=0.98*altitude previous_time 工程上对上面的公式，使用“递归”这个术语进行描述。通过递归前一时刻的值，不断计算当前值。最终我们递归到初始的“基本情况”，比如一个已知的高度。 试着移动上面的滑块，看看飞机针对不同百分比的高度变化。 Part 2:处理噪声 当然, 实际从传感器比如GPS或气压计获得测量高度时，传感器的数据或多或少有所偏差。如果传感器的偏移量为常数，我们可以简单地添加或减去这偏移量来确定我们的高度。不过通常情况下，传感器的偏移量是一个时变量，使得我们所观测到的传感器数据相当于实际高度加上噪声: observed_altitude current_time=altitude current_time+noise current_time 试着移动上面的滑块看到噪声对观察到的高度的影响。噪音被表示为可观测的海拔范围的百分比。
clear all v=150; %%目标速度 v_sensor=0;%%传感器速度 t=1; %%扫描周期 xradarpositon=0; %%传感器坐标yradarpositon=0; %% ppred=zeros(4,4); Pzz=zeros(2,2); Pxx=zeros(4,2); xpred=zeros(4,1); ypred=zeros(2,1); sumx=0; sumy=0; sumxukf=0; sumyukf=0; sumxekf=0; sumyekf=0; %%%统计的初值 L=4; alpha=1; kalpha=0; belta=2; ramda=3-L; azimutherror=0.015; %%方位均方误差rangeerror=100; %%距离均方误差processnoise=1; %%过程噪声均方差 tao=[t^3/3 t^2/2 0 0; t^2/2 t 0 0; 0 0 t^3/3 t^2/2; 0 0 t^2/2 t]; %% the input matrix of process G=[t^2/2 0 t 0 0 t^2/2 0 t ]; a=35*pi/180; a_v=5/100; a_sensor=45*pi/180; x(1)=8000; %%初始位置
y(1)=12000; for i=1:200 x(i+1)=x(i)+v*cos(a)*t; y(i+1)=y(i)+v*sin(a)*t; end for i=1:200 xradarpositon=0; yradarpositon=0; Zmeasure(1,i)=atan((y(i)-yradarpositon)/(x(i)-xradarpositon))+random('Normal',0,azimutherror,1,1); Zmeasure(2,i)=sqrt((y(i)-yradarpositon)^2+(x(i)-xradarpositon)^2)+random('Normal',0,rangeerror,1,1); xx(i)=Zmeasure(2,i)*cos(Zmeasure(1,i));%%观测值 yy(i)=Zmeasure(2,i)*sin(Zmeasure(1,i)); measureerror=[azimutherror^2 0;0 rangeerror^2]; processerror=tao*processnoise; vNoise = size(processerror,1); wNoise = size(measureerror,1); A=[1 t 0 0; 0 1 0 0; 0 0 1 t; 0 0 0 1]; Anoise=size(A,1); for j=1:2*L+1 Wm(j)=1/(2*(L+ramda)); Wc(j)=1/(2*(L+ramda)); end Wm(1)=ramda/(L+ramda); Wc(1)=ramda/(L+ramda);%+1-alpha^2+belta; %%%权值 if i==1 xerror=rangeerror^2*cos(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*sin(Zmeasure(1,i))^2; yerror=rangeerror^2*sin(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*cos(Zmeasure(1,i))^2; xyerror=(rangeerror^2-Zmeasure(2,i)^2*azimutherror^2)*sin(Zmeasure(1,i))*cos(Zmeasure(1,i)); P=[xerror xerror/t xyerror xyerror/t; xerror/t 2*xerror/(t^2) xyerror/t 2*xyerror/(t^2); xyerror xyerror/t yerror yerror/t;
卡尔曼滤波器综述 瞿伟军 G10074 1、卡尔曼滤波的起源 1960年，匈牙利数学家卡尔曼发表了一篇关于离散数据线性滤波递推算法的论文，这意味着卡尔曼滤波的诞生。斯坦利.施密特(Stanley Schmidt)首次实现了卡尔曼滤波器，卡尔曼在NASA埃姆斯研究中心访问时，发现他的方法对于解决阿波罗计划的轨道预测很有用，后来阿波罗飞船的导航电脑使用了这种滤波器。关于这种滤波器的论文由Swerling (1958)、Kalman (1960)与 Kalman and Bucy (1961)发表。 2、卡尔曼滤波的发展 卡尔曼滤波是一种有着相当广泛应用的滤波方法，但它既需要假定系统是线性的，又需要认为系统中的各个噪声与状态变量均呈高斯分布，而这两条并不总是确切的假设限制了卡尔曼滤波器在现实生活中的应用。扩展卡尔曼滤波器（EKF）极大地拓宽了卡尔曼滤波的适用范围。EKF的基本思路是，假定卡尔曼滤滤对当前系统状态估计值非常接近于其真实值，于是将非线性函数在当前状态估计值处进行台劳展开并实现线性化。另一种非线性卡尔曼滤波叫线性化卡尔曼滤波。它与EKF的主要区别是前者将非线函数在滤波器对当前系统状态的最优估计值处线性化，而后者因为预先知道非线性系统的实际运行状态大致按照所要求、希望的轨迹变化，所以这些非线性化函数在实际状态处的值可以表达为在希望的轨迹处的台劳展开式，从而完成线性化。 不敏卡尔曼滤波器（UKF）是针对非线性系统的一种改进型卡尔曼滤波器。UKF处理非线性系统的基本思路在于不敏变换，而不敏变换从根本上讲是一种描述高斯随机变量在非线性化变换后的概率分布情况的方法。不敏卡尔曼滤波认为，与其将一个非线性化变换线性化、近似化，还不如将高斯随机变量经非线性变换后的概率分布情况用高斯分布来近似那样简单，因而不敏卡尔曼滤波算法没
% 扩展卡尔曼滤波器估计单相电压幅值、相位、频率参数（含直流）function test2_EKF close all; clc; tic; %计时 %模型：y=A0+A1*cos(omega*t+phy1) %离散化：y(k)=A0(k)+A1(k)*cos(omega(k)*k*Ts+phy1(k)) %状态变量：x1(k)=A0(k),x2(k)=omega(k),x3(k)=A1(k)*cos(omega(k)*k*Ts+phy1(k) ),x4(k)=A1(k)*sin(omega(k)*k*Ts+phy1(k)) %下一时刻状态变量为(假设状态不突变)：A0(k+1)=A0(k),A1(k+1)=A1(k),omega(k+1)=omega(k),phy1(k+1)=phy1 (k); %则对应状态为：x1(k+1)=x1(k),x2(k+1)=x2(k),x3(k+1)=x3(k)*cos(x2(k)*Ts)- x4(k)*sin(x(2)*Ts),x4(k+1)=x3(k)*sin(x2(k)*Ts)+x4(k)*cos(x(2)*Ts); %状态空间描述：X(k+1)=f(X(k))+W(k);y(k)=H*X(k)+v(k) %f(X(k))=[x1(k);x2(k);x3(k)*cos(x2(k)*Ts)- x4(k)*sin(x(2)*Ts);x3(k)*sin(x2(k)*Ts)+x4(k)*cos(x(2)*Ts)] %偏导(只求了三个)：f`(X(k))=[1,0,0;0,1,0;0,-x3(k)*Ts*sin(x2(k)*Ts)-x4(k)*Ts*cos(x2(k)*Ts),cos(x2(k)*Ts);0,x3(k)*Ts*cos(x2(k)*Ts)- x4(k)*Ts*sin(x2(k)*Ts),sin(x2(k)*Ts)]
卡尔曼滤波简介及其算法实现代码 卡尔曼滤波算法实现代码（C，C＋＋分别实现） 卡尔曼滤波器简介 近来发现有些问题很多人都很感兴趣。所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。现在先讨论一下卡尔曼滤波器，如果时间和能力允许，我还希望能够写写其他的算法，例如遗传算法，傅立叶变换，数字滤波，神经网络，图像处理等等。 因为这里不能写复杂的数学公式，所以也只能形象的描述。希望如果哪位是这方面的专家，欢迎讨论更正。 卡尔曼滤波器– Kalman Filter 1．什么是卡尔曼滤波器 （What is the Kalman Filter?） 在学习卡尔曼滤波器之前，首先看看为什么叫“卡尔曼”。跟其他著名的理论（例如傅立叶变换，泰勒级数等等）一样，卡尔曼也是一个人的名字，而跟他们不同的是，他是个现代人！ 卡尔曼全名Rudolf Emil Kalman，匈牙利数学家，1930年出生于匈牙利首都布达佩斯。1953，1954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器，正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》（线性滤波与预测问题的新方法）。如果对这编论文有兴趣，可以到这里的地址下载： https://www.doczj.com/doc/093744275.html,/~welch/media/pdf/Kalman1960.pdf。 简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 2．卡尔曼滤波器的介绍 （Introduction to the Kalman Filter） 为了可以更加容易的理解卡尔曼滤波器，这里会应用形象的描述方法来讲解，而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。 在介绍他的5条公式之前，先让我们来根据下面的例子一步一步的探索。 假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就
南京航空航天大学 随机信号小论文题目扩展卡尔曼滤波 学生姓名梅晟 学号SX1504059 学院电子信息工程学院 专业通信与信息系统
扩展卡尔曼滤波 一、引言 20世纪60年代，在航空航天工程突飞猛进而电子计算机又方兴未艾之时，卡尔曼发表了论文《A New Approach to Linear Filtering and Prediction Problems》（一种关于线性滤波与预测问题的新方法），这让卡尔曼滤波成为了时域内有效的滤波方法，从此各种基于卡尔曼滤波的方法横空出世，在目标跟踪、故障诊断、计量经济学、惯导系统等方面得到了长足的发展。 二、卡尔曼滤波器 卡尔曼滤波是一种高效率的递归滤波器(自回归滤波器), 它能够从一系列的不完全及包含噪声的测量中，估计动态系统的状态。卡尔曼滤波的一个典型实例是从一组有限的，包含噪声的，对物体位置的观察序列（可能有偏差）预测出物体的位置的坐标及速度。 卡尔曼在NASA埃姆斯研究中心访问时，发现他的方法对于解决阿波罗计划的轨道预测很有用，后来阿波罗飞船的导航电脑便使用了这种滤波器。目前,卡尔曼滤波已经有很多不同的实现。卡尔曼最初提出的形式现在一般称为简单卡尔曼滤波器。除此以外，还有施密特扩展滤波器、信息滤波器以及很多Bierman, Thornton 开发的平方根滤波器的变种。也许最常见的卡尔曼滤波器是锁相环，它在收音机、计算机和几乎任何视频或通讯设备中广泛存在。 三、扩展卡尔曼滤波器 3.1 被估计的过程信号 卡尔曼最初提出的滤波理论只适用于线性系统,Bucy，Sunahara等人提出并研究了扩展卡尔曼滤波(Extended Kalman Filter，简称EKF)，将卡尔曼滤波理论进一步应用到非线性领域。EKF的基本思想是将非线性系统线性化，然后进行卡尔曼滤波，因此EKF是一种次优滤波。 同泰勒级数类似，面对非线性关系时，我们可以通过求过程和量测方程的偏导来线性化并计算当前估计。假设过程具有状态向量x∈?n，其状态方程为非线性随机差分方程的形式。 x k=f x k?1,u k?1,w k?1(1.1) 观测变量z∈?m为： z k=?(x k,v k)(1.2) 随机变量w k和v k代表过程激励噪声和观测噪声。它们为相互独立，服从正态分布的白色噪声：
设高度的测量误差是均值为0、方差为1的高斯白噪声随机序列，该物体的初始高度0h 和速度0V 也是高斯分布的随机变量，且0000019001000,var 10/02Eh h m P EV m s V ???????? ===? ??????? ???? ????。试求该物体高度和速度随时间变化的最优估计。（2/80.9s m g =） 解： 1. 令()()()h k X k v k ?? =? ??? t=1 R （k ）=1 Q(k)=0 根据离散时间卡尔曼滤波公式，则有： (1)(1,)()()X k k k X k U k φ+=++ (1)(1)(1)(1)Y k H k X k V k +=++++ (1,)k k φ+＝ 11t -?? ??? ? ()U k ＝ 20.5gt gt ??-???? (1)H k +＝[]10 滤波初值：^ 1900(0|0)(0)10X EX ?? ==???? 0100(0|0)var[(0)]2P X P ?? ===? ??? 一步预测：^^ (1|)(1,)(|)()X k k k k X k k U k φ+=++ (1|)(1,)(|)(1,)T P k k k k P k k k k φφ+=++ 滤波增益：1 (1)(1|)(1)[(1)(1|)(1)(1)]T T K k P k k H k H k P k k H k R k -+=+++++++ 滤波计算：^ ^ ^ (1|1)(1|)(1)[(1)(1)(1|)]X k k X k k K k Y k H k X k k ++=++++-++ (1|1)[(1)(1)](1|)P k k I K k H k P k k ++=-+++ 2. 实验结果
MATLAB 对卡尔曼滤波器的仿真实现 刘丹，朱毅，刘冰 武汉理工大学信息工程学院，武汉（430070） E-mail ：liudan_ina@https://www.doczj.com/doc/093744275.html, 摘 要：本文以卡尔曼滤波器原理为理论基础，用MATLAB 进行卡尔曼滤波器仿真、对比卡尔曼滤波器的预测效果，对影响滤波其效果的各方面原因进行讨论和比较，按照理论模型进行仿真编程，清晰地表述了编程过程。 关键词：数字信号处理；卡尔曼滤波器；MATLAB ；仿真过程 中图分类号: TN912.3 1. 引言 随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的学科和技术领域。数字信号处理已在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理中,数字滤波占有极其重要的地位,目前对数字滤波器的设计有多种方法,其中著名的MATLAB 软件包在多个研究领域都有着广泛的应用,它的频谱分析[1]和滤波器的分析设计功能很强,从而使数字信号处理变得十分简单、直观。本文分析了数字滤波器的设计方法,举出了基于MATLAB 软件的信号处理工具在数字滤波器设计中的应用。 2. 卡尔曼滤波基本原理 卡尔曼滤波过程实际上是获取维纳解的递推运算过程[2]。从维纳解导出的卡尔曼滤波器实际上是卡尔曼滤波过程结束后达到稳态的情况，这时Kalman Filtering 的结果与Wiener Solution 是相同的[3]。具体推导如下： )()1|1(?)|(?n Gy n n x f n n x +??= )|(?)()(n n x n x n e ?= 已知由此求c a cG a f F G n e E n ,)1(( ..min )]([)(2?=??→?==ε 由 f G f G ,0??????????=??εε ⑴ )]1|1(?)()[()1|1(?)|(????+??=n n x ac n y n G n n x a n n x 可以是时变的，非平稳的随机信号 ⑵ Q n a n P +?=)1()(2 ε均为正数。 ⑶ ) () ()(2n P C R n CP n G += ⑷ )()](1[)()(n P n CG n G C P n ??== ε )(n G 是个随时间变化的量，每次输入输出，)(n G 就调整一次，并逐渐逼近Kalman Filter 的增益G ，而)1()(? 扩展卡尔曼滤波原理： 在原有卡尔曼滤波的基础上，为了解决多目标值的跟踪与估计，形成了扩展卡尔曼滤波。起matlab主要程序如下： clear all v=150; %%目标速度 v_sensor=0;%%传感器速度 t=1; %%扫描周期 xradarpositon=0; %%传感器坐标 yradarpositon=0; %% ppred=zeros(4,4); Pzz=zeros(2,2); Pxx=zeros(4,2); xpred=zeros(4,1); ypred=zeros(2,1); sumx=0; sumy=0; sumxukf=0; sumyukf=0; sumxekf=0; sumyekf=0; %%%统计的初值 %滤波算法描述： L=4; alpha=1; kalpha=0; belta=2; ramda=3-L; azimutherror=0.015; %%方位均方误差 rangeerror=100; %%距离均方误差 processnoise=1; %%过程噪声均方差 tao=[t^3/3 t^2/2 0 0; t^2/2 t 0 0; 0 0 t^3/3 t^2/2; 0 0 t^2/2 t]; %% the input matrix of process G=[t^2/2 0 t 0 0 t^2/2 0 t ]; a=35*pi/180; a_v=5/100; a_sensor=45*pi/180; x(1)=8000; %%初始位置 y(1)=12000; for i=1:200 x(i+1)=x(i)+v*cos(a)*t; y(i+1)=y(i)+v*sin(a)*t; end for i=1:200 xradarpositon=0; yradarpositon=0; Zmeasure(1,i)=atan((y(i)-yradarpositon)/(x(i)-xradarpositon))+random('Normal',0,azimutherror,1,1); Zmeasure(2,i)=sqrt((y(i)-yradarpositon)^2+(x(i)-xradarpositon)^2)+random('Normal',0,rangeerror,1,1); xx(i)=Zmeasure(2,i)*cos(Zmeasure(1,i));%%观测值 yy(i)=Zmeasure(2,i)*sin(Zmeasure(1,i)); measureerror=[azimutherror^2 0;0 rangeerror^2]; processerror=tao*processnoise; vNoise = size(processerror,1); wNoise = size(measureerror,1); A=[1 t 0 0; 0 1 0 0; 0 0 1 t; 0 0 0 1]; Anoise=size(A,1); for j=1:2*L+1 Wm(j)=1/(2*(L+ramda)); Wc(j)=1/(2*(L+ramda)); end Wm(1)=ramda/(L+ramda); Wc(1)=ramda/(L+ramda);%+1-alpha^2+belta; %%%权值 if i==1 xerror=rangeerror^2*cos(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*sin(Zmeasure(1,i))^2; yerror=rangeerror^2*sin(Zmeasure(1,i))^2+Zmeasure(2,i)^2*azimutherror^2*cos(Zmeasure(1,i))^2; xyerror=(rangeerror^2-Zmeasure(2,i)^2*azimutherror^2)*sin(Zmeasure(1,i))*cos(Zmeasure(1,i)); 信息融合大作业 ——维纳最速下降法滤波器，卡尔曼滤波器设计及Matlab仿真 1.滤波问题浅谈 估计器或滤波器这一术语通常用来称呼一个系统，设计这样的系统是为了从含有噪声的数据中提取人们感兴趣的，接近规定质量的信息。由于这样一个宽目标，估计理论应用于诸如通信、雷达、声纳、导航、地震学、生物医学工程、 金融工程等众多不同的领域。例如，考虑一个数字通信系统，其基本形式由发 射机、信道和接收机连接组成。发射机的作用是把数字源(例如计算机)产生的0、1符号序列组成的消息信号变换成为适合于信道上传送的波形。而由于符号间干扰和噪声的存在，信道输出端收到的信号是含有噪声的或失真的发送信号。接收机的作用是，操作接收信号并把原消息信号的一个可靠估值传递给系统输出端的某个用户。随着通信系统复杂度的提高，对原消息信号的还原成为通信系统中最为重要的环节，而噪声是接收端需要排除的最主要的干扰，人们也设计出了针对各种不同条件应用的滤波器，其中最速下降算法是一种古老的最优化技术，而卡尔曼滤波器随着应用条件的精简成为了普适性的高效滤波器。2．维纳最速下降算法滤波器 2.1 最速下降算法的基本思想 考虑一个代价函数，它是某个未知向量的连续可微分函数。函数 将的元素映射为实数。这里，我们要寻找一个最优解。使它满足如下条件 (2.1) 这也是无约束最优化的数学表示。 特别适合于自适应滤波的一类无约束最优化算法基于局部迭代下降的算法： 从某一初始猜想出发，产生一系列权向量，使得代价函数在算法的每一次迭代都是下降的，即 其中是权向量的过去值，而是其更新值。 我们希望算法最终收敛到最优值。迭代下降的一种简单形式是最速下降法，该方法是沿最速下降方向连续调整权向量。为方便起见，我们将梯度向量表示为基于扩展卡尔曼滤波算法的matlab程序
维纳最速下降法滤波器卡尔曼滤波器设计及Matlab仿真
相关主题
文本预览