滤波的主要目的是减少噪声与干扰对数据的影响，让数据更加接近真实值。
比如常说的：
“高通滤波器”，就是高频通过，低频去掉的滤波器。
“低通滤波器”，就是低频通过，高频去掉的滤波器。
本文主要介绍： 高低通滤波器，三角滤波器，滤波器组 等。
# 准备 import matplotlib.pyplot as plt import numpy as np from scipy import signal # 定义一个方波 sr = 16_000 T = 4 x = np.linspace(0, T, T * sr, endpoint=False) y = 1 * np.sign(np.sin(2 * np.pi * 0.5 * x)) # 绘制 plt.figure(figsize=(8, 2)) plt.plot(y) plt.show()
高/低通滤波器
# 基准频率（HZ） cutoff_hz = 5 nyquist = 0.5 * sr normalized_cutoff = cutoff_hz / nyquist # 低通滤波，留下低于 cutoff_hz 的频率 sos_low = signal.butter( 8, normalized_cutoff, btype='low', output='sos' ) y_out_low = signal.sosfiltfilt(sos_low, y) print(f"sos_low.shape : {sos_low.shape}") print(f"y.shape : {y.shape}") print(f"y_out_low.shape : {y_out_low.shape}") # 高通滤波， 留下高于 cutoff_hz 的频率 sos_high = signal.butter( 8, normalized_cutoff, btype='high', output='sos' ) y_out_high = signal.sosfiltfilt(sos_high, y) print(f"sos_high.shape : {sos_high.shape}") print(f"y.shape : {y.shape}") print(f"y_out_high.shape : {y_out_high.shape}")
sos_low.shape : (4, 6)
y.shape : (64000,)
y_out_low.shape : (64000,)
sos_high.shape : (4, 6)
y.shape : (64000,)
y_out_high.shape : (64000,)
# 绘图对比 _, axs = plt.subplots(4, 1) axs[0].plot(y) axs[1].plot(y_out_low) axs[2].plot(y_out_high) axs[3].plot(y_out_low + y_out_high) plt.tight_layout() plt.show()
mse = np.square(y_out_low + y_out_high - y).mean() mse
4.5419663300030064e-07
如上对比可知： 高通滤波剩余信号 + 低通滤波剩余信号 = 原始信号
三角滤波器
这里的实验：
- 利用快速傅里叶变换提取方波主要频率
- 然后自定义一个三角滤波器
- 让方波主要频率通过三角滤波器
- 最后对比三角滤波器的滤波效果。
主要频率提取
# 信号长度 n = len(y) # 傅里叶变换 y_fft = np.fft.fft(y) y_freq = np.fft.fftfreq(n, 1 / sr) y_fft = y_fft[:n//2] y_freq = y_freq[:n//2] # 信号能量 magnitude = 2.0 * np.abs(y_fft) / n
# 提取主要频率 threshold = 0.1 mask = magnitude > threshold detected_freqs = y_freq[mask] detected_amps = magnitude[mask] # 绘制频率图 plt.figure(figsize=(5, 3)) plt.plot(detected_freqs, detected_amps) for xi, yi in zip(detected_freqs, detected_amps): plt.text(xi, yi, f"({xi},{yi:.2f})", ha='center', va='bottom') plt.grid(True) plt.show()
三角滤波器
# 三角滤波器生成方法 def triangular_filter(freqs, left, center, right): filter_response = np.zeros_like(freqs) # 上升部分 rising = (freqs >= left) & (freqs < center) filter_response[rising] = (freqs[rising] - left) / (center - left) # 下降部分 falling = (freqs >= center) & (freqs <= right) filter_response[falling] = (right - freqs[falling]) / (right - center) return filter_response # 创建一个三角滤波器 frequencies = np.linspace(0, max(detected_freqs), len(detected_freqs)) print(f"frequencies : {frequencies}") filter_response = triangular_filter(frequencies, 0, 2.2, 4.4) print(f"filter_response : {filter_response}") plt.plot(frequencies, filter_response) plt.show()
frequencies : [0. 1.1 2.2 3.3 4.4 5.5]
filter_response : [0. 0.5 1. 0.5 0. 0. ]
滤波效果对比
# 滤波操作：过滤频率 filter_freqs = detected_freqs * filter_response print(f"detected_freqs : {detected_freqs}") print(f"filter_response : {filter_response}") print(f"filter_freqs : {filter_freqs}")
detected_freqs : [0.5 1.5 2.5 3.5 4.5 5.5]
filter_response : [0. 0.5 1. 0.5 0. 0. ]
filter_freqs : [0. 0.75 2.5 1.75 0. 0. ]
# fft主要频率重建 re_wave = np.zeros_like(y) for freq, amp in zip(detected_freqs, detected_amps): re_wave += amp * np.sin(2 * np.pi * freq * x) # fft主要频率 进过三角滤波器之后，重建 re_filter_wave = np.zeros_like(y) for freq, amp in zip(filter_freqs, detected_amps): re_filter_wave += amp * np.sin(2 * np.pi * freq * x) # 绘图对比 plt.figure(figsize=(9, 3)) plt.plot(x, y, label='y') plt.plot(x, re_wave, label='Rebuilded Wave') plt.plot(x, re_filter_wave, label='Rebuilded Filter Wave') plt.tight_layout() plt.grid(True) plt.legend() plt.show()
好吧，这个图，似乎看不出滤波器什么的内容。
不过（滤波操作：过滤频率）的数字，可以仔细看看。
filter_freqs = detected_freqs * filter_response
滤波器组
单纯的 filter_response 是一个滤波器，一般只对一个小范围的频率进行滤波。
可以组建滤波器组，对所有范围的频率进行过滤。
# 基本参数 sample_rate = 16000 n_fft = 512 n_mels = 64 f_min = 0.0 f_max = sample_rate / 2.0 n_freqs = int(n_fft // 2 + 1) # 所有频率 all_freqs = np.linspace(0, sample_rate / 2.0, n_freqs) # 目标刻度 t_pts = np.linspace(f_min, f_max, n_mels + 2) # 创建三角滤波器组 def create_triangular_fbank(all_freqs, f_pts): # 滤波器 filter_fbank = np.zeros((len(all_freqs), len(f_pts) - 2)) # 实现三角滤波器 for i in range(n_mels): left = t_pts[i] center = t_pts[i+1] right = t_pts[i+2] # up up = (all_freqs >= left) & (all_freqs < center) filter_fbank[up, i] = (all_freqs[up] - left) / (center - left) # down down = (all_freqs >= center) & (all_freqs <= right) filter_fbank[down, i] = (right - all_freqs[down]) / (right - center) return filter_fbank # 创建滤波器组 filter_fbank = create_triangular_fbank(all_freqs, t_pts) print(filter_fbank.T.shape)
(64, 257)
# 绘制滤波器 plt.imshow(filter_fbank[:, :10].T, aspect='auto') plt.show()
奇怪，三角滤波器组，的三角去哪儿了？
好吧，图应该这样画：
# 绘制三角滤波器组 fig = plt.figure() ax = fig.add_subplot(111, projection='3d') for i in range(filter_fbank.shape[1]): x = np.full(filter_fbank.shape[0], i) y = np.arange(filter_fbank.shape[0]) z = filter_fbank[:, i] # 绘制 ax.plot(x, y, z) # 只绘制10条：filter_fbank[:, :10].T if i == 10: break ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') ax.view_init(elev=20, azim=-30) plt.show()
