如何用scipy.signal.butter实现多带通滤波器

9jyewag0  于 2023-11-19  发布在  其他
关注(0)|答案(2)|浏览(144)

基于带通滤波器here,我尝试使用以下代码制作一个多波段滤波器。然而,滤波后的信号接近零,这影响了频谱绘制的结果。每个波段的滤波器系数是否应该归一化?请有人建议我如何修复滤波器?

  1. from scipy.signal import butter, sosfreqz, sosfilt
  2. from scipy.signal import spectrogram
  3. import matplotlib
  4. import matplotlib.pyplot as plt
  5. from scipy.fft import fft
  6. import numpy as np
  7. def butter_bandpass(lowcut, highcut, fs, order=5):
  8. nyq = 0.5 * fs
  9. low = lowcut / nyq
  10. high = highcut / nyq
  11. sos = butter(order, [low, high], analog=False, btype='band', output='sos')
  12. return sos
  13. def multiband_filter(data, bands, fs, order=10):
  14. sos_list = []
  15. for lowcut, highcut in bands:
  16. sos = butter_bandpass(lowcut, highcut, fs, order=order)
  17. scalar = max(abs(fft(sos, 2000)))
  18. # sos = sos / scalar
  19. sos_list += [sos]
  20. # sos_list = [butter_bandpass(lowcut, highcut, fs, order=order) for lowcut, highcut in bands]
  21. # Combine filters into a single filter
  22. sos = np.vstack(sos_list)
  23. # Apply the multiband filter to the data
  24. y = sosfilt(sos, data)
  25. return y, sos_list
  26. def get_toy_signal():
  27. t = np.arange(0, 0.3, 1 / fs)
  28. fq = [-np.inf] + [x / 12 for x in range(-9, 3, 1)]
  29. 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]
  30. 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]
  31. toy_signal = np.array([])
  32. for kj in range(len(mel)):
  33. note_signal = np.sum([np.sin(2 * np.pi * 440 * 2 ** ff * t)
  34. for ff in [fq[acc[kj]] - 1, fq[acc[kj]], fq[mel[kj]] + 1]], axis=0)
  35. zeros = np.zeros(int(0.01 * fs))
  36. toy_signal = np.concatenate((toy_signal, note_signal, zeros))
  37. toy_signal += np.random.normal(0, 1, len(toy_signal))
  38. toy_signal = toy_signal / (np.max(np.abs(toy_signal)) + 0.1)
  39. t_toy_signal = np.arange(len(toy_signal)) / fs
  40. return t_toy_signal, toy_signal
  41. if __name__ == "__main__":
  42. fontsize = 12
  43. # Sample rate and desired cut_off frequencies (in Hz).
  44. fs = 3000
  45. f1, f2 = 100, 200
  46. f3, f4 = 470, 750
  47. f5, f6 = 800, 850
  48. f7, f8 = 1000, 1000.1
  49. cut_off = [(f1, f2), (f3, f4), (f5, f6), (f7, f8)]
  50. # cut_off = [(f1, f2), (f3, f4)]
  51. # cut_off = [(f1, f2)]
  52. # cut_off = [f1]
  53. t_toy_signal, toy_signal = get_toy_signal()
  54. # toy_signal -= np.mean(toy_signal)
  55. # t_toy_signal = wiener(t_toy_signal)
  56. fig, ax = plt.subplots(6, 1, figsize=(8, 12))
  57. fig.tight_layout()
  58. ax[0].plot(t_toy_signal, toy_signal)
  59. ax[0].set_title('Original toy_signal', fontsize=fontsize)
  60. ax[0].set_xlabel('Time (s)', fontsize=fontsize)
  61. ax[0].set_ylabel('Magnitude', fontsize=fontsize)
  62. ax[0].set_xlim(left=0, right=max(t_toy_signal))
  63. sos_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in cut_off]
  64. # Combine filters into a single filter
  65. sos = np.vstack(sos_list)
  66. # w *= 0.5 * fs / np.pi # Convert w to Hz.
  67. #####################################################################
  68. # First plot the desired ideal response as a green(ish) rectangle.
  69. #####################################################################
  70. # Plot the frequency response
  71. for i in range(len(cut_off)):
  72. w, h = sosfreqz(sos_list[i], worN=2000)
  73. ax[1].plot(0.5 * fs * w / np.pi, np.abs(h), label=f'Band {i + 1}: {cut_off[i]} Hz')
  74. ax[1].set_title('Multiband Filter Frequency Response')
  75. ax[1].set_xlabel('Frequency [Hz]')
  76. ax[1].set_ylabel('Gain')
  77. ax[1].legend()
  78. # ax[1].set_xlim(0, max(*cut_off) + 100)
  79. #####################################################################
  80. # Spectrogram of original signal
  81. #####################################################################
  82. f, t, Sxx = spectrogram(toy_signal, fs,
  83. nperseg=930, noverlap=0)
  84. ax[2].pcolormesh(t, f, np.abs(Sxx),
  85. norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx), vmax=np.max(Sxx)),
  86. )
  87. ax[2].set_title('Spectrogram of original toy_signal', fontsize=fontsize)
  88. ax[2].set_xlabel('Time (s)', fontsize=fontsize)
  89. ax[2].set_ylabel('Frequency (Hz)', fontsize=fontsize)
  90. #####################################################################
  91. # Compute filtered signal
  92. #####################################################################
  93. # Apply the multiband filter to the data
  94. # toy_signal_filtered = sosfilt(sos, toy_signal)
  95. toy_signal_filtered = np.sum([sosfilt(sos, toy_signal) for sos in sos_list], axis=0)
  96. #####################################################################
  97. # Spectrogram of filtered signal
  98. #####################################################################
  99. f, t, Sxx = spectrogram(toy_signal_filtered, fs,
  100. nperseg=930, noverlap=0)
  101. ax[3].pcolormesh(t, f, np.abs(Sxx),
  102. norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx),
  103. vmax=np.max(Sxx))
  104. )
  105. ax[3].set_title('Spectrogram of filtered toy_signal', fontsize=fontsize)
  106. ax[3].set_xlabel('Time (s)', fontsize=fontsize)
  107. ax[3].set_ylabel('Frequency (Hz)', fontsize=fontsize)
  108. ax[4].plot(t_toy_signal, toy_signal_filtered)
  109. ax[4].set_title('Filtered toy_signal', fontsize=fontsize)
  110. ax[4].set_xlim(left=0, right=max(t_toy_signal))
  111. ax[4].set_xlabel('Time (s)', fontsize=fontsize)
  112. ax[4].set_ylabel('Magnitude', fontsize=fontsize)
  113. N = 1512
  114. X = fft(toy_signal, n=N)
  115. Y = fft(toy_signal_filtered, n=N)
  116. # fig.set_size_inches((10, 4))
  117. ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
  118. ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), 'g-', label='FFT filtered signal')
  119. ax[5].set_xlim(xmax=fs / 2)
  120. ax[5].set_ylim(ymin=-20)
  121. ax[5].set_ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
  122. ax[5].set_xlabel("frequency (Hz)", fontsize=fontsize)
  123. ax[5].grid()
  124. ax[5].legend(loc='upper right')
  125. plt.tight_layout()
  126. plt.show()
  127. plt.figure()
  128. # fig.set_size_inches((10, 4))
  129. plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
  130. plt.plot(np.arange(N) / N * fs, 20 * np.log10(abs(Y)), 'g-', label='FFT filtered signal')
  131. plt.xlim(xmax=fs / 2)
  132. plt.ylim(ymin=-20)
  133. plt.ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
  134. plt.xlabel("frequency (Hz)", fontsize=fontsize)
  135. plt.grid()
  136. plt.legend(loc='upper right')
  137. plt.tight_layout()
  138. plt.show()

