Implementing this MATLAB function as a filter (Lagrange interpolation)

Below is some code I wrote that uses peicewise barycentric lagrange
interpolation to change the sampling rate of a signal from 9Hz to 30Hz. My
numerical implementation below works far better than any Farrow or variable
fractional delay filter built in MATLAB. How can implement my function
below as a filter? Or how can I build it from blocks in Simulink?

-------------


clear, clc

fun = inline('sin(2*pi*x)');

fs = 9; % original sampling rate
x = 0 : 1/fs : 10; % original sample points
y = fun(x); % original sample values

fsn = 30; % new sampling rate
xi = 2 : 1/fsn : 8; % new sample points
yr = fun(xi); % new sample values, reference
yi = zeros(1, length(xi)); % new sample values, interpolated

c = 6; % Number of points to use for Lagrange interpolation

% starting times and values
times = x(1:c); % buffer, sample times
values = y(1:c); % buffer, sample values

% Evaluate barycentric weights
w = zeros(1, c); 
for j = 1:c
    prod = 1;
    for k = 1:c
        if ( k ~= j )
            prod = prod * ( times(j) - times(k) );
        end
    end
    w(j) = 1/prod;
end

% Evaluate points
index = 0;
for i = 1:length(xi)

    t = xi(i); % current sample time
    
    % Advance the buffer to center the current sample time
    while t > times(c/2)
        index = index+1;
        times = (1+index:c+index)/fs - 1/fs;
        values = y(1+index:c+index);
    end

    % if the current sample time matches an old one, use the same y value
    if find( abs(xi(i)-times) < 1e-15 )
        yi(i) = values(  abs(xi(i)-times) < 1e-15  );
        continue;
    end
    
    % Calculate the points
    num = 0;
    dem = 0;
    for j = 1:c
        temp = ( w(j) / (xi(i)-times(j)) );
        num = num + temp * values(j);
        dem = dem + temp;
    end
    yi(i) = num/dem;
    
end

plot(xi, yr, 'b', xi, yi, 'r')

Hello,

For a given order, you can write out the equation for each coefficient
('barycentric weights' in your code).

Each is a polynomial with the fractional position of an output sample in
the input stream as variable. 
This allows to map it to a generic Farrow resampler.

My implementation is here
http://www.dsprelated.com/showcode/208.php
but it should do more or less the same as the one you've posted.

BTW, whether or not it really is the best solution depends very much on the
bandwidth of the signal. 

For any specific bandwidth, you should get even better results by using
coefficients from the function gPulseCoefs() here:
http://www.dsprelated.com/showcode/236.php
but their calculation is more involved.

I rewrote your implementation (didn't change polyCoeff()) so that it
functions as though the data was being streamed through the filter in the
time domain (like how a real filter would operate). I believe I have that
part working correctly, but I can't figure out how to apply this filter to
functions other than the impulse. 

-----

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 = 3;
   
    % *******************************************************************
    % Set up signals
    % *******************************************************************  
 
    nIn = order + 1;
    nOut = 20*pi;

    % 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);
   
    % impulse input to obtain impulse response
    inData = zeros(1, nIn);
    inData(1) = 1;
   
    % Get Lagrange polynomial coefficients
    c0 = lagrangeBasisPolynomial(order, 0, originDefinition);
    c1 = lagrangeBasisPolynomial(order, 1, originDefinition);
    c2 = lagrangeBasisPolynomial(order, 2, originDefinition);
    c3 = lagrangeBasisPolynomial(order, 3, 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 );
        
        % 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);
        
        % FIR operation: multiply FIR tap weight with input sample
        outData = d0*cTap0 + d1*cTap1 + d2*cTap2 + d3*cTap3;
        
        plot( i/nOut, outData, 'b.');
        pause(0.01)
    end

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

end

>> ...but I can't figure out how to apply this filter to functions other
than the impulse. 
 
Instead of 
inData = zeros(1, nIn);
inData(1) = 1;

put the samples into inData.

I tried that and it's seems so close to work. But when the number of input
samples matches the order+1, my results are 180deg out of phase. Changing
the number of input samples also changes the phase.

For example:
order = 5
inData = sin( [0:1/6:1-1/6] *2*pi );
plot: http://imgur.com/H9oo2
function: http://pastebin.com/fHN7zj39

