From help-octave-request at bevo dot che dot wisc dot edu Thu Jan 24 01:43:31 2002
Subject: Re: leasqr
From: Gert Van den Eynde <gvdeynde at sckcen dot be>
To: Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov>
Cc: "Christian T. Steigies" <cts at debian dot org>, help-octave@bevo.che.wisc.edu
Date: 24 Jan 2002 08:42:10 +0100


--=-ENQisBE2gKVpFObRt+LX
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

Hi Christian, Paul,


Another possible algorithm is Prony's algorithm. It's designed to fit
sums of exponential functions (c(i)*exp(a(i)*x)) but it also works for
imaginary a(i), and thus for sine/cosine pairs. One drawback is the
requirement that the data points need to be equidistant. My professor at
university always claimed this algorithm was more stable than the least
squares approach.... Some time ago, I implemented this for Octave but
never got around to post it to the mailing-list. So now, here it is...
Usual disclaimers applied. Hope it works. Comments and improvements
welcome.



Have a nice day,

Gert



*-*-*-*
Gert Van den Eynde
Reactor Physics & MYRRHA Dept.
SCK-CEN
Belgium
*-*-*-*

--=-ENQisBE2gKVpFObRt+LX
Content-Disposition: attachment; filename=prony.m
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain; charset=ISO-8859-1

## usage [alpha,c,rms] =3D 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 =3D 0; step =3D 0.05; xend =3D 5; x =3D x0:step:xend;
##    y =3D 2*exp(1.3*x)-0.5*exp(2*x);
##
##    error =3D (rand(1,length(y))-0.5)*1e-4;
##
##    [alpha,c,rms] =3D prony(2,x0,step,y+error)
##
##  alpha =3D
##
##    2.0000
##    1.3000
##
##  c =3D
##
##    -0.50000
##     2.00000
##
##    rms =3D 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
## na dot gvandeneynde at na-net dot ornl dot gov
##
## 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] =3D prony(deg,x1,h,y)


% Check input

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

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

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


   N =3D length(y);

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

% Solve for polynomial coeff

%  Compose matrix

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

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

%  Least squares solution
   s =3D A1\b1;

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

%  Compute roots
   u =3D roots(p);


% Compose second matrix

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

% Solve linear system

    f =3D A2\b2;

% Decompose all

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


% Compute rms

    rms =3D 0.0;

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

    rms =3D sqrt(rms);

endfunction

--=-ENQisBE2gKVpFObRt+LX--



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------