有理谱估计的参数化方法

非参数化方法对研究的信号除了平稳性假设外没有作任何假设。参数化或基于模型的谱估计方法均假设信号满足函数形式已知的模型，然后对假设中模型的参数进行估计。从估计模型中可以获得感兴趣的信号谱特性。显然，当假设模型与实际非常接近的情况下，参数化方法相对于非参数化方法能提供更为精确的谱估计值；但是，在研究信号的信息极少甚至没有的应用中，功率谱密度估计的非参数化方法仍然有用。

谱估计参数化方法分为两种——连续谱和离散谱。本章的有理谱参数化方法是连续谱类中研究的热点，我们主要讨论参数的估计方法，需要知道或选择谱模型的结构或阶数。

l 简约原理（Parsimony Principle）：Better estimate may be obtained by using an appropriate data model with fewer unknowns.

l 奥卡姆剃刀（Occam's razor）：one should not make more assumptions than the minimum needed

如果模型选错了，功率谱的估计将永远有偏。所以要想使用参数化方法，拥有对数据模型的正确的先验知识很有必要。

4.3.1 有理功率谱

有理谱密度为的一个有理函数，即两个关于的多项式比值：

其中，，。由微积分中的Weierstrass定理可知，当上式中阶数p、q选的足够大时，任意的连续功率谱密度都可以由其逼近，这也是有理谱成为连续谱类热点的原因，此模型引起很多研究者的兴趣。

由于，由上面所示模型建立功率谱密度可以进行如下形式分解

是一个标量，A和B是多项式：

模型也可以在Z域得到类似表示：

因此，可以把谱分解为：

其中，

l 的零点和极点均关于单位圆成对出现

l 以上内容被称为谱分解定理（Spectral Factorization）

4.3.2 三个经典的信号模型：ARMA，AR，MA

4.3.2.1 模型的基本结构

任意有理谱密度能和一个信号联系起来，该信号可由功率为的白噪声通过一个传输函数为

的有理滤波器得到。

利用。该滤波器在时域可以表示为

其中，是方差为的白噪声。因此，通过谱分解定理，谱的参数化模型可以变为信号自身的模型à谱估计问题可以简化信号的建模问题。

满足的信号x(n)称为自回归滑动平均模型（Autoregressive Moving Average），记为ARMA或信号。如果B(z)=1，即q=0，那么x(n)为自回归信号，记为AR或；如果A(z)=1，即p=0，那么x(n)为滑动平均信号，记为MA或。

l 假设对任意，均有限，那么没有任何零点恰好在单位圆上

l 若模型中A(z)的所有零点（模型的极点）都严格在单位圆内，则对应的模型是稳定的

l 若模型中B(z)的所有零点（模型的零点）都严格在单位圆内，则对应的模型具有最小相位。

l AR模型对有尖峰的建模较好；MA模型对有波谷的建模较好；ARMA能对既有尖峰又有波谷的进行建模。

由无穷级数的性质可知：

l 一个或模型可以等价于一个无穷阶数的AR模型

l 一个或模型可以等价于一个无穷阶数的MA模型

4.3.2.2 模型的协方差结构

我们期望推导由参数a、b和表示的ARMA过程协方差表达式。这种表达式使用由数据得到的协方差估值代替协方差的真值，为ARMA参数估计提供了简便方法。

式可以写为：

将上式与相乘并取数学期望得：

由于滤波器是渐近稳定和因果的，所以可得：

因此

因此ARMA模型求完期望之后右面变为：

因此协方差结构为：

通常，是关于系数a，b的非线性函数，由于s<0时，因此对于时，上式可化简为如下AR系数估计器：

4.3.3 AR信号

在ARMA类中，自回归（全极点）信号是应用中最常使用的一种类型，通过令中，A多项式的零点接近于单位圆，AR方程可以模拟具有尖锐峰的谱。这是一个非常重要的特性，因为窄带谱在现实中很常见。另外，AR模型参数估计问题可以通过解线性方程组得到，并能确保估计得到的AR多项式的稳定性。

AR谱有两种估计方法，第一种是基于上一节推出的协方差和AR参数线性关系的Yule-Walker方法，第二种是使用时域方程得到的AR参数的最小二乘解。在第二种方法中，我们将看到，“最小二乘方法”与线性预测问题的紧密联系。

4.3.3.1 Yule-Walker方法

对于AR信号，q=0，B(z)=1，因此当k>0时，有如下关系成立：

联立上式，将方程组写成矩阵形式：

称上式为Yule-Walker方程或正则方程，如果相关序列r已知，则可以求解待估计参数

对正则方程进行变换，去掉第一行

拆分

其解为，一旦求出，那么由被去掉的第一行就能求出。

l 当我们只有有限个数据时，只能求得的估计值，可以使用估计值获得对的估计值。

l 是正定矩阵则求解出来的由组成的多项式A(z)的根在单位圆之内

l 是Toeplitz&Hermitian矩阵，能利用结构寻找快速算法

