卡尔曼滤波器matlab代码

1、信息交融大作业维纳最速降落法滤波器，卡尔曼滤波器设计及Matlab仿真时间：2010-12-6专业：信息工程班级：09030702学号：姓名：马志强滤波问题浅谈预计器或滤波器这一术语往常用来称号一个系统，设计这样的系统是为了从含有噪声的数据中提取人们感兴趣的，靠近规定质量的信息。因为这样一个宽目标，预计理论应用于诸如通讯、雷达、声纳、导航、地震学、生物医学工程、金融工程等众多不一样的领域。比如，考虑一个数字通讯系统，其基本形式由发射机、信道和接收机连结构成。发射机的作用是把数字源(比如计算机)产生的0、符号序列构成的信息信号变换成为合适于信道上传递的波形。而因为符号间扰乱和噪声的存在，信道输出

2、端收到的信号是含有噪声的或失真的发送信号。接收机的作用是，操作接收信号并把原信息信号的一个靠谱估值传达给系统输出端的某个用户。跟着通讯系统复杂度的提升，对原信息信号的复原成为通讯系统中最为重要的环节，而噪声是接收端需要清除的最主要的扰乱，人们也设计出了针对各样不一样条件应用的滤波器，此中最速降落算法是一种古老的最优化技术，而卡尔曼滤波器跟着应用条件的精简成为了普适性的高效滤波器。2维纳最速降落算法滤波器2.1最速降落算法的基本思想考虑一个代价函数，它是某个未知向量的连续可微分函数。函数将的元素映照为实数。这里，我们要找寻一个最优解。使它知足以下条件(2.1)这也是无拘束最优化的数学表示。特别合

3、适于自适应滤波的一类无拘束最优化算法鉴于局部迭代降落的算法：从某一初始猜想出发，产生一系列权向量，使得代价函数在算法的每一次迭代都是降落的，即此中是权向量的过去值，而是其更新值。我们希望算法最后收敛到最优值。迭代降落的一种简单形式是最速降落法，该方法是沿最速降落方向连续调整权向量。为方便起见，我们将梯度向量表示为(2.2)所以，最速降落法能够表示为(2.3)此中代表进度，是正常数，称为步长参数，1/2因子的引入是为了数学上处理方便。在从到的迭代中，权向量的调整量为(2.4)为了证明最速降落算法知足式(2.1)，在处进行一阶泰勒睁开，获取(2.5)此式关于较小时是建立的。在式(2.4)中设为负值

4、向量，因此梯度向量负值向量，所以使用埃尔米特转置。将式(2.4)用到式(2.5)中，获取也为此式表示当为正数时，数减小，当时，代价函数趋于最小值。所以，跟着的增添，代价函。2.2最速降落算法应用于维纳滤波器考虑一个横向滤波器，其抽头输入为，对应的抽头权值为。抽头输入是来自零均值、有关矩阵为的广义安稳随机过程的抽样值。除了这些输入外，滤波器还要一个希望响应，以便为最优滤波供给一个参照。在时刻抽头输入向量表示为，滤波器输出端希望响应的预计值为,此中是由抽头输所张成的空间。空过比较希望响应及其预计值，能够获取一个预计偏差，即(2.6)这里是抽头权向量与抽头输入向量的内积。能够进一步表示为相同，抽头输

5、入向量可表示为假如抽头输入向量在时刻的代价函数和希望响应是结合安稳的，此时均方偏差或许是抽头权向量的二次函数，于是能够获取(2.7)此中，为目标函数的方差，抽头输入向量与希望响应关向量，及为抽头输入向量的有关矩阵。进而梯度向量能够写为的相互(2.8)此中在列向量中部和虚部矩阵和互有关向量和分别是代价函数对应第个抽头权值的偏导数。对最速降落算法应用而言，假定式(2.8)已知，则关于给定的抽头权向量为的实中有关(2.9)它描绘了为那滤波中最速降落法的数学表达式。3.卡尔曼滤波器3.1卡尔曼滤波器的基本思想卡尔曼滤波器是用状态空间观点描绘其数学公式的，此外新奇的特色是，他的解递归运算，能够不加改正地

