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

数字信号处理II
——随机信号的功率谱估计方法
一、实验目的
1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。

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

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

4．体会非参数化功率谱估计方法的优缺点。

二、实验原理与方法
假设信号()x n 为平稳随机过程，其自相关序列定义为:
{}*()()()m E x n x n m φ+@
(0.1)
其中{}E •表示取数学期望，{}*
•表示取共轭。

根据定义，()x n 的功率谱密度()P w 与自相关序列()m φ存在如下关系：
()()j m
m P m e
ωωφ+∞
-=-∞
=
∑ (0.2)
1
()()2j m
m P e
d π
ωπφωωπ
-
=
⎰
(0.3)
然而，实际中我们很难得到准确的自相关序列()m φ，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。

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

1．自相关函数法
假设已知随机信号()x n 的N 个观测样本，则其自相关序列可以用下式进行估计:
||1
*0
1
ˆ()()()||1||N m n m x n x n m m N N m φ
--==+≤--∑ (0.4)
当仅使用长度为2M -1的自相关序列时，对其进行傅立叶变换即可得到功率谱估计如下:
11
ˆˆ()()M j m m M P
m e ωωφ
--=-+=∑
(0.5)
其中M 为加窗长度，Re ()c
M
W m 为矩形窗函数，定义如下:
Re 1,||()0,||c
M
m M W
m m M <⎧=⎨≥⎩
(0.6)
因此， ˆ()P
w 在一定程度上可以看作是“真正的功率谱()P w ”与窗函数傅立叶变换的卷积。

矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣，旁瓣效应使那些处于旁瓣附近功率较小的频率分量被淹没掉。

为了降低旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为：
0.540.46cos ,||()0,
||Ham
M m m M W m M
m M
π⎧
+<⎪=⎨⎪≥⎩ (0.7)
这种窗函数可以有效的抑制旁瓣，但此时主瓣宽度增大，从而降低了谱估计的分辨率。

这种主瓣和旁瓣之间的矛盾在非参数化功率谱估计方法中是无法解决的。

2．周期图方法
假设已知随机信号()x n 的N 个样本，利用周期图方法，信号()x n 的功率谱估计为:
1
20
1ˆ()|()|N j m n P
x n e N ωω--==∑ (0.8)
利用上述方法得到的谱估计方差与信号()x n 的功率谱平方2[()]P w 成正比，为了减小它的方差，可以将信号序列进行分段处理，然后再求各分段结果的平均，这就是平均周期图方法，即Bartlett 方法,它特别适用于FFT 直接计算功率谱估值。

（1）Bartlett 平均周期图方法
将一个随机序列()x n (0)n N ≤≤分成K 段，每段长度为L ，各段之间互不重叠，因而N LK =。

可以想到，第i 段的信号序列可表示为：
()((1))
0n<,1i x n x n i L L i K =+-≤≤≤ (0.9)
每一段的周期图又可写成：
1
2
1ˆ()|()|1L j m i i n I x n e i K L ωω--==≤≤∑
(0.10)
于是，功率谱估计定义为：
1
1ˆˆ()()K
i i P
I K
ωω==∑
(0.11)
对于固定的记录长度来讲，分段数K 增大可使谱估计的方差减小，但是由于L 的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率也是矛盾的。

除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效应。

为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以得到较平滑的谱估计，当然，相应的分辨率也有所下降。

（2）平滑平均周期图方法
这是一种改进的Bartlett 周期图方法，将长度为 N 的平稳随机信号序列()x n 分成K 段，每段长度为L ，即/L N K =。

这里在计算周期图之前，先用窗函数
L W (n)给每段序列(n)i x 加权，K 个修正的周期图定义为：
2-1
-=0
1ˆ()=(n)W (n)e
1L W j m
i
i
L
n I x i K LU
ωω≤≤∑ (0.12)
其中U 表示窗函数序列每个样本的平均能量，即：
2
-1
=0
1=
W (n)
L L
n U L
∑ (0.13)
在这种情况下，功率谱估计可按下面表达式给出：
=1
1ˆˆ
()=()K
W i
i P
I K
ωω∑ (0.14)
本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦
信号进行谱估计：
(n+)=1(n)=+(n)s
i i N j i i x a e w ωϕ∑
(0.15)
其中N s 为正弦个数，i ω，i ϕ和i a 分别为第i 个正弦信号的数字频率、相位
和幅度，i ϕ随机分布在 [0,2)π之间，(n)w 为零均值，方差等于2
w σ的复高斯白
噪声。

