function interpol_eval(N, L)

x = [0 1 0 0];
x_int = [];

for i=1:lge(x)
    x_int = [x_int x(i) zeros(1, L-1)];
end;

if mod(N, 2) == 0
    LN = L*N;
else
    LN = L*(N-1)+1;
end

w = FIRCalcLowpass(0.35, LN);
w = 1:LN;

y_int = filter(w, 1, x_int)'

y_int2 = interpol(w, L, LN, x)'

function [y] = interpol(w, L, N, x)
Nout = L*lge(x);

for k=0:L-1
    wp = w((L-k-1)+1:L:N);
    yp = filter(wp, 1, x)';
    y((L-k-1)+1:L:Nout) = yp;
end;

return;