2
$\begingroup$

I need to design a filter bank of band pass filters (code for implementation is shown below). There is a problem in that if my center frequencies fc are much smaller than the sampling frequency fs then you get the behavior as shown in the first graph:

enter image description here

A temporary solution is to change the sampling frequency to fs = fh*4 where 'fh' is the hi cutoff frequency of each respective band. This gives

enter image description here

Which is what I want, however because the sampling frequency is no longer the true one (relative to the time series I then want to filter) I cant use it for filtering purposes using lfilter(b,a,x) where x is my time series because all of the bands will be identical since they have the same non-dimensional frequency! My question is what can I do to make this work as expected for a true fs? Presumably there is an issue with finite precision?

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter
from scipy.signal import freqz



def butter_bandpass(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def fbank(fs):
    """
    Calculate third octave filter bank.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    # ISO standard center frequencies (bins)
    F0 = [50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500,
          630, 800, 1000, 1250, 1600, 2000, 2500, 3150, 4000,
          5000, 6300, 8000, 10000, 12500, 16000, 20000]


    plt.figure(figsize=(10,6))
    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        fs = fh*4 # Temporary solution
        fs = 5e4
        b, a = butter_bandpass(fl,fh,fs)
        w, h = freqz(b, a, worN=2000)
        plt.loglog((fs * 0.5 / np.pi) * w, abs(h), label="fc = %d" % fc)

    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Gain")
    plt.legend(loc='best')
    plt.savefig("test_filter.png", bbox_inches='tight')
    plt.show()
    return 0
$\endgroup$
  • $\begingroup$ Can you upsample your data so that its sampling frequency matches your filters'? $\endgroup$ – MBaz Sep 29 '16 at 21:47
  • $\begingroup$ @MBaz It would be a bit messy - I would need to subsample my data differently for each frequency band and im also not sure about how that would affect the result in terms of overall power. $\endgroup$ – Dipole Sep 29 '16 at 21:50
  • $\begingroup$ Yeah, designing narrow bandpass filters is hard in general. Increasing the filter order helps, but only up to a point. One approach could be to design a wider filter and apply it multiple times. Also, you could try centering the filter at a more convenient center frequency, and then upconvert/downconvert your data to match the filter. (But I'm no expert in filter design -- hopefully someone else will chime in). $\endgroup$ – MBaz Sep 29 '16 at 22:12
1
$\begingroup$

Although Im not sure how I can display the frequency response using this method, I have found that it is workable to use second-order section method for filtering. So say we have a signal x with sampling frequency fs_real that needs to filtered using a constant percentage bandwidth (in this case third octaves) then the corresponding power spectrum will be given by the get_cpb function:

def butter_bandpass_sos(lowcut, highcut, fs, order=9):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = butter(order, [low, high], btype='band', output='sos')
    return sos

def get_cpb(x,fs_real):
    """
    Calculate third octave spectrum from pressure time series.
    An octave is the interval between two frequencies having
    a ratio of 2:1. Hence for 1/3 octave
    fh = fc*(sqrt(2)^(1/3))
    fl = fc/(sqrt(2)^(1/3))
    """

    cpb = []

    for fc in F0:
        fh = fc*(2**(1/6.0))
        fl = fc/(2**(1/6.0))
        sos = butter_bandpass_sos(fl,fh,fs_real)
        y = sosfilt(sos,x)
        cpb.append( 20*np.log10(np.std(y)/2e-5) )


    return F0, cpb
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.