IIR Butterworth型模拟低通滤波器设计原理

Butterworth技术指标及要求

filter_design_001.png

$\omega p:表示通带截止频率$
$\omega s：表示阻带截止频率$
$\delta p:表示通带波动$
$\delta s:表示阻带波动$
衰减函数(dB): A( $\omega$ ) = -10 $lg|H(j\omega)|^2$ = -20 $lg|H(j\omega)|$
通带最大衰减(dB): Ap = -10 $lg|1-\delta p|^2$ = -20 $lg|1-\delta p|$
阻带最小衰减(dB):As = -10 $lg|\delta s|^2$ = -20 $lg|\delta s|$
Butterworth增益和衰减特性曲线如下：

filter_design_002.png

Butterworth模拟低通滤波器的频域特性

$|H(j\omega)|^2 =\frac{1}{1+(\omega/\omega c)^{2N}}$

filter_design_003.png

N:表示滤波器阶
$\omega c$ : 表示3dB截频
$|H(j0)|=1,|H(j \infty) = 0|$
$A(w) = -10lg|H(\omega c)|^2 = 10lg2 \approx 3dB$
所以 $\omega c$ 称为3dB 截频。若 $\omega c$ =1,则为归一化的BWF。
幅度响应单调下降。 $\omega c$ 一般不是通带截频 $\omega p$
$|H(j \omega)|^2$ 在 $\omega$ =0点1到2N-1阶导数为零,称为最大平坦性。

Butterworth模拟低通滤波器设计基本步骤

步骤1:确定模拟滤波器的阶数N
步骤2:确定模拟滤波器的3dB截频 $\omega_c$
步骤3:计算模拟滤波器的系统函数极点
步骤4:得到模拟低通滤波器的系统函数HL(s)

步骤一：确定模拟滤波器的阶数N

已知给定的条件 $(\omega_p,Ap)$ 、 $(\omega_s,As)$ ，根据 $A(w) = -10lg|H(\omega_c)|^2$ ,再由 $|H(j\omega)|^2 =\frac{1}{1+(\omega/\omega c)^{2N}}$ 可得到 $A(w) = -10lg\frac{1}{1+(\omega/\omega_c)^{2N}}$ ,推出：
$N \geq \frac{\lg[(10^{0.1As}-1)/10^{0.1Ap}-1)]}{2\lg[\omega_s/ \omega_p]}$

步骤二：确定模拟滤波器的3dB截频 $\omega_c$

根据以上公式可得出

$\frac{\omega_p}{(10^{0.1Ap - 1})^{1/2N}} \leq \omega_c \leq \frac{\omega_s}{(10^{0.1As - 1})^{1/2N}}$

filter_design_004.png

当 $\omega_c$ 由 $\omega_p$ 决定的时候所得到的衰减图像表明在通带刚好满足要求，在阻带有裕量
当 $\omega_c$ 由 $\omega_s$ 决定的时候所得到的衰减图像表明在阻带刚好满足要求，在通带有裕量
当 $\omega_c$ 在 $\omega_s$ 和 $\omega_p$ 之间任意取值时，最后得到的结果是在通带和阻带都有裕量

步骤3:确定模拟滤波器的系统函数极点

利用实系数模拟系统频率响应的共轭对称性有
$|H(j\omega)|^2 = H(j\omega)H^*(j\omega)=H(j\omega)H(-j\omega)=H(s)H(-s)|_{s=j\omega}$
求出系统函数的极点(位于s左半平面)如下:
$S_k = \omega_c * e^{j\pi(\frac{1}{2}+\frac{2k-1}{2N})}; k=1,2,3 \ldots$

步骤4:确定模拟低通滤波器的系统函数 $H_L(s)$

$H_L(s) = \prod_{k=1}^N\frac{-S_k}{S-S_k}$

filter_design_005.png

使用matlab求模拟低通滤波器的系统函数

连续系统的复频域描述如下,若描述连续LTI系统的微分方程为：
$a_ny^{n}(t) + a_{n-1}y^{n-1}(t) +\ldots+a_1y^{'}(t) + a_0y(t) = \\ b_mx^{m}(t) + b_{m-1}x^{m-1}(t) +\ldots+b_1x^{'}(t) + b_0x(t)$
利用Laplace变换微分特性,可得描述该系统的复频域方程
$[a_nS^{n} + a_{n-1}S^{n-1}+\ldots+a_1S+a_0]Y_{ZS} = \\ [b_mS^{m} + b_{m-1}S^{m-1}+\ldots+b_1S+b_0]X(s)$
最后得到连续时间LTI系统的系统函数如下：

$H(s) = \frac{Y_{ZS}}{X(s)} = \frac{b_mS^{m} + b_{m-1}S^{m-1}+\ldots+b_1S+b_0}{a_nS^{n} + a_{n-1}S^{n-1}+\ldots+a_1S+a_0}$

