From 75902125d5e1659a35b33df735c3af501fc2a0ae Mon Sep 17 00:00:00 2001
From: Raph Levien <raph@google.com>
Date: Mon, 8 Jul 2013 00:55:31 -0700
Subject: [PATCH] Test program for comparing ladder filter quality

The sawtest.py script generates audio samples for various filter
implementations (including matrix exponential and TPT) with various
input signals and settings. It should be useful for comparing the
quality of the various implementations, even though it's pretty slow.
---
 lab/ladder/sawtest.py | 383 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 383 insertions(+)
 create mode 100644 lab/ladder/sawtest.py

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()