高阶谱分析

φx
(1)
=
N
∑
i=2
⎧⎪ ⎨ ⎪⎩
s
(
i
)
i
−
(i
s(i − −1)
1)
⎫⎪ ⎬ ⎪⎭
+
φx
(N
N
)
(6.7)
ω2 π
k k-i
0i kπ
ω1
图 6.1 在 Brillinger 算法中的双谱相位线
从(6.6)式可见，Fourier 相位φx (ω ) 是沿图 6.1 所示的对角双谱相位线求和计算的。
6.3.1 Brillinger 算法（1977）
Brillinger 提出的相位恢复算法是一种利用全部双谱值的递推方案。假设ω1 取离散
值 0,1,…,k，ω2 取 k，k-1,…,0，则利用(6.1)式可生成(k+1)个方程，将这些方程求和，得
k
k
k −1
∑
Ψ
x 3
(i,
k
−
i
)
=
2∑
φx
(i
)
下面给出两种双谱相位解模糊的方法。 方法Ⅰ[Marron 等 1990 年]
设 Ψˆ 3x (ω1,ω2 ) 表示模 2π 双谱相位，即有
复以上过程直到求出所有相位值。将 j = 2,3, , N 分别代入(6.10)式，有类似的相位关
系。因此，每个相位φx ( p) 有 ( p −1) 2 个独立的表示（对于 p 为奇数），或 p 2 个表示
（对于 p 为偶数）。 为了改善表示的信噪比（SNR），可以对这些表示取平均。注意，应对指数因子
⎢
⎥
A
=
⎢1 ⎢⎢0
0 2
0 0
0 −1
0 0
... ...
1⎥ 0⎥⎥
,
⎢0 1 1 0 −1 ... 0⎥
⎢ ⎢
.
.
.
.
.
...
.
⎥ ⎥
⎢ . . . . . ... . ⎥
⎢ ⎢
.
.
.
.
.
...
.
⎥ ⎥
⎢⎣ . . . . . ... 0⎥⎦
其中，矩阵的维数：若
N
为
偶
数
，
为
⎛ ⎜⎜⎝
⎛ ⎜⎝
N 2
的限制。因此，
{ } exp{ jφx ( p)} = (const.)⋅ ∑ exp
j
(φx
(
p
)
+
φx
(
p
−
q
)
−Hale Waihona Puke Baidu
Ψ
x 3
(
p
−
q,
q
)
)
q
(6.13)
B-L-W 算法虽然也是递推型的，但不要求解相位模糊方法，且利用了全部双谱值。
6.3.4 Matsuoka-Ulrych 算法（1984）
以上讨论的算法都是递推型的。Matsuoka 和 Ulrych 在 1984 年提出了一种非递推
k −1
φˆx (k ) −φx (k ) = −∑ 2π ⋅n (i,1)
(6.9)
i=0
其中 n (⋅) 是一整数函数。换句话说，φˆx (k ) 和φx (k ) 相差 2π 的整数倍。
授课教师：姬红兵教授 hbji@xidian.edu.cn
83
更新日期 2011 年 4 月 25 日
研究生课程：现代信号处理－高阶统计量分析
6.2 利用高阶谱估计信号的幅度和相位
对于一个线性非高斯信号，已知其高阶谱即可恢复它的 Fourier 幅度（在一个尺度 因子内）和相位（差一个线性相移）。本节讨论假设已知双谱幅度和相位时，恢复信号 Fourier 幅度和相位的算法。这些算法可以直接推广到如三谱等高阶谱的情况。
设信号 X (k) 的 FT 为 X (ω) ，则其双谱相位 Ψ3x (ω1,ω2 ) 和 Fourier 相位φx (ω ) 间有以
研究生课程：现代信号处理－高阶统计量分析
课程编号：0211007（博）0221023（硕） 西安电子科技大学
第六章 利用高阶谱恢复信号的非参数方法
6.1 引言
利用高阶谱恢复信号的非参数方法是指利用双谱或三谱恢复信号的 Fourier 幅度和 相位的方法，这些方法无需用参数模型（如 AR、MA、ARMA 等）适配数据来求解问 题。非参数方法需要信号双谱和三谱的先验信息。本章研究几种不同的信号恢复的非 参数算法，包括： （1） 相位恢复算法：递推和非递推； （2） 基于多倒谱的相位和幅度恢复算法； （3） 仅用双谱相位的信号恢复方法； （4） 双通道盲解卷积中的信号恢复。
容易证明，如果利用模
2π
相位
Ψˆ
x 3
(ω1,ω2
)
代替真实双谱相位，则
(6.6)
式和(6.7)
式将
产生有误差的 Fourier 相位φˆx (ω ) ，即φˆx (ω ) −φx (ω ) ≠ 2π k (ω ) 。
Brillinger 算法是递推的，因此对误差敏感，特别在初始相位值φx (1) 的估计中。
(6.10)
因此，对于 j = 1 ，有
φx
(
i
)
=
φx
(1)
+
φx
(
i
−1)
−
Ψ
x 3
(
i
−
1,1)
(6.11)
用初始条件φx (0) = 0 ，φx (1) 为任意值。根据(6.11)式，可恢复相位如下
φx (2) = φx (1) +φx (1) − Ψ3x (1,1)
(6.12.1)
φx
(3)
ω2 π
ωH
ωL
1
0 ωL ωH π
ω1
图 6.2 计算带限信号相位所能用的双谱值
6.3.3 Bartelt-Lohman-Wirnitzer 算法（1984）
Bartelt-Lohman-Wirnitzer 算法基本上是 Lii-Rosenblatt 算法的推广。由(6.1)式可得
φx (i) = φx ( j) +φx (i − j) − Ψ3x (i − j, j)
⎞2 ⎟⎠
×
(
N
−1)
⎞ ⎟⎟⎠
；若 N 为奇数，为
⎛ ⎜ ⎝
(
N
− 1) (
4
N
+
1)
×
(
N
−
1)
⎞ ⎟ ⎠
。
可以证明，A 是满秩矩阵。从(6.15)式可得最小二乘解
( ) Φx = ATA -1 ATΨ3x
(6.16)
该算法利用了所有的双谱值，且是非递推的。因此，它不受累积误差的影响。但
是，当使用模 2π 双谱相位代替真实双谱相位时，(6.16)式不能得到正确的φx (i) 。
课程编号：0211007（博）0221023（硕） 西安电子科技大学
该算法实质上仅限于利用沿直线
Ψ
x 3
(i,1)
，i
=
0,1,
, N 上的双谱值。因此，没有考
虑现有的全部双谱信息，这在带限信号的情况中会产生严重的问题。例如，考虑一带
限信号，其低频和高频截止频率分别为ωL 和ωH 。如果ωL > 1 ，则 Lii-Rosenblatt 算法 将失效，因为所用的双谱值全为零。
−
(
k
+
1)
φx
(
k
)
=
2∑
φx
(
i
)
−
(
k
−
1)
φx
(
k
)
i=0
i=0
i=0
(6.5)
可求得φx (k ) 为
∑ ∑ ∑ φx
(k)
=
k
2 −1
k −1
φx
i=0
(i)
−
k
1 −1
k i=0
Ψ
x 3
(i, k
−
i)
=
k
1 −1
⎡⎢⎣2
k −1 i=0
φx
(i)
−
s
(k )⎤⎥⎦
，k=2,3,…N
(6.6)
ln C3x (ω1,ω2 ) = ln X (ω1 ) + ln X (ω2 ) + ln X (ω1 + ω2 )
(6.3)
比较 (6.2) 式和(6.3) 式，可知两种问题可以用同样的算法，只要 X (ω ) ≠ 0 。
6.3 相位恢复算法
双谱最吸引人的特点是它保留了信号的 Fourier 相位信息，本节讨论几种相位恢复 算法，并说明它们的优点和不足。大部分算法基于(6.1)式。实际上，双谱相位仅可用 下式计算
6.3.2 Lii-Rosenblatt 算法（1982）
Lii-Rosenblatt 算法也是递推的，但它仅利用了一条双谱线上的值。将对应于 ω1 = 0,1, , k −1和ω2 = 1的双谱相位值求和，得
k −1
φx
(
k
)
=
−∑
Ψ
x 3
(
i,1)
+
kφx
(1)
+
φx
(
0
),
k = 2,3,
,N
因此，假设φx ( N ) = 0 将意味着信号将被时移。
利用关于φx ( N ) 的这一假设，φ (1) 可按如下方法计算
授课教师：姬红兵教授 hbji@xidian.edu.cn
82
更新日期 2011 年 4 月 25 日
研究生课程：现代信号处理－高阶统计量分析
课程编号：0211007（博）0221023（硕） 西安电子科技大学
双谱值。然而依据相位恢复算法，这将导致有误差的 Fourier 相位φˆx (ω ) 。对于将要讨
论的算法，从计算双谱相位导出的 Fourier 相位φˆx (ω ) 与真实相位仅差一线性相位项，
即φˆx (ω ) = φx (ω ) + 2π k (ω ) ，其中 k (ω ) 是一整数函数。
线性相位差一般不重要，因为它仅对应一个信号的时移。如果情况不是这样，那 么在应用相位恢复算法前，需要估计双谱相位的二维解相位模糊算法。
(1,
N
−
1)
=
φx
(1)
+
φx
(
N
−
1)
−
φx
(
N
)
Ψ3x (2, 2) = 2φx (2) −φx (4)
………
( ) ( ) Ψ3x
N 2,N 2
= 2φx
N 2
−φx ( N ) ，如果 N 是偶数。
取φx ( N ) = 0 ，将(6.14)式写为矩阵形式
(6.14)
AΦ x
=
Ψ
x 3
其中，
=
φx
(1)
+
φx
(
2
)
−
Ψ
x 3
(
2,1)
……..
φx ( N ) = φx (1) +φx ( N −1) − Ψ3x ( N −1,1)
（6.12.2) (6.12.3)
其中，将φx (1) 代入(6.12.1)式可得φx (2) 。同样，将φx (2) 代入(6.12.2)式可得φx (3) ，重
更新日期 2011 年 4 月 25 日
研究生课程：现代信号处理－高阶统计量分析
课程编号：0211007（博）0221023（硕） 西安电子科技大学
⎡2 −1 0 0 0 ... 0⎤
⎢⎢1 1 −1 0 0 ... 0⎥⎥
⎢ . . . . . ... . ⎥
⎢ ⎢
.
.
.
.
.
...
.
⎥ ⎥
⎢ . . . . . ... . ⎥
方法。如下：将ω1 = 1, 2,..., N 2 和ω2 = ω1,ω1 +1,..., N − ω1 代入(6.1)式，有下面方程组
Ψ
x 3
(1,1)
=
2φx
(1)
−
φx
(
2)
Ψ
x 3
(1,
2
)
=
φx
(1)
+
φx
(
2)
−
φx
(3)
Ψ
x 3
(1,
3)
=
φx
(
2
)
+
φx
(
3)
−
φx
(
5
)
………
Ψ
x 3
Φx = ⎡⎣φx (1),φx (2),...,φx ( N )⎤⎦T
( ) Ψ
x 3
=
⎡⎣Ψ
x 3
(1,1)
,
Ψ
x 3
(1,
2
)
,
...,
Ψ
x 3
(1,
N
−
1)
,
Ψ
x 3
(
2,
2
)
,
...,
Ψ
x 3
N 2,N 2
⎤T ⎦
和
(6.15)
授课教师：姬红兵教授 hbji@xidian.edu.cn
85
(6.8)
i=0
再设φx (0) = 0 ，φx (1) 为任意值，它建立了时间信号的位置。用这些初始条件，即
使用模 2π
双谱相位
Ψˆ
x 3
(i,1) 代替真实双谱相位 Ψ3x
(i,1) ，(6.8)式仍能给出正确的
Fourier
相位估计φˆx (k ) 。
设φx (1) = 0 ，由(6.8)式可得
授课教师：姬红兵教授 hbji@xidian.edu.cn
84
更新日期 2011 年 4 月 25 日
研究生课程：现代信号处理－高阶统计量分析
课程编号：0211007（博）0221023（硕） 西安电子科技大学
exp{ jφx ( p)} 取平均，而不是直接对相位求和，因为对φx ( p) 的每一个表示求和受模 2π
k
其中
s
(
k
)
=
∑
Ψ
x 3
(
i,
k
−
i
)
i=0
注意到，k=N 对应于ω = π 。为了实现(6.6)式中的递推，需要初始值φx (0) 和φx (1) 。
为了确定这两个值，首先假设φx (0) = 0, φx ( N ) = 0 。实际上，φx ( N ) = kπ ，k = 0, ±1, ±2,... 。
授课教师：姬红兵教授 hbji@xidian.edu.cn
81
更新日期 2011 年 4 月 25 日
研究生课程：现代信号处理－高阶统计量分析
课程编号：0211007（博）0221023（硕） 西安电子科技大学
Ψˆ 3x
(ω1,ω2
)
=
arctan
Im Re
⎡⎣C3x ⎡⎣C3x
(ω1,ω2 (ω1,ω2
)⎤⎦ )⎤⎦
(6.4)
其中 Im[i] ， Re[i]分别表示双谱 C3x (ω1,ω2 ) 的虚部和实部。
即使忽略双谱的估计误差，计算的相位
Ψˆ
x 3
(ω1,ω2
)
也与真实相位
Ψ3x
(ω1,ω2
)
差
2π k (ω1,ω2 ) ，其中 k (ω1,ω2 ) 仅取整数值。这一相位差对双谱无关紧要，因为它不影响
下关系
Ψ
x 3
(ω1
,
ω2
)
=
φx
(ω1
)
+
φx
(ω2
)
−
φx
(ω1
+
ω2
)
(6.1)
且幅度之间关系为
C3x (ω1,ω2 ) = X (ω1 ) ⋅ X (ω2 ) ⋅ X (ω1 + ω2 )
(6.2)
问题是如何从
C3x
(ω1,ω2 )
和
Ψ
x 3
(ω1,
ω2
)
恢复
X
(ω )
和φx (ω ) 。对(6.2)式取对数有