第四章 高阶谱估计的常规方法

{ } { } var
Re
⎡⎣
Mˆ
x n
(λ1
,"
,
λn
−1
)
⎤⎦
≅ var
Im
⎡⎣
Mˆ
x n
(λ1
,"
,
λn−1
)
⎤⎦
≈
N
⋅
M
x 2
(λ1
)"
M
x 2
(λn−1) ⋅
M
x 2
(λ1
+" +
λn −1 )
(4.18)
其中
M
x 2
(λ)
是{x(k
)}
的真实功率谱，
M
x n
(λ1,",
λn−1
)
是
{x(k
步骤 5. 对于随机信号可用累积量与矩的关系式（(2.5)式）生成累积量 cˆnx (τ1,",τ n−1) 。
如果每段去均值，则有
cˆ2x (τ 1 ) = mˆ 2x (τ 1 )
cˆ3x (τ1,τ 2 ) = mˆ 3x (τ1,τ 2 )
授课教师：姬红兵教授 hbji@xidian.edu.cn
无偏估计： lim[E[rˆ] − r] = 0 ， 即估计的数学期望等于其真值； N →∞
一致估计： lim E[(r − rˆ)2 ] = 0 ，即估计的方差随数据增长而趋于零的估计； N →∞ 实际中，遇到的情况是利用有限长度数据估计一过程的高阶累积量谱。一般有两
种主要方法： (1) 常规（Fourier 型）方法 (2) 参数方法：AR、MA、ARMA 或 Volterra 模型
Δ
其中， w(τ1,",τ n−1) 是有界支撑域的连续窗函数，且 Δn 是带宽，通常取 Δ n =1/ Ln 。注
意到，对于 K = 1，可选择 Ln 和 M ，以使当 Ln → ∞, M → ∞ 时，有 (L2n / M ) → 0 。
此外，利用高阶谱的对称性可减少（4.4）和（4.5）式的计算复杂度。
2．在估计的高阶统计量的支撑区外为零，即
w(τ1,",τ n−1) = 0 ， τi > Ln ,i = 1, 2,", n − 1
授课教师：姬红兵教授 hbji@xidian.edu.cn
59
更新日期 2011 年 4 月 1 日
研究生课程：现代信号处理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
)}
的真实
n
阶矩谱。
显然，高阶周期图估计是渐近无偏的，而其方差正比于功率谱的乘积。进而，当 N 增大时，估计方差也增大。因此，高阶周期图估计是非一致估计（inconsistent）。
通常有两种减小估计方差的途径： (1) 在邻域频率上平滑周期图；（频域平滑） (2) 平均不相邻时间段上的周期图估计。（Welch 法，时域平滑）
M
x n
(λ1
,"
,
λn
−1
)
=
1 N
⋅ Fx (λ1)" Fx (λn−1) ⋅ Fx∗ (λ1
+" + λn−1)
Brilliges 和 Rosenblatt 指出，对于大的 N，有
{ } E
Mˆ
x n
(λ1,"
,
λn−1
)
≈
M
x n
(λ1,",
λn −1 )
(4.16) (4.17)
并有渐近方差
4.3.3 窗函数
与常规功率谱估计类似，为了获得平滑的估计，必须选择合适的窗。用于高阶谱 估计的窗函数 w(τ1,",τ n−1) 应当满足下列性质：
1．具有高阶矩或累积量的对称性，例如，在实信号双谱估计中，窗函数应满足：
w(τ1,τ 2 ) = w(τ 2 ,τ1) = w(−τ1,τ 2 −τ1) = w(τ1 −τ 2 , −τ 2 )
图 4.1（见教材）示出一些常用的一维窗函数。图 4.2（见教材）示出（4.9）式最
优 MSE 窗，以及用（4.6）式产生的二维窗函数。这些窗函数可以用双谱偏差上确界（J）、
双谱方差（V）、以及估计双谱与真实值之间的均方误差（MSE）等来评价。Rao 和 Gabr
[1984]证明了 MSE 直接正比于一个效率指数 E，定义为 E = V ⋅ B ，其中
授课教师：姬红兵教授 hbji@xidian.edu.cn
59
更新日期 2011 年 4 月 1 日
研究生课程：现代信号处理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
4.2 间接法
该方法是先用有限长数据估计高阶统计量，然后用多维窗函数产生高阶谱。因此， 类似于用 Blackman-Tuckey 法估计功率谱。
研究生课程：现代信号处理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
第四章 高阶谱估计的常规方法
4.1 引言
功率谱估计回顾
谱估计的现代史是从 Tuckey 于 1949 年的突破开始的。Blackman-Tuckey 给出了用 Wiener 相关法从采样数据序列得到功率谱估计的实现方法——BT 法。
其中， n = 2,3,"， i = 1,2," K ，τ k = 0,±1,±2,"， s1 = max(0,−τ1,",−τ n−1 ) ，
s2 = min(M −1, M −1 −τ1,", M −1 −τ n−1 ) ， τ k ≤ Ln 。 Ln 为确定估计的 n 阶矩函数的支
撑区。 步骤 4. 平均所有段，即
∑ mˆ
x n
(τ
1
,"
,τ
n
−1
)
=
1 K
K
m (i ) n
(τ 1 ,",τ
n−1 )
，
n
=
2,3,",
τ
k
i =1
≤ Ln
（4.2）
mˆ
x n
(τ
1
,"
,τ
n−1
)
是
m
x n
(τ
1
,"
,τ
n
−1
)
的
一
般
估计。如
果
信号是确
定
性的，则
mˆ nx (τ 1,",τ n−1 ) = m1n (τ 1,",τ n−1 ) ，即 K=1，因此 N=M。
其中， τk ≤ Ln, k = 1, 2,3
Rosenblatt 和 Van Ness[1965 年]指出这种估计有两个基本要求： (1) 估计应是渐近无偏的； (2) 当 K → ∞ 或 K =1， N → ∞ 时，估计的方差应趋于零，即一致估计。
（4.3）
4.2.2 高阶谱估计
n 阶矩谱估计：
累积量谱估计：
Ln
Ln
∑ ∑ Mˆ
x n
(ω1 ," , ωn −1
)
=
"
mˆ nx (τ1,",τ n−1)
τ1 =− Ln τ n−1 =− Ln
⋅ w(τ1,",τ n−1) ⋅ exp{− j(ω1τ1 +" + ωn−1τ n−1)}
（4.4）
∑ ∑ Cˆnx (ω1,",ωn−1) = Ln " Ln cˆnx (τ1,",τ n−1) τ1 =− Ln τ n−1 =− Ln ⋅ w(τ1,",τ n−1) ⋅ exp{− j(ω1τ1 +" + ωn−1τ n−1)} （4.5）
{ } 步骤 3. 假设 x(i) (k), k = 0,1,", M −1 是第 i 段数据 (i = 1,", K ) ，则其高阶矩估计为
∑ mn(i) (τ 1 ,",τ n−1 )
=
1 M
s2
x (i) (k )x (i) (k
k =s1
+ τ1)" x(i) (k
+ τ n−1 )
（4.1）
授课教师：姬红兵教授 hbji@xidian.edu.cn
59
更新日期 2011 年 4 月 1 日
研究生课程：现代信号Hale Waihona Puke Baidu理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
4.3 直接法
由（2.27）或（3.47）多谱的定义提供了一种估计多谱的方法。如果有有限长观测 数据，估计它的离散傅立叶变换系数，作三重（四重）乘积并平均即得双谱（三谱 ） 的估计。类似于平均周期图法或 Welch 法。 4.3.1 高阶周期图（Higher-Order Periodogram）
59
更新日期 2011 年 4 月 1 日
研究生课程：现代信号处理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
cˆ4x (τ1,τ 2 ,τ 3 ) = mˆ 4x (τ1,τ 2 ,τ 3 ) − mˆ 2x (τ1 )mˆ 2x (τ 3 − τ 2 ) − mˆ 2x (τ 2 )mˆ 2x (τ 3 − τ1 ) − mˆ 2x (τ 3 )mˆ 3x (τ 2 − τ1 )
设{x(k)} ， k = 0,1,", N − 1，是一零均值实平稳时间序列，其 DFT 为
∑ Fx (λ)
=
N −1 x(k) exp ⎧⎨−
k =0
⎩
j
2π N
λk
⎫ ⎬
⎭
， λ = 0,1,", N −1
(4.15)
如果{x(k)} 没有去均值，则可置 Fx (0) = 0 。
高阶周期图定义为
4.2.1 高阶统计量估计
设 {x(1), x(2),"x(N)}是给定的数据集，它可表示一个严格平稳随机过程或确定性
序列的实现，则有以下估计方法： 步骤 1. 将数据分 K 段，每段 M 个样点，即 N=K·M
然而，如果数据样本对应一个确定性能量信号，则数据分段是不合适的。同样， 如果过程是确定性周期的，则 M 应等于信号的周期或周期的整数倍； 步骤 2. 每段数据去均值；
−π
表 4.2（见教材）列出了 MSE 最优窗和图 4.2 所示的乘积窗函数对应的 J，V，B 值。可见，MSE 最优窗具有最小的 E 值。Sasaki 窗具有最小的 J 值，这是因为它最初 是由使双谱偏差上确界最小化而导出的。另一方面，Parzen 窗具有最小的 V 值。 例 4.1（见教材） 例 4.2（见教材）
3．在原点等于 1，即 w(0,",0) = 1 (归一化条件)
4．具有实非负傅立叶变换，即W (ω1,",ωn−1 ) ≥ 0 ， ωi ≤ π ，i = 1,2,", n − 1 ，窗函数 也应具有有限能量。 利用标准一维滞后窗函数可以容易地生成一类满足这些性质的窗函数，即有
w(τ1,",τ n−1) = d (τ1)d (τ 2 )"d (τ n−1)d (τ1 +τ 2 +" +τ n−1)
本章讨论常规方法以及他们的统计特性和计算复杂性： 常规方法可分成下列三类： (1) 间接法（Indirect）：是定义（2.24）式或(3.45)式的近似； (2) 直接法(Direct)：是定义（2.27）式或(3.47)式的近似 (3) 复解调制法。（Complex demodulation） 虽然常规法直接且其实现可利用 FFT，但估计的统计方差和频率分辨率的限制对 其应用能力有严格的限制。
1965 年 FFT 的出现产生了周期图法； BT 法与周期图法的致命弱点是频率分辨率的限制。为了克服这一缺点，1967 年 Burg 受到他本人在地震应用研究中线性预测方法的启发，导出了最大熵谱估计法；E. Parzen 于 1968 年正式提出了 AR 谱估计方法；此后十几年来发展了许多高分辨的谱估 计方法，称为现代谱估计方法，而 BT 法与周期图法称为传统谱估计法。 无论哪种方法，每一种谱估计技术都可以认为是一种模型法。具体地说，就是根 据对过程的先验知识建立一个近似实际过程的模型。其次，利用观察数据或自相关函 数来估计假设的模型参数，最后做谱估计。常用的模型有：周期图与 BT（正弦谐波总 和模型）、AR、MA、ARMA、Prony、最大似然等。各种估计之间的性能变化，就是 由于假设的模型与实际过程拟合的好坏不同引起的。虽然不同模型可以产生类似的结 果，但是有些模型所需要的参量可能相对少一些。因此，从表征过程的角度来看，这 些方法就更为合适。
∫ ∫ B
=
−
1
( 2π
)2
π −π
π
ω1ω2W (ω1,ω2 ) dω1dω2
−π
∫ ∫ V
=
1
(2π )2
π −π
π
W (ω1,ω2 ) 2dω1dω2
−π
Sasaki 等[1975]证明了双谱近似正比于指数 V，双谱偏差上确界正比于指数 J，
∫ ∫ J
=
1
( 2π
)2
π −π
π
(ω1 − ω2 )2W (ω1,ω2 ) dω1dω2
4.3.2 矩谱的直接估计
设{x(1)" x(N )} 是平稳随机或确定序列的观测值。假设采样周期T
= 1 ，Δ n
Δ
=
1 Nn
是
授课教师：姬红兵教授 hbji@xidian.edu.cn
59
更新日期 2011 年 4 月 1 日
研究生课程：现代信号处理－高阶统计量分析 课程编号：0211007（博）0221023（硕） 西安电子科技大学
（4.6）
其中， d (τ ) 是一维函数，具有性质：
d (τ ) = d (−τ )
d (τ ) = 0 ， τ > Ln
d (0) = 1
D(ω) ≥ 0 ， ω ≤ π
（4.7）
然而，并不是所有的一维窗函数都满足条件，即对于任意ω ， D(ω) ≥ 0 。例如：
Hanning 窗在频域具有负旁瓣。