现代数字信号处理实验报告

现代数字信号处理实验报告
1、估计随机信号的样本自相关序列。

先以白噪声()x n 为例。

(a) 产生零均值单位方差高斯白噪声的1000个样点。

(b) 用公式：
999
1ˆ()()()1000x n r k x n x n k ==-∑
估计()x n 的前100个自相关序列值。

与真实的自相关序列()()x r k k δ=相比较，讨论你的估计的精确性。

(c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为()x n 自相关的估值，即：
999
00
1ˆ()(100)(100) , 0,1,...,991000x m n r k x n m x n k m k ===+-+=∑∑
与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列
()()x r k k δ=吗？
(d) 再将1000点的白噪声()x n 通过滤波器1
1
()10.9H z z
-=
-产生1000点的y (n )，试重复(b)的工作，估计y (n )的前100个自相关序列值，并与真实的自相关序列()y r k 相比较，讨论你的估计的精确性。

仿真结果：
(a)
图1.1 零均值单位方差高斯白噪声的1000个样本点
分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。

(b)
图1.2 ()
x n的前100个自相关序列值
分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k ≠0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列r x (k)=δ(k)比较接近，k ≠0时估计值非常接近0，说明了估计的结果是比较精确的。

(c)
图1.3基于Bartlett 法的前100个自相关序列值
与(b)的结果相比，同样在k=0时达到峰值，k ≠0时0值附近上下波动；估计值的方差比较小，随着k 的增大波动幅度逐渐变小，在k 较大时它更接近真实自相关序列()()x r k k δ=。

即采用分段方法得到的自相关序列的估计值更加接近r x (k)=δ(k)。

分析仿真图也可以看出：将样本数据分段，将所有子段的样本自相关的平均值作为()x n 自相关的估值时，可以有效的降低自相关估计的方差，而分段样本估计的优点在于，估计自相关序列与实际自相关序列的方差减小，且当分段数越大，估计值越趋向于无偏估计。

(d)
图1.4 y(n)的前100个自相关序列值与真实值的对比从图中可以看出在k=0时估计与真实的自相关序列之间有较小的误差，随着k的增大，估计得到的值有较大的波动，存在一定误差。

源程序
clc
clear
%%产生1000个高斯白噪声的样本点
x=randn(1,1000);
K=1000;
figure(1);
k=0:K-1;
stem(k,x,'.'); %绘制1000个高斯白噪声
title('零均值单位方差高斯宝噪声，1000个样本点');
xlabel('k');ylabel('x[k]');
mean_x=mean(x)
var_x=var(x)
%%
for k=0:99
for n=k+1:1000
y_ess(n)=x(n)*x(n-k);
end
r_ess(k+1)=sum(y_ess)/1000;
end
figure(2);
k=[0:99];
stem(k,r_ess,'.');
title('根据样本点估计出的前100自相关序列值');
xlabel('k');ylabel('r_ess[k]');
hold on;
realvalue=[1,zeros(1,99)];
stem(k,realvalue,'r','.');
legend('根据样本点估计出的前100自相关序列值','真实的自相关序列'); error1=r_ess-realvalue;
mean_error_b=mean(error1)
var_error_b=var(error1)
%%
for k=0:99
for m=0:9
for n=k+1:100
y_ess2(m+1,n)=x(n+100*m)*x(n-k+100*m);
end
end
r_ess2(k+1)=sum(sum(y_ess2))/1000;
end
figure(3);
k=0:99;
stem(k,r_ess2,'b.');
hold on;
realvalue2=[1,zeros(1,99)];
stem(k,realvalue2,'r.','.');
title('Bartlett法估计功率谱方法得出的前100个自相关序列值');
xlabel('k');ylabel('r_ess2[k]');
legend('Bartlett法估计功率谱方法得出的前100个自相关序列值','真实的自相关序列');
error2=r_ess2-realvalue2;
mean_error_c=mean(error2)
var_error_c=var(error2)
%%
y=zeros(1,1000);
B=[1];
A=[1,-0.9];
y=filter(B,A,x);
r_ess3=zeros(1,100);
for k=0:99
for n=(k+1):1000
r_ess3(k+1)=r_ess3(k+1)+y(n)*y(n-k);
end
r_ess3(k+1)=r_ess3(k+1)/1000;
end
figure(4);
stem(r_ess3,'.');
title('y[n]前100个自相关序列估计值');
xlabel('k'),ylabel('r_ess3(k)');
hold on;
p=[1,zeros(1,99)];
h=filter(B,A,p);
for i=1:100
h1(i)=h(101-i);
end
rh=conv(h,h1);
rh=rh(100:199);
realvalue3=conv(p,rh);
realvalue3=realvalue3(1:100);
stem(realvalue3,'r.','.');
legend('y[n]前100个自相关序列估计值','y[n]的真实自相关序列');
2、计算机练习2：AR过程的线性建模与功率谱估计。

考虑AR过程：
x n a x n a x n a x n a x n b v n =-+-+-+-+
()(1)(1)(2)(2)(3)(3)(4)(4)(0)() v n是单位方差白噪声。

()
(a) 取b(0)=1, a(1)=0.-7348, a(2)=-1.8820, a(3)=-0.7057, a(4)=-0.8851，产生x(n)的N=256个样点。

(b) 计算其自相关序列的估计ˆ()
r k，并与真实的自相关序列值相比较。

x
(c) 将ˆ()x r
k 的DTFT 作为x (n )的功率谱估计，即： 1
21
1
ˆˆ()()|()|N j jk j x
x k N P e r
k e X e N
ωωω--=-+==∑。

(d) 利用所估计的自相关值和Yule-Walker 法(自相关法)，估计
(1), (2), (3), (4)a a a a 和(0)b 的值，并讨论估计的精度。

(e) 用(d)中所估计的()a k 和(0)b 来估计功率谱为：2
2
4
1|(0)|ˆ()1()j x
jk k b P e a k e ωω
-==+∑。

(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。

(g) 重复上面的(d)~(f)，只是估计AR 参数分别采用如下方法：(1) 协方差法；(2) Burg 方法；(3) 修正协方差法。

试比较它们的功率谱估计精度。

仿真结果： （a ）
图2.1 x (n )的N =256个样点
（b ）
图2.2自相关序列的估计值与真实的对比
图2.2中估计的自相关序列与真实的自相关序列差异较大；在k>100后，真实的自相关序列就波动得很小，而估计的自相关序列则仍有较大的波动，但总体上来言两者都随着k的增大而逐渐衰减，当k>150时，真实自相关序列逐渐趋于0，而估计出的自相关序列却仍有较大的波动，这是因为估计的点数较少，使得估计精度不够。

（c）
图2.3 估计的功率谱与真实功率谱对比（d）Yule-Walker法(自相关法)
估计的参数值为：
b(0)= 1.1537
[a(1) a(2) a(3) a(4)]=[-0.7174 -1.7952 -0.6387 -0.8167]
图2.4 Yule-Walker法估计的功率谱与真实功率谱对比
Yule-Walker法(自相关法)估计的参数，其系数的符号正确，但数值大小相差较大，并且多次实验的相差值也很大，参数估计的精度远远不够。

因此从图2.4中也能看出，该方法得到功率谱与真实的谱相差很大
（e）协方差法
图2.5 协方差法估计的功率谱与真实功率谱对比
采用协方差法估计的参数，其系数与真实的系数非常接近，该方法能够较精确的估计出功率谱。

（f）修正协方差
图2.6 修正的协方差法估计的功率谱与真实功率谱对比采用修正的协方差法估计的参数，其系数虽然没有协方差法和burg法那么精确，但是估计误差也很小。

从图2.6中也能看出，该方法能够较精确的估计出功率谱。

（g）Burg算法
图2.7 burg法估计的功率谱与真实功率谱对比
采用burg估计的参数，其系数几乎等于真实的系数，分析图2.7中也能看
出，该方法非常精确的估计出功率谱，几乎与真实的功率谱曲线重合。

源程序：
clc;clear;
N=256;
NN=20000;
v1=normrnd(0,1,50,NN);
v=v1(:,1:N);
r=zeros(1,N);
r1=zeros(1,N);
w=0:2*pi/100:2*pi;
P=zeros(1,length(w));
PP1=zeros(1,length(w));
PP2=zeros(1,length(w));
PP3=zeros(1,length(w));
PP4=zeros(1,length(w));
for s=1:50
x1=filter([1],[1,0.7348,1.8820,0.7057,0.8851],v1(s,:));
x=x1(1:N);
for k=1:N
rx(k)=0;
for n=k:N
rx(k)=rx(k)+x(n)*x(n-(k-1));
end
rx(k)=rx(k)/(N);
end
r=r+rx;
for k=1:N
rx1(k)=0;
for n=k:NN
rx1(k)=rx1(k)+x1(n)*x1(n-(k-1));
end
rx1(k)=rx1(k)/(NN);
end
r1=r1+rx1;
for i=1:length(w)
P(i)=P(i)+rx(1);
for n=2:N
P(i)=P(i)+rx(n)*exp(-j*(n-1)*w(i))+rx(n)*exp(j*(n-1)*w(i));
end
end
A=toeplitz(rx(1:4)',rx(1:4));
B=-rx(2:5)';
Ap=A\B;
b0=rx(1);
for i=1:4
b0=b0+Ap(i)*rx(i+1);
end
b0=sqrt(b0);
for i=1:length(w)
P1(i)=1;
for k=1:4
P1(i)=P1(i)+Ap(k)*exp(-j*k*w(i));
end
P1(i)=b0^2/abs(P1(i))^2;
end
PP1=PP1+P1;
A=[];
for k=1:4
c=[];
for l=1:4
rr=0;
for n=5:N
rr=rr+x(n-l)*x(n-k);
end
c=[c;rr];
end
A=[A c];
end
B=[];
for l=1:4
rr=0;
for n=5:N
rr=rr+x(n-l)*x(n);
end
B=[B;rr];
end
B=B*(-1);
Ap=A\B;
for i=1:length(w)
P2(i)=1;
for k=1:4
P2(i)=P2(i)+Ap(k)*exp(-j*k*w(i));
end
P2(i)=x(1)^2/abs(P2(i))^2;
end
PP2=PP2+P2;
A=[];
for k=1:4
c=[];
for l=1:4
rr=0;
for n=5:N
rr=rr+x(n-l)*x(n-k)+x(n-4+l)*x(n-4+k);
end
c=[c;rr];
end
A=[A c];
end
B=[];
for l=1:4
rr=0;
for n=5:N
rr=rr+x(n-l)*x(n)+x(n-4+l)*x(n-4);
end
B=[B;rr];
end
B=B*(-1);
Ap=A\B;
for i=1:length(w)
P3(i)=1;
for k=1:4
P3(i)=P3(i)+Ap(k)*exp(-j*k*w(i));
end
P3(i)=x(1)^2/abs(P3(i))^2;
end
PP3=PP3+P3;
p=4;
ef=zeros(1+p,N);
eb=zeros(1+p,N);
ef(1,:)=x;
eb(1,:)=x;
km=[];
for m=2:p+1
mol=0;
den=0;
for n=m:N
mol=mol+(-2)*ef(m-1,n)*eb(m-1,n-1);
den=den+(ef(m-1,n))^2+(eb(m-1,n-1))^2;
end
km(m-1)=mol/den;
for n=m:N
ef(m,n)=ef(m-1,n)+km(m-1)*eb(m-1,n-1);
eb(m,n)=eb(m-1,n-1)+km(m-1)*ef(m-1,n);
end
end
aa=[1];
for i=1:4
aa=[aa;0];
bb=aa(length(aa):-1:1);
aa=aa+bb*km(i);
end
for i=1:length(w)
P4(i)=1;
for k=2:5
P4(i)=P4(i)+aa(k)*exp(-j*(k-1)*w(i));
end
P4(i)=1/abs(P4(i))^2;
end
PP4=PP4+P4;
end
figure(1)
stem(1:N,x,'.');
title('x[n]的256个样本点');
xlabel('n');ylabel('x[n]');
figure(2)
plot(0:N-1,r/50); hold on;
plot(0:N-1,r1/50,'r');
title('x[n]的256个样本点估计出的前256个自相关序列和真实值'); ylabel('Rx(k)');
xlabel('k');
legend('估计值','真实值');
grid on;
axis([0 256 -40 40]);
hold off;
figure(3)
plot(w/pi,10*log10(P/50)); hold on;
title('功率谱估计');
ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(PP1/50),'r');
plot(w/pi,10*log10(PP2/50),'g');
plot(w/pi,10*log10(PP3/50),'y');
plot(w/pi,10*log10(PP4/50),'k');
aap=[0.7348,1.8820,0.7057,0.8851];
for i=1:length(w)
P5(i)=1;
for k=1:4
P5(i)=P5(i)+aap(k)*exp(-j*k*w(i));
end
P5(i)=1/abs(P5(i))^2;
end
plot(w/pi,10*log10(P5),':');
legend('Rx(k)的DTFT','Yule-Walker'); grid on;
hold off;
figure(4)
plot(w/pi,10*log10(P/50)); hold on;
title('功率谱估计比较');
ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(P5),'r');
legend('Rx(k)的DTFT','真实频谱');
grid on;
hold off;
figure(5)
plot(w/pi,10*log10(PP1/50)); hold on; title('Yule-Walker法功率谱估计比较'); ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(P5),'r');
legend('Yule-Walke法','真实频谱');
grid on;
hold off;
figure(6)
plot(w/pi,10*log10(PP2/50)); hold on; title('协方差法功率谱估计比较'); ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(P5),'r');
legend('协方差法','真实频谱');
grid on;
hold off;
figure(7)
plot(w/pi,10*log10(PP3/50)); hold on; title('修正协方差法功率谱估计比较'); ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(P5),'r');
legend('修正协方差法','真实频谱'); grid on;
hold off;
figure(8)
plot(w/pi,10*log10(PP4/50)); hold on; title('Burg法功率谱估计比较');
ylabel('P(dB)');
xlabel('w/pi');
plot(w/pi,10*log10(P5),'r'); legend('Burg 法','真实谱'); grid on; hold off;
3、计算机练习3：维纳噪声抑制(例6.6的扩展实验)。

