clear, clc

    % Orders must be odd
    order = 3;

    % *******************************************************************
    % Set up signals
    % *******************************************************************

    fun = inline('sin(2*pi*x*3/40) + sin(2*pi*x*5/40)', 'x');

    % Input points and samples
    nIn = 401;
    xIn = 1:nIn;
    inData = fun(xIn);

    % Output points and reference samples
    ratio = 53/40;
    xOut = 1 : 1/ratio : nIn;
    interpolatorDelay = order/2-1;
    ref = fun(xOut - interpolatorDelay);

    nOut = numel(xOut);
    outData = zeros(1, nOut);

    % *******************************************************************
    % Apply Lagrange interpolation to input data
    % *******************************************************************

    % Lagrange polynomial coefficients for each tap
    for ixTap = 0:order
        c(:,ixTap+1) = lagrangeBasisPolynomial(order, ixTap);
    end

    % Calculate samples at the new points xOut
    outData = [];
    i=1;
    for point = xOut

        % Calculate the fractional delay of 'point'
        xOutFrac = point - floor(point);

        outTemp = 0;
        for ixTap = 0:order
            % inSample is the input sample that appears at FIR tap 'ixTap'
            % to contribute to the output sample

            % Use end samples for the out-of-range samples
            % inSample = mod(floor(point) + ixTap - order, nIn) + 1;

            % Just skip the out-of-range samples
            inSample = floor(point) + ixTap - order + 1;

            if inSample < 1 || inSample > nIn
                continue;
            end

            % The value of said input sample
            d = inData( inSample );

            % Evaluate the Lagrange polynomial, resulting in the
            % time-varying FIR tap weight
            cTap = polyval(c(:,ixTap+1), xOutFrac);

            % FIR operation: multiply FIR tap weight with input sample and
            % add to output sample
            outTemp = outTemp + d.*cTap;
        end

        outData(end+1) = outTemp;
    end

    % *******************************************************************
    % Trim ends to ignore filter response time
    % *******************************************************************
    trim = ceil(length(xOut) * 0.05); % trim first and last 5%
    xOut = xOut(trim+1:end-trim);
    outData = outData(trim+1:end-trim);
    ref = ref(trim+1:end-trim);

    % *******************************************************************
    % Plot
    % *******************************************************************
    figure(1), clf, hold on
    plot(xIn, inData, 'k+-');
    plot(xOut, outData, 'b.-');
    stem(xOut, ref, 'r');
    legend('input', 'interpolation result', 'ideal result');
    title('Lagrange interpolation, impulse response');
    axis([25 45 -2 2])

    % *******************************************************************
    % signal-to-noise ratio
    % *******************************************************************
    SNR_dB = snr(ref, outData)