function lagrangeResamplingDemo()
    originDefinition = 0; % see comment for LagrangeBasisPolynomial() below

    % Regarding order: From [1]
    % "Indeed, it is easy to show that whenever Q is odd, none of the
    % impulse responses corresponding to Lagrange interpolation can have
    % linear phase"
    % Here, order = Q+1 => odd orders are preferable
    order = 5;

    % *******************************************************************
    % Set up signals
    % *******************************************************************
    nIn = order + 1;
    nOut = 80;

    % Reference data: from [1] fig. 8, linear-phase type
    % ref = [-0.032, -0.056, -0.064, -0.048, 0, ...
    %       0.216, 0.448, 0.672, 0.864, 1, ...
    %       0.864, 0.672, 0.448, 0.216, 0, ...
    %       -0.048, -0.064, -0.056, -0.032, 0];
    % tRef = (1:size(ref, 2)) / size(ref, 2);
    tRef = [0:1/10:1-1/10];
    ref = sin( tRef*2*pi );

    % impulse input to obtain impulse response
    %inData = zeros(1, nIn);
    %inData(1) = 1;
    inData = sin( [0:1/6:1-1/6] *2*pi );
    nIn = length(inData);

    % Get Lagrange polynomial coefficients
    c0 = lagrangeBasisPolynomial(order, 0, originDefinition);
    c1 = lagrangeBasisPolynomial(order, 1, originDefinition);
    c2 = lagrangeBasisPolynomial(order, 2, originDefinition);
    c3 = lagrangeBasisPolynomial(order, 3, originDefinition);
    c4 = lagrangeBasisPolynomial(order, 4, originDefinition);
    c5 = lagrangeBasisPolynomial(order, 5, originDefinition);

    % *******************************************************************
    % Plot
    % *******************************************************************
    figure(2), clf, hold on
    stem(tRef, ref, 'r');

    % *******************************************************************
    % Apply Lagrange interpolation to input data
    % *******************************************************************
    for i=0:nOut-1

        % outTimeAtInputInteger
        temp = floor( (i/nOut)*nIn );

        % the value of the input sample that appears at FIR tap 'ixTap'
        d0 = inData( mod(temp + 0 - order, nIn) + 1 );
        d1 = inData( mod(temp + 1 - order, nIn) + 1 );
        d2 = inData( mod(temp + 2 - order, nIn) + 1 );
        d3 = inData( mod(temp + 3 - order, nIn) + 1 );
        d4 = inData( mod(temp + 4 - order, nIn) + 1 );
        d5 = inData( mod(temp + 5 - order, nIn) + 1 );

        % Evaluate the Lagrange polynomials, resulting in the time-varying FIR tap weight
        evalFracTime = -0.5 + mod(i,nOut/nIn) * (nIn/nOut);

        cTap0 = polyval(c0, evalFracTime);
        cTap1 = polyval(c1, evalFracTime);
        cTap2 = polyval(c2, evalFracTime);
        cTap3 = polyval(c3, evalFracTime);
        cTap4 = polyval(c4, evalFracTime);
        cTap5 = polyval(c5, evalFracTime);

        % FIR operation: multiply FIR tap weight with input sample
        outData(i+1) = d0*cTap0 + d1*cTap1 + d2*cTap2 + d3*cTap3 + d4*cTap4 + d5*cTap5;

        plot( i/nOut, outData(i+1), 'b.');
        pause(0.01)
    end

    legend('impulse response', 'reference result');
    title('Lagrange interpolation, impulse response');

end