function eval_blep()

nsin = 2048;
nharm_max = 1800;
L = 48000;
fs = 48000;
fstart = 110;
fend = 110;
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,
    bb(mm, :) = sin((mm-1)*xx*pi);
end;
aa = sin(xx*pi);

for mm=1:nharm
    for nn=1:length(xx)
        if (aa(nn) == 0)
            blitm(mm,nn) = (mm-1);
        else
            blitm(mm,nn) = bb(mm, nn)/aa(nn).*Hw(nn);
        end
    end
end;
for mm=1:nharm
    blepm(mm, :) = cumsum(blitm(mm, :))/nsin;
end;

%blep2 = blepm(nharm_max/2+1, :)'
%aa2 = aa(:)'
%bb2 = bb(nharm_max/2+1, :)'

x = 0.5;
z = 0;
f = fstart;
tri = 0;
for n = 1:L,

	if x >= 0.5,
		x = x - 1;
		p = fs/f;
		dx  = 1.0/p;
		m = min(fix((p/2.0)) + 1, nharm_max);
		z = ~z;
	end;
    nn = (x+0.5)*(nsin-1) + 1;
    ni = fix(nn);
    nf = nn - ni;

    h = lagrange(1, nf);
    
 %   blep = (blepm(m, min(ni+1, nsin)) - blepm(m, ni))*nf + blepm(m, ni);
 
    [blep, Zi] = filter(h, 1, blepm(m
    blep = h(1)*blepm(m, max(1, ni-1)) + h(2)*blepm(m, ni);
    
	saw = x - blep;
	if (z == 0)
        sqr = blep;
    else
        sqr = 1-blep;
    end
    tri = tri + 2*(sqr-0.5)*dx;
	vtri(n) = tri;
	vsaw(n) = saw+0.5;
	vsqr(n) = 0.5*(sqr-0.5);
	vblep(n) = 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