フィルタ

##### FILTER ##################################################################

def biquad_lpf(fc, samprate):
    """
    Biquad low-pass filter
    Usage:
       samprate = 44100
       b, a = biquad_lpf(1000, samprate)
       w, h = sig.freqz(b, a)
       plt.plot(w/np.pi*(samprate/2.0), amplitude_to_db(np.absolute(h)))
       plt.xlabel("Frequency (Hz)")
       plt.ylabel("Amplitude (dB)")
       plt.grid()
       plt.show()
    """
    K = np.tan(np.pi * fc / samprate)
    b0 =   K**2 / (1 + np.sqrt(2)*K + K**2)
    b1 = 2*K**2 / (1 + np.sqrt(2)*K + K**2)
    b2 =   K**2 / (1 + np.sqrt(2)*K + K**2)
    a0 = 1
    a1 = 2*(K**2-1) / (1 + np.sqrt(2)*K + K**2)
    a2 = (1 - np.sqrt(2)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
    b = [b0, b1, b2]
    a = [a0, a1, a2]
    return(b, a)

def biquad_hpf(fc, samprate):
    K = np.tan(np.pi * fc / samprate)
    b0 =  1 / (1 + np.sqrt(2)*K + K**2)
    b1 = -2 / (1 + np.sqrt(2)*K + K**2)
    b2 =  1 / (1 + np.sqrt(2)*K + K**2)
    a0 = 1
    a1 = 2*(K**2-1) / (1 + np.sqrt(2)*K + K**2)
    a2 = (1 - np.sqrt(2)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
    b = [b0, b1, b2]
    a = [a0, a1, a2]
    return(b, a)

def biquad_lowshelf(fc, G, samprate):
    K = np.tan(np.pi * fc / samprate)
    if (G >= 0.0):
        V0 = 10 ** (G/20)
        b0 = (1 + np.sqrt(2*V0)*K + V0*K**2) / (1 + np.sqrt(2)*K + K**2)
        b1 =                    2*(V0*K**2 - 1) / (1 + np.sqrt(2)*K + K**2)
        b2 = (1 - np.sqrt(2*V0)*K + V0*K**2) / (1 + np.sqrt(2)*K + K**2)
        a0 = 1
        a1 =                       2*(K**2 - 1) / (1 + np.sqrt(2)*K + K**2)
        a2 =       (1 - np.sqrt(2)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
    else:
        V0 = 10 ** (-G/20)
        b0 =       (1 + np.sqrt(2)*K + K**2) / (1 + np.sqrt(2*V0)*K + V0*K**2)
        b1 =                       2*(K**2 - 1) / (1 + np.sqrt(2*V0)*K + V0*K**2)
        b2 =       (1 - np.sqrt(2)*K + K**2) / (1 + np.sqrt(2*V0)*K + V0*K**2)
        a0 = 1
        a1 =                    2*(V0*K**2 - 1) / (1 + np.sqrt(2*V0)*K + V0*K**2)
        a2 = (1 - np.sqrt(2*V0)*K + V0*K**2) / (1 + np.sqrt(2*V0)*K + V0*K**2)
    b = [b0, b1, b2]
    a = [a0, a1, a2]
    return(b, a)

def biquad_highshelf(fc, G, samprate):
    K = np.tan(np.pi * fc / samprate)
    if (G >= 0.0):
        V0 = 10 ** (G/20)
        b0 =   (V0 + np.sqrt(2*V0)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
        b1 =                      2*(K**2 - V0) / (1 + np.sqrt(2)*K + K**2)
        b2 =   (V0 - np.sqrt(2*V0)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
        a0 = 1
        a1 =                       2*(K**2 - 1) / (1 + np.sqrt(2)*K + K**2)
        a2 =       (1 - np.sqrt(2)*K + K**2) / (1 + np.sqrt(2)*K + K**2)
    else:
        V0 = 10 ** (-G/20)
        b0 =       (1 + np.sqrt(2)*K + K**2) / (V0 + np.sqrt(2*V0)*K + K**2)
        b1 =                       2*(K**2 - 1) / (V0 + np.sqrt(2*V0)*K + K**2)
        b2 =       (1 - np.sqrt(2)*K + K**2) / (V0 + np.sqrt(2*V0)*K + K**2)
        a0 = 1
        a1 =              2*(K**2/V0 - 1) / (1 + np.sqrt(2/V0)*K + K**2/V0)
        a2 = (1 - np.sqrt(2/V0)*K + K**2/V0) / (1 + np.sqrt(2/V0)*K + K**2/V0)
    b = [b0, b1, b2]
    a = [a0, a1, a2]
    return(b, a)

def biquad_peak(fc, G, Q, samprate):
    K = np.tan(np.pi * fc / samprate)
    if (G >= 0.0):
        V0 = 10 ** (G/20)
        b0 = (1 + V0/Q*K + K**2) / (1 + 1/Q*K + K**2)
        b1 =        2*(K**2 - 1) / (1 + 1/Q*K + K**2)
        b2 = (1 - V0/Q*K + K**2) / (1 + 1/Q*K + K**2)
        a0 = 1
        a1 =        2*(K**2 - 1) / (1 + 1/Q*K + K**2)
        a2 =  (1 - 1/Q*K + K**2) / (1 + 1/Q*K + K**2)
    else:
        V0 = 10 ** (-G/20)
        b0 =  (1 + 1/Q*K + K**2) / (1 + V0/Q*K + K**2)
        b1 =        2*(K**2 - 1) / (1 + V0/Q*K + K**2)
        b2 =  (1 - 1/Q*K + K**2) / (1 + V0/Q*K + K**2)
        a0 = 1
        a1 =        2*(K**2 - 1) / (1 + V0/Q*K + K**2)
        a2 = (1 - V0/Q*K + K**2) / (1 + V0/Q*K + K**2)
    b = [b0, b1, b2]
    a = [a0, a1, a2]
    return(b, a)

def filterbank(x, samprate, oct=1):
    ”'
    FILTERBANK   octave filterbank using raised cosine frequency window
        out, cfreqs = filterbank(x, samprate) create filtered signals `out' from x.
        Columns of `out' is a filtered signal with corresponding center
        frequency in `cfreqs'.  By default, 1/1-octave band filters are used.
        Filtered signals at the lowest frequency is processed with low-pass,
        the highest frequency is processed with high-pass, and all frequencies
        in between are processed with band-pass filters.  This yields to
        (almost) perfect reconstruction of the original signal, such that:
            out, cfreqs = filterbank(x, samprate)
            z = sum(out, 1)
        creates z almost the same as x.

        out, cfreqs = filterbank(x, samprate, oct) creates and filters using
        1/oct-octave band filters.
    ”'
    # prepare center frequencies
    cfreqs = 1000 * 2**(np.arange(-100, 100)/oct)
    cfreqs = cfreqs[np.logical_and(20<=cfreqs, cfreqs<=samprate/2)]
    # create octave-band bandpass filter
    if len(x) % 2 == 1:
        x = np.concatenate((x, [0]))
    out = np.zeros((len(x), len(cfreqs)))
    xx = scipy.fft(x)
    for n in range(len(cfreqs)):
        fc = cfreqs[n]
        fl = fc * 2**(-1/oct)
        fh = fc * 2**(+1/oct)
        W = np.zeros(len(xx)//2)
        f = np.linspace(0, samprate/2, len(W))
        ff = np.logical_and(fl<=f, f<=fh)
        cc = (np.cos(np.log2(f[ff]/fc)*math.pi*oct) + 1) / 2
        W[ff] = cc
        if n==0:
            W[f < fc] = 1
        elif n==len(cfreqs)-1:
            W[fc < f] = 1
        W = np.concatenate((W, np.flipud(W)))
        W[0] = 0
        W[len(W)//2] = 0
        # apply filter
        y = np.real(scipy.ifft(xx * W))
        out[:,n] = y
    return(out, cfreqs)