Search code examples
matlabinterpolationhandlecurve-fittingpiecewise

Storing and accessing a piecewise interpolant using a single handle

I have a dataset containing two vectors of points, X and Y that represents measurements of an "exponential-like" phenomenon (i.e. Y = A*exp(b*x)). When fitting it with an exponential equation I'm getting a nice-looking fit, but when I'm using it to compute things it turns out that the fit is not quite as accurate as I would hope.

Currently my most promising idea is a piecewise exponential fit (taking about 6 (x,y) pairs each time), which seems to provide better results in cases I tested manually. Here's a diagram to illustrate what I mean to do:

// Assuming a window of size WS=4:
 - - - - - - - - - - - -  //the entire X span of the data
[- - - -]- - // the fit that should be evaluated for X(1)<= x < X(floor(WS/2)+1)
 -[- - - -]- // the fit that should be evaluated for X(3)<= x < X(4)
...
 - - - - - -[- - - -]- - // X(8)<= x < X(9)
... //and so on

Some Considerations:

  • I considered filtering the data before fitting, but this is tricky since I don't really know anything about the type of noise it contains.
  • I would like the piecewise fit (including all different cases) to be accessible using a single function handle. It's very similar to MATLAB's Shape-preserving "PCHIP" interpolant, only that it should use an exponential function instead.
  • The creation of the fit does not need to happen during the runtime of another code. I even prefer to create it separately.
  • I'm not worried about the potential unsmoothness of the final function.
  • The only way of going about this I could think of is defining an .m file as explained in Fit a Curve Defined by a File, but that would require manually writing conditions for almost as many cases as I have points (or obviously write a code that generates this code for me, which is also a considerable effort).

Relevant code snippets:

clear variables; clc;
%% // Definitions
CONSTS.N_PARAMETERS_IN_MODEL = 2; %// For the model y = A*exp(B*x);
CONSTS.WINDOW_SIZE = 4; 
assert(~mod(CONSTS.WINDOW_SIZE,2) && CONSTS.WINDOW_SIZE>0,...
    'WINDOW_SIZE should be a natural even number.');
%% // Example Data
X = [0.002723682418630,0.002687088539567,0.002634005004610,0.002582978173834,...
     0.002530684550171,0.002462144527884,0.002397219225698,0.002341097974950,...
     0.002287544321171,0.002238889510803]';
Y = [0.178923076435990,0.213320004768074,0.263918364216839,0.324208349386613,...
     0.394340431220509,0.511466688684182,0.671285738221314,0.843849959919278,...
     1.070756756433334,1.292800046096531]';
assert(CONSTS.WINDOW_SIZE <= length(X),...
    'WINDOW_SIZE cannot be larger than the amount of points.');
X = flipud(X); Y = flipud(Y); % ascending sorting is needed later for histc
%% // Initialization
nFits = length(X)-CONSTS.WINDOW_SIZE+1;
coeffMat(nFits,CONSTS.N_PARAMETERS_IN_MODEL) = 0; %// Preallocation
ft = fittype( 'exp1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%% // Populate coefficient matrix
for ind1 = 1:nFits
    [xData, yData] = prepareCurveData(...
        X(ind1:ind1+CONSTS.WINDOW_SIZE-1),Y(ind1:ind1+CONSTS.WINDOW_SIZE-1));
    %// Fit model to data:
    fitresult = fit( xData, yData, ft, opts );
    %// Save coefficients:
    coeffMat(ind1,:) = coeffvalues(fitresult);
end
clear ft opts ind1 xData yData fitresult ans
%% // Get the transition points
xTrans = [-inf; X(CONSTS.WINDOW_SIZE/2+1:end-CONSTS.WINDOW_SIZE/2); inf];

At this point, xTrans and coeffMat contain all the required information to evaluate the fits. To show this we can look at a vector of some test data: 

testPts = [X(1), X(1)/2, mean(X(4:5)), X(CONSTS.WINDOW_SIZE)*1.01,2*X(end)];

...and finally the following should roughly happen internally within the handle:

%% // Get the correct fit# to be evaluated:
if ~isempty(xTrans) %// The case of nFits==1
    rightFit = find(histc(testPts(3),xTrans));
else
    rightFit = 1;
end
%% // Actually evaluate the right fit
f = @(x,A,B)A*exp(B*x);
yy = f(testPts(3),coeffMat(rightFit,1),coeffMat(rightFit,2));

And so my problem is how to hold that last bit (along with all the fit coefficients) inside a single handle, that accepts an arbitrarily-sized input of points to interpolate on?

Related resources:

Solution

  • It's not all clear but why not to puts things into a class ?

    classdef Piecewise < handle

    methods

        % Construction
        function [this] = Piecewise(xmeas, ymeas)
            ... here compute xTrans and coeffMat...
        end

        % Interpolation
        function [yinterp] = Evaluate(xinterp)
            ... Here use previously computed xTrans and coeffMat ...
        end
    end

    properties(SetAccess=Private, GetAccess=Private)
        xTrans;
        coeffMat;
    end

end

    In this way you can prcompute xTrans vector and coeffMat matrix during construction and then later reuse these properties when you need to evaluate interpolant at xinterp values in Evaluate method.

    % The real measured data
xmeas = ...
ymeas = ...

% Build piecewise interpolation object
piecewise = Piecewise(x,meas, ymeas);

% Rebuild curve at any new positions
xinterp = ...
yinterp = piecewise.Evaluate(xinterp);

    Function like syntax

    If you truly prefer to have function-handle like syntax, you can still internally use above Piecewise object and add extra static method to return it as a function handle: 

    classdef Piecewise < handle

    ... see code above...

    methods(Static=true)
        function [f] = MakeHandle(xmeas, ymeas)
        %[
            obj = Piecewise(xmeas, ymeas);
            f = @(x)obj.Evaluate(x);
        %]
    end

end

    This can be used like this:

    f = Piecewise.MakeHandle(xmeas, ymeas);
yinterp = f(xinterp);

    PS1: You can later put Evaluate and Piecewise constructor methods as private if you absolutely wanna to force this syntax.

    PS2: To fully hide object oriented design, you can turn MakeHandle into a fully classic routine (will work the same as if static and users won't have to type Piecewise. in front of MakeHandle).

    A last solution without oop

    Put everything in a single .m file : 

    function [f] = MakeHandle(xmeas, ymeas)
    ... Here compute xTrans and coeffMat ...
    f = @(x)compute(x, xTrans, coeffMat);% Passing xTrans/coeffMatt as hidden parameters
end
function [yinterp] = compute(xinterp, xTrans, coeffMat)
    ... do interpolation here...
 end