字符串


的数据
以下是使用@Warren Weckesser后的注解:

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


的数据
以下是使用@Warren Weckesser后的评论:

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



下面是一个使用窄带的示例:


9ceoxa92

9ceoxa921#

Warren在评论中写的是一种更好的推荐方法,只需计算单独带通滤波信号的总和。
话虽如此,对于想要创建和应用 * 单 * 多波段滤波器的人来说,他可以尝试通过combining filters来实现这一点:

  • 低通(以削减最后一个通滤波器以上的所有内容),
  • 高通(削减一切低于第一通滤波器),
  • N-1个带阻滤波器,其中N是通带的数量(以切割滤波器之间的部分)。

但要使其稳定可能很困难(小心滤波器阶数),使其陡峭也更难。
我觉得这很有趣,我自己也试过:

  1. from scipy.signal import butter, lfilter
  2. import matplotlib.pyplot as plt
  3. from scipy.fft import fft
  4. import numpy as np
  5. def multi_band_filter(bands, subfilter_order=5):
  6. # high-pass filter
  7. nyq = 0.5 * fs
  8. normal_cutoff = bands[0][0] / nyq
  9. b, a = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False)
  10. all_b = [b]
  11. all_a = [a]
  12. # band-stop filters
  13. for idx in range(len(bands) - 1):
  14. normal_cutoff1 = bands[idx][1] / nyq
  15. normal_cutoff2 = bands[idx+1][0] / nyq
  16. b, a = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False)
  17. all_b.append(b)
  18. all_a.append(a)
  19. # low-pass filter
  20. normal_cutoff = bands[-1][1] / nyq
  21. b, a = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False)
  22. all_b.append(b)
  23. all_a.append(a)
  24. # combine filters:
  25. combined_a = all_a[0]
  26. for a in all_a[1:]:
  27. combined_a = np.convolve(a, combined_a)
  28. combined_b = all_b[0]
  29. for b in all_b[1:]:
  30. combined_b = np.convolve(b, combined_b)
  31. return combined_b, combined_a
  32. bands = [[400, 700], [1000, 1500]]
  33. fs = 8000
  34. time = np.arange(0, 1 - 0.5/fs, 1/fs)
  35. 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)
  36. b, a = multi_band_filter(bands)
  37. filtered_signal = lfilter(b, a, signal_to_filter)
  38. original_spectrum = fft(signal_to_filter)
  39. filtered_signal_spectrum = fft(filtered_signal)
  40. plt.figure(figsize=(16, 10))
  41. plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
  42. plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
  43. plt.xlim([0, 4000])
  44. plt.show()

