high_shelf2.py

import math
import numpy
import matplotlib.pyplot as plt

def IdealH(f):
    T = 1. / 441000.          # we use the tenfold sampling frequency
    OmegaU = 1. / 15e-6
    OmegaL = 15. / 50. * OmegaU
    # Calculate filter coefficients
    V0 = OmegaL / OmegaU
    H0 = V0 - 1.
    B = V0 * numpy.tan(OmegaU * T / 2.)
    A1 = (B - 1.) / (B + 1.)
    B0 = (1. + (1. - A1) * H0 / 2.)
    B1 = (A1 + (A1 - 1.) * H0 / 2.)
    # helper variables
    D = B1 / B0
    O = 2 * math.pi * T
    # Ideal transfer function
    return B0*numpy.sqrt((1 + 2*numpy.cos(f*O)*D + D*D)/(1 + 2*numpy.cos(f*O)*A1 + A1*A1))


def plot(a0, a1, a2, b1, b2):
    sampF = 44100.
    o = 2 * math.pi * (1. / sampF)
    f = numpy.arange(0, sampF / 2, sampF / 1000)
    #
    H = numpy.sqrt((a0*a0 + a1*a1 + a2*a2 + 2.*(a0*a1 + a1*a2)*numpy.cos(f*o) + 2.*(a0*a2)*numpy.cos(2.*f*o)) / (1. + b1*b1 + b2*b2 + 2.*(b1 + b1*b2)*numpy.cos(f*o) + 2.*b2*numpy.cos(2.*f*o)))
    IH = IdealH(f)
    #
    H = -20. * numpy.log10(H)
    IH = -20. * numpy.log10(IH)
    #
    plt.semilogx(f, H, 'r--')
    plt.semilogx(f, IH, 'g--')
    plt.show()


def foo(gain, cutoffFrequency, Q):
    samplingFrequency = 44100.
    #
    v = math.pow(10., abs(gain) / 20.0)
    k = math.tan(math.pi * cutoffFrequency / samplingFrequency)
    n = 0.0
    #
    if gain >= 0:
        n = 1. / (1. + k/Q + k * k)
        a0 = (v + math.sqrt(v)/Q * k + k * k) * n
        a1 = 2 * (k * k - v) * n
        a2 = (v - math.sqrt(v)/Q * k + k * k) * n
        b1 = 2 * (k * k - 1) * n
        b2 = (1 - k/Q + k * k) * n
    else:
        n = 1. / (v + math.sqrt(v)/Q * k + k * k)
        a0 = (1 + k/Q + k * k) * n
        a1 = 2 * (k * k - 1) * n
        a2 = (1 - k/Q + k * k) * n
        b1 = 2 * (k * k - v) * n
        b2 = (v - math.sqrt(v)/Q * k + k * k) * n
    #
    plot(a0, a1, a2, b1, b2)

cutoff = (10000. + 3100.)/2.
gain = -9.477
A = math.pow(10., gain / 2. / 20.)
slope = .5
Q = 1 / math.sqrt((A + 1. / A) * (1. / slope - 1.) + 2.)
foo(gain, cutoff, Q)