三、实验内容及步骤
1、仔细阅读有关功率谱估计的内容，按照图1给出的框图编制自相关函数法谱
估计程序。

运行程序，输入2
=100,=10,=1,=0.6,=0,=0,s i i w
N M N ωπϕσ选择矩形窗。

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

图 1 实验流程框图 图2 M=10时，矩形窗的谱估计
从图中可以看出谱峰位置在0.6*pi 处，这个位置是正确的。

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

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

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

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

（一般来讲，与矩形窗相比，加 Hamming 窗，主瓣宽度增加一倍，但旁瓣减小许多）。

图4 M=10时， hamming 窗的谱估计
与矩形窗相比，选用汉明窗，功率谱的主瓣宽度增加了近一倍，旁瓣也减少了许多，分辨率降低了。

(3) 改变1=/4ϕπ，其它输入同步骤 1，观察初始相位的变化对谱估计性能的影响。

图5 4/1πϕ=，矩形窗的谱估计
由图可知改变初始相位，对频谱估计的结果没有影响。

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

图6
2
1w σ=，矩形窗的谱估计
当噪声增加，信噪比减小时，会使谱估计得到的结果的精确性减小，而且增加的噪声越大，对谱估计影响就越大，
(5) 改变2
=10,=1w
N σ，其它输入同步骤 1，结合（4）的内容，观察数据长度及信噪比对谱估计性能的影响。

图7 N=10，2
1w σ=，矩形窗的谱估计
由图知，在N 较大时，改变信噪比对谱估计影响比较小；而当N 较小时，改变信噪比对谱估计的影响较大。

在上图中我们可以看出此时频谱峰值已偏离了0.6pi ，同时旁瓣的幅值增大了许多。

3、运行程序，输入 2
121212=100,=2,=0.6,=0.8,==0,==1,=0,s w
N N a a ωπωπϕϕσ选择矩形窗，调整自相关序列长度M ，使得两个正弦频率分量临界分辨出来，记录此时M 值，并绘出功率谱图。

同样，在加 Hamming 窗的情况下，记录使两个正弦频率分量临界分辨出来的M 值，并绘出此时的功率谱图。

图8 当M=7时，矩形窗的谱估计
图9 当M=8时，矩形窗的谱估计
由上面两图可知，在M=8时，加矩形窗可以将两个正弦频率分量临界分辨出来。

图10 当M=9时，汉明窗的谱估计
图11 当M=10时，汉明窗的谱估计
由上两图可知，在M=10时，加汉明窗可以将两个正弦频率临界分辨出来。

4、运行自相关函数法谱估计程序，输入2
=100,=0,=1,M=10,s w
N N σ选择矩形窗，观察利用自相关函数法得到白噪声信号谱估计。

改变M =3、20，观察M 的变化对白噪声谱估计的影响。

图12 当M=3时，加矩形窗的白噪声谱估计 图13 当M=10时，加矩形窗的白噪声谱估计 图14 当M=20时，加矩形窗的白噪声谱估计
由上面三幅图知，M 越大，加矩形窗的谱估计各个分量变得越清晰，更容易分辨。

5． 根据图1给出的框图，编制周期图法谱估计程序，子程序FFT 可以直接调用Matlab 中的函数。

运行程序，输入=100,=128,=10,F N N K 选择窗函数为矩
形窗，输入2
111=1,=0.6,=0,=1,=0,s w
N a ωπϕσ观察谱峰位置是否正确，并与步骤 1 结果比较。

图15 N/K=10时的周期图法谱估计
此图峰值在0.6*pi 处，位置正确，并且由对比可知，主瓣峰值明显下降了，而且其宽充增加，旁瓣数目有所减少。

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

（1）重复步骤2(1)，选矩形窗可得如下图所示的周期图法谱估计
图16 N/K=5时，矩形窗的Bartlett 周期图谱估计
对比每2(1)步中结果可知，用Bartlett 周期图谱估计的主瓣在N/K 减小时，主瓣宽度变大，分辨率降低，峰值下降，而且旁瓣数目减少，其规律和自相关函数法中减小M 时，规律一致。

