## usage [alpha,c,rms] = prony(deg,x1,h,y)
##
## Implementation of Prony's method to perform
## non-linear# Implementation of Prony's method to perform
## non-linear exponential fitting on equidistant
## data.
##
## Fit function: \sum_1^{deg} c(i)*exp(alpha(i)*x)
##
## Data is equidistant, starts at x1 with a stepsize h
##
## Function returns linear coefficients c, non-linear
## coefficients alpha and root mean square error rms
## This method is known to be more stable than 'brute-
## force' non-linear least squares fitting.
##
## The non-linear coefficients alpha can be imaginary,
## fitting function becomes a sum of sines and cosines.
##
## Example
##    x0 = 0; step = 0.05; xend = 5; x = x0:step:xend;
##    y = 2*exp(1.3*x)-0.5*exp(2*x);
##
##    error = (rand(1,length(y))-0.5)*1e-4;
##
##    [alpha,c,rms] = prony(2,x0,step,y+error)
##
##  alpha =
##
##    2.0000
##    1.3000
##
##  c =
##
##    -0.50000
##     2.00000
##
##    rms = 0.00028461
##
## The fitting is very sensitive to the number of data points.
## It doesn't perform very well for small data sets (theoretically,
## you need at least 2*deg data points, but if there are errors on
## the data, you certainly need more.
##
## Be aware that this is a very (very,very) ill-posed problem.
## By the way, this algorithm relies heavily on computing the roots
## of a polynomial. I used 'roots.m', if there is something better
## please use that code.
##
## Copyright (C) 2000: Gert Van den Eynde
## SCK-CEN (Nuclear Energy Research Centre)
## Boeretang 200
## 2400 Mol
## Belgium
## address@hidden
##
## This code is under the Gnu Public License (GPL) version 2 or later.
## I hope that it is useful, but it is WITHOUT ANY WARRANTY;
## without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE.


function [alpha,c,rms] = prony(deg,x1,h,y)


% Check input

if (deg < 1) error ('deg must be >= 1');end

if (h <= 0) error ('Stepsize h must be > 0');end;

%if (length(y) <= 2*deg)
% error ('Number of data points must be larger than 2*deg');end;


   N = length(y);

   alpha = zeros(deg,1);
   c = zeros(deg,1);
   rms = 0.0;

% Solve for polynomial coeff

%  Compose matrix

   A1 = zeros(N-deg,deg);
   b1 = zeros(N-deg,1);

   for k=1:N-deg,
     for l=1:deg
       A1(k,l) = y(k+l-1);
     end
     b1(k,1) = -y(k+deg);
   end

%  Least squares solution
   s = A1\b1;

%  Compose polynomial
   p = [1,fliplr(s')];

%  Compute roots
   u = roots(p);


% Compose second matrix

   for k=1:N,
     for l=1:deg,
       A2(k,l) = u(l)^(k-1);
     end
     b2(k) = y(k);
   end

% Solve linear system

    f = A2\b2;

% Decompose all

    for l=1:deg,
      alpha(l) = log(u(l))/h;
      c(l) = f(l)/exp(alpha(l)*x1);
    end


% Compute rms

    rms = 0.0;

    for k=1:N,
      approx = 0.0;
        for l=1:deg,
          approx = approx + c(l)*exp(alpha(l)*(x1+h*(k-1)));
        end
      rms = rms + (approx - y(k))^2;
    end

    rms = sqrt(rms);

endfunction