本篇文章转载自博客园MrLouis

IIR数字滤波器的设计

冲激响应不变法

  • 冲激响应不变法：就是用其单位冲激响应序列模仿模拟滤波器的单位冲激响应的抽样值

  • 设计的具体步骤及方法

    首先要设计一个响应的模拟滤波器。
    模拟滤波器的单位冲激响应非周期$\Rightarrow$频率离散$(-\infty,+\infty)$，频谱范围可能大于折叠频率，即奈奎斯特频率 $\Rightarrow$取样，频率周期延拓$\Rightarrow$频谱可能产生混叠。
    数字域中极点在单位圆内$\Leftrightarrow$模拟域中极点在左半平面。
    可以通过提高抽样频率来减少混叠，但设计指标若以数字域频率给定时，不能通过提高抽样频率改善混叠现象。因为$\Omega=\frac{\omega}{T}$，$\omega$是数字域频率，$\Omega$是模拟域频率，增大抽样频率$\Rightarrow$减小采样周期$\Rightarrow$增大截止频率。

    模拟域 ：$\sum_{k=1}^N \frac{A_k}{s-s_k}$

    数字域：$\sum_{k=1}^N \frac{A_k}{1-e^{s_kT}z^{-1}}$
    s平面：$s=s_k$

    z平面：$e^{s_kT}$两者系数相同

    • 当采样周期很小时，数字滤波器增益很大，容易溢出，可以修正，令$h(n)=h_a(nT)$，因为$\frac{1}{T}H_a(j\frac{\omega}{T})$
    • 该方法对应的MATLAB程序：$[BZ, AZ] = impinvar(B,A,Fs)$，其中BZ，AZ为数字滤波器单位冲击响应的分子和分母系数；B，A为模拟滤波器的分子和分母的系数，Fs为抽样频率。

  • 优点

    • 时域逼近良好
    • 保持线性关系：$w=\Omega{T}$

  • 缺点

    • 频率响应混叠，只适用低通、带通滤波器。

阶跃响应不变法

  • 变换原理：数字滤波器的阶跃响应模仿模拟滤波器的阶跃响应，即$g(n)=g_a(t)|t=nt=g_a(nT)$.
  • $g(n)=u(n)*h(n)→G(z)=\frac{z}{z-1}H(z)$ ,$g_a(t)=u(t)h_a(t)→G_a(s)=\frac{1}{s}H_a(s)$
    变换步骤：$H_a(s)\rightarrow G_a(s)\rightarrow g_a(t)\rightarrow g(n)\rightarrow G(z)\rightarrow H(z)$
  • 优缺点同冲激响应不变法

双线性不变法

  • 为什么需要双线性变换法？—避免频谱混叠

    在推导过程中有一个中间域，最终的关系式为$z=\frac{1+s}{1-s}$ ,为了使模拟滤波器与数字滤波器的某一频率有对应关系，引入一个常数$c$，则有$s=c\frac{1-\frac{1}{z}}{1+\frac{1}{z}}$ ,即$\frac{c+s}{c-s}$。若低频处有确切的对应关系应为$c=\frac{2}{T}$。
    s平面虚轴与z平面单位圆的对应关系$\Omega=ctan\frac{w}{2}\Omega$为$s$域，$w$为数字域

  • 优点：无频率混叠，避免了多值映射。

  • 缺点：除零频率附近，严重非线性关系，在临界频率点（通带截止频率和阻带截止频率—这个是技术指标)产生畸变。
  • 解决临界频率点的畸变—预畸变
    $\Omega=ctan\frac{w}{2}\Omega$ ,$c$一般取$\frac{2}{T}$，在临界频率点根据$w$去设计$Ω$ ,得到模拟滤波器对应的技术指标，再去设计模拟滤波器，再用双线性变换法设计数字滤波器。
  • 该方法对应的MATLAB方法：$[BZ, AZ] = bilinear(B,A,Fs)$

常用模拟滤波器设计

  • 由幅度平方函数确定模拟滤波器的系统函数
    $|H_a(j\Omega)|^2=H_a(j\Omega)H_a^*(j\Omega)=H_a(j\Omega)H_a(-j\Omega)=H_a(s)H_a(-s)|(s=j\Omega)，H_a(jΩ)$可以看做模拟滤波器单位冲激响应的傅里叶变化。
    $H_a(s)$和$H_a(−s)$的零极点对称，成对出现，且其零极点分别分布在左半平面（只有分布在左半平面系统才是稳定的）和右半平面。

