Skip to content

Instantly share code, notes, and snippets.

@shane5ul
Last active November 27, 2024 23:51
Show Gist options
  • Save shane5ul/8b360bb605baa9e29e5e6ede364f4d7d to your computer and use it in GitHub Desktop.
Save shane5ul/8b360bb605baa9e29e5e6ede364f4d7d to your computer and use it in GitHub Desktop.
This program uses the Prony method to fit data to F(t) = \sum_{i=1}^{m} a_i e^{b_i t}.
def prony(t, F, m):
"""Input : real arrays t, F of the same size (ti, Fi)
: integer m - the number of modes in the exponential fit
Output : arrays a and b such that F(t) ~ sum ai exp(bi*t)"""
import numpy as np
import numpy.polynomial.polynomial as poly
# Solve LLS problem in step 1
# Amat is (N-m)*m and bmat is N-m*1
N = len(t)
Amat = np.zeros((N-m, m))
bmat = F[m:N]
for jcol in range(m):
Amat[:, jcol] = F[m-jcol-1:N-1-jcol]
sol = np.linalg.lstsq(Amat, bmat)
d = sol[0]
# Solve the roots of the polynomial in step 2
# first, form the polynomial coefficients
c = np.zeros(m+1)
c[m] = 1.
for i in range(1,m+1):
c[m-i] = -d[i-1]
u = poly.polyroots(c)
b_est = np.log(u)/(t[1] - t[0])
# Set up LLS problem to find the "a"s in step 3
Amat = np.zeros((N, m))
bmat = F
for irow in range(N):
Amat[irow, :] = u**irow
sol = np.linalg.lstsq(Amat, bmat)
a_est = sol[0]
return a_est, b_est
@rizky-r-lab
Copy link

Hi @shane5ul ,
Thank you for providing simple prony method code. I am using it for analysis for power system small signal analysis.
Thanks

@Vhvj971
Copy link

Vhvj971 commented Nov 27, 2024

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment