added matlab

matlab/osc/blep/eval_minblep.asv

@@ -0,0 +1,116 @@
+function eval_minblep()
+
+nsin = 4096;
+nharm_max = 1800;
+L = 48000/2;
+fs = 48000;
+fstart = 440;
+fend = 440;
+f2start = 200;
+f2end = 600;
+df = (fend-fstart)/L
+df2 = (f2end-f2start)/L
+saw = 0;
+
+% 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;
+
+w = [1; 2*ones(nsin/2-1,1); ones(1 - rem(nsin,2),1); zeros(nsin/2-1,1)]';
+for mm=1:nharm
+    x_rc = real(ifft(log(abs(fft(blitm(mm, :))))));
+    y = real(ifft(exp(fft(w.*x_rc))));
+    minblepm(mm, :) = cumsum(y)/nsin;
+end;
+
+x = 0;
+x2 = 0;
+z = 0;
+f = fstart;
+f2 = f2start;
+tri = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+startup = 1;
+blep = 1;
+for n = 1:L,
+
+    if x >= 1 || startup,
+        x = 0;
+        p = fs/f;
+        dx  = 1.0/p;
+        m = min(fix((p/2.0)) + 1, nharm_max);
+        blep_offset = 0;
+        blep_gain = 1;
+        z = ~z;
+    end;
+    if x2 >= 1 || startup,
+        x2 =0;
+        dx2  = f2/fs;
+%        blep_offset = -(1-x);
+%        blep_gain = x;
+%        x = 0;
+%        z = ~z;
+    end;
+    startup = 0;
+    
+    nn = (x+0.0)*(nsin-1) + 1;
+    ni = fix(nn);
+    nf = nn - ni;
+    h = lagrange(NLG, nf);
+    blep = h(NLG+1)*minblepm(m, max(1, ni));
+    j = 1;
+    for i = NLG:-1:1
+        blep = blep + h(i)*minblepm(m, max(1, ni-j));
+        j = j + 1;
+    end;
+
+    saw = x - blep_gain*blep + blep_offset;
+	if (z == 0)
+        sqr = blep;
+    else
+        sqr = 1-blep;
+    end
+    tri = tri + 2*(sqr-0.5)*dx;
+	vsqr(n) = 0.5*(sqr-0.5);
+	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;
+    x2 = x2 + dx2;
+	f2 = f2 + df2;
+end;
+		
+
+close all
+plot(1:nsin, blepm(fix(nharm/2), :), 1:nsin, minblepm(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

matlab/osc/blep/eval_minblep.m

@@ -0,0 +1,116 @@
+function eval_minblep()
+
+nsin = 4096;
+nharm_max = 1800;
+L = 48000/2;
+fs = 48000;
+fstart = 440;
+fend = 440;
+df = (fend-fstart)/L
+f2start = 200;
+f2end = 600;
+df2 = (f2end-f2start)/L
+saw = 0;
+
+% 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;
+
+w = [1; 2*ones(nsin/2-1,1); ones(1 - rem(nsin,2),1); zeros(nsin/2-1,1)]';
+for mm=1:nharm
+    x_rc = real(ifft(log(abs(fft(blitm(mm, :))))));
+    y = real(ifft(exp(fft(w.*x_rc))));
+    minblepm(mm, :) = cumsum(y)/nsin;
+end;
+
+x = 0;
+x2 = 0;
+z = 0;
+f = fstart;
+f2 = f2start;
+tri = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+startup = 1;
+blep = 1;
+for n = 1:L,
+
+    if x >= 1 || startup,
+        x = 0;
+        p = fs/f;
+        dx  = 1.0/p;
+        m = min(fix((p/2.0)) + 1, nharm_max);
+        blep_offset = 0;
+        blep_gain = 1;
+        z = ~z;
+    end;
+    if x2 >= 1 || startup,
+        x2 =0;
+        dx2  = f2/fs;
+%        blep_offset = -(1-x);
+%        blep_gain = x;
+%        x = 0;
+%        z = ~z;
+    end;
+    startup = 0;
+    
+    nn = (x+0.0)*(nsin-1) + 1;
+    ni = fix(nn);
+    nf = nn - ni;
+    h = lagrange(NLG, nf);
+    blep = h(NLG+1)*minblepm(m, max(1, ni));
+    j = 1;
+    for i = NLG:-1:1
+        blep = blep + h(i)*minblepm(m, max(1, ni-j));
+        j = j + 1;
+    end;
+
+    saw = x - blep_gain*blep + blep_offset;
+	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;
+    x2 = x2 + dx2;
+	f2 = f2 + df2;
+end;
+		
+
+close all
+plot(1:nsin, blepm(fix(nharm/2), :), 1:nsin, minblepm(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

matlab/osc/blep/eval_minblep_pre.asv

@@ -0,0 +1,33 @@
+function eval_minblep_pre()
+
+f = 440;
+N = 4000;
+t = (-N/2:N/2-1)/N;
+
+x = sin(t*f*pi)./sin(t*pi); .* kaiser(N, 8)';
+x(find(isnan(x))) = f;
+%x = x / f;
+
+% Real cepstrum
+x_rc = real(ifft(log(abs(fft(x)))));
+w = [1; 2*ones(N/2-1,1); ones(1 - rem(N,2),1); zeros(N/2-1,1)]';
+y = real(ifft(exp(fft(w.*x_rc))));
+
+xi = cumsum(x)/N;
+yi = cumsum(y)/N;
+
+close all;
+subplot(2, 1, 1)
+plot(t, x/f, t, xi, 'r'); grid;
+
+subplot(2, 1, 2)
+plot(0:N-1, y/f, 0:N-1, yi, 'r'); grid;
+
+figure;
+plot(t, xi); grid;
+
+figure;
+plot(t, yi); grid;
+
+figure;
+plot(w); grid;

matlab/osc/blep/eval_minblep_pre.m

@@ -0,0 +1,33 @@
+function eval_minblep_pre()
+
+f = 440;
+N = 4000;
+t = (-N/2:N/2-1)/N;
+
+x = sin(t*f*pi)./sin(t*pi);% .* kaiser(N, 8)';
+x(find(isnan(x))) = f;
+%x = x / f;
+
+% Real cepstrum
+x_rc = real(ifft(log(abs(fft(x)))));
+w = [1; 2*ones(N/2-1,1); ones(1 - rem(N,2),1); zeros(N/2-1,1)]';
+y = real(ifft(exp(fft(w.*x_rc))));
+
+xi = cumsum(x)/N;
+yi = cumsum(y)/N;
+
+close all;
+subplot(2, 1, 1)
+plot(t, x/f, t, xi, 'r'); grid;
+
+subplot(2, 1, 2)
+plot(0:N-1, y/f, 0:N-1, yi, 'r'); grid;
+
+figure;
+plot(t, xi); grid;
+
+figure;
+plot(t, yi); grid;
+
+figure;
+plot(w); grid;

matlab/osc/blip/eval_blep.asv

@@ -0,0 +1,96 @@
+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

matlab/osc/blip/eval_blep.m

@@ -0,0 +1,90 @@
+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

matlab/osc/blip/eval_blep_hs.asv

@@ -0,0 +1,105 @@
+function eval_blep_hs()
+
+nsin = 1024;
+nharm_max = 1800;
+L = 48000;
+fs = 48000;
+fstart = 44;
+fend = 44;
+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;
+x2 = 0.5;
+z = 0;
+f = fstart;
+
+tri = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+
+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(NLG, nf);
+    
+ %   blep = (blepm(m, min(ni+1, nsin)) - blepm(m, ni))*nf + blepm(m, ni);
+ 
+%    blep = h(1)*blepm(m, max(1, ni-3)) +  h(2)*blepm(m, max(1, ni-2)) +  h(3)*blepm(m, max(1, ni-1)) +  h(4)*blepm(m, max(1, ni));
+    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
+    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

matlab/osc/blip/eval_blep_hs.m

@@ -0,0 +1,145 @@
+function eval_blep_hs()
+
+nsin = 1024;
+nharm_max = 256;
+L = 4800;
+fs = 48000;
+fstart = 44;
+fend = 44;
+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.0;
+x2 = -0.5;
+z = 0;
+f = fstart;
+f2 = 0.6*f;
+dx2 = f2/fs;
+tri = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+sync = 0;
+restart = 0;
+startup = 1;
+minblep_active = 0;
+minblep_count = 0;
+minblep_length = 200;
+minblep = 1-(0:minblep_length-1)/minblep_length;
+for n = 1:L,
+
+    sync = 0;
+    if startup,
+        x = 0.0;
+        startup = 0;
+		p = fs/f;
+		dx  = 1.0/p;
+		m = min(fix((p/2.0)) + 1, nharm_max);
+		z = 0;
+    end;
+    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;
+	if (x2 >= 0.5),
+		x2 = x2 - 1;
+        sync = 1;
+        startup = 1;
+        minblep_count = minblep_length;
+	end;
+    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;
+    if (minblep_count > 0)
+        saw = saw_last - minblep(minblep_count)/(1-2*saw_last+minblep_length*dx);
+        minblep_count = minblep_count - 1;
+        minblep_active = 1;
+    else
+        minblep_active = 0;
+        saw = x - blep;
+        saw_last = saw;
+    end
+	if (z == 0)
+        sqr = blep;
+    else
+        sqr = 1-blep;
+    end
+    tri = tri + 2*(sqr-0.5)*dx;
+	vtri(n) = tri;
+
+	vsaw(n) = 2*(saw+0.5);
+	vsqr(n) = 0.5*(sqr-0.5);
+	vblep(n) = blep;
+    vx(n) = x;
+    vsync(n) = sync;
+    vminblep_active(n) = minblep_active;
+	x2 = x2 + dx2;
+	x = x + dx;
+	f = f + df;
+
+end;
+		
+
+close all
+plot(1:nsin, blepm(fix(nharm/2), :)); grid;
+
+figure;
+subplot(3, 1, 1)
+plot(1:L, vx, 1:L, vsync, 'r'); grid;
+subplot(3, 1, 2)
+plot(1:L, vblep); grid;
+subplot(3, 1, 3)
+plot(1:L, vsaw, 1:L, vminblep_active); 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

matlab/osc/blip/eval_blit.asv

@@ -0,0 +1,103 @@
+function eval_blit()
+
+sinc_use_lut = 1
+
+L = 40000;
+fs = 48000;
+fstart = 110;
+fend = 1100;
+df = (fend-fstart)/L;
+
+nhw = 4096;
+Hw = kaiser(nhw, 8)';
+
+if sinc_use_lut == 1
+    
+    nharm = 512;
+    nsin = 4096;
+    xx = (0:nsin)/nsin - 0.5;
+    mmm=0;
+    for mm=1:nharm,
+        mmm=mmm+1;
+        bb(mm, :) = sin(mmm*xx*pi);
+    end;
+    aa = sin(xx*pi);
+end;
+
+sinc_m(mm, find(aa == 0)) = 1;
+
+x = 0.5;
+z = -1;
+f = fstart;
+sqr = 0.0;
+c4 = 0;
+a = 1;
+b = 1;
+for n = 1:L,
+
+	if x >= 0.5,
+		x = x - 1;
+		p = fs/f;
+		fraq  = 1.0/p;
+		m = fix((p/2.0)) + 1;
+		saw = 0.0;
+		c3 = 0.0;
+		z = -z;
+	end;
+    if sinc_use_lut == 1
+        nn = (x+0.5)*nsin + 1;
+        ni = fix(nn);
+        nf = nn - ni;
+        bbi = (bb(m, ni+1) - bb(m, ni))*nf + bb(m, ni);     % Fractional delay (linear interpolation)
+        aai = (aa(ni+1) - aa(ni))*nf + aa(ni);              % Fractional delay (linear interpolation)
+        if (aai ~= 0)
+            blit = fraq* bbi / aai  * Hw(fix((x+0.5)*nhw)+1);
+        else
+            blit = fraq;
+        end
+    else
+        b = sin(m*x*pi);
+        a = sin(x*pi);
+        if (a ~= 0)
+            blit = fraq * b/a * Hw(fix((x+0.5)*nhw)+1);
+        else
+            blit = fraq;
+        end
+    end
+    
+	saw = saw + blit;
+	sqr = sqr + z*blit;
+	vsaw(n) = 2*(saw - c3);
+	vsqr(n) = 2*(sqr - 0.5);
+	vblit(n) = blit;
+    vc3(n) = b;
+	x = x + fraq;
+	c3 = c3 + fraq;
+	f = f + df;
+    c4 = c4 + fraq - blit;
+    vc4(n) = a;
+end;
+		
+
+close all
+
+figure;
+plot(1:L, vsaw, 1:L, vc4, 1:L, vblit); grid;
+
+figure;
+plot(vc4); grid;
+
+figure;
+plot(vc3); grid;
+
+if sinc_use_lut == 1
+    figure;
+    plot(sinc_m(nharm, :)); grid;
+end;
+
+wavwrite(0.5*vc4, fs, 'c4.wav');
+wavwrite(0.5*vc3, fs, 'c3.wav');
+wavwrite(0.5*vsqr, fs, 'sqr.wav');
+wavwrite(0.5*vsaw, fs, 'saw.wav');
+wavwrite(0.5*vblit, fs, 'blit.wav');
+

matlab/osc/blip/eval_blit.bak.m

@@ -0,0 +1,116 @@
+function eval_blit()
+Nb = 33;
+Nos = 1000;
+fs = 48000;
+f = 440;
+fcut = 16000;
+P = fs/f
+Pi = floor(P)
+Pf = fix(Nos*(fs/f - Pi))
+
+x = (0:Nb*Nos-1) - (Nb*Nos-1)/2;
+H = sinc(2*fcut/fs.*x/Nos).*Kaiser(Nos*Nb, 8)';
+for ii=0:Nos-1,
+    for jj=0:Nb-1
+        Hs(ii+1,jj+1) = H(Nos-ii+Nos*jj);
+    end
+end;
+
+%Hs = init_steps(Nb, Nos);
+
+k = 1./sum(Hs');
+%k = ones(1, Nb);
+y0 = 0;
+y1 = 0;
+y2 = .5;
+ks = 0.002*f/55
+ii = 0;
+ff = 0;
+nn = 1;
+jj = Nb;
+kk = 1;
+for i = 1:32000,
+    if ii == 0
+        jj = 1;
+        kk = ff + 1;
+        ii = Pi;
+        ff = ff + Pf;
+        if ff > (Nos-1)
+            ff = ff - Nos;
+            ii = ii + 1;
+        end
+    end;
+    ii = ii - 1;
+    if jj <= Nb
+        y0 = y0 + k(kk)*Hs(kk,jj); % - ks*y0;
+        jj = jj + 1;
+    end;
+    y0 = y0 - ks*y0;
+    y1 = 1 - y0;
+    y2 = 0.001*y1 + 0.999*y2;
+    blit(nn) = y1-y2;
+    nn = nn + 1;
+end;
+wavwrite(0.7*blit, fs, 'blit.wav');
+
+for ii = 1:Nos
+    step(ii,:) = k(ii).*filter(Hs(ii,:), 1, [ones(1,Nb)]);
+end
+
+close all
+plot(1:length(blit), blit, '-o'); grid;
+
+figure
+plot(abs(fft(blit))); grid;
+
+figure
+plot(step', '-'); grid;
+
+figure
+plot(Hs(1,:)', '-+'); grid;
+
+function steps = init_steps(step_width, phase_count)
+
+low_pass  = 0.999; % lower values filter more high frequency
+high_pass = 0.990; % lower values filter more low frequency
+
+%phase_count = 32; % number of phase offsets to sample band-limited step at
+%step_width  = 16; % number of samples in each final band-limited step
+
+%steps [phase_count] [step_width]; // would use short for speed in a real program
+
+	% Generate master band-limited step by adding sine components of a square wave
+	master_size = step_width * phase_count;
+	% master [master_size]; // large; might want to malloc() instead
+	for i = 0:master_size-1
+		master(i+1) = 0.5;
+    end;
+    
+    gain = 0.5 / 0.777; % adjust normal square wave's amplitude of ~0.777 to 0.5
+	sine_size = 256 * phase_count + 2;
+	max_harmonic = sine_size / 2 / phase_count;
+	for h = 1:2:max_harmonic
+		amplitude = gain / h;
+		to_angle = 3.14159265358979323846 * 2 / sine_size * h;
+    	for i = 0:master_size-1
+			master(i+1) = master(i+1) + sin( (i - master_size / 2) * to_angle ) * amplitude;
+        end
+		gain = gain * low_pass;
+    end
+	
+	% Sample master step at several phases
+	for phase = 0:phase_count-1
+		error = 1.0;
+		prev = 0.0;
+    	for i = 0:step_width-1
+			cur = master (i * phase_count + (phase_count - 1 - phase)+1);
+			delta = cur - prev;
+			error = error - delta;
+			prev = cur;
+			steps (phase+1, i+1) = delta;
+        end
+		
+		% each delta should total 1.0
+		steps (phase+1, step_width / 2) = steps (phase+1, step_width / 2) + error * 0.5;
+		steps (phase+1, step_width / 2 + 1) = steps (phase+1, step_width / 2 + 1) + error * 0.5;
+    end

matlab/osc/blip/eval_blit.m

@@ -0,0 +1,89 @@
+function eval_blit()
+
+sinc_use_lut = 0
+
+L = 40000;
+fs = 48000;
+fstart = 8*440;
+fend = 8*440;
+df = (fend-fstart)/L
+
+nhw = 4096;
+Hw = kaiser(nhw, 8)';
+
+max_m = fix((fs/fstart/2.0)) + 1
+
+sinc_table_calc = sinc_use_lut;
+
+if sinc_table_calc == 1
+    nharm = max_m
+    nsin = 1024;
+    xx = (0:nsin)/nsin - 0.5;
+    mmm=0;
+    for mm=1:nharm,
+        mmm=mmm+1;
+        bb(mm, :) = sin(mmm*xx*pi);
+    end;
+    aa = sin(xx*pi);
+end;
+
+x = 0.5;
+z = -1;
+f = fstart;
+sqr = 0.0;
+saw = 0.0;
+tri = 0.0;
+a = 1;
+b = 1;
+for n = 1:L,
+
+	if x >= 0.5,
+		x = x - 1;
+		p = fs/f;
+		fraq  = 1.0/p;
+		m = fix((p/2.0)) + 1;
+		z = -z;
+	end;
+    if sinc_use_lut == 1
+        nn = (x+0.5)*nsin + 1;
+        ni = fix(nn);
+        nf = nn - ni;
+        bbi = (bb(m, ni+1) - bb(m, ni))*nf + bb(m, ni);     % Fractional delay (linear interpolation)
+        aai = (aa(ni+1) - aa(ni))*nf + aa(ni);              % Fractional delay (linear interpolation)
+        if (aai ~= 0)
+            blit = fraq* bbi / aai  * Hw(fix((x+0.5)*nhw)+1);
+        else
+            blit = fraq;
+        end
+    else
+        b = sin(m*x*pi);
+        a = sin(x*pi);
+        if (a ~= 0)
+            blit = fraq * b/a * Hw(fix((x+0.5)*nhw)+1);
+        else
+            blit = fraq;
+        end
+    end
+    
+	saw = saw + fraq - blit;
+	sqr = sqr + z*blit;
+    tri = tri + 2*(sqr-0.5)*fraq;
+	vtri(n) = tri;
+	vsaw(n) = saw;
+	vsqr(n) = 2*(sqr - 0.5);
+	vblit(n) = blit;
+	x = x + fraq;
+	f = f + df;
+end;
+		
+
+close all
+
+figure;
+plot(1:L, vsaw, 1:L, vblit); 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*vblit, fs, 'blit.wav');
+

matlab/osc/blip/eval_blit2.m

@@ -0,0 +1,128 @@
+function eval_blit()
+L = 64000;
+Nb = 65;
+Nos = 8;
+fs = 48000;
+fstart = 440;
+fend = 1760;
+df = (fend-fstart)/L
+fcut = 16000;
+
+x = (0:Nb*Nos-1) - (Nb*Nos-1)/2;
+H = sinc(2*fcut/fs.*x/Nos).*Kaiser(Nos*Nb, 8)';
+for ii=0:Nos-1,
+    for jj=0:Nb-1
+        Hs(ii+1,jj+1) = H(Nos-ii+Nos*jj);
+    end
+end;
+
+%Hs = init_steps(Nb, Nos);
+
+f = fstart;
+k = 1./sum(Hs')
+%k = ones(1, Nb);
+y0 = 0;
+y1 = 0;
+y2 = .5;
+ks = 0.002*f/55
+ii = 0;
+ff = 0;
+nn = 1;
+jj = Nb;
+kk = 1;
+for i = 1:L,
+	P = fs/f;
+	Pi = floor(P);
+	Pf = fix(Nos*(fs/f - Pi));
+	
+    if ii == 0
+        jj = 1;
+        kk = ff + 1;
+        ii = Pi;
+        ff = ff + Pf;
+        if ff > (Nos-1)
+            ff = ff - Nos;
+            ii = ii + 1;
+        end
+    end;
+    ii = ii - 1;
+    if jj <= Nb
+        y0 = k(kk)*Hs(kk,jj); % - ks*y0;
+        jj = jj + 1;
+    end;
+    y0 = y0 - ks*y0;
+    y1 = 1 - y0;
+    y2 = 0.001*y1 + 0.999*y2;
+    blit(nn) = y1-y2;			% DC removal
+    nn = nn + 1;
+	f = f + df;
+    
+end;
+wavwrite(0.7*blit, fs, 'blit.wav');
+
+for ii = 1:Nos
+    step(ii,:) = k(ii).*filter(Hs(ii,:), 1, [ones(1,Nb)]);
+end
+
+DC = y2
+
+close all
+plot(H); grid;
+
+figure
+plot(1:length(blit), blit, '-o'); grid;
+
+figure
+plot(abs(fft(blit))); grid;
+
+figure
+plot(step', '-'); grid;
+
+figure
+plot(Hs(1,:)', '-+'); grid;
+
+function steps = init_steps(step_width, phase_count)
+
+low_pass  = 0.999; % lower values filter more high frequency
+high_pass = 0.990; % lower values filter more low frequency
+
+%phase_count = 32; % number of phase offsets to sample band-limited step at
+%step_width  = 16; % number of samples in each final band-limited step
+
+%steps [phase_count] [step_width]; // would use short for speed in a real program
+
+	% Generate master band-limited step by adding sine components of a square wave
+	master_size = step_width * phase_count;
+	% master [master_size]; // large; might want to malloc() instead
+	for i = 0:master_size-1
+		master(i+1) = 0.5;
+    end;
+    
+    gain = 0.5 / 0.777; % adjust normal square wave's amplitude of ~0.777 to 0.5
+	sine_size = 256 * phase_count + 2;
+	max_harmonic = sine_size / 2 / phase_count;
+	for h = 1:2:max_harmonic
+		amplitude = gain / h;
+		to_angle = 3.14159265358979323846 * 2 / sine_size * h;
+    	for i = 0:master_size-1
+			master(i+1) = master(i+1) + sin( (i - master_size / 2) * to_angle ) * amplitude;
+        end
+		gain = gain * low_pass;
+    end
+	
+	% Sample master step at several phases
+	for phase = 0:phase_count-1
+		error = 1.0;
+		prev = 0.0;
+    	for i = 0:step_width-1
+			cur = master (i * phase_count + (phase_count - 1 - phase)+1);
+			delta = cur - prev;
+			error = error - delta;
+			prev = cur;
+			steps (phase+1, i+1) = delta;
+        end
+		
+		% each delta should total 1.0
+		steps (phase+1, step_width / 2) = steps (phase+1, step_width / 2) + error * 0.5;
+		steps (phase+1, step_width / 2 + 1) = steps (phase+1, step_width / 2 + 1) + error * 0.5;
+    end

matlab/osc/blip/eval_blit3.m

@@ -0,0 +1,68 @@
+function eval_blit()
+L = 64000;
+Nb = 137;
+Nphases = 256;
+fs = 48000;
+fstart = 110*4;
+fend = 110*4;
+df = (fend-fstart)/L;
+fcut = 16000;
+
+x = (0:Nb*Nphases-1) - (Nb*Nphases-1)/2;
+H = sinc(2*fcut/fs.*x/Nphases);%.*Kaiser(Nphases*Nb, 8)';
+Hw = Kaiser(Nb, 18)';
+
+for ii=0:Nphases-1,
+    for jj=0:Nb-1
+        Hs(ii+1,jj+1) = H(Nphases-ii+Nphases*jj);
+    end
+end;
+
+
+f = fstart;
+Pic = 1;
+blit(1:L) = zeros(1, L);
+blit_bp(1:L) = zeros(1, L);
+pol = 1;
+for n = Nb:L,
+	P = fs/f;
+	Pi = floor(P);
+	Pf = 1+fix(Nphases*(P - Pi));
+
+	if Pic >= Pi,
+		blit(n-Nb+1:n) = blit(n-Nb+1:n) + Hs(Pf,:).*Hw;
+		blit_bp(n-Nb+1:n) = blit_bp(n-Nb+1:n) + pol*Hs(Pf,:).*Hw;
+		Pic = 0;
+		pol = -pol;
+	end;
+	Pic = Pic + 1;
+	f = f + df;
+end;
+
+mean_blit = 0.0034;
+mean_saw =  0.00;
+
+saw(1:L) = zeros(1, L);
+sqr(1:L) = zeros(1, L);
+y_saw = 0;
+y_sqr = 0;
+for n = 2:L,
+	y_saw = y_saw + blit(n)-mean_blit;
+	y_sqr = y_sqr + blit_bp(n);
+	saw(n) = -y_saw -1;
+	sqr(n) = y_sqr;
+end;
+
+saw_nodc = filter([1 -1], [1 -0.995], saw);
+sqr_nodc = filter([1 -1], [1 -0.995], sqr);
+
+close all;
+plot(saw); grid;
+
+figure;
+plot(saw_nodc); grid;
+
+wavwrite(0.8*blit, fs, 'blit.wav');
+wavwrite(0.7*saw_nodc, fs, 'saw.wav');
+wavwrite(0.5*sqr_nodc, fs, 'sqr.wav');
+

matlab/osc/blip/eval_blit4.asv

@@ -0,0 +1,92 @@
+function eval_blit()
+
+sinc_use_lut = 0
+
+L = 40000;
+fs = 48000;
+fstart = 110;
+fend = 1550;
+df = (fend-fstart)/L;
+
+nhw = 4096;
+Hw = kaiser(nhw, 8)';
+
+ya = zeros(3,1);
+yb = zeros(3,1);
+
+x = 0.5;
+z = -1;
+f = fstart;
+sqr = 0.0;
+c4 = 0;
+for n = 1:L,
+
+    
+	if x >= 0.5,
+		x = x - 1;
+		p = fs/f;
+		fraq  = 1.0/p;
+		m = fix((p/2.0)) + 1;
+		saw = 0.0;
+		c3 = 0.0;
+		z = -z;
+        ip = 0; % phase of the first output sample in radians
+        w = f*pi / fs
+        ba1 = 2.0 * cos(w)
+        bb1 = 2.0 * cos(w)
+        
+	end;
+    
+    if sinc_use_lut == 1
+        xi = fix((x+0.5)*nsin)+1
+        if  aa(xi) ~= 0
+            blit = fraq*sinc_m(m, xi)  * Hw(fix((x+0.5)*nhw)+1);;
+        else
+            blit = fraq;
+        end
+    else
+        b = fraq*sin(m*x*pi);
+        a = sin(x*pi);
+        if (a ~= 0)
+            blit = b/a * Hw(fix((x+0.5)*nhw)+1);
+        else
+            blit = fraq;
+        end
+    end
+    
+	saw = saw + blit;
+	sqr = sqr + z*blit;
+	vsaw(n) = 2*(saw - c3);
+	vsqr(n) = 2*(sqr - 0.5);
+	vblit(n) = blit;
+    vc3(n) = c3;
+	x = x + fraq;
+	c3 = c3 + fraq;
+	f = f + df;
+    c4 = c4 + fraq - blit;
+    vc4(n) = c4;
+end;
+		
+
+close all
+
+figure;
+plot(1:L, vsaw, 1:L, vc4, 1:L, vblit); grid;
+
+figure;
+plot(vblit); grid;
+
+figure;
+plot(vc3); grid;
+
+if sinc_use_lut == 1
+    figure;
+    plot(sinc_m(nharm, :)); grid;
+end;
+
+wavwrite(0.5*vc4, fs, 'c4.wav');
+wavwrite(0.5*vc3, fs, 'c3.wav');
+wavwrite(0.5*vsqr, fs, 'sqr.wav');
+wavwrite(0.5*vsaw, fs, 'saw.wav');
+wavwrite(0.5*vblit, fs, 'blit.wav');
+

matlab/osc/blip/eval_blit4.m

@@ -0,0 +1,88 @@
+function eval_blit4()
+
+L = 4000;
+fs = 48000;
+fstart = 110;
+fend = 110;
+df = (fend-fstart)/L;
+
+nhw = 4096;
+Hw = kaiser(nhw, 8)';
+
+ya = zeros(3,1);
+yb = zeros(3,1);
+
+x = 0.5;
+z = -1;
+f = fstart;
+sqr = 0.0;
+c4 = 0;
+
+for n = 1:L,
+
+    
+	if x >= 0.5,
+		x = x - 1;
+		p = fs/f;
+		fraq  = 1.0/p;
+		m = fix((p/2.0)) + 1
+		saw = 0.0;
+		c3 = 0.0;
+		z = -z;
+        ip = -pi/2; % phase of the first output sample in radians
+        w = f*pi / fs;
+        ba1 = 2.0 * cos(w);
+        bb1 = 2.0 * cos(m*w);
+        ya(2)=sin(ip-w);
+        ya(3)=sin(ip-2*w);
+        yb(2)=sin(ip-m*w);
+        yb(3)=sin(ip-2*m*w);
+	end;
+    
+    ya(1) = ba1*ya(2) - ya(3);
+    ya(3) = ya(2);
+    ya(2) = ya(1);
+
+    yb(1) = bb1*yb(2) - yb(3);
+    yb(3) = yb(2);
+    yb(2) = yb(1);
+
+    b = fraq*yb(1);
+    a = ya(1);
+    if (a ~= 0)
+        blit = b/a * Hw(fix((x+0.5)*nhw)+1);
+    else
+        blit = fraq;
+    end
+    
+	saw = saw + blit;
+	sqr = sqr + z*blit;
+	vsaw(n) = 2*(saw - c3);
+	vsqr(n) = 2*(sqr - 0.5);
+	vblit(n) = blit;
+    vc3(n) = b;
+	x = x + fraq;
+	c3 = c3 + fraq;
+	f = f + df;
+    c4 = c4 + fraq - blit;
+    vc4(n) = a;
+end;
+		
+
+close all
+
+figure;
+plot(1:L, vsaw, 1:L, vc4, 1:L, vblit); grid;
+
+figure;
+plot(vc4); grid;
+
+figure;
+plot(vc3); grid;
+
+wavwrite(0.5*vc4, fs, 'c4.wav');
+wavwrite(0.5*vc3, fs, 'c3.wav');
+wavwrite(0.5*vsqr, fs, 'sqr.wav');
+wavwrite(0.5*vsaw, fs, 'saw.wav');
+wavwrite(0.5*vblit, fs, 'blit.wav');
+

matlab/osc/eval_hardsync.asv

@@ -0,0 +1,125 @@
+function eval_hardsync()
+
+nsin = 4096;
+nharm_max = 1800;
+L = 48000/2;
+fs = 48000;
+fstart = 440;
+fend = 440;
+df = (fend-fstart)/L
+f2start = 400;
+f2end = 400;
+df2 = (f2end-f2start)/L
+
+% Calc blep table
+nharm = min(fix((fs/min(fstart, f2start)/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;
+x2 = 0.0;
+z = 0;
+f = fstart;
+f2 = f2start;
+tri = 0;
+saw = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+startup = 1;
+
+for n = 1:L,
+
+	if x >= 1 || startup;
+		x = 0;
+		p = fs/f;
+		dx  = 1.0/p;
+		m = min(fix((p/2.0)) + 1, nharm_max);
+        blep_offset = 0;
+        blep_gain = 1;
+		z = ~z;
+	end;
+    if x2 >= 1 || startup,
+        x2 =0;
+		p2 = fs/f2;
+		dx2  = 1.0/p2;
+		m2 = min(fix((p2/2.0)) + 1, nharm_max);
+%        blep_offset = -(1-(x+0.5));
+%        blep_gain = x+0.5;
+%        x = 0;
+%        z = ~z;
+    end;
+    startup = 0;
+
+    % BLEP 1
+    nn = (x+0.0)*(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;
+    
+    % BLEP 2
+    nn = (x2+0.0)*(nsin-1) + 1;
+    ni = fix(nn);
+    nf = nn - ni;
+    h = lagrange(NLG, nf);
+    blep2 = h(NLG+1)*blepm(m2, max(1, ni));
+    j = 1;
+    for i = NLG:-1:1
+        blep2 = blep2 + h(i)*blepm(m2, 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,1) = sqr;
+    vtri(n,2) = tri;
+	vsaw(n,1) = x;
+	vsaw(n,2) = tri;
+	vblep(n,1) = x;
+	vblep(n,2) = blep;
+	x = x + dx;
+	f = f + df;
+    x2 = x2 + dx2;
+	f2 = f2 + df2;
+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

matlab/osc/eval_hardsync.m

@@ -0,0 +1,107 @@
+function eval_hardsync()
+
+nsin = 4096;
+nharm_max = 1800;
+L = 48000/2;
+fs = 48000;
+fstart = 440;
+fend = 440;
+df = (fend-fstart)/L
+f2start = 200;
+f2end = 600;
+df2 = (f2end-f2start)/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;
+x2 = 0.0;
+z = 0;
+f = fstart;
+f2 = f2start;
+tri = 0;
+NLG = 2;
+h = lagrange(1, 0.35)
+startup = 1;
+
+for n = 1:L,
+
+	if x >= 1 || startup;
+		x = x - 1;
+		p = fs/f;
+		dx  = 1.0/p;
+		m = min(fix((p/2.0)) + 1, nharm_max);
+		z = ~z;
+	end;
+    if x2 >= 1 || startup,
+        x2 =0;
+		p2 = fs/f2;
+		dx2  = 1.0/p2;
+		m2 = min(fix((p2/2.0)) + 1, nharm_max);
+%        blep_offset = -(1-x);
+%        blep_gain = x;
+%        x = 0;
+%        z = ~z;
+    end;
+    startup = 0;
+    nn = (x+0.0)*(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;
+    x2 = x2 + dx2;
+	f2 = f2 + df2;
+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

matlab/osc/eval_harmonics.asv

@@ -0,0 +1,11 @@
+function eval_harmonics(name)
+
+base = exp(1);
+ke = 10;
+
+[x, fs, nbits] = wavread(name);
+
+y = (1 - base.^(-ke*fabs(x))) * sign(x);
+
+
+wavrite(y, fs, nbits, 

matlab/osc/eval_harmonics.m

@@ -0,0 +1,11 @@
+function eval_harmonics(name)
+
+base = exp(1);
+ke = 2;
+
+[x, fs, nbits] = wavread(name);
+
+y = (1 - base.^(-ke*abs(x))) .* sign(x);
+
+yname = sprintf('out.wav', name);
+wavwrite(y*0.8, fs, nbits, yname);