第3章 功率谱估计和信号频率估计方法

− jω m
2 1 = U N (ω ) N
而当 M
ˆ ( m) N − 1 时，这相当于对长度为 2 N − 1 的 r
做截断处理，也即施加了一个矩形窗，即
ˆM (m ) = w2 M +1 (m) r ˆ (m) r
UESTC 何子述，夏威
(R )
所以，BT法实际上是对周期图法的平滑。
12
3.1.3 经典功率谱估计性能讨论
ˆ ( m) 另外，还有一种常用的 r (m) 的估计 r
1 ˆ ( m) = r N− m
∗ u n u ∑ N ( ) N (n − m), n=0
N −1
m ≤ N −1
其均值为
ˆ (m)} = r (m) E {r
UESTC 何子述，夏威 9
ˆ (m)的方 若信号 u (n)是零均值的实高斯随机信号， 则r
典的功率谱估计应该是有偏的。但是，当 N → ∞ ，
W2(N)−1 (ω )趋向于冲激函数， 因此该估计又是渐近无偏的。
T
⑵ 方差
2 σ 假定 u (n)是零均值的实高斯白噪声，方差为 u ，
ˆ (ω ) 和 S ˆ (ω ) 的协方差为 S 1 2
2 2⎫ ⎧ ⎪ ⎪ ⎡ ⎤ ⎡ ⎤ + − N ω ω N ω ω ( ) ( ) ⎪ ⎪ 1 2 1 2 ⎢ sin ⎥ ⎢ sin ⎥ ⎪ 4 ⎪ ⎪ σu ⎪ ⎢ ⎥ ⎢ ⎥ ⎪ ⎪ 2 2 ˆ ˆ cov S (ω1 ), S (ω2 ) = 2 ⎨⎢ ⎥ +⎢ ⎥ ⎬ N ⎪ ⎢ sin ω1 + ω2 ⎥ ⎢ sin ω1 − ω2 ⎥ ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ ⎥ ⎢ ⎥ ⎪ ⎪ 2 2 ⎢ ⎥ ⎢ ⎥ UESTC 何子述，夏威 14 ⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭
1 ˆ (m) = r N
UESTC 何子述，夏威
∗ u n u ( ) ∑ N N (n − m), n=0
N −1
m ≤ N −1
(1)
根据维纳-辛钦（Wiener-Khintchine）定理，对上式
3
求傅里叶变换，即
ˆ (ω ) = S BT
m =− M
∑
M
ˆ (m) e− jω m , r
PER
}= 1 ∫ π N 2 ω ( )}
⎡W (ω )⎤ d ω = 1 ⎦ −π ⎣ N
π 2
m=− M
∑
M
w2 ( m ) < 1
ˆ (ω ) 的方差小于 S ˆ (ω ) 的方差，这正是W ω S 这说明， BT ( ) PER ˆ (ω ) 平滑的结果。 对S PER ˆ (ω ) 谱的平滑（即方差减小） •由以上讨论可知， S BT 是以牺牲分辨率为代价的。由于 W (ω )主瓣比三角窗 的主瓣宽，因而使其分辨率下降。谱的平滑同时也
第3章 功率谱估计和信号频率估计方法
UESTC 何子述，夏威
1
•本章要回答的问题是，怎样利用随机过程u (n) 的N 个观测数据 uN (0), uN (1), , uN ( N −1)估计出随机过程 的功率谱 S (ω ) ？ •经典功率谱估计 •参数模型法估计 •基于相关矩阵特征分解的信号频率估计
差为
1 ˆ var {r (m)} = N− m ⎡ ⎤ l 2 ⎢ ⎥ ⎡ ⎤ 1 − + + − r l r l m r l m ( ) ( ) ( ) ∑ ⎢ ⎥⎦ ⎢ ⎥⎣ N m − l =−( N −1−| m|) ⎢ ⎥⎦ ⎣
N −1−| m|
ˆ (m) 为无偏估计， 当 m 接近于 N 时， 由以上两式得， r
SPER (ω ) 的均值为
L
∑
L
i =1 n= 0
i − jω n u n e ( ) ∑ N
M −1
2
1 (T ) W E {SPER (ω )} = 2 M −1 (ω − λ ) S (λ ) d λ ∫ 2 π 2π
SPER (ω ) 的方差为
1 var {S PER (ω )} = 2 L
UESTC 何子述，夏威
2
3.1 经典功率谱估计方法
•经典功率谱估计是基于传统傅里叶变换的思想，其 中的典型代表有Blackman和Tukey提出的自相关谱估 计（简称为BT法），和周期图法。
3.1.1 BT法
1.自相关函数的估计与傅里叶变换 设 u N (n) 为 u (n)的 N 个观测值，则 u (n) 的自相关函数 估计值
M ≤ N −1
(2)
以此结果作为对真实功率谱 S (ω ) 的估计，因为式(2)是 通过自相关函数间接得到的，称为间接法。
ˆ (m) 的运算量就很大。这时 r 如果m 和 N 都比较大， ˆ (m) 进行计算。 可以采用FFT对 r
对式（1）求傅里叶变换，并整理得
UESTC 何子述，夏威 4
m=−( N −1)
( ) w2 的乘积， 。 N −1 ( m) 的长度为 2 N − 1
T
(2) 方差 ˆ(m) 的方差为 r
ˆ (m)} = E r ˆ(m) − E {r ˆ(m)} var {r
2
{ } ˆ (m) } − E {r ˆ(m)} = E{ r
2
2
假定信号 u (n) 是零均值的实高斯随机信号，得
{
}
{
}
1 π W (ω ) d ω = w (0) = 1 ∫ 2π − π
约束下，它是渐近无偏估计。不过由于W (ω ) 的影响， 其偏差趋于零的速度要小于周期图法，因此对周 期图作平滑的结果是使偏差变大。
UESTC 何子述，夏威 17
⑵ 方差
Kr
{ = ˆ var {S
ˆ (ω ) var S BT
导致估计的偏差变大。由此可以看出，在方差，偏
UESTC 何子述，夏威
差和分辨率之间存在着矛盾，在实际应用中，只能 18
根据需要做出折衷的选择。
3.1.4 经典功率谱估计的改进
周期图法估计出的谱性能不好，当观测数据长度 太大时，谱的曲线起伏加剧；而数据太短时，谱的分 辨率又不好。因此需要加以改进。 1 Bartlett法 Bartlett法的基本步骤是：将 N 点的观测数据 u N (n) 分为 L 段，每段的长度为 M 即
(4)
其中， U N (ω ) = ∑ u N (n) e− jωn
n=0
N −1
因为这种功率谱估计方法是直接通过观察数据的
UESTC 何子述，夏威 11
傅里叶变换求得的，所以人们习惯上称之为直接法。 当 M = N − 1 时，周期图法和BT法是相同的，即
m=−( N −1)
∑
N −1
ˆ ( m) e r
UESTC 何子述，夏威 5
步骤2 求 U 2 N (k )的幅度平方，然后除以 N ，得 N 步骤3
2 1 ˆ0 (m) U k 对 N 2 N ( ) 进行IFFT，得r
1
U 2 N (k )
2
m = 0,1,
, 2 N −1
ˆ (m) 的关系为 ˆ0 (m) 与 r r
⎧ ˆ0 (m) ⎪ r ⎪ ˆ (m) = ⎨ r ⎪r ⎪ ⎩ ˆ0 (m + 2 N )
{
}
若取 ω1 = 2k π
N
和 ω2 = 2l π N ，则上式变为
2 2⎫ ⎧ ⎪ ⎪ ⎡ ⎤ ⎡ ⎤ ⎪ ⎪ sin k l sin k l + − π π ( ) ( ) 4 ⎪⎢ ⎪ ⎥ ⎢ ⎥ ˆ ˆ cov S (ω1 ) , S (ω2 ) = σu ⎨⎢ + ⎬ ⎥ ⎢ N sin k − l π N ⎥ ⎪ ⎪ sin + N k l N π ( ) ( ) (( ) )⎥⎦ ⎪ ⎢⎣ ⎥⎦ ⎢⎣ ⎪ ⎪ ⎪ ⎩ ⎭
ˆ (m) 的方差很大，但当 N 估计 r
ˆ ( m) 是 r ( m) m 时， r
的渐近一致估计。
UESTC 何子述，夏威 10
3.1.2 周期图法
ˆ (ω ) 周期图（Periodogram）法又称直接法。以 S PER
表示周期图法估计出的功率谱，则
2 1 ˆ S PER (ω ) = U N (ω ) N
1 M = N − 1时的估计性能 在这种情况下，周期图法和BT法的性能是一致的。 ⑴ 均值 BT法的均值为
1 (T ) ˆ S (ω ) ∗ W2 N −1 (ω ) E SBT (ω ) = 2π
{
}
由上式可知，功率谱估计的均值可以表示为信号的
UESTC 何子述，夏威 13
T) 真实功率谱 S (ω )和窗函数 W2(N −1 (ω )的卷积，因此，经
UESTC 何子述，夏威
1 i ˆ ˆ i (ω ) S ω = S var var ( ) ∑ PER PER L i =1
L
{
}
{
}
ቤተ መጻሕፍቲ ባይዱ
21
Bartlett功率谱估计频率分辨率下降 由以上讨论可知， 为原来的 1 L ，方差也减小为周期图法的 1 L ，因此 Bartlett功率谱估计较周期图法的结果更为平滑。 2 Welch法 这种方法是Welch在1967年提出的，又称修正平均 周期图法，是应用较广的一种方法。它是对Bartlett法 的改进。
{
}
当 N 增大时，会使互不相关的点增多，这就加剧了 估计的功率谱曲线的起伏。 若取 ω1 = ω 2 = ω ，功率谱估计的方差为
⎡ ⎛ sin N ω ⎞2 ⎤ ⎥ ˆ (ω ) = σ 4 ⎢1 + ⎜ ⎟ var S ⎟ ⎜ u ⎢ ⎟⎥ ⎜ ⎝ ⎠ ω N sin ⎢⎣ ⎥⎦
{
}
UESTC 何子述，夏威