"""
Model MCE 4 pole butterworth filters.
"""

from pylab import *

class MCEButterworth:
    def __init__(self, params):
        self.params = params
        self.accum = zeros(4)
    
    def spectrum(self, f):
        K = 1./2**14
        scalars = [K, K, K, K, 1., 1.]
        b11, b12, b21, b22, k1, k2 = [s*p for s,p in zip(scalars, self.params)]
        z = exp(-1j*f)
        H = (1. + z)**4 / (1. - b11*z + b12*z**2) / (1. - b21*z + b22*z**2)
        return H  / 2**(k1+k2)

    def gain(self):
        return self.spectrum(0)


# wb rca fltr_coeff b11 b12 b21 b22 k1 k2

coeffs = {
    'type1': [ 32092, 15750, 31238, 14895, 0, 11],
    'type2': [ 32295, 15915, 32568, 16188, 3, 14],
}

f0 = 15.4e3
fs = f0 / 150

n = 5000
f = f0/2 * arange(n) / n

b = MCEButterworth(coeffs['type2'])
y = b.spectrum(2.*pi * f / f0)

f_wrapped = abs(((f + fs/2) % fs) - fs/2)

subplot(211)
Y = (abs(y) / b.gain())
plot(f, Y, 'k')
plot(f_wrapped[f>fs/2], Y[f>fs/2], 'k', ls='dashed')
xlim(0., fs/2)
ylabel('Gain')

subplot(212)
plot(f, unwrap(angle(y)))
xlim(0., fs/2)
ylim(-2*pi, 0.)

ylabel('Phase (rads)')
xlabel('Freq (Hz)')

gcf().set_size_inches(5., 4.)
savefig('transfer.png')
show()