Prony算法步骤及示例代码

---

0 前言

下文描述的Prony算法中，N为数据的个数，形式为(t₀,x_t₀),(t₁,x_t₁),…,(t_N-1,x_t_N-1)，另外x_t_i=x(i)，i=0,1,…,N-1。

---

1 算法步骤

下面是Prony算法的步骤：

1、构造矩阵R，其中p_e为扩展后的阶数，且有p_e>>p，p为Prony算法的真实阶数，一般取p_e=N/2：

其中：

上式中的x*(n-i)是x(n-i)的共轭，如果使用的数据中无虚数而是实数的话就有x(n-i)=x*(n-i)。

2、下面使用SVD分解来确定矩阵R的有效秩p：首先，R=U∑V^H，Σ矩阵的主对角线上的元素σ_ii是矩阵R的奇异值，并按此顺序排列：σ₁₁≥σ₂₂≥…≥σ_hh≥0。考虑下面的归一化比值：

上图中的h为矩阵R的奇异值总个数，通常选择一个接近于1的数作为门限值(例如0.995)，并把v(k)大于该门限值的最小k值定为矩阵R的有效秩p。

3、定义(p+1)×(p+1)维矩阵S^(p)：

上式中的σ_ii是矩阵R的奇异值，v_jⁱ定义为：

v(i,j)是对矩阵R进行SVD分解后得到的矩阵V的第i行第j列上的因素。根据

求解a=(a₁,a₂,…,a_p)的值，S^(-p)为S^(p)的逆矩阵，而S^(-p)(i,j)指其第i行第j列上的因素。

4、求式：

的根z_i，其中i=1,…,p。然后需要通过

与

两个式子递推出x^(n)，这里n=0,1,…,N-1。

5、看下图：

上图的矩阵Φ是N×p维的范德蒙特矩阵。其中：

于是参数向量b=(b₁,b₂,…,b_p)由

给出。

6、最后利用下图计算振幅A、相位φ、频率f和衰减系数a(又称为吸收系数或衰减因子)，其中k=1,2,…,p：

其中|x|是复数x的模长，Im(x)为复数x的虚部，Re(x)为复数x的实部，Δt是采样时间的间隔。

---

2 示例代码

希望按我的思路浏览下文。

首先，下面是模拟信号的表达式及其图像：

下面的代码框对上述信号进行采样，采样点数为N=500，采样间隔Δt=0.001s：

import numpy as np

t=np.arange(0,0.5,0.001)
x_t=[]
for i in range(len(t)):
    x=t[i]
    y=200*np.exp(-6*x)*np.sin(2*np.pi*60*x+np.pi/4)+50*np.sin(2*np.pi*40*x+np.pi/5)
    x_t.append(y)

下面是显示500个离散采样点对应的信号的代码及图像：

from matplotlib import pyplot as plt

plt.plot(t,x_t)
plt.show()

下面是Prony算法主程序：

import numpy as np
import math
from matplotlib import pyplot as plt

N=len(x_t) #样本数量
R=[] #初始化样本函数矩阵R为空列表
pe=int(np.floor(N/2)) #矩阵R的阶数为(1+pe)×(1+pe)

for i in range(0,pe+1): #i的范围为[0,pe]
    temp_list=[]
    for j in range(0,pe+1): #j的范围为[0,pe]
        temp_sum=0
        for n in range(pe,N): #n的范围为[pe,N-1]
            temp_sum+=x_t[n-j]*x_t[n-i]
        temp_list.append(temp_sum)
    R.append(temp_list)

R=np.array(R) #将列表转换为数组形式

u,sigma,vh=np.linalg.svd(R) #对矩阵R进行SVD分解，R=u·sigma·vh
#上面的sigma是一维数组，存储的元素为非零奇异值
#直接进行u·sigma·vh是不能运算的，需要将sigma变换为相应阶的矩阵
V=vh.T
limit=0.995 #门限值limit

sums=0
for sigm in sigma:
    sums+=sigm*sigm
    
