巴特沃斯数字低通滤波器的设计步骤

Hello，欢迎来到我的博客~~~

1 低通滤波器性能指标

低通滤波器的主要性能指标有四个：

通带截止频率：通带内幅值下降至一定程度时对应的频率，单位为：Hz或角频率
阻带截止频率：阻带内幅值高于一定程度时对应的频率，单位为：Hz或角频率
通带纹波（衰减）：通带中最大幅值和最小幅值之间的差值，单位为：dB
阻带衰减：阻带中最大幅值和通带最大幅值的差值，单位为：dB

幅值单位：dB，计算方式：20log⁡10A20\log_{10}A20log10A
示例：

假如幅值响应是100，那么相当于20log⁡10100=4020\log_{10}100=4020log10100=40 dB

假如通带最大幅值是100，阻带最大幅值是10，那么阻带衰减就是20log⁡1010−20log⁡10100=−2020\log_{10}10 - 20\log_{10}100=-2020log1010−20log10100=−20 dB

2 巴特沃斯模拟低通滤波器设计步骤

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

极点形式：
H(s)=G0∏k=1n1s−ωcsk,sk=exp⁡j(2k+n−1)2n,k=1,2,3,...,nH(s)=G_0\prod_{k=1}^n \frac{1}{s-\omega_cs_k},\ s_k=\exp\frac{j(2k+n-1)}{2n}, \ k=1,2,3,...,n H(s)=G0k=1∏ns−ωcsk1, sk=exp2nj(2k+n−1), k=1,2,3,...,n
其中，ωc\omega_cωc就是-3dB截止频率，nnn是滤波器阶数，G0G_0G0是直流增益，一般取1。
多项式形式：
H(s)=G0∑k=0nak(s/ωc)k,ak=∏μ=1kcos⁡((μ−1)γ)sin⁡(μγ),a0=1,γ=π2n,k=1,2,3,...,nH(s)=\frac{G_0}{\sum_{k=0}^n a_k({s}/{\omega_c})^k} ,\ a_k=\prod_{\mu=1}^{k}\frac{\cos((\mu-1)\gamma)}{\sin(\mu\gamma)},\ a_0=1, \ \gamma=\frac{\pi}{2n}, \ k=1,2,3,...,n H(s)=∑k=0nak(s/ωc)kG0, ak=μ=1∏ksin(μγ)cos((μ−1)γ), a0=1, γ=2nπ, k=1,2,3,...,n
对应于不同阶数的滤波器，sks_ksk和aka_kak的值都是可以事先计算好，然后查表得到的，如下示例：

2.2 频率响应

根据巴特沃斯滤波器的传递函数，可以推得其频率响应为：
∣H(jω)∣2=G021+(ωωc)2n|H(j\omega)|^2=\frac{G_0^2}{1+(\frac{\omega}{\omega_c})^{2n}} ∣H(jω)∣2=1+(ωcω)2nG02

根据上述的频率响应，很容易分析得到当ω\omegaω趋于0时，频率响应趋于G02G_0^2G02，当ω\omegaω趋于无穷时，频率响应趋于0。

2.3 设计步骤

1. 给定通带截止频率ωp\omega_pωp，阻带截止频率ωs\omega_sωs，通带纹波αs\alpha_sαs和阻带衰减αp\alpha_pαp，计算滤波器阶数N

计算两个辅助变量：ksp=100.1αs−1100.1αp−1k_{sp}=\sqrt{\frac{10^{0.1\alpha_s}-1}{10^{0.1\alpha_p}-1}}ksp=100.1αp−1100.1αs−1, λsp=ωsωp\lambda_{sp}=\frac{\omega_s}{\omega_p}λsp=ωpωs

计算阶数（向上取整）：N=lg⁡ksplg⁡λspN=\frac{\lg{k_sp}}{\lg{\lambda_{sp}}}N=lgλsplgksp

2. 根据滤波器阶数，通过查表得到滤波器传递函数的系数

