# 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, 2016 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. Sep 29, 2016 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, 2016 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