设计一个满足下列指标BW型模拟低通滤波器fp=1KHz , fs=2KHz,Ap ≤1dB, As ≥40dB,求出其系统函数，并画出其频率响应特性图

Wp=2*pi*1000;%通带截止频率,将数字频率化成模拟频率
Ws=2*pi*2000;%阻带截止频率,将数字频率化成模拟频率
Ap=1;%通带最大衰耗
As=40; %阻带最小衰耗
% 求得阶数和3db截频wc
[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %Computer the order of analog filter,s表示模拟
fprintf('Order of the filter=%.0f\n',N)
%num 表示系统函数分子系数,den表示系统函数分母系数
[num,den] = butter(N,Wc,'s'); %Compute AF coefficients
disp('Numerator polynomial'); fprintf('%.4e\n',num);
disp('Denominator polynomial'); fprintf('%.4e\n',den);

omega=[Wp Ws]; % 求出通带频率点和截止频率点对应的频率响应值
h = freqs(num,den,omega); %Compute Ap and As of AF

ap = -20*log10(abs(h(1)));% 将通带频率响应转换成衰耗值，
as = -20*log10(abs(h(2)));% 将阻带频率响应转换成衰耗值，

hp = abs(h(1));
hs = abs(h(2));

gain_p = 20*log10(abs(h(1))); 
gain_s = 20*log10(abs(h(2))); 

fprintf('Ap= %.4f\n',ap);%打印ap
fprintf('As= %.4f\n',as);%打印as

% 看0~6000pi之间，步进为200 的频率响应值
omega = [0: 200: 3000*2*pi];
h = freqs(num,den,omega); %Compute the frequency response of the AF
subplot(211);
plot(omega/(2*pi),abs(h));
hold on;%在画完函数虚线之后保持曲线图
plot([0 Wp/(2*pi)],[hp hp],'r--'); %画两个虚线
plot([Wp/(2*pi) Wp/(2*pi)],[0 hp],'r--');

plot([0 Ws/(2*pi)],[hs hs],'r--'); %画两个虚线
plot([Ws/(2*pi) Ws/(2*pi)],[0 hs],'r--');

hold off;
xlabel('Frequency in Hz'); 
ylabel('H(jw)');

gain=20*log10(abs(h)); 
subplot(212);
plot(omega/(2*pi),gain);
hold on;%在画完函数虚线之后保持曲线图
plot([0 Wp/(2*pi)],[gain_p gain_p],'r--'); %画两个虚线
plot([Wp/(2*pi) Wp/(2*pi)],[0 gain_p],'r--');

plot([0 Ws/(2*pi)],[gain_s gain_s],'r--'); %画两个虚线
plot([Ws/(2*pi) Ws/(2*pi)],[0 gain_s],'r--');

hold off;

xlabel('Frequency in Hz'); 
ylabel('Gain in dB');

最后得出其频率响应以及增益曲线图如下：

filter_design_006.png

其系数为

Order of the filter=8 %8阶滤波器
Numerator polynomial % 分之系数
0.0000e+00
0.0000e+00
0.0000e+00
0.0000e+00
0.0000e+00
0.0000e+00
0.0000e+00
0.0000e+00
6.2187e+30
Denominator polynomial %分母系数
1.0000e+00
3.6222e+04
6.5603e+08
7.7093e+12
6.4060e+16
3.8498e+20
1.6360e+24
4.5108e+27
6.2187e+30
Ap= 0.6167
As= 40.0000

$H(s) = \frac{6.2187*10^{30}S^{8}}{1+3.6222*10^{4}S^{1}+6.5603*10^{8}S^{2}+7.7093*10^{12}S^{3}+6.4060*10^{16}S^{4}+3.8498*10^{20}S^{5}+1.6360*10^{24}S^{6}+4.5108*10^{27}S^{7}+6.2187*10^{30}S^{8}}$

使用matlab分析该系统的零极点分布以及稳定性

num=[6.2187e+30];
den=[1.0000e+00 3.6222e+04 6.5603e+08 7.7093e+12 6.4060e+16 3.8498e+20 1.6360e+24 4.5108e+27 6.2187e+30];
subplot(211);
sys=tf(num,den);%求出系统的零极点
pzmap(sys);%画出系统零极点
subplot(212);
omega = [0: 200: 3000*2*pi];
H=freqs(num,den,omega);%得到频率响应
plot(omega/(2*pi),abs(H));%画出频率响应
xlabel('Frequency Hz');
title('Magnitude Respone');

filter_design_007.png

由图可知极点的分布都在S平面的左半平面，说明是系统稳定
N=8为偶数，一共有8个极点，和上面理论分析互相匹配