diff --git a/lab/ladder/sawtest.py b/lab/ladder/sawtest.py
new file mode 100644
index 0000000..1224ade
--- /dev/null
+++ b/lab/ladder/sawtest.py
@@ -0,0 +1,383 @@
+import argparse
+import struct
+import numpy as np
+from math import *
+
+samplerate = 44100 * 1
+
+# generate a saw matching the parameters of "fun with sigmoids"
+def sawperiod():
+	result = []
+	f0 = 63. * 2 * pi / samplerate
+	for i in range(samplerate / 63):
+		y = 0
+		x0 = i * f0
+		for partial in range(1, 351):
+			gain = 1.0 / partial
+			if partial > 300:
+				gain *= (351 - partial) * (1.0 / (351 - 300))
+			y += gain * sin(partial * x0)
+		result.append(y)
+	return result
+
+def saw(n):
+	return sawperiod() * (n / (samplerate / 63))
+
+def sinsweep(n):
+	fmax = 22000 * pi * 2 / samplerate
+	scale = fmax * .5 / n
+	return [sin(i * i * scale) for i in range(n)]
+
+def sweep(n):
+	n2 = n / 2
+	lamin = log(20 * 2 * pi / samplerate)
+	lamax = log(14000 * 2 * pi / samplerate)
+	result = []
+	slope = (lamax - lamin) / n2
+	for i in range(n2):
+		a = exp(lamin + slope * i)
+		result.append(a)
+	return result + result[-1::-1]
+
+# Based on code by mystran (Teemu Voipio)
+# See http://www.kvraudio.com/forum/viewtopic.php?t=349859
+
+def tanhXdx(x):
+	a = x * x
+	return ((a + 105.0)* a + 945.0) / ((15.0 * a + 420.0) * a + 945.0)
+
+def tpt_nl(xs, aas, k):
+	z1 = 0
+	s = [0] * 4
+	r = k
+	result = []
+	for i in range(len(xs)):
+		xin = xs[i]
+		a = aas[i]
+		f = tan(a * .5)
+		ih = .5 * (xin + z1)
+		z1 = xin
+
+		t0 = tanhXdx(ih - r * s[3])
+		t1 = tanhXdx(s[0])
+		t2 = tanhXdx(s[1])
+		t3 = tanhXdx(s[2])
+		t4 = tanhXdx(s[3])
+		g0 = 1 / (1 + f * t1)
+		g1 = 1 / (1 + f * t2)
+		g2 = 1 / (1 + f * t3)
+		g3 = 1 / (1 + f * t4)
+		f3 = f * t3 * g3
+		f2 = f * t2 * g2 * f3
+		f1 = f * t1 * g1 * f2
+		f0 = f * t0 * g0 * f1
+		y3 = (g3*s[3] + f3*g2*s[2] + f2*g1*s[1] + f1*g0*s[0] + f0*xin)/(1 + r * f0)
+		xx = t0 * (xin - r * y3)
+		y0 = t1 * g0 * (s[0] + f * xx)
+		y1 = t2 * g1 * (s[1] + f * y0)
+		y2 = t3 * g2 * (s[2] + f * y1)
+		s[0] += 2 * f * (xx - y0)
+		s[1] += 2 * f * (y0 - y1)
+		s[2] += 2 * f * (y1 - y2)
+		s[3] += 2 * f * (y2 - t4 * y3)
+		result.append(y3)
+	return result
+
+def antti_nl(xs, aas, k, tanhfunc = tanh):
+	y0 = 0
+	y1 = 0
+	y2 = 0
+	y3 = 0
+	ty0 = 0
+	ty1 = 0
+	ty2 = 0
+	ty3 = 0
+	yy = 0
+	result = []
+	for i in range(len(xs)):
+		xin = xs[i]
+		a = aas[i]
+		tx = tanhfunc(xin - k * (y3 + yy) * .5)
+		y0 += a * (tx - ty0)
+		yy = y3
+		ty0 = tanhfunc(y0)
+		y1 += a * (ty0 - ty1)
+		ty1 = tanhfunc(y1)
+		y2 += a * (ty1 - ty2)
+		ty2 = tanhfunc(y2)
+		y3 += a * (ty2 - ty3)
+		ty3 = tanhfunc(y3)
+		x = 0
+		result.append(y3)
+	return result
+
+def expm_series(A, n = 16):
+	B = np.identity(len(A))
+	C = np.identity(len(A))
+	for i in range(1, n):
+		C = A * C / i
+		B = B + C
+	return B
+
+def expm_hyb(A, n1 = 8, n2 = 8):
+	A = expm_series(A / (1 << n2), n1)
+	for i in range(n2):
+		A = A * A
+	return A	
+
+def mkjacobian(a, k):
+	return np.matrix([[0, 0, 0, 0, 0],
+		[a, -a, 0, 0, -k * a],
+		[0, a, -a, 0, 0],
+		[0, 0, a, -a, 0],
+		[0, 0, 0, a, -a]])
+
+def expm_nl(xs, aas, k, tanhfunc = tanh, every = 64):
+	kk = k
+	y = np.zeros([4, 1])
+	ty = y
+	result = []
+	for i in range(len(xs)):
+		if i % every == 0:
+			a = aas[i]
+			A = expm_hyb(mkjacobian(a, k))
+
+			B = A[1:, 0]
+			A = A[1:, 1:]
+			AM = A - np.identity(4)
+
+			for j in range(4):
+				AM[j, 3] += B[j, 0] * kk
+
+		x = xs[i]
+		tx = tanhfunc(x - kk * y[3, 0])
+		y += B * tx + AM * ty
+		ty = np.matrix([tanhfunc(x[0]) for x in y]).T
+
+		result.append(y[3, 0])
+	return result
+
+def invsqrt(x):
+	return x / sqrt(1 + x ** 2)
+
+def clip(x):
+	return min(1, max(-1, x))
+
+fir = map(float, '''-0.00000152158394097619
+-0.00000932875737718674
+-0.00001008290833020705
+ 0.00000728628960094510
+ 0.00002272556429291851
+-0.00000017886444625648
+-0.00004038828866646377
+-0.00002127235603321839
+ 0.00005623314602326363
+ 0.00006233931357689778
+-0.00005884569707596701
+-0.00012410666372671377
+ 0.00003231781361429016
+ 0.00019973378481126099
+ 0.00004097428652472217
+-0.00027125558280978066
+-0.00017513777576291543
+ 0.00030833229435342778
+ 0.00037363561338823163
+-0.00027028256701905416
+-0.00062153272838533420
+ 0.00011242294792251808
+ 0.00087909017443692499
+ 0.00020309303892565388
+-0.00107928469499717137
+-0.00069305675670084180
+ 0.00113148992142813524
+ 0.00133794323705832747
+-0.00093286682198106608
+-0.00206892788224323455
+ 0.00038777937772768273
+ 0.00276141532345176117
+ 0.00056670689248167661
+-0.00323852878331680298
+-0.00193259906597622396
+ 0.00328731491441583683
+ 0.00362542589345449945
+-0.00268785554870503091
+-0.00545702740033716938
+ 0.00125317559909794512
+ 0.00713205061850966104
+ 0.00112492149537589880
+-0.00826148096867550946
+-0.00443085208574918715
+ 0.00839336065985300112
+ 0.00848839486761648020
+-0.00705662837862476577
+-0.01293781936476604867
+ 0.00380913247729417143
+ 0.01722862947876550518
+ 0.00172513205576642695
+-0.02062091118580804822
+-0.00985374287540894019
+ 0.02217056516274141381
+ 0.02089533076058044253
+-0.02062760756367006815
+-0.03546365160613443313
+ 0.01401630243169032369
+ 0.05541063386809867708
+ 0.00212045427486363437
+-0.08794052708207862612
+-0.04526389505656687462
+ 0.18221762668014460096
+ 0.41075960732889965632
+ 0.41075960732889965632
+ 0.18221762668014460096
+-0.04526389505656687462
+-0.08794052708207862612
+ 0.00212045427486363437
+ 0.05541063386809867708
+ 0.01401630243169032369
+-0.03546365160613443313
+-0.02062760756367006815
+ 0.02089533076058044253
+ 0.02217056516274141381
+-0.00985374287540894019
+-0.02062091118580804822
+ 0.00172513205576642695
+ 0.01722862947876550518
+ 0.00380913247729417143
+-0.01293781936476604867
+-0.00705662837862476577
+ 0.00848839486761648020
+ 0.00839336065985300112
+-0.00443085208574918715
+-0.00826148096867550946
+ 0.00112492149537589880
+ 0.00713205061850966104
+ 0.00125317559909794512
+-0.00545702740033716938
+-0.00268785554870503091
+ 0.00362542589345449945
+ 0.00328731491441583683
+-0.00193259906597622396
+-0.00323852878331680298
+ 0.00056670689248167661
+ 0.00276141532345176117
+ 0.00038777937772768273
+-0.00206892788224323455
+-0.00093286682198106608
+ 0.00133794323705832747
+ 0.00113148992142813524
+-0.00069305675670084180
+-0.00107928469499717137
+ 0.00020309303892565388
+ 0.00087909017443692499
+ 0.00011242294792251808
+-0.00062153272838533420
+-0.00027028256701905416
+ 0.00037363561338823163
+ 0.00030833229435342778
+-0.00017513777576291543
+-0.00027125558280978066
+ 0.00004097428652472217
+ 0.00019973378481126099
+ 0.00003231781361429016
+-0.00012410666372671377
+-0.00005884569707596701
+ 0.00006233931357689778
+ 0.00005623314602326363
+-0.00002127235603321839
+-0.00004038828866646377
+-0.00000017886444625648
+ 0.00002272556429291851
+ 0.00000728628960094510
+-0.00001008290833020705
+-0.00000932875737718674
+-0.00000152158394097619
+'''.split())
+
+fir_n = len(fir)
+
+# would probably be faster to use numpy.convolve and slice, but oh well
+def downsample(data):
+	result = []
+	buf = [0] * fir_n
+	i = 0
+	for x in data:
+		buf[i] = x
+		i = (i + 1) & (fir_n - 1)
+		if (i & 1) == 0:
+			y = 0
+			for j in range(len(fir)):
+				y += fir[j] * buf[(i + j) & (fir_n - 1)]
+			result.append(y)
+	return result
+
+def wavwrite(seq, fn, sr = 44100):
+  f = file(fn, 'wb')
+  n_samples = len(seq)
+  f.write(struct.pack('<4sI4s4sIHHIIHH4sI',
+        'RIFF',
+        36 + 2 * n_samples,
+        'WAVE',
+        'fmt ',
+        16,
+        1, 1,
+        sr,
+        2 * sr,
+        2, 16,
+        'data',
+        2 * n_samples))
+  for x in seq:
+    f.write(struct.pack('<h', min(32767, max(-32767, int(16384 * x)))))
+
+def main():
+	parser = argparse.ArgumentParser()
+	parser.add_argument("--oversample", help="oversample factor")
+	parser.add_argument("--signal", help="saw or sinsweep")
+	parser.add_argument("--k", help="resonance, 4 = oscillate")
+	parser.add_argument("--filter", help="tpt, antti, or expm")
+	parser.add_argument("--tanhfunc", help="tanh or invsqrt")
+	parser.add_argument("--cutoff")
+	parser.add_argument("--gain", help="gain in dB")
+	parser.add_argument("--ogain", help="outputgain in dB")
+	parser.add_argument("--out", help="output wav file")
+	args = parser.parse_args()
+	oversample = 1
+	if args.oversample:
+		oversample = int(args.oversample)
+	k = 0
+	if args.k:
+		k = float(args.k)
+	signal = 'saw'
+	if args.signal:
+		signal = args.signal
+	global samplerate
+	samplerate *= oversample
+	n = 6 * samplerate
+	if signal == 'saw':
+		input = saw(n)
+		aas = sweep(n)
+	elif signal == 'sinsweep':
+		input = sinsweep(n)
+		aas = [1000./samplerate * 2 * pi] * n
+	gain = 1
+	if args.gain:
+		gain = 10 ** (float(args.gain)/20)
+	input = [y * gain for y in input]
+	tanhfunc = tanh
+	if args.tanhfunc == 'invsqrt':
+		tanhfunc = invsqrt
+	if args.filter == 'tpt':
+		result = tpt_nl(input, aas, k)
+	elif args.filter == 'antti':
+		result = antti_nl(input, aas, k, tanhfunc = tanhfunc)
+	else:
+		result = expm_nl(input, aas, k, tanhfunc = tanhfunc)
+	while oversample > 1:
+		result = downsample(result)
+		oversample /= 2
+	ogain = 1
+	if args.ogain:
+		ogain = 10 ** (float(args.ogain)/20)
+	if args.out:
+		wavwrite(result, args.out)
+
+main()