数字信号处理实验报告姓名：任伟平 专业: 通信与信息系统 学号: 2015111806 日期：2015.11实验内容任务一：一连续平稳的随机信号()t x ，自相关函数()tx er -=τ，信号()t x 为加性噪声所干扰，噪声是白噪声，测量值的离散值()k z 为已知，s T s 02.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，自编卡尔曼滤波递推程序，估计信号()t x 的波形。

任务二：设计一维纳滤波器。

（1）产生三组观测数据：首先根据()()()n w n as n s +-=1产生信号()n s ，将其加噪（信噪比分别为20dB ，10dB ，6dB ），得到观测数据() n x 1，() n x 2，() n x 3。

（2）估计() n x i ， 1=i ，2，3的AR 模型参数。

假设信号长度为L ，AR 模型阶数为N ，分析实验结果，并讨论改变L ，N 对实验结果的影响。

实验任务一 1. 卡尔曼滤波原理1.1 卡尔曼滤波简介早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推最优估计理论。

它用状态空间法描述系统，由状态方程和量测方程所组成，即知道前一个状态的估计值和最近一个观测数据，采用递推的算法估计当前的状态值。

由于卡尔曼滤波采用递推法，适合于计算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很多领域，并取得了很好的结果。

卡尔曼滤波一经出现，就受到人们的很大重视，并 在实践中不断丰富和完善，其中一个成功的应用是设计运载体的高精度组合导航系统。

卡尔曼滤波具有以下的特点:（1）算法是递推的，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼算法适用于计算机处理。

（2）用递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号可以是平稳的，也可以是非平稳的，即卡尔曼滤波适用于非平稳过程。

（3）卡尔曼滤波采取的误差准则仍为估计误差的均方值最小。

1.2 卡尔曼滤波的状态方程和测量方程假设某系统k 时刻的状态变量为k x ，状态方程和量测方程（输出方程）表示为kk k k k k k k v x C y w x A x +=+=---111其中，k x 是状态变量；1-k w 表示输入信号是白噪声；k v 是观测噪声；k y 是观测数据。

为了推导简单，假设状态变量的增益矩阵A 不随时间发生变化，k w ，k v 都是均值为零的正态白噪声，方差分别是k Q 和k R ，并且初始状态与k w ，k v 都不相关，γ表示相关系数。

即：[][]kjk v v k vk k kjk w w k w k k R R v E v Q Q w E w j k j k δγσδγσ======,2,2,,0:,,0:其中⎩⎨⎧≠==jk jk kj 01 δ1.3 卡尔曼滤波的递推算法卡尔曼滤波采用递推算法来实现，其基本思想是先不考虑输入信号k w 和观测噪声k v 的影响，得到状态变量和输出信号（即观测数据）的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。

因此，卡尔曼滤波器的关键是计算出加权矩阵的最佳值。

当不考虑观测噪声和输入信号时，状态方程和量测方程为1'1'-∧∧∧-∧∧===k k k k k k k k k x A C x C y x A x显然，由于不考虑观测噪声的影响，输出信号的估计值与实际值是有误差的，用k y ~表示'~k k k y y y ∧-=为了提高状态估计的质量，用输出信号的估计误差k y ~来校正状态变量⎪⎭⎫ ⎝⎛-+=⎪⎪⎭⎫ ⎝⎛-+=-∧-∧∧-∧∧11'1C k k k k k k k k kk k k k x A y H x A y y H x A x 其中，k H 为增益矩阵，即加权矩阵。

经过校正后的状态变量的估计误差及其均方值分别用k x ~和k P 表示，把未经校正的状态变量的估计误差的均方值用'k P 表示⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-=-=∧∧∧∧∧Tk k k k k T k k k k k kk k x x x x E P x x x x E P x x x ''~卡尔曼滤波要求状态变量的估计误差的均方值k P 为最小，因此卡尔曼滤波的关键即为通过选择合适的k H ，使得k P 取得最小值。

首先推导状态变量的估计值k x ∧和状态变量的的估计误差k x ~，然后计算~x 的均方值k P ，通过化简k P ，得到一组卡尔曼滤波的递推公式：()()⎪⎪⎪⎩⎪⎪⎪⎨⎧-=+=+=⎪⎭⎫ ⎝⎛-+=----∧-∧∧'11'1''11kk k k k T k k k kk Tk k k T k k k k k k k k k k k P C H I P QA P A P R C P C C P H x A C y H x A x 假设初始条件k A ，k C ，k Q ，k R ，k y ，1-k x ∧，1-k P 已知，其中[]00x E x =∧，[]00var x P =，那么递推流程如下：1-k x ∧，1-k P 'kPk H k x ∧，k P2. 卡尔曼滤波递推程序编程思想题目分析（1）由于信号()t x 为加性噪声所干扰，可知k k k v x y +=，所以1=k C（2）又因为噪声为白噪声，所以()102===vv v k S R σ（3）因为()tx er -=τ，所以()()()()()111211110111111--------∞==--=-∞=-∞=-∞=--∞=-∞=----=-+-=+===∑∑∑∑z e z e e z e ez z e z e ze zezm r z S m m mm m m mm m m mmm m mxx xx由此可知，1111----=z e z B z ，即()()()111-=---n w n x e n x ，可得到：1-=e A k，因为抽样间隔s T s 02.0=，所以：02.0--==e eA sT k 。

（4）因此()()()n x en x n w sT --+=1，所以[]()()[]T ww w k e n w n w E r Q 2210--====σ因此04.01--=e Q k编程分析由上面的分析可知初始条件k A ，k C ，k Q ，k R ，k y 已知，在仿真中假设00=x ，则00=∧x，10=P ，由以上参数可得卡尔曼实际递推公式()()⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-+=+=⎪⎭⎫ ⎝⎛-+=-----∧--∧-∧'04.0104.0'1''102.0102.011kk k k kkk k k k k k k P H I P e P e P P P H x e y H x e x将得到的公式代入前面分析的递推公式，即可进行迭代得到结果k x 。

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.9 0.8 10 0.2 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]; %数据长度N=length(Y);for k=1:Nif k==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%k>1时，开始递推%递推公式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. 实验结果0.10.20.30.40.50.60.70.80.9-20-15-10-5051015ty (t )卡尔曼滤波-测量信号y(t)波形0.10.20.30.40.50.60.70.80.9tx (t )卡尔曼滤波-估计信号x(t)波形实验任务二 1. 维纳滤波器原理维纳-霍夫方程()()()()()k r k h m k r m h k r xx m xx xd *0=-=∑+∞=当()n h 是一个长度为M 的因果序列（即一个长度为M 的FIR 滤波器）时，维纳-霍夫方程表述为()()()()() ，，，210*1==-=∑-=k k r k h m k r m h k r xx M m xx xd定义0.10.20.30.40.50.60.70.80.9-20-15-10-551015()()()()()()()()()()()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=02120111011021xx xx xx xx xx xx xx xx xx xx xd xd xd xdM r M r M r M r r r M r r r M r r r h h h R R h则可写成矩阵的形式，即h R Rxx xd=对上式求逆，得到R R hxd xx 1-=由以上式子可知：若已知期望信号与观测数据的互相关函数及观测数据的自相关函数，则可以通过矩阵求逆运算，得到维纳滤波器的最佳解。