字符串

SOS版本

  1. from scipy.signal import butter, sosfilt, freqz
  2. import matplotlib.pyplot as plt
  3. from scipy.fft import fft
  4. import numpy as np
  5. def multi_band_filter(bands, subfilter_order=5):
  6. # high-pass filter
  7. nyq = 0.5 * fs
  8. normal_cutoff = bands[0][0] / nyq
  9. sos = butter(subfilter_order, normal_cutoff, btype='highpass', analog=False, output='sos')
  10. all_sos = [sos]
  11. # band-stop filters
  12. for idx in range(len(bands) - 1):
  13. normal_cutoff1 = bands[idx][1] / nyq
  14. normal_cutoff2 = bands[idx+1][0] / nyq
  15. sos = butter(subfilter_order, [normal_cutoff1, normal_cutoff2], btype='bandstop', analog=False, output='sos')
  16. all_sos.append(sos)
  17. # low-pass filter
  18. normal_cutoff = bands[-1][1] / nyq
  19. sos = butter(subfilter_order, normal_cutoff, btype='lowpass', analog=False, output='sos')
  20. all_sos.append(sos)
  21. # combine filters:
  22. combined_sos = np.vstack(all_sos)
  23. return combined_sos
  24. bands = [[400, 700], [1000, 1500]]
  25. fs = 8000
  26. time = np.arange(0, 1 - 0.5/fs, 1/fs)
  27. 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)
  28. sos = multi_band_filter(bands)
  29. filtered_signal = sosfilt(sos, signal_to_filter)
  30. original_spectrum = fft(signal_to_filter)
  31. filtered_signal_spectrum = fft(filtered_signal)
  32. plt.figure(figsize=(16, 10))
  33. plt.plot(np.linspace(0, fs, len(original_spectrum)), np.abs(original_spectrum), color='b')
  34. plt.plot(np.linspace(0, fs, len(filtered_signal_spectrum)), np.abs(filtered_signal_spectrum), color='orange')
  35. plt.xlim([0, 4000])
  36. plt.show()
  37. w, h = freqz(b, a)
  38. freq_domain = np.linspace(0, fs/2, len(w))
  39. plt.figure(figsize=(16, 10))
  40. plt.plot(freq_domain, 20 * np.log10(abs(h)), 'b')
  41. plt.show()


