【Prony算法】之公式推导及算法步骤

---

1 公式推导过程

Prony方法假定的模型是一组p个具有任意振幅、相位、频率和衰减因子的指数函数，其离散时间的函数形式如下(记为式为4-7-1)：

在上式中，N为采样数据的个数；m为模型阶数，且N≥2m。

式4-7-1用作近似测量数据x(0),…,x(N-1)的模型。更一般地，z_m和b_m假定是复数，其公式如下(记为4-7-2)：

其中的符号为：幅值A_m、频率f_m、相位(单位是弧度)θ_m和衰减因子α_m，Δt代表采样间隔。

使平方误差(记为式4-7-3)：

最小，便可以求出{A_m,θ_m,α_m,f_m}，但解上述这样一个非线性的最小二乘问题是困难的，通常这种求解是一个迭代过程，本文不打算详细介绍。

下面仅讨论{A_m,θ_m,α_m,f_m}的线性估计问题。

首先定义多项式(记为4-7-4)：

根据式4-7-1可知，表示：

的一种方法是：

用a_m乘上式，然后对p+1个乘积求和，即有：

又因为：

所以：

上式记为4-7-5，该式最后等于零是因为第二项求和恰好是式(4-7-4)的位于根z_l处的多项式φ(z_l)，而φ(z_l)=0。【我并不理解本段】

式(4-7-5)意味着

满足下面递推的差分方程式(记为4-7-6)，确实是这样：

现在定义真实测量的数据x(n)和其近似值x^(n)之间的差为e(n)，即：

上式记为4-7-7，将式(4-7-6)代入式(4-7-7)得：

上式记为4-7-8。参数的最小二乘估计应通过使

最小化求出，但这将得到一组难于求解的非线性方程，下面不这样做。

现在定义(记为4-7-9)：

使得式4-7-8成为下式(记为4-7-10)：

现在需要使

最小，而不是让

最小。

由式4-7-10可得：

上图中的表达式对应于下面的矩阵方程：

Prony方法需要求解上面的矩阵方程。

为使：

最小，可先令上式的偏导数∂ε_p/∂a_m=0得到式(记为4-7-12)：

式中x(n-i)为样本数据；x*(n-i)为x(n-i)的共轭，下文再出现此表达式也是同样的解释，不再重复。

对应的最小误差能量为下式(记为4-7-13)：

现在定义下式(记为4-7-14)：

则式4-7-12可写作(记为4-7-15)：

求解上述方程即可得到参数a₁,a₂,…,a_p和最小误差能量ε_p的估计值。

根据刚刚得到的参数a₁,a₂,…,a_p与

可求出其根z_i，其中i=1,…,p。

现在假设z_i已被求出，则式4-7-1就简化为关于未知参数b_m的线性方程，其矩阵形式(记为4-7-16)为：

其中：

上式的矩阵Φ是N×p维的范德蒙特矩阵。因此参数向量b的解由：

给出。其中矩阵Φ^HΦ可表示为下式(记为4-7-18)：

其中，有下式(记为4-7-19)：

---

2 算法步骤

下面是Prony算法步骤：

1、构造矩阵R，其中p_e为扩展后的阶数，且有p_e>>p，p_e的最大值=N/2：

2、用SVD-TLS方法确定R的有效秩p；

3、根据式：

求解上述方程即可得到参数a₁,a₂,…,a_p和最小误差能量ε_p的估计值。

4、求式：

的根z_i，其中i=1,…,p。然后需要通过下式递推出x^(n)，其中x^(0)=x(0)：

除此之外，还有：

5、根据式4-7-16~4-7-19计算参数b₁,b₂,…,b_p；

6、利用下式计算振幅、相位、频率和衰减系数，其中k=1,2,…,p(请使用i替换下式的k)：

---

3 参考文献

1、《现代信号处理》-张贤达-清华大学出版社

2、Prony算法分析低频振荡的有效性研究_王铁强。

3、一种基于综合DFT和Prony算法的谐波与间谐波分析方法

---

END