如上文所示，从表格中找出对应阶数滤波器的系数值：aka_kak
3. 计算3 dB截止频率ωc\omega_cωc
ωc=ωp(100.1ap−1)−12N\omega_c=\omega_p(10^{0.1a_p}-1)^{-\frac{1}{2N}}ωc=ωp(100.1ap−1)−2N1

4. 代入系数和ωc\omega_cωc，得到最终的滤波器传递函数

H(s)=G0∑k=0nak(s/ωc)kH(s)=\frac{G_0}{\sum_{k=0}^n a_k({s}/{\omega_c})^k}H(s)=∑k=0nak(s/ωc)kG0

3 巴特沃斯数字低通滤波器设计步骤（IIR实现）

选择一个归一化的模拟滤波器（确定巴特沃斯低通滤波器的阶数）
确定数字滤波器的3 dB截止频率
利用公式计算模拟滤波器的3 dB截止频率

fa=fsπtan⁡πfdfsf_a=\frac{f_s}{\pi}\tan{\frac{\pi f_d}{f_s}}fa=πfstanfsπfd
fsf_sfs是采样频率，fdf_dfd是数字滤波器截止频率

将模拟截止频率ωc=2πfa\omega_c=2\pi f_aωc=2πfa带入模拟滤波器传递函数H(s)H(s)H(s)
用双线性变换，把模拟滤波器传递函数中的sss替换为zzz，得到H(z)H(z)H(z)

双线性变换：s=2fs(z−1z+1)s=2f_s(\frac{z-1}{z+1})s=2fs(z+1z−1)

4 巴特沃斯高通、带通、带阻数字滤波器的设计

要设计高通、带通、带阻等数字滤波器，有两种思路。

低通模拟滤波器 =》=》高通、带通、带阻模拟滤波器 =》=》高通、带通、带阻数字滤波器
低通模拟滤波器 =》=》高通、带通、带阻数字滤波器

这里主要介绍的是第二种思想，方法如下图所示：

4.1 变量说明

ω=2πfpfs\omega=\frac{2\pi f_p}{f_s}ω=fs2πfp，其中，fpf_pfp是数字通带截止频率
Ω=2πFpfs\Omega=\frac{2\pi F_p}{f_s}Ω=fs2πFp，其中，FpF_pFp是模拟通带截止频率（注意这里的符号和上文有差异，不要混淆）
ωp1=2πfp1fs\omega_{p1}=\frac{2\pi f_{p1}}{f_s}ωp1=fs2πfp1，其中，fp1f_{p1}fp1是数字下通带截止频率（带通滤波器）
ωp2=2πfp2fs\omega_{p2}=\frac{2\pi f_{p2}}{f_s}ωp2=fs2πfp2，其中，fp2f_{p2}fp2是数字上通带截止频率（带通滤波器）
ωst1=2πfst1fs\omega_{st1}=\frac{2\pi f_{st1}}{f_s}ωst1=fs2πfst1，其中，fst1f_{st1}fst1是数字下阻带截止频率（带阻滤波器）
ωst2=2πfst2fs\omega_{st2}=\frac{2\pi f_{st2}}{f_s}ωst2=fs2πfst2，其中，fst2f_{st2}fst2是数字上阻带截止频率（带阻滤波器）

4.2 设计步骤

根据上图，设计步骤可以描述如下：

选定巴特沃斯滤波器的阶数，可以得到一个归一化的巴特沃斯低通滤波器，形式如下：
H(p)=1∑k=0NakpkH(p)=\frac{1}{\sum_{k=0}^N a_kp^k} H(p)=∑k=0Nakpk1
针对不同的滤波器形式，利用图中公式计算出Ω\OmegaΩ，并在H(p)H(p)H(p)中带入p=sΩp=\frac{s}{\Omega}p=Ωs，得到H(s)H(s)H(s)
针对不同的滤波器形式，利用图中公式，将H(s)H(s)H(s)公式中的sss用zzz变量替换，得到H(z)H(z)H(z)

注意事项：
上图中对应低通、高通、带通、带阻都有sss和Ω\OmegaΩ的计算方法
需要注意的是，
对于低通和高通而言，ω\omegaω一般指的都是通带333 dB截止频率
对于带通和带阻而言，ω\omegaω一般指的都是上通带或上阻带333 dB截止频率

