image thumbnail

Nonlinear fitting n-dimensional data with arbitrary functions

by

 

01 Jul 2008 (Updated )

Demonstrates the abilities of Matlab functions lsqcurvefit, fmincon and fminsearch to fit complex mo

fit_func
function y = fit_func(fun, x, xdata)
%Either calls a function pointer or uses a built-in function.
%
% fit_func is called by fit_nl and fit_nl_ex and either simply forward
% everything and call y=fun(x, xdata) or calls a built-in function
% depending on the keyword given in fun.
% Built-in are 1D-3D Gaussian/Lorentzian models with a variable number of
% Gaussian/Lorentzian peaks.
%
% Syntax
%   function y = fit_func(fun, x, xdata)
%
% Input parameters
%   fun     Either a keyword or a function pointer for a function of type
%           y = myfun(x, xdata)
%           Possible keywords are 'xd_gaussian' or 'xd_lorentzian'
%           with x=1,2,3
%   x       Parameter vector. Number of parameters for built-in functions
%           are given below.
%   xdata   The input data for the model function. For built-in functions
%           depending on the dimension this will be a struct with fields x,
%           y and z which are arrays created by ndgrid(), i.d. they contain
%           the x, y or z values of an (x,y,z) position
%
% Output parameters
%   y       Output value of function fun. Can be multi-dimensional.
%
%
% Built-in functions
% ------------------ (Read this if You want to use them!)
%
% Gaussian/Exponential function and sums of them:
%
% The basic form of one(!) 3D Gaussian peak is:
%   p1 * 2^(-(x-p2)^2/(p5/2)^2 -(y-p3)^2/(p6/2)^2 -(z-p4)^2/(p8/2)^2)
% where p1 is the amplitude (p2, p3, p4) defines the center of the peak and
% (p5, p6, p7) are the FWHMs in direction x, y and z. 7 different
% parameters are needed.
%
% For one 2D Gaussian 5 parameters are needed: p1 - amplitude, (p2, p3)
% center position and (p4, p5) FWHMs.
%   [p1 * 2^(-(x-p2)^2/(p4/2)^2 -(y-p3)^2/(p5/2)^2)]
%
% For one 1D Gaussian 3 parameters are enough: p1 - amplitude, p2 - center,
% p3 - FWHM. [p1 * 2^(-(x-p2)^2/(p3/2)^2)]
%
% The general model of the built-in functions is that we have a constant
% background plus the sum of m different gaussian peaks. Therefore for 1D
% depending on the number of Gaussians we will need 1+3m parameters,
% for 2D 1+5m and for 3D 1+7m parameters.
%
% That means that when the dimension of the gaussian is given by the
% keyword ('1/2/3d_gaussian') and the number of parameters is given by
% length(x) than the built-in function knows how many peaks are wanted.
%
% The order of the parameter is first constant background and then the
% parameters for the first gaussian peak, than next gaussian peaks and so
% on...
%
% Example: Two 1D Gaussian peaks
%   keyword is '1d_gaussian', length(x) must be 7
%   function to be computed will be (pi=x(i)=parameter):
%   p1 + p2*2^(-(x-p3)^2/(p4/2)^2)] + p5*2^(-(x-p6)^2/(p7/2)^2)]
%
%
% Lorentzian functions and sums of them:
%
% This is completely equivalent to the system above except the function
% definitions are different.
% The basic definition of a 3D Lorentzian is:
%   p1 / ((x-p2)^2/(p5/2)^2 + (y-p3)^2/(p6/2)^2 + (z-p4)^2/(p7/2)^2 + 1)
%   with amplitude p1, center (p2, p3, p4) and FWHMS (p5, p6, p7)
% And for a 2D Lorentzian:
%   p1 / ((x-p2)^2/(p4/2)^2 + (y-p3)^2/(p5/2)^2 + 1)
% And for a 1D Lorentzian:
%   p1 / ((x-p2)^2/(p3/2)^2 + 1)
%
% Example: Two 1D Lorentzian peaks
%   keyword is '1d_lorentzian', length(x) must be 7
%   function to be computed will be (pi=x(i)=parameter):
%   p1 + p2 / ((x-p3)^2/(p4/2)^2 + 1) + p5 / ((x-p6)^2/(p7/2)^2 + 1)

if nargin < 3
    error('Not enough arguments given!');
end

