# Problem with scipy.signal.butter for designing filter bank

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: 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 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

• Can you upsample your data so that its sampling frequency matches your filters'? – MBaz Sep 29 '16 at 21:47
• @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. – Dipole Sep 29 '16 at 21:50
• 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). – MBaz Sep 29 '16 at 22:12

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