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

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

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

(k)=δ(k)

x

比较接近，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.4y(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 过程：

()(1)(1)(2)(2)(3)(3)(4)(4)(0)()x n a x n a x n a x n a x n b v n =-+-+-+-+

()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)计算其自相关序列的估计?()x r

k ，并与真实的自相关序列值相比较。 (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.1x (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.12()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

220

1?()()()N xv n r k x n V n k N -==-∑ 求解Wiener-Hopf 22v xv R w r =方程可以得到最优FIR 维纳滤波器。

均方误差的样本估计可以用1

2110

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时的估计值与真实值对比