% check fun for keywords and select function pointers for built-in
% functions
m = [];
if ischar(fun)
    switch fun
        case '1d_gaussian'
            m.fun = @f_gaussian;
            m.dim  = 1;
        case '2d_gaussian'
            m.fun = @f_gaussian;
            m.dim  = 2;
        case '3d_gaussian'
            m.fun = @f_gaussian;
            m.dim  = 3;
        case '1d_lorentzian'
            m.fun = @f_lorentzian;
            m.dim  = 1;
        case '2d_lorentzian'
            m.fun = @f_lorentzian;
            m.dim  = 2;
        case '3d_lorentzian'
            m.fun = @f_lorentzian;
            m.dim  = 3;
        otherwise
            error('Unknown built-in function!');
    end
    % check if the number of parameters is consistent with chosen built-in
    % function
    if mod(length(x) - 1 , 2 * m.dim + 1) ~= 0
        % not all numbers of parameters are allowed
        error('Size of x is not 1+(2*d+1)m with d-dimensionality and m-number of peaks!');
    end
    % check if x_data is a struct with fields x, y, z (depending on
    % dimensionality)
    if ~isstruct(xdata)
        error('Argument xdata must be struct with fields x(,y,z) for built-in functions!');
    end
    if ~isfield(xdata, 'x')
        error('Field x missing in struct xdata!');
    end
    if m.dim > 1 & ~isfield(xdata, 'y')
        error('Field y missing in struct xdata!');
    end
    if m.dim == 3 & ~isfield(xdata, 'z');
        error('Field z missing in struct xdata!');
    end
    % check if size of x,y and z arrays is the same
    if m.dim == 2 & ~isequal(size(xdata.x), size(xdata.y))
        error('Field x and y in xdata must have same size!');
    end
    if m.dim == 3 & ~isequal(size(xdata.x), size(xdata.y), size(xdata.z))
        error('Field x, y and z in xdata must have same size!');
    end
else
    % not a built-in function just forward
    m.fun = fun;
end

% the rest is simple, just call the function
y = m.fun(x, xdata);

% builtin-functions -------------------------------------------------------

    function y = f_gaussian(x, xdata)
        % we can assume that the number of parameters in x and the
        % structure in xdata is correct because this was checked above

        % the number of different gaussians in the sum is n
        n = (length(x) - 1) / (2 * m.dim + 1);

        % set background
        y = zeros(size(xdata.x)) + x(1);

        % for all peaks
        for ki = 1 : n
            % get the right parameters for this peak out of x
            si = 2 + (ki - 1) * (2 * m.dim + 1);
            p = x(si:si+2*m.dim);
            % depending on dimensionality add a peak
            switch m.dim
                case 1
                    y = y + p(1) * power(2., -(xdata.x - p(2)).^2 / (p(3) / 2.)^2);
                case 2
                    y = y + p(1) * power(2., -(xdata.x - p(2)).^2 / (p(4) / 2.)^2 -(xdata.y - p(3)).^2 / (p(5) / 2.)^2);
                case 3
                    y = y + p(1) * power(2., -(xdata.x - p(2)).^2 / (p(5) / 2.)^2 -(xdata.y - p(3)).^2 / (p(6) / 2.)^2 -(xdata.z - p(4)).^2 / (p(7) / 2.)^2);
            end % no error handling necessary m.dim was set by us
        end
    end

    function y = f_lorentzian(x, xdata)
        % we can assume that the number of parameters in x and the
        % structure in xdata is correct because this was checked above

        % the number of different lorentzians in the sum is n
        n = (length(x) - 1) / (2 * m.dim + 1);

        % set background
        y = zeros(size(xdata.x)) + x(1);

        % for all peaks
        for ki = 1 : n
            % get the right parameters for this peak out of x
            si = 2 + (ki - 1) * (2 * m.dim + 1);
            p = x(si:si+2*m.dim);
            % depending on dimensionality add a peak
            switch m.dim
                case 1
                    y = y + p(1) ./ ((xdata.x - p(2)).^2 / (p(3) / 2.)^2 + 1);
                case 2
                    y = y + p(1) ./ ((xdata.x - p(2)).^2 / (p(4) / 2.)^2 + (xdata.y - p(3)).^2 / (p(5) / 2.)^2 + 1);
                case 3
                    y = y + p(1) ./ ((xdata.x - p(2)).^2 / (p(5) / 2.)^2 + (xdata.y - p(3)).^2 / (p(6) / 2.)^2 + (xdata.z - p(4)).^2 / (p(7) / 2.)^2 + 1);
            end % no error handling necessary m.dim was set by us
        end
    end

end
Error using ==> fit_func at 84
Not enough arguments given!

Error in ==> fit_func at 84
    error('Not enough arguments given!');

Contact us