How to implement multi-band-pass filter with scipy.signal.butter

Question

Based on the band-pass filter here, I am trying to make a multi-band filter using the code bellow. However, the filtered signal is close to zero which affects the result when the spectrum is plotted. Should the coefficients of the filter of each band be normalized? Can you please someone suggest how I can fix the filter?

from scipy.signal import butter, sosfreqz, sosfilt
from scipy.signal import spectrogram
import matplotlib
import matplotlib.pyplot as plt
from scipy.fft import fft
import numpy as np


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


def multiband_filter(data, bands, fs, order=10):
    sos_list = []
    for lowcut, highcut in bands:
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        scalar = max(abs(fft(sos, 2000)))
        # sos = sos / scalar
        sos_list += [sos]

    # sos_list = [butter_bandpass(lowcut, highcut, fs, order=order) for lowcut, highcut in bands]

    # Combine filters into a single filter
    sos = np.vstack(sos_list)

    # Apply the multiband filter to the data
    y = sosfilt(sos, data)

    return y, sos_list


def get_toy_signal():
    t = np.arange(0, 0.3, 1 / fs)

    fq = [-np.inf] + [x / 12 for x in range(-9, 3, 1)]

    mel = [5, 3, 1, 3, 5, 5, 5, 0, 3, 3, 3, 0, 5, 8, 8, 0, 5, 3, 1, 3, 5, 5, 5, 5, 3, 3, 5, 3, 1]
    acc = [5, 0, 8, 0, 5, 0, 5, 5, 3, 0, 3, 3, 5, 0, 8, 8, 5, 0, 8, 0, 5, 5, 5, 0, 3, 3, 5, 0, 1]

    toy_signal = np.array([])

    for kj in range(len(mel)):
        note_signal = np.sum([np.sin(2 * np.pi * 440 * 2 ** ff * t)
                              for ff in [fq[acc[kj]] - 1, fq[acc[kj]], fq[mel[kj]] + 1]], axis=0)

        zeros = np.zeros(int(0.01 * fs))
        toy_signal = np.concatenate((toy_signal, note_signal, zeros))

    toy_signal += np.random.normal(0, 1, len(toy_signal))

    toy_signal = toy_signal / (np.max(np.abs(toy_signal)) + 0.1)
    t_toy_signal = np.arange(len(toy_signal)) / fs

    return t_toy_signal, toy_signal


