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