function eval_blep()

nsin = 4096;
nharm_max = 1800;
L = 48000/2;
fs = 48000;
fstart = 440;
fend = 440;
df = (fend-fstart)/L

% Calc blep table
nharm = min(fix((fs/fstart/2.0)) + 1, nharm_max)
Hw = kaiser(nsin, 8);
xx = (0:nsin-1)/nsin - 0.5;

for mm=1:nharm
    blitm(mm,:) = sin(xx*(mm-1)*pi)./sin(xx*pi).*Hw';
    blitm(mm,(find(isnan(blitm(mm,:))))) = (mm-1);
end;

for mm=1:nharm
    blepm(mm, :) = cumsum(blitm(mm, :))/nsin;
end;

x = 0.0;
z = 0;
f = fstart;
tri = 0;
NLG = 2;
h = lagrange(1, 0.35)
startup = 1;

for n = 1:L,

	if x >= 0.5 || startup;
		x = x - 1;
		p = fs/f;
		dx  = 1.0/p;
		m = min(fix((p/2.0)) + 1, nharm_max);
		z = ~z;
	end;
    startup = 0;
    nn = (x+0.5)*(nsin-1) + 1;
    ni = fix(nn);
    nf = nn - ni;
    
    h = lagrange(NLG, nf);
    blep = h(NLG+1)*blepm(m, max(1, ni));
    j = 1;
    for i = NLG:-1:1
        blep = blep + h(i)*blepm(m, max(1, ni-j));
        j = j + 1;
    end;
	saw = x - blep;
	if (z == 0)
        sqr = blep;
    else
        sqr = 1-blep;
    end
	vsqr(n) = 0.5*(sqr-0.5);
    tri = tri + 2*(sqr-0.5)*dx;
	vtri(n) = tri;
	vsaw(n,1) = x;
	vsaw(n,2) = saw+0.5;
	vblep(n,1) = x;
	vblep(n,2) = blep;
	x = x + dx;
	f = f + df;
end;
		

close all
plot(1:nsin, blepm(fix(nharm/2), :)); grid;

wavwrite(0.5*vtri, fs, 'tri.wav');
wavwrite(0.5*vsqr, fs, 'sqr.wav');
wavwrite(0.5*vsaw, fs, 'saw.wav');
wavwrite(0.5*vblep, fs, 'blep.wav');

function h = lagrange(N, delay)
    %LAGRANGE h=lagrange(N,delay) returns order N FIR
    % filter h which implements given delay
    % (in samples). For best results,
    % delay should be near N/2 +/- 1.
    n = 0:N;
    h = ones(1,N+1);
    for k = 0:N
        index = find(n ~= k);
        h(index) = h(index) * (delay-k)./ (n(index)-k);
    end