仿真原理Word文档格式.docx

1、 x（n - 3） x（n - 4） 状态变量增益矩阵为：1.352 -1.338 0.662 0.240 1000 A= 0100 0010 仅对 x 方向进行估计，即考虑一维的情况下，有： 状态变量与输出信号之间的增益矩阵 Ck =1 0 0 0;量测噪声协方差阵为 Rk =1;1 0 0 0 0 1 0 0 噪声的均方误差阵为 Pk = ； 0 0 1 0 0 0 0 1 3、仿真要求分别利用 Wiener 滤波器和 Kalman 滤波器通过测量信号估计 x（n） 的波形4、仿真原理（1） 维纳滤波原理维纳滤波的原理是根据全部过去的和当前的观测数据 x（n）,x（n-1）,来估计信号的当

2、前值。以均方误差最小条件下求解系统的传递函数 H（z）或单位冲激响应 h（n）。维纳滤波的信号模型是从信号与噪声的相关函数得到。Wiener 滤波器的一般结构如下：-线性滤波器y（n） = d（n）x（n）e（n）d （n）有一期望响应 d（n）,滤波器系数的设计准则是使得滤波器的输出 y（n）是均方意义上对期望响应的最优线性估计。线性系统输出为：y（n） = s（n）= h（m）x（n - m）m2 均方误差为： E e2 （n） = E s（n）- h（m）x（n - m） m=0 E 维纳滤波器的设计，实际上就是在最小均方误差条件下，即hj=0 ，确定滤波器的冲激响应 h（n）或系统函数

3、 H（z），可等效于求解维纳-霍夫方程：h=xx-1xs。（2） 卡尔曼滤波原理不需要全部过去的观察数据,只根据前一个估计值和最近一个观察数据来估计信号的当前值。它是用状态空间法描述系统，即由状态方程和量测方程组成。解是以估计值（是状态变量的估计值）的形式给出的，其算法是递推。卡尔曼滤波的信号模型则是从状态方程和量测方程得到。卡尔曼滤波的信号模型（一维）如下：A（z）wkskukyk =sk +uk离散系统的 n 维状态方程： x（k ） = Ak xk -1 + wk -1离散系统的 m 维量测方程： ykck xkuk令 Ak 表示状态变量之间的增益矩阵， Ck 为状态变量与输出信号之间的

4、增益矩阵，不随时间发生变化，动态噪声wk 与观测噪声uk 都是零均值的正态噪声，kkk k且两者互不相关， R var u = E uuT 为量测噪声协方差矩阵，中Qvarw =E w wT 为动态噪声协方差矩阵。kkk k系统初始条件为：E x0 = m00000varx = E （x - m）（x - m）T = pkcovx , w = E x wT = 0covx ,u = E xuT 卡尔曼滤波的基本思想是先不考虑激励噪声wk 和观测噪声uk ，得到状态的估计值 x 和观测数据的估计值 y ,再用观测数据的估计误差 y% =y - y去修正状kkkkk态的估计值 xk ，通过选择修正

5、矩阵 H 使得状态估计误差的均方值 Pk 最小。卡尔曼滤波的递推公式如下：xk= Ak xk -1 + Hk （ yk - Ck Ak xk -1 ）-1kk kk k kk H = PCT （C PCT + R ）P = A P AT + Qkk k -1 kk -1P = （I - H C ）Pkk kk假设初始条件 Ak，Ck ，Q，k ，R，kykxk -1Pk -1 已知，其中x0 = E x0 ，P0 = varx0 ，则卡尔曼滤波的递推流程如下：x，PPk -1k -1P =A P AT +Qkkk k -1kk -1H = PCT + Rxk = Ak xk -1 + Hk （

6、 yk - Ck Ak xk -1 ）Hk）-1PkP= （I - H C）Pkk kk二、仿真流程（1） 维纳滤波仿真流程如下：初始条件+状态方程观测信号与期望信号的互相关函数 yx观测样本 y=x+v，v 表示加性高斯白噪声真实值 x输出滤波估计值 xk_s维纳滤波器系数：h=yy-1yx观测信号的自相关函数 yy误差=真实值-滤波估计值_（2） 卡尔曼滤波仿真流程如下：初始值x0k=k+1计算估计值：观测值 yk yk=x+vk kk）（P = I - H CP计算滤波后的均方误差阵：kk kk k k TH = P CC P C + R计算增益矩阵：P =A P AT +Q计算未考虑噪

