仿真原理Word文档格式.docx
《仿真原理Word文档格式.docx》由会员分享，可在线阅读，更多相关《仿真原理Word文档格式.docx（22页珍藏版）》请在冰点文库上搜索。
x（n-3）÷
è
ø
è
ø
x（n-4）÷
状态变量增益矩阵为：
1.352-1.3380.6620.240ö
1 0 0 0÷
A=ç
0 1 0 0÷
0 0 1 0÷
仅对x方向进行估计，即考虑一维的情况下，有：
状态变量与输出信号之间的增益矩阵Ck=[1000];
量测噪声协方差阵为Rk=[1];
1000ö
0100÷
噪声的均方误差阵为Pk=ç
；
0010÷
0001÷
3、仿真要求
分别利用Wiener滤波器和Kalman滤波器通过测量信号估计x（n）的波形
4、仿真原理
（1）维纳滤波原理
维纳滤波的原理是根据全部过去的和当前的观测数据x（n）,x（n-1）, …来估计信号的当前值。
以均方误差最小条件下求解系统的传递函数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）
m
ì
ï
é
¥
ù
2ü
均方误差为：
Eé
ë
e2（n）ù
û
=Eí
ê
s（n）-å
h（m）x（n-m）ú
ý
î
m=0
û
ï
þ
¶
Eé
维纳滤波器的设计，实际上就是在最小均方误差条件下，即
¶
hj
=0，
确定滤波器的冲激响应h（n）或系统函数H（z），可等效于求解维纳-霍夫方程：
[h]=[φxx]-1[φxs]。
（2）卡尔曼滤波原理
不需要全部过去的观察数据,只根据前一个估计值和最近一个观察数据来估计信号的当前值。
它是用状态空间法描述系统，即由状态方程和量测方程组成。
解是以估计值（是状态变量的估计值）的形式给出的，其算法是递推。
卡尔曼滤波的信号模型则是从状态方程和量测方程得到。
卡尔曼滤波的信号模型（一维）如下：
A（z）
wk sk
uk yk=sk+uk
离散系统的n维状态方程：
x（k）=Akxk-1+wk-1
离散系统的m维量测方程：
yk＝c＋kxk uk
令Ak表示状态变量之间的增益矩阵，Ck为状态变量与输出信号之间的增益矩阵，不随时间发生变化，动态噪声wk与观测噪声uk都是零均值的正态噪声，
k k kk
且两者互不相关，R＝var[u]=Eé
uuTù
为量测噪声协方差矩阵，
中Q＝ var[w]=Eé
wwTù
为动态噪声协方差矩阵。
k k kk
系统初始条件为：
E[x0]=m0
0 0 0 0
var[x]=Eé
（x-m）（x-m）Tù
=p
k
cov[x,w]=Eé
xwTù
=0
cov[x,u]=Eé
xuTù
卡尔曼滤波的基本思想是先不考虑激励噪声wk和观测噪声uk，得到状态的
估计值xˆ'
和观测数据的估计值yˆ'
再用观测数据的估计误差y%=y-yˆ
去修正状
k k k k k
态的估计值xˆk，通过选择修正矩阵H使得状态估计误差的均方值Pk最小。
卡尔曼滤波的递推公式如下：
xˆk
=Akxˆk-1+Hk（yk-CkAkxˆk-1）
-1
k kk kkk k
H=P'
CT（CP'
CT+R）
í
P'
=APAT+Q
k kk-1k k-1
P=（I-HC）P'
k kk k
假设初始条件Ak，，Ck，Q，k，R，k
yk xˆk-1
Pk-1已知，其中
xˆ0=E[x0]，P0=var[x0]，则卡尔曼滤波的递推流程如下：
xˆ ，P P'
k-1
k-1
P'
=APAT+Q k
k kk-1
k k-1
H=P'
CT+R
xˆk=Akxˆk-1+Hk（yk-CkAkxˆk-1）
Hk
）-1
Pk P
=（I-HC
）P'
k kk k
二、 仿真流程
（1）维纳滤波仿真流程如下：
初始条件+状态方程
观测信号与期望信号的互相关函数φyx
观测样本y=x+v，
v表示加性高斯白噪声
真实值x
输出滤波估计值xk_s
维纳滤波器系数：
[h]=[φyy]-1[φyx]
观测信号的自相关函数φyy
误差=真实值-滤波估计值
_
（2）卡尔曼滤波仿真流程如下：
初始值xˆ0
k=k+1
计算估计值：
观测值ykyk=x+v
kk k
'
）
（
P=I-HC P
计算滤波后的均方误差阵：
k kk kkk
T
H=PC CPC+R
计算增益矩阵：
P'
=APAT+Q
计算未考虑噪声时的均方误差阵：
给定Ak、Ck、Rk、Pk的初始值
三、仿真结果及分析
1、利用维纳滤波器进行估计：
（1）Wiener滤波器的阶数取101阶，观测点数取100，u（n）的方差为1时，
通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图1所示：
□估估估
估估估估
6
4
2
-2
-4
-6
0 10 20 30 40 50 60 70 80 90 100
图1采用维纳滤波，噪声方差为1时，真实轨迹、观测样本及估计轨迹的比较从图1可以看出，利用Wiener滤波器通过测量信号估计x（n）的波形，得
到的估计轨迹接近于真实轨迹。
当噪声方差为1时，估计轨迹与真实轨迹间的误差很小。
（2）Wiener滤波器的阶数取101阶，观测点数取100，u（n）的方差为1
时，通过仿真得到平均误差如下图2所示：
1.5
0.5
-0.5
-1.5
图2采用维纳滤波，噪声方差为1时的平均误差
从图2可以看出，利用Wiener滤波器通过测量信号估计x（n）的波形，当
噪声方差为1时，得到的估计轨迹与真实轨迹的平均误差较小，而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。
（3）Wiener滤波器的阶数取101阶，观测点数取100，u（n）的方差为4时，
通过仿真得到真实轨迹、观测样本及估计轨迹的比较图形如下图3所示：
15
10
5
-5
-10
-15
图3采用维纳滤波，噪声方差为4时，真实轨迹、观测样本及估计轨迹的比较从图3可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的误差变
大了。
利用Wiener滤波器通过测量信号估计x（n）的波形，得到的估计轨迹接近
于真实轨迹。
（4）Wiener滤波器的阶数取101阶，观测点数取100，u（n）的方差为4
时，通过仿真得到平均误差如下图4所示：
3
-3
图4采用维纳滤波，噪声方差为4时的平均误差
从图4可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的平均误差变大了。
而且刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。
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）的波形，当
噪声方差为1时，得到的估计轨迹与真实轨迹的平均误差较小。
（3）采用卡尔曼滤波，观测点数取100，u（n）的方差为4时，通过仿真得
到真实轨迹、观测样本及估计轨迹的比较图形如下图7所示：
-20
图7采用卡尔曼滤波，噪声方差为4时，真实轨迹、观测样本及估计轨迹的比较
从图7可以看出，当噪声方差为4时，估计轨迹与真实轨迹间的误差变大了。
利用卡尔曼滤波器通过测量信号估计x（n）的波形，得到的估计轨迹接近
（4）采用卡尔曼滤波，观测点数取100，u（n）的方差为4时，通过仿真
得到平均误差如下图8所示：
8
-8
图8采用卡尔曼滤波，噪声方差为4时的平均误差
从图8可以看出，利用卡尔曼滤波器通过测量信号估计x（n）的波形，当噪
声方差为4时，估计轨迹与真实轨迹间的平均误差变大了。
四、结论
通过以上的仿真结果及其分析可以得到如下结论：
（1）利用Wiener滤波器和Kalman滤波器都可以通过测量信号对x（n）的波
形进行估计，估计轨迹接近于真实轨迹；
（2）随着噪声方差的变大，通过维纳滤波和卡尔曼滤波得到的估计轨迹与真实轨迹间的误差也随之变大；
（3）相比较而言，利用维纳滤波器进行估计的误差稍小于通过卡尔曼滤波器进行估计的误差；
（4）在利用维纳滤波器和卡尔曼滤波器进行估计时，刚开始的时候平均误差相对较大，随着时间的推移，平均误差有逐渐减小的趋势。
五、附件：
1、维纳滤波仿真源程序（噪声方差取1）
clc;
clearall;
maxlag=100;
%%%%%%%%%%%%%%%%%%%%%%%%
N=100;
%观测点数取100x=zeros（N,1）;
y=zeros（N,1）;
var=1;
%%%%%%%%%%%%%%%%%%%%%%%%%列出状态方程x
（1）=randn（1,1）;
%令x（-1）=x（-2）=x（-3）=x（-4）=0x
（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）;
forn=5:
N
x（n）=1.352*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=y'
xk_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））;
%对称化自相关函数矩阵使之成为方阵，滤波器的阶数为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;
%计算滤波的均值,计算滤波误差的均值
fori=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、维纳滤波仿真源程序（噪声方差取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];
%状态变量之间的
增益矩阵Ak
Ck=[1000];
%状态变量与输出信号之间的增益矩阵Ck,（一维:
仅仅对x方向进行估计）
Rk=[1];
%量测噪声协方差阵Pk=[1000
0100
0010
0001];
%噪声的均方误差阵xk=[y
（1）;
Qk=[1];
%%%%%%%%%%%%%%%%%%%%%%%%Kalman滤波开始,估计循环forr=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为估计值
%ende_x=0;
4、卡尔曼滤波仿真源程序（噪声方差取4）
cle
升级会员 