20220905-五阶巴特沃斯低通滤波器设计及实验-巴特沃思

1 巴特沃斯低通滤波器原理

巴特沃斯低通滤波器的传递函数如下

$H(p)=1b0+b1p+b2p2+\dots+bN-1pN-1+pNp=sωc$ H(p)=\frac{1}{b_0+b_1p+b_2p^2+…+b_{N-1}p^{N-1}+p^N}\\ p=\frac{s}{\omega_c}

其中，

N——滤波器的阶数

$ωc$ \omega_c ——滤波器的截止频率

$bi$ b_i ——滤波器的系数，可以查表确定

图 1.1 巴特沃斯低通滤波器系数表

2 五阶巴特沃斯低通滤波器

2.1 传递函数推导

根据第一节的知识，查表可以得到五阶巴特沃斯低通滤波器的传递函数如下

$G(p)=11+3.2361p+5.2361p2+5.2361p3+3.2361p4+p5$ G(p)=\frac{1}{1+3.2361p+5.2361p^2+5.2361p^3+3.2361p^4+p^5}

将 $p=sωc$ p=\frac{s}{\omega_c} 代入上式可得到模拟滤波器的传递函数

$G(s)=ω5ω5+3.2361ω4s+5.2361ω3s2+5.2361ω2s3+3.2361ωs4+s5$ G(s)=\frac{\omega^5}{\omega^5+3.2361\omega ^4s+5.2361\omega ^3s^2+5.2361\omega^2 s^3+3.2361\omega s^4+s^5}

2.2 频率特性分析

取截止频率为 1000Hz，则角频率为6283.2 rad/s编程绘制伯德图如下。由图 2.1 可知-3db 带宽为 6260 rad/s 左右，与设计结果基本一致。由图2.2 可知，在10倍频程处，衰减-100db，这说明五阶巴特沃斯低通滤波器的衰减速度是-100 / dB/dec ，频率每增大10 倍，幅值衰减100db，即幅值衰减为原来的10 万分之一，可见衰减速度还是很快的。

图 2.1 五阶巴特沃斯低通滤波器伯德图-3dB点

图 2.2 五阶巴特沃斯低通滤波器伯德图-100dB点

3 使用低通滤波器处理音频信号

使用上一节设计的低通滤波器对上一篇文章提到的音频进行滤波，滤波结果如下所示。图 3.1 中可以看到滤波后的信号幅值有所下降；图 3.2 可以看出滤波之后的信号的频谱强度从 1000Hz 开始衰减，到 2000Hz 的时候已经衰减到接近零，这说明滤波器上一节设计的低通滤波器是有效的。

图 3.1 滤波前后的音频信号曲线

图 3.2 滤波前后信号频谱

4 总结

巴特沃斯低通滤波器其实就是多个一阶低通滤波器的串联，通过查表的方法可以方便地得到归一化的表达式。

2. 根据要求的截止频率，可以算出模拟形式的低通滤波器的传递函数。

3. 将该滤波器通过双线性变换方法进行离散化，可以得到离散化的滤波器。

4. 使用离散化的滤波器。可以对采集的音频信号进行滤波。

5. 滤波之后，高频信号被衰减到几乎为零，而低频信号得以保留。

6. 播放音频发现滤掉高频信号之后，音乐中的背景噪声减弱了很多，但是音乐中的高频成分也没有了，感觉有点闷。

5 代码

% 绘制二阶巴特沃斯低通滤波器的频谱 % 20220904 clc;clear;close all; % ————<五阶巴特沃斯低通滤波器>————- w_c = 2*pi*1000; % 截止频率 rad/s num = [w_c^5]; % 滤波器系数 den = [1, 3.2361*w_c, 5.2361*w_c^2, 5.2361*w_c^3, 3.2361*w_c^4, w_c^5]; sys = tf(num, den) figure(1); % 绘制bode图 bode(sys); hold on; grid on; % ————-<读取音频信号>————— [x, f_s] = audioread(“1.m4a”); x = x(:,1); % 只取一个声道 len = length(x); t = (1:1:len)/f_s; % ——-<模拟滤波器的离散化>——— sys_2 =c2d(sys, 1/f_s, “tustin”); um_2 = sys_2.num{1}; den_2 = sys_2.den{1}; % ——–<对音频文件进行滤波>——- x_2 = filter(um_2, den_2, x); % sound(x, f_s); % 播放滤波前的音频 sound(x_2, f_s); % 播放滤波后的音频 % —-<绘制滤波前后的波形曲线进行对比>———- figure(2); plot(t, x, .-); grid on; hold on; plot(t, x_2, .-); xlabel(“t/s”); ylabel(“y”); set(gca, ylim, [-2,2]); legend(“原始信号”, “滤波后”); % ————<绘制滤波前后的频谱曲线进行对比>————- fft_x = fft(x, len); fft_x_2 = fft(x_2, len); f_x = (1:1:len)*f_s/len; figure(3); plot(f_x, abs(fft_x), .-); grid on; hold on; plot(f_x, abs(fft_x_2), .-); xlabel(“f/Hz”); ylabel(“强度”); legend(“原始信号”, “滤波后”); set(gca, “xlim”, [0,4000]);