octave_packages/optim-1.2.0/expfit.m

   1 ## Copyright (C) 2000 Gert Van den Eynde <na.gvandeneynde@na-net.ornl.gov>
   2 ## Copyright (C) 2002 Rolf Fabian <fabian@tu-cottbus.de>
   3 ##
   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
   7 ## version.
   8 ##
   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
  12 ## details.
  13 ##
  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/>.
  16
  17 ## USAGE  [alpha,c,rms] = expfit( deg, x1, h, y )
  18 ##
  19 ## Prony's method for non-linear exponential fitting
  20 ##
  21 ## Fit function:   \sum_1^{deg} c(i)*exp(alpha(i)*x)
  22 ##
  23 ## Elements of data vector y must correspond to
  24 ## equidistant x-values starting at x1 with stepsize h
  25 ##
  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.
  30 ##
  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.
  35 ##
  36 ## Example
  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)
  41 ##
  42 ##  alpha =
  43 ##    2.0000
  44 ##    1.3000
  45 ##  c =
  46 ##    -0.50000
  47 ##     2.00000
  48 ##  rms = 0.00028461
  49 ##
  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.
  54 ##
  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.
  59 ##
  60 ## Demo for a complex fit-function:
  61 ## deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
  62 ## h = x(2) - x(1)
  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 )
  66
  67 function [a,c,rms] = expfit(deg,x1,h,y)
  68
  69 % Check input
  70 if deg<1, error('expfit: deg must be >= 1');     end
  71 if  ~h,   error('expfit: vanishing stepsize h'); end
  72
  73 if ( N=length( y=y(:) ) ) < 2*deg
  74    error('expfit: less than %d samples',2*deg);
  75 end
  76
  77 % Solve for polynomial coefficients
  78 A = hankel( y(1:N-deg),y(N-deg:N) );
  79 s = - A(:,1:deg) \ A(:,deg+1);
  80
  81 % Compose polynomial
  82 p = flipud([s;1]);
  83
  84 % Compute roots
  85 u = roots(p);
  86
  87 % nonlinear coefficients
  88 a = log(u)/h;
  89
  90 % Compose second matrix A(i,j) = u(j)^(i-1)
  91 A = ( ones(N,1) * u(:).' ) .^ ( [0:N-1]' * ones(1,deg) );
  92
  93 % Solve linear system
  94 f = A\y;
  95
  96 % linear coefficients
  97 c = f./exp( a*x1 );
  98
  99 % Compute rms of y - approx
 100 % where approx(i) = sum_k ( c(k) * exp(x(i)*a(k)) )
 101 if nargout > 2
 102    x = x1+h*[0:N-1];
 103    approx = exp( x(:) * a(:).' ) * c(:);
 104    rms = sqrt( sumsq(approx - y) );
 105 end
 106
 107 endfunction
 108
 109 %!demo   % same as in help - part
 110 %! deg= 2; N= 20; x1= -(1+i), x= linspace(x1,1+i/2,N).';
 111 %! h = x(2) - x(1)
 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 )
 115
 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).';
 118 %! h = x(2) - x(1)
 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 )
 122
 123 %!demo
 124 %! x0 = 1.5; step = 0.05; xend = 5;
 125 %! a = [1.3, 2]';
 126 %! c = [2, -0.5]';
 127 %! v = 1e-4;
 128 %!
 129 %! x = x0:step:xend;
 130 %! y = exp (x(:) * a(:).') * c(:);
 131 %! err = randn (size (y)) * v;
 132 %! plot (x, y + err);
 133 %!
 134 %! [a_out, c_out, rms] = expfit (2, x0, step, y+err)
 135
 136