数字信号处理2实验---随机信号的线性谱估计方法

—随机信号的线性谱估计方法

一、实验目的

1. 利用自相关函数法和周期图法实现随机信号的功率谱估计。

2. 观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣宽度的影响。

3. 学习使用FFT提高谱估计的运算速度。

4. 体会线型谱估计方法的优缺点。

二、实验原理和方法

假设信号x(n)为平稳随机过程，其自相关序列定义为

(m)E[x*(n)x(nm)] (3-1)

其中E表示取数学期望，*表示共轭运算。根据定义，x(n)的功率谱密度P()与自相关序列(m)存在下面关系：

P()m(m)ejm (3-2)

1(m)2P()ejmd (3-3)

但是，实际中我们很难得到准确的自相关序列(m)，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。

目前，常用的线性谱估计方法有两种，即相关函数法和周期图方法，本实验对这两种方法分别予以讨论。

1． 自相关函数法

假设我们已知随机信号x(n)的M长的自相关序列{(m)}，利用自相关函数法可以得到x(n)的功率谱估计：

1(m)LmLmx(i)x(im)*i1

ˆ()PM1mM1ˆ(m)ejm (3-4)

利用窗函数，上式又可表达为

ˆ()PmWRMˆ(m)ejm(m) (3-5)

其中，

RWM(m)为矩形窗函数，定义为

1RWM(m)0mMmM (3-6)

Rˆ()WMP()P因此，实际上是真正功率谱与窗函数(m)傅立叶变换的卷积。矩形窗函数

不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣。为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为

0.540.46cosmW(m)M0HMmMmM (3-7)

这种窗函数可以有效的抑制旁瓣，但是，此时主瓣宽度增大，从而降低了谱估计的分辨率，这种主瓣和旁瓣之间的矛盾在线性谱估计方法中是无法解决的。

2． 周期图方法

假设已知随机信号x(n)的N个样本，利用周期图方法，信号x(n)的功率谱估计为

ˆ()1PNx(n)en0N12jn (3-8)

2利用上述方法得到的谱估计方差与信号x(n)的功率谱平方P()成正比，为了减小它的

方差，可以将信号序列进行分段处理，然后再求各分段结果的平均，这就是平均周期图方法，即Bartlett方法。

(1) Bartlett 平均周期图方法