假定图6.8中所需的信号()d n 是一个正弦序列0()sin()d n n ωφ=+, 00.02ωπ=, 噪声序列1()v n 和2()v n 都是AR(1) 过程，分别由如下的一阶差分方程产生：
11()0.8(1)()v n v n g n =-+ 22()0.6(1)()v n v n g n =--+
其中()g n 是零均值、单位方差的白噪声，与()d n 不相关。

(a) 试用Matlab 程序产生x (n )和2()v n 的500个样点，画出波形图。

(b) 基于x (n )和2()v n 的500个样点，设计p 阶的最优FIR 维纳滤波器，由2()v n 估计1()v n ，
进而估计出()d n ，其中阶数p 分别取为p=3，6，9，12，试计算各种情况下估计
1()v n 时的平均平方误差(均方误差的样本估计，要叙述估计方案)，并画出对d (n )估计
的结果。

(c) 有时辅助观测数据中也会漏入一些d (n )信号，即辅助观测信号不仅是2()v n ，而是
02()()()v n v n d n α=+
试针对p=12的情况，分别取几个不同的α值(如0.1, 0.5, 1.0)，研究这时的噪声抑制性
能。

(d) 若只有一路观测1()()()x n d n v n =+的1000个样点，你能想办法近似完成对噪声1()
v n 的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。