7、声时的均方误差阵：给定 Ak、Ck、Rk、Pk 的初始值三、仿真结果及分析1、利用维纳滤波器进行估计：（1） Wiener 滤波器的阶数取 101 阶，观测点数取 100，u（n） 的方差为 1 时，通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图 1 所示： 估 估 估 估 估 估估642-2-4-60102030405060708090100图 1 采用维纳滤波，噪声方差为 1 时，真实轨迹、观测样本及估计轨迹的比较从图 1 可以看出，利用 Wiener 滤波器通过测量信号估计 x（n） 的波形，得到的估计轨迹接近于真实轨迹。当噪声方差为 1 时，估计轨迹与真实轨迹间的误差很小。（

8、2） Wiener 滤波器的阶数取 101 阶，观测点数取 100，u（n） 的方差为 1时，通过仿真得到平均误差如下图 2 所示：1.50.5-0.5-1.5图 2 采用维纳滤波，噪声方差为 1 时的平均误差从图 2 可以看出，利用 Wiener 滤波器通过测量信号估计 x（n） 的波形，当噪声方差为 1 时，得到的估计轨迹与真实轨迹的平均误差较小，而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。（3） Wiener 滤波器的阶数取 101 阶，观测点数取 100，u（n） 的方差为 4 时，通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图 3 所示：15

9、105-5-10-15图 3 采用维纳滤波，噪声方差为 4 时，真实轨迹、观测样本及估计轨迹的比较从图 3 可以看出，当噪声方差为 4 时，估计轨迹与真实轨迹间的误差变大了。利用 Wiener 滤波器通过测量信号估计 x（n） 的波形，得到的估计轨迹接近于真实轨迹。（4） Wiener 滤波器的阶数取 101 阶，观测点数取 100，u（n） 的方差为 4时，通过仿真得到平均误差如下图 4 所示：3-3图 4 采用维纳滤波，噪声方差为 4 时的平均误差从图 4 可以看出，当噪声方差为 4 时，估计轨迹与真实轨迹间的平均误差变大了。而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐

10、减小的趋势。2、利用卡尔曼滤波器进行估计：（1） 采用卡尔曼滤波，观测点数取 100，u（n） 的方差为 1 时，通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图 5 所示：图 5采用卡尔曼滤波，噪声方差为 1 时，真实轨迹、观测样本及估计轨迹的比较从图 5 可以看出，利用卡尔曼滤波器，通过测量信号估计 x（n） 的波形，得到的估计轨迹接近于真实轨迹。（2） 采用卡尔曼滤波，观测点数取 100，u（n） 的方差为 1 时，通过仿真得到平均误差如下图 6 所示：2.5-2.5图 6 采用卡尔曼滤波，噪声方差为 1 时的平均误差从图 6 可以看出，利用卡尔曼滤波器通过测量信号估计 x（n）

11、 的波形，当噪声方差为 1 时，得到的估计轨迹与真实轨迹的平均误差较小。（3） 采用卡尔曼滤波，观测点数取 100，u（n） 的方差为 4 时，通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图 7 所示：-20图 7 采用卡尔曼滤波，噪声方差为 4 时，真实轨迹、观测样本及估计轨迹的比较从图 7 可以看出，当噪声方差为 4 时，估计轨迹与真实轨迹间的误差变大了。利用卡尔曼滤波器通过测量信号估计 x（n） 的波形，得到的估计轨迹接近（4） 采用卡尔曼滤波，观测点数取 100，u（n） 的方差为 4 时，通过仿真得到平均误差如下图 8 所示：8-8图 8 采用卡尔曼滤波，噪声方差为 4 时

