diff --git a/matlab/env/eval_adsr.asv b/matlab/env/eval_adsr.asv
new file mode 100755
index 0000000..639a5c2
--- /dev/null
+++ b/matlab/env/eval_adsr.asv
@@ -0,0 +1,53 @@
+function eval_adsr(Ta, Td, Vs, Tr, ks)
+
+% eval_adsr(0.1, 0.2, 0.2, 0.2)
+
+fs = 1000;
+
+a_a = -log(1-ks)/(Ta*fs);
+a_d = 5/(Td*fs);
+a_r = 5/(Tr*fs);
+state = 0;
+s = 0;
+tol = 1E-4;
+n = 1;
+p = 0;
+k = 1/ks;
+while (1),
+    
+    if state == 0       % Attack
+        s = a_a + (1-a_a)*s;
+        y(n) = k*s;
+        if s >= (ks-tol)
+            state = 1;
+            s = 1-ks;
+        end;
+    end;
+    if state == 1       % Decay
+        s = (1-a_d)*s;
+        y(n) = s/(1-ks);
+        if s <= (Vs+tol)
+            state = 2;
+            p = 0;
+        end;
+    end;
+    if state == 2       % sustain
+        y(n) = Vs;
+        if p >= 1*fs
+            state = 3;
+            s = Vs;
+        end;
+        p = p + 1;
+    end;
+    if state == 3       % Release
+        s = (1-a_r)*s;
+        y(n) = s;
+        if s <= (tol)
+            break;
+        end;
+    end;
+    n = n + 1;
+end;
+
+plot((0:n-1)/fs, y); grid;
+set(gca,'xtick',[0:0.1:(n-1)/fs])
diff --git a/matlab/env/eval_adsr.m b/matlab/env/eval_adsr.m
new file mode 100755
index 0000000..f1da472
--- /dev/null
+++ b/matlab/env/eval_adsr.m
@@ -0,0 +1,52 @@
+function eval_adsr(Ta, Td, Vs, Tr, ks)
+
+% eval_adsr(0.1, 0.2, 0.2, 0.2)
+
+fs = 1000;
+
+a_a = -log(1-ks)/(Ta*fs);
+a_d = 5/(Td*fs);
+a_r = 5/(Tr*fs);
+state = 0;
+s = 0;
+tol = 1E-4;
+n = 1;
+p = 0;
+k = 1/ks;
+while (1),
+    
+    if state == 0       % Attack
+        s = a_a + (1-a_a)*s;
+        y(n) = s/ks;
+        if s >= (ks-tol)
+            state = 1;
+        end;
+    end;
+    if state == 1       % Decay
+        s = (1-a_d)*s;
+        y(n) = s/(ks);
+        if s <= (Vs*ks+tol)
+            state = 2;
+            p = 0;
+        end;
+    end;
+    if state == 2       % sustain
+        y(n) = Vs;
+        if p >= 1*fs
+            state = 3;
+            s = Vs;
+        end;
+        p = p + 1;
+    end;
+    if state == 3       % Release
+        s = (1-a_r)*s;
+        y(n) = s;
+        if s <= (tol)
+            break;
+        end;
+    end;
+    n = n + 1;
+end;
+
+plot((0:n-1)/fs, y); grid;
+set(gca,'xtick',[0:0.1:(n-1)/fs])
diff --git a/matlab/env/eval_adsr.m.bak b/matlab/env/eval_adsr.m.bak
new file mode 100755
index 0000000..292415e
--- /dev/null
+++ b/matlab/env/eval_adsr.m.bak
@@ -0,0 +1,56 @@
+function eval_adsr(Ta, Td, Vs, Tr)
+
+% eval_adsr(0.1, 0.2, 0.2, 0.2)
+
+fs = 10000;
+
+aa = 1/(Ta*fs)
+ba = (1-aa)
+a_d = 5/(Td*fs)
+a_r = 5/(Tr*fs)
+state = 0;
+s = 0;
+tol = 1E-4;
+ke = 1-exp(-1)
+n = 0;
+p = 0;
+k = 1/ke;
+kk = 0;
+s
+while (1),
+    
+    if state == 0
+        s = aa + ba*s;
+        if s >= (ke-tol)
+            state = 1;
+            k = (1-Vs)/ke;
+            kk = Vs;
+        end;
+    end;
+    if state == 1
+        s = (1-a_d)*s;
+        if s <= (tol)
+            state = 2;
+            s = 1;
+            k = Vs;
+            kk = 0;
+            p = 0;
+        end;
+    end;
+    if state == 2
+        if p >= 0.1*fs
+            state = 3;
+        end;
+        p = p + 1;
+    end;
+    if state == 3
+        s = (1-a_r)*s;
+        if s <= (tol)
+            break;
+        end;
+    end;
+    n = n + 1;
+    y(n) = k*s + kk;
+end;
+
+plot((0:n-1)/fs, 0.8-y); grid;
\ No newline at end of file
diff --git a/matlab/env/eval_adsr2.asv b/matlab/env/eval_adsr2.asv
new file mode 100755
index 0000000..7f35ca1
--- /dev/null
+++ b/matlab/env/eval_adsr2.asv
@@ -0,0 +1,52 @@
+function eval_adsr(Ta, Td, Vs, Tr)
+
+% eval_adsr(0.1, 0.2, 0.2, 0.2)
+
+fs = 1000;
+
+a_a = 1/(Ta*fs);
+a_d = -1/(Td*fs/log(Vs))
+a_d = 5/(Td*fs)
+a_r = 5/(Tr*fs)
+state = 0;
+s = 0;
+tol = 1E-4;
+target = 1;
+n = 0;
+p = 0;
+k = 1/0.63;
+while (1),
+    
+    if state == 0       % Attack
+        s = a_a + (1-a_a)*s;
+        if s >= (0.63-tol)
+            state = 1
+            target = Vs;
+        end;
+    end;
+    if state == 1       % Decay
+        s = (1-a_d)*s;
+        if s <= (0.63*Vs+tol)
+            state = 2
+            target = 0;
+            p = 0;
+        end;
+    end;
+    if state == 2       % sustain
+        s = Vs/k;
+        if p >= 1*fs
+            state = 3
+        end;
+        p = p + 1;
+    end;
+    if state == 3       % Release
+        s = (1-a_r)*s;
+        if s <= (target+tol)
+            break;
+        end;
+    end;
+    n = n + 1;
+    y(n) = k*s;
+end;
+
+plot((0:n-1)/fs, y); grid;
\ No newline at end of file
diff --git a/matlab/filter/eval_vcf_lut.asv b/matlab/filter/eval_vcf_lut.asv
new file mode 100755
index 0000000..e408946
--- /dev/null
+++ b/matlab/filter/eval_vcf_lut.asv
@@ -0,0 +1,18 @@
+function eval_vcf_lut()
+N = 1000;
+fs = 48000;
+fmin = 18;
+fmax = 18000;
+sweepbase = 10;
+
+filterParam = struct('sweepbase', 10, 'omega_min', 2*fmin/fs, 'omega_max', 2*fmax/fs, 'num_octaves', log(fmax/fmin)/log(sweepbase));
+filterParam
+
+cv = (0:N-1)/N;
+plot (c
+function cv = omega2cv(filterParam, omega)
+	cv = log(omega/filterParam.omega_min)/(filterParam.num_octaves * log(filterParam.sweepbase));
+	cv = min(1, max(0, cv));
+
+function omega = cv2omega(filterParam, cv)
+	omega = filterParam.omega_min*filterParam.sweepbase ^ (filterParam.num_octaves*min(1, max(0, cv))); 
diff --git a/matlab/filter/eval_vcf_lut.m b/matlab/filter/eval_vcf_lut.m
new file mode 100755
index 0000000..149bd36
--- /dev/null
+++ b/matlab/filter/eval_vcf_lut.m
@@ -0,0 +1,23 @@
+function eval_vcf_lut()
+N = 1000;
+fs = 48000;
+fmin = 18;
+fmax = 18000;
+sweepbase = 10;
+
+filterParam = struct('sweepbase', 10, 'omega_min', fmin/fs, 'omega_max', fmax/fs, 'num_octaves', log(fmax/fmin)/log(sweepbase));
+filterParam
+
+cv = (0:N-1)/N;
+
+subplot(2, 1, 1)
+plot (cv, fs*cv2omega(filterParam, cv)); grid; legend('F'); xlabel('cv');
+subplot(2, 1, 2)
+plot (cv, cos(2*pi*cv2omega(filterParam, cv)), cv, sin(2*pi*cv2omega(filterParam, cv))); grid; legend('kc', 'ks'); xlabel('cv');
+
+function cv = omega2cv(filterParam, omega)
+	cv = log(omega/filterParam.omega_min)/(filterParam.num_octaves * log(filterParam.sweepbase));
+	cv = min(1, max(0, cv));
+
+function omega = cv2omega(filterParam, cv)
+	omega = filterParam.omega_min*filterParam.sweepbase .^ (filterParam.num_octaves*min(1, max(0, cv))); 
diff --git a/matlab/lfo/eval_lfo.asv b/matlab/lfo/eval_lfo.asv
new file mode 100755
index 0000000..0dcbf9d
--- /dev/null
+++ b/matlab/lfo/eval_lfo.asv
@@ -0,0 +1,34 @@
+function eval_lfo()
+
+L = 10000;
+fs = 44100;
+fstart = 40;
+fend = 400;
+df = (fend-fstart)/L;
+
+ip = 0 % phase of the first output sample in radians
+w = freq*pi / samplerate
+b1 = 2.0 * cos(w)
+
+% Init
+y1=sin(ip-w)
+y2=sin(ip-2*w)
+
+% Loop
+for n=1:L,
+    
+y0 = b1*y1 - y2
+y2 = y1
+y1 = y0
+
+
+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');
+
diff --git a/matlab/lfo/eval_lfo.m b/matlab/lfo/eval_lfo.m
new file mode 100755
index 0000000..27d9f17
--- /dev/null
+++ b/matlab/lfo/eval_lfo.m
@@ -0,0 +1,44 @@
+function eval_lfo()
+
+L = 10000;
+fs = 44100;
+fstart = 44;
+fend = 400;
+df = (fend-fstart)/L;
+
+ip = 0 % phase of the first output sample in radians
+
+% Init
+
+% Loop
+f = fstart;
+fchg = 1;
+
+y1=sin(ip-f*pi / fs);
+y2=sin(ip-2*f*pi / fs);
+
+for n=1:L,
+
+if fchg == 1,
+    b1 = 2.0 * cos(f*pi / fs);
+    fchg = 1000;
+end;
+
+y0 = b1*y1 - y2;
+y2 = y1;
+y1 = y0;
+
+fchg = fchg - 1;
+
+lfo(n) = y0;
+f = f + df;
+
+close all
+
+end;
+
+figure;
+plot(1:L, lfo); grid;
+
+wavwrite(0.5*lfo, fs, 'lfo.wav');
+
diff --git a/matlab/lfo/eval_lfo2.asv b/matlab/lfo/eval_lfo2.asv
new file mode 100755
index 0000000..15168cf
--- /dev/null
+++ b/matlab/lfo/eval_lfo2.asv
@@ -0,0 +1,47 @@
+function eval_lfo2()
+
+L = 80000;
+fs = 48000;
+fstart = 10;
+fend = 1000;
+df = (fend-fstart)/L;
+
+ip = 0.0 % phase of the first output sample in radians
+
+% Init
+
+% Loop
+f = fstart;
+fchg = 0;
+
+ylast = 0;
+a = 0.5;
+b = 2.0 * sin(f*pi / fs);
+y0 = a*cos(2*pi*ip);
+y1 = a*sin(2*pi*ip);
+
+for n=1:L,
+
+y0 = y0 - b*y1;
+y1 = y1 + b*y0;
+
+fchg = (mod(n, 1000) == 0);
+
+if fchg
+    b = 2.0 * sin(f*pi / fs);
+end
+
+ylast = y1;
+
+lfo(n) = y0;
+f = f + df;
+
+close all
+
+end;
+
+figure;
+plot(1:L, lfo); grid;
+
+wavwrite(0.5*lfo, fs, 'lfo.wav');
+
diff --git a/matlab/lfo/eval_lfo2.m b/matlab/lfo/eval_lfo2.m
new file mode 100755
index 0000000..b15e13b
--- /dev/null
+++ b/matlab/lfo/eval_lfo2.m
@@ -0,0 +1,47 @@
+function eval_lfo2()
+
+L = 80000;
+fs = 48000;
+fstart = 10;
+fend = 5000;
+df = (fend-fstart)/L;
+
+ip = 0.0 % phase of the first output sample in radians
+
+% Init
+
+% Loop
+f = fstart;
+fchg = 0;
+
+ylast = 0;
+a = 0.5;
+b = 2.0 * sin(f*pi / fs);
+y0 = a*cos(2*pi*ip);
+y1 = a*sin(2*pi*ip);
+
+for n=1:L,
+
+y0 = y0 - b*y1;
+y1 = y1 + b*y0;
+
+fchg = (mod(n, 1000) == 0);
+
+if fchg
+    b = 2.0 * sin(f*pi / fs);
+end
+
+ylast = y1;
+
+lfo(n) = y0;
+f = f + df;
+
+close all
+
+end;
+
+figure;
+plot(1:L, lfo); grid;
+
+wavwrite(0.5*lfo, fs, 'lfo.wav');
+
diff --git a/matlab/lfo/eval_smooth.asv b/matlab/lfo/eval_smooth.asv
new file mode 100755
index 0000000..bde255e
--- /dev/null
+++ b/matlab/lfo/eval_smooth.asv
@@ -0,0 +1,21 @@
+function eval_smooth(Ta)
+
+fs = 10000;
+
+b = 1;
+if (abs(Ta) > 0)
+    b = 1/(abs(Ta)*fs);
+end;
+
+a = 1 - b;
+
+x = [ zeros(1, 1000) ones(1, 1000) zeros(1, 1000)];
+
+y = 0;
+for i=1:length(x)
+    if (Ta < 0)
+    y = b*x(i) + a*y;
+    yn(i) = y;
+end;
+
+plot(yn); grid;
\ No newline at end of file
diff --git a/matlab/lfo/eval_smooth.m b/matlab/lfo/eval_smooth.m
new file mode 100755
index 0000000..a6999f3
--- /dev/null
+++ b/matlab/lfo/eval_smooth.m
@@ -0,0 +1,34 @@
+function eval_smooth(Ta)
+
+fs = 10000;
+
+b = 1;
+if (abs(Ta) > 0)
+    b = 1/(abs(Ta)*fs);
+end;
+
+a = 1 - b;
+
+x = [ zeros(1, 1000) ones(1, 1000) zeros(1, 1000)];
+N = length(x);
+
+y = 0;
+y1 = 0;
+for i=1:length(x)
+    y1 = b*x(i) + a*y1;
+%    if (Ta > 0)
+        y1n(i) = y1;
+%    else
+        y2n(i) = 2*x(i) - y1;
+%    end
+end;
+
+subplot(3, 1, 1)
+plot(1:N, x, 'r'); grid;
+xlabel('Original LFO Wellenform (Square)')
+subplot(3, 1, 2)
+plot(1:N, y1n, 'r'); grid;
+xlabel('Positive smoothed')
+subplot(3, 1, 3)
+plot(1:N, y2n, 'r'); grid;
+xlabel('Negative smoothed')
diff --git a/matlab/lfo/lfo.wav b/matlab/lfo/lfo.wav
new file mode 100755
index 0000000..8a6ba35
Binary files /dev/null and b/matlab/lfo/lfo.wav differ
diff --git a/matlab/limiter/eval_limiter.asv b/matlab/limiter/eval_limiter.asv
new file mode 100755
index 0000000..8c70d6c
--- /dev/null
+++ b/matlab/limiter/eval_limiter.asv
@@ -0,0 +1,38 @@
+function eval_limiter(name)
+
+Tr = 1;     % ms
+Tf = 100;  % ms
+input_gain = 2;
+output_gain = 0.5;
+
+limit_thresh = 0.25;
+
+[x, fs, nbits, opts] = wavread(name);
+x = x.*input_gain;
+
+N = fix(length(x)/);
+
+a_r = 1000/(Tr*fs);
+a_f = 1000/(Tf*fs);
+
+y = 0;
+for n=1:N
+    if (y < abs(x(n)))
+		y = y + a_r*(abs(x(n)) - y);
+	else
+		y = y + a_f*(abs(x(n)) - y);
+    end
+    env(n) = y;
+end;
+
+limiter_gain = limit_thresh./max(limit_thresh, env)';
+x_limited = x(1:N).*limiter_gain.*output_gain;
+
+wavwrite(x_limited,fs,'limiter_out.wav');
+
+subplot(2, 1, 1)
+plot(0:N-1, x(1:N), 0:N-1, env, 'r-'); grid;
+subplot(2, 1, 2)
+plot(0:N-1, x_limited, 0:N-1, limiter_gain, 'r-'); grid;
+
+
diff --git a/matlab/limiter/eval_limiter.m b/matlab/limiter/eval_limiter.m
new file mode 100755
index 0000000..8c92020
--- /dev/null
+++ b/matlab/limiter/eval_limiter.m
@@ -0,0 +1,38 @@
+function eval_limiter(name)
+
+Tr = 1;     % ms
+Tf = 100;  % ms
+input_gain = 2;
+output_gain = 0.5;
+
+limit_thresh = 0.25;
+
+[x, fs, nbits, opts] = wavread(name);
+x = x.*input_gain;
+
+N = fix(length(x));
+
+a_r = 1000/(Tr*fs);
+a_f = 1000/(Tf*fs);
+
+y = 0;
+for n=1:N
+    if (y < abs(x(n)))
+		y = y + a_r*(abs(x(n)) - y);
+	else
+		y = y + a_f*(abs(x(n)) - y);
+    end
+    env(n) = y;
+end;
+
+limiter_gain = limit_thresh./max(limit_thresh, env)';
+x_limited = x(1:N).*limiter_gain.*output_gain;
+
+wavwrite(x_limited,fs,'limiter_out.wav');
+
+subplot(2, 1, 1)
+plot(0:N-1, x(1:N), 0:N-1, env, 'r-'); grid;
+subplot(2, 1, 2)
+plot(0:N-1, x_limited, 0:N-1, limiter_gain, 'r-'); grid;
+
+
diff --git a/matlab/limiter/limit_test1.wav b/matlab/limiter/limit_test1.wav
new file mode 100755
index 0000000..cadc0e7
Binary files /dev/null and b/matlab/limiter/limit_test1.wav differ
diff --git a/matlab/limiter/limit_test2.wav b/matlab/limiter/limit_test2.wav
new file mode 100755
index 0000000..9fb11e4
Binary files /dev/null and b/matlab/limiter/limit_test2.wav differ
diff --git a/matlab/limiter/limiter_out.wav b/matlab/limiter/limiter_out.wav
new file mode 100755
index 0000000..747188d
Binary files /dev/null and b/matlab/limiter/limiter_out.wav differ
diff --git a/matlab/midi/eval_midi_sync.asv b/matlab/midi/eval_midi_sync.asv
new file mode 100755
index 0000000..a6da3c9
--- /dev/null
+++ b/matlab/midi/eval_midi_sync.asv
@@ -0,0 +1,77 @@
+function eval_midi_sync(pe, fe)
+N = 5000;
+t = (0:N-1)/N;
+
+x1 = mod(20*t, 1) >=0.5;
+x2 = mod((20+fe)*t+pe, 1) >=0.5;
+
+dx1 = 50/N
+dx2 = (50+fe)/N
+
+qi1 = 0;
+qi2 = 0;
+xi1 = 0;
+xi2 = 0;
+xx1 = 0;
+xx2 = 0;
+
+lead = 0;
+lag = 0;
+ddx = 0;
+for i=1:N,
+%[qi1, qi2, xi1, xi2] = pd(x1(i), x2(i), qi1, qi2, xi1, xi2);
+[qi1, qi2, xi1, xi2] = pd(xx1 >=0.5, xx2 >=0.5, qi1, qi2, xi1, xi2);
+err = (qi1-qi2);
+lead = 0.*err;
+lag = lag + 0.000001*err;
+ddx = (lead+lag);
+xx1 = mod(xx1 + dx1, 1);
+xx2 = mod(xx2 + dx2 + ddx, 1);
+
+
+q1(i) = qi1;
+q2(i) = qi2;
+
+end;
+
+dx2 = dx2 + ddx
+dx1 = dx1
+
+subplot(3, 1, 1)
+plot(t, x1, t, x2); grid; axis([0 1 -0.2 1.2])
+subplot(3, 1, 2)
+plot(t, q1); grid; axis([0 1 -0.2 1.2])
+subplot(3, 1, 3)
+plot(t, q2); grid; axis([0 1 -0.2 1.2])
+
+
+function [Q1, Q2, xo1, xo2] = pd(x1, x2, qi1, qi2, xi1, xi2)
+xo1 = xi1;
+xo2 = xi2;
+qo1 = qi1;
+qo2 = qi2;
+for i=1:length(x1)
+   % Q1 
+   % Detect edge
+   if (x1(i) - xo1) > 0.5
+       qo1 = 1;
+   end;
+   xo1 = x1(i);
+   
+   % Q2 
+   % Detect edge
+   if (x2(i) - xo2) > 0.5
+       qo2 = 1;
+   end;   
+   xo2 = x2(i);
+   
+   % Asynchronous reset
+   if (qo1 * qo2) > 0.5
+       qo1 = 0;
+       qo2 = 0;
+   end;
+
+   Q1(i) = qo1;
+   Q2(i) = qo2;
+   
+end;
\ No newline at end of file
diff --git a/matlab/midi/eval_midi_sync.m b/matlab/midi/eval_midi_sync.m
new file mode 100755
index 0000000..fa4234f
--- /dev/null
+++ b/matlab/midi/eval_midi_sync.m
@@ -0,0 +1,75 @@
+function eval_midi_sync(pe, fe)
+N = 5000;
+t = (0:N-1)/N;
+
+dx1 = 40/N
+dx2 = (40+fe)/N
+
+qi1 = 0;
+qi2 = 0;
+xi1 = 0;
+xi2 = 0;
+xx1 = 0;
+xx2 = pe;
+
+lead = 0;
+lag = 0;
+ddx = 0;
+for i=1:N,
+    
+    [qi1, qi2, xi1, xi2] = pd(xx1 >=0.5, xx2 >=0.5, qi1, qi2, xi1, xi2);
+    err = (qi1-qi2);
+    lead = 0.002*err;
+    lag = lag + 0.000005*err;
+    ddx(i) = (lead+lag);
+    xx1 = mod(xx1 + dx1, 1);
+    xx2 = mod(xx2 + dx2 + ddx(i), 1);
+
+
+    q1(i) = qi1;
+    q2(i) = qi2;
+    x1(i) = xx1  >=0.5;
+    x2(i) = xx2  >=0.5;
+
+end;
+
+subplot(4, 1, 1)
+plot(t, x1, t, x2); grid; axis([0 1 -0.2 1.2])
+subplot(4, 1, 2)
+plot(t, q1); grid; axis([0 1 -0.2 1.2])
+subplot(4, 1, 3)
+plot(t, q2); grid; axis([0 1 -0.2 1.2])
+subplot(4, 1, 4)
+plot(t, ddx); grid;
+
+
+function [Q1, Q2, xo1, xo2] = pd(x1, x2, qi1, qi2, xi1, xi2)
+xo1 = xi1;
+xo2 = xi2;
+qo1 = qi1;
+qo2 = qi2;
+for i=1:length(x1)
+   % Q1 
+   % Detect edge
+   if (x1(i) - xo1) > 0.5
+       qo1 = 1;
+   end;
+   xo1 = x1(i);
+   
+   % Q2 
+   % Detect edge
+   if (x2(i) - xo2) > 0.5
+       qo2 = 1;
+   end;   
+   xo2 = x2(i);
+   
+   % Asynchronous reset
+   if (qo1 * qo2) > 0.5
+       qo1 = 0;
+       qo2 = 0;
+   end;
+
+   Q1(i) = qo1;
+   Q2(i) = qo2;
+   
+end;
\ No newline at end of file
diff --git a/matlab/osc/blep/eval_minblep.asv b/matlab/osc/blep/eval_minblep.asv
new file mode 100755
index 0000000..f27af9e
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blep/eval_minblep.m b/matlab/osc/blep/eval_minblep.m
new file mode 100755
index 0000000..16922f0
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blep/eval_minblep_pre.asv b/matlab/osc/blep/eval_minblep_pre.asv
new file mode 100755
index 0000000..72fb350
--- /dev/null
+++ b/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;
\ No newline at end of file
diff --git a/matlab/osc/blep/eval_minblep_pre.m b/matlab/osc/blep/eval_minblep_pre.m
new file mode 100755
index 0000000..fc4bb3f
--- /dev/null
+++ b/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;
\ No newline at end of file
diff --git a/matlab/osc/blip/blep.wav b/matlab/osc/blip/blep.wav
new file mode 100755
index 0000000..ca8ef57
Binary files /dev/null and b/matlab/osc/blip/blep.wav differ
diff --git a/matlab/osc/blip/blit.wav b/matlab/osc/blip/blit.wav
new file mode 100755
index 0000000..5475c84
Binary files /dev/null and b/matlab/osc/blip/blit.wav differ
diff --git a/matlab/osc/blip/c3.wav b/matlab/osc/blip/c3.wav
new file mode 100755
index 0000000..9d0baee
Binary files /dev/null and b/matlab/osc/blip/c3.wav differ
diff --git a/matlab/osc/blip/c4.wav b/matlab/osc/blip/c4.wav
new file mode 100755
index 0000000..ad9faaa
Binary files /dev/null and b/matlab/osc/blip/c4.wav differ
diff --git a/matlab/osc/blip/eval_blep.asv b/matlab/osc/blip/eval_blep.asv
new file mode 100755
index 0000000..d429d5e
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blep.m b/matlab/osc/blip/eval_blep.m
new file mode 100755
index 0000000..ac48451
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blep_hs.asv b/matlab/osc/blip/eval_blep_hs.asv
new file mode 100755
index 0000000..02ef73c
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blep_hs.m b/matlab/osc/blip/eval_blep_hs.m
new file mode 100755
index 0000000..374063e
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blit.asv b/matlab/osc/blip/eval_blit.asv
new file mode 100755
index 0000000..f6acdd0
--- /dev/null
+++ b/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');
+
diff --git a/matlab/osc/blip/eval_blit.bak.m b/matlab/osc/blip/eval_blit.bak.m
new file mode 100755
index 0000000..bd1e326
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blit.m b/matlab/osc/blip/eval_blit.m
new file mode 100755
index 0000000..fa02623
--- /dev/null
+++ b/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');
+
diff --git a/matlab/osc/blip/eval_blit2.m b/matlab/osc/blip/eval_blit2.m
new file mode 100755
index 0000000..cf44271
--- /dev/null
+++ b/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
diff --git a/matlab/osc/blip/eval_blit3.m b/matlab/osc/blip/eval_blit3.m
new file mode 100755
index 0000000..ffddd4f
--- /dev/null
+++ b/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');
+
diff --git a/matlab/osc/blip/eval_blit4.asv b/matlab/osc/blip/eval_blit4.asv
new file mode 100755
index 0000000..4c2cc16
--- /dev/null
+++ b/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');
+
diff --git a/matlab/osc/blip/eval_blit4.m b/matlab/osc/blip/eval_blit4.m
new file mode 100755
index 0000000..e3b4582
--- /dev/null
+++ b/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');
+
diff --git a/matlab/osc/eval_hardsync.asv b/matlab/osc/eval_hardsync.asv
new file mode 100755
index 0000000..6f77f0b
--- /dev/null
+++ b/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
diff --git a/matlab/osc/eval_hardsync.m b/matlab/osc/eval_hardsync.m
new file mode 100755
index 0000000..61dd297
--- /dev/null
+++ b/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
diff --git a/matlab/osc/eval_harmonics.asv b/matlab/osc/eval_harmonics.asv
new file mode 100755
index 0000000..d035720
--- /dev/null
+++ b/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, 
\ No newline at end of file
diff --git a/matlab/osc/eval_harmonics.m b/matlab/osc/eval_harmonics.m
new file mode 100755
index 0000000..0e23438
--- /dev/null
+++ b/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);
diff --git a/matlab/osc/saw.wav b/matlab/osc/saw.wav
new file mode 100755
index 0000000..203f47c
Binary files /dev/null and b/matlab/osc/saw.wav differ
diff --git a/matlab/osc/sqr.wav b/matlab/osc/sqr.wav
new file mode 100755
index 0000000..91f0c31
Binary files /dev/null and b/matlab/osc/sqr.wav differ
diff --git a/matlab/osc/tri.wav b/matlab/osc/tri.wav
new file mode 100755
index 0000000..92b1c4d
Binary files /dev/null and b/matlab/osc/tri.wav differ
diff --git a/matlab/params/paramscale.asv b/matlab/params/paramscale.asv
new file mode 100755
index 0000000..5b05d95
--- /dev/null
+++ b/matlab/params/paramscale.asv
@@ -0,0 +1,66 @@
+function paramscale(base, kexp, scenter, pcenter, pmin, pmax)
+
+N = 1000;
+slider = (0:N)/N;
+
+p = toParam(base, kexp, scenter, pcenter, pmin, pmax, slider);
+s = toSlider(base, kexp, scenter, pcenter, pmin, pmax, p);
+
+subplot (2, 1, 1)
+plot (0:N, p); grid; xlabel('param');
+subplot (2, 1, 2)
+plot (0:N, s); grid; xlabel('slider');
+
+function param = toParam(base, kexp, scenter, pcenter, pmin, pmax, slider)
+	
+%    pcenter = pmax*scenter + pmin*(1-scenter)
+
+    for i=1:length(slider),
+    	s = min(1, max(0, slider(i)));
+        
+        if (s < scenter)
+            if (base == 1)
+                p = pcenter-(pcenter-pmin)*(scenter-s)/scenter;
+                % p = pcenter-(pcenter-pmin)*(scenter-s)/scenter;
+            else
+                p = pcenter-(pcenter-pmin)*(base^(kexp*(scenter-s)/scenter)-1)/(base^kexp-1);
+                % p = pcenter-(pcenter-pmin)*(pow(base,kexp*(scenter-s)/scenter)-1)/(pow(base,kexp)-1);
+            end
+        else
+            if (base == 1)
+                p = pcenter+(pmax-pcenter)*(-scenter+s)/(1-scenter);
+                % p = pcenter+(pmax-pcenter)*(-scenter+s)/(1-scenter);
+            else
+                p = pcenter+(pmax-pcenter)*(base^(kexp*(-scenter+s)/(1-scenter))-1)/(base^kexp-1);
+                % p = pcenter+(pmax-pcenter)*(pow(base,kexp*(-scenter+s)/(1-scenter))-1)/(pow(base,kexp)-1);
+            end;
+        end
+        param(i) = p;
+    end
+    
+function slider = toSlider(base, kexp, scenter, pcenter, pmin, pmax, param)
+
+%pcenter = pmax*scenter + pmin*(1-scenter);
+
+    for i=1:length(param),
+        p = param(i);
+        
+        if (p < pcenter)
+            if (base == 1)
+                s = (-pmin+p)*scenter/(pcenter-pmin);
+                % s = (-pmin+p)*scenter/(pcenter-pmin);
+            else
+                s = scenter*(log(base)*kexp-log(-(-pcenter*base^kexp+pmin+p*base^kexp-p)/(pcenter-pmin)))/log(base)/kexp;
+                % s = scenter*(log(base)*kexp-log(-(-pcenter*pow(base,kexp)+pmin+p*pow(base,kexp)-p)/(pcenter-pmin)))/log(base)/kexp;
+            end
+        else
+            if (base == 1)
+                s = (pcenter-scenter*pmax-p+scenter*p)/(-pmax+pcenter);
+                % s = (pcenter-scenter*pmax-p+scenter*p)/(-pmax+pcenter);
+            else
+                s = (kexp*log(base)*scenter+log(-(-pcenter*base^kexp+pmax+p*base^kexp-p)/(-pmax+pcenter))-log(-(-pcenter*base^kexp+pmax+p*base^kexp-p)/(-pmax+pcenter))*scenter)/log(base)/kexp;
+                % s = (kexp*log(base)*scenter+log(-(-pcenter*pow(base,kexp)+pmax+p*pow(base,kexp)-p)/(-pmax+pcenter))-log(-(-pcenter*pow(base,kexp)+pmax+p*pow(base,kexp)-p)/(-pmax+pcenter))*scenter)/log(base)/kexp;
+            end
+        end
+        slider(i) = s;
+    end
diff --git a/matlab/params/paramscale.m b/matlab/params/paramscale.m
new file mode 100755
index 0000000..162649d
--- /dev/null
+++ b/matlab/params/paramscale.m
@@ -0,0 +1,13 @@
+function paramscale(base, kexp, scenter, pcenter, pmin, pmax)
+
+N = 1000;
+slider = (0:N)/N;
+
+p = toParam(base, kexp, scenter, pcenter, pmin, pmax, slider);
+s = toSlider(base, kexp, scenter, pcenter, pmin, pmax, p);
+
+subplot (2, 1, 1)
+plot (slider, p); grid; xlabel('toParam(s)');
+subplot (2, 1, 2)
+plot (slider, s); grid; xlabel('toSlider(p)');
+
diff --git a/matlab/params/toParam.m b/matlab/params/toParam.m
new file mode 100755
index 0000000..ba76cbe
--- /dev/null
+++ b/matlab/params/toParam.m
@@ -0,0 +1,28 @@
+function param = toParam(base, kexp, scenter, pcenter, pmin, pmax, slider)
+	
+	if ((pcenter < pmin) || (pcenter > pmax))
+	    pcenter = pmax*scenter + pmin*(1-scenter);
+    end
+
+    for i=1:length(slider),
+    	s = min(1, max(0, slider(i)));
+        
+        if (s < scenter)
+            if (base == 1)
+                p = pcenter-(pcenter-pmin)*(scenter-s)/scenter;
+                % p = pcenter-(pcenter-pmin)*(scenter-s)/scenter;
+            else
+                p = pcenter-(pcenter-pmin)*(base^(kexp*(scenter-s)/scenter)-1)/(base^kexp-1);
+                % p = pcenter-(pcenter-pmin)*(pow(base,kexp*(scenter-s)/scenter)-1)/(pow(base,kexp)-1);
+            end
+        else
+            if (base == 1)
+                p = pcenter+(pmax-pcenter)*(-scenter+s)/(1-scenter);
+                % p = pcenter+(pmax-pcenter)*(-scenter+s)/(1-scenter);
+            else
+                p = pcenter+(pmax-pcenter)*(base^(kexp*(-scenter+s)/(1-scenter))-1)/(base^kexp-1);
+                % p = pcenter+(pmax-pcenter)*(pow(base,kexp*(-scenter+s)/(1-scenter))-1)/(pow(base,kexp)-1);
+            end;
+        end
+        param(i) = p;
+    end
diff --git a/matlab/params/toSlider.m b/matlab/params/toSlider.m
new file mode 100755
index 0000000..d031aa0
--- /dev/null
+++ b/matlab/params/toSlider.m
@@ -0,0 +1,28 @@
+function slider = toSlider(base, kexp, scenter, pcenter, pmin, pmax, param)
+
+	if ((pcenter < pmin) || (pcenter > pmax))
+	    pcenter = pmax*scenter + pmin*(1-scenter);
+    end
+
+    for i=1:length(param),
+        p = param(i);
+        
+        if (p < pcenter)
+            if (base == 1)
+                s = (-pmin+p)*scenter/(pcenter-pmin);
+                % s = (-pmin+p)*scenter/(pcenter-pmin);
+            else
+                s = scenter*(log(base)*kexp-log(-(-pcenter*base^kexp+pmin+p*base^kexp-p)/(pcenter-pmin)))/log(base)/kexp;
+                % s = scenter*(log(base)*kexp-log(-(-pcenter*pow(base,kexp)+pmin+p*pow(base,kexp)-p)/(pcenter-pmin)))/log(base)/kexp;
+            end
+        else
+            if (base == 1)
+                s = (pcenter-scenter*pmax-p+scenter*p)/(-pmax+pcenter);
+                % s = (pcenter-scenter*pmax-p+scenter*p)/(-pmax+pcenter);
+            else
+                s = (kexp*log(base)*scenter+log(-(-pcenter*base^kexp+pmax+p*base^kexp-p)/(-pmax+pcenter))-log(-(-pcenter*base^kexp+pmax+p*base^kexp-p)/(-pmax+pcenter))*scenter)/log(base)/kexp;
+                % s = (kexp*log(base)*scenter+log(-(-pcenter*pow(base,kexp)+pmax+p*pow(base,kexp)-p)/(-pmax+pcenter))-log(-(-pcenter*pow(base,kexp)+pmax+p*pow(base,kexp)-p)/(-pmax+pcenter))*scenter)/log(base)/kexp;
+            end
+        end
+        slider(i) = s;
+    end