Sorry, that was a stupid question. It's obviously just the group delay from
the filter. But I am concerned with the quality of the filter after
accounting for the group delay.

Here is the code from my OP
http://pastebin.com/TZaz804c

Here is my modification of your code
http://pastebin.com/FKghzXk0

The function snr() is defined as

function s = snr(y_ref, y)
    s = 20*log10(norm(y_ref(:))/norm(y_ref(:)-y(:)));
end


When using my original code (for the signal, order, and sampling rate
chosen in the above pasted examples), I get an SNR of 62.9dB. When using
this filter method, I get 17.25dB. It looks like it's performing on par
with most of the Farrow filter implementations I've seen. If I'm not
mistaken, your implementation is just a standard Lagrange Farrow filter? In
my experience, Lagrange Farrow filters don't perform nearly as well as the
their numerical counterparts. In the code from my OP, a 4x oversampled
signal that is interpolated by 9th-order Lagrange has an SNR of over 100dB.
I've never seen anything close to same level of accuracy from a Farrow
filter operating under similar conditions. I'm hoping that it's possible
and I've just been implementing my filters incorrectly.

Hello,

The waveforms aren't lined up accurately.
Use 

xr = [0:10/nOut:10-10/nOut]+0.00501; % fsn = 50Hz

The 0.00501 is found by trial-and-error and not perfect, but 'close enough'
for a good visual match.

There is another problem:
In addition to the above line, put the following 

    figure()
    hold on;
    plot(yr(1:end-4), '+-')
    plot(outData(5:end), 'k*-')

to the end. You'll notice that the interpolation is now in general very
accurate, but there are a few major errors. What exactly causes them, I
can't say right now.

Hello, again.

I put your example into my original code. It works fine, and I get a SNR of
66.94 dB with order = 3. Probably some error happened during the rewrite.

The new version can be found here:
http://www.dsprelated.com/blogimages/MarkusNentwig/sn_lagrangeResampling/lagrange2.m

Cheers

Markus

I rewrote your new version to work in the time domain (I'd like to
implement in hardware at some point), and this time is actually works.
http://pastebin.com/TXqgZdpS

I have a question about the choice of input samples to use in the
calculation of each output sample. It isn't always a linear situation.
Depending on the new sampling rate, you might use the input samples in
order:

1 2 3 4
2 3 4 5
3 4 5 6
4 5 6 7

or you might skip a step

1 2 3 4
2 3 4 5
4 5 6 7
5 6 7 8

or you might hold a step

1 2 3 4
2 3 4 5
2 3 4 5
3 4 5 6

I'm not sure how to implement this in hardware (or Simulink first). I was
looking at your diagram in "Polyphase filter / Farrows interpolation"
http://www.dsprelated.com/blogimages/MarkusNentwig/dspblog/polyphase/img103.gif

There isn't really a mechanism to decide when to hold a set of samples for
an extra calculation or when to skip a set. It's easy enough in MATLAB
because you can do random access on any of the input samples, but that's
not the case in the real world.

I was trying to simulate holding a set by doing something like this (for
3rd order)

samples = [1 2 3 4];
if xOutFrac > 0.2
    samples = circshift(samplearr,[1 -1]);
    samples(end) = samplearr(end-1)+1;
end

But the value of the fractional delay which causes a hold appears to be
dependent on the new sampling rate (though I'm not sure, I have no theory
to back this up).

Also thank you for help so far! I'm much closer to getting this to work in
Simulink, if only we can figure out the sample selection. A circular buffer
is easy to implement, but when to advance it is the trick, I think.

Hi,

take an accumulator with a sufficient number of fractional bits (to the
right side of the decimal point).
For each output sample, add an increment that depends on the sample rate:
less than one for interpolation, greater than one for decimation.

The lower accumulator bits (on the right side of the decimal point) goes
into the Lagrange interpolation: A value 0..1 that tells, where to
interpolate between samples.

And the higher bits tell the position in the input stream: If the value
changes by one, load one new input sample. If the value changes by two,
shift two new input samples into the delay chain etc. 
In hardware, you can possibly limit the implementation to 0 and 1, unless
you plan to less than half the original rate.

BTW, you can read my code exactly as explained above. The only difference
is that it doesn't use a for-next loop, similar to hardware, but handles
all output samples in parallel.