巴特沃斯低通逼近

幅度平方函数为$|H_a(j\Omega)|^2 = \frac{1}{1+(\frac{\Omega}{\Omega_c})^{2N}}$，其中N为滤波器阶数，$Ω_c$为3dB截止频率，因为$\Omega=\Omega_c$时，幅度为原来的$\frac{1}{\sqrt{2}}$，

img

  • 极点分布

$|H_a(j\Omega)|{\Omega=\frac{s}{j}}^2 = H_a(s)H_a(-s) = \frac{1}{1+(\frac{s}{j\Omega})^{2N}}$，令$1+(\frac{s}{j\Omega})^{2N}=0$即可解得极点分布。

极点在s平面呈象限对称，个数为2N，极点间角度间隔为$\frac{\pi}{N}$，极点不落在虚轴上，若N为奇数，实轴上有极点；若为偶数，实轴上无极点。

  • 滤波器的系统函数
  • 如何确定N和$\Omega_c$

已知技术指标：$Ω_c$(通带截止频率)、$δ1$(通带内允许的最大衰减)、$Ω_s$(阻带截止频率)、$δ_2$(阻带允许的最小衰减).
$\delta_1=-20lg|H_a(j\Omega)|,|H_a(j\Omega)|^2=\frac{1}{1+(\frac{\Omega_p}{\Omega_c})^{2N}}\rightarrow1+(\frac{\Omega_p}{\Omega_c})^{2N}=10^{0.1\delta_1}\rightarrow N$①
同理,有$1+(\frac{\Omega_s}{\Omega_c})^{2N}=10^{0.1\delta_2}$② 由这两个方程可以解得N和$\Omega_c$ 。

$Ω_c$既可以若由第一个方程解得，那么等同于将通带的技术指标卡死，表示阻带指标有富余，若由第二个方程解得，则与前一种情况相反。

  • MATLAB仿真
1
2
3
4
5
[N, Wc]=buttord(wp,ws,Rp,Rs,'s') %根据性能指标得到巴特沃斯滤波器的阶数及截止频率
[z,p,k]=buttap(N) %返回N阶归一化原型模拟滤波器的零极点以及增益
[B,A]=zp2tf(z,p,k) %根据零极点增益得到原型多项式的系数（归一化）
[Bt,At]=lp2lp(B,A,Wo) %得到模拟低通滤波器的多项式系数（去归一化）
[BZ,AZ]=bilinear(Bt,At,Fs) %使用双线性变换法求得数字低通滤波器额的多项式系数

切比雪夫低通逼近

  • 类型
    I型：通带范围内波动；II型：阻带范围内波动。
    相比于巴特沃斯：巴特沃斯在整个频带内都是下降的，卡死通带或阻带中其中一个，都将引起另一个的快速变化。

  • 幅度平方函数（以I型为例）
    $|H_a(j\Omega)|^2=\frac{1}{1+\varepsilon^2C_N^2(\frac{\Omega}{\Omega_c})}$
    $ε$：通带波纹大小，$0<ε<1$值越大，波纹越大
    $Ω_c$：通带截止频率，不一定为3dB带宽
    $N$：滤波器阶数
    $C_N(x)$：N阶切比雪夫多项式

$C_N(x)=\begin{cases}\cos(N\cos^{-1}x)& \text{|x|<=1}\ ch(Nch^{-1}x)&|x|>1\end{cases}$

该多项式的图像如下：

从图像可以看出，n为奇数时，$|C_N(0)|=0$；n为偶数时，$|C_N(0)|=1$。$\Omega=\Omega_c时，|C_N(x)|=1，|H_a(j\Omega)|=\frac{1}{\sqrt{1+\varepsilon^2}}$$在|x|<=1时，|C_N(x)|在0-1$间反复变化，导致$|H_a(j\Omega)|在\frac{1}{\sqrt{1+\varepsilon^2}}$之间反复变化，从而形成波纹$ 。$而$\delta_1$表示通带内允许的最大衰减，故应有$\delta_1=20lg\sqrt{1+\varepsilon^2}，$从而可以解得参数$\varepsilon。$ 再将阻带允许的最小衰减和阻带截止频率带入即可得到阶数N。

