1 ## Copyright (C) 2000 Gert Van den Eynde <na.gvandeneynde@na-net.ornl.gov>
2 ## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de>
4 ## This program is free software; you can redistribute it and/or modify it under
5 ## the terms of the GNU General Public License as published by the Free Software
6 ## Foundation; either version 3 of the License, or (at your option) any later
9 ## This program is distributed in the hope that it will be useful, but WITHOUT
10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 ## You should have received a copy of the GNU General Public License along with
15 ## this program; if not, see <http://www.gnu.org/licenses/>.
17 ## USAGE [alpha,c,rms] = expfit( deg, x1, h, y )
19 ## Prony's method for non-linear exponential fitting
21 ## Fit function: \sum_1^{deg} c(i)*exp(alpha(i)*x)
23 ## Elements of data vector y must correspond to
24 ## equidistant x-values starting at x1 with stepsize h
26 ## The method is fully compatible with complex linear
27 ## coefficients c, complex nonlinear coefficients alpha
28 ## and complex input arguments y, x1, non-zero h .
29 ## Fit-order deg must be a real positive integer.
31 ## Returns linear coefficients c, nonlinear coefficients
32 ## alpha and root mean square error rms. This method is
33 ## known to be more stable than 'brute-force' non-linear
34 ## least squares fitting.
37 ## x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
38 ## y = 2*exp(1.3*x)-0.5*exp(2*x);
39 ## error = (rand(1,length(y))-0.5)*1e-4;
40 ## [alpha,c,rms] = expfit(2,x0,step,y+error)
50 ## The fit is very sensitive to the number of data points.
51 ## It doesn't perform very well for small data sets.
52 ## Theoretically, you need at least 2*deg data points, but
53 ## if there are errors on the data, you certainly need more.
55 ## Be aware that this is a very (very,very) ill-posed problem.
56 ## By the way, this algorithm relies heavily on computing the
57 ## roots of a polynomial. I used 'roots.m', if there is
58 ## something better please use that code.
60 ## Demo for a complex fit-function:
61 ## deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
63 ## y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
64 ## A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
65 ## [alpha,c,rms]= expfit( deg, x1, h, y )
67 function [a,c,rms] = expfit(deg,x1,h,y)
70 if deg<1, error('expfit: deg must be >= 1'); end
71 if ~h, error('expfit: vanishing stepsize h'); end
73 if ( N=length( y=y(:) ) ) < 2*deg
74 error('expfit: less than %d samples',2*deg);
77 % Solve for polynomial coefficients
78 A = hankel( y(1:N-deg),y(N-deg:N) );
79 s = - A(:,1:deg) \ A(:,deg+1);
87 % nonlinear coefficients
90 % Compose second matrix A(i,j) = u(j)^(i-1)
91 A = ( ones(N,1) * u(:).' ) .^ ( [0:N-1]' * ones(1,deg) );
99 % Compute rms of y - approx
100 % where approx(i) = sum_k ( c(k) * exp(x(i)*a(k)) )
103 approx = exp( x(:) * a(:).' ) * c(:);
104 rms = sqrt( sumsq(approx - y) );
109 %!demo % same as in help - part
110 %! deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
112 %! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
113 %! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
114 %! [alpha,c,rms]= expfit( deg, x1, h, y )
116 %!demo % demo for stepsize with negative real part
117 %! deg= 2; N= 20; x1= +3*(1+i), x= linspace(x1,1+i/3,N).';
119 %! y= (2+i)*exp( (-1-2i)*x ) + (-1+3i)*exp( (2+3i)*x );
120 %! A= 5e-2; y+= A*(randn(N,1)+randn(N,1)*i); % add complex noise
121 %! [alpha,c,rms]= expfit( deg, x1, h, y )
124 %! x0 = 1.5; step = 0.05; xend = 5;
130 %! y = exp (x(:) * a(:).') * c(:);
131 %! err = randn (size (y)) * v;
132 %! plot (x, y + err);
134 %! [a_out, c_out, rms] = expfit (2, x0, step, y+err)