第四章 功率谱估计

1
2
-
WB (e j )Pxx (e j() )d
式中
Pxx (e j ) WB (e j )
FT[rxx (m)] FT[wB (m)]
1 N
sin(N / 2) sin( / 2)
2
WB (ej )称为三角谱窗函数。 式表明，周期图的统计平均值 等于它的真值卷积三角谱窗函数，因此周期图是有偏估计，
BT法：1958年，R.Blackmant和J.Tukey提出， 经典谱估计 先估计自相关函数，再计算功率谱。
周期图法：1898年，Schuster利用傅里叶级数
去拟合待分析的信号，提出周期图的术语，但
直到FFT出现，周期图法才受到人们的重视。
这种方法直接对观测数据进行FFT，取模平方，
除以N得到功率谱。
第四章 功率谱估计
4.1 引言
对信号和系统进行分析研究、处理有两类方法：一类是 时域进行，维纳-卡尔曼滤波和自适应滤波都属于时域处理 方法；另一类是频域方法。
对于确定性信号，傅里叶变换是在频率分析研究的理论 基础，但对于随机信号，其傅里叶变换不存在，因此研究它 的功率谱。
在实际应用中，通常只能采集或观测到平稳随机过程的 一个抽样序列的一段(有限个)数据，如果根据这有限个已知 数据来估计随机过程的功率谱的问题，简称谱估计(谱分析) 问题。
当N 时，var[rˆxx (m)] 0 ，并且
var[rˆxx (m)] var[rˆxx (m)]
rˆxx (m) 虽然是有偏估计，但是渐近一致估计，估计量的方差 小于 rˆxx (m)的方差。实际应用中多用这种有偏自相关估计。 也用符号 rˆxx (m) 表示。
12
BT法功率谱估计采用有偏自相关函数估计法，
X N (e j1 ) 2
X N (e j2 ) 2 ]
1 2
N n k p q
E[x(k)x(n)x( p)x(q)]
RN (k)RN (n)RN ( p)RN (q)e j1(nk e) j2 ( pq) 20
利用正态白噪声、多元正态随机变量多阶矩公式，有 E[x(k)x(n)x( p)x(q)] E[x(k)x(n)]E[x( p)x(q)]
功率谱的计算公式如下：
PˆBT (e j ) L1 Sxx (m)e jm
m0
PˆBT (k) FFT[Sxx (m)]
k 0,1, , L 1
rˆxx (m)w(m) Sxx (m) 0
rˆxx (m L)w(m L)
0 m M 1 M m L M 1 L M m L 1
N m 1 n0
x(n)x(n m)
1 N m N 1
估计性能分析：
估计量的偏差：
Nm rˆxx (m) N rˆxx (m)
因为 rˆxx (m) 是无偏估计，两边取均值，得
Nm
Nm
E[rˆxx (m)] N E[rˆxx (m)] N rxx (m) rxx (m)
有偏估计 10
图b 最大熵谱估计法 (自回归PSD法)
图c Pisarenko谐波分解法
6
4.2 经典谱估计
4.2.1BT法
BT法是先估计自相关函数，再进行傅里叶变换得到功率谱。 下面介绍自相关函数的估计。
利用观测到的实随机序列 x(n)(n 0,1, , N 1)，估计自相关 函数的两种方法是：无偏自相关函数估计和有偏自相关函数 估计。 1. 无偏自相关函数的估计
Pxx
(e
j
)
2 w
H (e j ) 2
w(n)
H (z)
x(n)
如果由观测数据能够估计出信号模型的参数，信号的功率 谱可以计算出来。谱估计问题变成了由观测数据估计信号 模型参数的问题。
模型种类很多，如AR模型、MA模型等。
合适地选择模型，功率谱估计质量比经典谱估计的估计质 量有很大提高。
5
图a BT法
1
按照Weiner-Khintchine定理，信号的功率谱和其自相
关函数服从一对傅里叶变换关系，称为功率谱的定义。
Pxx (e j )
rxx (m)e jm
m
rxx (m)
1
2
Pxx (e j )e jmd
(4.1.1) (4.1.2)
rxx (m) E[x*(n)x(n m)]
(4.1.3)
N m 1 k 1 m N
[rx2x (k) rxx (k m)rxx (k m)][N m k ]
一般观测数据量N很大， N m k
mk
N m k N (1
)N
N
var[rˆxx (m)]
(N
N m )2
N m 1 k 1 m N
[rx2x (k) rxx (k m)rxx (k m)]
Nm rˆxx (m) N rˆxx (m) 知估计量 rˆxx (m) 的方差为
11
var[rˆxx
(m)]
(
N
N
m
)2
var[rˆxx
(m)]
2
N m
var[rˆxx (m)]
N
N
N m 1
(N m )2 k 1 m N
[rx2x (k) rxx (k m)rxx (k m)]
19
式中
E[
x(k
)
x(n)]
rxx
(n
k
)
ห้องสมุดไป่ตู้
2
x
(n
k
)
(n
k
)
1 0
nk nk
var[IN ()]
1 N
n
R2 N
(n)
2 x
2 x
上式说明周期图是无偏估计，这是由于假设信号是实白噪声。
求均方值时，先求两个频率1 和 2处的均方值，最后令 1 2
E[I
N
(1
)I
N
(2
)]
E[
1 N2
1
N 1
N 1
x(k )e jk
x* (n)e jn
N k0
n0
1 N 1 N 1 x(k )x* (n)e j(k n)
N k0 n0 令 m k n，即 k m n，则
Pˆxx (e j ) N1
m( N 1)
1 N 1m N n0
x*
(n)
x(m
n)
e
jm
上式中的方括号部分正是有偏自相关函数的计算公式，则
E[x(n)x( p)]E[x(k)x(q)] E[x(n)x(q)]E[x(k)x( p)]
E[
x(k
)
x(n)
x(
p)
x(q)]
4 x
0
k n, p q; p n, q k;q n, p k 其它
将上式代入周期图的均方值公式中，得
E[I ( )I ( )]= N [N e e ] N 1 N 2
偏移量为
m B rxx (m) E[rˆxx (m)] N rxx (m)
由
Nm E[rˆxx (m)] N rxx (m)
可见，只有当 m 0 时，rˆxx (m) 才是
无偏的，其它 m 都是有偏的，但当N 时，B 0 ，因此
rˆxx (m) 是渐近无偏的。
估计量的方差：
由两种估计的关系式
m
lim
1
N
x* (n)x(n m)e jm
2N 1 m N
n N
lim
1
N
[ x* (n)e jn ][ x(n m)e j(nm) ]
N 2N 1 nN
m
令 l n m ，则
Pxx
(e
j
)
lim
N
1
N
[
2N 1 nN
x(n)e jn ]*[
l
x(l)e jl ]
lim 1
便得到周期图法的定义：
Pˆxx (e j )
1 N
N 1 n0
2
x(n)e jn
上式的绝对值符号内的部分可以用FFT计算，这样周期图法
的计算框图如下图所示。
15
下面证明周期图和BT法的等价关系，再分析周期图的估计质量。 1.周期图和BT法的等价关系
Pˆxx (e j )
1 N
N 1 n0
2
x(n)e jn
E[rˆxx (m)]
N
1
m
N m 1
n0
E[x(n)x(n m)] rxx (m)
B 0 ，这是一种无偏估计。
估计量的方差：
var[rˆxx
(m)]
E{rˆxx
(m)
E[rˆxx
(m)])2}
E[rˆx2x(m)]
r2 xx
(m)
8
可以证明，
var[rˆxx (m)]
(N
1 m )2
这里自相关函数是渐近一致估计，但经过傅里叶变换得到功
率谱的估计却不一定是渐近一致估计，可以证明它是非一致
估计，是一种不好的估计方法。
14
4.2.2周期图法 功率谱的另一种定义式重写如下：
Pxx
(e
j
)
lim
N
E
1 2N 1
N
n N
x(n)e
jn
2
忽略上式中求统计平均的运算，观测数据为 x(n)(0 n N 1)，
对于平稳随机信号，服从各态历经性，集合平均可以用 时间平均代替，可推出功率谱的另一个定义。
将(4.1.3)式中的集合平均用时间平均代替，得
1N
rxx
(m)
lim
N
2N
1
n N
x*(n)x(n m)
(4.1.4)
2
将(4.1.4)式代入(4.1.1)式，得到
Pxx (e j )
rxx (m)e jm
但当 N 时，wB (m) 1，三角谱窗函数趋近于 函数，周期 图的统计平均值趋于它的真值，因此周期图属于渐近无偏估
计。
18
2)周期图的方差
由于周期图方差的精确表示式很复杂，为分析简单，通
常假设
x(n)
是实的零均值的正态白噪声信号，方差是
2，即
x
功率谱是常数
2 x
，其周期图用
IN
() 表示，N表示观测数据
w(m) w(n) 0
(M 其它
1)
m
(M
1)，M
N
要求加窗后的功率谱仍是非负的，这样窗函数的选择必须满足
一个原则，即它的傅里叶变换必须是非负的，如巴特利特窗 13
为了采用FFT计算PˆBT (e j )
M 1
rˆxx (m)w(m)e jm ，设FFT的变换
m ( M 1)
域为(0 ~ L 1) ，必须将求和域 ( M 1, M 1) 移到 (0 ~ L 1)，
rˆxx (m)
1 N
N m 1 n0
x(n)x(n m)
对上式进行傅里叶变换，得到BT法的功率谱估计值为
PˆBT (e j )
rˆxx (m)e jm
m
为了减少谱估计的方差，经常用窗函数 w(n)对自相关函数进行
加权，此时谱估计公式为
式中
PˆBT (e j ) M 1 rˆxx (m)w(m)e jm 加权协方差谱估计 m ( M 1)
此式表明，只有当 N m, N 时，估计量的方差才趋
于0，此估计是渐近一致估计。但是当 m N 时，方差将很大，
因此，这种方法在一般情况下不是一种好的估计方法。虽然是
无偏的，但不能算是一致的。
9
2. 有偏自相关函数的估计
有偏自相关函数用 rˆxx (m) 表示，估计器为
rˆxx (m)
1 N
估计器为
rˆxx (m)
N
1 m
N m1 n0
x(n)x(n m)
0 m N 1
rˆxx (m) rˆxx (m)
1 N m 0
7
也可写成一个表达式
rˆxx (m)
N
1
m
N m 1
n0
x(n)x(n m)
估计性能分析：
估计量的偏差：
1 N m N 1
B E[rˆxx (m)] rxx (m)
的长度，则
IN
()
1 N
X (e j ) 2 1 N 1 N 1 N n0 k0
x(k)x(n)e jke jn
var[I
N
()]
E[
I
2 N
()]
E
2
[I
N
()]
下面先求周期图的均值，再求其均方值：
E[IN ()]
1 N
X (e j ) 2
1
N n k
E[x(k )x(n)]RN (k )RN (n)e j(nk )
经典谱估计致命的缺点是频率分辨率低，原因是傅里叶变
换域是无限大，而观测数据是有限长，观测不到的数据被
认为是0。这相当于将信号在时域加了矩形窗，在频域使真
正的功率谱卷积一个sinc函数。
4
现代谱估计 参数模型法：AR、MA、ARMA模型 非参数模型法：Pisarenko、MUSIC、ESPRIT法
参数模型法以信号模型为基础。
Pˆxx (e j ) N 1 rˆxx (m)e jm
m( N 1)
Pˆxx (e j ) PˆBT (e j )
等价
16
2.周期图法谱估计质量分析 1)周期图的偏移
已知自相关函数的估计值 rˆxx (m)，(N 1) m N 1，求功率
谱的统计平均值，
Pˆxx (e j ) N 1 rˆxx (m)e jm
N
2
x(n)e jn
N 2N 1 nN
(4.1.5)
Pxx (ej )是随机变量，必须对 Pxx (ej ) 取统计平均值，得到
Pxx
(e
j
)
lim
N
E
1 2N 1
N
n N
x(n)e jn
2
上式被认为是功率谱的另一定义。
(4.1.6)
3
经典谱估计，也称为线性谱估计 谱估计方法
现代谱估计，也称为非线性谱估计
m( N 1)
E[Pˆxx (e j )]
N 1
E[rˆxx (m)]e jm
m( N 1)
N 1
m( N 1)
N
m N
rxx (m)e jm
wB (m)rxx (m)e jm
m
17
式中
N m
wB
(m)
N
m N
0
其它
式 中乘积的傅里叶变换，在频域服从卷积关系，得到
E[Pˆxx
(e
j
)]