4.3 设计示例

4.3.1 巴特沃斯数字高通滤波器

给定条件：滤波器阶数为1，数字滤波器通带截止频率为30 Hz，采样频率为100 Hz
第一步：查表，得到归一化巴特沃斯低通滤波器形式：
H(p)=1p+1H(p)=\frac{1}{p+1} H(p)=p+11
第二步：计算ω\omegaω和Ω\OmegaΩ，在H(p)H(p)H(p)中带入p=sΩp=\frac{s}{\Omega}p=Ωs
ω=2π∗30100=0.6π,Ω=cot⁡ω2=0.7265\omega=\frac{2\pi*30}{100}=0.6\pi,\ \Omega=\cot{\frac{\omega}{2}}=0.7265 ω=1002π∗30=0.6π, Ω=cot2ω=0.7265
H(s)=0.7265s+0.7265H(s)=\frac{0.7265}{s+0.7265} H(s)=s+0.72650.7265
第三步：将H(s)H(s)H(s)公式中的sss用zzz变量替换
H(z)=0.7265z−1z+1+0.7265=0.7265z−0.72651.7265z+0.2735=0.4208z−0.4208z+0.1584H(z)=\frac{0.7265}{\frac{z-1}{z+1}+0.7265}=\frac{0.7265z-0.7265}{1.7265z+0.2735}=\frac{0.4208z-0.4208}{z+0.1584} H(z)=z+1z−1+0.72650.7265=1.7265z+0.27350.7265z−0.7265=z+0.15840.4208z−0.4208
以上是设计得到的巴特沃斯数字高通滤波器，利用Matlab可以验算结果，和公式计算的完全一致。

Matlab命令：
fc = 30; fs = 100; [a,b]=butter(1, fc/(fs/2), ‘high’)
结果：
a = [0.4208, -0.4208] % 分子多项式系数
b = [1.0000, 0.1584] % 分母多项式系数

4.3.2 巴特沃斯数字带阻滤波器

