clear, clc

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

fs = 40;
x = 0 : 1/fs : 10;
y = fun(x);

fsn = 50;
xi = 2 : 1/fsn : 8;
yr = fun(xi);

c = 6;

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

% 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;
yi = zeros(1, length(xi));
for i = 1:length(xi)

    t = xi(i);

    while t > times(c/2)
        index = index+1;
        times = (1+index:c+index)/fs - 1/fs;
        values = y(1+index:c+index);
    end

    if find( abs(xi(i)-times) < 1e-15 )
        yi(i) = values(  abs(xi(i)-times) < 1e-15  );
        continue;
    end

    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')

snr(yr,yi)