极点个数为2N个，分布在椭圆上。

  • 设计步骤

  • MATLAB仿真

    1
    2
    3
    4
    5
    [N,Wc]=cheblord(wp,ws,Rp,Rs,'s') %根据性能指标求出切比雪夫滤波器的阶数和截止频率
    [z,p,k]=cheb1ap(N,rp) %返回N阶归一化原型模拟滤波器的零极点及增益
    [B,A]=zp2tf(z,p,k) %求出归一化低通原型多项式系数
    [Bt,At]=lp2lp(B,A,Wo) %求出模拟低通滤波器的系数
    [BZ,AZ]=impinvar(Bt,At,Fs) %冲激响应不变法求数字低通滤波器的系数

频率变换法

  • 为什么需要频率变换法？
    之间的设计方法都是基于低通滤波器进行设计的，为了得到其他类型的滤波器，需要进行频率变换

模拟域频带变换法

  • 变换步骤
    归一化模拟低通滤波器$⇒$模拟域频带变换$⇒$模拟低通、高通、带通、带阻滤波器$⇒$双线性变换$⇒$数字低通、高通、带通、带阻滤波器。
    总结来讲，就是在模拟域就进行频带变换。

  • 具体方法

    归一化低通滤波器的系统函数：$H_a(s),s=\sigma+j\Omega$

    频带变换后的系统函数：$\overline{H}_s(\overline{s}),\overline{s}=\overline{\sigma}+j\overline{\Omega}$

    • 模拟低通→→模拟低通
      $s=\frac{\Omega_c\overline{s}}{\overline{\Omega_c}}$，当$\Omega_c=1$时，相当于去归一化。
      对应的MATLAB函数为[Bt,At]=lp2lp(B,A,Wo)
    • 模拟低通→→模拟带通
      $s=\frac{\overline{s}^2+\overline{\Omega}_0^2}{\overline{s}}$，$\overline{\Omega}_0$为带通滤波器的中心频率$\overline{\Omega}_0=\overline{\Omega}_1\overline{\Omega}_2$。$\overline{\Omega}_1、\overline{\Omega}_2$分别为带通滤波器的上通带截止频率和下通带截止频率。为了得到归一化的带通滤波器，可令$s=\frac{\overline{s}^2+\overline{\Omega}_0^2}{\overline{s}}$两边同时除以$\Omega_c(B)$.
      对应的MATLAB函数为[Bt,At]=lp2bp(B,A,Wo,Bw)，Bw为带宽。
      数字带通滤波器性能指标→→模拟带通滤波器性能能指标$→$归一化$→$模拟低通滤波器归一化性能指标→→模拟低通滤波器。
      带通滤波器归一化频率计算方法：$\eta=\frac{\overline{\Omega}}{B}$
      低通滤波器归一化频率计算方法：$\lambda=\frac{\Omega}{\Omega_c},\Omega_c$为通带截止频率。
      两者间关系：$\lambda=\frac{\eta^2-\eta_0^2}{\eta}$。

数字域频带变换法

  • 变换步骤
    归一化模拟低通滤波器$⇒$数字低通滤波器$⇒$数字域频带变换$⇒$数字低通、高通、带通、带阻滤波器。
    总结来讲，就是在数字域进行频带变换。

  • 具体方法
    低通数字滤波器的系统函数：$HL(z)$频带变化后的数字滤波器的系统函数：$Hd(Z)$要求：(1)两个系统函数的单位圆要对应；(2)都是因果稳定；关系：$G(Z^{-1})=\pm\prod{k=1}^N\frac{Z^{-1}-a^*}{1-a_kZ^{-1}}$（全通函数）应为有理函数)

    • 数字低通$→$数字低通
      $N=1,z^{-1}=G(Z^{-1})=\frac{Z^{-1}-\alpha}{1-\alpha Z^{-1}}$,根据两个低通滤波器的对应关系（技术指标）,最终得到$\alpha=\frac{\sin\frac{\theta_c-\omega_c}{2}}{\sin\frac{\theta_c+\omega_c}{2}}$img

    • 数字低通$→$数字高通
      低通频率响应在单位圆上旋转$180^∘$,即得到高通频率响应：Z=-Z。

      $G(Z^{-1})=\frac{-Z^{-1}-\alpha}{1+\alpha Z^{-1}}$

      数字低通滤波器与高通滤波器之间的映射图