给定条件：滤波器阶数为2，数字滤波器上阻带截止频率为15 Hz，下阻带截止频率为10 Hz，采样频率为100 Hz
第一步：查表，得到归一化巴特沃斯低通滤波器形式：
H(p)=1p2+1.4142p+1H(p)=\frac{1}{p^2+1.4142p+1} H(p)=p2+1.4142p+11
第二步：计算ω\omegaω，ωst1\omega_{st1}ωst1，ωst2\omega_{st2}ωst2和Ω\OmegaΩ，在H(p)H(p)H(p)中带入p=sΩp=\frac{s}{\Omega}p=Ωs
ωst1=2π∗10100=0.2π,ωst2=2π∗15100=0.3π,cos⁡ω0=cos⁡ωst2+ωst12cos⁡ωst2−ωst12=cos⁡0.3π+0.2π2cos⁡0.3π−0.2π2=0.7159,Ωst1=sin⁡ωst1cos⁡ωst1−cos⁡ω0=sin⁡0.2πcos⁡0.2π−0.7159=6.3123,Ωst2=sin⁡ωst2cos⁡ωst2−cos⁡ω0=sin⁡0.3πcos⁡0.3π−0.7159=−6.3148,Ω=Ωst1\begin{aligned} &\omega_{st1}=\frac{2\pi*10}{100}=0.2\pi, \ \omega_{st2}=\frac{2\pi*15}{100}=0.3\pi,\\ &\cos\omega_0=\frac{\cos\frac{\omega_{st2}+\omega_{st1}}{2}}{\cos\frac{\omega_{st2}-\omega_{st1}}{2}}=\frac{\cos\frac{0.3\pi+0.2\pi}{2}}{\cos\frac{0.3\pi-0.2\pi}{2}}=0.7159,\\ &\Omega_{st1}=\frac{\sin\omega_{st1}}{\cos\omega_{st1}-\cos\omega_0}=\frac{\sin0.2\pi}{\cos0.2\pi-0.7159}=6.3123,\\ &\Omega_{st2}=\frac{\sin\omega_{st2}}{\cos\omega_{st2}-\cos\omega_0}=\frac{\sin0.3\pi}{\cos0.3\pi-0.7159}=-6.3148,\\ &\Omega = \Omega_{st1} \end{aligned} ωst1=1002π∗10=0.2π, ωst2=1002π∗15=0.3π,cosω0=cos2ωst2−ωst1cos2ωst2+ωst1=cos20.3π−0.2πcos20.3π+0.2π=0.7159,Ωst1=cosωst1−cosω0sinωst1=cos0.2π−0.7159sin0.2π=6.3123,Ωst2=cosωst2−cosω0sinωst2=cos0.3π−0.7159sin0.3π=−6.3148,Ω=Ωst1
H(s)=6.31232s2+6.3123s+6.31232=39.8451s2+6.3123s+39.8451H(s)=\frac{6.3123^2}{s^2+6.3123s+6.3123^2}=\frac{39.8451}{s^2+6.3123s+39.8451} H(s)=s2+6.3123s+6.312326.31232=s2+6.3123s+39.845139.8451
第三步：将H(s)H(s)H(s)公式中的sss用zzz变量替换
H(z)=39.8451(z2−1z2−2cos⁡ω0z+1)2+6.3123z2−1z2−2cos⁡ω0z+1+39.8451=39.8451(z2−1.4318z+1)2(z2−1)2+6.3123(z2−1)(z2−1.4318z+1)+39.8451(z2−1.4318z+1)2=39.8451(z4−2.8636z3+4.0501z2−2.8636z+1)47.1574z4−123.1384z3+159.3747z2−105.0625z+34.5328=0.8449z4−2.4196z3+3.4221z2−2.4196z+0.8449z4−2.6112z3+3.3796z2−2.2279z+0.7323\begin{aligned} H(z)&=\frac{39.8451}{(\frac{z^2-1}{z^2-2\cos\omega_0z+1})^2+6.3123\frac{z^2-1}{z^2-2\cos\omega_0z+1}+39.8451}\\ &=\frac{39.8451(z^2-1.4318z+1)^2}{(z^2-1)^2+6.3123(z^2-1)(z^2-1.4318z+1)+39.8451(z^2-1.4318z+1)^2}\\ &=\frac{39.8451(z^4 - 2.8636z^3 + 4.0501z^2 - 2.8636z + 1)}{47.1574z^4 - 123.1384z^3 + 159.3747z^2 - 105.0625z+ 34.5328}\\ &=\frac{0.8449z^4 - 2.4196z^3 + 3.4221z^2 - 2.4196z + 0.8449}{z^4 - 2.6112z^3 + 3.3796z^2 - 2.2279z+ 0.7323}\\ \end{aligned} H(z)=(z2−2cosω0z+1z2−1)2+6.3123z2−2cosω0z+1z2−1+39.845139.8451=(z2−1)2+6.3123(z2−1)(z2−1.4318z+1)+39.8451(z2−1.4318z+1)239.8451(z2−1.4318z+1)2=47.1574z4−123.1384z3+159.3747z2−105.0625z+34.532839.8451(z4−2.8636z3+4.0501z2−2.8636z+1)=z4−2.6112z3+3.3796z2−2.2279z+0.73230.8449z4−2.4196z3+3.4221z2−2.4196z+0.8449
以上是设计得到的巴特沃斯数字高通滤波器，利用Matlab可以验算结果，和公式计算的基本一致。

Matlab命令：
fst1 = 10; fst2 = 15; fs = 100; [a,b]=butter(2,[fst1/(fs/2),fst2/(fs/2)],‘stop’)
结果：
a = [0.8006, -2.2926, 3.2425, -2.2926, 0.8006] % 分子多项式系数
b = [1.0000, -2.5494, 3.2024, -2.0359, 0.6414] % 分母多项式系数

5 参考资料

参考链接：从模拟滤波器到数字滤波器
数字信号处理公式变程序（四）——巴特沃斯滤波器（上）
数字信号处理教程（超浓缩版）