实验三、随机信号的功率谱估计方法.docx

随机信号的功率谱估计方法

一、 实验目的
1、 利用自相关函数法和周期图法实现对随机信号的功率谱估计
2、 观察数据长度、自相关序列长度、信噪比、窗函数、平均次 数等谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。
3、 学习使用FFT提高谱估计的运算速度。
4、 体会非参数化功率谱估计方法的优缺点。
二、 实验原理
假设信号x(n)为平稳随机过程，其自相关函数定义为
∅m≜E{x*(n)x(n+m)} (3-1)
其中E表示取数学期望，*表示共轭运算。根据定义，x(n)的功率谱密度与自相关序列存在下面关系：
(3-2)
(3-3)
但是，实际中我们很难得到准确的自相关序列，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。
目前，常用的线性谱估计方法有两种，即相关函数法和周期图方法，本实验对这两种方法分别予以讨论。
1． 自相关函数法
假设我们已知随机信号x(n)的M长的自相关序列{}，利用自相关函数法可以得到x(n)的功率谱估计：
(3-4)
利用窗函数，上式又可表达为
(3-5)
其中，为矩形窗函数，定义为
(3-6)
因此，实际上是真正功率谱与窗函数傅立叶变换的卷积。矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣。为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为
(3-7)
这种窗函数可以有效的抑制旁瓣，但是，此时主瓣宽度增大，从而降低了谱估计的分辨率，这种主瓣和旁瓣之间的矛盾在线性谱估计方法中是无法解决的。
2． 周期图方法
假设已知随机信号的N个样本，利用周期图方法，信号x(n)的功率谱估计为
(3-8)
利用上述方法得到的谱估计方差与信号的功率谱平方成正比，为了减小它的方差，可以将信号序列进行分段处理，然后再求各分段结果的平均，这就是平均周期图方法，即Bartlett方法。
（1）Bartlett 平均周期图方法
将一个随机序列 (0≤n≤N）分成K段，每段长度为L，各段之间互不重迭，因而N=KL，可以想到，第i段的信号序列可表示为
(3-9)
对于每一段的周期图又可写成
, (3-10)
于是，功率谱估计定义为
(3-11)
因此，对于固定的记录长度来讲，分段数K增大可使谱估计的方差减小，但是由于L的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率也是矛盾的。
除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以得到较平滑的谱估计，当然，相应的分辨率也有所下降。
（2）平滑平均周期图方法
这时一种改进的Bartlett周期图方法，它特别适用于FFT直接计算功率谱估值。将长度为N的平稳随机信号序列x(n)分成K段，每段长度为L，即L=N/K。但这里在计算周期图之前，先用窗函数给每段序列加权，K个修正的周期图定义为
, (3-12)
其中U表示窗函数序列的能量，
(3-13)
在这种情况下，功率谱估计可按下面表达式给出：
(3-14)
本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计：
(3-15)
其中NS为正弦个数，ωi,φi和ai分别为第i个正弦信号的数字频率、相位和幅度，φi随机的分布在（0, 2π）之间，w(n)为零均值方差等于σw2的复高斯白噪声。
三、实验内容和步骤
1． 仔细阅读有关线性谱估计的内容，根据给出的框图3-1编制自相关函数法谱估计的程序。运行程序，输入N=100,M=10, ω1=0.6π,φ1=0,a1=1,σw2=0选择矩形窗。观察谱峰位置是否正确（注意：由于窗效应可能引起谱估计的非正定）。
如下图所示，

图（1）
由上图（1）的仿真结果看出，谱峰的位置处于0.6π处，故估计较精确。

开始
输入参数：数据长度N，自相关函数个数M，平均次数K
信号产生：输入正弦个数Ns，每个正弦信号的数字频率、相位和幅度，白噪声信号的方差σw2，按照公式（3.16）产生复正弦加白噪声信号的N个采样
结束
自相关函数法
由N个x(n)估计出自相关序列(M长),并对此自相关序列加矩形窗或Hamming窗，利用公式(3.5)计算[0，2π）之间的128个功率谱抽样点
周期图方法
输入 FFT 点数 ，按照公式(3.11 )、(3.12 )和(3.13)计算NF点功率谱

图3-1

2． 观察并记录参数变化对谱估计性能的影响。
（1）改变M=5，其它输入同步骤1，观察功率谱估计的主瓣宽度和旁瓣大小随自相关序列长度的变化情况。

图（2）
上图（2）为M=5是，矩形窗函数的功率谱估计的仿真图。比较M=5与M=10时的图形可以看出：M越小，主瓣宽度越大，分辨率越低，谱峰高度越低低。

（2）选择窗函数为Hamming窗，其它输入同步骤1，观察不同的窗函数对谱估计性能的影响。

图（3）

比较图（1）图（3）可以看出，图（3）的主瓣宽度增加了近一倍，旁瓣相应也减少了许多，同时发现Hamming窗的分辨率降低了。

（3）改变，其它输入同步骤1，观察初始相位的变化对谱估计性能的影响。

图（4）
比较图（1）、（4）可以看出，说明相位变化对谱估计的性能没有影响。

（4）改变σw2=1，其它输入同步骤1，观察信噪比变化对谱估计性能的影响。
（5）改变N=10，，其它输入同步骤1，结合(4)的内容，观察数据长度及信噪比对谱估计性能的影响。

图（5） 图（6）
如上图（5）为σw2=1时，加矩形窗普估计的仿真图，对比图（1）可以看出，当增大高斯白噪声的方差时，对谱估计的峰值没有产生影响，但对旁瓣的宽度和幅值有影响。
如上图（6）所示，在N较大时，改变信噪比对谱估计影响比较小，而当N较小时，改变信噪比对谱估计的影响较大。有图（6）可以看出此时频谱峰值已偏离了0.6π，同时旁瓣的幅值增大了许多。

3．运行程序，输入，，选择矩形窗，调整自相关序列长度M，使得两个正弦频率分量临界分辨出来，纪录此时的M值，并绘制此时的功率谱图。同样，在加Hamming窗的情况下，记录使两个正弦频率分量临界分开的M值，并绘制此时的功率谱图。

图（7） 图（8）
图（7）、（8）分别是M=7和M=8时，加矩形窗的普估计仿真图，可以看出M=8时可以使两个正弦频率分量临界分辨出来。

图（9） 图（10）
图（9）、（10）分别是M=9和M=10时，加Hamming窗的普估计仿真图，可以看出M=10时可以使两个正弦频率分量临界分辨出来。

4．运行自相关函数法谱估计程序，输入，选择矩形窗，观察利用自相关函数法得到的白噪声信号谱估计。改变M=3，20，观察M的变化对白噪声谱估计的影响。
如下图（11）、（12）、（13）分别为M=3、10、20时，加矩形窗的普估计的仿真图，可以看出，M越大，谱估计的各个分量越越清晰。

图（11）

图（12） 图（13）
5. 根据框图，编制周期图法谱估计程序。自程序FFT可以直接调用Matlab中的函数。运行程序，输入。选择窗函数为矩形窗，，观察谱峰位置是否正确，并与步骤1结果比较。
图（14）
对比图（1），上图中的谱峰位置同样也出现在0.6π处，谱的幅值明显降低，主瓣宽度增加，旁瓣数量减少。

6．利用周期图方法重复步骤2、3、4的内容，这里L=N/K相当于自相关函数中的M，观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。
1）、

图（15） 图（16）
如图（15）、（16）分别为N=100，K=20时，加矩形窗和汉明窗后的普估计图。N/K减小，周期图法谱估计的主瓣宽度越大，分辨率降低，幅值降低，而且旁瓣数量减少。这和自相关函数法中减小M时，谱估计出现的规律是一致的。
加汉明窗时的普估计，对比于图（16）和图（3）可以看出，周期图法得到的谱的主瓣宽度加倍，旁瓣减小许多。周期图法谱估计和自相关法谱估计出现的规律相似。

2）、改变，其它输入同步骤1，观察初始相位的变化对谱估计性能的影响。

图（17）
对比图（14）以看出，改变初始相位的值对谱估计没影响。

改变σw2=1，其它输入同步骤1，观察信噪比变化对谱估计性能的影响。

图（18）
如上图所示，增大高斯白噪声的方差时，对谱估计的峰值没有影响，但却使旁瓣增大。

3）、改变N=50，，其它输入同步骤1，结合(4)的内容，观察数据长度及信噪比对谱估计性能的影响。

图（19）
由图可知，当N减小，同时增大噪声方差时，矩形窗的周期图谱估计的主瓣宽度增加。

4）、运行程序，输入，，选择矩形窗，调整自相关序列长度M，使得两个正弦频率分量临界分辨出来，纪录此时的M值，并绘制此时的功率谱图。同样，在加Hamming窗的情况下，记录使两个正弦频率分量临界分开的M值，并绘制此时的功率谱图。
固定K=10，分别取N=30和40，运行程序可以得到下图（20）、（21）所示的加矩形窗普估计的图形。

图（20） 图（21）

图（20）、（21）分别是N/K=3、4时，加矩形窗的普估计仿真图，可以看出N/K=4时可以使两个正弦频率分量临界分辨出来。只不过此时谱峰的位置已不在0.6π和0.8π处了，普估计不准确了。

5）、固定K=10，分别取N=80和90，运行程序可以得到下图（22）、（23）所示的加汉明窗普估计的图形。

图（22） 图（23）

图（22）、（23）分别是N/K=8、9时，加汉明窗时的普估计仿真图，可以看出N/K=9时可以使两个正弦频率分量临界分辨出来。此时估计仍发生了错误。

6）、运行自相关函数法谱估计程序，输入，选择矩形窗，观察利用自相关函数法得到的白噪声信号谱估计。改变M=3，20，观察M的变化对白噪声谱估计的影响。
固定N=100，分别取K=50，10，5，即N/K=2、10、20，运行程序，对应的有图（24）、（25）、（26）的仿真结果。

图（24）

图（25） 图（26）
对比以上三图可以看出，N/K的值越大，谱估计的各个分量就越清晰，也就越容易分辨了。
四、实验主要结论
1、使用不同的窗函数谱估计的质量是不一样的，矩形窗的主瓣较窄，分辨率较好，但方差较大，噪声水平较高；而Hamming窗的主瓣较宽，分辨率较低，但方差较小，噪声水平较低。
2、周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计：增大 M(或 N/K)值，会使主瓣宽度减小，即分辨率增加；反之，M 越小，主瓣越宽，分辨率越小。另一方面，对于周期图法，仍存在着频率分辨率低、方差性能不好的问题，原因是谱估计时对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，当数据很短时，这个问题更为突出。
3、初始相位对谱估计没有影响。

五、思考题

1. 证明：式（3-4）可以表示为：

其中，Real[·]表示取实部。
证明：
Ρw=m=-M+1M-1ϕme-jwm
=m=0M-1ϕme-jwm+m=-M+10ϕme-jwm-ϕ(0)

=m=0M-1ϕme-jwm+m=0M-1ϕ-mejwm-ϕ(0)
又自相关函数为偶函数，则上式可化简为：
Ρw=m=0M-1ϕme-jwm+n=0M-1ϕmejwm-ϕ(0)=m=0M-1ϕm（e-jwm+ejwm）-ϕ(0)
即：

2、 已知实信号的 N 个样本，可以定义的估计如下：
ΡAw=n=-N+1N-1ϕme-jwm

其中
试证明：

证明：
ΡAw=m=-N+1N-1ϕme-jwm

=m=-N+1N-11Nn=0N-m-1x*(n)xn+me-jwm
令 m+n=k,上式为：
=1Nk=0n-1xk n=0N-1x*ne-jwk-n
=1N[k=0n-1x(k)e-jwk]· [n=0N-1xnejwn]*
即有：

附录一
clear
clc
N=100;
M=10;
Ns=1;
w1=0.6*pi;
pha1=0;
a1=1;
sigma_w2=0;
L=256;
wn=sqrt(sigma_w2)*(randn(N,1));
xn=zeros(N,1);
for n=1:N
xn(n)=a1*exp(j*(w1*n+pha1))+wn(n);
end
r_xx=xcorr(xn,'unbiased');

win_p=2*M-1;
W_rect=rectwin(win_p);
rect1=r_xx(N-M+1:N+M-1).*W_rect;
rect2=fft(rect1,L);
rect3=abs(rect2(1:L/2));
figure(1),clf
stem((0:L/2-1)/(L/2),rect3,'b')%归一化
title('矩形窗功率谱估计');

W_han=hamming(win_p);
han1=r_xx(N-M+1:N+M-1).* W_han;
han2=fft(han1,L);
han3=abs(han2(1:L/2));
figure(2),clf
stem((0:L/2-1)/(L/2),han3,'b')
title('hamming窗功率谱估计');

附录二
clear
clc
N=100;
M=10;
Ns=2;
w1=0.6*pi;
w2=0.8*pi;
pha1=0;
pha2=0;
a1=1;
a2=1;
sigma_w2=0;
L=256;
wn=sqrt(sigma_w2)*(randn(N,1));
xn=zeros(N,1);
for n=1:N
xn(n)=a1*exp(j*(w1*n+pha1))+a2*exp(j*(w2*n+pha2))+wn(n);
end
r_xx=xcorr(xn,'unbiased');

win_p=2*M-1;
W_rect=rectwin(win_p);
rect1=r_xx(N-M+1:N+M-1).*W_rect;
rect2=fft(rect1,L);
rect3=abs(rect2(1:L/2));
figure(1),clf
stem((0:L/2-1)/(L/2),rect3,'b')%归一化
title('矩形窗功率谱估计');

W_han=hamming(win_p);
han1=r_xx(N-M+1:N+M-1).* W_han;
han2=fft(ham1,L);
han3=abs(ham2(1:L/2));
figure(2),clf
stem((0:L/2-1)/(L/2),han3,'b')
title('hamming窗功率谱估计');

附录三
clear
clc
N=100;
Nf=128;
K=30;
Ns=1;
w1=0.6*pi;
pha1=0;
a1=1;
sigma_w2=1;
L=N/K;
wn=sqrt(sigma_w2)*(randn(N,1));
xn=zeros(N,1);
for i=1:N
xn(i)=a1*exp(j*(w1*i+pha1))+wn(i);
end

%加矩形窗
w_rect=rectwin(L);
for i=1:K
x_w_rect=xn((i-1)*L+1:i*L,:).*w_rect;
x_w_dft=fft(x_w_rect,Nf);
ixrec(i,:)=x_w_dft.^2/u/L;
end
pw_rect=abs(sum(ixrec)/K);

%加hamming窗
w_han=hamming(L);
for i=1:K
x_w_han=xn((i-1)*L+1:i*L,:).*w_han;
x_han_dft=fft(x_w_han,Nf);
ixhamm(i,:)=x_han_dft.^2/Nf/L;
end
pw_han=abs(sum(ixhamm)/K);

figure(1),clf
stem((0:Nf/2-1)/(Nf/2-1),pw_rect(1:Nf/2),'b')
title('加矩形窗谱估计');
figure(2),clf
stem((0:Nf/2-1)/(Nf/2-1),pw_han(1:Nf/2),'b')
title('加hamming窗谱估计');