第9章 信号建模和参量谱估计

功率谱
j 2πmf k
R x (ω ) =
∑
i =1
P
2 π | α i | 2 δ (ω − ω i ) + ρ w
自相关矩阵
R x = E{x(n)x H (n)} = R s + R w
2 2 = ∑ | α p | 2 v( f p ) v H ( f p ) + σ w I = VAV H + σ w I p =1 P
Q w = [q P +1 L q M ]
信号特征向量生成的信号子空间 噪声特征向量生成的噪声子空间
信号子空间和噪声子空间
从M维空间到一组向量Z=[z1 z2 … zL]所生成的 L(L<M)维子空间的投影矩阵
P = Z(ZH Z)−1 Z H
Q Qs = I
H s
Ps = Q s Q
H s
Q Qw = I
rank( Rs)<=P
秩为P的Rs性质： 如果M>P， Rs奇异； 如果M<=P， Rs正定
自相关矩阵的特征值分解
H R x = ∑ λmq mq m = QΛ ⋅ Q H m =1 M
其中，λm是以降序排列的特征值， λ1≥ λ2 ≥ …≥ λm qm是它们相应的特征向量， Λ是由降序排列的特征值构成的对角阵，Q的列为对应的特征向量；
第9章 信号建模和参量谱 估计
林耀荣
9.1 建模过程：理论和实际
输出信号的功率谱密度PSD
R (e ) = σ
jω 2 w
1 + ∑ d k e − jωk 1 + ∑ a k e − j ωk
k =1 k =1 Q
Q
2
| D ( e − jω ) |2 =σ | A(e − jω ) |2
2 w
m = P +1 λm
∑
M
1
v ( f )q m
H
2
=
1
m = P +1
∑
M
1
λm
| Qm (e j 2πf ) |2
根-MUSIC算法
分母的z变换
1 Pmusic ( z ) = ∑ Qm ( z )Q m ( * ) z m = P +1
* M
每个噪声特征向量伪谱z变换的和
2M-1阶多项式有M-1对根； 其中，离单位圆最近的P个根对应复指数， 根的相位为估计的频率
| αl | e
− j [ 2π ( n − m ) f l +ψ l ]
}
2 Rx (m) = ∑ | α k |2 e j 2πmf k + σ wδ (m) k =1
P
1 = ∑| αk | 2π k =1 = ∑| αk | e
2 k =1 P
π
e j 2πmf k dψ k
Fourier变换
9.6 谐波模型和频率估计技术
白噪声激励的线性时不变系统输出-零极点模 型 白噪声中的正弦函数或复指数信号-谐波模型 谐波模型信号的实例
语音处理中的共振峰频率 雷达中的移动目标回波 阵列处理中的空间传播信号
9.6.1 谐波模型
谐波模型
实值信号：复共轭对（正弦）形式的复指数信号 复值信号：单一频率点
H 2
=
1 | QM (e j 2πf ) |2
P个峰值对应估计的频率，称为伪谱
9.6.2 Pisarenko谐波分解
QM (e
j 2 πf
) = v ( f )q M = ∑ qM ( k )e j 2πf ( k −1)
H k =1
M
对应最小特征值的第M个特征向量的傅里叶变换
伪谱可以通过计算qM的FFT有效实现 或者，通过计算第M个特征向量多项式的零点实现：
=0 P<m≤M
每个噪声特征向量的伪谱
Rm (e
j 2πf
)=
1 v ( f )q m
H 2
1 = | Qm (e j 2πf ) |2
伪谱有M-1个峰值； 其中P个对应复指数的频率； 其余M-1-P个为噪声引起的伪峰值； 不同的噪声特征向量的伪峰值位置不同；
9.6.3 MUSIC算法
减小伪谱中的伪峰值-平均M-P个噪声特征向量的伪谱
为时间窗频率向量
自相关矩阵
R x = E{x( n)x H ( n)} = R s + R w
Rs (m) = E ( x(n) x ∗ (n − m)) = E{∑ | α k | e
k =1 P j ( 2πnf k +ψ k )
白噪声的自相关矩阵
∑| α
l =1
P
l
|e
− j [ 2π ( n − m ) f l +ψ l ]
时间窗口的长度：M=P+1 因此，噪声子空间只包含一个特征值：QM=qM ，对应最小的特征值λM 由于信号和噪声子空间正交
v ( f p )q M = ∑ qM ( k )e
H k =1 M j 2πf p ( k −1)
=0 m≤P
可知：
R phd (e j 2πf ) =
1 v ( f )q M
频率估计技术
估计复指数信号的频率
功率谱估计
估计的频率-频谱中的峰值频率 非参数方法估计功率谱 全极点模型谱估计方法 没有考虑噪声中复指数信号的模型
基于谐波模型的频率估计方法-高分辨率方法
Pisarenko谐波分解 MUSIC 最小范数 ESPRIT算法
9.6.1 谐波模型
长度为M的确定时间窗 x(n) = [ x(n) x(n + 1) ... x(n + M − 1)]T
MATLAB
[Q0,D]=eig(R ); %特征值分解 [lambda, index] = sort (abs(diag(D) ) ); %特征值升序排列 lambda = lambda(M:-1:1); Q = Q0(:,index(M:-1:1)); Rphd = 1./abs(fftshift(Q(:,M), NFFT))).^2
2 因为 R x = VAV H + σ wI
特征值可以分为2部分：
λm = M | α m | +σ , m ≤ P
2 2 w
P个最大的特征值对 应复指数信号和噪声
2 λm = σ w , m > P
剩下的特征值对应噪声
信号子空间和噪声子空间
将相关矩阵按信号和噪声的特征向量分成两部分
R x = ∑ ( M | α m | +σ )q mq +
谐波模型 噪声中的P个复指数的信号模型 其中
fp =
x ( n) = ∑ A p cos(ω p n + φ p )
p =1
P
x ( n) = ∑ α p e
p =1
P
j 2πnf p
+ w( n )
ω p Fp = 2π Fs
α p =| α p | e
jψ p
ψp是[0, 2π]内均匀分布的随机变量 w(n)为白噪声 |αp|、fp为确定量，待估计 复指数的功率谱为线性谱
L L M L M L L
⎤ ⎥ x( M ) ⎥ ⎥ M ⎥ x(n + M − 1) ⎥ ⎥ M ⎥ x( N + M − 3) ⎥ x( N + M − 2)⎥ ⎦
x( M − 1)
相关矩阵的估计-采样相关矩阵：
ˆ = 1 XH X Rx N
9.6.2 Pisarenko谐波分解(M=P+1)
}
2 R w = σ wI
1 = 2π
∫π
− P
π
| α k | e j ( 2πnf k +ψ k ) | α k | e − j[ 2π ( n − m ) f k +ψ k ]dψ k ∑
k =1 P
P
+ E{∑
P
k =1 l =1,l ≠ k 2
∑| αk | e ∫π
−
j ( 2πnf k +ψ k )
其中 V = [ v ( f1 ) v ( f 2 ) L v ( f p )] MxP矩阵，列是频率fp的复指数信号的时间窗频率向量
v H ( fi ) v( f j ) , i ≠ j
⎡| α1 |2 0 L 0 ⎤ ⎢ ⎥ 2 | α2 | O M ⎥ ⎢ 0 A= ⎢ M O O 0 ⎥ ⎢ 2⎥ L 0 | αP | ⎥ ⎢ 0 ⎣ ⎦
Rmusic (e j 2πf ) = MUSIC伪谱 1
m = P +1
v H ( f )q m ∑
M
2
=
1
m= P +1
| Qm (e j 2πf ) |2 ∑
M
MUSIC伪谱假设：所有噪声特征值相等 为了解决可能存在不相等噪声特征值的问题： 特征值归一化每个特征向量的伪谱-特征向量方法
Rev (e j 2πf ) = 1
x ( n ) = ∑ α p v ( f p )e
p =1 P j 2πnf p
+ w ( n) = s( n) + w ( n)
其中
w (n) = [ w(n) w(n + 1) ... w(n + M − 1)]T 为白噪声的时间窗向量
v ( f ) = [1 e j 2πf
L e j 2π ( M −1) f ]T
H w
H Pw = Q wQ w
特征向量是正交的
PwQ s = 0
Pw v ( f p ) = 0
Ps Q w = 0
时间窗口频率向量完全在信号子空间内
Ps v ( f p ) = v ( f p )
采样相关矩阵
x(1) ⎡ xT (0) ⎤ ⎡ x(0) ⎥ ⎢ ⎢ T x(2) ⎢ x (1) ⎥ ⎢ x(1) ⎥ ⎢ ⎢ M M M ⎥ ⎢ ⎢ T X = ⎢ x ( n) ⎥ = ⎢ x ( n) x(n + 1) ⎥ ⎢ ⎢ M M M ⎥ ⎢ ⎢ T ⎢x ( N − 2)⎥ ⎢ x( N − 2) x( N − 1) ⎢xT ( N − 2)⎥ ⎢ x( N − 1) x( N ) ⎦ ⎣ ⎣
9.6.3 MUSIC算法 (M>P+1)
时间窗口的长度：M>P+1 噪声子空间的维数大于1 对较大维数的噪声子空间平均，更强壮的频率 估计方法
9.6.3 MUSIC算法
Pw v ( f p ) = 0
v ( f p )q m = ∑ qm (k )e
H k =1 M j 2πf p ( k −1)
2 m =1 2 w H m H 2 H = Q s Λ s Q s + σ wQ wQ w P m = P +1 2 H σ wq m q m ∑ M
Q = [q1 q 2 L q P q P +1 L q M ]
Qs ⊥ Qw
H Qs Qw = 0
Qs
Q s = [q1 q 2 L q P ]
Qw
QM ( z ) = ∑ qM (k ) z
k =1 M −k M −1 k =1
= ∏ (1 − e j 2πf k z −1 )
多项式的P=M-1个根的相位对应P个复指数的频率fk
Pisarenko谐波分解的缺点ຫໍສະໝຸດ Baidu
对噪声敏感
仅用最小特征值对应的特征向量估计复指数频率 对这个噪声特征向量估计中的任何误差都很敏感