temp=0
for i in range(len(sigma)):
    temp+=sigma[i]*sigma[i]
    if np.sqrt(temp/sums)>limit:
        break

p=i+1 #最终确定的阶数

S_p=[[0 for i in range(p+1)] for i in range(p+1)]

for j in range(1,p+1):
    s_sigma=sigma[j-1]
    for i in range(1,pe-p+1+1):
        v_i_j=[]
        for x in range(p+1):
            v_i_j.append([V[i+x-1][j-1]])
        v_i_j=np.array(v_i_j)
        temp=s_sigma*s_sigma*(np.matmul(v_i_j,v_i_j.T))
        S_p+=temp

S_p_inv=np.linalg.inv(S_p)
a_i=[]
for i in range(1,p+1):
    a_i.append(S_p_inv[i][0]/S_p_inv[0][0])
    
a=np.array(a_i)
a_i.insert(0,1)
z=np.roots(a_i)

x_t_hat=[]
for i in range(p):
    x_t_hat.append(x_t[i])
for n in range(p,N):
    temp_sum=0
    for m in range(1,p+1):
        temp_sum-=a[m-1]*x_t_hat[n-m]
    x_t_hat.append(temp_sum)
x_t_hat=np.array(x_t_hat)
x_t_hat=x_t_hat.reshape(N,1)

phi_T=[]
for z_i in z:
    temp=[]
    for i in range(N):
        temp.append(pow(z_i,i))
    phi_T.append(temp)
phi_T=np.array(phi_T)
phi_H=phi_T.conjugate()
phi=phi_T.T

phi_H_phi=np.matmul(phi_H,phi)
temp1=np.linalg.inv(phi_H_phi)
temp2=np.matmul(temp1,phi_H)
b=np.matmul(temp2,x_t_hat)

#最后一步，求四个参数
A_i=[] #振幅
phi_i=[] #相位
f_i=[] #频率
alpha_i=[] #吸收系数

for i in range(p):
    A_i.append(abs(b[i]))
    phi_i.append(np.arctan(b[i].imag/b[i].real))
    f_i.append(np.arctan(z[i].imag/z[i].real)/(2*np.pi*(t[1]-t[0])))
    alpha_i.append(np.log(abs(z[i]))/(t[1]-t[0]))

print('最终计算得到的【'+str(p)+'】组余弦分量的四个参数分别如下：\n振幅\t频率\t相位\t吸收系数')

for i in range(p):
    print(str(A_i[i])+'\t'+str(f_i[i])+'\t'+str(phi_i[i])+'\t'+str(alpha_i[i]))

输出为：

最终计算得到的【4】组余弦分量的四个参数分别如下：
振幅	频率	相位	吸收系数
[99.99999996]	60.00000000207298	[-0.78539816]	-6.0000000656132055
[99.99999996]	-60.00000000207298	[0.78539816]	-6.0000000656132055
[25.00000002]	40.00000000784647	[-0.94247779]	-3.761068124316161e-07
[25.00000002]	-40.00000000784647	[0.94247779]	-3.761068124316161e-07

下面是最后一个代码块，将原信号(蓝色)与通过Prony算法得到的信号(红色)放在同一个坐标系中进行对比：

import numpy as np
from matplotlib import pyplot as plt

xp_t=0
for i in range(len(A_i)):
    A_k,xita_k,f_k,a_k=A_i[i],phi_i[i],f_i[i],alpha_i[i]
    xp_t+=A_k*np.exp(a_k*t)*np.cos(2*np.pi*f_k*t+xita_k)

plt.title('The Graph of the Signal')
plt.xlabel('Times(s)')
plt.ylabel('Amplitude')
plt.plot(t+0.002,x_t,c='b') #原信号，蓝色，为加以区别，特将原信号信号水平平移一段距离
plt.plot(t,xp_t,c='r') #Prony算法得到的信号，红色
plt.show()

肉眼可见，图像完全重合。

---

3 算法评价

本文的算法步骤及代码对有噪声的数据可能效果不好，还需要改进。

---

END