现代谱估计计算机仿真实验报告
胡敏
在许多工程应用中，利用观测到的一组样本数据估计并分析一个平稳随机信号的功率谱密度是十分重要的。例如，在雷达信号处理中，由回波信号的功率谱密度、谱峰的宽度、高度和位置，可以确定目标的位距离和运动速度；在阵列信号处理中，空间功率谱描述了信号功率随空间角度的分布情况。在许多信号处理应用中，谐波过程经常会遇到，它对应的功率谱为线谱，谐波过程的功率谱估计就是要确定谐波的个数，频率和功率（合称谐波恢复）。为了更好的学习现代信号处理中该部分的内容，我们做了相应的计算机仿真实验。
1 实验目的
1、深入理解现代谱估计和谐波恢复的基本理论，包括ARMA 模型，ARMA 谱估计，ARMA 模型识别，Pisarenko 谐波分解，信号子空间和噪声子空间，旋转不变技术(ESPRIT)；
2、熟悉与上述谱估计和谐波恢复理论相关的数学方法以及各自的特点，包括最小二乘估计（LS ），奇异值分解（SVD ），总体最小二乘估计（TLS ），特征值分解和广义特征值分解；
3、体会ARMA 功率谱估计中的Cadzow 谱估计子和Kaveh 谱估计子，ARMA 模型的识别方法，Pisarenko 谐波恢复方法，ARMA 建模谐波恢复方法，MUSIC 方法进行谐波恢复，两种ESPRIT 方法（LS-ESPRIT 和TLS-ESPRIT 进行谐波恢复；
2 实验原理
2.1 ARMA 谱估计
相当多的平稳随机过程都可以通过用白噪声激励线性时不变系统来产生，而线性系统又可以用线性差分方程进行描述，这种差分模型就是自回归—滑动平均（ARMA ）模型。而且，任何一个有理式的功率谱密度都可以用一个ARMA 随机过程的功率谱密度精确逼近。ARMA 随机过程定义为服足下列线性差分方程的离散随机过程{})(n x ：
∑∑==-+=-+q
j j p i i j n e b n e i n x a n x 1
1
)()()()( （1）
式中)(n e 是一离散白噪声；式（1）所示的差分方程称为ARMA 模型，系统p a a Λ,1和q b b ,,1Λ分别称为自回归（AR ）参数和滑动平均（MA ）参数，而p 和q 分别叫做AR 阶数和MA 阶数。ARMA 过程{})(n x 可以简记为ARMA ),(q p ，使用移位算子1
-z 可以把它写作如下形式：
)()()()(n e z B n x z A = （2）
其中： p p z a z
a z A --+++=Λ1
11)(
q q z b z b z B --+++=Λ111)( 若),0(~)(2
σN n e ，则平稳ARMA ),(q p 过程的功率谱密度为：
2
22
222
)
()()
()()(ω
ωσσωω
j j e z x e A e B z A z B P j === （3）
用（3）式进行谱估计必须事先辨识ARMA 模型和激励噪声的方差2
σ，而MA 参数的估计需要求解非线性方程。为了避开非线性运算，Cadzow 和Kaveh 分别提出了一种线性谱估计子：
1、Cadzow 谱估计子
ω
ω
σ
ωj j e z e z x z A z N z A z N z A z A z B z B P =--=--+
==)
()()()()
()()
()()(11112
（4）
∑=-=p
i i i z n z N 0
)( (5)
∑=-=p
i i k i k a n 0
)(ρ p k Λ,1,0= (6)
⎩⎨
⎧≠==0
)(0
2/)()(k k C k k C k x x ρ (7)
其中)(k C x 为{})(n x 的协方差函数。因此用Cadzow 谱估计子只需要确定AR 阶数和AR 参数就能进行ARMA 谱估计。
2、Kaveh 谱估计子
ω
ω
σ
ωj j e z q
q
k k
k
e z x z c z A z A z A z A z B z B P =-=--=--∑==)()(1)
()()
()()(1112
（8）
∑∑
=+-=p
j x
j
i
p
i
k j i k C
a a c 0
)( q k Λ,1,0= （9）
Kaveh 谱估计子需要确定AR 阶数，AR 参数和MA 阶数来进行ARMA 谱估计。
3、AR 阶数和参数的确定
对于一个ARMA ),(q p 过程，可以推出其相关函数满足如下方程：
∑∞
=+=0
2
)()()(i x i h i h R τσ
τ (10)
其中)(n h 为系统的冲激响应，根据其定义可以得到：
n
q
k k
p i i
b
k n b i n h a =-=-∑∑==0
)()(δ （11）
由（10）式和（11）式就能推导得到著名的修正Yule-Walker 方程，简称MYW 方程：
0)()(1
=-=∑=p
i x i x i l R a l R ，q l >∀ （12）
若已知AR 阶数p ，就能通过求解p 个MYW 方程来求解AR 参数：
r Ra -= （13）
其中：⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡-+-+-++-+-=)()2()1()2()()1()1()1()
(q R p q R p q R p q R q R q R p q R q R q R x x x x x x x x x ΛM M M M ΛΛR []
T
p a a a ,,,21Λ=a
[]T
x x x p q R q R q R )(,),2(),1(+++=Λr
该方程可以用最小二乘法估计a 的值：
r R R R a T T 1^
)(-= （14）
而实际问题中，AR 阶数往往是未知的，此时可用奇异值分解法确定AR 阶数，总体最小二乘估计AR 参
数，合称SVD-TLS 算法。考虑超定方程：
0a R e e = （15）
其中：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-+-++-+++-++=)()1()()2()1()2()1()()
1(e e x e x e
x e e x e x e x e e x e x e x e p M q R M q R M q R p q R q R q R p q R q R q R ΛM M M M ΛΛR []
T
p p p e e
a a a a a ,,,,,,,1121ΛΛ+=a
e p M >>，p p e ≥，q q e ≥。若p q p q e e -≥-，就可以通过对e R 奇异值分解：
H e ΣR V U = （16）
∑中包含1+e p 个奇异值，将其归一化：
11/σσσkk def
kk =-，11+≤≤e p k （17）
选择一个接近于零的数作为阈值，把-
kk σ大于此值得最大整数k 作为有效秩p ，它就是AR 阶数。 根据总体最小二乘方法可以得到矩阵：
H i
j p j p
n i i j jj p v v )(1
11
2)
(∑
∑=-+==σS
（18）
其中：
v
i j
是酉矩阵V 第j 列的一个加窗段，定义为：
T
i
j j p k v j i v j i v v )],(,),,1(),,([++=Λ （19）