12、的平均误差从图 8 可以看出，利用卡尔曼滤波器通过测量信号估计 x（n） 的波形，当噪声方差为 4 时，估计轨迹与真实轨迹间的平均误差变大了。四、结论通过以上的仿真结果及其分析可以得到如下结论：（1） 利用 Wiener 滤波器和 Kalman 滤波器都可以通过测量信号对 x（n） 的波形进行估计，估计轨迹接近于真实轨迹；（2） 随着噪声方差的变大，通过维纳滤波和卡尔曼滤波得到的估计轨迹与真实轨迹间的误差也随之变大；（3） 相比较而言，利用维纳滤波器进行估计的误差稍小于通过卡尔曼滤波器进行估计的误差；（4） 在利用维纳滤波器和卡尔曼滤波器进行估计时，刚开始的时候平均误差相对较大，随着时间的推移

13、，平均误差有逐渐减小的趋势。五、附件：1、维纳滤波仿真源程序（噪声方差取 1）clc;clear all; maxlag=100;%N=100;%观测点数取 100 x=zeros（N,1）;y=zeros（N,1）;var=1;%列出状态方程x（1）=randn（1,1）;%令 x（-1）=x（-2）=x（-3）=x（-4）=0 x（2）=randn（1,1）+1.352*x（1）; x（3）=randn（1,1）+1.352*x（2）-1.338*x（1）;x（4）=randn（1,1）+1.352*x（3）-1.338*x（2）+0.602*x（1）;for n=5:Nx（n）=1.35

14、2*x（n-1）-1.338*x（n-2）+0.602*x（n-3）-0.24*x（n-4）+randn（1,1）; %x 为真实值end;v=randn（N,1）;y=x+v;%z_x 为观测样本值=真值+噪声%滤波x = x;y = yxk_s（1）=y（1）;%赋初值xk_s（2）=y（2）;xk_s（3）=y（3）;xk_s（4）=y（4）; xk=y（1）;y（2）;y（3）;y（4）;%维纳滤波器的生成rx,lags=xcorr（y,maxlag,biased）;%观测信号的自相关函数rx1=toeplitz（rx（101:end）;%对称化自相关函数矩阵使之成为方阵，滤波器的阶数

15、为 101 阶rx2=xcorr（x,y,maxlag,%观测信号与期望信号的互相关函数rx2=rx2（101:end）;h=inv（rx1）*rx2%维纳-霍夫方程xk_s=filter（h,1,y）;%加噪信号通过滤波器后的输出%e_x=0; eq_x=0; e_x1=N:1;%计算滤波的均值,计算滤波误差的均值for i=1:e_x（i）=x（i）-xk_s（i）; %误差=真实值-滤波估计值end%作图t=1:N;figure（1）;plot（t,x,r-,t,y,g:,t,xk_s,b-.legend（真实轨迹,观测样本估计轨迹 figure（2）;plot（e_x）;平均误差2、维

16、纳滤波仿真源程序（噪声方差取 4）var=2;v=4*randn（N,1）;3、卡尔曼滤波仿真源程序（噪声方差取 1）clear;I=1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;%I 为四阶单位阵Ak=1.352,-1.338,0.662,0.240;1,0,0,0;0,0,1,0;%状态变量之间的增益矩阵 AkCk=1 0 0 0;%状态变量与输出信号之间的增益矩阵 Ck,（一维:仅仅对 x 方向进行估计）Rk=1;%量测噪声协方差阵Pk=1 0 0 00 1 0 00 0 1 00 0 0 1 ;%噪声的均方误差阵xk=y（1）;Qk = 1;%Kalman 滤波开始,估计循环for r=5:yk=y（r）;Pk1=Ak*Pk*Ak+Qk;%（未考虑噪声）k 时刻滤波的均方误差矩阵Hk=Pk1*Ck*inv（Ck*Pk1*Ck+Rk）; %增益方程xk=Ak*xk+Hk*（yk-Ck*Ak*xk）;%递推公式Pk=（I-Hk*Ck）*Pk1;%滤波后的均方误差矩阵xk_s（r）=xk（1,1）;%xk_s 为估计值%end e_x=0;4、卡尔曼滤波仿真源程序（噪声方差取 4） cle