Using Prony's method to fit a sum of exponentials

Submitted by
Tillsten on 10 November 2011
Download
File
Update history
Revision 4 of 4: previous
Updated by
joseA on 17 September 2013
Tags
# License: Creative Commons Zero (almost public domain) http://scpyce.org/cc0 import scipy.linalg as la import numpy.linalg as nla import numpy as np def expfit(deg,x1,h,y): #Build matrix A=la.hankel(y[:-deg],y[-deg-1:]) a=-A[:,:deg] b= A[:,deg] #Solve it s=nla.lstsq(a,b)[0] #Solve polynomial p=np.flipud(np.hstack((s,1))) u=np.roots(p) #Calc exponential factors a=np.log(u)/h #Build power matrix A=np.power((np.ones((len(y),1))*u[:,None].T),np.arange(len(y))[:,None]*np.ones((1,deg))) #solve it f=nla.lstsq(A,y)[0] #calc amplitudes c=f/np.exp(a*x1) #build x, approx and calc rms x=x1+h*np.arange(len(y)) approx=0*x for ii in range(len(a)): approx += np.exp(x*a[ii]).dot(c[ii]) rms=np.sqrt(((approx-y)**2).sum()/len(y)) return a, c,rms if __name__=='__main__': import matplotlib.pyplot as plt x0 = 0 step = 0.005 xend = 3 deg = 2 a = [2, -0.5, 0.8] # b = [1.3+0.5*1j, 2.0, 1] b = [1.3, 2.0, 1] x = np.arange(x0,xend+step,step) y = 0.0 for ii in range(len(a)): y += a[ii]*np.exp(b[ii]*x) error = (np.random.rand(len(y))-0.5)*1e-4 beta,alpha,rms = expfit(deg,x0,step,y+y*error) prony = 0*x for ii in range(len(alpha)): prony += np.exp(x*beta[ii]).dot(alpha[ii]) plt.plot(x,y, linewidth=2.0,label='Function') plt.plot(x,prony,'r--', linewidth=2.0,label='Prony') plt.legend() plt.show() print 'alfa = ' + str(alpha) print 'beta = ' + str(beta) print 'RMS = ' + str(rms)

This a basic implementation of Prony’s method. In this form, it is very susceptible to noise.

Added by Jose:

I had some problems using “lstsq” from scipy.linalg that were solved when I used instead the version from numpy.linalg.

I have also changed the way the RMS is computed.

Please Sign in or Register to leave a comment

• Creative Commons Zero. No rights reserved.
Users have permission to do anything with the code and other material on this page. (More details)
• X 17
• X 16
• X 12
• X 9
• X 8
• X 7
• X 7
• X 6
• X 5
• X 4
• X 4
• X 4
• X 3
• X 3
• X 3
• Trademarks are property of their respective owners. Code and comments are owned by their respective posters. © 2013 All Rights Reserved