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

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

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

极点形式： H ( s ) = G 0 ∏ k = 1 n 1 s − ω c s k ,   s k = exp ⁡ j ( 2 k + n − 1 ) 2 n ,   k = 1 , 2 , 3 , . . . , n H(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)=G0​k=1∏n​s−ωc​sk​1​, sk​=exp2nj(2k+n−1)​, k=1,2,3,...,n 其中， ω c \omega_c ωc​就是-3dB截止频率， n n n是滤波器阶数， G 0 G_0 G0​是直流增益，一般取1。 多项式形式： H ( s ) = G 0 ∑ k = 0 n a k ( s / ω c ) k ,   a k = ∏ μ = 1 k cos ⁡ ( ( μ − 1 ) γ ) sin ⁡ ( μ γ ) ,   a 0 = 1 ,   γ = π 2 n ,   k = 1 , 2 , 3 , . . . , n H(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=0n​ak​(s/ωc​)kG0​​, ak​=μ=1∏k​sin(μγ)cos((μ−1)γ)​, a0​=1, γ=2nπ​, k=1,2,3,...,n 对应于不同阶数的滤波器， s k s_k sk​和 a k a_k ak​的值都是可以事先计算好，然后查表得到的，如下示例：

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

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

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

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

计算阶数（向上取整）： N = lg ⁡ k s p lg ⁡ λ s p N=\frac{\lg{k_sp}}{\lg{\lambda_{sp}}} N=lgλsp​lgks​p​

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

如上文所示，从表格中找出对应阶数滤波器的系数值： a k a_k ak​ 3. 计算3 dB截止频率 ω c \omega_c ωc​ ω c = ω p ( 1 0 0.1 a p − 1 ) − 1 2 N \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 ) = G 0 ∑ k = 0 n a k ( s / ω c ) k H(s)=\frac{G_0}{\sum_{k=0}^n a_k({s}/{\omega_c})^k} H(s)=∑k=0n​ak​(s/ωc​)kG0​​

f a = f s π tan ⁡ π f d f s f_a=\frac{f_s}{\pi}\tan{\frac{\pi f_d}{f_s}} fa​=πfs​​tanfs​πfd​​ f s f_s fs​是采样频率， f d f_d fd​是数字滤波器截止频率

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

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

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

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

选定巴特沃斯滤波器的阶数，可以得到一个归一化的巴特沃斯低通滤波器，形式如下： H ( p ) = 1 ∑ k = 0 N a k p k H(p)=\frac{1}{\sum_{k=0}^N a_kp^k} H(p)=∑k=0N​ak​pk1​

针对不同的滤波器形式，利用图中公式计算出 Ω \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)公式中的 s s s用 z z z变量替换，得到 H ( z ) H(z) H(z)

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

给定条件：滤波器阶数为1，数字滤波器通带截止频率为30 Hz，采样频率为100 Hz 第一步：查表，得到归一化巴特沃斯低通滤波器形式： H ( p ) = 1 p + 1 H(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 π ∗ 30 100 = 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.7265 s + 0.7265 H(s)=\frac{0.7265}{s+0.7265} H(s)=s+0.72650.7265​ 第三步：将 H ( s ) H(s) H(s)公式中的 s s s用 z z z变量替换 H ( z ) = 0.7265 z − 1 z + 1 + 0.7265 = 0.7265 z − 0.7265 1.7265 z + 0.2735 = 0.4208 z − 0.4208 z + 0.1584 H(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] % 分母多项式系数

给定条件：滤波器阶数为2，数字滤波器上阻带截止频率为15 Hz，下阻带截止频率为10 Hz，采样频率为100 Hz 第一步：查表，得到归一化巴特沃斯低通滤波器形式： H ( p ) = 1 p 2 + 1.4142 p + 1 H(p)=\frac{1}{p^2+1.4142p+1} H(p)=p2+1.4142p+11​ 第二步：计算 ω \omega ω， ω s t 1 \omega_{st1} ωst1​， ω s t 2 \omega_{st2} ωst2​和 Ω \Omega Ω，在 H ( p ) H(p) H(p)中带入 p = s Ω p=\frac{s}{\Omega} p=Ωs​ ω s t 1 = 2 π ∗ 10 100 = 0.2 π ,   ω s t 2 = 2 π ∗ 15 100 = 0.3 π , cos ⁡ ω 0 = cos ⁡ ω s t 2 + ω s t 1 2 cos ⁡ ω s t 2 − ω s t 1 2 = cos ⁡ 0.3 π + 0.2 π 2 cos ⁡ 0.3 π − 0.2 π 2 = 0.7159 , Ω s t 1 = sin ⁡ ω s t 1 cos ⁡ ω s t 1 − cos ⁡ ω 0 = sin ⁡ 0.2 π cos ⁡ 0.2 π − 0.7159 = 6.3123 , Ω s t 2 = sin ⁡ ω s t 2 cos ⁡ ω s t 2 − cos ⁡ ω 0 = sin ⁡ 0.3 π cos ⁡ 0.3 π − 0.7159 = − 6.3148 , Ω = Ω s t 1 \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​−ωst1​​cos2ωst2​+ωst1​​​=cos20.3π−0.2π​cos20.3π+0.2π​​=0.7159,Ωst1​=cosωst1​−cosω0​sinωst1​​=cos0.2π−0.7159sin0.2π​=6.3123,Ωst2​=cosωst2​−cosω0​sinωst2​​=cos0.3π−0.7159sin0.3π​=−6.3148,Ω=Ωst1​​ H ( s ) = 6.312 3 2 s 2 + 6.3123 s + 6.312 3 2 = 39.8451 s 2 + 6.3123 s + 39.8451 H(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)公式中的 s s s用 z z z变量替换 H ( z ) = 39.8451 ( z 2 − 1 z 2 − 2 cos ⁡ ω 0 z + 1 ) 2 + 6.3123 z 2 − 1 z 2 − 2 cos ⁡ ω 0 z + 1 + 39.8451 = 39.8451 ( z 2 − 1.4318 z + 1 ) 2 ( z 2 − 1 ) 2 + 6.3123 ( z 2 − 1 ) ( z 2 − 1.4318 z + 1 ) + 39.8451 ( z 2 − 1.4318 z + 1 ) 2 = 39.8451 ( z 4 − 2.8636 z 3 + 4.0501 z 2 − 2.8636 z + 1 ) 47.1574 z 4 − 123.1384 z 3 + 159.3747 z 2 − 105.0625 z + 34.5328 = 0.8449 z 4 − 2.4196 z 3 + 3.4221 z 2 − 2.4196 z + 0.8449 z 4 − 2.6112 z 3 + 3.3796 z 2 − 2.2279 z + 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ω0​z+1z2−1​)2+6.3123z2−2cosω0​z+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] % 分母多项式系数

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