6、应用于安稳和非安稳环境。特别是，其状态的每一次更新预计都由前一次预计和新的输入数据计算获取，所以只要储存前一次预计。除了不需要储存过去的全部观察数据外，卡尔曼滤波计算比直接根据滤波过程中每一步全部过去数据进行估值的方法都更为有效。+图3.1线性动向失散时间系统的信号流图表示“状态”的观点是这类表示的基础。状态向量，简单地说状态，定义为数据的最小会合，这组数据足以独一地描绘系统的自然动向行为。换句话说，状态由展望系统将来特征时所素要的，与系统的过去行为有关的最少的数据构成。典型地，比较有代表性的状况是，状态是未知的。为了预计它，我们使用一组观测数据，在途顶用向量表示。成为观察向量或许简称观察值，

7、并假定它是维的。在数学上，图3.1表示的信号流图隐含着一下两个方程：过程方程(3.1)式中，向量表示噪声过程，可建模为零均值的白噪声过程，且其有关矩阵定义为丈量方程(3.2)此中是已知的丈量矩阵。向量称为丈量噪声，建模为零均值的白噪声过程，其有关矩阵为(3.3)丈量方程(3.2)确定了可测系统输出与状态之间的关系，如图3.1所示。3.2新息过程为了求解卡尔曼滤波问题，我们将应用鉴于新息过程的方法。依据以前所述，用向量表示时刻到时刻全部观察数据过去值给定的状况下，你时刻观察数据的最小均方预计。过去的值用观察值表示，他们张成的向量空间用表示。进而能够定义新息过程以下：(3.4)此中向量表示观察数据

8、的新息。3.3应用新息过程进行状态预计下边，我们依据信息过程导出状态的最小均方预计。依据推导，这个估计能够表示成为新息过程序列的线性组合，即(3.5)此中是一组待定的矩阵。依据正交性原理，展望状态偏差向量与新息过程正交，即(3.6)将式(3.5)代入式(3.6)，并利用新息过程的正交性质，即得(3.7)所以，式(3.7)两边同时右乘逆矩阵，可得的表达式为(3.8)最后，将式(3.8)带入式(3.5)，可得最小军方差预计(3.9)故关于，有(3.10)但是，时刻的状态与时刻的状态的关系式由式能够推导出对于，有(3.11)此中只与观察数占有关。所以可知，与相互正交(此中)。利用式(3.11)以及当

9、时的计算公式，可将式(3.10)右侧的乞降项改写为(3.12)为了进一步议论，引入以下基本定义。3.4卡尔曼增益定义矩阵(3.13)此中是状态向量和新息过程的互有关矩阵。利用这必定义和式(3.12)的结果，能够将式(3.10)简单重写为(3.14)式(3.14)拥有明确的物理意义。它注明：线性动向系统状态的最小均方预计能够由前一个预计求得。为了表示对卡尔曼创始性贡献的认同，将矩阵称为卡尔曼增益。此刻剩下独一要解决的问题是，如何以一种便于计算的形式来表示卡尔曼增益。为此，第一将与乘积的希望表示为(3.15)式中利用了状态与噪声向量互不有关这一事实。其次，因为展望状态偏差向量与预计正交，所以与趁机

10、的希望为零。这样，用展望状态偏差向量取代相乘因子，将不会引起式(3.15)变化，故有(3.16)由此，可将上式进一步变化为(3.17)此刻我们从头定义卡尔曼增益。为此，将式(3.17)代入式(3.13)得(3.18)此刻我们已经认识了卡尔曼滤波的整个过程和相应的参数设置，为了能够更为方便利用计算机仿真切现，特将此中参数变量进行小结。卡尔曼变量和参数小结变量定义维数时刻状态时刻状态值从时刻到时刻的转移矩阵时刻的丈量矩阵过程噪声的有关矩阵过程噪声的有关矩阵给定观察值在时刻状态的展望预计给定观察值在时刻状态的滤波预计时刻卡尔曼增益矩阵时刻新息向量新息向量的有关矩阵中偏差有关矩阵中偏差有关矩阵鉴于单步

11、展望的卡尔曼滤波器的小结观察值=转移矩阵=丈量矩阵=过程噪声的有关矩阵=丈量噪声的有关矩阵=4Matlab仿真为了简化，这里只议论简单的一维单输入单输出线性系统模型，此中加入白噪声作为系统的扰动，详细仿真结果能够获取以下4.1维纳最速降落法滤波器仿真结果以上为最速降落法中不一样的递归步长所致使的追踪成效变化，关于最速降落法中的步长是影响其算法稳固的重点，最速降落算法稳固的充分必需条件是条件步长因子为小于输入自有关矩阵的最大特色值倒数的2倍。上边的序列分别从有关矩阵的随大特色之2倍的0.4倍开始变化至其1倍，最后一幅图象能够看出其已经不再收敛，下边是大于输入有关矩阵的最大特色值2倍步长时所表现的