图6.8 有辅观测数据的维纳噪声抑制器的原理图
仿真结果：
（a ）
图3.1 2()v n 的波形
图3.2 x (n )的500个样点的波形
（b ）
基于()X n 和2()V n 的500个样点，可以得到
122201ˆ()()()N v n r
k V n V n k N -==-∑、1
2
20
1ˆ()()()N xv n r k x n V n k N -==-∑ 求解Wiener-Hopf 22v xv R w r =方程可以得到最优FIR 维纳滤波器。

均方误差的样本估计可以用1
2
1
1
1ˆˆ(()())N n v
n v n N
ξ-==-∑计算得当p=3、6、9、12时，
估计1()v n 时的平均平方误差分别为0.7849、0.2173、0.0747、0.0453。

图3.3 滤波器阶数p=3时的估计值与真实值对比
图3 .4滤波器阶数p=6时的估计值与真实值对比
图3.5 滤波器阶数p=9时的估计值与真实值对比
图3.6 滤波器阶数p=12时的估计值与真实值对比
分析以上4个图：红色曲线代表真实值。

蓝色曲线代表估计值。

由结果来看，随着滤波器阶数的提高，误差越来越小，对)(n
d的估计也更加精确了，这是因为阶数越高，使用的自相关序列的值的个数就越多，估计的值也就越准确了。

（c）
图3.7 漏入的参数a=0.1时的估计值与真实值对比
图3.8 漏入的参数a=0.5时的估计值与真实值对比
图3.9 漏入的参数a=1.0时的估计值与真实值对比
漏泄的参数为0.1、0.5、1.0时，估计误差依次为 0.2312、 1.8892、 2.4204，当辅助观测信号)(2n v 受到)(n d 干扰时，由仿真结果可以看出，干扰的程度越深，即α的值越大，估计的性能越差。

当漏泄系数5.0<α时，还可分辨是正弦波形；当漏泄系数5.0>α时，波形已经基本失真，不能起到分辨效果。

（d ）原理框图如下：
如图，采用“自适应滤波”的方法近似完成对噪声)(1n v 的抑制；假设)(n d 和
)(1n v 都是实值的、零均值过程，且互不相关，另外假设)(n d 是窄带过程，
)(1n v 是宽带过程，对0k k ≥，可以近似认为)(1n v 与)(1k n v -自相关序列为0，此
时：
}
)]()({[)}({)]}()()[({2})]()({[)}({}
|)()()({||})({|2211221212n y n d E n v E n y n d n v E n y n d E n v E n y n v n d E n e E -+=-+-+=-+= 将观测信号进行延迟产生“参考信号”，将这个参考信号通过自适应滤波器来估计出)(n d 。

因为)(n d 与)(1n v 不相关，则有
)}()({2)]}()()[({211n y n v E n y n d n v E -=-
另外，由于到自适应滤波器的输入是)(0k n x -，因此其输出
∑∑==--+--=--=p
k p
k n n k k n v k k n d k w k k n x k w n y 0
0100)]()()[()()()(
从而
∑=--+--=p
k n k k n v n v E k k n d n v E k w n y n v E 0011011)}]()({)}()({)[()}()({
因为)(n d 和)(1n v 是互不相关，与)(01k k n v --也是近似不相关的，所以
0)}()({1=n y n v E ，从而最后的均方误差变为
})]()({[)}({|})({|2212n y n d E n v E n e E -+=
可见，最小化|})({|2n e E 等效于最小化})]()({[2n y n d E -，所以自适应滤波器的输出)(n y 就是)(n d 的最小均方估计，即实现了噪声抑制。

源程序： clc;clear; N=500;
g=normrnd(0,1,1,N); v1=filter([1],[1,-0.8],g); v2=filter([1],[1,0.6],g); figure(1)
plot(0:N-1,v2);
title('辅助噪声观测v2(n)');
ylabel('v2(n)');
xlabel('n');
grid on;
d=sin([0:N-1]*0.02*pi-pi/2);
x=d+v1;
figure(2)
plot(0:N-1,x,'r'); hold on;
plot(0:N-1,d,'-.');
title('观测信号x(n)和正弦信号d(n)'); ylabel('x(n)');
xlabel('n');
legend('观测信号x(n)','正弦信号d(n)'); grid on;
hold off;
%d(n),p=3,6,9,12
for k=1:13
rv2(k)=0;
for n=k:N
rv2(k)=rv2(k)+v2(n)*v2(n-(k-1));
end
rv2(k)=rv2(k)/N;
end
%Rxv2(k)
for k=1:13
rxv2(k)=0;
for n=k:N
rxv2(k)=rxv2(k)+x(n)*v2(n-(k-1));
end
rxv2(k)=rxv2(k)/N;
end
%
p=zeros(1,4);
p(1)=3;p(2)=6;p(3)=9;p(4)=12;
A=toeplitz(rv2(1:p(1)+1)',rv2(1:p(1)+1)); B=rxv2(1:p(1)+1)';
w1=A\B;
A2=toeplitz(rv2(1:p(2)+1)',rv2(1:p(2)+1)); B2=rxv2(1:p(2)+1)';
w2=A2\B2;
A3=toeplitz(rv2(1:p(3)+1)',rv2(1:p(3)+1));
B3=rxv2(1:p(3)+1)';
w3=A3\B3;
A4=toeplitz(rv2(1:p(4)+1)',rv2(1:p(4)+1)); B4=rxv2(1:p(4)+1)';
w4=A4\B4;
%d(n)
v11=filter(w1',[1],v2);
v12=filter(w2',[1],v2);
v13=filter(w3',[1],v2);
v14=filter(w4',[1],v2);
d1=x-v11;
d2=x-v12;
d3=x-v13;
d4=x-v14;
(v1-v11)*(v1-v11)'/N
(v1-v12)*(v1-v12)'/N
(v1-v13)*(v1-v13)'/N
(v1-v14)*(v1-v14)'/N
figure(3)
plot(0:N-1,d1); hold on;
plot(0:N-1,d,'r');
title('FIR维纳滤波器阶数p=3输出结果'); grid on;
hold off;
figure(4)
plot(0:N-1,d2); hold on;
plot(0:N-1,d,'r');
title('FIR维纳滤波器阶数p=6输出结果'); grid on;
hold off;
figure(5)
plot(0:N-1,d3); hold on;
plot(0:N-1,d,'r');
title('FIR维纳滤波器阶数p=9输出结果'); grid on;
hold off;
figure(6)
plot(0:N-1,d4); hold on;
plot(0:N-1,d,'r');
title('FIR维纳滤波器阶数p=12输出结果'); grid on;
hold off;
w2'
w3'
w4'
N=500;
g=normrnd(0,1,1,N);
v1=filter([1],[1,-0.8],g);
v2=filter([1],[1,0.6],g);
d=sin([0:N-1]*0.05*pi-pi/2);
x=d+v1;
d1=[];
a=[0.1 0.5 1.0];
for l=1:length(a)
v0=v2+a(l)*d;
for k=1:13
rv2(k)=0;
for n=k:N
rv2(k)=rv2(k)+v0(n)*v0(n-(k-1));
end
rv2(k)=rv2(k)/N;
%
for k=1:13
rxv2(k)=0;
for n=k:N
rxv2(k)=rxv2(k)+x(n)*v0(n-(k-1));
end
rxv2(k)=rxv2(k)/N;
end
%
pp=12;
A=toeplitz(rv2(1:pp+1)',rv2(1:pp+1));
B=rxv2(1:pp+1)';
w1=A\B;
%
v11=filter(w1',[1],v0);
d1=[d1 x-v11]; (v1-v11)*(v1-v11)'/N end
figure(7)
plot(0:N-1,d1(1:N)); hold on;
plot(0:N-1,d,'r');
title('辅助观测数据中漏入d[n]的0.1'); legend('估计结果','期望值');
axis([0 N -4 4]);
grid on;
hold off;
figure(8)
plot(0:N-1,d1(N+1:2*N)); hold on; plot(0:N-1,d,'r');
title('辅助观测数据中漏入d[n]的0.5'); legend('估计结果','期望值');
axis([0 N -4 4]);
grid on;
hold off;
figure(9)
plot(0:N-1,d1(2*N+1:3*N)); hold on; plot(0:N-1,d,'r');
title('辅助观测数据中漏入d[n]的1.0'); legend('估计结果','期望值');
axis([0 N -4 4]);
grid on;
hold off;。