将一个随机序列x(n) (0≤n≤N）分成K段，每段长度为L，各段之间互不重迭，因而N=KL，可以想到，第i段的信号序列可表示为

xi(n)x(n(i1)L) 0nL,i1,...,K (3-9)

对于每一段的周期图又可写成

2jniIiˆ()x(n)en0L1,i1,...,K (3-10)

于是，功率谱估计定义为

ˆ()1PKIˆ()ii1K (3-11)

因此，对于固定的记录长度来讲，分段数K增大可使谱估计的方差减小，但是由于L的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率也是矛盾的。

除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效应。为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以得到较平滑的谱估计，当然，相应的分辨率也有所下降。

(2) 平滑平均周期图方法

这时一种改进的Bartlett周期图方法，它特别适用于FFT直接计算功率谱估值。将长度为N的平稳随机信号序列x(n)分成K段，每段长度为L，即L=N/K。但这里在计算周期图之前，先用窗函数WL(n)给每段序列xi(n)加权，K个修正的周期图定义为

1Iiw()LUx(n)W(n)eiLn0L12jn,i1,...,K (3-12)

其中U表示窗函数序列WL(n)的能量，

1LUWL2(n)Ln0 (3-13)

在这种情况下，功率谱估计可按下面表达式给出：

ˆ()1PKIi1Kwi() (3-14)

本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦信号进行谱估计：

x(n)aiej(ini)w(n)i1NS (3-15)

其中NS 为正弦个数，i，i和ai分别为第i个正弦信号的数字频率、相位和幅度，i随机的分布在（0，2）之间，w(n)为零均值方差等于

2w的复高斯白噪声。

三、实验内容和步骤

1． 仔细阅读有关线性谱估计的内容，根据给出的框图编制自相关函数法谱估计的程序。运行程序，输入

2N100,M10,10.6,10,a11,w0,选择矩形窗。观察谱峰位置

是否正确（注意：由于窗效应可能引起谱估计的非正定）。

图表3.1

2． 观察并记录参数变化对谱估计性能的影响。

（1）改变M=5，其它输入同步骤1，观察功率谱估计的主瓣宽度和旁瓣大小随自相关序列长度的变化情况。

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

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

2（4）改变w1，其它输入同步骤1，观察信噪比变化对谱估计性能的影响。

（5）改变N=10，

2w1，其它输入同步骤1，结合(4)的内容，观察数据长度及信噪比

对谱估计性能的影响。

2N100,NS2,0.6,0.8,0,aa11212123．运行程序，输入，w0，选择

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

4．运行自相关函数法谱估计程序，输入

2N100,M10,NS0,w1，选择矩形窗，观

察利用自相关函数法得到的白噪声信号谱估计。改变M=3，20，观察M的变化对白噪声谱估计的影响。

5. 根据框图，编制周期图法谱估计程序。自程序FFT可以直接调用Matlab中的函数。运行程序，输入N100,NF128,K10。选择窗函数为矩形窗NS1，

2w10.6,10,a10,w0，观察谱峰位置是否正确，并与步骤1结果比较。

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

四. 实验结果及其分析 1. 运行程序，输入

2N100,M10,10.6,10,a11,w0,选择矩形窗。得到的实验

仿真结果如图4.1所示。

Rectangle Window2018161412108642000.10.20.30.40.50.60.70.80.91

图4.1 M=10时，矩形窗的谱估计

由图4.1的仿真结果看出，谱峰的位置处于0.6pi处，说明估计正确。

2．观察并记录参数变化对谱估计性能的影响

（1）M=5时，其它条件不变，其仿真结果如图4.2所示。

Rectangle Window987654321000.10.20.30.40.50.60.70.80.91

图4.2 M=5时，矩形窗的谱估计

比较M=10和M=5的仿真结果，可以看出M越小，则主瓣宽度越大，分辨率越低，谱峰高度也降低了；旁瓣越少，旁瓣最大幅度也减小；主瓣与旁瓣的幅值比、宽度比基本没有改变。

（2）M=10，选用hamming窗，其仿真结果如图4.3所示

Hamming Window10987654321000.10.20.30.40.50.60.70.80.91

图4.3 M=10，选用hamming窗的谱估计

与图4.1中矩形窗的谱估计相比，图4.3中我们选用的Hamming窗，其功率谱的主瓣宽度增加了近一倍，旁瓣相应也减少了许多，同时发现Hamming窗的分辨率降低了。

（3）1/4，其它输入同步骤1，仿真结果如图4.4所示

Rectangle Window2018161412108642000.10.20.30.40.50.60.70.80.91

图4.4 1/4，矩形窗的谱估计

谱估计的峰值位置仍处于0.6处，说明相位变化对谱估计的性能没有影响。

（4）改变

2w1，其它输入同步骤1，仿真结果如图4.5所示

Rectangle Window2018161412108642000.10.20.30.40.50.60.70.80.91

2图4.5 w1，矩形窗的谱估计

当增大高斯白噪声的方差时，对谱估计的峰值没有产生影响，但是影响了旁瓣的宽度和幅值。改变信噪比对频谱的影响较小。

（5）改变N=10，

2w1，其它输入同步骤1，其仿真结果如图4.6所示

Rectangle Window252015105000.10.20.30.40.50.60.70.80.91

图4.6 N=10，

2w1，矩形窗的谱估计

结合图4.1和图4.6，当N较大时，改变信噪比对谱估计影响很小；当N较小时，改变信噪比对谱估计的影响较大。图4.6中，频谱峰值偏离了0.6pi，同时旁瓣的幅值增大了许多。

2N100,NS2,0.6,0.8,0,aa1w1212123．运行程序，输入，0，加矩

形窗，调整M，其仿真结果如图4.7、4.8所示；加Hamming窗时分别如图4.9、4.10所示

Rectangle Window12108642000.10.20.30.40.50.60.70.80.91

图4.7 M=7时，矩形窗的谱估计

Rectangle Window12108642000.10.20.30.40.50.60.70.80.91

图4.8 M=8时，矩形窗的谱估计

由图4.7和图4.8可以看出，加矩形窗时，当M=8可将两个正弦频率分量临界分辨出

来。

Hamming Window10987654321000.10.20.30.40.50.60.70.80.91

图4.9 M=9时，Hamming窗的谱估计

Hamming Window12108642000.10.20.30.40.50.60.70.80.91

图4.10 M=10时，Hamming窗的谱估计

由图4.9和图4.10可以看出，加Hamming窗时，当M=10可将两个正弦频率分量临界分辨出来。

4．加矩形窗，4.13所示

2N100,NS0,w1，当M=10、3、20，仿真结果分别如图4.11、4.12、

Rectangle Window1.41.210.80.60.40.2000.10.20.30.40.50.60.70.80.91 图4.11 M=3时，加矩形窗的白噪声谱估计

Rectangle Window21.81.61.41.210.80.60.40.2000.10.20.30.40.50.60.70.80.91

图4.12 M=10时，加矩形窗的白噪声谱估计

Rectangle Window2.521.510.5000.10.20.30.40.50.60.70.80.91 图4.13 M=20时，加矩形窗的白噪声谱估计

图4.11、4.12、4.13对比发现：M越大，谱估计的各个分量变得越清晰，也就越容易分辨了。

5. 运行周期图法的谱估计程序，输入N100,NF128,K10，仿真结果如图4.14所示

Barlett + Rectangle0.40.350.30.250.20.150.10.05000.10.20.30.40.50.60.70.80.91

图4.14 N/K=10时的周期图法谱估计

对比步骤1所得的图4.1，图4.14中的谱峰位置同样也出现在0.6pi处，但主瓣幅值明显降低，主瓣宽度增加，另外旁瓣数量减少。

6．利用周期图方法重复步骤2(1)-(4)、3、4的内容，依次得到如下的仿真结果

（1）重复步骤2（1），选矩形窗可得如图4.15的周期图法谱估计

Barlett + Rectangle0.10.090.080.070.060.050.040.030.020.01000.10.20.30.40.50.60.70.80.91

图4.15 N/K=5时，选择矩形窗的周期图法谱估计

对比图4.14，N/K减小，周期图法谱估计的主瓣宽度越大，分辨率降低，幅值降低，而且旁瓣数量减少。这和自相关函数法中减小M时，谱估计出现的规律是一致的。

（2）重复步骤2（2），选Hamming窗可得如图4.16的周期图法谱估计

Bartlett + Hamming0.070.060.050.040.030.020.01000.10.20.30.40.50.60.70.80.91 图4.16 N/K=5时，选择Hamming窗的周期图法谱估计

对比图4.14，与矩形窗相比，选择Hamming窗，主瓣宽度加倍，旁瓣减小许多。周期图法谱估计和自相关法谱估计出现的规律相似。

（3）重复步骤2（3）得如图4.17的周期图法谱估计

Barlett + Rectangle0.40.350.30.250.20.150.10.05000.10.20.30.40.50.60.70.80.91 图4.17 取1/4，N/K=10时的周期图法谱估计

对比图4.14，取1/4和10时，两者的周期图法谱估计没有发生变化，所以初始相位的变化对谱估计的性能没有影响。

（4）重复步骤2（4）得如图4.18的周期图法谱估计

Barlett + Rectangle0.40.350.30.250.20.150.10.05000.10.20.30.40.50.60.70.80.91 图4.18 取

2w1，N/K=10时的周期图法谱估计

对比图4.14，当增大高斯白噪声的方差时，对谱估计的峰值没有产生影响，但是影响了旁瓣的宽度和幅值。改变信噪比对频谱的影响较小。

（5）用周期图法重复步骤3，运行程序，输入10.6,20.8,120,a1a21，

2w0，加矩形窗，调整段数K，其仿真结果如图4.19、4.20所示；加Hamming窗时，

分别如图4.21、4.22所示

Barlett + Rectangle0.080.070.060.050.040.030.020.01000.10.20.30.40.50.60.70.80.91 图4.19 K=20(N/K=5)时，加矩形窗的周期图法谱估计

-179876543210x 10Barlett + Rectangle00.10.20.30.40.50.60.70.80.91

图4.20 K=25(N/K=4)时，加矩形窗的周期图法谱估计

对比图4.19和图4.20，可以看出，加矩形窗时，当K=20时可将两个正弦频率分量临界分辨出来。我们仔细观察发现，图4.29中，谱峰的位置已不再是0.6pi和0.8pi了；在发现图4.20中，谱估计也已经发生了比较严重的错误。

Bartlett + Hamming0.040.0350.030.0250.020.0150.010.005000.10.20.30.40.50.60.70.80.91

图4.21 K=20(N/K=5)时，加Hamming窗的周期图法谱估计

6x 10-17Bartlett + Hamming54321000.10.20.30.40.50.60.70.80.91

图4.22 K=25(N/K=4)时，加Hamming窗的周期图法谱估计

对比图4.21和图4.22，可以看出，加矩形窗时，当K=20时可将两个正弦频率分量临界分辨出来。同样也出现了图4.19和图4.20中的现象。

2N100,NS0,w1（6）利用周期图方法重复步骤4，加矩形窗，20，仿真结果分别如图4.23、4.24、4.25所示

，当N/K=10、4、

Barlett + Rectangle0.0140.0120.010.0080.0060.0040.002000.10.20.30.40.50.60.70.80.91 图4.23 N/K=4时，加矩形窗的白噪声谱估计

Barlett + Rectangle0.060.050.040.030.020.01000.10.20.30.40.50.60.70.80.91 图4.24 N/K=10时，加矩形窗的白噪声谱估计

Barlett + Rectangle0.080.070.060.050.040.030.020.01000.10.20.30.40.50.60.70.80.91 图4.25 N/K=20时，加矩形窗的白噪声谱估计

图4.23、4.24、4.25对比发现：N/K越大，谱估计的各个分量变得越清晰，也就越容易分辨了。

综上，对比自相关函数法的谱估计，由步骤5、6的仿真结果，我们发现周期图法谱估计减小了峰值的幅度，加宽主瓣宽度，但是周期图使方差减小，从而得到一致的谱估计。周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计。

五． 实验主要结论

通过仿真，我们得到以下实验结论：

1． 用矩形窗进行谱估计时，主瓣宽度小，但是旁瓣的幅值较大。若条件不变，我们选用矩形窗时，主瓣宽度增大，分辨率降低，但是良好的抑制了旁瓣。

2． 第三步中，分开两个频谱成分，矩形窗在M=8即可分开，而Hamming窗在M=10时才能刚好分开，这也验证了Hamming窗的主瓣宽度大、分辨率低这一事实。

3． 增大M(或N/K)值，会使主瓣宽度减小，即分辨率增加，但是主瓣和旁瓣的幅值比不会发生变化。

4． 初始相位对谱估计没有影响。

5． 当信号序列较长时，能够采样到全部(或者大部分)的频率成分，若增大信号噪声方差，对谱估计的影响是很小的；而当信号序列较小时，只能采到少量的频率成分，对频谱估计的影响较大，甚至无法得到正确的谱估计。

6． 对于白噪声的谱估计，M(或N/K)值越大，则谱估计的越准确。这是由于增大了信号量，则会增强频谱的分辨率。

六． 思考题

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

ˆ(m)ejwm]ˆ(0)ˆ(w)2Real[P0M1，其中

证明： 先将(3-14)改写为三项相加的形式

ˆ(w)PM1M1ˆ(m)ejwmM10ˆ(m)ejwmM10ˆ(m)eM10jwmˆ(0)

ˆ(m)e0M1jwmˆ(m)ejwmˆ(0)

由于

ˆ(m)1L|m|*x(i)x(im)L|m|i1，

L|m|1ˆ(m)ˆ(m)x(i)x*(im)L|m|i1对其求共轭

*M1于是

ˆ(m)e0jwmˆ(m)e*0M1jwmˆ(m)ejwm)*(0M1

所以得

ˆ(m)eˆ(w)P0M1jwmˆ(m)e(0M1jwm*ˆ(0)2Real[)ˆ(m)ejwm]ˆ(0)0M1

2 .已知实信号x(n)的N 个样本x(0),x(1),…,x (N-1),可以定义PA()的估计：

ˆˆ(w)PAm(N1)N1ˆ(m)ejwm

ˆ(m)1N其中

N|m|1i0x(i)x*(im)

ˆ(w)1PAN试推导：

x(n)ei0N12jwm

解：

ˆ(w)PAm(N1)N1ˆ(m)ejwm

1m(N1)NN1Nm1n0x(n)x*(nm)ejwm

1m0NN1x(n)x(nm)e*n0N1jwm

N1N1*jw(nm)1jwnx(nm)ex(i)em0Nn0

令 lnm

N1N11*jwlˆ(w)x(l)ex(n)ejwnPAm0Ni0* N1x(l)ejwlN11m0Nx(n)ejwni0 N121jwnNx(n)en0

附录：

1. 自相关函数法谱估计程序

clear all; close all; clc;

N = 100;

M = 20;

NS = 2;

omega1 = [0.6*pi 0.8*pi];

phi1 = [0 0];



a1=[1 1];

sigma_2 = 1;

L = 256;

%-------------------------------------------

%%式(3-15)产生的信号

xxn = zeros(N,1);

for k = 1 : NS

xn1 = a1(k)*exp(j*(omega1(k)*(1:N)'+phi1(k)));

xxn = xxn + xn1;

end

wn = sqrt(sigma_2).*(randn(N,1));

xn = xxn + wn;

%-------------------------------------------

%%Rectangle Window

phi_x = xcorr(xn,'unbiased');

Lr = 2*M - 1;

Rectangle = window(@rectwin, Lr);

phi_xm = phi_x(N-M+1 : N+M-1).*Rectangle;

pw = fft(phi_xm,L);

pw_a = abs(pw(1:L/2));

%-------------------------------------------

%% Hamming Window

Lh = 2*M - 1;

Hammwin = window(@hamming, Lh);

phi_xh = phi_x(N-M+1 : N+M-1).*Hammwin;

pw_hamm = fft(phi_xh, L);

pw_hamm_a = abs(pw_hamm(1:L/2));

%-------------------------------------------

%%绘图

figure(1);

stem((0:L/2-1)/(L/2),pw_a,'.');grid on;

title('Rectangle Window');

figure(2);

stem((0:L/2-1)/(L/2),pw_hamm_a,'.');grid on;

title('Hamming Window');

2. 周期图法谱估计程序

clear all; close all; clc;

N = 100;

NS = 0;

omega1 = [0.6*pi 0.8*pi];

phi1 = [0 0];

a1=[1 1];

sigma_2 = 1;

K = 5; %分成K段

Lk = N/K; %每段的长度

NF = 128;

%-------------------------------------------

%%式(3-15)产生的信号

xxn = zeros(N,1);

for k = 1 : NS

xn1 = a1(k)*exp(j*(omega1(k)*(1:N)'+phi1(k)));

xxn = xxn + xn1;

end

wn = sqrt(sigma_2).*(randn(N,1));

xn = xxn + wn;

%-------------------------------------------

%%% Bartlett平均周期图方法

%%加矩形窗

Rectangle = window(@rectwin, Lk);

U = sum(Rectangle.^2)/Lk; %计算窗函数序列的能量

for i = 1:K

xinR = xn((i-1)*Lk+1 : i*Lk , :).*Rectangle;

xinR_dft = fft(xinR,NF); %用FFT计算NF点功率

IiwR(i,:) = xinR_dft.^2/NF/U;

end

pwR_peri = abs(sum(IiwR)/K); %功率谱估计

%-------------------------------------------

%%加Hamming窗

Hamming = window(@hamming, Lk);

U = sum(Hamming.^2)/Lk;

for i = 1:K

xinH = xn((i-1)*Lk+1 : i*Lk , :).*Hamming;

xinH_dft = fft(xinH,NF);

IiwH(i,:) = xinH_dft.^2/NF/U;

end

pwH_peri = abs(sum(IiwH)/K);

%-------------------------------------------

%%绘图

figure(1);

stem((0:NF/2-1)/(NF/2-1),pwR_peri(1:NF/2),'.');grid on;

title('Barlett + Rectangle');

figure(2);

stem((0:NF/2-1)/(NF/2-1),pwH_peri(1:NF/2),'.'); grid on;

title('Bartlett + Hamming');