12、追踪结果能够看出其已经显然发散，不再是我们所希望的滤波算法。所以能够总结出，对于最速降落法来说，步长的选用是很重要的，依据不一样条件的需求，选用正确的步长，能为算法的迅速高效供给基础。4.2卡尔曼滤波器仿真结果从图中能够发现，卡尔曼滤波器能够特别有效地在比较大的扰乱下比较正确地反应真切值，假如观察端加入扰乱较大时，卡尔曼滤波器能够较为有效地进行滤波，可是当状态端的扰乱增大时，卡尔曼滤波器的滤波成效也会随之降落。以下列图，是加大了状态端的扰乱，所体现的滤波成效。如上图所示，状态端的扰乱致使状态不稳固，卡尔曼滤波器的预计值也出现了比较大的颠簸。假如将状态端的扰乱再增大，则会出现更为严重的滤波考验，

13、滤波成效以下。这是的状态已经很牵强了，所以，研究更为有效的多方法卡尔曼滤波器也显得十分必需了。4.3一种不需初始化的卡尔曼滤波器仿真这类滤波器不过实现了无需对部分变量进行初始化的设计，没有特别意义上的改良经典卡尔曼滤波器自己性能的特色。仿真图以下4.4后联光滑滤波的卡尔曼滤波器仿真不过在经典卡尔曼滤波器后端联接了光滑滤波器，对性能改良的成效其实不特别显然，仿真图以下如图中所表示，即便光滑过的估值与观察值之间的差异也不是特别令人满意，所以，关于经典卡尔曼滤波的研究还需要更深一步进行，因为时间和能力有限，本次的作业关于卡尔曼及其余滤波器的研究只好达到这类程度，希望在此后的学习中，能发现更好的对经典

14、卡尔曼滤波器的改良方法。5Matlab源代码(部分参照自互联网)5.1经典卡尔曼滤波器clearN=200;w(1)=0;x(1)=5;a=1;c=1;Q1=randn(1,N)*1;%过程噪声Q2=randn(1,N);%丈量噪声fork=2:N;x(k)=a*x(k-1)+Q1(k-1);end%状态矩阵fork=1:N;Y(k)=c*x(k)+Q2(k);endp(1)=10;s(1)=1;fort=2:N;Rww=cov(Q1(1:t);Rvv=cov(Q2(1:t);p1(t)=a.2*p(t-1)+Rww;b(t)=c*p1(t)/(c.2*p1(t)+Rvv);%kalman增益

15、s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1);p(t)=p1(t)-c*b(t)*p1(t);endt=1:N;plot(t,s,r,t,Y,g,t,x,b);%红色卡尔曼，绿色观察值，蓝色状态值legend(kalmanestimate,ovservations,truth);5.2最速降落法clcclearallN=30;q=2.1;%q1&q2/Ryx最大特色值hn=zeros(1,N);hn(:)=5;vg=0;Rxx=xcorr(1);Ryx=min(min(corrcoef(1,1+randn);echoofffori=1:N-1;%vg=2*Rxx*hn

16、(:,i)-2*Ryx;%hn(:,i+1)=hn(:,i)-1/2*q*vg;vg=2*Rxx*hn(i)-2*Ryx;hn(i+1)=hn(i)-1/2*q*vg;m(i)=1;endt=1:N-1;plot(t,hn(t),r-,t,m(t),b-);5.3后联光滑滤波器的卡尔曼滤波器clearclc;N=300;CON=5;x=zeros(1,N);x(1)=1;p=10;Q=randn(1,N)*0.2;%过程噪声协方差R=randn(1,N);%观察噪声协方差y=R+CON;%加过程噪声的状态输出fork=2:NQ1=cov(Q(1:k-1);%过程噪声协方差Q2=cov(R(1:

17、k-1);x(k)=x(k-1);%预预计k时刻状态变量的值p=p+Q1;%对应于预估值的协方差kg=p/(p+Q2);%kalmangainx(k)=x(k)+kg*(y(k)-x(k);p=(1-kg)*p;endFilter_Wid=10;smooth_res=zeros(1,N);kalman_p=zeros(1,N);fori=Filter_Wid+1:Ntempsum=0;kalman_m=0;forj=i-Filter_Wid:i-1tempsum=tempsum+y(j);kalman_m=kalman_m+x(j);endkalman_p(i)=kalman_m/Filter_Wid;smooth_res(i)=tempsum/Filter_Wid;%光滑滤波endfigure(1);hist(