l 为强调和对阶数p的依赖关系可以将正则方程写成

4.3.3.2 Yule-Walker方程的阶递推解法

因为在多数应用中，缺少真正关于阶数n的先验信息，所以可以选择遍历不同阶数的AR模型——对求解Yule-Walker方程组，使用一般求解的方法，计算量要达到，要减小计算量，需要充分利用上一节讨论的特殊代数结构。

Yule-Walker方程组中的是高度结构化的，可以使用Levinson-Durbin算法(LDA)以复杂度完成的求解计算；此外Delsarte-Genin算法（DGA），也称Split-Levinson算法在实信号下比LDA快一倍。

两种算法都是以阶数p的递推形式求解，且仅要求其中的矩阵是正定、Hermitian和Toeplitz的。

4.3.3.3 最小二乘方法

最小二乘方法是基于时域关系的最小二乘（Least-squares，LS）最小化准则，可以通过与AR估计器关系密切的线性预测问题来推导最小二乘估计器。

l 最小二乘思想：预测值和真实值在向量空间中欧式距离最小，且最小时误差向量与真实值正交。

首先把Yule-Walker方程和线性预测问题联系起来，假设是p阶AR过程，那么满足：

其中，我们认为是由前p个样本的线性预测，为相应预测误差。

在非理想情况下可能无法获知参数的真值，此时满足：

由最小二乘LS的思想可知，使预测误差方差最小化的矢量即为AR系数向量，此时满足，方差仅由白噪声产生，而无人为预测误差。

由之前在中的定义

对上式中求导，并令求导结果等于0，

得最小值对应的为：

与之相对应的最小预测误差为：

上面两个式子其实就是之前所示得Yule-Walker方程对应的第一行和拆出来得低一阶得正则方程：

由此可见，Yule-Walker方程可以看作是利用之前最近的p个样本来求其最佳线性预测问题的解。为此AR模型有时也被称作线性预测模型。

最小二乘AR估计方法就是基于上述最小化问题的有限样本近似解。已知一有限的观测集，我们通过有限样本代价函数近似最小化

假设当n<1和n>N时x(n)=0，最小化的矢量为：

从代价函数可以看出，和的定义依赖于的选择。

l 若得到所谓的自相关法

l 若即去除了式中所有的零值，得到所谓协方差法，又叫协方差LS法或LS法

l 若，叫做前加窗方法

l 若，叫做后加窗方法

AR估计中最小二乘LS和Yule-Walker的区别与联系

注意到，将中和稍作处理就能得到和有限样本估计值，由此最小二乘方法可以认为是Yule-Walker方程的近似解。

l 有偏估计值和，是Toeplitz矩阵

l 无偏估计值和，不是Toeplitz矩阵，不可以用LDA和DGA求解

因此，LS自相关法与Yule-Walker法等效，LS协方差法与Yule-Walker法的区别随着N的增加而减少。

l 对于大样本来说AR参数的Yule-Walker方法与LS估计几乎一致；对于中小样本来说，YW与协方差LS表现略有不同，由YW方法估计的AR模型总能保证是稳定的，而估计的LS模型可能不稳定。

l LS方法比YW方法更精确，即LS比YW估计出的参数更接近真实值，这个原因暂时还在研究中，一种假说是在之外和中对应的零元素导致YW估计有偏，同时N不远大于n时，偏差可能很大。

4.3.4 求解超定方程的LS计算方法

4.3.4.1 Frobenius norm

Let denote the vector space of all matrices of size (with rows and columns) with entries in the field .

For all matrix in ,

where are the singular values of .

4.3.4.2 线性最小二乘（LLS）

直接求解正则方程，得

4.3.4.3 数值解对噪声更鲁棒的QR Decomposition LS Mehtod

使用辅助变换（householder transformation），我们能找到一个正交的酉矩阵Q，使得，

因此

其中，R是一个上三角方阵,最小二乘解为：

4.3.4.4 总体最小二乘（Total Least Square）

l 在LLS中，得到最小二乘解得方法是使得中

l 在TLS中，得到最小二乘解得方法是使得中

低信噪比条件下，TLS比LS估计的准确；高信噪比下，二者结果接近。

TLS和RLS的图解

4.3.4.5 使用奇异值分解（SVD）计算最小二乘LS

对A进行SVD分解：

0矩阵使得矩阵A可以非满秩，若A是满秩的，则，其中

叫做矩阵A的特征值，r是矩阵A的秩。最小二乘解如下：

4.3.4.6 使用奇异值分解（SVD）计算总体最小二乘TLS

表示矩阵的奇异值分解

将分块表示为

4.3.5 MA信号

MA模型基于以下几点原因在谱估计实际工程中使用受限：

l MA模型具有全零点结构，阶数低时很难模拟窄带尖峰谱

l MA适合模拟的平坦峰和尖锐的谷（零点），但实际应用中这种谱不常见

l MA参数估计问题通常是非线性的，求解起来比线性的AR参数估计要困难

