实验报告-卡尔曼滤波

数字信号处理实验报告姓名：专业:通信与信息系统学号:日期：2015.11实验内容任务一：一连续平稳的随机信号tx，自相关函数txer，信号tx为加性噪声所干扰，噪声是白噪声，测量值的离散值kz为已知，sTs02.0，-3.2，-0.8，-14，-16，-17，-18，-3.3，-2.4，-18，-0.3，-0.4，-0.8，-19，-2.0，-1.2，-11，-14，-0.9，-0.8，10，0.2，0.5，-0.5，2.4，-0.5，0.5，-13，0.5，10，-12，0.5，-0.6，-15，-0.7，15，0.5，-0.7，-2.0，-19，-17，-11，-14，自编卡尔曼滤波递推程序，估计信号tx的波形。任务二：设计一维纳滤波器。（1）产生三组观测数据：首先根据nwnasns1产生信号ns，将其加噪（信噪比分别为20dB，10dB，6dB），得到观测数据nx1，nx2，nx3。（2）估计nxi，1i，2，3的AR模型参数。假设信号长度为L，AR模型阶数为N，分析实验结果，并讨论改变L，N对实验结果的影响。实验任务一1.卡尔曼滤波原理1.1卡尔曼滤波简介早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。它用状态空间法描述系统，由状态方程和量测方程所组成，即知道前一个状态的估计值和最近一个观测数据，采用递推的算法估计当前的状态值。由于卡尔曼滤波采用递推法，适合于计算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很多领域，并取得了很好的结果。卡尔曼滤波一经出现，就受到人们的很大重视，并在实践中不断丰富和完善，其中一个成功的应用是设计运载体的高精度组合导航系统。卡尔曼滤波具有以下的特点:（1）算法是递推的，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼算法适用于计算机处理。（2）用递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号可以是平稳的，也可以是非平稳的，即卡尔曼滤波适用于非平稳过程。（3）卡尔曼滤波采取的误差准则仍为估计误差的均方值最小。1.2卡尔曼滤波的状态方程和测量方程假设某系统k时刻的状态变量为kx，状态方程和量测方程（输出方程）表示为kkkkkkkkvxCywxAx111其中，kx是状态变量；1kw表示输入信号是白噪声；kv是观测噪声；ky是观测数据。为了推导简单，假设状态变量的增益矩阵A不随时间发生变化，kw，kv都是均值为零的正态白噪声，方差分别是kQ和kR，并且初始状态与kw，kv都不相关，表示相关系数。即：kjkvvkvkkkjkwwkwkkRRvEvQQwEwjkjk,2,2,,0:,,0:其中jkjkkj011.3卡尔曼滤波的递推算法卡尔曼滤波采用递推算法来实现，其基本思想是先不考虑输入信号kw和观测噪声kv的影响，得到状态变量和输出信号（即观测数据）的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。因此，卡尔曼滤波器的关键是计算出加权矩阵的最佳值。当不考虑观测噪声和输入信号时，状态方程和量测方程为1'1'kkkkkkkkkxACxCyxAx显然，由于不考虑观测噪声的影响，输出信号的估计值与实际值是有误差的，用ky~表示'~kkkyyy为了提高状态估计的质量，用输出信号的估计误差ky~来校正状态变量11'1CkkkkkkkkkkkkkxAyHxAyyHxAx其中，kH为增益矩阵，即加权矩阵。经过校正后的状态变量的估计误差及其均方值分别用kx~和kP表示，把未经校正的状态变量的估计误差的均方值用'kP表示TkkkkkTkkkkkkkkxxxxEPxxxxEPxxx''~卡尔曼滤波要求状态变量的估计误差的均方值kP为最小，因此卡尔曼滤波的关键即为通过选择合适的kH，使得kP取得最小值。首先推导状态变量的估计值kx和状态变量的的估计误差kx~，然后计算~x的均方值kP，通过化简kP，得到一组卡尔曼滤波的递推公式：'11'1''11kkkkkTkkkkkTkkkTkkkkkkkkkkkPCHIPQAPAPRCPCCPHxACyHxAx假设初始条件kA，kC，kQ，kR，ky，1-kx，1-kP已知，其中00xEx，00varxP，那么递推流程如下：1-kx，1-kP'kPkHkx，kP2.卡尔曼滤波递推程序编程思想题目分析（1）由于信号tx为加性噪声所干扰，可知kkkvxy，所以1kC（2）又因为噪声为白噪声，所以102vvvkSR（3）因为txer，所以1112111100111111zezeezeezzezezezezmrzSmmmmmmmmmmmmmmmxxxx由此可知，1111zezBz，即111nwnxenx，可得到：1eAk，因为抽样间隔sTs02.0，所以：02.0eeAsTk。（4）因此nxenxnwsT1，所以T因此04.01eQk编程分析由上面的分析可知初始条件kA，kC，kQ，kR，ky已知，在仿真中假设00x，则00x，10P，由以上参数可得卡尔曼实际递推公式'04.0104.0'1''102.0102.011kkkkkkkkkkkkkPHIPePePPPHxeyHxex将得到的公式代入前面分析的递推公式，即可进行迭代得到结果kx。3.MATLAB源代码根据以上分析，编写matlab程序如下：%%%---------------卡尔曼滤波-----------------%-----说明%X(k+1)=Ak*X(k)+W(k);%Y(k)=Ck*X(k)+V(k)%%clear;clc;%基本参数值Ak=exp(-0.02);Ck=1;Qk=1-exp(-0.04);Rk=1;%初始值设置X0=0;P0=1;%观测值y(k)Y=[-3.2-0.8-14-16-17-18-3.3-2.4-18-0.3-0.4-0.8-19-2.0-1.2...-11-14-0.90.8100.20.52.4-0.50.5-130.510-120.5-0.6-15-0.715...0.5-0.7-2.0-19-17-11-14];%数据长度N=length(Y);fork=1:Nifk==1%k=1时由初值开始计算P_(k)=Ak*P0*Ak'+Qk;H(k)=P_(k)*Ck'*inv(Ck*P_(k)*Ck'+Rk);X(k)=Ak*X0+H(k)*(Y(k)-Ck*Ak*X0);I=eye(size(H(k)));P(k)=(I-H(k)*Ck)*P_(k);else%k1时，开始递推%递推公式P_(k)=Ak*P(k-1)*Ak'+Qk;H(k)=P_(k)*Ck'*inv(Ck*P_(k)*Ck'+Rk);X(k)=Ak*X(k-1)+H(k)*(Y(k)-Ck*Ak*X(k-1));I=eye(size(H(k)));P(k)=(I-H(k)*Ck)*P_(k);endendM=1:N;T=0.02*M;%作图，画出x(t)的波形figure(1)plot(T,Y,'r','LineWidth',1);xlabel('t');ylabel('y(t)');title('卡尔曼滤波-测量信号y(t)波形');grid;figure(2)plot(T,X,'b','LineWidth',1);xlabel('t');ylabel('x(t)');title('卡尔曼滤波-估计信号x(t)波形');grid;4.实验结果00.10.20.30.40.50.60.70.80.9-20-15-10-5051015ty(t)卡尔曼滤波-测量信号y(t)波形00.10.20.30.40.50.60.70.80.9-12-10-8-6-4-202tx(t)卡尔曼滤波-估计信号x(t)波形实验任务二1.维纳滤波器原理维纳-霍夫方程krkhmkrmhkrxxmxxxd*0当nh是一个长度为M的因果序列（即一个长度为M的FIR滤波器）时，维纳-霍夫方程表述为，，，210*10kkrkhmkrmhkrxxMmxxxd定义00.10.20.30.40.50.60.70.80.9-20-15-10-5051015卡尔曼滤波测量信号z(t)估计信号x(t)02120111011021xxxxxxxxxxxxxxxxxxxxxdxdxdxdMrMrMrMrrrMrrrMrrrhhhRRh则可写成矩阵的形式，即hRRxxxd对上式求逆，得到RRhxdxx1由以上式子可知：若已知期望信号与观测数据的互相关函数及观测数据的自相关函数，则可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度较大时，计算工作量很大，并且需要计算xxR的逆矩阵，从而要求的存储量也很大预测是根据观测到对的过去数据来估计当前或将来的信号值。维纳预测是已知以前时刻的p个数据pnxnxnx，，，21，估计当前时刻n，或者未来Nn时刻的信号值，即估计0^NNns，，估计得到的结果仍然要求满足均方误差最小的准则。信号可以预测是由于信号内部存在着关联性。预测是利用数据前后的关联性，根据其中一部分推知其余部分。一步线性预测的时域解已知1nx，2nx，…，pnx，预测nx，假设噪声0nv，这样的预测成为一步线性预测。设定系统的单位脉冲响应为nh。根据现行系统的基本理论，输出信号pkknxkhnxny1令khapk，则pkpkknxanx1-预测误差pkpkpkpkknxaknxanxnxnxne01其中10pa212pkpkknxanxEneE要使均方误差为最小值，要求102laneEpl，2，...，p又因为，我们可以得到10*llnxneE，2，...，p所以01pkxxpkxxlkralr1l，2，...，p（1）由于预测器的输出nx是输入信号的线性组合，所以可得：0*nxneE以上说明误差信号与输入信号满足正交性原理，预测误差与预测信号值同样满足正交性原理。预测误差的最小均方值pkxxpkxxpkpkkrarnxknxanxEnxneEnxnxneEneE11****min20（2）由（1）（2）联立方程组