if __name__ == "__main__":

    fontsize = 12
    # Sample rate and desired cut_off frequencies (in Hz).
    fs = 3000

    f1, f2 = 100, 200
    f3, f4 = 470, 750
    f5, f6 = 800, 850
    f7, f8 = 1000, 1000.1
    cut_off = [(f1, f2), (f3, f4), (f5, f6), (f7, f8)]
    # cut_off = [(f1, f2), (f3, f4)]
    # cut_off = [(f1, f2)]
    # cut_off = [f1]

    t_toy_signal, toy_signal = get_toy_signal()
    # toy_signal -= np.mean(toy_signal)
    # t_toy_signal = wiener(t_toy_signal)

    fig, ax = plt.subplots(6, 1, figsize=(8, 12))
    fig.tight_layout()

    ax[0].plot(t_toy_signal, toy_signal)
    ax[0].set_title('Original toy_signal', fontsize=fontsize)
    ax[0].set_xlabel('Time (s)', fontsize=fontsize)
    ax[0].set_ylabel('Magnitude', fontsize=fontsize)
    ax[0].set_xlim(left=0, right=max(t_toy_signal))

    sos_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in cut_off]

    # Combine filters into a single filter
    sos = np.vstack(sos_list)

    # w *= 0.5 * fs / np.pi  # Convert w to Hz.
    #####################################################################
    # First plot the desired ideal response as a green(ish) rectangle.
    #####################################################################

    # Plot the frequency response
    for i in range(len(cut_off)):
        w, h = sosfreqz(sos_list[i], worN=2000)
        ax[1].plot(0.5 * fs * w / np.pi, np.abs(h), label=f'Band {i + 1}: {cut_off[i]} Hz')

    ax[1].set_title('Multiband Filter Frequency Response')
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Gain')
    ax[1].legend()
    # ax[1].set_xlim(0, max(*cut_off) + 100)

    #####################################################################
    # Spectrogram of original signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal, fs,
                            nperseg=930, noverlap=0)
    ax[2].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx), vmax=np.max(Sxx)),
                     )

    ax[2].set_title('Spectrogram of original toy_signal', fontsize=fontsize)
    ax[2].set_xlabel('Time (s)', fontsize=fontsize)
    ax[2].set_ylabel('Frequency (Hz)', fontsize=fontsize)

    #####################################################################
    # Compute filtered signal
    #####################################################################

    # Apply the multiband filter to the data
    # toy_signal_filtered = sosfilt(sos, toy_signal)
    toy_signal_filtered = np.sum([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)

    #####################################################################
    # Spectrogram of filtered signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal_filtered, fs,
                            nperseg=930, noverlap=0)

    ax[3].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx),
                                                    vmax=np.max(Sxx))
                     )

    ax[3].set_title('Spectrogram of filtered toy_signal', fontsize=fontsize)
    ax[3].set_xlabel('Time (s)', fontsize=fontsize)
    ax[3].set_ylabel('Frequency (Hz)', fontsize=fontsize)

    ax[4].plot(t_toy_signal, toy_signal_filtered)
    ax[4].set_title('Filtered toy_signal', fontsize=fontsize)
    ax[4].set_xlim(left=0, right=max(t_toy_signal))
    ax[4].set_xlabel('Time (s)', fontsize=fontsize)
    ax[4].set_ylabel('Magnitude', fontsize=fontsize)

    N = 1512
    X = fft(toy_signal, n=N)
    Y = fft(toy_signal_filtered, n=N)

    # fig.set_size_inches((10, 4))
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), 'g-', label='FFT filtered signal')
    ax[5].set_xlim(xmax=fs / 2)
    ax[5].set_ylim(ymin=-20)
    ax[5].set_ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
    ax[5].set_xlabel("frequency (Hz)", fontsize=fontsize)
    ax[5].grid()
    ax[5].legend(loc='upper right')

    plt.tight_layout()
    plt.show()

    plt.figure()
    # fig.set_size_inches((10, 4))
    plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), 'g-', label='FFT filtered signal')
    plt.xlim(xmax=fs / 2)
    plt.ylim(ymin=-20)
    plt.ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
    plt.xlabel("frequency (Hz)", fontsize=fontsize)
    plt.grid()
    plt.legend(loc='upper right')
    plt.tight_layout()
    plt.show()

The following is after using @Warren Weckesser comment:

toy_signal_filtered = np.mean([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)

The following is after using @Warren Weckesser comment:

toy_signal_filtered = np.sum([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)

Here is an example where a narrow band is used:

Answer 1

Easier and recommended method is what Warren wrote in comments. Just calculate sum of separately band-pass filtered signals.

That being said, for someone who wants to create and apply single multi-band filter, he can try to achieve this by combining filters:

• lowpass (to cut everything above last pass-filter),
• highpass (to cut everything below first pass-filter),
• N-1 band-stop filters, where N being number of pass-bands (to cut parts in-between filters).

It may be difficult though to make it stable (be careful with filter orders) and harder to make it steep.

Found it interesting and tried myself:

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


def multi_band_filter(bands, subfilter_order=5):
    # high-pass filter
    nyq = 0.5 * fs
    normal_cutoff = bands[0][0] / nyq
    b, a = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False)
    all_b = [b]
    all_a = [a]

    # band-stop filters
    for idx in range(len(bands) - 1):
        normal_cutoff1 = bands[idx][1] / nyq
        normal_cutoff2 = bands[idx+1][0] / nyq
        b, a = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False)
        all_b.append(b)
        all_a.append(a)

    # low-pass filter
    normal_cutoff = bands[-1][1] / nyq
    b, a = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False)
    all_b.append(b)
    all_a.append(a)

    # combine filters:
    combined_a = all_a[0]
    for a in all_a[1:]:
        combined_a = np.convolve(a, combined_a)
    combined_b = all_b[0]
    for b in all_b[1:]:
        combined_b = np.convolve(b, combined_b)

    return combined_b, combined_a


bands = [[400, 700], [1000, 1500]]

fs = 8000
time = np.arange(0, 1 - 0.5/fs, 1/fs)
signal_to_filter = np.sum([np.sin(2 * np.pi * (freq + 0.01 * np.random.random()) * time + np.pi*np.random.random()) for freq in range(10, 3800)], axis=0)

b, a = multi_band_filter(bands)
filtered_signal = lfilter(b, a, signal_to_filter)
original_spectrum = fft(signal_to_filter)
filtered_signal_spectrum = fft(filtered_signal)