（2）重复步骤2(2)，选Hamming 窗可得如下图所示的周期图法谱估计
图17 N/K=5时，汉明窗的Bartlett 周期图谱估计
对比第2(2)步中结果可知与矩形窗相比，汉明窗知主瓣宽度加倍，旁瓣减小许多。

Bartlett 周期图谱估计和自相关谱出现的规律相似。

（3）重复步骤2(3)得到如下图所示的周期图法谱估计
图18 4/1πϕ=，N/K=10时矩形窗的Bartlett 周期图谱估计
对比每2(3)可知，φ变化对谱估计没有太大影响。

（4）重复步骤2(4)得如下图所示的周期图法谱估计
图19 取12
=w σ，N/K=10时矩形窗的Bartlett 周期图谱估计
对比第2(4)步可知，增大高斯白噪声的方差时，对谱估计的峰值没有影响，但却使旁瓣增大。

（5）重复步骤2(5)得如下图所示的周期图法谱估计
图20 N=50，12
=w σ时 矩形窗的Bartlett 周期图谱估计
由图可知，当N 减小，同时增大噪声方差时，矩形窗的Bartlett 周期图谱估计的
主瓣宽度增加。

（5）用周期图法重复步骤3，输入100,2,N NS ==10.6,ωπ=20.8,
ωπ=12120,1a a φφ====，02
=w σ，加矩形窗，调整段数K ，其仿真结果如下图所示；
加Hamming 窗时，仿真如下图所示
图21 K=20时，矩形窗的Bartlett 周期图谱估计 图22 K=25时，矩形窗的Bartlett 周期图谱估计
由上两幅图可知，加矩形窗时，K=20时可以将两个正弦频率分量临界分辨出来，但峰值位置有所变化。

图23 K=20时，汉明窗的Bartlett 周期图谱估计 图24 K=25时，汉明窗的Bartlett 周期图谱估计
同样由上面两图可知，加汉明窗时，当K=20时，可以将两个正弦频率分量临界分辨出来。

同样其频率中心也发生了变化。

（6）利用周期图方法重复步骤4，加矩形窗，2
100,0,1w
N NS σ===，当
N/K=10、4、20，仿真结果分别如下图所示
图25 N/K=4时，矩形窗的Bartlett 周期图谱估计 图26 N/K=10时，矩形窗的Bartlett 周期图谱估计 图27 N/K=20时，矩形窗的Bartlett 周期图谱估计
由上图可知，N/K 越大，谱估计的各个分量变得越清晰，更容易分辨。

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

周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计。

四、思考题
1. 证明式（1.5）可以表示为：
-1
-=0
ˆˆˆ()=2Real[(m)e ]-(0),M j m m P
ωωφφ∑其中Real[]•表示取实部。

证明：10
11
1
ˆˆˆˆ()()()()M M j m j m j m
m M m M m P
m e m e m e
ωωωωφ
φ
φ-----=-+=-+===+∑
∑
∑Q
由自相关函数的对称性得：
即是： 10
ˆˆ()()(+)M j k
j m m P
m e e ω
ωωφ--==∑ˆ-(0)φ
即得： -1
-=0
ˆˆˆ()=2Real[(m)e ]-(0)M j m m P
ωωφφ∑ 2. 已知实信号(n)x 的 N 个样本(0),(1),...,(1),x x x N -，可以定义()A P ω的估计如下：
11
ˆˆ()(),N j m A
m N P m e ωωφ--=-+=∑其中||1
*0
1ˆ()()()N m i m x i x i m N φ--==+∑
试证明：21
1ˆ()()N j n
A
n P x n e
N
ωω--==∑。

