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