plt.figure(figsize=(16, 10))
plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
plt.xlim([0, 4000])
plt.show()

SOS version

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


def multi_band_filter(bands, subfilter_order=5):
    # high-pass filter
    nyq = 0.5 * fs
    normal_cutoff = bands[0][0] / nyq
    sos = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False, output='sos')
    all_sos = [sos]

    # band-stop filters
    for idx in range(len(bands) - 1):
        normal_cutoff1 = bands[idx][1] / nyq
        normal_cutoff2 = bands[idx+1][0] / nyq
        sos = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False, output='sos')
        all_sos.append(sos)

    # low-pass filter
    normal_cutoff = bands[-1][1] / nyq
    sos = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False, output='sos')
    all_sos.append(sos)

    # combine filters:
    combined_sos = np.vstack(all_sos)
    return combined_sos


bands = [[400, 700], [1000, 1500]]

fs = 8000
time = np.arange(0, 1 - 0.5/fs, 1/fs)
signal_to_filter = np.sum([np.sin(2 * np.pi * (freq + 0.01 * np.random.random()) * time + np.pi*np.random.random()) for freq in range(10, 3800)], axis=0)

sos = multi_band_filter(bands)
filtered_signal = sosfilt(sos, signal_to_filter)
original_spectrum = fft(signal_to_filter)
filtered_signal_spectrum = fft(filtered_signal)

plt.figure(figsize=(16, 10))
plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
plt.xlim([0, 4000])
plt.show()

w, h = freqz(b, a)
freq_domain = np.linspace(0, fs/2, len(w))
plt.figure(figsize=(16, 10))
plt.plot(freq_domain, 20 * np.log10(abs(h)), 'b')
plt.show()

As you can see, the slope of the filter is not very steep.

Answer 2

Your Filtered Signals get "multiplicated" by each other because you calculate them in Series. In my code below i changed your code so it calculates your filtered_signals one by one and then plots them one by one on your graphs. I hope thats what you tried to achieve:

from scipy.signal import butter, sosfreqz, sosfilt
from scipy.signal import spectrogram
import matplotlib
import matplotlib.pyplot as plt
from scipy.fft import fft
import numpy as np
import scipy


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    print(high)
    sos = butter(order, [low, high], analog=False, btype='band', output='sos')
    return sos


def multiband_filter(data, bands, fs, order=10):
    sos_list = []
    for lowcut, highcut in bands:
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        scalar = max(abs(fft(sos, 2000)))
        sos = sos / scalar
        sos_list += [sos]

    # sos_list = [butter_bandpass(lowcut, highcut, fs, order=order) for lowcut, highcut in bands]

    # Combine filters into a single filter
    sos = np.vstack(sos_list)

    # Apply the multiband filter to the data
    y = sosfilt(sos, data)

    return y, sos_list


def get_toy_signal():
    t = np.arange(0, 0.3, 1 / fs)

    fq = [-np.inf] + [x / 12 for x in range(-9, 3, 1)]

    mel = [5, 3, 1, 3, 5, 5, 5, 0, 3, 3, 3, 0, 5, 8, 8, 0, 5, 3, 1, 3, 5, 5, 5, 5, 3, 3, 5, 3, 1]
    acc = [5, 0, 8, 0, 5, 0, 5, 5, 3, 0, 3, 3, 5, 0, 8, 8, 5, 0, 8, 0, 5, 5, 5, 0, 3, 3, 5, 0, 1]

    toy_signal = np.array([])

    for kj in range(len(mel)):
        note_signal = np.sum([np.sin(2 * np.pi * 440 * 2 ** ff * t)
                              for ff in [fq[acc[kj]] - 1, fq[acc[kj]], fq[mel[kj]] + 1]], axis=0)

        zeros = np.zeros(int(0.01 * fs))
        toy_signal = np.concatenate((toy_signal, note_signal, zeros))

    toy_signal += np.random.normal(0, 1, len(toy_signal))

    toy_signal = toy_signal / (np.max(np.abs(toy_signal)) + 0.1)
    t_toy_signal = np.arange(len(toy_signal)) / fs

    return t_toy_signal, toy_signal


