Problem with scipy.signal.butter for designing filter bank
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