l MA和ARMA估计问题面临的问题很相似，因此更常采用较为一般的ARMA模型代替MA

4.3.6 ARMA信号

4.3.6.1 修正的Yule-Walker方法

用于估计ARMA谱密度的修正YW方法由两部分组成：

第一步

利用估计AR系数

对于，将写成矩阵形式：

如果Q=p，则上式为由n个方程构成的含有n个未知变量的方程组，这就是在AR情况下成立的广义YW方程组，称这些方程构成了修正的Yule-Walker方程组。在这之中若使用样本估值代替理论协方差，可以得到：

由这些方程组可以解出称它们是的修正Yule-Walker估值，上式的方阵非奇异，具有快速算法。

l 增加Q值，充分利用信号更长时间的协方差的意义：若ARMA模型中的零点严格位于单位圆内，，那么由Modified Yule-Walker估计的参数具有合理的精确度，然而，当ARMA模型中的零点和极点间隔很近且在单位圆附近时，由MYW得到的估计值非常不准确，这种具有几乎重合的零极点且零极点模值接近1的ARMA模型与窄带信号相对应。由于信号在频域越集中，时域就越发散，因此窄带信号的协方差序列衰减很慢。这意味着我们增加延迟序列协方差矩阵的信息能提高系数估值的精确度。

通过选择Q>p就可以利用附加信息，并且由此得到超定方程组。由于超定方程组一般不存在精确阶。因此此时应用LS或者TLS求解方程：

这里的‘hat’符号表示由ACS估计值代替真实值后的ACS矩阵和矢量。上式的加权最小二乘解可由下式给出：

为阶正定加权阵，Q>p的情况下，由导出的AR估值称为超定MYW估值。

加权阵的适当选择也能够增加AR系数估值的精确度，当时可以得到常规最小二乘估值，精确度的提高可以通过选择对角元素为正且递减的对角阵来得到，这样的反映了高ACS延迟股指中降低的置信度。

STOICA等人已经推导出了最优加权阵，在所有可能的中该最优加权使的协方差最小。可惜，该最优加权依赖于ARMA的参数，因此，为了使用最优加权法，需要使用两步“bootstrap”。

l 第一步，选择一个固定的来得到初始参数估计值

l 第二步，使用这些初始估值来获得最优化

最优加权的收益相比计算最优加权矩阵的计算开销而言较小，因此，精确度的提高主要通过时增加Q值来实现。

第二步

用估计的AR系数和ACS估计值在中估计系数

一旦得到了AR估值，我们就可以转到ARMA谱中MA部分的估计问题，因为MA部分功率谱密度可由下式给出：

因此估计出就足以描述MA的谱。

令当时

表示MA部分的协方差，考虑

可以得到：

将上一步计算出的和的估值代入上式，得到估计器：

最后ARMA的谱公式如下：

由于不能保证上式中的分子对所有全为正，因此这种方法可能导致负的ARMA谱估值

4.3.6.2 两步最小二乘方法（Two Stage Least Square Method）

如果噪声序列已知，那么ARMA模型中的参数估计问题将是一个简单的输入输出系统参数估计问题（系统辨识）。最简单的方法就是最小二乘方法

第一步：得到噪声序列估值

使用阶数L较高的来拟合，可以使用YW方程来求解估计

因此噪声序列可以估计为：

第二部：系统辨识

令

，，

观测方程：

则最小二乘解为：

QRLS也同样适用于此最小二乘解法，这个方法在ARMA的零点在单位圆附近时，产生不稳定解。

4.3.7 超分辨谱分析——MUSIC算法和Pisarenko算法

核心思想：接收到的信号中要么是真实信号，要么是噪声，那么，我们只要能找到噪声并将之排出，剩下的就是信号了。

MUSIC算法假定输入序列是由一组角频率为的正弦信号组成的，那么真实角频率必是如下方程的解

其中，

是由天线阵结构决定的，阵元数为M，阵元数也可以看作是延迟单元的个数。假设阵元接收的是长度为M的序列，MUSIC算法先计算其的自相关矩阵

n从M开始是为了保证中没有0值。

接下来对自相关矩阵R作特征值分解（Eigendecomposition），分解完后由R的特征向量组成如下矩阵的列，并将之划分为四块：

左上角大小的方阵对应K个较大特征值，是信号；右下角大小的仿真对应（M-K）个较小的特征值，是噪声。

就是由特征向量为列组成的方阵的右下角矩阵。

相当于将信号与噪声作互相关，即是互相关的功率谱。

由于要利用接收序列中的噪声，因此正弦信号频率分量个数K必须要比阵元数M少才可以，不能等于，否则无法分解出；当阵元数仅比信号数多一个的时候，即时，此时是MUSIC算法的特例，叫做Pisarenko算法（Pisarenko harmonic decomposition）。

4.3.8 仿真对比各种参数化谱估计方法

l 原始信号的周期图

l Yule-Walker方法

l Levinson-Durbin递推方法

l Burg方法

l 协方差法

l 改进的协方差法

l 超分辨MUSIC方法