证明：1||1
1
*1
1
1ˆˆ()()()()N m N N j m j m A
m N i m N P m e x i x i m e N ωωωφ------=-+==-+==+=
∑∑∑
Q 1||1
*()
1
1
[()]()N m N j i j m i i m N x i e x i m e N
ωω-----+==-++∑∑
211
*00
11[()]()()N N j n j n j n
n n x n e x n e x n e
N N
ωωω-----====∑∑
五、实验程序
1、自相关函数法（信号不为零）
clc;
clear all; n=100; m=10; ns=2;
w1=0.6*pi; w2=0.8*pi; fi1=0; fi2=0; a1=1; a2=1; sigma2=0; l=256;
sn1=zeros(n,1); sn2=zeros(n,1); xn=zeros();
sn1=a1*exp(i*(w1*(1:n)'+fi1)); sn2=a2*exp(i*(w2*(1:n)'+fi2)); wn=sqrt(sigma2).*(randn(n,1)); xn=sn1+sn2+wn;
xxn=xcorr(xn,'unbiased'); %产生信号 long=2*m -1; %矩形窗 Rectangle=window(@rectwin,long); p1=xxn(n -m+1:n+m -1).*Rectangle; p2=fft(p1,l);
p3=abs(p2(1:l/2));
long1=2*m -1; %汉明窗 Hammwin=window(@hamming,long1); p11=xxn(n -m+1:n+m -1).*Hammwin; p22=fft(p11,l);
p33=abs(p22(1:l/2)); figure(1);
stem((0:l/2-1)/(l/2),p3,'.'); title('矩形窗功率谱估计') grid on; figure(2)
stem((0:l/2-1)/(l/2),p33,'.'); title('汉明窗功率谱估计')
grid on;
2、自相关函数法（白噪声）
clear all;
n=100;
m=20;
w=0.6*pi;
sigma2=1;
l=256;
xn=zeros();
xn=sqrt(sigma2).*(randn(n,1));
xxn=xcorr(xn,'unbiased'); %产生信号long=2*m-1; %矩形窗Rectangle=window(@rectwin,long);
p1=xxn(n-m+1:n+m-1).*Rectangle;
p2=fft(p1,l);
p3=abs(p2(1:l/2));
long1=2*m-1; %汉明窗Hammwin=window(@hamming,long1);
p11=xxn(n-m+1:n+m-1).*Hammwin;
p22=fft(p11,l);
p33=abs(p22(1:l/2));
figure(1);
stem((0:l/2-1)/(l/2),p3,'.');
title('矩形窗功率谱估计')
grid on;
3、周期图法（信号不为零）
clc;
clear all;
n=50;
nf=128;
ns=2;
k=10; %所分段数w1=0.6*pi;
w2=0.8*pi;
fi1=0;
fi2=0;
a1=1;
a2=0;
sigma2=1;
l=n/k; %每段长度sn1=zeros(n,1);
sn2=zeros(n,1);
xn=zeros();
sn1=a1*exp(i*(w1*(1:n)'+fi1));
sn2=a2*exp(i*(w2*(1:n)'+fi2));
wn=sqrt(sigma2).*(randn(n,1));
xn=sn1+sn2+wn; %产生信号Rectangle=window(@rectwin,l); %Bartlett周期图方法u=sum(Rectangle.^2)/l; %计算窗函数序列的能量
for i=1:k
xnrectwin=xn((i-1)*l+1:i*l,:).*Rectangle;
xnrectwindft=fft(xnrectwin,nf);
ixnrec(i,:)=xnrectwindft.^2/nf/u;
end
pwrec=abs(sum(ixnrec)/k);
Hamming=window(@hamming,l); %加Hamming窗
u=sum(Hamming.^2)/l; %窗函数序列能量for i=1:k
xnhamm=xn((i-1)*l+1:i*l,:).*Hamming;
xnhammdft=fft(xnhamm,nf);
ixnhamm(i,:)=xnhammdft.^2/nf/u;
end
pwhamm=abs(sum(ixnhamm)/k);
figure(1)
stem((0:nf/2-1)/(nf/2-1),pwrec(1:nf/2),'.');
title('矩形窗Barlett周期图法')
grid on;
figure(2)
stem((0:nf/2-1)/(nf/2-1),pwhamm(1:nf/2),'.');
title('汉明窗Barlett周期图法')
grid on;
4、周期图法（白噪声）
clc;
clear all;
n=100;
nf=128;
k=5; %所分段数
sigma2=1;
l=n/k; %每段长度
xn=sqrt(sigma2).*(randn(n,1)); %产生信号
Rectangle=window(@rectwin,l); %Bartlett周期图方法
u=sum(Rectangle.^2)/l; %计算窗函数序列的能量for i=1:k
xnrectwin=xn((i-1)*l+1:i*l,:).*Rectangle;
xnrectwindft=fft(xnrectwin,nf);
ixnrec(i,:)=xnrectwindft.^2/nf/u;
end
pwrec=abs(sum(ixnrec)/k);
figure(1)
stem((0:nf/2-1)/(nf/2-1),pwrec(1:nf/2),'.'); title('矩形窗Barlett周期图法')
grid on;。