第四章 功 率 谱 估 计第四章 功 率 谱 估 计
4.1 引言
4.2 经典谱估计
4.3 现代谱估计中的参数建模
4.4 AR模型谱估计的性质
4.5 AR谱估计的方法
4.6 最大熵谱估计与最大似然谱估计
4.7 特征分解法谱估计第四章 功 率 谱 估 计
4.1 引 言我们知道,对信号和系统进行分析研究,处理有两类方法：一类是在时域进行,前面我们学习的维纳 —卡尔曼滤波和自适应滤波都属于这种方法；本章则是在频率域进行研究的另一类方法 。
这两类方法都是信号处理的重要方法 。 对确定性信号傅里叶变换是在频率域分析研究的理论基础,但对于随机信号,其傅里叶变换并不存在,因此转向研究它的功率谱 。 按照 Weiner-
Khintchine定理,信号的功率谱和其自相关函数服从一对傅里叶变换关系,公式如下
mj
m
xx
j
xx mrP
e)()e(
(4.1.1)
第四章 功 率 谱 估 计
dPmr njjxxxx e)e(π2 1)(
)]()([)( * mnxnxEmr xx
(4.1.2)
(4.1.3)
（ 4.1.1） 式被称做功率谱的定义,对于平稳随机信号,服从各态历经定理,集合平均可以用时间平均代替,由 （ 4.1.1）
式还可以推出功率谱的另一个定义,推导如下：
将 （ 4.1.3） 式中的集合平均用时间平均代替,得到
)()(12 1li m)( * mnxnxNmr
N
NnN
xx
(4.1.4)
第四章 功 率 谱 估 计将（ 4.1.4）式代入 (4.1.1)式,得到
m
mnjnj
N
Nn
N
nj
m
N
m
j
xx
mnxnx
N
mnxnx
N
P
)(*
*
e)(e)(
12
1
lim
e)()(
12
1
lim)e(
令 l=n+m,则
2
*
j
e)(
12
1
lim
e)()(
12
1
lim)e(
j ωω
N
Nn
N
lj
l
N
Nn
nj
N
xx
nx
N
lxenx
N
P
(4.1.5)
第四章 功 率 谱 估 计上式中 x(n)是观测数据,Pxx(ejω)是随机变量,必须对 Pxx(ejω)
取统计平均值,得到
2
j ω )(
12
1lim)e( N
Nn
nj
Nxx
enx
N
EP?(4.1.6)
上式被认为是功率谱的另一定义 。
（ 4.1.1） 式表明功率谱是无限多个自相关函数的函数,但观测数据只有有限个,只能得到有限个自相关函数 。 按照
（ 4.1.6） 式求功率谱,也需要无限个观测数据 。 因此根据有限个样本数据,分析计算随机序列的真正功率谱,是求功率谱的中心问题,毫无疑问,这是一个功率谱的估计问题 。 在第一章已介绍了统计估计的一般估计准则,主要有偏移,估计量方差和估计量的均方误差 （ 有效性 ）,这里不再重复,下面直接用它们分析估计质量 。
第四章 功 率 谱 估 计现代谱估计以信号模型为基础,图 4.1.1表示的是 x(n)的信号模型,输入白噪声 w(n)均值为 0,方差为 σ 2w,x(n)的功率谱由下式计算：
2j2j |)e(|)e( HP wxx?
（ 4.1.7）
如果由观测数据能够估计出信号模型的参数,信号的功率谱可以按照 （ 4.1.7） 式计算出来,这样,估计功率谱的问题变成了由观测数据估计信号模型参数的问题 。 模型有很多种类,例如
AR模型,MA 模型等等,针对不同的情况,合适地选择模型,
功率谱估计质量比较经典谱估计的估计质量有很大的提高 。 遗憾的是,尚无任何理论能指导我们选择一个合适的模型,我们只能根据功率谱的一些先验知识,或者说一些重要的谱特性,
选择模型 。
第四章 功 率 谱 估 计
H ( z )
w ( n ) x ( n )
图 4.1.1 平稳随机序列的信号模型第四章 功 率 谱 估 计
4.2 经 典 谱 估 计
4.2.1 BT法
BT法是先估计自相关函数,然后按照 (4.1.1)式进行傅里叶变换得到功率谱 。 设对随机信号 x(n),只观测到一段样本数据,
n=0,1,2,…,N-1。 关于如何根据这一段样本数据估计自相关函数,第一章已经作了详细介绍,结果是共有两种估计方法,即有偏自相关函数估计和无偏自相关函数估计 。 有偏自相关函数估计的误差相对较小,这种估计是一种渐近一致估计,将该估计公式重写如下：
1||
0
* )()(1)(?
mN
n
xx mnxnxNmr
（ 4.2.1）
第四章 功 率 谱 估 计对上式进行傅里叶变换，得到 BT法的功率估计值为
m
j ωω
xx
j mrP -
BT e)(?)e(
(4.2.2)
为了减少谱估计的方差,经常用窗函数 w(m)对自相关函数进行加权,此时谱估计公式为
j ωω
xx
Mm
j mrP -
)1(
BT e)(?)e(
（ 4.2.3）
式中
0
)()( mwmw
-(M-1)≤ m≤ (M-1)
其它,M≤ N
（ 4.2.4）
第四章 功 率 谱 估 计有时称 (4.2.3)式为加权协方差谱估计 。 它要求加窗后的功率谱仍是非负的,这样窗函数 w(m)的选择必须满足一个原则,即它的傅里叶变换必须是非负的,例如巴特利特窗就满足这一条件 。
为了采用 FFT计算 （ 4.2.3） 式,设 FFT的变换域为 (0~L-1)，
必须将求和域 (-M+1,M-1)移到 （ 0~L-1）,功率谱的计算公式如下：
)()(?
e)()e(?
BT
1
0
j-j ω
BT
mSFFTkP
mSP
xx
L
m
ωm
xx
k=0,1,2,…,L-1
第四章 功 率 谱 估 计
(4.2.7)
)()(?
0
)()(?
)(
LmwLmr
mwmr
mS
xx
xx
xx
0≤ m≤ M-1
M≤ m≤ L-M-1
L-M≤ m≤ L-1
按照 （ 4.2.1） 式估计自相关函数,我们已经证明这是渐近一致估计,但经过傅里叶变换得到功率谱的估计,功率谱估计却不一定仍是渐近一致估计,可以证明它是非一致估计,是一种不好的估计方法 。 下面我们将证明,BT法中用有偏自相关函数进行估计时,它和用周期图法估计功率谱是等价的,因此 BT
法估计质量和周期图法的估计质量是一样的 。
第四章 功 率 谱 估 计
4.2.2 周期图法将功率谱的另一定义（ 4.1.6）式重写如下,
2
jj e)(
12
1
lim)e(
N
Nn
n
Nxx
nx
N
EP
如果忽略上式中求统计平均的运算,观测数据为,x(n) 0≤n≤N-1，
便得到周期图法的定义,
21
0
j-j e)(1)e(
N
n
n
xx nxNP
(4.2.8)
第四章 功 率 谱 估 计图 4.2.1 用周期图法计算功率谱框图观测数据 x ( n ) FFT
取模的平方
1/ N
P xx (e j )
^
第四章 功 率 谱 估 计
1,周期图与 BT法的等价关系周期图法的功率谱估计公式用 （ 4.2.8） 式表示,下面由该公式出发推导它们的等价关系 。
1
0
1
0
)(*
1
0
*
1
0
2
1
0
e)()(
1
e)(e)(
1
e)(
1
)e(?
N
k
N
n
nkj
N
n
nj
N
k
nj
N
n
njj
xx
nxkx
N
nxkx
N
nx
N
P
令 m=k-n,即 k=m+n，则
)(
1
0
*j e)()(1)e(? nkj
mN
n
xx nmxnxNP
第四章 功 率 谱 估 计上式中的方括号部分正是有偏自相关函数的计算公式,因此得到
mj
xx
N
Nm
xx mrP
)e(?)e(?
1
)1(
j
)()(1)(? *
||1
0
nmxnxNmr
mN
n
xx
)e(?)e(? BT jjxx PP?
因此证明了利用有偏自相关函数的 BT法和周期图法的等价关系。
第四章 功 率 谱 估 计
2,周期图法谱估计质量分析
1)
已知自相关函数的估计值,m=-(N-1),-N,-N+1,…,0,
1,2,…,N-1,按照 （ 4.2.2） 式求功率谱的统计平均值,得到
)(? mrxx
nxxN
Nm
xx mrEPE
j-
1
)1(
j e)(?)]e(?[?
有偏自相关函数统计平均值已由第一章 (1.3.30)式确定,将该式代入上式,得到
n
xx
N
Nm
xx mrN
mNPE j-1
)1(
j )e(||)]e(?[
第四章 功 率 谱 估 计
m
n
xxBxx mrmwPE
j-j )e()()]e(?[ （ 4.2.9）
式中
其它0
||
||
)(
Nm
N
mN
mw B （ 4.2.10）
(4.2.9)式中,两序列乘积的傅里叶变换,在频域服从卷积关系，
得到
1)]e(?[ jxxPE
ΩWPE jBjxx d)e(π2 1)]e(?[ )(（ 4.2.11）
式中
)]([)e( mrFTP xxjxx
第四章 功 率 谱 估 计
2
)2/s i n (
)2/s i n (1)]([)e(
N
N
mwFTW BjB
(4.2.12)
WB(ejω)称为三角谱窗函数 。 （ 4.2.11） 式表明,周期图的统计平均值等于它的真值卷积三角谱窗函数,因此周期图是有偏估计,但当 N→∞ 时,wB(m)→ 1,三角谱窗函数趋近于 δ 函数,周期图的统计平均值趋于它的真值,因此周期图属于渐近无偏估计 。
第四章 功 率 谱 估 计
2)
由于周期图的方差的精确表示式很繁冗,为分析简单起见,
通常假设 x(n)是实的零均值的正态白噪声信号,方差是 σ x2，
即功率谱是常数 σ x2,其周期图用 IN(ω )表示,N表示观测数据的长度 。 按照周期图的定义,周期图表示为
)]([)]([)](v a r [
ee)()(
1
|)e(|
1
)(
22
j-
1
0
1
0
j2j
NNN
n
N
n
N
k
k
N
IEIEI
nxkx
N
X
N
I
下面先求周期图的均值，再求其均方值：
)(j-e)()()]()([1)]([ n - k
N
k
N
n
N kRnRkxnxENIE
第四章 功 率 谱 估 计式中
kn
kn
kn
knknrkxnxE xxx
0
1
)(
)()()]()([
2
222 )(1)]([
x
n
xNN nRNIE
（ 4.2.13）
上式说明周期图是无偏估计,但前面已推导出周期图是有偏估计
（ 一般情况 ）,这里由于对信号作了实白噪声的假设,才有无偏估计的结果 。 在求均方值时,先求两个频率 ω 1和 ω2处的均方值,
最后令 ω =ω 1=ω 2。
第四章 功 率 谱 估 计
)(j-)(j-
2
2
j
2
j
221
21
21
ee)]()()()([
)()()()(
1
)e()e(
1
)]()([
qpkn
NNN
q
N
pkn
NNNN
qxpxkxnxE
qRpRkRnR
N
XX
N
EIIE
利用正态白噪声、多元正态随机变量的多阶矩公式，有
)]()([)]()([
)]()([)]()([)]()([)]()([)]()()()([
pxkxEqxnxE
qxkxEqxpxEqxpxEkxnxEqxpxkxnxE
其它0
,;,,,
)]()()([
4 kpnqkqnpqpnk
qxkxnxE x
第四章 功 率 谱 估 计将上式代入周期图的均方值公式中,得到
2
21
21
2
21
214
1
0
1
0
1
0
))((
1
0
))((2
2
4
21
2/)s i n (
2/)s i n (
2/)s i n (
2/)s i n (
1
ee)]()([ 2121
N
N
N
N
N
N
IIE
x
N
n
N
n
N
k
n - kj
N
k
n - kjx
NN
(4.2.14) 将 ω =ω 1=ω 2代入上式，得到
2
4
42
)s i n (
)s i n (
1)](v a r [
)s i n (
)s i n (
2)]([
N
N
I
N
N
IE
xN
xN
(4.2.15)
第四章 功 率 谱 估 计显然,当 N趋于无限大时,周期图的方差并不趋于 0,而趋于功率谱真值的平方,即
4)](v a r [ x
NNI
(4.2.16)
这里无论怎样选择 N,周期图的方差总是和 σ 4x同一个数量级 。 我们知道,信号的功率谱真值是 σ 2x,说明周期图的方差很大,周期图的均方误差也是非常大 。 用这种方法估计的功率谱在 σ 2x附近起伏很大,故周期图是非一致估计,是一种很差的功率谱估计方法 。
第四章 功 率 谱 估 计为了进一步说明数据长度 N对功率谱估计的影响，下面求两个频率处的协方差函数。
)]()([)]()([)](),(c o v [ 212121 NNNNNN IIEIIIII
将（ 4.2.13）式和 (4.2.14)式代入上式,得到
2
21
21
2
21
214
21 2/)s in (
2/)(s in
2/)s in (
2/)(s in)]()(c o v [
N
N
N
NII
xNN
令,ω 1=2π k/N,ω 2=2π l/N,式中 k,l均是整数，得到
224
21 ]/)(πs in [
)](πs in [
]/)(πs in [
)](πs in [)]()(c o v [
NlkN
lk
NlkN
lkII
xNN
(4.2.17)
第四章 功 率 谱 估 计图 4.2.2 白噪声的周期图
0 1 2 3
5
4
3
2
1
0
0 1 2 3
5
4
3
2
1
0
P
xx
(e
j
)
P
xx
(e
j
)
5
4
3
2
1
0
0 1 2 3
0 1 2 3
5
4
3
2
1
0
P
xx
(e
j
)
P
xx
(e
j
)
/?
/?
( a ) ( b )
/ /?
( c ) ( d )
N ＝ 32 N ＝ 64
N ＝ 256
N ＝ 128
第四章 功 率 谱 估 计
4.2.3 经典谱估计方法改进
1,平均周期图法平均周期图法是基于这样的思想：对一个随机变量进行观测,得到 L组独立记录数据,用每一组数据求其均值,然后将 L
个均值加起来求平均 。 这样得到的均值,其方差将是用一组数据得到的均值的方差的 1/L。
假设随机信号 x(n)的观测数据区间为,0≤n≤M-1,共进行了
L次独立观测,得到 L组记录数据,每一组记录数据用 xi(n),i=1,
2,3,…,L表示,第 i组的周期图用下式表示：
21
0
)(1)(?
M
n
nj
ii enxMI
(4.2.18)
第四章 功 率 谱 估 计将得到的 L个周期图进行平均，作为信号 x(n)的功率谱估计，
公式如下：
L
i
ixx ILP
1
j )(1)e(
(4.2.19)
为了分析偏移，对上式求统计平均，得到
2
j
)-j(j
π
π-
j
)2/s i n (
)2/s i n (1
)]([)e(
)e()e(
π2
1
)]e(?[
M
M
mwFTW
dPWPE
BB
xxBxx
(4.2.20)
(4.2.21)
周期图的统计平均值已经求出,如（ 4.2.11）,(4.2.12)式所示，
重写如下：
第四章 功 率 谱 估 计上式表明,平均周期图仍然是有偏估计,偏移和每一段的数据个数 M有关；由于 M<N,平均周期图的偏移比周期图的偏移大,表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽 。 由于三角谱窗主瓣的宽度变宽,分辨率更加降低,因此也可以说,
偏移的大小反映分辨率的低与高 。
按照 （ 4.2.19） 式求方差,由于是 L次独立观测,L个周期图相互独立,因此平均周期图的方差为
)](v a r [1)]e(?v a r [ ijxx ILP? （ 4.2.22）
即平均周期图的估计方差是周期图的方差的 1/L。 显然,是以分辨率的降低换取了估计方差的减少,当然,估计的均方误差也减少 。
第四章 功 率 谱 估 计图 4.2.3 平均周期图法
P
xx
(e
j
)
/?
0
3
2
1
0 1 2 3
0
3
2
1
0 1 2 3
0
3
2
1
0 1 2 3
P
xx
(e
j
)
/?
P
xx
(e
j
)
/?
( a ) ( b )
( c )
N ＝ 256,L ＝ 2 N ＝ 256,L ＝ 4
N ＝ 256,L ＝ 8
第四章 功 率 谱 估 计
2.
这种方法是用一适当的功率谱窗函数 W(ejω)与周期图进行卷积,来达到使周期图平滑的目的的 。
d)e()(π2 1)e(? )(π jNlwxx WIP （ 4.2.23）
式中
m
xx
N
Nm
N mrI
j-
1
)1(
e)(?)(?
是有偏自相关函数)(? mr
xx
de)e(π2 1)( jjπ
π
nWnw?
-(M-1)≤ n≤ M-1
（ 4.2.24）
第四章 功 率 谱 估 计那么
m
xx
M
Mm
j
xx emwmrP
j
1
)1(
)()(?)e(
(4.2.25)
将 (4.2.25)式和 (4.2.3)式进行对比,它们是一样的,说明周期图的窗函数法就是前面提到的 BT法的加权协方差谱估计 。
在 （ 4.2.23） 式中,周期图和谱窗函数卷积得到功率谱,等效于在频域对周期图进行修正,使周期图通过一个线性非频变系统,滤除掉周期图中的快变成分,谱窗函数需具有低通特性 。
第四章 功 率 谱 估 计对（ 4.2.25）式求统计平均,得到
m
xx
M
Mm
xx mwmrEPE
j-
1
)1(
j e)()](?[)]e(?[?
将 (1.3.30)式和 (1.3.32)式代入上式，得到
m
Bxx
M
Mm
xx mwmwmrePE
j-
1
)1(
j e)()()()](?[?
(4.2.26)
式中
Nm
N
mmw
B ||
||1)(
第四章 功 率 谱 估 计上式表明,周期图的窗函数法仍然是有偏估计,其偏移和 wB(m)、
w(m)两个窗函数有关,如果 w(m)窗的宽度比较窄,M比 N小得多,这样 |m|<<N,则 wB(m)~1,（ 4.2.26） 式可以近似写成下式：
d)e()e(
π2
1
e)()()]e(?[
)-j(j
π
π
j-
1
)1(
j
WP
mwmrPE
xx
m
xx
M
Mm
xx
(4.2.27)
第四章 功 率 谱 估 计
3,修正的周期图求平均法这种方法和平均周期图法一样,首先把数据长度为 N的信号
x(n)分成 L段,每一段数据长度为 M,N=LM。 然后把窗函数 w(n)
加到每一个数据段上,求出每一段的周期图,形成修正的周期图,再对每一个修正的周期图进行平均 。 第 i段的修正周期图为
21
0
j-e)()(1)(?
M
n
ωn
ii nwnxUI?
i=1,2,3,…,M-1
(4.2.28)
式中
1
0
2 )(1
M
n
nw
M
U
(4.2.29)
第四章 功 率 谱 估 计同样,将每一段的修正的周期图之间近似看成互不相关,最后功率谱估计为
L
i
ixx ILP
1
j )(1)e( (4.2.30)
对上式求统计平均,得到
d)e()e(π2 1)]e(?[ )-j(j
π
j WPPE
xxxx
(4.2.31)
式中
21
0
jj e)(1)e(?
M
n
nw
MU
W (4.2.32)
第四章 功 率 谱 估 计
U称为归一化因子 。 (4.2.31)式的证明可参考文献 ［ 2］ 。 为使功率谱估计渐近无偏,（ 4.2.28） 式和 (4.2.32)式中的归一化因子是不可缺少的 。 这种在计算周期图之前,先对各数据段加窗函数的方法,使平均周期图方法的估计方差减少,当然分辨率同样减少 。 但这种方法对窗函数没有限制；此外,分段时,相邻的两段可以有重叠,进一步使方差减少,可以重叠
50%［ 6］ 。 总之,传统的功率谱估计方法无论采取哪一种改进方法,总是以减少分辨率为代价,换取估计方差的减少,提高分辨率的问题无法根本解决 。
第四章 功 率 谱 估 计对于由白噪声和正弦信号或者窄带信号组成的信号,在计算周期图之前,给数据加窗是有利的,如果不加数据窗,相近的低电平信号可能被高电平信号的旁瓣淹没掉 。 在图 4.2.4（ a）
中,ω/π =0.12处的正弦信号的旁瓣几乎掩盖了在 ω/π =0的信号,经过给数据加哈明窗处理后,功率谱如图 4.2.4(b)所示 。
由于加了哈明窗,大大压低了旁瓣,使低电平信号清晰可见,
但由于主瓣加宽,功率谱波峰变宽了,却降低了信号的分辨率 。
一般来说,两个等幅的正弦信号的频率相隔很近,可以不加数据窗,频率间隔应该大于 2π /N才能分辨 。
第四章 功 率 谱 估 计图 4.2.4
(a) 没有数据窗； (b) 加哈明数据窗
4 0,0 0
2 0,0 0
0,0 0
－ 2 0,0 0
－ 0,5 0 － 0,3 0 － 0,1 0 0,1 0 0,3 0 0,5 0
P
xx
(e
j
)
/
d
B
( a )
4 0,0 0
2 0,0 0
0,0 0
－ 2 0,0 0
－ 0,5 0 － 0,3 0 － 0,1 0 0,1 0 0,3 0 0,5 0
P
xx
(e
j
)
/
d
B
( b )
/ /?
第四章 功 率 谱 估 计
4.3 现代谱估计中的参数建模在第一章中我们介绍了用参数模型来描述随机信号的方法,
如果能确定信号 x(n)的信号模型,根据信号观测数据求出模型参数,系统函数用 H(z)表示,模型输入白噪声方差为 σ w2,信号的功率谱用下式求出：
2j2j |)e(|)e( HP wxx?
按照这种求功率谱的思路,功率谱估计可分成三个步骤：
（ 1） 选择合适的信号模型；
（ 2） 根据 x(n)有限的观测数据,或者它的有限个自相关函数估计值,估计模型的参数 ;
（ 3） 计算模型输出功率谱。
第四章 功 率 谱 估 计
4.3.1 模型选择选择模型的主要考虑是模型能够表示谱峰,谱谷和滚降的能力 。 对于具有尖峰的谱,应该选用具有极点的模型,如 AR和
ARMA模型；对于具有平坦的谱峰和深谷的信号,可以选用 MA模型；既有极点也有零点的谱应选用 ARMA模型,相对地说,ARMA
模型适用范围较宽 。 对于滚降太快的谱,没有一种模型可以准确地表示功率谱,可以选用高阶的 AR模型近似表示 。 如果选择不合适,例如,选用 MA模型去估计具有尖峰的功率谱,估计效果会很差 。 图 4.3.1表示的是用 MA模型估计二阶 AR信号功率谱的例子,图中 (a),（ b）,（ c） 中的 AR信号谱峰较平坦,用二阶 MA模型功率谱拟和真实谱时,差别较大,随着阶数的提高,
估计的谱愈来愈近似于真实的谱 。 但是对于 （ d）,（ e）,
（ f） 中 AR信号谱峰很窄,再用 MA信号模型拟和时,直到 MA模型阶数提高到 10阶,其效果仍很差 。
第四章 功 率 谱 估 计实际信号中一般都含有和信号不相关的噪声,对这一类带有噪声的信号,如果信号是 AR模型,由于有噪声需要用 ARMA模型 ; 如果用 AR模型,则需要阶数更高 ［ 5］ 。
选择模型的另一个考虑是尽量减少模型参数 。 当然这和选择模型是否合适有关系,前面介绍过三种信号模型均有普遍应用价值,但模型选择合适,功率谱估计和实际的谱拟合得好,
如果不合适,只有提高阶数才能得到较近似的谱,这样需要估计的参数增多,同样也会降低谱估计的质量 。 因此应该在选择模型合适的基础上,尽量减少模型的参数 。
第四章 功 率 谱 估 计
4.3.2 模型参数和自相关函数之间的关系根据观测数据如何来确定模型的参数对各种功率谱估计方法不尽相同 。 这里先介绍模型参数和信号自相关函数之间的关系,这些关系在功率谱估计中起着很重要的作用 。
假设模型的差分方程和系统函数分别用下式表示：
q
i
i
P
i
i inwbinxanx
01
)()()(
)(
)(
1
)(
1
0
zA
zB
za
zb
zH
P
i
i
i
q
i
i
i
(4.3.1)
(4.3.2)
第四章 功 率 谱 估 计图 4.3.1 对 AR(2)信号的模型选择
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
15
20
25
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
/?
/?
/?
/?
( a ) ( b )
( c )
( d )
dB
dB
dB
dB
真实 A R( 2 ) 的 P S D
MA (2 ) 模型真实 A R( 2 ) 的 P S D
MA (5 ) 模型真实 A R( 2 ) 的 P S D
MA (2 ) 模型真实 A R( 2 ) 的 P S D
MA (1 0 ) 模型第四章 功 率 谱 估 计
10
5
0
－ 5
－ 10
15
20
25
10
5
0
－ 5
－ 10
15
20
25
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,50,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
/ /?
( e ) ( f )
dBdB
真实 A R(2 ) 的 PS D
MA (1 0 ) 模型真实 A R(2 ) 的 PS D
MA (5 ) 模型图 4.3.1 对 AR(2)信号的模型选择第四章 功 率 谱 估 计
1,ARMA模型的系数和信号自相关函数之间的关系在第一章我们已推导出系统输出功率谱与输入功率谱之间的关系,用 （ 1.4.7） 式表示,重写如下：
*
* 1)()()(
zHzHzPzP wwxx
式中,Pxx(z)和 Pww(z)分别表示模型输出信号 x(n） 和输入信号 w(n)
功率谱的 Z变换形式 ; H（ z） 表示模型的系统函数 。 将上式两边同乘以 A(z),得到
)(
1
)(
1
)()()(
*
*2
*
*
zB
z
H
zB
z
HzPzAzP
w
wwxx
(4.3.3)
第四章 功 率 谱 估 计将上式进行 Z反变换，公式左边的反变换为
P
l
xxAAxxxx lmrlhmhmrzAzPZ
0
1 )()()()()]()([T
公式右边的反变换为
)()()(1T
0
2
*
*21 lmqlhzB
zHZ
q
l
Bww
式中
)(1T)( ***1 mhzHZmq
第四章 功 率 谱 估 计代入上式,得到
q
l
Bww lmhlhzBzHZ
0
*2
*
*21 )()()(1T
这样得到
q
l
Bwxx
p
l
A lmhlhlmrlh
0
*2
0
)()()()(?
由于模型 H(z)是因果的,h(n)=0,n＜ 0,可以得到
p
l
xxA
p
l
mq
l
BwxxA
xx
lmrlh
mlhlhlmrlh
mr
1
1 0
*2
)()(
)()()()(
)(
m=0,1,2,3,…,q
m≥ q+1
第四章 功 率 谱 估 计此式即是 ARMA模型输出自相关函数与模型参数之间的关系式 。
如果能由信号的观测数据估计出信号的自相关函数,可以按照
（ 4.3.4） 式求出 ARMA模型的参数,再由 （ 1.4.8） 式求出信号的功率谱 。 但我们看到,在公式中由于 项的存在,模型输出自相关函数和模型参数之间的关系是非线性的,
从而增加了估计功率谱的困难 。 但是当 m>q时,(4.3.4)式却是一个线性方程,用矩阵方程表示如下：
)()(
0
* mlhlh
mq
l
B
第四章 功 率 谱 估 计
)(
)2(
)1(
)(
)2(
)1(
)()2()1(
)2()()1(
)1()1()(
pqr
qr
qr
pqh
h
h
qrpqrpqr
pqrqrqr
pqrqrqr
xx
xx
xx
A
A
A
xxxxxx
xxxxxx
xxxxxx
上式共有 p个方程 。 可以用该方程首先计算出 AR部分的 p个系数
hA(i),i=1,2,3,…,p; 然后代入 (4.3.4)式,设法求出 MA部分的系数 。
第四章 功 率 谱 估 计
2.AR模型的系数和信号自相关函数之间的关系
AR模型的系统函数为,H(z)=1/A(z),相当于 ARMA模型中
B(z)=1的情况,这样在公式中
ml
ml
mlmlh
nnh
B
B
0
1
)()(
)()(
因为 h(n)是因果性的,因此 m＞ 0时 h*(-m)=0,将上面公式代入到
（ 4.3.4） 式中,得到
P
l
wxxA
P
l
xxA
xx
lrlh
mrlh
mr
1
2
1
)()(
)1()(
)(
m≥ 1
m=0
(4.3.6)
第四章 功 率 谱 估 计也可以将上式中 m≥ 1的情况写成矩阵形式：
)(
)2(
)1(
)(
)2(
)1(
)0()2()1(
)2()0()1(
)1()1()0(
pr
r
r
ph
h
h
rprpr
prrr
prrr
xx
xx
xx
A
A
A
xxxxxx
xxxxxx
xxxxxx
(4.3.7) （ 4.3.6）式的矩阵形式如下式：
0
0
1
)0()1()(
)2()0()1(
)1()1()0( 2
1
w
p
xxxxxx
xxxxxx
xxxxxx
h
a
rprpr
prrr
prrr?
(4.3.8)
第四章 功 率 谱 估 计或者用模型参数表示：
0
0
1
)0()1()(
)1()0()1(
)()1()0( 2
1
w
p
xxxxxx
xxxxxx
xxxxxx
h
a
rprpr
prrr
prrr?
(4.3.9) 令
)0()1()(
)1()0()1(
)()1()0(
)1(
xxxxxx
xxxxxx
xxxxxx
p
xx
rprpr
prrr
prrr
R
(4.3.10)
第四章 功 率 谱 估 计称为自相关矩阵,它满足 （ H表示共轭转置 ）,
是一个埃尔米特 （ Hermitian） 矩阵,且沿任一对角线的元素相等,是一个托布列斯 （ Toeplitz） 矩阵 。 也是正定矩阵 。
上面推导出的 （ 4.3.6） 式或 (4.3.9)式确定了 AR模型参数
（ 包括模型输入噪声方差 ） 和信号自相关函数之间的关系 。
我们注意到这是一个线性方程,如果能够由信号的观测数据求出其自相关函数,可以按照 （ 4.3.7） 式,通过解一组线性方程得到模型参数,相对 ARMA模型,这是 AR模型的优点 。
)1(?pxxR H
xxxx RR?
)1(?pxxR
第四章 功 率 谱 估 计
3.MA
MA模型的系统函数 H(z)=B(z),相当于 ARMA模型中 A(z)=1,
hA(n)=δ (n)的情况,此时 h(n)=hB(n),由 （ 4.3.4）
MA模型系数和信号自相关函数的关系为
0
)()(
)( 0
*2
mq
l
BBw
xx
mlhlh
mr
m=0,1,…,q
m≥ q+1 (4.3.11)
上式表明,MA模型的参数和信号自相关函数之间也是非线性关系。
第四章 功 率 谱 估 计上面我们分别推导了三种信号模型的参数和信号自相关函数之间的关系 。 这些关系式为我们提供了一种估计功率谱的方法,即首先根据信号观测数据估计信号自相关函数,再按照所选择信号模型,解上面相应的方程,求出模型参数,最后按照下式求出信号的功率谱,
2
1
022j ω2j ω
e1
e
|)e(|)e(
P
i
i
i
q
i
i
i
wwxx
a
b
HP
(4.3.12)
第四章 功 率 谱 估 计
4.4 AR模型谱估计的性质
4.4.1 AR模型的线性预测在第二章中,我们已推导出维纳线性一步预测器系数和信号自相关函数之间的关系式 (称为 Yule-Walker方程 ),重写如下：
0
0
)]([1
)0()1()(
)1()0()1(
)()1()0(
m i n
2
1
neE
a
a
rprpr
prrr
prrr
pp
p
xxxxxx
xxxxxx
xxxxxx
(4.4.1)
第四章 功 率 谱 估 计式中,e(n)表示线性一步预测误差，其公式为
p
i
pi inxanxnxnxne
1
)()()(?)()(
(4.4.2)
E［ e2(n)］ min表示 e(n)的均方差最小值； api(i=1,2,3,…,p)表示预测器的系数,它和线性预测器单位脉冲响应 h(n)之间差一符号,
即
pnanh pn,,3,2,1)(
对（ 4.4.2）式进行 Z变换,得到
p
i
i
pi zazXzXzE
1
)()()(
(4.4.3)
第四章 功 率 谱 估 计令 He(z)=E(z)/X(z)，由上式，得到
p
i
i
pie zazH
1
1)(
(4.4.4)
称 He(z)为线性一步预测误差滤波器,其作用是将信号 x(n)转换成预测误差 e(n),如图 4.4.1所示 。 一般认为 e(n)具有白噪声的性质,因此 He(z)也称为白化滤波器 。
第四章 功 率 谱 估 计图 4.4.1 预测误差滤波器
H
e
( z )
x ( n ) e ( n )
x ( n )
e ( n )
z － 1 z － 1 z － 1
a
p 1
a
p 2
a
p,p － 1
a
p,p
第四章 功 率 谱 估 计我们知道 AR模型的系统函数为
P
i
i
i za
zH
1
1
1)( (4.4.5)
对比 （ 4.4.4）,(4.4.5)两式,当 api=ai(i=1,2,3,…,p)时,
He(z)和 H(z)互为逆滤波器,He(z)=1/H(z),基于以上分析,也可以将 AR模型定义为
)()(?)( nwnxnx
(4.4.6)
式中,是基于信号前 p个样本的最佳一步线性预测，公式为)(? nx
p
i
i
pi zanx
1
)(?
第四章 功 率 谱 估 计
w(n)是模型输入白噪声 。 将 AR模型参数和信号自相关函数之间的关系式 （ 4.3.9 ） 和 (4.4.1) 式进行对比,得到,api=ai,
w(n)=e(n),E［ e2(n)］ min=σ2w,信号自相关函数和它的 AR模型参数之间的关系服从 Yule-Walker方程 。 注意,只有当 AR模型的阶数与线性预测器的阶数相同时,以上结论才是正确的 。 由于
AR模型具有这种特性,因而 AR模型法也称为线性预测 AR模型法 。
第四章 功 率 谱 估 计
AR模型与线性预测之间的关系,可以被用来解卷积,假设信号 s(n)通过一个 AR系统,系统单位脉冲响应为 h(n),响应是
x(n),即
x(n)=s(n) * h(n)
如果 s(n)具有白噪声性质,可以利用 AR模型与预测滤波器之间的关系对上式进行解卷积,得到 s(n)信号 。 方法是：先由
x(n)的观测数据估计 AR模型的参数,得到 AR模型的系统函数
H(z),再让 x(n)通过 H(z)的逆滤波器 H-1(z),便得到信号 s(n)。
第四章 功 率 谱 估 计
4.4.2 预测误差滤波器的最小相位特性我们知道,AR模型必须因果稳定,即极点均在单位圆内,
才能保证信号 x(n)是平稳随机信号,于是 AR模型 H(z)和预测误差滤波器 He(z)互为逆滤波,那么 He(z)应为最小相位系统 。 但是由解 Yule-Walker方程得到 AR模型的参数,其极点不一定在单位圆内 。 下面将证明当最佳 P阶线性预测系数与 AR模型参数相同时,由此得到的极点保证在单位圆内,AR滤波器稳定,
预测误差滤波器 He(z)或者 A(z)是最小相位系统 。 这里自相关矩阵 是正定的 。)1(?p
xxR
第四章 功 率 谱 估 计解 Yule-Walker方程得到的是最佳线性预测滤波器的系数,
此时预测误差滤波器输出功率 Pe达到最小,用 Pe min表示,即
)d(e|)(e| jππ- 2jmi n xxe PAP
(4.4.7)
式中,,ai是最佳预测系数 。 下面先用反证法证明 A(z)的全部零点不在单位圆外部 （ 即全部零点在单位圆上或者单位圆内部 ） 。 设 A(z)的第 i个零点 zi在单位圆外部,即 |zi|
＞ 1,用 1/z*i代替 zi,这时 A(z)的幅度函数不变,按照 （ 4.4.7） 式计算出的预测误差滤波器输出功率 Pe应不变,仍是最小的 。 但是下面将用公式证明预测误差滤波器输出功率 Pe不是最小的,
这一矛盾的结论只能说明 A(z)的零点不可能在单位圆外部 。
p
i
i
i zazA
1
1)(
第四章 功 率 谱 估 计
)(')1()1()1(
)1(1)(
1
111
1
1
1
zAzzzzzz
zzzazA
p
ij
j
jji
p
j
j
p
i
i
i
(4.4.8)
式中
p
ij
j
j zzzA
1
1 )1()('
将 (4.4.8)式代入（ 4.4.7）式,得到
d)e(|)e('|e1 j2j2j-mi n xxie PAzP
(4.4.9)
第四章 功 率 谱 估 计式中
2
*
2
2
*
j2
2
j22j 11||1||e1|||e1| j
i
i
i
i
i
ii ezzzezzzz
因为
1?iz
因此
2
j-
*
2j- e11|e1|
i
i zz
(4.4.10)
第四章 功 率 谱 估 计
4.4.3 AR模型隐含自相关函数延拓特性
AR模型的自相关函数和模型系数之间的关系服从 Yule-
Walker方程,重写如下：
p
l
wxxA
p
l
xxA
xx
lrlh
lmrlh
mr
1
2
1
)()(
)()(
)(
m≥ 1
m=0
上式中,对于 m≥1的情况,公式本身就是一个递推方程,如果已由观测数据计算出 p+1个自相关函数,用,m=0,1,2,,…,p
表示,对于 m＞ p的情况,可以用该公式外推得到,公式如下：
)(? mrxx
第四章 功 率 谱 估 计
P
l
xxA
xx
xx lmrlh
mr
mr
1
'
'
)(?)(
)(?
)(?
0≤ m≤ p
0>p
(4.4.11)
上式中,系数 hA(l)需用 （ 4.3.6） 式求出 。 因此 AR模型隐含着自相关函数外推的特性 。 我们知道,经典谱估计 BT法中,自相关函数只能限于由观测数据计算出的有限个自相关函数,其它的认为是 0,造成了谱估计分辨率低,模糊 。 也正是 AR模型具有自相关函数外推特性,使它具有高分辨率的优点 。
第四章 功 率 谱 估 计
4.5 AR谱估计的方法
4.5.1 自相关法 ——列文森 （ Levenson） 递推法自相关法的出发点是选择 AR模型的参数使预测误差功率最小，预测误差功率为
n
p
i
pi
n
inxanx
N
ne
N 1
22 )()(|1|)(|1?
(4.5.1a)
假设信号 x(n)的数据区在 0≤n≤N-1范围,有 p个预测系数,N个数据经过冲激响应为 api(i=0,1,2,…,p)的滤波器,输出预测误差
e(n)的长度为 N+P,因此应用下式计算,
第四章 功 率 谱 估 计
2
1
1
0
2
1
0
|)()(|1|)(|1
p
i
pi
PN
n
PN
n
inxanxNneN?
(4.5.1b)
显然,e(n)的长度长于数据的长度,上式中数据 x（ n） 的两端需补充零点,这相当于无穷长的信号经过加窗处理,得到长度为
N的数据 。 用 （ 4.5.1b） 式对系数 api的实部和虚部求微分的方法使预测误差功率最小,得到
)(?
)2(?
)1(?
)0(?)2(?)1(?
)2(?)0(?)1(?
)1(?)1(?)0(?
2
1
pr
r
r
a
a
a
rprpr
prrr
prrr
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
(4.5.2)
第四章 功 率 谱 估 计式中自相关函数采用有偏自相关估计,即
)(?
)()(
1
)(?
*
*
1
0
mr
mnxnx
Nmr
xx
mN
nxx
m=0,1,2,…,p
m=-p+1,-p+2,…,-1 (4.5.3)
对比 （ 4.3.7） 式,（ 4.5.2） 式就是已推导出的 Yule-Walker方程,因此自相关法也是基于解 Yule-Walker方程的一种方法 。
首先由信号的观测数据估计出其自相关函数,再解该方程,得到模型参数,便可求出信号的功率谱 。
Yule-Walker法 。 但是直接解该方程,需要计算逆在矩阵,不方便 。 在第三章自适应滤波器中,曾介绍了基于 Yule-Walker方程中自相关矩阵的性质,导出 Levenson-Durbin递推法,这是一种高效的解方程方法 。 Levenson[CD*2]Durbin
递推法简化重写如下：
第四章 功 率 谱 估 计
2
1
1
1
,1 /)()(?
k
k
l
xxlkxxkk lkrakra?
（ 4.5.4）
*,1,1 ikkkkikki aaaa i=1,2,3,…,k-1 （ 4.5.5）
2 122 )||1( kkkk a （ 4.5.6）
由 k=1开始递推,递推到 k=p,依次得到 {a11,σ 21},{a21,a22,
σ 22},…,{ap1,ap2,…,app,σ 2p}。 AR模型的各个系数以及模型输入白噪声方差求出后,信号功率谱用下式计算,2
1
j-
22j ω2j ω
e1
1
|)e(|)e(
P
i
i
i
wwxx
a
HP
（ 4.5.7）
第四章 功 率 谱 估 计
（ 4.5.6） 式表明：,说明随着阶数增加,
预测误差功率将减少或者不变,为此要求 |akk|≤1,akk称为反射系数 。 另外,递推公式提供了一种确定模型阶数的实验方法,如模型的阶数不知道,由低阶开始递推,当递推到 M阶时,预测误差满足允许的值,停止递推,选 AR模型的阶数为 M。 这种递推法效率高,且当阶数变化时,无需从头计算 。
如果知道信号的 N个观测数据 （ x(n),0≤n≤N-1 ）,利用列文森递推法计算功率谱的计算流程图如图 4.5.1所示 。 图中采用有偏自相关估计 （ 4.5.3） 式计算,σ2是要求的方差 。
2232221 p
)(? mrxx
第四章 功 率 谱 估 计图 4.5.1 利用列文森递推法计算功率谱的流程图输入,x ( n ),n ＝ 0,1,2,3,…,N － 1,?
2
计算：
)0(?1),0(?/)1(?
2
1,1
2
11,1 xxxxxx
rarra?
)(? mr
xx
p =2
2
1
2
2
*
,1,1
2
1
1
1
,1
1
1,,2,1,0
)(?)(?
pppp
kppppkppk
p
p
i
xxipxxpp
a
pkaaaa
iprapra
22
p
p
i
i
pi
wxx
a
P
1
j
2j
e1
1
)(e
结束
p = p + 1
N
Y
第四章 功 率 谱 估 计
4.5.2 协方差法与修正协方差法
1,协方差法这种方法和自相关法一样,仍然利用使预测误差功率最小的方法求模型参数,但由观测数据求预测误差功率的公式如下式：
1 1
2
1
2 |)()(|1|)(|1
N
pn
N
pn
p
k
pk knxanxpNnepN? （ 4.5.8）
将该式对比自相关法中求预测误差功率的公式 （ 4.5.1b）,不同的是求和限不同 。 该公式中使用的观测数据均已得到,不需要在数据两端补充零点,因此比较自相关法去掉了加窗处理的不合理假设 。 为求得模型参数仍然应用复梯度法使 （ 4.5.8） 式达到最小,公式如下：
第四章 功 率 谱 估 计
)0,(
)0,2(
)0,1(
),()2,()1,(
),2()2,2()1,2(
),1()2,1()1,1(
2
1
pc
c
c
a
a
a
ppcpcpc
pccc
pccc
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
（ 4.5.9）
式中
1
* )()(1),(
N
pn
xx knxjnxpNkjc
（ 4.5.10）
白噪声的方差为
1
m i n
2 ),0()0,0(
N
pk
xxpkxxw kcac （ 4.5.11）
第四章 功 率 谱 估 计由观测数据 x(n)(n=0,1,2,…,N-1),利用上面三个公式可以求出模型的参数,{api(i=1,2,3,…,p); σ2w}。 按照定义,
（ 4.5.9） 式中的 cxx(j,k)可以称作协方差函数,它有两个变量,
因此也适合于非平稳随机信号 。 （ 4.5.9） 式中的协方差矩阵是埃尔米特 （ Hermitian） 矩阵,cxx(k,j)=c*xx(j,k),是半正定的 。
这种方法近似于自相关法 。 一些实验结果说明它的分辨率优于自相关法 ［ 5］,另外对于纯正弦信号数据,可以有效地估计正弦信号的频率 。 进一步的内容可参考文献 ［ 14］,［ 15］
第四章 功 率 谱 估 计
2,修正协方差法修正协方差法使用前向和后向预测误差平均值最小的方法,
估计 AR模型的参数,进而估计信号的功率谱 。
Nuttall 在 1976年提出,称为前向 —后向法；同年,Ulrych
和 Ctayton也独立提出,称之为最小二乘法,因此该法也称为前后向线性预测误差功率最小的最小二乘法 。
前面我们已推导出信号的前向和后向预测,分别重写如下：
p
k
pk
p
k
pk
knxanx
knxanx
1
*
1
)()(?
)()(? (4.5.12)
(4.5.13)
第四章 功 率 谱 估 计式中 apk 是 AR模型的参数 。 最小预测误差平均功率是模型输入白噪声的方差,即 ρ p=σ 2w,前后向预测误差平均功率为
)(5.0 pbpep (4.5.14)
式中前向和后向预测误差功率 ρ pe,ρ pb分别用下式表示：
2
1
0 1
*
2
1
1
)()(
1
)()(
1
pN
n
p
k
pkpb
N
pn
p
k
pkpe
knxanx
pN
knxanx
pN
(4.5.15)
(4.5.16)
第四章 功 率 谱 估 计和协方差法一样,上式中仅对那些用到的观测数据预测误差求和 。 为了使预测误差平均功率最小,求 ρp对 apk(k=1,2,3,…,p)
的实部和虚部的微分,或者用复梯度法求,得到
pl
lnxknxanx
lnxknxanx
pNa
p
k
pk
pN
n
N
pn
p
k
pk
pl
p
,,3,2,1
0)]()()(
)()()([
1
1
**
1
0
1
*
1
(4.5.17)
第四章 功 率 谱 估 计经过简化，得到
1
*
1
0
*
1
1
*
1
0
*
)()()()(
)()()()(
N
pn
pN
n
p
k
N
pn
pN
n
pk
lnxnxlnxnx
lnxknxlnxknxa
令
1
*
1
0
* )()()()(
)(2
1),( N
pn
pN
n
xx knxjnxknxjnxpNkjc
(4.5.18)
第四章 功 率 谱 估 计将上式写成矩阵形式：
)0,(
)0,2(
)0,1(
),()2,()1,(
),2()2,2()1,2(
),1()2,1()1,1(
2
1
pc
c
c
a
a
a
ppcpcpc
pccc
pccc
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
(4.5.19) 白噪声的方差估计值为
p
k
xxpkxx
pN
n
p
k
pk
N
pn
p
k
pkpw
kcac
nxknxanx
nxknxanx
pN
1
1
0 1
**
1
*
1
m i n,
2
),0()0,0(
)]()()((
)()()(([
)(2
1
(4.5.20)
第四章 功 率 谱 估 计例 4.5.1 已知信号的四个观察数据为 x(n)={x(0),x(1),x(2),
x(3)}={2,4,1,3},AR （ 1）
模型参数 。
解 （ 1） 自相关法：
5.0
09)3()41(4)24(2
0)1()(
2
1)(
)(
2
1
)1()()(|)(|
4
1
11
11111111
4
011
4
011
11
4
0
2
a
aaaa
nxne
a
ne
ne
a
nxanxnene
nn
n
第四章 功 率 谱 估 计
(2) 协方差法,
7 1 4.0
0)3()41(4)24(2
0)1()(
3
2)(
)(
3
2
)1()()(|)(|
3
1
|)(|
1
11
111111
3
111
3
111
11
3
1
2
1
2
a
aaa
nxne
a
ne
ne
a
nxanxnenene
pN
nn
n
N
pn
第四章 功 率 谱 估 计
4.5.3 伯格（ Burg）递推法设信号 x(n)观测数据区间为,0≤ n≤ N-1,前向,后向预测误差功率分别用 ρ p,e和 ρ p,b表示,预测误差平均功率用 ρ p表示,公式分别为
)(
2
1
|)(|
1
|)(|
1
,,
1
2
,
1
2
,
bpepp
N
pn
b
pbp
N
pn
f
pep
ne
pN
ne
pN
(4.5.21)
(4.5.22)
(4.5.23)
第四章 功 率 谱 估 计下面解释 （ 4.5.21）,(4.5.22)式中的求和限 。 我们已经知道,
前向,后向预测误差公式分别为
p
k
pk
b
p
p
k
pk
f
p
kpnxapnxne
knxanxne
1
*
1
)()()(
)()()( (4.5.24)
(4.5.25)
在上面两个公式中,信号项的自变量最大的是 n,最小的是
n-p。 为了保证计算范围不超出给定的数据范围,在 （ 4.5.21）
和 (4.5.22)式中,选择求和范围为,p≤n≤N-1。
第四章 功 率 谱 估 计为求预测误差平均功率 ρ p最小时的反射系数 kp，令
0?
p
p
k
(4.5.26)
将第三章中已推导出的前向、后向预测误差递推公式重写如下：
)()1()(
)1()()(
1
*
1
11
neknene
neknene
f
pp
b
p
b
p
b
pp
f
p
f
p
(4.5.27)
(4.5.28)
将上面两个公式带入（ 4.5.21）,(4.5.22)、（ 4.5.23）式中，得到
1
2
1
*
1
2
11 ]|)()(||)1()([|)(2
1 N
pn
f
pp
f
p
b
pp
f
pp nekneneknepN?
第四章 功 率 谱 估 计利用复梯度法，求 ρ p对 kp的实部和虚部的微分，并令其结果等于 0，得到
0)]()()1(
)]1()1()([
)(
1
1
*
1
*
1
1
*
111
nenekne
nenekne
pNk
f
p
f
pp
b
p
N
pn
b
p
b
pp
f
p
p
p?
解出 kp，得到
)|)1(||)((|
)1()(2
1
2
1
2
1
1
*
11
N
pn
b
p
f
p
N
pn
b
p
f
p
p
nene
nene
k (4.5.29)
第四章 功 率 谱 估 计上式就是利用伯格递推法求第 p个反射系数的公式 。 将伯格递推法求 AR模型参数的递推公式总结如下：
)|)1(||)((|
)1()(2
1,,3,2,1,0)()(
1,,3,2,1,0)()(
)0(?
|)(|
1
)0(?
2
1
1
2
1
1
*
11
0
0
0
1
0
2
N
pn
b
p
f
p
N
pn
b
p
f
p
p
b
f
xx
N
n
xx
nene
nene
k
Nnnxne
Nnnxne
r
nx
N
r
(4.5.30)
(4.5.31)
(4.5.32)
(4.5.33)
(4.5.34)
第四章 功 率 谱 估 计
2,,1,)()1()(
1,,2,1)()()(
1,,3,2,1
)||1(
1
*
11
11
,
*
,1,1
1
2
Nppnneknene
Nppnneknene
ka
piakaa
k
f
pp
b
p
b
p
b
pp
f
p
f
p
ppp
ipppippi
ppp
(4.5.36)
(4.5.35)
(4.5.37)
(4.5.38)
(4.5.39)
对于 p阶 AR模型的输入白噪声方差 σ 2w=ρ p。 利用伯格递推法求
AR模型参数的流程图如图 4.5.2所示 。
第四章 功 率 谱 估 计图 4.5.2 伯格递推法流程图输入,x ( n ),n ＝ 0,1,2,3,…,N － 1
阶数 IP
1)(
1
)()()()(
2
1
0
0
00
pnx
N
nxnenxne
N
n
bf
p = IP
结束
1 pp
1
2
1
2
1
*
1
1
1
)()(
)1()(2
N
pn
b
p
f
p
b
p
N
pn
f
p
p
nene
nene
k
pwppp
k
2
1
2
1
pppipppippi
kaakaa
*
,1,1
1,,2,1),()1()(
1,,2,1),1()()(
1
*
1
11
Nppnneknene
Nppnneknene
f
pp
b
p
b
p
b
ppp
f
p
2
,;,,3,2,1:
wip
pia输出
2
1
j
2
j
e1
)(e
p
i
i
pi
w
xx
a
P
N
Y
第四章 功 率 谱 估 计
4.5.4 关于 AR模型阶次的选择在 AR模型谱估计中,模型阶次的选择也是一个关键问题 。
一般模型的最好选择是先验未知的,实际中需预先选定模型阶次 。 如果是纯 P阶 AR信号,选择模型阶次 k＜ P时,将产生对谱的平滑作用,降低谱的分辨率,如图 4.5.3所示 。 图中,AR信号 P=4,选择模型的阶次 k=2,产生的平滑作用使两个峰变成一个峰,分辨率明显降低 。 如果选择 k≥ P,且假定观测的数据没有误差 （ 没有干扰 ）,估计的参数应是：
kppi
pia
a ipki
,,2,10
,,3,2,1,
第四章 功 率 谱 估 计图 4.5.3 AR模型阶次太小时的平滑作用
5 0,0 0
3 0,0 0
1 0,0 0
－ 1 0,0 0
－ 3 0,0 0
0,0 0 0,1 0 0,2 0 0,3 0 0,4 0 0,5 0
真实 A R(4 ) 的 PS D
A R(2 ) 的 PS D 模型
P
S
D
/
d
B
/?
第四章 功 率 谱 估 计因此,对于纯 P阶 AR信号,应选择阶次 k≥P。 一般对于纯 AR信号,
只要记录数据的长度不十分短,模型阶次的估计结果比较满意
［ 21］ 。 但是如果是白噪声中的 AR信号 （ 观测数据有误差或者信号中含有白噪声 ）,此时选择 ARMA模型合适,如选择了 AR模型,其阶次应加大,较低的阶次会使谱估计产生偏移,降低分辨率 。 当然,这也和信噪比有关,信噪比愈低,平滑作用愈严重,愈需要高的阶次,因此信噪比低应选高的阶次 。 一般来说,
阶次愈高,分辨率愈高；但阶次太高,会使估计误差加大,谱峰分裂,因此,对于白噪中的 AR信号,其阶次的选择应折衷考虑 ［ 5］ 。 下面介绍三种模型阶次的估计方法 。
第四章 功 率 谱 估 计为了确定模型的阶次,所提出的各种方法都是基于估计出的预测误差功率,即选择模型的阶次,使预测误差功率极小化,
这是因为所讨论的 AR模型参数估计方法,预测误差功率都是随模型的阶次的增加而减少,或者不变的 。 但是不能仅仅考虑预测误差功率的减少,也要考虑因模型阶次增高引起模型参数加多,引起估计误差的加大 。 根据这一原理,阿凯克提出两种估计方法,第一种方法称为最终预测误差 （ FPE） 准则,用下式计算后再取最小,估计模型阶次：
kkN
kNkFPE)(
（ 4.5.40）
第四章 功 率 谱 估 计式中,表示 k阶 AR模型的白噪声方差 （ 预测误差功率 ） 的估计值 ;N是观测数据的长度 ; k是模型阶次 。 上式中,随 k增加而减少,但 (N+k)/(N-k)却随 k增加而加大 。 第二种方法应用更广泛,称为阿凯克信息论准则：
kp?
kp?
knk k 2?ln)(A I C （ 4.5.41）
选择使上式最小的 k值作为模型的阶次 。 但这种方法有过高估计模型阶次的趋势,例如真实的阶次是 2,而 AIC(k)最小值对应的是 4,为改进,提出用 k lnN代替公式中的 2k的方法 ［ 22］,［ 23］ 。
AIC和 FPE估计阶次的性能是相近的,但对于短记录数据,
建议使用 AIC准则 ［ 24］,对于长记录 （ N→∞ ）,两种方法得到相同的模型阶次 。
第四章 功 率 谱 估 计第三种方法由 Parzen在 1976年提出,称为自回归传递函数准则 （ CAT） ［ 25］,定义是：
k
k
i iN
k1~11)(C A T
1
(4.5.42)
式中
ii iN
N~
同样,模型的阶次取使上式最小的 k值 。 以上三种估计模型阶次的方法性能类似 ［ 26］ 。 还有其它的模型阶次估计方法,
可参考文献 ［ 27］ ～ ［ 29］ 。
第四章 功 率 谱 估 计
4.6 最大熵谱估计与最大似然谱估计
4.6.1 最大熵谱估计
1,利用最大熵的原则外推自相关函数按照 Shannon对熵的定义,当随机变量 X取离散值时，熵的定义为
i
ii ppH ln
(4.6.1)
式中 pi是出现状态 i的概率。当 X取连续值时，熵的定义为
)]([ l nd)(ln)( xpExxpxpH
(4.6.2)
第四章 功 率 谱 估 计式中,p(x)是 X的概率密度函数,对于离散随机序列,概率密度函数用联合概率密度函数代替 。 显然,熵代表一种不确定性,最大熵代表最大的不确定性,或者说最大的随机性 。 下面我们研究对于有限的自相关函数值不作任何改变,对于未知自相关函数用最大熵原则外推,即不作任何附加条件的外推方法 。 假设 x(n)
是零均值正态分布的平稳随机序列,它的 N维高斯概率密度函数为
XNRXNRxxxp
xx
H
xx
N
N
12/12/
21 ))((2
1e x p))(( d e t)π2(),,,(?
式中
H21 ],,,[ NxxxX
第四章 功 率 谱 估 计
)0()1()(
)1()0()1(
)()1()0(
)(
xxxxxx
xxxxxx
xxxxxx
xx
rNrNr
Nrrr
Nrrr
NR
按照（ 4.6.2）式,x(n)信号的熵为
])))(( d e t ()πe2l n [ ( 2/12/ NRH xxN? （ 4.6.3）
式中 det(Rxx(N))表示矩阵 Rxx(N)的行列式,由上式表明为使熵最大,要求 det(Rxx(N)最大 。
第四章 功 率 谱 估 计若已知 N+1个自相关函数值 rxx(0),rxx(1),…,rxx(N),下面用最大熵方法外推 rxx(N+1)。 设 rxx(N+1)确实是信号自相关函数的第
N+2个值,根据自相关函数的性质,由 N+2个自相关函数组成的矩阵为
)0()1()()1(
)()1()0()1(
)1()()1()0(
)1(
xxxxxxxx
xxxxxxxx
xxxxxxxx
xx
rrNrNr
NrNrrr
NrNrrr
NR
（ 4.6.4）
第四章 功 率 谱 估 计它必须是非负定的矩阵,即
0)]1(d et [NR xx （ 4.6.5）
0)](d e t [,),1(d e t [),0(d e t [?NRRR xxxxxx?
将行列式展开,det(Rxx(N+1))是 rxx(N+1)的二次函数,该二次函数系数的符号是,(-1)1+N+2(-1)1+N+1=-1,且 det(Rxx(N+1)) 对
rxx(N+1)的二次导数是 -2det［ Rxx(N-1)］,它是负值,负值表示
det (Rxx(N+1))对 rxx(N+1)的一次导数是减函数,det(Rxx·(N+1))作为 rxx(N+1)的函数,凹口向下,那么只有一个最大值 。 为选择
rxx(N+1)使 det(Rxx(N+1)最大,解下列方程：
0)]1(d e t [)1(d d NRNr xx
xx
(4.6.6)
第四章 功 率 谱 估 计用数学归纳法，得到
0
)1()()1(
)2()1()2(
)1()0()1(
xxxxxx
xxxxxx
xxxxxx
rNrNr
Nrrr
Nrrr
（ 4.6.7）
上式是 rxx(N+1) 的一次函数,可以解出 rxx(N+1) 。 继续再将
rxx(N+1) 代入 Rxx(N+2) 和 det(Rxx(N+2)) 中,求 det(Rxx(N+2)) 对
rxx(N+2)的最大值,得到 rxx(N+2); 以此类推,可推出任意多个其它自相关函数值,而不必假设它们为零,这就是最大熵谱估计的基本思想 。
第四章 功 率 谱 估 计
2,最大熵谱估计与 AR模型谱估计的等价性我们已经知道 AR
Yule-Walker方程,即
N
l
wxxA
N
l
xxA
xx
lrlh
lmrlh
mr
1
2
1
)()(
)()(
)(
m≥ 1
m=0
将 m≥ 1的情况写成矩阵形式：
)(
)2(
)1(
)0()2()1(
)2()0()1(
)1()1()0(
2
1
Nr
r
r
a
a
a
rNrNr
Nrrr
Nrrr
xx
xx
xx
N
xxxxxx
xxxxxx
xxxxxx
第四章 功 率 谱 估 计式中,ai=hA(i),i=1,2,3,…,N,ai是 AR模型系数 。 解该方程,可以得到模型系数 ai,即
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
)0()2()(
)2()0()2(
)1()1()1(
1
（ 4.6.8）
第四章 功 率 谱 估 计
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
)0()1()(
)2()1()2(
)1()0()1(
2
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
N
N
)1()1()(
)3()1()2(
)2()0()1(
)1(
（ 4.6.9）
（ 4.6.10）
第四章 功 率 谱 估 计
)0()()1(
)2()2()1(
)1()1()0(
xxxxxx
xxxxxx
xxxxxx
rNrNr
Nrrr
Nrrr
D
(4.6.11)
在（ 4.3.6）式中，令 m=N+1，得到
N
l
xxlxx lNraNr
1
0)1()1(
(4.6.12)
将以上求出的系数 a1,a2,…,aN代入上式,求出 rxx(N+1)。 而最大熵外推自相关函数的公式是 (4.6.7)式,按照该公式的最后一行展开,得到第四章 功 率 谱 估 计
0)1()1(
0)1()1()()1(
0)1()1()()1(
1
21
21
kNraNr
araNraNrNr
DarDaNrDaNrDNr
xx
N
k
kxx
Nxxxxxxxx
Nxxxxxxxx
（ 4.6.13）
上式即是最大熵外推自相关函数的公式,对比 （ 4.6.12） 式,
两公式完全一样,证明了 AR模型功率谱估计和最大熵谱估计的等价性 。 这里最大熵外推自相关函数等价于已知 N+1个自相关函数,匹配一个 N阶 AR信号模型的系数 。 一旦通过解 Yule-Walker
方程,解出模型参数,最大熵谱估计用下式计算信号功率谱：
N
k
kj
k
wj
xx
a
P
1
2
e1
)e(
（ 4.6.14）
第四章 功 率 谱 估 计
4.6.2 最大似然谱估计 ——最小方差谱估计最大似然谱估计是用一个 FIR滤波器实现,该滤波器对所关心频率的正弦信号,可以无失真地通过,而对于其它频率的信号,让其频响尽可能地小,亦即将它们尽可能地滤除 。 此时,
滤波器输出的均方值,就作为信号的功率谱估计 。
设实信号用 x(n)表示,FIR滤波器系统函数用 A(z)表示：
p
n
nznazA
0
)()(
p
k
pAXknxkananxny
0
T)()()()()(
输出 y(n)为
(4.6.15)
第四章 功 率 谱 估 计式中
T
T
)](,),1(),([
)](,),1(),0([
pnxnxnxX
paaaA p
输出信号的均方值为
ppHppTHppTHpy ARAAXXEAAXXAEnyEP ][][]|)([|? 2
(4.6.16)
上式中 T表示转置,H表示共轭转置,Rp=E［ XXT］ 是 Toeplith
自相关矩阵,为求,必须先求 FIR滤波器的系数 。 求这些系数的原则是：在所关心频率 ωi处,信号 x(n)无失真地通过,
即在 ω i处的传输函数为 1：
yP?
第四章 功 率 谱 估 计
pi
H
p
n
p
n
AEA
naA
i
ii
)()(e
1e)()e(
j
j-
0
j
式中
Hjj2j ]eee1[ pωωωH
p iiiE
(4.6.17)
另外一个原则是在 ωi附近的频率分量尽量衰减掉,即 ω≠ωi处,滤波器输出 y(n)的均方差 最小,即 (4.6.16)式最小,此时 作为信号 x(n)的功率谱估计 。 因此,最大似然谱估计称为最小方差谱估计更为合适,但由于习惯也可以仍称为最大似然谱估计 。 在以上原则下,使方差 最小的滤波器系数和 分别为
［ 30］,［ 31］
yP? yP?
xP?
yP? x
P?
第四章 功 率 谱 估 计
)()(
1
)()(
)(
1
1
1
pi
H
p
yx
ipi
H
p
ipp
p
RE
PP
RE
ER
A
应该指出,此时 并不是真正意义上的信号功率谱,只是描述了信号功率谱的相对强度 。
yP?
第四章 功 率 谱 估 计下面分析最小方差谱估计与 AR模型谱估计之间的关系：
1)2()1()(
01)1()2(
001)1(
0001
1
1
papapa
aa
a
A
ppp
pp
p
pp
（ 4.6.20）
第四章 功 率 谱 估 计
2
2
1
2
0
0
0
p
w
P
(4.6.21)
第四章 功 率 谱 估 计
4.7 特征分解法谱估计
4.7.1 正弦波用退化 AR模型表示无论是实正弦波还是复正弦波,都可以用一个退化 AR模型表示,设 P个实正弦波组成的信号用下式表示：
)(s in)(
1
ii
P
i
i nqnx
（ 4.7.1）
式中,初相位 θ i是在区间 （ -π,π ） 均匀分布的随机变量,
首先分析下面的三角恒等式：
])1(s i n [c o s2])2(s i n [)s i n ( nnn
-π ＜ ω ＜ π
第四章 功 率 谱 估 计令 x(n)=sin(ω n+θ ),则上式变为
)2()1(c o s2)( nxnxnx? （ 4.7.2）
将上式进行 Z变换，得到
0)co s21)(( 21 zzzX? （ 4.7.3）
这样 (4.7.2)式的特征多项式为
0c o s21 21 zz? （ 4.7.4）
第四章 功 率 谱 估 计上式的两个根分别是,z1=ejω,z2=e-jω,它们共轭成对,且模为 1。
由这两个根可以确定正弦波的频率 。 对比 AR模型的系统函数,
可以把正弦波信号用一个特殊的 AR（ 2） 模型表示,括弧中的 2
表示模型是二阶的 。 该 AR模型的激励白噪声方差趋于 0,极点趋于单位圆 。 通常称为退化的 AR模型 。 这一模型系数有两个,即
2 cosω 和 1,（ 4.7.2） 式是模型的差分方程 。
第四章 功 率 谱 估 计对于 P个实正弦波,特征多项式是
0)c o s21(
1
21
p
i
i zz?
（ 4.7.5）
上式是 z-1的 2P阶多项式，可以表示为
10 0
2
0
aza
P
k
k
k
（ 4.7.6）
注意上式中的系数 ak(k=1,2,3,…,2P),必须保证它的根共轭成对 。
考虑到根共轭成对,也可表示为
P
i
ii
P
k
k
k zzzzza
1
*
2
0
))((
（ 4.7.7）
第四章 功 率 谱 估 计这样由 (4.7.6)式,P个正弦波组合的模型用下面 2P阶差分方程描述
P
k
k knxanx
2
1
)()(
（ 4.7.8）
对于复正弦波情况,P个复正弦波组成的信号是
P
i
n
i
iiqnx
1
)j(e)(
（ 4.7.9）
用一个退化的 AR（ p）模型表示的差分方程为
p
i
k n - kxanx
1
)()(
（ 4.7.10）
第四章 功 率 谱 估 计其特征多项式为
10 0
0
aza
P
k
k
k
(4.7.11)
其根为
ijiz?e?
1≤ i≤ P
注意这里的根不是共轭成对出现的 。
总结以上 P个正弦波组合是一个退化的 AR(2P)过程,独立参量个数为 P个； P个复正弦波的组合是退化的 AR（ P） 过程,独立参量个数仍为 P个 。 实正弦过程相应的退化 AR过程的阶数比复正弦情况的阶数高 1倍 。
第四章 功 率 谱 估 计
4.7.2 白噪声中正弦波组合用一特殊的 ARMA模型表示白噪声中正弦波组合的信号为
)()s in ()()()(
1
nwnqnwnxny
P
i
iii
（ 4.7.12）
式中,w(n)为白噪声,且
0)]()([,)]()([,0)]([,2 lwnxElwnwEnwE lnw
将 (4.7.12)式中 x(n)的用 AR(2P)表示，即将 (4.7.8)式带入 (4.7.12)
式中,得到
)()()(
2
1
nwinxany
P
i
i
（ 4.7.13）
第四章 功 率 谱 估 计将（ 4.7.12）式中的 n用 n-i代替,
x(n-i)=y(n-i)-w(n-i)
再将上式带入 (4.7.13)式,得到
P
i
i
P
i
i inwanwinyany
2
1
2
1
)()()()(
(4.7.14)
上式可以看成一个特殊的 ARMA(2P,2P)模型,括弧中的两个 2P分别表示 ARMA模型系统函数分子和分母的阶次 。 它与一般的
ARMA模型比较,有三方面不同：
第四章 功 率 谱 估 计
(1) 它的 AR部分和 MA部分具有相同的参数,它们存在共同的因子；
(2) 由于特征多项式 (4.7.6)式的根的模为 1,故 AR部分特征多项式不满足平稳性条件,MA部分特征多项式也不满足可逆性条件；
（ 3） AR部分的 y(n)=x(n)+w(n),y(n)是含白噪声的观测值,
而通常为信号的 x(n)不含白噪声 。
第四章 功 率 谱 估 计
4.7.3 特征分解法谱估计这种特殊的 ARMA模型结构,不能用一般的 ARMA模型结构求解 。 下面介绍特征分解技术 。
将（ 4.7.14）式写成矩阵形式：
YTA=WTA (4.7.15)
式中
T
T
221
T
)]2(,),1(),([
],,,,1[
)]2(,),1(),([
pnwnwnwW
aaaA
pnynynyY
p
第四章 功 率 谱 估 计用向量 Y左乘 (4.7.15)式,并取数学期望,得到
E［ YY T］ A=E［ YW T］ A (4.7.16)
式中
yy
yyyyyy
yyyyyy
yyyyyy
R
rprpr
prrr
prrr
YYE?
)0()12()2(
)12()0()1(
)2()1()0(
][
T
IWWEWWXEYWE w2TTT ][])[(][
第四章 功 率 谱 估 计将上面关系式带入 (4.7.16)式,得到
AAR wyy 2
(4.7.17)
式中,σ w2是自相关函数 Ryy的特征值 ; A是对应 σ w2的特征矢量 。
由于 x(n)与 w(n)互不相关,由 （ 4.7.12） 式求 y(n)的自相关函数,
得到
IRR wxxyy 2
式中 Rxx是 x(n)的自相关函数,可以推导出在特征方程 (4.7.17)
式中,σ w2是 Ryy的最小特征值,且 Ryy的阶数超过 (2P+1)(2P+1)时,
σ w2就是其多重最小特征值 。 这一结论为我们寻求向量 A提供了重要的依据 。
(4.7.18)
第四章 功 率 谱 估 计当特征向量 A求出后,就可以通过解特征方程,解出各个根,求出各正弦波的频率值,特征方程为
01 222211
2
0
P
P
P
k
k
k zazazaza?
（ 4.7.19）
该方程有 2P个根，这些根在 z平面的单位圆上，它们是
iωiz je?
i=1,2,3,…,2P （ 4.7.20）
式中的 ω i即是正弦波的频率 。 当各正弦频率由 (4.7.19)式求出后,各正弦波的功率也可以求出,对于白噪声中 P个实正弦组合信号,它的自相关函数为
P
i
iwyy Pr
1
2)0(? （ 4.7.21）
第四章 功 率 谱 估 计
mPmr i
P
i
iyy?c o s)(
1
m≠ 0 (4.7.22)
式中
2
2
i
i
qP?
qi是频率为 ωi的正弦波的幅度,Pi是其功率。
将 (4.7.22)式写成矩阵形式,
FP=r (4.7.23)
式中
P
P
P
PPP
F
c o sc o sc o s
2c o s2c o s2c o s
c o sc o sc o s
21
21
21
(4.7.24)
第四章 功 率 谱 估 计
PT=［ P1,P2,…,PP］
)](,),2(),1([T Prrrr yyyyyy
(4.7.25)
(4.7.26)
由于各正弦波的频率已求出,矩阵 F,r已知,由 (4.7.23)式解出各正弦波的功率或幅度 。 最后噪声功率由下式求出,
p
i
iyyw Pr
1
2 )0(?
(4.7.27)
以上就是皮萨论科谱分解法的全过程 。 其中关键的一步是求自相关函数 Ryy的最小特征值及对应的特征向量 。 求最小特征值对应的特征向量经典方法是幂法 ［ 32］,可以利用 Ryy的布利斯结构,
采用升幂法迭代公式求解 A,记第 k次迭代值为 A(k),迭代公式为
)()1( 1 kARkA yy
(4.7.28)
第四章 功 率 谱 估 计上面公式将收敛到 Ryy的最小特征值对应的特征矢量 A,A的理想初值是
T]1,,1,1[)0(A
(4.7.29)
通常只需迭代几次就能求得所需的稳定特征矢量 A（ ∞ ）,一旦求得 A,由 （ 4.7.17） 式求出 Ryy的最小特征值 σ 2w:
AA
ARA
T
yy
T
w?
2?
(4.7.30)
第四章 功 率 谱 估 计此外,皮萨论科方法需要精确知道自相关函数矩阵,但在实际中只能由 y(n)的观测数据进行估计,估计误差会影响估计频率的精确性,为此,玛坡尔 (Marple)提出了改进方法,这方面内容可参考文献 ［ 32］ 。 在处理实际问题时,通常不知道正弦波的个数,即模型的阶数 2p是不知道的 。 如果 Ryy的维数是
(2p+1)(2p+1),最小特征值是 σ 2w,假如增大 Ryy的维数,它的最小特征值重复数目增加,但都是 σ 2w 。 这样可以逐步增加阶数,直到估计出相邻的最小特征值相等时,就认为得到了模型的正确阶数 。 为了清楚地了解皮萨论科方法的计算过程,
具体计算步骤归纳如下：
第四章 功 率 谱 估 计
(1) 由 y(n)的观测数据估计出 2p+1个相关函数 ryy(m),m=0,
1,2,…,2p值,给定阶数 p初始值 。
(2) 求 Ryy的最小特征值和特征向量 。 可以采用 （ 4.7.28）
~(4.7.30)式进行迭代计算 。
(3) 进行最小特征值判断 。 如此次求得的最小特征值和上一个 p值时求得的相同,则认为上次的阶数 p即为所需的阶数,否则继续增加阶数进行计算 。
(4) 由 （ 4.7.19） 式求特征多项式的根 。
(5) 由求得的特征多项式的根确定各正弦波的频率 。
(6) 由方程（ 4.7.23）和 (4.7.27)式,求各正弦波的功率和噪声功率。
4.1 引言
4.2 经典谱估计
4.3 现代谱估计中的参数建模
4.4 AR模型谱估计的性质
4.5 AR谱估计的方法
4.6 最大熵谱估计与最大似然谱估计
4.7 特征分解法谱估计第四章 功 率 谱 估 计
4.1 引 言我们知道,对信号和系统进行分析研究,处理有两类方法：一类是在时域进行,前面我们学习的维纳 —卡尔曼滤波和自适应滤波都属于这种方法；本章则是在频率域进行研究的另一类方法 。
这两类方法都是信号处理的重要方法 。 对确定性信号傅里叶变换是在频率域分析研究的理论基础,但对于随机信号,其傅里叶变换并不存在,因此转向研究它的功率谱 。 按照 Weiner-
Khintchine定理,信号的功率谱和其自相关函数服从一对傅里叶变换关系,公式如下
mj
m
xx
j
xx mrP
e)()e(
(4.1.1)
第四章 功 率 谱 估 计
dPmr njjxxxx e)e(π2 1)(
)]()([)( * mnxnxEmr xx
(4.1.2)
(4.1.3)
（ 4.1.1） 式被称做功率谱的定义,对于平稳随机信号,服从各态历经定理,集合平均可以用时间平均代替,由 （ 4.1.1）
式还可以推出功率谱的另一个定义,推导如下：
将 （ 4.1.3） 式中的集合平均用时间平均代替,得到
)()(12 1li m)( * mnxnxNmr
N
NnN
xx
(4.1.4)
第四章 功 率 谱 估 计将（ 4.1.4）式代入 (4.1.1)式,得到
m
mnjnj
N
Nn
N
nj
m
N
m
j
xx
mnxnx
N
mnxnx
N
P
)(*
*
e)(e)(
12
1
lim
e)()(
12
1
lim)e(
令 l=n+m,则
2
*
j
e)(
12
1
lim
e)()(
12
1
lim)e(
j ωω
N
Nn
N
lj
l
N
Nn
nj
N
xx
nx
N
lxenx
N
P
(4.1.5)
第四章 功 率 谱 估 计上式中 x(n)是观测数据,Pxx(ejω)是随机变量,必须对 Pxx(ejω)
取统计平均值,得到
2
j ω )(
12
1lim)e( N
Nn
nj
Nxx
enx
N
EP?(4.1.6)
上式被认为是功率谱的另一定义 。
（ 4.1.1） 式表明功率谱是无限多个自相关函数的函数,但观测数据只有有限个,只能得到有限个自相关函数 。 按照
（ 4.1.6） 式求功率谱,也需要无限个观测数据 。 因此根据有限个样本数据,分析计算随机序列的真正功率谱,是求功率谱的中心问题,毫无疑问,这是一个功率谱的估计问题 。 在第一章已介绍了统计估计的一般估计准则,主要有偏移,估计量方差和估计量的均方误差 （ 有效性 ）,这里不再重复,下面直接用它们分析估计质量 。
第四章 功 率 谱 估 计现代谱估计以信号模型为基础,图 4.1.1表示的是 x(n)的信号模型,输入白噪声 w(n)均值为 0,方差为 σ 2w,x(n)的功率谱由下式计算：
2j2j |)e(|)e( HP wxx?
（ 4.1.7）
如果由观测数据能够估计出信号模型的参数,信号的功率谱可以按照 （ 4.1.7） 式计算出来,这样,估计功率谱的问题变成了由观测数据估计信号模型参数的问题 。 模型有很多种类,例如
AR模型,MA 模型等等,针对不同的情况,合适地选择模型,
功率谱估计质量比较经典谱估计的估计质量有很大的提高 。 遗憾的是,尚无任何理论能指导我们选择一个合适的模型,我们只能根据功率谱的一些先验知识,或者说一些重要的谱特性,
选择模型 。
第四章 功 率 谱 估 计
H ( z )
w ( n ) x ( n )
图 4.1.1 平稳随机序列的信号模型第四章 功 率 谱 估 计
4.2 经 典 谱 估 计
4.2.1 BT法
BT法是先估计自相关函数,然后按照 (4.1.1)式进行傅里叶变换得到功率谱 。 设对随机信号 x(n),只观测到一段样本数据,
n=0,1,2,…,N-1。 关于如何根据这一段样本数据估计自相关函数,第一章已经作了详细介绍,结果是共有两种估计方法,即有偏自相关函数估计和无偏自相关函数估计 。 有偏自相关函数估计的误差相对较小,这种估计是一种渐近一致估计,将该估计公式重写如下：
1||
0
* )()(1)(?
mN
n
xx mnxnxNmr
（ 4.2.1）
第四章 功 率 谱 估 计对上式进行傅里叶变换，得到 BT法的功率估计值为
m
j ωω
xx
j mrP -
BT e)(?)e(
(4.2.2)
为了减少谱估计的方差,经常用窗函数 w(m)对自相关函数进行加权,此时谱估计公式为
j ωω
xx
Mm
j mrP -
)1(
BT e)(?)e(
（ 4.2.3）
式中
0
)()( mwmw
-(M-1)≤ m≤ (M-1)
其它,M≤ N
（ 4.2.4）
第四章 功 率 谱 估 计有时称 (4.2.3)式为加权协方差谱估计 。 它要求加窗后的功率谱仍是非负的,这样窗函数 w(m)的选择必须满足一个原则,即它的傅里叶变换必须是非负的,例如巴特利特窗就满足这一条件 。
为了采用 FFT计算 （ 4.2.3） 式,设 FFT的变换域为 (0~L-1)，
必须将求和域 (-M+1,M-1)移到 （ 0~L-1）,功率谱的计算公式如下：
)()(?
e)()e(?
BT
1
0
j-j ω
BT
mSFFTkP
mSP
xx
L
m
ωm
xx
k=0,1,2,…,L-1
第四章 功 率 谱 估 计
(4.2.7)
)()(?
0
)()(?
)(
LmwLmr
mwmr
mS
xx
xx
xx
0≤ m≤ M-1
M≤ m≤ L-M-1
L-M≤ m≤ L-1
按照 （ 4.2.1） 式估计自相关函数,我们已经证明这是渐近一致估计,但经过傅里叶变换得到功率谱的估计,功率谱估计却不一定仍是渐近一致估计,可以证明它是非一致估计,是一种不好的估计方法 。 下面我们将证明,BT法中用有偏自相关函数进行估计时,它和用周期图法估计功率谱是等价的,因此 BT
法估计质量和周期图法的估计质量是一样的 。
第四章 功 率 谱 估 计
4.2.2 周期图法将功率谱的另一定义（ 4.1.6）式重写如下,
2
jj e)(
12
1
lim)e(
N
Nn
n
Nxx
nx
N
EP
如果忽略上式中求统计平均的运算,观测数据为,x(n) 0≤n≤N-1，
便得到周期图法的定义,
21
0
j-j e)(1)e(
N
n
n
xx nxNP
(4.2.8)
第四章 功 率 谱 估 计图 4.2.1 用周期图法计算功率谱框图观测数据 x ( n ) FFT
取模的平方
1/ N
P xx (e j )
^
第四章 功 率 谱 估 计
1,周期图与 BT法的等价关系周期图法的功率谱估计公式用 （ 4.2.8） 式表示,下面由该公式出发推导它们的等价关系 。
1
0
1
0
)(*
1
0
*
1
0
2
1
0
e)()(
1
e)(e)(
1
e)(
1
)e(?
N
k
N
n
nkj
N
n
nj
N
k
nj
N
n
njj
xx
nxkx
N
nxkx
N
nx
N
P
令 m=k-n,即 k=m+n，则
)(
1
0
*j e)()(1)e(? nkj
mN
n
xx nmxnxNP
第四章 功 率 谱 估 计上式中的方括号部分正是有偏自相关函数的计算公式,因此得到
mj
xx
N
Nm
xx mrP
)e(?)e(?
1
)1(
j
)()(1)(? *
||1
0
nmxnxNmr
mN
n
xx
)e(?)e(? BT jjxx PP?
因此证明了利用有偏自相关函数的 BT法和周期图法的等价关系。
第四章 功 率 谱 估 计
2,周期图法谱估计质量分析
1)
已知自相关函数的估计值,m=-(N-1),-N,-N+1,…,0,
1,2,…,N-1,按照 （ 4.2.2） 式求功率谱的统计平均值,得到
)(? mrxx
nxxN
Nm
xx mrEPE
j-
1
)1(
j e)(?)]e(?[?
有偏自相关函数统计平均值已由第一章 (1.3.30)式确定,将该式代入上式,得到
n
xx
N
Nm
xx mrN
mNPE j-1
)1(
j )e(||)]e(?[
第四章 功 率 谱 估 计
m
n
xxBxx mrmwPE
j-j )e()()]e(?[ （ 4.2.9）
式中
其它0
||
||
)(
Nm
N
mN
mw B （ 4.2.10）
(4.2.9)式中,两序列乘积的傅里叶变换,在频域服从卷积关系，
得到
1)]e(?[ jxxPE
ΩWPE jBjxx d)e(π2 1)]e(?[ )(（ 4.2.11）
式中
)]([)e( mrFTP xxjxx
第四章 功 率 谱 估 计
2
)2/s i n (
)2/s i n (1)]([)e(
N
N
mwFTW BjB
(4.2.12)
WB(ejω)称为三角谱窗函数 。 （ 4.2.11） 式表明,周期图的统计平均值等于它的真值卷积三角谱窗函数,因此周期图是有偏估计,但当 N→∞ 时,wB(m)→ 1,三角谱窗函数趋近于 δ 函数,周期图的统计平均值趋于它的真值,因此周期图属于渐近无偏估计 。
第四章 功 率 谱 估 计
2)
由于周期图的方差的精确表示式很繁冗,为分析简单起见,
通常假设 x(n)是实的零均值的正态白噪声信号,方差是 σ x2，
即功率谱是常数 σ x2,其周期图用 IN(ω )表示,N表示观测数据的长度 。 按照周期图的定义,周期图表示为
)]([)]([)](v a r [
ee)()(
1
|)e(|
1
)(
22
j-
1
0
1
0
j2j
NNN
n
N
n
N
k
k
N
IEIEI
nxkx
N
X
N
I
下面先求周期图的均值，再求其均方值：
)(j-e)()()]()([1)]([ n - k
N
k
N
n
N kRnRkxnxENIE
第四章 功 率 谱 估 计式中
kn
kn
kn
knknrkxnxE xxx
0
1
)(
)()()]()([
2
222 )(1)]([
x
n
xNN nRNIE
（ 4.2.13）
上式说明周期图是无偏估计,但前面已推导出周期图是有偏估计
（ 一般情况 ）,这里由于对信号作了实白噪声的假设,才有无偏估计的结果 。 在求均方值时,先求两个频率 ω 1和 ω2处的均方值,
最后令 ω =ω 1=ω 2。
第四章 功 率 谱 估 计
)(j-)(j-
2
2
j
2
j
221
21
21
ee)]()()()([
)()()()(
1
)e()e(
1
)]()([
qpkn
NNN
q
N
pkn
NNNN
qxpxkxnxE
qRpRkRnR
N
XX
N
EIIE
利用正态白噪声、多元正态随机变量的多阶矩公式，有
)]()([)]()([
)]()([)]()([)]()([)]()([)]()()()([
pxkxEqxnxE
qxkxEqxpxEqxpxEkxnxEqxpxkxnxE
其它0
,;,,,
)]()()([
4 kpnqkqnpqpnk
qxkxnxE x
第四章 功 率 谱 估 计将上式代入周期图的均方值公式中,得到
2
21
21
2
21
214
1
0
1
0
1
0
))((
1
0
))((2
2
4
21
2/)s i n (
2/)s i n (
2/)s i n (
2/)s i n (
1
ee)]()([ 2121
N
N
N
N
N
N
IIE
x
N
n
N
n
N
k
n - kj
N
k
n - kjx
NN
(4.2.14) 将 ω =ω 1=ω 2代入上式，得到
2
4
42
)s i n (
)s i n (
1)](v a r [
)s i n (
)s i n (
2)]([
N
N
I
N
N
IE
xN
xN
(4.2.15)
第四章 功 率 谱 估 计显然,当 N趋于无限大时,周期图的方差并不趋于 0,而趋于功率谱真值的平方,即
4)](v a r [ x
NNI
(4.2.16)
这里无论怎样选择 N,周期图的方差总是和 σ 4x同一个数量级 。 我们知道,信号的功率谱真值是 σ 2x,说明周期图的方差很大,周期图的均方误差也是非常大 。 用这种方法估计的功率谱在 σ 2x附近起伏很大,故周期图是非一致估计,是一种很差的功率谱估计方法 。
第四章 功 率 谱 估 计为了进一步说明数据长度 N对功率谱估计的影响，下面求两个频率处的协方差函数。
)]()([)]()([)](),(c o v [ 212121 NNNNNN IIEIIIII
将（ 4.2.13）式和 (4.2.14)式代入上式,得到
2
21
21
2
21
214
21 2/)s in (
2/)(s in
2/)s in (
2/)(s in)]()(c o v [
N
N
N
NII
xNN
令,ω 1=2π k/N,ω 2=2π l/N,式中 k,l均是整数，得到
224
21 ]/)(πs in [
)](πs in [
]/)(πs in [
)](πs in [)]()(c o v [
NlkN
lk
NlkN
lkII
xNN
(4.2.17)
第四章 功 率 谱 估 计图 4.2.2 白噪声的周期图
0 1 2 3
5
4
3
2
1
0
0 1 2 3
5
4
3
2
1
0
P
xx
(e
j
)
P
xx
(e
j
)
5
4
3
2
1
0
0 1 2 3
0 1 2 3
5
4
3
2
1
0
P
xx
(e
j
)
P
xx
(e
j
)
/?
/?
( a ) ( b )
/ /?
( c ) ( d )
N ＝ 32 N ＝ 64
N ＝ 256
N ＝ 128
第四章 功 率 谱 估 计
4.2.3 经典谱估计方法改进
1,平均周期图法平均周期图法是基于这样的思想：对一个随机变量进行观测,得到 L组独立记录数据,用每一组数据求其均值,然后将 L
个均值加起来求平均 。 这样得到的均值,其方差将是用一组数据得到的均值的方差的 1/L。
假设随机信号 x(n)的观测数据区间为,0≤n≤M-1,共进行了
L次独立观测,得到 L组记录数据,每一组记录数据用 xi(n),i=1,
2,3,…,L表示,第 i组的周期图用下式表示：
21
0
)(1)(?
M
n
nj
ii enxMI
(4.2.18)
第四章 功 率 谱 估 计将得到的 L个周期图进行平均，作为信号 x(n)的功率谱估计，
公式如下：
L
i
ixx ILP
1
j )(1)e(
(4.2.19)
为了分析偏移，对上式求统计平均，得到
2
j
)-j(j
π
π-
j
)2/s i n (
)2/s i n (1
)]([)e(
)e()e(
π2
1
)]e(?[
M
M
mwFTW
dPWPE
BB
xxBxx
(4.2.20)
(4.2.21)
周期图的统计平均值已经求出,如（ 4.2.11）,(4.2.12)式所示，
重写如下：
第四章 功 率 谱 估 计上式表明,平均周期图仍然是有偏估计,偏移和每一段的数据个数 M有关；由于 M<N,平均周期图的偏移比周期图的偏移大,表现在三角谱窗主瓣的宽度比周期图主瓣的宽度宽 。 由于三角谱窗主瓣的宽度变宽,分辨率更加降低,因此也可以说,
偏移的大小反映分辨率的低与高 。
按照 （ 4.2.19） 式求方差,由于是 L次独立观测,L个周期图相互独立,因此平均周期图的方差为
)](v a r [1)]e(?v a r [ ijxx ILP? （ 4.2.22）
即平均周期图的估计方差是周期图的方差的 1/L。 显然,是以分辨率的降低换取了估计方差的减少,当然,估计的均方误差也减少 。
第四章 功 率 谱 估 计图 4.2.3 平均周期图法
P
xx
(e
j
)
/?
0
3
2
1
0 1 2 3
0
3
2
1
0 1 2 3
0
3
2
1
0 1 2 3
P
xx
(e
j
)
/?
P
xx
(e
j
)
/?
( a ) ( b )
( c )
N ＝ 256,L ＝ 2 N ＝ 256,L ＝ 4
N ＝ 256,L ＝ 8
第四章 功 率 谱 估 计
2.
这种方法是用一适当的功率谱窗函数 W(ejω)与周期图进行卷积,来达到使周期图平滑的目的的 。
d)e()(π2 1)e(? )(π jNlwxx WIP （ 4.2.23）
式中
m
xx
N
Nm
N mrI
j-
1
)1(
e)(?)(?
是有偏自相关函数)(? mr
xx
de)e(π2 1)( jjπ
π
nWnw?
-(M-1)≤ n≤ M-1
（ 4.2.24）
第四章 功 率 谱 估 计那么
m
xx
M
Mm
j
xx emwmrP
j
1
)1(
)()(?)e(
(4.2.25)
将 (4.2.25)式和 (4.2.3)式进行对比,它们是一样的,说明周期图的窗函数法就是前面提到的 BT法的加权协方差谱估计 。
在 （ 4.2.23） 式中,周期图和谱窗函数卷积得到功率谱,等效于在频域对周期图进行修正,使周期图通过一个线性非频变系统,滤除掉周期图中的快变成分,谱窗函数需具有低通特性 。
第四章 功 率 谱 估 计对（ 4.2.25）式求统计平均,得到
m
xx
M
Mm
xx mwmrEPE
j-
1
)1(
j e)()](?[)]e(?[?
将 (1.3.30)式和 (1.3.32)式代入上式，得到
m
Bxx
M
Mm
xx mwmwmrePE
j-
1
)1(
j e)()()()](?[?
(4.2.26)
式中
Nm
N
mmw
B ||
||1)(
第四章 功 率 谱 估 计上式表明,周期图的窗函数法仍然是有偏估计,其偏移和 wB(m)、
w(m)两个窗函数有关,如果 w(m)窗的宽度比较窄,M比 N小得多,这样 |m|<<N,则 wB(m)~1,（ 4.2.26） 式可以近似写成下式：
d)e()e(
π2
1
e)()()]e(?[
)-j(j
π
π
j-
1
)1(
j
WP
mwmrPE
xx
m
xx
M
Mm
xx
(4.2.27)
第四章 功 率 谱 估 计
3,修正的周期图求平均法这种方法和平均周期图法一样,首先把数据长度为 N的信号
x(n)分成 L段,每一段数据长度为 M,N=LM。 然后把窗函数 w(n)
加到每一个数据段上,求出每一段的周期图,形成修正的周期图,再对每一个修正的周期图进行平均 。 第 i段的修正周期图为
21
0
j-e)()(1)(?
M
n
ωn
ii nwnxUI?
i=1,2,3,…,M-1
(4.2.28)
式中
1
0
2 )(1
M
n
nw
M
U
(4.2.29)
第四章 功 率 谱 估 计同样,将每一段的修正的周期图之间近似看成互不相关,最后功率谱估计为
L
i
ixx ILP
1
j )(1)e( (4.2.30)
对上式求统计平均,得到
d)e()e(π2 1)]e(?[ )-j(j
π
j WPPE
xxxx
(4.2.31)
式中
21
0
jj e)(1)e(?
M
n
nw
MU
W (4.2.32)
第四章 功 率 谱 估 计
U称为归一化因子 。 (4.2.31)式的证明可参考文献 ［ 2］ 。 为使功率谱估计渐近无偏,（ 4.2.28） 式和 (4.2.32)式中的归一化因子是不可缺少的 。 这种在计算周期图之前,先对各数据段加窗函数的方法,使平均周期图方法的估计方差减少,当然分辨率同样减少 。 但这种方法对窗函数没有限制；此外,分段时,相邻的两段可以有重叠,进一步使方差减少,可以重叠
50%［ 6］ 。 总之,传统的功率谱估计方法无论采取哪一种改进方法,总是以减少分辨率为代价,换取估计方差的减少,提高分辨率的问题无法根本解决 。
第四章 功 率 谱 估 计对于由白噪声和正弦信号或者窄带信号组成的信号,在计算周期图之前,给数据加窗是有利的,如果不加数据窗,相近的低电平信号可能被高电平信号的旁瓣淹没掉 。 在图 4.2.4（ a）
中,ω/π =0.12处的正弦信号的旁瓣几乎掩盖了在 ω/π =0的信号,经过给数据加哈明窗处理后,功率谱如图 4.2.4(b)所示 。
由于加了哈明窗,大大压低了旁瓣,使低电平信号清晰可见,
但由于主瓣加宽,功率谱波峰变宽了,却降低了信号的分辨率 。
一般来说,两个等幅的正弦信号的频率相隔很近,可以不加数据窗,频率间隔应该大于 2π /N才能分辨 。
第四章 功 率 谱 估 计图 4.2.4
(a) 没有数据窗； (b) 加哈明数据窗
4 0,0 0
2 0,0 0
0,0 0
－ 2 0,0 0
－ 0,5 0 － 0,3 0 － 0,1 0 0,1 0 0,3 0 0,5 0
P
xx
(e
j
)
/
d
B
( a )
4 0,0 0
2 0,0 0
0,0 0
－ 2 0,0 0
－ 0,5 0 － 0,3 0 － 0,1 0 0,1 0 0,3 0 0,5 0
P
xx
(e
j
)
/
d
B
( b )
/ /?
第四章 功 率 谱 估 计
4.3 现代谱估计中的参数建模在第一章中我们介绍了用参数模型来描述随机信号的方法,
如果能确定信号 x(n)的信号模型,根据信号观测数据求出模型参数,系统函数用 H(z)表示,模型输入白噪声方差为 σ w2,信号的功率谱用下式求出：
2j2j |)e(|)e( HP wxx?
按照这种求功率谱的思路,功率谱估计可分成三个步骤：
（ 1） 选择合适的信号模型；
（ 2） 根据 x(n)有限的观测数据,或者它的有限个自相关函数估计值,估计模型的参数 ;
（ 3） 计算模型输出功率谱。
第四章 功 率 谱 估 计
4.3.1 模型选择选择模型的主要考虑是模型能够表示谱峰,谱谷和滚降的能力 。 对于具有尖峰的谱,应该选用具有极点的模型,如 AR和
ARMA模型；对于具有平坦的谱峰和深谷的信号,可以选用 MA模型；既有极点也有零点的谱应选用 ARMA模型,相对地说,ARMA
模型适用范围较宽 。 对于滚降太快的谱,没有一种模型可以准确地表示功率谱,可以选用高阶的 AR模型近似表示 。 如果选择不合适,例如,选用 MA模型去估计具有尖峰的功率谱,估计效果会很差 。 图 4.3.1表示的是用 MA模型估计二阶 AR信号功率谱的例子,图中 (a),（ b）,（ c） 中的 AR信号谱峰较平坦,用二阶 MA模型功率谱拟和真实谱时,差别较大,随着阶数的提高,
估计的谱愈来愈近似于真实的谱 。 但是对于 （ d）,（ e）,
（ f） 中 AR信号谱峰很窄,再用 MA信号模型拟和时,直到 MA模型阶数提高到 10阶,其效果仍很差 。
第四章 功 率 谱 估 计实际信号中一般都含有和信号不相关的噪声,对这一类带有噪声的信号,如果信号是 AR模型,由于有噪声需要用 ARMA模型 ; 如果用 AR模型,则需要阶数更高 ［ 5］ 。
选择模型的另一个考虑是尽量减少模型参数 。 当然这和选择模型是否合适有关系,前面介绍过三种信号模型均有普遍应用价值,但模型选择合适,功率谱估计和实际的谱拟合得好,
如果不合适,只有提高阶数才能得到较近似的谱,这样需要估计的参数增多,同样也会降低谱估计的质量 。 因此应该在选择模型合适的基础上,尽量减少模型的参数 。
第四章 功 率 谱 估 计
4.3.2 模型参数和自相关函数之间的关系根据观测数据如何来确定模型的参数对各种功率谱估计方法不尽相同 。 这里先介绍模型参数和信号自相关函数之间的关系,这些关系在功率谱估计中起着很重要的作用 。
假设模型的差分方程和系统函数分别用下式表示：
q
i
i
P
i
i inwbinxanx
01
)()()(
)(
)(
1
)(
1
0
zA
zB
za
zb
zH
P
i
i
i
q
i
i
i
(4.3.1)
(4.3.2)
第四章 功 率 谱 估 计图 4.3.1 对 AR(2)信号的模型选择
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
10
5
0
－ 5
－ 10
15
20
25
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
/?
/?
/?
/?
( a ) ( b )
( c )
( d )
dB
dB
dB
dB
真实 A R( 2 ) 的 P S D
MA (2 ) 模型真实 A R( 2 ) 的 P S D
MA (5 ) 模型真实 A R( 2 ) 的 P S D
MA (2 ) 模型真实 A R( 2 ) 的 P S D
MA (1 0 ) 模型第四章 功 率 谱 估 计
10
5
0
－ 5
－ 10
15
20
25
10
5
0
－ 5
－ 10
15
20
25
0,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,50,0－ 0,1－ 0,2－ 0,3－ 0,4－ 0,5 0,1 0,2 0,3 0,4 0,5
/ /?
( e ) ( f )
dBdB
真实 A R(2 ) 的 PS D
MA (1 0 ) 模型真实 A R(2 ) 的 PS D
MA (5 ) 模型图 4.3.1 对 AR(2)信号的模型选择第四章 功 率 谱 估 计
1,ARMA模型的系数和信号自相关函数之间的关系在第一章我们已推导出系统输出功率谱与输入功率谱之间的关系,用 （ 1.4.7） 式表示,重写如下：
*
* 1)()()(
zHzHzPzP wwxx
式中,Pxx(z)和 Pww(z)分别表示模型输出信号 x(n） 和输入信号 w(n)
功率谱的 Z变换形式 ; H（ z） 表示模型的系统函数 。 将上式两边同乘以 A(z),得到
)(
1
)(
1
)()()(
*
*2
*
*
zB
z
H
zB
z
HzPzAzP
w
wwxx
(4.3.3)
第四章 功 率 谱 估 计将上式进行 Z反变换，公式左边的反变换为
P
l
xxAAxxxx lmrlhmhmrzAzPZ
0
1 )()()()()]()([T
公式右边的反变换为
)()()(1T
0
2
*
*21 lmqlhzB
zHZ
q
l
Bww
式中
)(1T)( ***1 mhzHZmq
第四章 功 率 谱 估 计代入上式,得到
q
l
Bww lmhlhzBzHZ
0
*2
*
*21 )()()(1T
这样得到
q
l
Bwxx
p
l
A lmhlhlmrlh
0
*2
0
)()()()(?
由于模型 H(z)是因果的,h(n)=0,n＜ 0,可以得到
p
l
xxA
p
l
mq
l
BwxxA
xx
lmrlh
mlhlhlmrlh
mr
1
1 0
*2
)()(
)()()()(
)(
m=0,1,2,3,…,q
m≥ q+1
第四章 功 率 谱 估 计此式即是 ARMA模型输出自相关函数与模型参数之间的关系式 。
如果能由信号的观测数据估计出信号的自相关函数,可以按照
（ 4.3.4） 式求出 ARMA模型的参数,再由 （ 1.4.8） 式求出信号的功率谱 。 但我们看到,在公式中由于 项的存在,模型输出自相关函数和模型参数之间的关系是非线性的,
从而增加了估计功率谱的困难 。 但是当 m>q时,(4.3.4)式却是一个线性方程,用矩阵方程表示如下：
)()(
0
* mlhlh
mq
l
B
第四章 功 率 谱 估 计
)(
)2(
)1(
)(
)2(
)1(
)()2()1(
)2()()1(
)1()1()(
pqr
qr
qr
pqh
h
h
qrpqrpqr
pqrqrqr
pqrqrqr
xx
xx
xx
A
A
A
xxxxxx
xxxxxx
xxxxxx
上式共有 p个方程 。 可以用该方程首先计算出 AR部分的 p个系数
hA(i),i=1,2,3,…,p; 然后代入 (4.3.4)式,设法求出 MA部分的系数 。
第四章 功 率 谱 估 计
2.AR模型的系数和信号自相关函数之间的关系
AR模型的系统函数为,H(z)=1/A(z),相当于 ARMA模型中
B(z)=1的情况,这样在公式中
ml
ml
mlmlh
nnh
B
B
0
1
)()(
)()(
因为 h(n)是因果性的,因此 m＞ 0时 h*(-m)=0,将上面公式代入到
（ 4.3.4） 式中,得到
P
l
wxxA
P
l
xxA
xx
lrlh
mrlh
mr
1
2
1
)()(
)1()(
)(
m≥ 1
m=0
(4.3.6)
第四章 功 率 谱 估 计也可以将上式中 m≥ 1的情况写成矩阵形式：
)(
)2(
)1(
)(
)2(
)1(
)0()2()1(
)2()0()1(
)1()1()0(
pr
r
r
ph
h
h
rprpr
prrr
prrr
xx
xx
xx
A
A
A
xxxxxx
xxxxxx
xxxxxx
(4.3.7) （ 4.3.6）式的矩阵形式如下式：
0
0
1
)0()1()(
)2()0()1(
)1()1()0( 2
1
w
p
xxxxxx
xxxxxx
xxxxxx
h
a
rprpr
prrr
prrr?
(4.3.8)
第四章 功 率 谱 估 计或者用模型参数表示：
0
0
1
)0()1()(
)1()0()1(
)()1()0( 2
1
w
p
xxxxxx
xxxxxx
xxxxxx
h
a
rprpr
prrr
prrr?
(4.3.9) 令
)0()1()(
)1()0()1(
)()1()0(
)1(
xxxxxx
xxxxxx
xxxxxx
p
xx
rprpr
prrr
prrr
R
(4.3.10)
第四章 功 率 谱 估 计称为自相关矩阵,它满足 （ H表示共轭转置 ）,
是一个埃尔米特 （ Hermitian） 矩阵,且沿任一对角线的元素相等,是一个托布列斯 （ Toeplitz） 矩阵 。 也是正定矩阵 。
上面推导出的 （ 4.3.6） 式或 (4.3.9)式确定了 AR模型参数
（ 包括模型输入噪声方差 ） 和信号自相关函数之间的关系 。
我们注意到这是一个线性方程,如果能够由信号的观测数据求出其自相关函数,可以按照 （ 4.3.7） 式,通过解一组线性方程得到模型参数,相对 ARMA模型,这是 AR模型的优点 。
)1(?pxxR H
xxxx RR?
)1(?pxxR
第四章 功 率 谱 估 计
3.MA
MA模型的系统函数 H(z)=B(z),相当于 ARMA模型中 A(z)=1,
hA(n)=δ (n)的情况,此时 h(n)=hB(n),由 （ 4.3.4）
MA模型系数和信号自相关函数的关系为
0
)()(
)( 0
*2
mq
l
BBw
xx
mlhlh
mr
m=0,1,…,q
m≥ q+1 (4.3.11)
上式表明,MA模型的参数和信号自相关函数之间也是非线性关系。
第四章 功 率 谱 估 计上面我们分别推导了三种信号模型的参数和信号自相关函数之间的关系 。 这些关系式为我们提供了一种估计功率谱的方法,即首先根据信号观测数据估计信号自相关函数,再按照所选择信号模型,解上面相应的方程,求出模型参数,最后按照下式求出信号的功率谱,
2
1
022j ω2j ω
e1
e
|)e(|)e(
P
i
i
i
q
i
i
i
wwxx
a
b
HP
(4.3.12)
第四章 功 率 谱 估 计
4.4 AR模型谱估计的性质
4.4.1 AR模型的线性预测在第二章中,我们已推导出维纳线性一步预测器系数和信号自相关函数之间的关系式 (称为 Yule-Walker方程 ),重写如下：
0
0
)]([1
)0()1()(
)1()0()1(
)()1()0(
m i n
2
1
neE
a
a
rprpr
prrr
prrr
pp
p
xxxxxx
xxxxxx
xxxxxx
(4.4.1)
第四章 功 率 谱 估 计式中,e(n)表示线性一步预测误差，其公式为
p
i
pi inxanxnxnxne
1
)()()(?)()(
(4.4.2)
E［ e2(n)］ min表示 e(n)的均方差最小值； api(i=1,2,3,…,p)表示预测器的系数,它和线性预测器单位脉冲响应 h(n)之间差一符号,
即
pnanh pn,,3,2,1)(
对（ 4.4.2）式进行 Z变换,得到
p
i
i
pi zazXzXzE
1
)()()(
(4.4.3)
第四章 功 率 谱 估 计令 He(z)=E(z)/X(z)，由上式，得到
p
i
i
pie zazH
1
1)(
(4.4.4)
称 He(z)为线性一步预测误差滤波器,其作用是将信号 x(n)转换成预测误差 e(n),如图 4.4.1所示 。 一般认为 e(n)具有白噪声的性质,因此 He(z)也称为白化滤波器 。
第四章 功 率 谱 估 计图 4.4.1 预测误差滤波器
H
e
( z )
x ( n ) e ( n )
x ( n )
e ( n )
z － 1 z － 1 z － 1
a
p 1
a
p 2
a
p,p － 1
a
p,p
第四章 功 率 谱 估 计我们知道 AR模型的系统函数为
P
i
i
i za
zH
1
1
1)( (4.4.5)
对比 （ 4.4.4）,(4.4.5)两式,当 api=ai(i=1,2,3,…,p)时,
He(z)和 H(z)互为逆滤波器,He(z)=1/H(z),基于以上分析,也可以将 AR模型定义为
)()(?)( nwnxnx
(4.4.6)
式中,是基于信号前 p个样本的最佳一步线性预测，公式为)(? nx
p
i
i
pi zanx
1
)(?
第四章 功 率 谱 估 计
w(n)是模型输入白噪声 。 将 AR模型参数和信号自相关函数之间的关系式 （ 4.3.9 ） 和 (4.4.1) 式进行对比,得到,api=ai,
w(n)=e(n),E［ e2(n)］ min=σ2w,信号自相关函数和它的 AR模型参数之间的关系服从 Yule-Walker方程 。 注意,只有当 AR模型的阶数与线性预测器的阶数相同时,以上结论才是正确的 。 由于
AR模型具有这种特性,因而 AR模型法也称为线性预测 AR模型法 。
第四章 功 率 谱 估 计
AR模型与线性预测之间的关系,可以被用来解卷积,假设信号 s(n)通过一个 AR系统,系统单位脉冲响应为 h(n),响应是
x(n),即
x(n)=s(n) * h(n)
如果 s(n)具有白噪声性质,可以利用 AR模型与预测滤波器之间的关系对上式进行解卷积,得到 s(n)信号 。 方法是：先由
x(n)的观测数据估计 AR模型的参数,得到 AR模型的系统函数
H(z),再让 x(n)通过 H(z)的逆滤波器 H-1(z),便得到信号 s(n)。
第四章 功 率 谱 估 计
4.4.2 预测误差滤波器的最小相位特性我们知道,AR模型必须因果稳定,即极点均在单位圆内,
才能保证信号 x(n)是平稳随机信号,于是 AR模型 H(z)和预测误差滤波器 He(z)互为逆滤波,那么 He(z)应为最小相位系统 。 但是由解 Yule-Walker方程得到 AR模型的参数,其极点不一定在单位圆内 。 下面将证明当最佳 P阶线性预测系数与 AR模型参数相同时,由此得到的极点保证在单位圆内,AR滤波器稳定,
预测误差滤波器 He(z)或者 A(z)是最小相位系统 。 这里自相关矩阵 是正定的 。)1(?p
xxR
第四章 功 率 谱 估 计解 Yule-Walker方程得到的是最佳线性预测滤波器的系数,
此时预测误差滤波器输出功率 Pe达到最小,用 Pe min表示,即
)d(e|)(e| jππ- 2jmi n xxe PAP
(4.4.7)
式中,,ai是最佳预测系数 。 下面先用反证法证明 A(z)的全部零点不在单位圆外部 （ 即全部零点在单位圆上或者单位圆内部 ） 。 设 A(z)的第 i个零点 zi在单位圆外部,即 |zi|
＞ 1,用 1/z*i代替 zi,这时 A(z)的幅度函数不变,按照 （ 4.4.7） 式计算出的预测误差滤波器输出功率 Pe应不变,仍是最小的 。 但是下面将用公式证明预测误差滤波器输出功率 Pe不是最小的,
这一矛盾的结论只能说明 A(z)的零点不可能在单位圆外部 。
p
i
i
i zazA
1
1)(
第四章 功 率 谱 估 计
)(')1()1()1(
)1(1)(
1
111
1
1
1
zAzzzzzz
zzzazA
p
ij
j
jji
p
j
j
p
i
i
i
(4.4.8)
式中
p
ij
j
j zzzA
1
1 )1()('
将 (4.4.8)式代入（ 4.4.7）式,得到
d)e(|)e('|e1 j2j2j-mi n xxie PAzP
(4.4.9)
第四章 功 率 谱 估 计式中
2
*
2
2
*
j2
2
j22j 11||1||e1|||e1| j
i
i
i
i
i
ii ezzzezzzz
因为
1?iz
因此
2
j-
*
2j- e11|e1|
i
i zz
(4.4.10)
第四章 功 率 谱 估 计
4.4.3 AR模型隐含自相关函数延拓特性
AR模型的自相关函数和模型系数之间的关系服从 Yule-
Walker方程,重写如下：
p
l
wxxA
p
l
xxA
xx
lrlh
lmrlh
mr
1
2
1
)()(
)()(
)(
m≥ 1
m=0
上式中,对于 m≥1的情况,公式本身就是一个递推方程,如果已由观测数据计算出 p+1个自相关函数,用,m=0,1,2,,…,p
表示,对于 m＞ p的情况,可以用该公式外推得到,公式如下：
)(? mrxx
第四章 功 率 谱 估 计
P
l
xxA
xx
xx lmrlh
mr
mr
1
'
'
)(?)(
)(?
)(?
0≤ m≤ p
0>p
(4.4.11)
上式中,系数 hA(l)需用 （ 4.3.6） 式求出 。 因此 AR模型隐含着自相关函数外推的特性 。 我们知道,经典谱估计 BT法中,自相关函数只能限于由观测数据计算出的有限个自相关函数,其它的认为是 0,造成了谱估计分辨率低,模糊 。 也正是 AR模型具有自相关函数外推特性,使它具有高分辨率的优点 。
第四章 功 率 谱 估 计
4.5 AR谱估计的方法
4.5.1 自相关法 ——列文森 （ Levenson） 递推法自相关法的出发点是选择 AR模型的参数使预测误差功率最小，预测误差功率为
n
p
i
pi
n
inxanx
N
ne
N 1
22 )()(|1|)(|1?
(4.5.1a)
假设信号 x(n)的数据区在 0≤n≤N-1范围,有 p个预测系数,N个数据经过冲激响应为 api(i=0,1,2,…,p)的滤波器,输出预测误差
e(n)的长度为 N+P,因此应用下式计算,
第四章 功 率 谱 估 计
2
1
1
0
2
1
0
|)()(|1|)(|1
p
i
pi
PN
n
PN
n
inxanxNneN?
(4.5.1b)
显然,e(n)的长度长于数据的长度,上式中数据 x（ n） 的两端需补充零点,这相当于无穷长的信号经过加窗处理,得到长度为
N的数据 。 用 （ 4.5.1b） 式对系数 api的实部和虚部求微分的方法使预测误差功率最小,得到
)(?
)2(?
)1(?
)0(?)2(?)1(?
)2(?)0(?)1(?
)1(?)1(?)0(?
2
1
pr
r
r
a
a
a
rprpr
prrr
prrr
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
(4.5.2)
第四章 功 率 谱 估 计式中自相关函数采用有偏自相关估计,即
)(?
)()(
1
)(?
*
*
1
0
mr
mnxnx
Nmr
xx
mN
nxx
m=0,1,2,…,p
m=-p+1,-p+2,…,-1 (4.5.3)
对比 （ 4.3.7） 式,（ 4.5.2） 式就是已推导出的 Yule-Walker方程,因此自相关法也是基于解 Yule-Walker方程的一种方法 。
首先由信号的观测数据估计出其自相关函数,再解该方程,得到模型参数,便可求出信号的功率谱 。
Yule-Walker法 。 但是直接解该方程,需要计算逆在矩阵,不方便 。 在第三章自适应滤波器中,曾介绍了基于 Yule-Walker方程中自相关矩阵的性质,导出 Levenson-Durbin递推法,这是一种高效的解方程方法 。 Levenson[CD*2]Durbin
递推法简化重写如下：
第四章 功 率 谱 估 计
2
1
1
1
,1 /)()(?
k
k
l
xxlkxxkk lkrakra?
（ 4.5.4）
*,1,1 ikkkkikki aaaa i=1,2,3,…,k-1 （ 4.5.5）
2 122 )||1( kkkk a （ 4.5.6）
由 k=1开始递推,递推到 k=p,依次得到 {a11,σ 21},{a21,a22,
σ 22},…,{ap1,ap2,…,app,σ 2p}。 AR模型的各个系数以及模型输入白噪声方差求出后,信号功率谱用下式计算,2
1
j-
22j ω2j ω
e1
1
|)e(|)e(
P
i
i
i
wwxx
a
HP
（ 4.5.7）
第四章 功 率 谱 估 计
（ 4.5.6） 式表明：,说明随着阶数增加,
预测误差功率将减少或者不变,为此要求 |akk|≤1,akk称为反射系数 。 另外,递推公式提供了一种确定模型阶数的实验方法,如模型的阶数不知道,由低阶开始递推,当递推到 M阶时,预测误差满足允许的值,停止递推,选 AR模型的阶数为 M。 这种递推法效率高,且当阶数变化时,无需从头计算 。
如果知道信号的 N个观测数据 （ x(n),0≤n≤N-1 ）,利用列文森递推法计算功率谱的计算流程图如图 4.5.1所示 。 图中采用有偏自相关估计 （ 4.5.3） 式计算,σ2是要求的方差 。
2232221 p
)(? mrxx
第四章 功 率 谱 估 计图 4.5.1 利用列文森递推法计算功率谱的流程图输入,x ( n ),n ＝ 0,1,2,3,…,N － 1,?
2
计算：
)0(?1),0(?/)1(?
2
1,1
2
11,1 xxxxxx
rarra?
)(? mr
xx
p =2
2
1
2
2
*
,1,1
2
1
1
1
,1
1
1,,2,1,0
)(?)(?
pppp
kppppkppk
p
p
i
xxipxxpp
a
pkaaaa
iprapra
22
p
p
i
i
pi
wxx
a
P
1
j
2j
e1
1
)(e
结束
p = p + 1
N
Y
第四章 功 率 谱 估 计
4.5.2 协方差法与修正协方差法
1,协方差法这种方法和自相关法一样,仍然利用使预测误差功率最小的方法求模型参数,但由观测数据求预测误差功率的公式如下式：
1 1
2
1
2 |)()(|1|)(|1
N
pn
N
pn
p
k
pk knxanxpNnepN? （ 4.5.8）
将该式对比自相关法中求预测误差功率的公式 （ 4.5.1b）,不同的是求和限不同 。 该公式中使用的观测数据均已得到,不需要在数据两端补充零点,因此比较自相关法去掉了加窗处理的不合理假设 。 为求得模型参数仍然应用复梯度法使 （ 4.5.8） 式达到最小,公式如下：
第四章 功 率 谱 估 计
)0,(
)0,2(
)0,1(
),()2,()1,(
),2()2,2()1,2(
),1()2,1()1,1(
2
1
pc
c
c
a
a
a
ppcpcpc
pccc
pccc
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
（ 4.5.9）
式中
1
* )()(1),(
N
pn
xx knxjnxpNkjc
（ 4.5.10）
白噪声的方差为
1
m i n
2 ),0()0,0(
N
pk
xxpkxxw kcac （ 4.5.11）
第四章 功 率 谱 估 计由观测数据 x(n)(n=0,1,2,…,N-1),利用上面三个公式可以求出模型的参数,{api(i=1,2,3,…,p); σ2w}。 按照定义,
（ 4.5.9） 式中的 cxx(j,k)可以称作协方差函数,它有两个变量,
因此也适合于非平稳随机信号 。 （ 4.5.9） 式中的协方差矩阵是埃尔米特 （ Hermitian） 矩阵,cxx(k,j)=c*xx(j,k),是半正定的 。
这种方法近似于自相关法 。 一些实验结果说明它的分辨率优于自相关法 ［ 5］,另外对于纯正弦信号数据,可以有效地估计正弦信号的频率 。 进一步的内容可参考文献 ［ 14］,［ 15］
第四章 功 率 谱 估 计
2,修正协方差法修正协方差法使用前向和后向预测误差平均值最小的方法,
估计 AR模型的参数,进而估计信号的功率谱 。
Nuttall 在 1976年提出,称为前向 —后向法；同年,Ulrych
和 Ctayton也独立提出,称之为最小二乘法,因此该法也称为前后向线性预测误差功率最小的最小二乘法 。
前面我们已推导出信号的前向和后向预测,分别重写如下：
p
k
pk
p
k
pk
knxanx
knxanx
1
*
1
)()(?
)()(? (4.5.12)
(4.5.13)
第四章 功 率 谱 估 计式中 apk 是 AR模型的参数 。 最小预测误差平均功率是模型输入白噪声的方差,即 ρ p=σ 2w,前后向预测误差平均功率为
)(5.0 pbpep (4.5.14)
式中前向和后向预测误差功率 ρ pe,ρ pb分别用下式表示：
2
1
0 1
*
2
1
1
)()(
1
)()(
1
pN
n
p
k
pkpb
N
pn
p
k
pkpe
knxanx
pN
knxanx
pN
(4.5.15)
(4.5.16)
第四章 功 率 谱 估 计和协方差法一样,上式中仅对那些用到的观测数据预测误差求和 。 为了使预测误差平均功率最小,求 ρp对 apk(k=1,2,3,…,p)
的实部和虚部的微分,或者用复梯度法求,得到
pl
lnxknxanx
lnxknxanx
pNa
p
k
pk
pN
n
N
pn
p
k
pk
pl
p
,,3,2,1
0)]()()(
)()()([
1
1
**
1
0
1
*
1
(4.5.17)
第四章 功 率 谱 估 计经过简化，得到
1
*
1
0
*
1
1
*
1
0
*
)()()()(
)()()()(
N
pn
pN
n
p
k
N
pn
pN
n
pk
lnxnxlnxnx
lnxknxlnxknxa
令
1
*
1
0
* )()()()(
)(2
1),( N
pn
pN
n
xx knxjnxknxjnxpNkjc
(4.5.18)
第四章 功 率 谱 估 计将上式写成矩阵形式：
)0,(
)0,2(
)0,1(
),()2,()1,(
),2()2,2()1,2(
),1()2,1()1,1(
2
1
pc
c
c
a
a
a
ppcpcpc
pccc
pccc
xx
xx
xx
pp
p
p
xxxxxx
xxxxxx
xxxxxx
(4.5.19) 白噪声的方差估计值为
p
k
xxpkxx
pN
n
p
k
pk
N
pn
p
k
pkpw
kcac
nxknxanx
nxknxanx
pN
1
1
0 1
**
1
*
1
m i n,
2
),0()0,0(
)]()()((
)()()(([
)(2
1
(4.5.20)
第四章 功 率 谱 估 计例 4.5.1 已知信号的四个观察数据为 x(n)={x(0),x(1),x(2),
x(3)}={2,4,1,3},AR （ 1）
模型参数 。
解 （ 1） 自相关法：
5.0
09)3()41(4)24(2
0)1()(
2
1)(
)(
2
1
)1()()(|)(|
4
1
11
11111111
4
011
4
011
11
4
0
2
a
aaaa
nxne
a
ne
ne
a
nxanxnene
nn
n
第四章 功 率 谱 估 计
(2) 协方差法,
7 1 4.0
0)3()41(4)24(2
0)1()(
3
2)(
)(
3
2
)1()()(|)(|
3
1
|)(|
1
11
111111
3
111
3
111
11
3
1
2
1
2
a
aaa
nxne
a
ne
ne
a
nxanxnenene
pN
nn
n
N
pn
第四章 功 率 谱 估 计
4.5.3 伯格（ Burg）递推法设信号 x(n)观测数据区间为,0≤ n≤ N-1,前向,后向预测误差功率分别用 ρ p,e和 ρ p,b表示,预测误差平均功率用 ρ p表示,公式分别为
)(
2
1
|)(|
1
|)(|
1
,,
1
2
,
1
2
,
bpepp
N
pn
b
pbp
N
pn
f
pep
ne
pN
ne
pN
(4.5.21)
(4.5.22)
(4.5.23)
第四章 功 率 谱 估 计下面解释 （ 4.5.21）,(4.5.22)式中的求和限 。 我们已经知道,
前向,后向预测误差公式分别为
p
k
pk
b
p
p
k
pk
f
p
kpnxapnxne
knxanxne
1
*
1
)()()(
)()()( (4.5.24)
(4.5.25)
在上面两个公式中,信号项的自变量最大的是 n,最小的是
n-p。 为了保证计算范围不超出给定的数据范围,在 （ 4.5.21）
和 (4.5.22)式中,选择求和范围为,p≤n≤N-1。
第四章 功 率 谱 估 计为求预测误差平均功率 ρ p最小时的反射系数 kp，令
0?
p
p
k
(4.5.26)
将第三章中已推导出的前向、后向预测误差递推公式重写如下：
)()1()(
)1()()(
1
*
1
11
neknene
neknene
f
pp
b
p
b
p
b
pp
f
p
f
p
(4.5.27)
(4.5.28)
将上面两个公式带入（ 4.5.21）,(4.5.22)、（ 4.5.23）式中，得到
1
2
1
*
1
2
11 ]|)()(||)1()([|)(2
1 N
pn
f
pp
f
p
b
pp
f
pp nekneneknepN?
第四章 功 率 谱 估 计利用复梯度法，求 ρ p对 kp的实部和虚部的微分，并令其结果等于 0，得到
0)]()()1(
)]1()1()([
)(
1
1
*
1
*
1
1
*
111
nenekne
nenekne
pNk
f
p
f
pp
b
p
N
pn
b
p
b
pp
f
p
p
p?
解出 kp，得到
)|)1(||)((|
)1()(2
1
2
1
2
1
1
*
11
N
pn
b
p
f
p
N
pn
b
p
f
p
p
nene
nene
k (4.5.29)
第四章 功 率 谱 估 计上式就是利用伯格递推法求第 p个反射系数的公式 。 将伯格递推法求 AR模型参数的递推公式总结如下：
)|)1(||)((|
)1()(2
1,,3,2,1,0)()(
1,,3,2,1,0)()(
)0(?
|)(|
1
)0(?
2
1
1
2
1
1
*
11
0
0
0
1
0
2
N
pn
b
p
f
p
N
pn
b
p
f
p
p
b
f
xx
N
n
xx
nene
nene
k
Nnnxne
Nnnxne
r
nx
N
r
(4.5.30)
(4.5.31)
(4.5.32)
(4.5.33)
(4.5.34)
第四章 功 率 谱 估 计
2,,1,)()1()(
1,,2,1)()()(
1,,3,2,1
)||1(
1
*
11
11
,
*
,1,1
1
2
Nppnneknene
Nppnneknene
ka
piakaa
k
f
pp
b
p
b
p
b
pp
f
p
f
p
ppp
ipppippi
ppp
(4.5.36)
(4.5.35)
(4.5.37)
(4.5.38)
(4.5.39)
对于 p阶 AR模型的输入白噪声方差 σ 2w=ρ p。 利用伯格递推法求
AR模型参数的流程图如图 4.5.2所示 。
第四章 功 率 谱 估 计图 4.5.2 伯格递推法流程图输入,x ( n ),n ＝ 0,1,2,3,…,N － 1
阶数 IP
1)(
1
)()()()(
2
1
0
0
00
pnx
N
nxnenxne
N
n
bf
p = IP
结束
1 pp
1
2
1
2
1
*
1
1
1
)()(
)1()(2
N
pn
b
p
f
p
b
p
N
pn
f
p
p
nene
nene
k
pwppp
k
2
1
2
1
pppipppippi
kaakaa
*
,1,1
1,,2,1),()1()(
1,,2,1),1()()(
1
*
1
11
Nppnneknene
Nppnneknene
f
pp
b
p
b
p
b
ppp
f
p
2
,;,,3,2,1:
wip
pia输出
2
1
j
2
j
e1
)(e
p
i
i
pi
w
xx
a
P
N
Y
第四章 功 率 谱 估 计
4.5.4 关于 AR模型阶次的选择在 AR模型谱估计中,模型阶次的选择也是一个关键问题 。
一般模型的最好选择是先验未知的,实际中需预先选定模型阶次 。 如果是纯 P阶 AR信号,选择模型阶次 k＜ P时,将产生对谱的平滑作用,降低谱的分辨率,如图 4.5.3所示 。 图中,AR信号 P=4,选择模型的阶次 k=2,产生的平滑作用使两个峰变成一个峰,分辨率明显降低 。 如果选择 k≥ P,且假定观测的数据没有误差 （ 没有干扰 ）,估计的参数应是：
kppi
pia
a ipki
,,2,10
,,3,2,1,
第四章 功 率 谱 估 计图 4.5.3 AR模型阶次太小时的平滑作用
5 0,0 0
3 0,0 0
1 0,0 0
－ 1 0,0 0
－ 3 0,0 0
0,0 0 0,1 0 0,2 0 0,3 0 0,4 0 0,5 0
真实 A R(4 ) 的 PS D
A R(2 ) 的 PS D 模型
P
S
D
/
d
B
/?
第四章 功 率 谱 估 计因此,对于纯 P阶 AR信号,应选择阶次 k≥P。 一般对于纯 AR信号,
只要记录数据的长度不十分短,模型阶次的估计结果比较满意
［ 21］ 。 但是如果是白噪声中的 AR信号 （ 观测数据有误差或者信号中含有白噪声 ）,此时选择 ARMA模型合适,如选择了 AR模型,其阶次应加大,较低的阶次会使谱估计产生偏移,降低分辨率 。 当然,这也和信噪比有关,信噪比愈低,平滑作用愈严重,愈需要高的阶次,因此信噪比低应选高的阶次 。 一般来说,
阶次愈高,分辨率愈高；但阶次太高,会使估计误差加大,谱峰分裂,因此,对于白噪中的 AR信号,其阶次的选择应折衷考虑 ［ 5］ 。 下面介绍三种模型阶次的估计方法 。
第四章 功 率 谱 估 计为了确定模型的阶次,所提出的各种方法都是基于估计出的预测误差功率,即选择模型的阶次,使预测误差功率极小化,
这是因为所讨论的 AR模型参数估计方法,预测误差功率都是随模型的阶次的增加而减少,或者不变的 。 但是不能仅仅考虑预测误差功率的减少,也要考虑因模型阶次增高引起模型参数加多,引起估计误差的加大 。 根据这一原理,阿凯克提出两种估计方法,第一种方法称为最终预测误差 （ FPE） 准则,用下式计算后再取最小,估计模型阶次：
kkN
kNkFPE)(
（ 4.5.40）
第四章 功 率 谱 估 计式中,表示 k阶 AR模型的白噪声方差 （ 预测误差功率 ） 的估计值 ;N是观测数据的长度 ; k是模型阶次 。 上式中,随 k增加而减少,但 (N+k)/(N-k)却随 k增加而加大 。 第二种方法应用更广泛,称为阿凯克信息论准则：
kp?
kp?
knk k 2?ln)(A I C （ 4.5.41）
选择使上式最小的 k值作为模型的阶次 。 但这种方法有过高估计模型阶次的趋势,例如真实的阶次是 2,而 AIC(k)最小值对应的是 4,为改进,提出用 k lnN代替公式中的 2k的方法 ［ 22］,［ 23］ 。
AIC和 FPE估计阶次的性能是相近的,但对于短记录数据,
建议使用 AIC准则 ［ 24］,对于长记录 （ N→∞ ）,两种方法得到相同的模型阶次 。
第四章 功 率 谱 估 计第三种方法由 Parzen在 1976年提出,称为自回归传递函数准则 （ CAT） ［ 25］,定义是：
k
k
i iN
k1~11)(C A T
1
(4.5.42)
式中
ii iN
N~
同样,模型的阶次取使上式最小的 k值 。 以上三种估计模型阶次的方法性能类似 ［ 26］ 。 还有其它的模型阶次估计方法,
可参考文献 ［ 27］ ～ ［ 29］ 。
第四章 功 率 谱 估 计
4.6 最大熵谱估计与最大似然谱估计
4.6.1 最大熵谱估计
1,利用最大熵的原则外推自相关函数按照 Shannon对熵的定义,当随机变量 X取离散值时，熵的定义为
i
ii ppH ln
(4.6.1)
式中 pi是出现状态 i的概率。当 X取连续值时，熵的定义为
)]([ l nd)(ln)( xpExxpxpH
(4.6.2)
第四章 功 率 谱 估 计式中,p(x)是 X的概率密度函数,对于离散随机序列,概率密度函数用联合概率密度函数代替 。 显然,熵代表一种不确定性,最大熵代表最大的不确定性,或者说最大的随机性 。 下面我们研究对于有限的自相关函数值不作任何改变,对于未知自相关函数用最大熵原则外推,即不作任何附加条件的外推方法 。 假设 x(n)
是零均值正态分布的平稳随机序列,它的 N维高斯概率密度函数为
XNRXNRxxxp
xx
H
xx
N
N
12/12/
21 ))((2
1e x p))(( d e t)π2(),,,(?
式中
H21 ],,,[ NxxxX
第四章 功 率 谱 估 计
)0()1()(
)1()0()1(
)()1()0(
)(
xxxxxx
xxxxxx
xxxxxx
xx
rNrNr
Nrrr
Nrrr
NR
按照（ 4.6.2）式,x(n)信号的熵为
])))(( d e t ()πe2l n [ ( 2/12/ NRH xxN? （ 4.6.3）
式中 det(Rxx(N))表示矩阵 Rxx(N)的行列式,由上式表明为使熵最大,要求 det(Rxx(N)最大 。
第四章 功 率 谱 估 计若已知 N+1个自相关函数值 rxx(0),rxx(1),…,rxx(N),下面用最大熵方法外推 rxx(N+1)。 设 rxx(N+1)确实是信号自相关函数的第
N+2个值,根据自相关函数的性质,由 N+2个自相关函数组成的矩阵为
)0()1()()1(
)()1()0()1(
)1()()1()0(
)1(
xxxxxxxx
xxxxxxxx
xxxxxxxx
xx
rrNrNr
NrNrrr
NrNrrr
NR
（ 4.6.4）
第四章 功 率 谱 估 计它必须是非负定的矩阵,即
0)]1(d et [NR xx （ 4.6.5）
0)](d e t [,),1(d e t [),0(d e t [?NRRR xxxxxx?
将行列式展开,det(Rxx(N+1))是 rxx(N+1)的二次函数,该二次函数系数的符号是,(-1)1+N+2(-1)1+N+1=-1,且 det(Rxx(N+1)) 对
rxx(N+1)的二次导数是 -2det［ Rxx(N-1)］,它是负值,负值表示
det (Rxx(N+1))对 rxx(N+1)的一次导数是减函数,det(Rxx·(N+1))作为 rxx(N+1)的函数,凹口向下,那么只有一个最大值 。 为选择
rxx(N+1)使 det(Rxx(N+1)最大,解下列方程：
0)]1(d e t [)1(d d NRNr xx
xx
(4.6.6)
第四章 功 率 谱 估 计用数学归纳法，得到
0
)1()()1(
)2()1()2(
)1()0()1(
xxxxxx
xxxxxx
xxxxxx
rNrNr
Nrrr
Nrrr
（ 4.6.7）
上式是 rxx(N+1) 的一次函数,可以解出 rxx(N+1) 。 继续再将
rxx(N+1) 代入 Rxx(N+2) 和 det(Rxx(N+2)) 中,求 det(Rxx(N+2)) 对
rxx(N+2)的最大值,得到 rxx(N+2); 以此类推,可推出任意多个其它自相关函数值,而不必假设它们为零,这就是最大熵谱估计的基本思想 。
第四章 功 率 谱 估 计
2,最大熵谱估计与 AR模型谱估计的等价性我们已经知道 AR
Yule-Walker方程,即
N
l
wxxA
N
l
xxA
xx
lrlh
lmrlh
mr
1
2
1
)()(
)()(
)(
m≥ 1
m=0
将 m≥ 1的情况写成矩阵形式：
)(
)2(
)1(
)0()2()1(
)2()0()1(
)1()1()0(
2
1
Nr
r
r
a
a
a
rNrNr
Nrrr
Nrrr
xx
xx
xx
N
xxxxxx
xxxxxx
xxxxxx
第四章 功 率 谱 估 计式中,ai=hA(i),i=1,2,3,…,N,ai是 AR模型系数 。 解该方程,可以得到模型系数 ai,即
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
)0()2()(
)2()0()2(
)1()1()1(
1
（ 4.6.8）
第四章 功 率 谱 估 计
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
)0()1()(
)2()1()2(
)1()0()1(
2
D
rNrNr
Nrrr
Nrrr
a
xxxxxx
xxxxxx
xxxxxx
N
N
)1()1()(
)3()1()2(
)2()0()1(
)1(
（ 4.6.9）
（ 4.6.10）
第四章 功 率 谱 估 计
)0()()1(
)2()2()1(
)1()1()0(
xxxxxx
xxxxxx
xxxxxx
rNrNr
Nrrr
Nrrr
D
(4.6.11)
在（ 4.3.6）式中，令 m=N+1，得到
N
l
xxlxx lNraNr
1
0)1()1(
(4.6.12)
将以上求出的系数 a1,a2,…,aN代入上式,求出 rxx(N+1)。 而最大熵外推自相关函数的公式是 (4.6.7)式,按照该公式的最后一行展开,得到第四章 功 率 谱 估 计
0)1()1(
0)1()1()()1(
0)1()1()()1(
1
21
21
kNraNr
araNraNrNr
DarDaNrDaNrDNr
xx
N
k
kxx
Nxxxxxxxx
Nxxxxxxxx
（ 4.6.13）
上式即是最大熵外推自相关函数的公式,对比 （ 4.6.12） 式,
两公式完全一样,证明了 AR模型功率谱估计和最大熵谱估计的等价性 。 这里最大熵外推自相关函数等价于已知 N+1个自相关函数,匹配一个 N阶 AR信号模型的系数 。 一旦通过解 Yule-Walker
方程,解出模型参数,最大熵谱估计用下式计算信号功率谱：
N
k
kj
k
wj
xx
a
P
1
2
e1
)e(
（ 4.6.14）
第四章 功 率 谱 估 计
4.6.2 最大似然谱估计 ——最小方差谱估计最大似然谱估计是用一个 FIR滤波器实现,该滤波器对所关心频率的正弦信号,可以无失真地通过,而对于其它频率的信号,让其频响尽可能地小,亦即将它们尽可能地滤除 。 此时,
滤波器输出的均方值,就作为信号的功率谱估计 。
设实信号用 x(n)表示,FIR滤波器系统函数用 A(z)表示：
p
n
nznazA
0
)()(
p
k
pAXknxkananxny
0
T)()()()()(
输出 y(n)为
(4.6.15)
第四章 功 率 谱 估 计式中
T
T
)](,),1(),([
)](,),1(),0([
pnxnxnxX
paaaA p
输出信号的均方值为
ppHppTHppTHpy ARAAXXEAAXXAEnyEP ][][]|)([|? 2
(4.6.16)
上式中 T表示转置,H表示共轭转置,Rp=E［ XXT］ 是 Toeplith
自相关矩阵,为求,必须先求 FIR滤波器的系数 。 求这些系数的原则是：在所关心频率 ωi处,信号 x(n)无失真地通过,
即在 ω i处的传输函数为 1：
yP?
第四章 功 率 谱 估 计
pi
H
p
n
p
n
AEA
naA
i
ii
)()(e
1e)()e(
j
j-
0
j
式中
Hjj2j ]eee1[ pωωωH
p iiiE
(4.6.17)
另外一个原则是在 ωi附近的频率分量尽量衰减掉,即 ω≠ωi处,滤波器输出 y(n)的均方差 最小,即 (4.6.16)式最小,此时 作为信号 x(n)的功率谱估计 。 因此,最大似然谱估计称为最小方差谱估计更为合适,但由于习惯也可以仍称为最大似然谱估计 。 在以上原则下,使方差 最小的滤波器系数和 分别为
［ 30］,［ 31］
yP? yP?
xP?
yP? x
P?
第四章 功 率 谱 估 计
)()(
1
)()(
)(
1
1
1
pi
H
p
yx
ipi
H
p
ipp
p
RE
PP
RE
ER
A
应该指出,此时 并不是真正意义上的信号功率谱,只是描述了信号功率谱的相对强度 。
yP?
第四章 功 率 谱 估 计下面分析最小方差谱估计与 AR模型谱估计之间的关系：
1)2()1()(
01)1()2(
001)1(
0001
1
1
papapa
aa
a
A
ppp
pp
p
pp
（ 4.6.20）
第四章 功 率 谱 估 计
2
2
1
2
0
0
0
p
w
P
(4.6.21)
第四章 功 率 谱 估 计
4.7 特征分解法谱估计
4.7.1 正弦波用退化 AR模型表示无论是实正弦波还是复正弦波,都可以用一个退化 AR模型表示,设 P个实正弦波组成的信号用下式表示：
)(s in)(
1
ii
P
i
i nqnx
（ 4.7.1）
式中,初相位 θ i是在区间 （ -π,π ） 均匀分布的随机变量,
首先分析下面的三角恒等式：
])1(s i n [c o s2])2(s i n [)s i n ( nnn
-π ＜ ω ＜ π
第四章 功 率 谱 估 计令 x(n)=sin(ω n+θ ),则上式变为
)2()1(c o s2)( nxnxnx? （ 4.7.2）
将上式进行 Z变换，得到
0)co s21)(( 21 zzzX? （ 4.7.3）
这样 (4.7.2)式的特征多项式为
0c o s21 21 zz? （ 4.7.4）
第四章 功 率 谱 估 计上式的两个根分别是,z1=ejω,z2=e-jω,它们共轭成对,且模为 1。
由这两个根可以确定正弦波的频率 。 对比 AR模型的系统函数,
可以把正弦波信号用一个特殊的 AR（ 2） 模型表示,括弧中的 2
表示模型是二阶的 。 该 AR模型的激励白噪声方差趋于 0,极点趋于单位圆 。 通常称为退化的 AR模型 。 这一模型系数有两个,即
2 cosω 和 1,（ 4.7.2） 式是模型的差分方程 。
第四章 功 率 谱 估 计对于 P个实正弦波,特征多项式是
0)c o s21(
1
21
p
i
i zz?
（ 4.7.5）
上式是 z-1的 2P阶多项式，可以表示为
10 0
2
0
aza
P
k
k
k
（ 4.7.6）
注意上式中的系数 ak(k=1,2,3,…,2P),必须保证它的根共轭成对 。
考虑到根共轭成对,也可表示为
P
i
ii
P
k
k
k zzzzza
1
*
2
0
))((
（ 4.7.7）
第四章 功 率 谱 估 计这样由 (4.7.6)式,P个正弦波组合的模型用下面 2P阶差分方程描述
P
k
k knxanx
2
1
)()(
（ 4.7.8）
对于复正弦波情况,P个复正弦波组成的信号是
P
i
n
i
iiqnx
1
)j(e)(
（ 4.7.9）
用一个退化的 AR（ p）模型表示的差分方程为
p
i
k n - kxanx
1
)()(
（ 4.7.10）
第四章 功 率 谱 估 计其特征多项式为
10 0
0
aza
P
k
k
k
(4.7.11)
其根为
ijiz?e?
1≤ i≤ P
注意这里的根不是共轭成对出现的 。
总结以上 P个正弦波组合是一个退化的 AR(2P)过程,独立参量个数为 P个； P个复正弦波的组合是退化的 AR（ P） 过程,独立参量个数仍为 P个 。 实正弦过程相应的退化 AR过程的阶数比复正弦情况的阶数高 1倍 。
第四章 功 率 谱 估 计
4.7.2 白噪声中正弦波组合用一特殊的 ARMA模型表示白噪声中正弦波组合的信号为
)()s in ()()()(
1
nwnqnwnxny
P
i
iii
（ 4.7.12）
式中,w(n)为白噪声,且
0)]()([,)]()([,0)]([,2 lwnxElwnwEnwE lnw
将 (4.7.12)式中 x(n)的用 AR(2P)表示，即将 (4.7.8)式带入 (4.7.12)
式中,得到
)()()(
2
1
nwinxany
P
i
i
（ 4.7.13）
第四章 功 率 谱 估 计将（ 4.7.12）式中的 n用 n-i代替,
x(n-i)=y(n-i)-w(n-i)
再将上式带入 (4.7.13)式,得到
P
i
i
P
i
i inwanwinyany
2
1
2
1
)()()()(
(4.7.14)
上式可以看成一个特殊的 ARMA(2P,2P)模型,括弧中的两个 2P分别表示 ARMA模型系统函数分子和分母的阶次 。 它与一般的
ARMA模型比较,有三方面不同：
第四章 功 率 谱 估 计
(1) 它的 AR部分和 MA部分具有相同的参数,它们存在共同的因子；
(2) 由于特征多项式 (4.7.6)式的根的模为 1,故 AR部分特征多项式不满足平稳性条件,MA部分特征多项式也不满足可逆性条件；
（ 3） AR部分的 y(n)=x(n)+w(n),y(n)是含白噪声的观测值,
而通常为信号的 x(n)不含白噪声 。
第四章 功 率 谱 估 计
4.7.3 特征分解法谱估计这种特殊的 ARMA模型结构,不能用一般的 ARMA模型结构求解 。 下面介绍特征分解技术 。
将（ 4.7.14）式写成矩阵形式：
YTA=WTA (4.7.15)
式中
T
T
221
T
)]2(,),1(),([
],,,,1[
)]2(,),1(),([
pnwnwnwW
aaaA
pnynynyY
p
第四章 功 率 谱 估 计用向量 Y左乘 (4.7.15)式,并取数学期望,得到
E［ YY T］ A=E［ YW T］ A (4.7.16)
式中
yy
yyyyyy
yyyyyy
yyyyyy
R
rprpr
prrr
prrr
YYE?
)0()12()2(
)12()0()1(
)2()1()0(
][
T
IWWEWWXEYWE w2TTT ][])[(][
第四章 功 率 谱 估 计将上面关系式带入 (4.7.16)式,得到
AAR wyy 2
(4.7.17)
式中,σ w2是自相关函数 Ryy的特征值 ; A是对应 σ w2的特征矢量 。
由于 x(n)与 w(n)互不相关,由 （ 4.7.12） 式求 y(n)的自相关函数,
得到
IRR wxxyy 2
式中 Rxx是 x(n)的自相关函数,可以推导出在特征方程 (4.7.17)
式中,σ w2是 Ryy的最小特征值,且 Ryy的阶数超过 (2P+1)(2P+1)时,
σ w2就是其多重最小特征值 。 这一结论为我们寻求向量 A提供了重要的依据 。
(4.7.18)
第四章 功 率 谱 估 计当特征向量 A求出后,就可以通过解特征方程,解出各个根,求出各正弦波的频率值,特征方程为
01 222211
2
0
P
P
P
k
k
k zazazaza?
（ 4.7.19）
该方程有 2P个根，这些根在 z平面的单位圆上，它们是
iωiz je?
i=1,2,3,…,2P （ 4.7.20）
式中的 ω i即是正弦波的频率 。 当各正弦频率由 (4.7.19)式求出后,各正弦波的功率也可以求出,对于白噪声中 P个实正弦组合信号,它的自相关函数为
P
i
iwyy Pr
1
2)0(? （ 4.7.21）
第四章 功 率 谱 估 计
mPmr i
P
i
iyy?c o s)(
1
m≠ 0 (4.7.22)
式中
2
2
i
i
qP?
qi是频率为 ωi的正弦波的幅度,Pi是其功率。
将 (4.7.22)式写成矩阵形式,
FP=r (4.7.23)
式中
P
P
P
PPP
F
c o sc o sc o s
2c o s2c o s2c o s
c o sc o sc o s
21
21
21
(4.7.24)
第四章 功 率 谱 估 计
PT=［ P1,P2,…,PP］
)](,),2(),1([T Prrrr yyyyyy
(4.7.25)
(4.7.26)
由于各正弦波的频率已求出,矩阵 F,r已知,由 (4.7.23)式解出各正弦波的功率或幅度 。 最后噪声功率由下式求出,
p
i
iyyw Pr
1
2 )0(?
(4.7.27)
以上就是皮萨论科谱分解法的全过程 。 其中关键的一步是求自相关函数 Ryy的最小特征值及对应的特征向量 。 求最小特征值对应的特征向量经典方法是幂法 ［ 32］,可以利用 Ryy的布利斯结构,
采用升幂法迭代公式求解 A,记第 k次迭代值为 A(k),迭代公式为
)()1( 1 kARkA yy
(4.7.28)
第四章 功 率 谱 估 计上面公式将收敛到 Ryy的最小特征值对应的特征矢量 A,A的理想初值是
T]1,,1,1[)0(A
(4.7.29)
通常只需迭代几次就能求得所需的稳定特征矢量 A（ ∞ ）,一旦求得 A,由 （ 4.7.17） 式求出 Ryy的最小特征值 σ 2w:
AA
ARA
T
yy
T
w?
2?
(4.7.30)
第四章 功 率 谱 估 计此外,皮萨论科方法需要精确知道自相关函数矩阵,但在实际中只能由 y(n)的观测数据进行估计,估计误差会影响估计频率的精确性,为此,玛坡尔 (Marple)提出了改进方法,这方面内容可参考文献 ［ 32］ 。 在处理实际问题时,通常不知道正弦波的个数,即模型的阶数 2p是不知道的 。 如果 Ryy的维数是
(2p+1)(2p+1),最小特征值是 σ 2w,假如增大 Ryy的维数,它的最小特征值重复数目增加,但都是 σ 2w 。 这样可以逐步增加阶数,直到估计出相邻的最小特征值相等时,就认为得到了模型的正确阶数 。 为了清楚地了解皮萨论科方法的计算过程,
具体计算步骤归纳如下：
第四章 功 率 谱 估 计
(1) 由 y(n)的观测数据估计出 2p+1个相关函数 ryy(m),m=0,
1,2,…,2p值,给定阶数 p初始值 。
(2) 求 Ryy的最小特征值和特征向量 。 可以采用 （ 4.7.28）
~(4.7.30)式进行迭代计算 。
(3) 进行最小特征值判断 。 如此次求得的最小特征值和上一个 p值时求得的相同,则认为上次的阶数 p即为所需的阶数,否则继续增加阶数进行计算 。
(4) 由 （ 4.7.19） 式求特征多项式的根 。
(5) 由求得的特征多项式的根确定各正弦波的频率 。
(6) 由方程（ 4.7.23）和 (4.7.27)式,求各正弦波的功率和噪声功率。
