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

1 N
UN (w)2
26
归一化功率谱(dB) 归一化功率谱(dB)
0 -5 -10 -15 -20
-25 -30 -35 -40
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w 2p
（a） N = 32
0 -5 -10 -15 -20 -25 -30 -35 -40 -45 -50
当 M = N - 1 时，周期图法和BT法是相同的，即
åN- 1
rˆ(m)e-
m= - (N - 1)
jwm =
1 N
U N (w) 2
而当 M = N - 1时，这相当于对长度为 2N - 1的 rˆ(m)
做截断处理，也即施加了一个矩形窗，即
rˆM (m) = w2(RM)+ 1 (m)rˆ(m)
的渐近一致估计。
另外，还有一种常用的 r(m) 的估计 rˆ(m)
å rˆ (m) =
1 N- m
N- 1
uN (n)uN* (n -
n= 0
m),
其均值为
E {rˆ(m)}= r (m)
m? N 1
9
若信号 u(n)是零均值的实高斯随机信号，则 rˆ(m)的方
差为
å var {rˆ(m)}=
N
1 -
|m|
N,
| m |? N 1 其它
7
的乘积，w2(TN)- 1(m) 的长度为 2N - 1。 (2) 方差
rˆ(m) 的方差为
{ } var {rˆ(m)}= E rˆ(m) - E{rˆ(m)}2 { } = E rˆ(m) 2 - E{rˆ(m)}2
假定信号 u(n) 是零均值的实高斯随机信号，得
SNR1 = 30dB,SNR2 = 30dB,SNR3 = 27dB
25
步骤3 将三个实正弦信号和高斯白噪声进行叠加， 得观测信号 uN (n) 。
(1) 周期图法
步骤1 计算信号的离散傅立叶变换：
N- 1
å UN (w) = uN (n)e- jwn
n= 0
步骤2 计算信号的功率谱：
SˆPER (w)=
这种方法是Welch在1967年提出的，又称修正平均 周期图法，是应用较广的一种方法。它是对Bartlett法 的改进。
22
Welch法也对 N 点的信号 uN (n)进行分段，只是分段 时允许每段的信号有所交叠，通常取相邻两段的信号
交叠一半，若每段的信号长度仍为M ，信号被分为 L
段，则
L= N- M /2 M /2
将每段信号uNi (n)和窗函数w(n)相乘，然后按式(5)得到
每段信号的功率谱估计
å SˆPi ER (w) =
1 MU
M- 1
uNi (n)w(n)e-
n= 0
2 jwn
23
修正的周期图为
邋 S%PER (w) =
1 LMU
L i= 1
M- 1
2
uNi (n)w(n)e- jwn
n= 0
Welch方法允许分段数据样本的重叠，于是可以得到 更多的周期图估计，从而进一步减小估计的功率谱 密度的方差。通过窗函数加权，可以减小了相邻样 本段之间的相关性。所以，Welch方法可以更好地控 制功率谱密度估计的方差特性。
å SˆPi ER (w) =
1 M
M- 1
2
uNi (n)e- jwn ,
n= 0
1＃i L (5)
然后对每段功率谱估计结果作平均，得到平均周期图
20
邋 ? ( ) ( ) ( ) SPER
w
=
1 L
L i= 1
SˆPi ER
w
=
1 LM
L M- 1
u
i N
i= 1 n= 0
2
n e- jwn
SPER (w)的均值为
周期图（Periodogram）法又称直接法。以 SˆPER (w) 表示周期图法估计出的功率谱，则
SˆPER (w)=
1 N
UN (w) 2
(4)
N- 1
å 其中， UN (w)= uN (n)e- jwn
n= 0
因为这种功率谱估计方法是直接通过观察数据的
11
傅里叶变换求得的，所以人们习惯上称之为直接法。
所以，BT法实际上是对周期图法的平滑。
12
3.1.3 经典功率谱估计性能讨论
1 M = N - 1时的估计性能 在这种情况下，周期图法和BT法的性能是一致的。
⑴ 均值
BT法的均值为
{ } E
SˆBT (w)
=
1 2p
S (w)*W2(NT)-
1
(w)
由上式可知，功率谱估计的均值可以表示为信号的
13
是一对傅立叶变换。
算法3.1（用FFT计算自相关函数的方法）
步骤1 对 uN (n)补 N个零，得 u2N (n)，对 u2N (n)做快速 傅立叶变换（FFT）得U2N (k) k = 0,1,L ,2N - 1
5
步骤2
求 U2N
(k)的幅度平方，然后除以N，得
1 N
U2N (k) 2
步骤3
对
1 N
{ } ò E
SPER (w)
=
1 2p
( W (T)
2p 2M- 1
w-
l )S(l )dl
SPER (w) 的方差为
å { } { } { } var
SPER (w)
=
1 L2
L
var
i= 1
SˆPi ER (w)
=
1 var L
SˆPi ER (w)
21
由以上讨论可知，Bartlett功率谱估计频率分辨率下降 为原来的 1 L ，方差也减小为周期图法的 1 L ，因此 Bartlett功率谱估计较周期图法的结果更为平滑。 2 Welch法
第3章 功率谱估计和信号频率估计方法
1
•本章要回答的问题是，怎样利用随机过程u(n) 的N 个观测数据 uN (0),uN (1),L ,uN (N - 1)估计出随机过程 的功率谱 S (w) ？ •经典功率谱估计 •参数模型法估计 •基于相关矩阵特征分解的信号频率估计
2
3.1 经典功率谱估计方法
å var {rˆ(m)}=
1 N
N- 1- |m| l= - ( N- 1- |m|)
轾犏犏臌1-
|m|+ |l| N
轾犏臌r2 (l)+ r(l + m)r(l -
m)
8
由于r(m) 是有限的，显然当 N 时，rˆ(m) 的方差 将趋近于零。所以，对于固定的延时 m ，rˆ(m) 是 r(m)
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 w 2p
2p
p
W (w)dw
-p
ò 1
p
W (w)d w = w(0)= 1
2p - p
约束下，它是渐近无偏估计。不过由于W (w) 的影响，
其偏差趋于零的速度要小于周期图法，因此对周
期图作平滑的结果是使偏差变大。
17
⑵ 方差
{ } var SˆBT (w)
1
{ } ò å Kr = var SˆPER (w) = 2p N
是周期图所存在的固有矛盾。
2 M = N - 1时的估计性能
在这种情况下，两种方法不一致，BT法是对周 期图法的平滑。
16
⑴ 均值
{ } E SˆBT (w) = S (w)*W2(NT)- 1 (w)*W (w)
BT法也是一种有偏估计，当 N 很大，且在下面两式
{ } ò E
SˆBT (w)
=
S(w) 1
真实功率谱 S(w)和窗函数 W2(NT)- 1 (w)的卷积，因此，经
典的功率谱估计应该是有偏的。但是，当 N
，
W2(NT)- 1 (w)趋向于冲激函数，因此该估计又是渐近无偏的。
⑵ 方差
假定
u(n)是零均值的实高斯白噪声，方差为
s
2 u
，
Sˆ (w1 )和 Sˆ (w2 ) 的协方差为
{ } cov
导致估计的偏差变大。由此可以看出，在方差，偏
差和分辨率之间存在着矛盾，在实际应用中，只能 18
根据需要做出折衷的选择。
3.1.4 经典功率谱估计的改进
周期图法估计出的谱性能不好，当观测数据长度 太大时，谱的曲线起伏加剧；而数据太短时，谱的分 辨率又不好。因此需要加以改进。
1 Bartlett法
Bartlett法的基本步骤是：将N点的观测数据 uN (n) 分为 L段，每段的长度为 M 即
N
从式（3）可以看出， 对于固定的延时 m ，lim E{rˆ(m)}= r(m)。即 rˆ(m)是对
N
r(m) 的渐近无偏估计；
对于固定的 N，当 m 越接近于 N 时，估计的偏差
越大； 由式(3)可知，rˆ(m) 的均值是真值r(m)和三角窗函数
w2(TN)- 1(m) =
ìïïíïïî
10,
m
l
=
N - 1- (N-
|m| 1- |m|)
轾 犏 犏 犏 臌1-
l N- m
轾 犏 臌r2 (l)+ r (l + m)r (l -
m)
由以上两式得，rˆ(m) 为无偏估计，当 m 接近于 N 时，
估计 rˆ(m) 的方差很大，但当 N ? m 时，rˆ(m)是 r(m)
的渐近一致估计。
10
3.1.2 周期图法
19
L= N M
第 i(1＃i L)段数据加矩形窗后，变为
uNi (n) = u 轾 臌n + (i - 1)M wM(R) (n), 0 ＃n M - 1,1＃i L 其中，wM(R) (n)是长度为 M 的矩形窗。
对于每段数据
u
i N
(n
)，先利用周期图法求得其功率
谱 SˆPi ER (w) ，即
-
p p
轾 臌W
(w)
2
d
w
=
1 M w2 (m)< 1
N m= - M
这说明，SˆBT (w)的方差小于 SˆPER (w) 的方差，这正是W (w)
对SˆPER (w) 平滑的结果。
•由以上讨论可知，SˆBT (w) 谱的平滑（即方差减小）
是以牺牲分辨率为代价的。由于 W (w)主瓣比三角窗
的主瓣宽，因而使其分辨率下降。谱的平滑同时也
若取 w1 = w2 = w ，功率谱估计的方差为
{ } var
Sˆ (w)
=
s
4 u
轾 犏 犏 犏 臌1+
骣 ççç桫sNinsiNn ww ÷÷÷2
15
由上式可知，当 N
时，功率谱估计的方差不趋
近于零，而趋近于
s
4 u
，因此，经典功率谱估计不是一
致估计。
由上面的讨论可知，为了保证 Sˆ(w) 的渐近无偏性， 希望 N 要大，但是 N 增大时又使Sˆ(w)起伏加剧，这
U2N
(k) 2
进行IFFT，得rˆ0
(m)
m
=
0,1,L
,2N -
1
rˆ0 (m) 与 rˆ(m) 的关系为
rˆ(m) =
ìïïíïïî
rˆ0 (m) rˆ0 (m +
2N)
0 ＃m N - 1 - N + 1＃m - 1
2 自相关函数的估计性能 (1) 均值 rˆ(m)的均值为
6
E{rˆ(m)}= N- | m | r(m) (3)
如果m 和 N 都比较大，rˆ(m) 的运算量就很大。这时 可以采用FFT对 rˆ(m) 进行计算。
对式（1）求傅里叶变换，并整理得
4
åN- 1
rˆ(m)e- jwm =
m= - (N - 1)
1 N
UN
(w) 2
其中，UN (w)2 是 uN (n) 的能量谱，除以N 后即为功率
谱。这说明，由式(1)估计出的 rˆ(m) 和 uN (n) 的功率谱
•经典功率谱估计是基于传统傅里叶变换的思想，其 中的典型代表有Blackman和Tukey提出的自相关谱估 计（简称为BT法），和周期图法。
3.1.1 BT法
1.自相关函数的估计与傅里叶变换
设 uN (n) 为 u(n)的 N 个观测值，则 u(n)的自相关函数
估计值
å rˆ (m) =
1 N
N- 1
uN (n)uN* (n -
n= 0
m),
m ? N 1 (1)
根据维纳-辛钦（Wiener-Khintchine）定理，对上式 3
求傅里叶变换，即
M
å SˆBT (w)=
rˆ(m)e- jwm ,
m= - M
M ? N 1 (2)
以此结果作为对真实功率谱 S (w)的估计，因为式(2)是
通过自相关函数间接得到的，称为间接法。
Sˆ(w1), Sˆ(w2 )
=
s
4 u
N2
禳 镲 镲 镲 镲 镲 睚 镲 镲 镲 镲 镲 铪轾 犏 犏 犏 犏 犏 犏 臌sinsiNn
(w1 + w2
2 w1 + w2
2
)
2
+
轾 犏 犏 犏 犏 犏 犏 臌sinsiNn
(w1 -
2 w1 -
2
w2 w2
)
2
14
若取 w1 = 2kp N 和 w2 = 2lp N ，则上式变为
24
3.1.5 经典谱估计算法仿真实例及性能比较 例3.1 按如下步骤产生一组实验数据。 步骤1 产生长度为 N 的一段均值为0，方差为1的实高
斯白噪声序列 v(n)
步骤2 产生长度为 N 的三个实正弦信号 ，其归一化 频率分别是 f1 = 0.1, f2 = 0.25, f3 = 0.27 ， 信噪比分别为
{ } cov
ห้องสมุดไป่ตู้
Sˆ(w1), Sˆ(w2 )
=
s
4 u
禳 镲 镲 镲 睚 镲 镲 镲 铪轾 犏 犏 犏 臌N
sin(k +
sin((k +
l)p l)p
2
N) +
轾 犏 犏 犏 臌N
sin(k -
sin((k -
l)p l)p
2
N)
当 N 增大时，会使互不相关的点增多，这就加剧了
估计的功率谱曲线的起伏。