Monday, September 25, 2017

Prony Method

Given N equispaced data-points \(F_i = F(t = i \Delta t)\), where \(i = 0, 1, ..., N-1\), the Prony method can be used to fit a sum of m decaying exponenitals: \[F(t) = \sum_{i=1}^{m} a_i e^{b_i t}. \] The 2m unknowns are \(a_i\) and \(b_i\).

In the Prony method, the number of modes in the exponential (m) is pre-specified. There are other methods, which are more general.

Here is a python subprogram which implements the Prony method.
If you have arrays t and F, it can be called as:

a_est, b_est = prony(t, F, m)

4 comments:

Wave energy ridder said...

Dear Clueless Fundatma,
Can you please post more information about the code or give a reference?
The following line seems to be missing in the code
N = len(t)
probably after the beginning of the function Prony()
Kind regards,
João Henriques

Sachin Shanbhag said...

Thanks for pointing that error out, Joao. There is a link to the method in the writeup above. I just implemented it.

Unknown said...

Dear Clueless Fundatma,
The method produces an error when you are computing a_est and b_est for a signal like:
y[k]=math.exp(-3.0*t[k])+math.exp(-2.0*t[k])*math.cos(50.0*2.0*math.pi*t[k]);

the correction is to assume a complex matrix for the least square solver.

I modified the line:

# Set up LLS problem to find the "a"s in step 3
Amat = np.zeros((N, m),dtype=complex)

This solves the problem when you need to compute the oscillation frequency of the sine or cosine waves.

Sachin Shanbhag said...

Thanks for that Santiago.