姓名： 专业：通信与信息系统 学号： 口期：2015.11
实验内容
任务一：
一连续平稳的随机信号兀（/）,自相关函数rv（r） = ^W,信号曲）为加性噪声所干扰，噪 声是白噪声，测量值的离散值揪）为已知，7；=0.02s，-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, H编卡尔曼滤波递推稈序，估计信号兀（\"
的波形。
任务二
设计一维纳滤波器。
（1） 产生三纽•观测数据：首先根据£（/?）=处（刃T）+W（〃）产生信号S（〃），
将其加噪（信 噪比分别为20dB, 10dB, 6dB）,得到观测数据可⑺），X2（z?）,心⑺）。
（2） 佔计兀“），i = 1 , 2 , 3的AR模型参数。假设信号长度为L, AR模型阶数为N , 分
析实验结果，并讨论改变厶，N对实验结果的影响。
实验任务一
1.卡尔曼滤波原理 1.1卡尔曼滤波简介
早在20世纪40年代，开始有人用状态变量模型来研究随机过程，到60年代初，由于 空间技术的发展，为了解决对非平稳、多输入输出随机序列的估计问题，卡尔曼提出了递推 最优估计理论。它用状态空间法描述系统，由状态方稈和量测方程所组成，即知道前一个状 态的估计值和最近一个观测数据，采用递推的算法估计当前的状态值。由于卡尔曼滤波采用 递推法，适合于计
算机处理，并且可以用来处理多维和非平稳随机信号，现已广泛应用于很 多领域，并取得了很好的结果。卡尔曼滤波一经出现，就受到人们的很大重视，并在实践 中不断丰富和完善，其小一个成功的应用是设计运载体的高精度组合导航系统。卡尔曼滤波 具有以下的特点：
（1） 算法是递推的，口状态空间法采用在时域内设计滤波器的方法，因而适用于多维
随机过程的估计；离散型卡尔曼算法适用于计算机处理。
（2） 川递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化
规律，因此信号可以是平稳的，也可以是非平稳的，即卡尔曼滤波适用于非平稳过程。
（3） 卡尔曼滤波采取的谋差准则仍为估计谋差的均方值最小。
1.2卡尔曼滤波的状态方程和测量方程
假设某系统£时刻的状态变量为无，状态方稈和量测方稈（输出方稈）表示为
X
k ~ \\-\\Xk-\\ + Wk-\\
其中，无是状态变量；
为了推导简单，假设状态变量的增益矩阵4不随时间发生变化，叫，叫都是均值为零
的正态白噪声，方差分别是Q和傀，并且初始状态与％，\"都不相关，厂表示相关系数。
即:
叫-表示输入信号是白噪声；叫是观测噪声；儿是观测数据。
叭:
其屮
E
M = 0,云=2,久“ =QQkj * :团坯］== R, Y
k
Vk 比=Rk %
13卡尔曼滤波的递推算法
卡尔曼滤波采用递推算法来实现，其基木思想是先不考虑输入信号叫和观测噪声叫的 影响，得到状态变量和输出信号（即观测数据）的估计值，再用输岀信号的佔计误力加权后 校正状态变呈的佔计值，使状态变量估计误弟的均方值最小。因此，卡尔曼滤波器的关键是 计算出加权矩阵的最佳值。
当不考虑观测噪声和输入信号时，状态方稈和星测方程为
A
A
忌=
AR Xk-\\
显然，由于不考虑观测噪声的影响，输岀信号的估计值与实际值是有误羌的，用儿表
— A
儿=九一儿
为了提高状态估计的质量，用输出信号的估计误差儿来校正状态变量
/ • \\
A
A
A
A
(
A
\\
Xk = Ak Xk-\\+Hk
儿—儿 兀R-i+ H* 儿-c/严
< 丿
=AR\\ 丿
其屮，//攵为增益矩阵，即加权矩阵。经过校正肩的状态变量的估计误差及其均方值分用忌和P,表示，把未经校正的状态变最的估计误差的均方值用P,表示
Xk =Xk-Xk
/
A V
A Y_
P林 -Xk X k 一忌
严\\ 八
/
/ • \\z
.、A
A
P严耳 -Xk x t
一双
E
\\
/
卡尔曼滤波要求状态变量的佔计误并的均方值人为最小，因此卡尔曼滤波的关键即为 人
通过选择合适的使得人取得最小值。首先推导状态变量的估计值忌和状态变量的的
估计误羌忌，然后计算兀的均方值人，通过化简人，得到一组卡尔曼滤波的递推公式：
A
A
(
人、
Xk = Ak Xk-\\ + WJ yk 一CkAk Xk-\\
■ H严 RC：aRC；+Ry
R =佔一+2一1 Pk=(I-HkCk)Pk
假设初始条件 Ak, Ck, Qk, Rk, yk, Xk-\\, PkA 已知，其中 xo = E[XQ] , PQ = var[x0], 那么递推流程如下:
别 2. 卡尔曼滤波递推程序编程思想
题目分析
(1) 由于信号述)为加性噪声所干扰，可知儿=母+«,所以q=l (2) 又因为噪声为白噪声,所以& =q2=S」0)=l (3) 因为= 所以
/H=OC
〃J=YO
sM)= 5X(\"=工戶L =工严严+工严八
〃〃加=0
2=-QC
J=YC
/H=OC W=0 /N=qc
_ 1 _ l-e~2
由此可知，B.=—
，-,,即兀(司一不bG — l)=wS — l),可得到：4•二不】，因为抽样 '1-e z
7-1
间隔 T = 0.025 ,所以：4•之％ 二 e'002。
(4) 因此yv(n) = x(n + l)-e~rxx(n)f 所以
Qk = 6： = rww\\0]= E[w(n)vv(n)] = 1-严
因此Q =1-严04
编程分析
A
由上面的分析可知初始条件4 , C, , Qk , /?,，儿已知,在仿真屮假设兀0 = 0,则
xo=O,
佗=1,由以上参数可得卡尔曼实际递推公式
A
A
(
忌=严°2柏+比儿一严°?
\\
Pk=e-^Pk_^\\-e-^ £=(/一比加
将得到的公式代入前面分析的递推公式，即可进行迭代得到结果母。
3. MATLAB源代码
根据以上分析，编matlab稈序如下：
%%
% ------------------ 卡尔 曼滤波 --------------------
% ----- 说明
%X(k + l)=Ak*X(k)+W(k); %Y(k)=Ck*X(k)+V(k) %%
clear;clc;
%基本参数值
Ak=exp (-0.02) ;Ck=l; Qk=l-exp(-0•04);Rk=l;
%初始值设置
X0=0;P0=l;
%观测值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=l:N
if k==l %k=l时山初值开始计算
P_(k)=Ak*P0*Akf +Qk;
H (k)=P_(k)*Ck,*inv(Ck*P_(k)*Ckf+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>l时，开始递推
%递推公式
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);
end end
M=l:N;T=0.02*M;
%作图，画出x (t)的波形
figure (1)
plot (Tz Yf * r * f * LineWidth * z1); xlabel('t');ylabel ( 'y (t)');
title ( Tv尔曼滤波-测量信号y比)波形 ); grid; figure(2)
plot(T,X, * b * A 'LineWidth * A1); xlabel('t');ylabel ('x (t)');
title ( 尔曼滤波-估计信号x (t)波形，); grid;
1
4.实验结果
t
卡尔曼滤波•估计信号X(t)波形
实验任务二
1. 维纳滤波器原理
维纳■霍夫方程
G/ （k）=》咖.（k-m） =/?（£）* q（k）
〃】=0
当力S）是一个长度为M的因果序列（即一个长度为M的FIR滤波器）时，维纳•霍夫 方程表述为
M-I
匚/依）二工力（加）G依-加）“&）%•依）k=0, 1, 2,...
m=0
定义
h\\ h2
• • ■ \"(o)] s(i) RXd = •
■ ■
鼻(1)
L(o) 2(1)
G(M-1 厂 L(M-2)
L(M-1)
则可写成矩阵的形式，即
2)
RXd Rxxh
对上式求逆，得到
=
h Rxx Rxd
由以上式了可知：若已知期望信号与观测数据的互相关函数及观测数据的白相关函数, 则可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，真接从时域求解因果 的维纳滤波器，当选择的滤波器的长度较大时，计算工作量很大，并且需要计算的逆矩 阵，从而要求的存储量也很人
=
预测是根据观测到对的过去数据来佔计半前或将来的信号值。维纳预测是已知以前时刻 的“个数据兀(〃一1),
x{n - p),估计肖前时刻死，或者未来〃+ N时刻的信
号值，即估计如 N\\ Nno,估计得到的结果仍然要求满足均方误差最小的准则。信号
可以预测是由于信号内部存在着关联性。预测是利用数据前后的关联性，根据其屮一部分推 知其余部分。
一步线性预测的时域解
已知x(n -1), A(H-2), x(n- p),预测,假设噪声v(n) = O ,这样的预测
成为一步线性预测。设定系统的单位脉冲响应为力⑺)。根据现行系统的基本理论，输出信
A
+〃
k=\\
y(〃) = x(〃) =为力(k)x(“ 一 R)
令 cipk =_h(k),则
A
J/； 心)
k=l
预测谋并
A
+\"
+\"
心)=x(n)-血)=心)+ 工％血-k)=工apkx(n -k) k=\\ k=Q
其屮％) = 1
E I如2
=E\\ \\x(n)+y\\a kx{n-k
I |
k=[ |
要使均方误差为最小值，要求
0
/ = 1, 2 , ..., p
又因为，我们可以得到
E^(/7)x(n-/)] = 0
/ = I , 2 ,…，p
所以
G(，)+工％-/) = 0
/ = 1，2, ...» p (1)
k=\\
由于预测器的输出;S)是输入信号的线性组合，所以可得：
E e\\n)x{n) =0
以上说明误差信号与输入信号满足止交性原理，预测误差与预测信号值同样满足止交性理。预测误并的最小均方值
E e\\n
E e (n)x(n)]
=E
*{n)+^apkx*{n-k) x{n)
(2)
= G(0)+
丈％M)
k=\\
由(1) (2)联立方程组，写成矩阵形式可得
原 ru(°) G(l)…「(\") L(o)…匚(卩-1) …
r
「1 ■
a
E
nin
L(l) …
rp\\
0
…G(。)
_M G(P-1)…u(°) _
这就是有名的Yule-Walker (维纳•霍夫)方程。
• • a
L \"」
• • •
0
2. 实验编程思想
在木实验中，首先根据要求产生加噪不同的观测数据州S), x2(n), x3(n),然后可 利用已知条件代入Yule-Walker方程，即可求解AR模型参数。
在木实验屮，假设d = 0.2,信号心)的初值$(1) = 0。
3. MATLAB 代码
function Wiener_predict(Lr N) %clc;clear;
%信噪比
SN1 = 6; SN2 = 10;SN3
20;
%产生信号s(n)
a=0.2;W = random(1 norm19
0, 1, L, 1) ; S(l) = 0;
S (n) = a * S(n-1) +
end
&产生观测信号
Am = sum(abs (S) .A2)
/ L;
P2 = Am / (10A (SN2/20));P3 =
Am /
Pl = Am / (10A (SN1/20));
(10A(SN3/20));
VI = random ( T norm1 z 0, Pl, L, 1) ; V2 = random ( 1 norm1 r 0f P2Z Lz 1) ;V3 =random(1 norm1f 0z P3‘ L, 1);
for n=l:L
XI (n) = S (n) + VI (n); X2 (n) = S (n) + V2 (n);
X3 (n) = S (n) + V3 (n);
!
11111end subplot(2z 2z1);plot (Sz b);title( {h S(n) );ylabel( 帕丿变 ; grid
end
for i = 1 : N + 1 on;
subplot(2,2,2);plot (XI,
1);grid on;
•b* ) ; title『观测信号 XI (n) 1 ) ; ylabel (1 幅度
subplot(2,2,3);plot (X2,
B1) ; title (f 观测借号 X2 (n) 1); ylabel (1 幅
1);grid on;
subplot (2f 2f 4);plot(X3Z •b1)
观测信号 X3(n) 1 ) ; ylabel 1);grid on;
fprintf ('\\n对XI信号来说 阶模型参数和误差分别为:\\nf ) AR(X1,N); fprintf ( 1 \\n对X2信号来说 阶模型参数和氏差分别为:\\n-)AR(X2,N); fprintf ( 1 \\n对X3信号来说 阶模型参数和误差分别为:\\n>)AR(X3zN);
function AR(X’N) L = length (X);rx zeros (lz N + 1);R = zeros(N + 1, N + 1); for i = 1 : (N +
1)
for j = 1 : (L -
i + 1)；
sum = sum
+ X(j) * X( j + i - 1)；
end
rx(i) = sum /
(L 一 i + 1); R(i, 1: (i-1)) = rx((i-1) :-l:l); R(i, i: (N+l)) = rx(1: (N - i + 2)); end
zx = rx (2: (N+l));ap = inv(R(1:1:N)) * (-zx) 1; a = [1f ap1];e = rx (1) + zx * ap;
disp (『AR 系数：1
r num2str (a) ] ) ; disp ( [ 1 均方误左：1
f num2str (e)]);
f 幅度
(
function Wiener_newl(L,N)
%%产生三纽•观测数据 $信噪比(dB)
SNR1 = 20; SNR2 = 10; SNR3 = 6;
%产生信号s (n)
a=0.2;W = random('norm', 0, 1, L, 1) ; S (1) = 0; for n = 2 : L
S (n) = a * S (n-1) + W (n);
end
%加噪声产生观测限号
Xl= awgn(SfSNR1f 1 measured 1f 1 linear 1); X2= awgn(SfSNR2 f 1 measured1f 1 linear 1); X3= awgn(S,SNR3仔 1 measured1r 1 linear 1);
subplot (2,2,1); plot (Sz 1 b 1 ) ; title ( 1 {H V S (n) 1 ) ; ylabel ( 1 幅卩2 1 ) ; grid on;
subplot(2,2,2);plot(XI, on;
subplot(2f 2f 3);plot(X2, on;
subplot(2,2,4);plot(X3Z on;
1 b1 ) ; title (1 1 b 1 ) ; title ('
观测信 v XI (n) 1 ) ; ylabel ( 1 幅度,);grid
f) ;t itle(f 观测信号 X2(n) 1 ) ; ylabel (!幅度 T ;grid
观测信 v X3 (n) 1 ) ; ylabel ( 1 幅度 ) ; grid
1
%%估计AR模型参数
fprintf ( 1 \\n对XI信号来说 fprintf (1 \\n对X2信号来说 fprintf (1 \\n对X3疽号来说
阶模型参数和谋差分别为: \\n!);AR(X1,N); 阶模型参数和误差分别为: \\nf);AR(X2,N); 阶模型参数和谋差分别为: \\nf);AR(X3,N);
function AR(X’N)
L = length(X);rx = zeros(lz N + 1);R = zeros(N + lz N + 1); for i = 1 : (N + 1)
sum = 0; for j = 1 :
(L-i+l);
sum = sum + X(j)
end
* X(j + i - 1);
rx(i) = sum / (L - i + 1); end
for i = 1 : N + 1
R(i, 1: (i-1)) = rx ( (i-1) ：-l:1); R(i, i: (N+l)) = rx (1: (N - i + 2)); end
zx = rx (2: (N+l)); ap = inv(R(1:N, 1:N))
* (-zx) * ;
a = [1, ap1]；e = rx (1) + zx * ap;
disp ( [ 1 AR 系数：1 , num2str (a) ] ) ; disp ( [ 1 均方课差：1 / num2str (e)]);
4. 实验结果与分析
(1) 观测数据产生
原始信号S(n)
观测信号X3(n)
图1•原始信号与观测信号(L=50)
（2）模型阶数N对实验结果的影响
N=1
对XI信号来说N阶模型参数和误茅分别为： 对XI信号来说N阶模型参数和误差分别为：
AR 系数：1
均方误差:1.1289
-0.27766
对X2信号來说N阶模型参数和误并分别为：
AR 系数：1 -0.29326
均方误差:0.97283
对X3信号来说N阶模型参数和误差分别为：
AR 系数：1
均方误差:1.0531
-0.26441
N=2
对XI信号来说N阶模型参数和误差分别为：
AR 系数:1 -0.34494 0.2854
均方误差：1.0958
对X2信号来说N阶模型参数和误差分别为：
AR 系数：1
均方误差:1-1639
-0.1696 0.10742
对X3信号来说N阶模型参数和误差分别为：
AR 系数：1 -0.19532 0.17033
均方误差:0.92331
N=3
对XI信号来说N阶模型参数和误差分别为：
AR 系数：1 -0.10673 0.050928 -0.19364
均方误差：1.4197
对X2信号來说N阶模型参数和误并分别为：
AR 系数：1 -0.35451 0.62013 -0.75585
均方谋差：0.95739 对X3信号来说N阶模型参数和误差分别为:
AR 系数：1 -0.12221 0.14428 -0.34185
均方误旁：0.99317
N=5
对XI信号来说N阶模型参数和误差分别为:
AR系数：1 -0.35515 均方误差：1.2405
0.56619 -0.54005 0.65254 -0.51327
对X2信号来说N阶模世参数和误并分别为:
AR系数：1
-0.27343 0.10227 0.028944 0.21289 -0.2508
均方误差:1.3557
对X3信号来说N阶模型参数和误羌分别为:
AR系数：1
-0.36594 0.41414 -0.41665 0.66894 -0.60712
均方误差:1.1025
分抓 由以上实验结果可知: 在数据的长度一定的条件下, 改变AR模型的阶数，均方误差会
改变，当阶数在某个值时, 均方误差的值最小，因此滤波器的阶数对实验结果有很人影响。 在木次实验屮，仿真情况有限，在以上仿真中我们可以看到当模型阶数N为某一I古I定值时, 均方误芳明显较小。
（3）信号长度L対实验结果的影响
L=100
对XI信号来说N阶模型参数和误差分别为:
AR 系数：1 0.022914 0.0078698
均方谋差：1.2033
対X2信号来说N阶模型参数和误差分别为:
AR 系数：1 0.017414 0.19629
均方误差：1.1607
对X3信号来说N阶模型参数和误差分别为:
AR 系数：1 0012789 0.13086
均方误差:1.1483
L=200
对XI.信号来说N阶模型参数和误差分别为:
AR 系数：1
均方误差:13371
017679 0.073726
対X2信号来说N阶模型参数和误差分别为：
AR 系数：1 -0.2649 0.16751
均方误差：0.99844
对X3信号來说N阶模型参数和误養分别为：
AR 系数：1 027145 0.17666
均方误差：0.99289 分析：
由以上仿真结果可知，实验中存在误差，但仍然可以看出，随着信号长度的增加，均方 误并减小，预测更准确。
L=100, N=1
观测信号 X2(n),SNR=10dB
0 50 100
XI信号:N阶模型参数和预测误芳的最小均方值分别为: AR 系数：1 -0.15954
预测误弟的最小均方值：1.0612
X2信号:N阶模世参数和预测谋并的最小均方值分别为: AR 系数：1
-0.1682
预测谋茅的最小均方值：1.1551
X3信号:N阶模型参数和预测误差的最小均方值分别为: AR 系数：1
-0.1161
预测误养的最小均方值：1.2883
N=2
XI ffi号:N阶模型参数和预测误養的最小均方值分别为: AR 系数：1
-0.20658
0.25733
预测误羌的最小均方值：0.98824
X2信号:N阶模型参数和预测误并的最小均方值分别为: AR 系数：1
-0.10574
0.15188
预测误差的最小均方值:1.0349
X3信号:N阶模型参数和预测误差的最小均方值分别为: AR 系数：1
-0.17376
0.23089
预测误羌的最小均方值：1.2323
N=5
XI信V： N阶模型参数和预测误差的最小均方值分别为： AR 系数:1
-0.21587
0.22397
-0.24306
0.24469
-0.13453
预测误弟的最小均方值：0.88869
X2信号:N阶模樂参数和预测误并的最小均方值分别为： AR 系数：1
-0.2537
0.31482
-0.19014
0.12243
0.040983
预测误旁的最小均方值：0.89859
X3信号:N阶模型参数和预测误差的最小均方值分别为： AR 系数：1
-0.27812
0.33384
-0.27881
0.25752
-0.11447
预测谋养的最小均方值：1.0384
N=10
XI信号:N阶模型参数和预测谋羌的最小均方值分别为: AR 系数：1
-0.085746
0.17042
-0.24887
03049
-0.29013
0.34871
-0.33129
0.48704 -0.39942 0.37265
预测误并的最小均方值：0.96799
X2信空N阶模型参数和预测误差的最小均方值分别为： AR 系数：1 010183 0.20689 -0.18967 0.023448 0.0093937
0.18436
-0.1808
0.13523
预测误差的最小均方值：0.98478
X3信号:N阶模型参数和预测误差的最小均方值分别为： AR 系数：1 021104 0.51543 -0.60156 -0.90162
0.95174
-0.74381
0.6246
预测误差的最小均方值:0.96821
L=1000
XI信号:N阶模型参数和预测误差的最小均方值分别为： AR 系数：1
-0.29706
0.2975
-0.25122
预测误差的最小均方值:0.99801
X2信号:N阶模型参数和预测误弟的最小均方值分别为： AR 系数：1
-0.23764
0.2173
-0.15246
预测误差的最小均方值：1.0381
X3信号:N阶模型参数和预测误差的最小均方值分别为： AR 系数：1
-0.29651
0.28487
-0.21534
预测误差的最小均方值:1.0984
L=500
XI信号:N阶模型参数和预测误差的最小均方值分别为: AR系数： 1
-0.34838
0.37975
-0.33385
预测误养的最小均方值： 1.1584
X2信号:N 阶模熨参数和预测误并的最小均方值分别为:AR系数： 1 -0.32197 0.34402 -0.24965
预测误差的最小均方值:1.1777
0.17418
0.76905
0.30972
0.2091
0.26745
0.3446
0.26156
--0.82715
-0.26171
-0.16053
-0.21993
-0.2554-0.26372-0.011321
0.92023
X3信号:N阶模型参数和预测误芳的最小均方值分别为： AR 系数：1
-0.29053
0.28967
-0.19091
0.17643
-0.14395
预测课茅的最小均方值:1.405
因篇幅问题不能全部显示，请点此查看更多更全内容