fitting sum of two decaying exponentials (over damped oscillator)

by

 

five-parameter fit to the equation A*exp(-t/T1) + B*exp( -t/T2 ) + offset

fit_overdamping( varargin )
function varargout = fit_overdamping( varargin )
%[A, B, T1, T2, offset, Ssq] = fit_overdamping( t, y, options )
%fitting overdamped harmonic oscillator of the form
%y(t) = A exp( -t/T1 ) + B exp( -t/T2) + offset
%where A and B are the amplitudes (assumed positive), and T1, T2 are the
%decay times. Ssq is the sum of residuals.
%
%options are 'notify' to show progress and 'plot' to give a plot
%you can also supply a guess or series of guesses via a row vector or
%matrix, m, where size( m, 2 ) = 5;
%
%if you just run program without arguments it will do an example
%
%examples are
%fit_overdamping( t, y, 'notify', 'plot')
%[A, B, T1, T2, offset] = fit_overdamping( t, y, [1 2 3 4 5])
%
%Note: this script uses LMFnlsq from MATLAB central file exchange
%sak@wpi.edu 2/8/12

fnc = @(A, B, T1, T2, offset, t ) A*exp(-t/T1) + B*exp( -t/T2 ) + offset;

if nargin == 0
    disp( 'doing sample fitting' );
    t = linspace( -2, 3, 100 );
    [A, B, T1, T2, offset] = deal( 1, 1, 10, .1, 1 );    
    y = fnc( A, B, T1, T2, offset, t )+randn(size( t ) )*.01;
    varargin = { t, y, 'plot', 'Xnotify' };
end;
[t, y] = deal( shiftdim( varargin{1} ), shiftdim( varargin{2} ) );
varargin(1:2) = [];

plot_flag = strcmpi( varargin, 'plot' );
varargin( plot_flag ) = [];

if any( strcmpi( varargin, 'notify' ) ); 
    notify = @(str) fprintf( str );
else
    notify = @(str) '';
end
varargin( strcmpi( varargin, 'notify' ) ) = [];

if any( plot_flag )
    figure( sum( mfilename ) );
    clf;
    semilogy( t, y, '.r' );
end

t_ = t - t(1);
if any( cellfun( 'isclass', varargin, 'double' ) )
    guess = varargin{cellfun( 'isclass', varargin, 'double' )};    
else
    [A, T1, T2] = ndgrid( linspace( 0, y(1)-y(end), 6), t_(end)*10.^[-3:1], t_(end)*10.^[-3:1] );
    guess = [reshape( A, [], 1), reshape( y(1)-y(end)-A, [], 1), reshape( T1, [], 1), reshape( T2, [], 1), repmat( y(end), numel( A ), 1 ) ];
end

J = @(param) [ exp(-t_/param(3)), exp(-t_/param(4)), param(1)/param(3)^2*t_.*exp(-t_/param(3)), param(2)/param(4)^2*t_.*exp(-t_/param(4)), ...
    ones( size( t_) )];
warning off;
param = nan( size( guess ) );
[Ssq, CNT] = deal( nan( size( guess, 1), 1) );
tic;
for i=1:size(guess, 1);
    notify( sprintf( 'trying param A = %g, T1 = %g, T2 = %g, phase = %g, offset = %g ', guess(i,:) ) );        
    [param(i,:), Ssq(i), CNT(i)] = LMFnlsq( @(param) fnc( param(1), param(2), param(3), param(4), param(5), t_ ) - y, guess(i,:), 'Jacobian', J, 'MaxIter', 20 );
    notify( sprintf( 'error %g \n', Ssq(i) ) );
end
notify(sprintf( '%.2f sec \n', toc ) );
warning on;
[best, i] = min( Ssq );
%refining solution
[param(i,:), Ssq(i), CNT(i), Res, XY] = LMFnlsq( @(param) fnc( param(1), param(2), param(3), param(4), param(5), t_ ) - y, param(i, :), 'Jacobian', J );
A = param(i, 1)*exp(t(1)/param(i, 3));
T1 = param(i,3);
B = param(i, 2)*exp(t(1)/param(i, 4));
T2 = param(i,4);
offset = param(i, 5 );
if any( plot_flag)
    hold on;
    plot( t, fnc(A, B, T1, T2, offset, t ), '-k' );
    title( sprintf( 'A = %g, T1 = %g, B = %g, T2 = %g, offset = %g ', A, T1, B, T2, offset ) );
end

varargout = {A, B, T1, T2, offset, Ssq};
varargout = varargout(1:nargout);

Contact us