自功率谱估计的经典方法

5．自功率谱估计的经典方法 1） 周期图法（直接法）
对于时间序列)(n x N ，其傅里叶变换（DTFT ——离散时间信号的傅里叶变换）为
∑-=-=1
)()(N n n
j N j N e
n x e X ωω
，⎰-
=
π
π
ωωωπ
d e e X n x n j j N N )(21)(
记为
)()(ωj N D TFT
N e X n x −−→←
)(n x N 的离散傅里叶变换（DFT ）为
∑-=-=1
02)()(N n kn N
j N N e
n x k X π
，∑-==1
2)(1)(N k kn N
j N
N e k X N n x π
记为
)()(k X n x N D FT
N −−→←
若)(n x N 是信号)(n x 在时间域截断的结果，即
)()()(n d n x n x N N ⋅= (5-58)
其中，)(n d N 是单边矩形窗，其表达式为
⎩⎨
⎧-≤≤=其它
,01
0,1)(N n n d N 而)(n x 是确定性功率信号（或随机信号的一个样本序列），则根据第三章的讨论结果知，
=)(ωj x e S 2,)(1
)(lim
lim ωωj N N j x N N e X N
e P ∞
→∞
→= (5-59) 反映了信号)(n x 的平均功率在频域的分布情况，称为平均功率谱密度。

因此，估计量
2,,)(1)()(ˆωωωj N j x N j PER x e X N
e P e S == (5-60) 为信号)(n x 的功率谱的一个估计。

此估计方法称为直接法或周期图法。

在)(ˆ,ωj PER x e
S 的实际运算中采用DFT ，ω在单位园上均匀取值。

当取N
π
ω2=∆时，（5-60）改写为
2,,)(1)()(ˆk X N
k P k S N
x N N PER x ==，1,,1,0-=N k (5-61) 其中，
∑-=-=1
2)()(N n nk N
j
N N e
n x k X π，1,,1,0-=N k
当取N
22π
ω=
∆时，需对)(n x N 补N 个零后再作DFT ，此时（5-60）改写为 22,22,)(1)()(ˆk X N
k P k S N
x N N PER x ==，12,,1,0-=N k (5-62) 其中，)(2k X N 参见(5-42)、(5-33)式。

∑-=-=
1
20
2222)()(N n nk N
j
N
N e
n x
k X π，12,,1,0-=N k
在)(2k X N 的自变量取偶数的点，有
∑-=-=
1
20
22222)()2(N n nk N
j
N
N e
n x
k X π∑-=-=1
2)(N n nk N
j
N e
n x π，
因此，有
)2(2k X N )(k X N = (5-63-1)
)2(2,k S N PER x )(,k S N PER x = (5-63-2)
)12(2＋k X N 是)(k X N 和)1(＋k X N 之间的频率细化；)12(2,＋k S N PER x 是)(,k S N PER x 和)1(,+k S N
PER x 之间的频率细化。

)(2k X N 和)(2,k S N PER x 的频率分辨率分别是)(k X N 和)(,k S N PER x 的二倍。

2） 自相关法(间接法)
根据维纳—辛钦定理，自相关函数与自功率谱是一个傅里叶变换对。

因此，基于序列)(n x N 的功率谱估计可以表示为
∑-=-=
M
M
m m j x
N j M BT
x e m R
e S
ωω
)(ˆ)(ˆ,2,，1-≤N M (5-64-1)
在)(ˆ2,ωj M BT
x e S 的实际运算中采样DFT 。

当（5-63）式中取1-=N M 时，若ω在单位园上均匀取N 2个点，则(5-63）式可改写为
∑-+-=-=
1
1
22,2,)(ˆ)(ˆN N m km N
j
x
N N
BT
x e
m R
k S
π
，12,,1,0-=N k (5-64-2)
若ω在单位园上均匀取N 个点，则(5-65）式可改写为
∑-+-=-=
1
1
2,,)(ˆ)(ˆN N m km N
j
x
N N BT
x e
m R
k S
π
，1,,1,0-=N k (5-64-3)
显然，)(ˆ2,k S N BT x 和)(ˆ,k S N BT x 之间有关系： =)2(ˆ2,k S N BT x )(ˆ,k S N BT
x ，1,,1,0-=N k (5-65) )12(ˆ2,＋k S N BT x 是)(ˆ,k S N BT
x 和)1(ˆ,＋k S N BT
x 之间的频率细化。

3） 直接法与间接法的关系
当1-=N M 时，有
1）)(ˆ)(ˆ2,2,k S k S N PER
x N BT x =，12,,1,0-=N k (5-66-1) 2）)(ˆ)(ˆ,,k S k S N PER
x N BT x =，1,,1,0-=N k (5-66-2) 证明：1）根据（5-45），有
)(ˆ0,m R x
N −−→←DFT =)(ˆ2,k S N PER x )(,2k P x N 又
∑-+-=-=
1
122,2,)(ˆ)(ˆN N m km N
j
x
N N
BT
x e m R
k S
π
∑∑-=--+-=-+=1
22,1
1
22,)(ˆ)(ˆN m km N
j
x
N N m km N
j
x
N e m R e
m R
π
π
由km N
j e
22π-的周期性，有
km N
j
N m k N
j
e e
22)2(22ππ
---= 再根据（5-41）式（参考图.5－3），并考虑到0)(ˆ0,=N R x
N ，得 ∑-+-=-=
1
1
22,2,)(ˆ)(ˆN N m km N
j
x
N N BT
x e
m R
k S
π
∑∑-=---=-++=
1
0220,1
1220,)(ˆ)2(ˆN m km N
j
x
N N m km N
j x
N e m R e
m N R
ππ＋
∑∑-=--+=--+=
1
220,1
21
)2(220,)(ˆ)(ˆN m km N
j
x
N N N m N m k N
j
x
N e
m R
e
m R
π
π
∑-=-=1
20
220,)(ˆN m km N
j
x
N e
m R
π
因此，有
)(ˆ0,m R x
N −−→←DFT )(ˆ2,k S N BT x 综合上述，关系1）得证。

2）根据（5-63-2）、（5-65）和（5-66-1）式，有
＝)(,k S N PER x )2(2,k S N PER x =)2(ˆ2,k S N BT x ＝)(ˆ,k S N BT x
即(5-66-2)式成立。

4） 质量
当1-=N M 时，周期图法与相关法等价。

因此，我们直接利用相关法的估计公式，分析估计
质量，并将)(ˆ,ωj BT
x e S 简写为)(ˆωj x e S 。

估计均值：
{}
⎭
⎬⎫⎩⎨⎧=∑-+-=-1
1,)(ˆ)(ˆN N m m j N x j x e m R E e S E ωω
{}
∑-+-=-=
1
1,)(ˆN N m m
j N
x e m R E ω∑-+-=--=
1
1
)(|
|N N m m j x e m R N m N ω ∑-+-=-=1
1
)()(N N m m j x N e m R m ων)()(ω
ωj N j x e V e S *=φπππφωφd e V e S j j x ⋅⋅⎰--)()(21)(＝
【∑-+-=-11)()(N N m m
j x N e m R m ων∑⎰-+-=--=1
1
])(21)[(N N m m j m j j x N e d e e S m v ωππφφφπ
φπππφωφd e m v e S N N m m j N j x ⎰∑--+-=--⎥⎦
⎤⎢⎣⎡=1
1)()()(21φπππφωφd e V e S j N j x ⎰--=)()(21)(】 式中，
⎪⎩⎪⎨⎧≤-=其它,0||,||1)(N m N m m v N −−→←DTFT 2
2sin 2sin 1)(⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣
⎡=ωωN N e V j N
结论：1）)(ˆωj x
e S 是有偏估计。

从时域上看，是由于真实自相关函数)(m R x 被乘以三角窗函数)(m v N ；从频域上看，真实功率谱被三角窗函数谱)(ωj N e V 卷积。

2）)(ˆωj x e S 是渐近无偏估计，即{}
)()(ˆωωj x N j x
e S e S E −−→−∞→ 证明：因为
⎰-
=
π
π
ωωωπ
d e e V m v m j j N N )(21
)(，⎰-
=
π
π
ωωπ
d e V v j N N )(21)0(
当∞→N 时，
0)(0,−−−→−≠∞→ωωN j N e V ，N e V j N −−→−→0
)(ωω∞−−→−∞→N
1)0()(21−−→−∞→-⎰N N
j N v d e V ＝ππ
ω
ωπ 因此，
)()(ωδω−−→−∞
→N j N e V
{}
)()(*)()()()(ˆωωωωωωδj x j x N j N
i x j x e S e S e V e S e S E =−−→−*=∞→ [ 矩形窗与三角窗的关系：
)()()(m d m d m v N N N *-=，2|)(|1
)(ωωj N j N e D N
e V =
]
估计方差：
当序列)(n x 服从高斯分布且为白噪声序列时，可以求得功率谱估计的方差
])sin sin (1[)}(ˆ{2
42ˆω
ωσσωN N e S Cov x
j x S x
+==
其中x σ是序列)(n x 的标准差。

显然，当N 趋于无穷时，)}(ˆ{ωj x e S Cov 不趋于零。

因此，)(ˆωj x
e S 不是一致估计。

5） 存在的主要问题
)(ˆωj x e S 不是一致估计，其结果是)(ˆωj x
e S 的谱曲线起伏剧烈；并且随N 增大，起伏程度往往并不会减弱。

上述结论对一般的序列仍然成立。

我们也可以从另一个角度来说明这个问题：
测量信号中，不可避免地存在噪声。

由第三章的内容知，对于平稳的各态遍历的随机信号，可以通过时间域求平均获得信号的特征量，即
=)(m R x ∑-=∞→∞
→+=1
,)()(1)(ˆlim lim
N n N N x N m n x n x N m R
成立。

但是，这一结论在频域未必成立。

也就是
≠=∞→∞
→2,)(1)(ˆlim lim
ωωj N N j N PER x N e X N e S ⎭⎬⎫⎩⎨⎧=∞
→2)(1)(lim ωωj N N j x e X N E e S 6） 改进
以减少估计量的方差为主要目标。

1）分段平均
基本思路：若K 个随机变量K X X X ,,,21 相互独立并具有相同的均值μ和方差2
σ，则它们的平均
][1
21K X X X K
X +++=
的均值仍然是μ，但方差是K /2
σ，减小了K 倍。

具体做法：将数据)(n x N 分为互不重叠的K 段，各段长度为M (K N /=)，即
])1([)(M i n x n x N i M -+=，1,,1,0-=M n ，K i ,,2,1 =
对各段数据分别求功率谱)(ˆωj i x
e S
，即 ∑
-+-=-=1
1
,)(1)(ˆM M n jn i x M j i x
e
n R M
e S
ω
ω
，K i ,,2,1 = )(ˆωj i x e S 的平均就是)(ˆωj x
e S 的估计，即 )(ˆωj x e S ∑==K i j i x e S K 1
)(ˆ1ω
此时，对于高斯分布且为白噪声序列，估计方差为
])sin sin (
1[2
4
2ˆω
ωσσ
M M K
x
S x
+=
它随K 值增加而减小。

然而，估计均值为
＝)(ˆωj x
e S )()(ωωj M j x e V e S * 其中，三角窗函数谱)(ω
j M e
V 的主瓣宽度随K 的减小而增宽，进
而造成估计的偏差增大。

)(ωj e V
2）加数据窗平滑
由于数据截断，即加数据窗，不可避免。

因此，希望寻找到好的数据窗，它具有： 1） 3dB 带宽B 小；
2） 最大边瓣峰值A 小；
3） 边瓣谱峰渐近衰减速度D 大。

窗函数的主要性质：窗函数具有低通，有平滑作用。

3）分段平均与窗函数结合
6．互功率谱估计的经典方法 与自功率谱估计相似
)()(1)(ˆ22ωωωj N
j N j xy e Y e X N
e S *= 7．经典谱估计方法的特点
优点：可以利用FFT 算法，物理概念明确；
缺点：1）谱的分辨率较低，s s NT T /2/πω=∆=∆Ω
2）不可避免的窗函数，使真正谱)(ω
j e
S 在窗口主瓣内的功率向边瓣部分“泄漏”
，从而降低了分辨率；较大的边瓣又可能掩盖)(ω
j e S 中较弱成分，或产生假的峰值；尤其是，当N 较小时。

3）)(ˆωj e S
不是)(ωj e S 的一致估计，且随N 的增大，谱线的起伏加剧； 4）没有一个窗函数（方法）能使估计谱在方差、偏差和分辨率各方面都得到改善。