1、估计随机信号的样本自相关序列。先以白噪声为例。
(a) 产生零均值单位方差高斯白噪声的1000个样点。
(b) 用公式：
估计的前100个自相关序列值。与真实的自相关序列相比较，讨论你的估计的精确性。
(c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为自相关的估值，即：
与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列吗？
(d) 再将1000点的白噪声通过滤波器产生1000点的y(n)，试重复(b)的工作，估计y(n)的前100个自相关序列值，并与真实的自相关序列相比较，讨论你的估计的精确性。
仿真结果：
(a)
图1.1 零均值单位方差高斯白噪声的1000个样本点
分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。
(b)
图1.2 的前100个自相关序列值
分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k≠0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列rx(k)=δ(k)比较接近，k≠0时估计值非常接近0，说明了估计的结果是比较精确的。
(c)
图1.3基于Bartlett法的前100个自相关序列值
与(b)的结果相比，同样在k=0时达到峰值，k≠0时0值附近上下波动；估计值的方差比较小，随着k的增大波动幅度逐渐变小，在k较大时它更接近真实自相关序列。即采用分段方法得到的自相关序列的估计值更加接近rx(k)=δ(k)。分析仿真图也可以看出：将样本数据分段，将所有子段的样本自相关的平均值作为自相关的估值时，可以有效的降低自相关估计的方差，而分段样本估计的优点在于，估计自相关序列与实际自相关序列的方差减小，且当分段数越大，估计值越趋向于无偏估计。
(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过程：
是单位方差白噪声。
(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) 计算其自相关序列的估计，并与真实的自相关序列值相比较。
(c) 将的DTFT作为x(n)的功率谱估计，即：
。
(d) 利用所估计的自相关值和Yule-Walker法(自相关法)，估计和的值，并讨论估计的精度。
(e) 用(d)中所估计的和来估计功率谱为：。
(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中所需的信号是一个正弦序列, , 噪声序列和都是AR(1) 过程，分别由如下的一阶差分方程产生：
其中是零均值、单位方差的白噪声，与不相关。
(a) 试用Matlab程序产生x(n)和的500个样点，画出波形图。
(b) 基于x(n)和的500个样点，设计p阶的最优FIR维纳滤波器，由估计，进而估计出，其中阶数p分别取为p=3，6，9，12，试计算各种情况下估计时的平均平方误差(均方误差的样本估计，要叙述估计方案)，并画出对d(n)估计的结果。
(c) 有时辅助观测数据中也会漏入一些d(n)信号，即辅助观测信号不仅是，而是
试针对p=12的情况，分别取几个不同的 值(如0.1, 0.5, 1.0)，研究这时的噪声抑制性能。
(d) 若只有一路观测的1000个样点，你能想办法近似完成对噪声的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。
图6.8 有辅观测数据的维纳噪声抑制器的原理图
仿真结果：
（a）
图3.1 的波形
图3.2 x(n)的500个样点的波形
（b）
基于和的500个样点，可以得到
、
求解Wiener-Hopf方程可以得到最优FIR维纳滤波器。
均方误差的样本估计可以用计算得当p=3、6、9、12时，估计时的平均平方误差分别为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个图：红色曲线代表真实值。蓝色曲线代表估计值。由结果来看，随着滤波器阶数的提高，误差越来越小，对的估计也更加精确了，这是因为阶数越高，使用的自相关序列的值的个数就越多，估计的值也就越准确了。
（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.82、 2.4204，当辅助观测信号受到干扰时，由仿真结果可以看出，干扰的程度越深，即的值越大，估计的性能越差。当漏泄系数时，还可分辨是正弦波形；当漏泄系数时，波形已经基本失真，不能起到分辨效果。
（d）原理框图如下：
如图，采用“自适应滤波”的方法近似完成对噪声的抑制；假设和都是实值的、零均值过程，且互不相关，另外假设是窄带过程，是宽带过程，对，可以近似认为与自相关序列为0，此时：
将观测信号进行延迟产生“参考信号”，将这个参考信号通过自适应滤波器来估计出。因为与不相关，则有
另外，由于到自适应滤波器的输入是，因此其输出
从而
因为和是互不相关，与也是近似不相关的，所以
，从而最后的均方误差变为
可见，最小化等效于最小化，所以自适应滤波器的输出就是的最小均方估计，即实现了噪声抑制。
源程序：
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;
w1'
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;
end
%
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;下载本文
