1import numpy 2import pylab 3 4def dumb(f): 5 return f * numpy.pi 6 7def pade(f): 8 return f * -0.90585 / (-0.28833 + f * f) 9 10def poly3taylor(f): 11 fsq = f * f 12 r = 0.3333333 * numpy.pi ** 3 13 r *= fsq 14 r += numpy.pi 15 r *= f 16 return r 17 18def poly3gradient(f, a=3.736e-01): 19 fsq = f * f 20 r = a * numpy.pi ** 3 21 r *= fsq 22 r += numpy.pi 23 r *= f 24 return r 25 26 27def poly5mdsp(f, a=3.1755e-01, b=2.033e-01): 28 fsq = f * f 29 r = b * numpy.pi ** 5 30 r *= fsq 31 r += a * numpy.pi ** 3 32 r *= fsq 33 r += numpy.pi 34 r *= f 35 return r 36 37 38def poly5gradient(f, a=3.260e-01, b=1.823e-01): 39 f = f * numpy.pi 40 fsq = f * f 41 r = b 42 r *= fsq 43 r += a 44 r *= fsq 45 r += 1.0 46 r *= f 47 return r 48 49 50def poly11mdsp(f): 51 fsq = f * f 52 r = 9.5168091e-03 * numpy.pi ** 11 53 r *= fsq 54 r += 2.900525e-03 * numpy.pi ** 9 55 r *= fsq 56 r += 5.33740603e-02 * numpy.pi **7 57 r *= fsq 58 r += 1.333923995e-01 * numpy.pi **5 59 r *= fsq 60 r += 3.333314036e-01 * numpy.pi **3 61 r *= fsq 62 r += numpy.pi 63 r *= f 64 return r 65 66 67def compute_filter_settings(cutoff, resonance): 68 g = numpy.tan(numpy.pi * cutoff) + resonance * 0 69 r = 1.0 / resonance + cutoff * 0 70 h = 1 / (1 + r * g + g * g) 71 return g, r, h 72 73 74def evaluate(groundtruth_f, approximate_g): 75 approximate_f = numpy.arctan(approximate_g) / numpy.pi 76 return numpy.log2(approximate_f / groundtruth_f) * 1200.0 77 78 79f = numpy.exp(numpy.linspace(numpy.log(16), numpy.log(10000), 1000.0)) 80f /= 48000.0 81g, _, _ = compute_filter_settings(f, 0.5) 82approximations = [pade, poly3gradient, poly5mdsp, poly5gradient, poly11mdsp] 83# 84# a_ = numpy.linspace(3.259e-01, 3.261e-01, 100) 85# b_ = numpy.linspace(1.822e-01, 1.823e-01, 100) 86# error = numpy.zeros((100, 100)) 87# best = 1e8 88# arg_best = None 89# for i, a in enumerate(a_): 90# for j, b in enumerate(b_): 91# error[i, j] = (evaluate(f, fast0(f, a, b)) ** 2).sum() 92# if error[i, j] < best: 93# best = error[i, j] 94# arg_best = (a, b) 95# 96# print arg_best 97# pylab.plot(coefficient, error) 98# pylab.show() 99 100pylab.figure(figsize=(15,10)) 101for i, a in enumerate(approximations): 102 n = len(approximations) 103 # pylab.subplot(n * 100 + 10 + i + 1) 104 pylab.plot(f * 48000, evaluate(f, a(f))) 105 pylab.xlabel('Hz') 106 pylab.ylabel('$\delta$ cents') 107 108pylab.legend(map(lambda x: x.__name__, approximations)) 109pylab.tight_layout() 110#pylab.savefig('plot.pdf') 111pylab.show() 112