100d 1xx 1c 1d 1x的字符串
正如你所看到的,滤波器的斜率不是很陡。

展开查看全部
swvgeqrz

swvgeqrz2#

你的滤波信号得到“乘”由对方,因为你计算他们在系列.在我下面的代码中,我改变了你的代码,所以它计算你的filtered_signals一个接一个,然后把它们一个接一个地绘制在你的图.我希望这就是你试图实现:


的数据

  1. from scipy.signal import butter, sosfreqz, sosfilt
  2. from scipy.signal import spectrogram
  3. import matplotlib
  4. import matplotlib.pyplot as plt
  5. from scipy.fft import fft
  6. import numpy as np
  7. import scipy
  8. def butter_bandpass(lowcut, highcut, fs, order=5):
  9. nyq = 0.5 * fs
  10. low = lowcut / nyq
  11. high = highcut / nyq
  12. print(high)
  13. sos = butter(order, [low, high], analog=False, btype='band', output='sos')
  14. return sos
  15. def multiband_filter(data, bands, fs, order=10):
  16. sos_list = []
  17. for lowcut, highcut in bands:
  18. sos = butter_bandpass(lowcut, highcut, fs, order=order)
  19. scalar = max(abs(fft(sos, 2000)))
  20. sos = sos / scalar
  21. sos_list += [sos]
  22. # sos_list = [butter_bandpass(lowcut, highcut, fs, order=order) for lowcut, highcut in bands]
  23. # Combine filters into a single filter
  24. sos = np.vstack(sos_list)
  25. # Apply the multiband filter to the data
  26. y = sosfilt(sos, data)
  27. return y, sos_list
  28. def get_toy_signal():
  29. t = np.arange(0, 0.3, 1 / fs)
  30. fq = [-np.inf] + [x / 12 for x in range(-9, 3, 1)]
  31. 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]
  32. 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]
  33. toy_signal = np.array([])
  34. for kj in range(len(mel)):
  35. note_signal = np.sum([np.sin(2 * np.pi * 440 * 2 ** ff * t)
  36. for ff in [fq[acc[kj]] - 1, fq[acc[kj]], fq[mel[kj]] + 1]], axis=0)
  37. zeros = np.zeros(int(0.01 * fs))
  38. toy_signal = np.concatenate((toy_signal, note_signal, zeros))
  39. toy_signal += np.random.normal(0, 1, len(toy_signal))
  40. toy_signal = toy_signal / (np.max(np.abs(toy_signal)) + 0.1)
  41. t_toy_signal = np.arange(len(toy_signal)) / fs
  42. return t_toy_signal, toy_signal
  43. if __name__ == "__main__":
  44. fontsize = 12
  45. # Sample rate and desired cut_off frequencies (in Hz).
  46. fs = 3000
  47. f1, f2 = 100, 200
  48. f3, f4 = 470, 750
  49. f5, f6 = 800, 850
  50. cut_off = [(f1, f2), (f3, f4), (f5, f6)]
  51. #cut_off = [(f1, f2), (f3, f4)]
  52. #cut_off = [(f1, f2)]
  53. # cut_off = [f1]
  54. t_toy_signal, toy_signal = get_toy_signal()
  55. # toy_signal -= np.mean(toy_signal)
  56. # t_toy_signal = wiener(t_toy_signal)
  57. fig, ax = plt.subplots(6, 1, figsize=(8, 12))
  58. fig.tight_layout()
  59. ax[0].plot(t_toy_signal, toy_signal)
  60. ax[0].set_title('Original toy_signal', fontsize=fontsize)
  61. ax[0].set_xlabel('Time (s)', fontsize=fontsize)
  62. ax[0].set_ylabel('Magnitude', fontsize=fontsize)
  63. ax[0].set_xlim(left=0, right=max(t_toy_signal))
  64. sos_lists=[]
  65. for cut in cut_off:
  66. print(cut)
  67. sos_list = [butter_bandpass(lowcut, highcut, fs, order=10) for lowcut, highcut in [cut]]
  68. sos_lists.append(np.vstack(sos_list))
  69. # Combine filters into a single filter
  70. sos = np.vstack(sos_list)
  71. # w *= 0.5 * fs / np.pi # Convert w to Hz.
  72. #####################################################################
  73. # First plot the desired ideal response as a green(ish) rectangle.
  74. #####################################################################
  75. # Plot the frequency response
  76. for i in range(len(cut_off)):
  77. w, h = sosfreqz(sos_lists[i], worN=2000)
  78. ax[1].plot(0.5 * fs * w / np.pi, np.abs(h), label=f'Band {i + 1}: {cut_off[i]} Hz')
  79. ax[1].set_title('Multiband Filter Frequency Response')
  80. ax[1].set_xlabel('Frequency [Hz]')
  81. ax[1].set_ylabel('Gain')
  82. ax[1].legend()
  83. # ax[1].set_xlim(0, max(*cut_off) + 100)
  84. #####################################################################
  85. # Spectrogram of original signal
  86. #####################################################################
  87. f, t, Sxx = spectrogram(toy_signal, fs,
  88. nperseg=930, noverlap=0)
  89. ax[2].pcolormesh(t, f, np.abs(Sxx),
  90. norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx), vmax=np.max(Sxx)),
  91. )
  92. ax[2].set_title('Spectrogram of original toy_signal', fontsize=fontsize)
  93. ax[2].set_xlabel('Time (s)', fontsize=fontsize)
  94. ax[2].set_ylabel('Frequency (Hz)', fontsize=fontsize)
  95. #####################################################################
  96. # Compute filtered signal
  97. #####################################################################
  98. # Apply the multiband filter to the data
  99. toy_signal_filtereds=[]
  100. for sos in sos_lists:
  101. toy_signal_filtered = sosfilt(sos, toy_signal, )
  102. toy_signal_filtereds.append([toy_signal_filtered])
  103. #####################################################################
  104. # Spectrogram of filtered signal
  105. #####################################################################
  106. f, t, Sxx = spectrogram(toy_signal_filtered, fs,
  107. nperseg=930, noverlap=0)
  108. ax[3].pcolormesh(t, f, np.abs(Sxx),
  109. norm=matplotlib.colors.LogNorm(vmin=np.min(Sxx),
  110. vmax=np.max(Sxx))
  111. )
  112. ax[3].set_title('Spectrogram of filtered toy_signal', fontsize=fontsize)
  113. ax[3].set_xlabel('Time (s)', fontsize=fontsize)
  114. ax[3].set_ylabel('Frequency (Hz)', fontsize=fontsize)
  115. for toy_signal_filtered in toy_signal_filtereds:
  116. ax[4].plot(t_toy_signal, toy_signal_filtered[0])
  117. ax[4].set_title('Filtered toy_signal', fontsize=fontsize)
  118. ax[4].set_xlim(left=0, right=max(t_toy_signal))
  119. ax[4].set_xlabel('Time (s)', fontsize=fontsize)
  120. ax[4].set_ylabel('Magnitude', fontsize=fontsize)
  121. N = 1512
  122. X = fft(toy_signal, n=N)
  123. Ys=[]
  124. for toy_signal_filtered in toy_signal_filtereds:
  125. Y = fft(toy_signal_filtered[0], n=N)
  126. Ys.append(Y)
  127. # fig.set_size_inches((10, 4))
  128. ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(X)), 'r-', label='FFT original signal')
  129. for bigY in Ys:
  130. ax[5].plot(np.arange(N) / N * fs, 20 * np.log10(abs(bigY)), 'g-', label='FFT filtered signal')
  131. ax[5].set_xlim(xmax=fs / 2)
  132. ax[5].set_ylim(ymin=-20)
  133. ax[5].set_ylabel(r'Power Spectrum (dB)', fontsize=fontsize)
  134. ax[5].set_xlabel("frequency (Hz)", fontsize=fontsize)
  135. ax[5].grid()
  136. ax[5].legend(loc=0)
  137. plt.tight_layout()
  138. plt.show()

字符串
对于频率响应,在你的情况下,你可能不应该使用对数标度图。但是ax[1].semilogy在y轴上使用对数标度图。


的数据

展开查看全部

相关问题