if __name__ == "__main__":

    fontsize = 12
    # Sample rate and desired cut_off frequencies (in Hz).
    fs = 3000

    f1, f2 = 100, 200
    f3, f4 = 470, 750
    f5, f6 = 800, 850
    cut_off = [(f1, f2), (f3, f4), (f5, f6)]
    #cut_off = [(f1, f2), (f3, f4)]
    #cut_off = [(f1, f2)]
    # cut_off = [f1]

    t_toy_signal, toy_signal = get_toy_signal()
    # toy_signal -= np.mean(toy_signal)
    # t_toy_signal = wiener(t_toy_signal)

    fig, ax = plt.subplots(6, 1, figsize=(8, 12))
    fig.tight_layout()

    ax[0].plot(t_toy_signal, toy_signal)
    ax[0].set_title('Original toy_signal', fontsize=fontsize)
    ax[0].set_xlabel('Time (s)', fontsize=fontsize)
    ax[0].set_ylabel('Magnitude', fontsize=fontsize)
    ax[0].set_xlim(left=0, right=max(t_toy_signal))
    sos_lists=[]
    for cut in cut_off:
        print(cut)
        sos_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in [cut]]
        sos_lists.append(np.vstack(sos_list))
    # Combine filters into a single filter
    sos = np.vstack(sos_list)



    # w *= 0.5 * fs / np.pi  # Convert w to Hz.
    #####################################################################
    # First plot the desired ideal response as a green(ish) rectangle.
    #####################################################################

    # Plot the frequency response
    for i in range(len(cut_off)):
        w, h = sosfreqz(sos_lists[i], worN=2000)
        ax[1].plot(0.5 * fs * w / np.pi, np.abs(h), label=f'Band {i + 1}: {cut_off[i]} Hz')

    ax[1].set_title('Multiband Filter Frequency Response')
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Gain')
    ax[1].legend()
    # ax[1].set_xlim(0, max(*cut_off) + 100)

    #####################################################################
    # Spectrogram of original signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal, fs,
                            nperseg=930, noverlap=0)
    ax[2].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx), vmax=np.max(Sxx)),
                     )

    ax[2].set_title('Spectrogram of original toy_signal', fontsize=fontsize)
    ax[2].set_xlabel('Time (s)', fontsize=fontsize)
    ax[2].set_ylabel('Frequency (Hz)', fontsize=fontsize)

    #####################################################################
    # Compute filtered signal
    #####################################################################

    # Apply the multiband filter to the data
    toy_signal_filtereds=[]
    for sos in sos_lists:
        toy_signal_filtered = sosfilt(sos, toy_signal, )
        toy_signal_filtereds.append([toy_signal_filtered])


    #####################################################################
    # Spectrogram of filtered signal
    #####################################################################

    f, t, Sxx = spectrogram(toy_signal_filtered, fs,
                            nperseg=930, noverlap=0)

    ax[3].pcolormesh(t, f, np.abs(Sxx),
                     norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx),
                                                    vmax=np.max(Sxx))
                     )

    ax[3].set_title('Spectrogram of filtered toy_signal', fontsize=fontsize)
    ax[3].set_xlabel('Time (s)', fontsize=fontsize)
    ax[3].set_ylabel('Frequency (Hz)', fontsize=fontsize)
    for toy_signal_filtered in toy_signal_filtereds:
        ax[4].plot(t_toy_signal, toy_signal_filtered[0])
    ax[4].set_title('Filtered toy_signal', fontsize=fontsize)
    ax[4].set_xlim(left=0, right=max(t_toy_signal))
    ax[4].set_xlabel('Time (s)', fontsize=fontsize)
    ax[4].set_ylabel('Magnitude', fontsize=fontsize)

    N = 1512
    X = fft(toy_signal, n=N)
    Ys=[]
    for toy_signal_filtered in toy_signal_filtereds:
        Y = fft(toy_signal_filtered[0], n=N)
        Ys.append(Y)

    # fig.set_size_inches((10, 4))
    ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
    for bigY in Ys:
        ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(bigY)), 'g-', label='FFT filtered signal')
    ax[5].set_xlim(xmax=fs / 2)
    ax[5].set_ylim(ymin=-20)
    ax[5].set_ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
    ax[5].set_xlabel("frequency (Hz)", fontsize=fontsize)
    ax[5].grid()
    ax[5].legend(loc=0)

    plt.tight_layout()
    plt.show()

For the Frequency response in your case you should probably dont use a logscaled plot. ax[1].semilogy however makes a plot with log scaling on the y-axis.