Bandpass filters with python for low frequencies

I have a signal, sampled at 1 kHz, that I would like to analyze in third octave bands. For that purpose I defined Butterworth filters with proportional bandwidth. This works well for higher frequencies, but produces very weird results for low frequencies. What am I missing? Plotted are the frequency responses of the Butterworth filters. I find similar behaviour in Hanning filters, and think that I am missing something fundamental.

Code used to define filters:

import scipy
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

def butter_bandpass(lowcut, highcut, fs, order=5, label=None):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = signal.butter(order, [low, high], btype='band')
    w, h = signal.freqz(b,a,worN=2000)
    plt.plot((fs * 0.5 / np.pi) * w, abs(h), label=label)
    return b, a

center_freqs = np.array([0.5* 2 ** (n * 1 / 3.) for n in range(0, 29)])

center_freqs.sort()
lower_freqs = 2. ** (-1 / 6.) * center_freqs
for lf in lower_freqs:
    butter_bandpass(lf, lf*2**(1/3.), fs=1000, order=5)

plt.show()

Answer

I would bet this is just numerical error in the transfer function. Try using butter_sos = butter(..., output='sos') instead of ba format, and sosfreqz(butter_sos, ...) instead of freqz. Does that fix it?

https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.sosfreqz.html

Edit: I just tried it and it does. :D

import scipy
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
def butter_bandpass(lowcut, highcut, fs, order=5, label=None):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = signal.butter(order, [low, high], btype='band', output='sos')
    w, h = signal.sosfreqz(sos,worN=20000)

    plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label=label)
    return sos

center_freqs = np.array([0.5* 2 ** (n * 1 / 3.) for n in range(0, 29)])
center_freqs.sort()
lower_freqs = 2. ** (-1 / 6.) * center_freqs
for lf in lower_freqs:
    butter_bandpass(lf, lf*2**(1/3.), fs=1000, order=5)

plt.show()