记有滤波器\(H_a(j\omega)\), 衡量滤波器的性能参数一般由通带边缘频率\(\omega_p\), 阻带起始频率\(\omega_s\), 通带波纹\(\delta_p\)(\(\alpha_p\)), 阻带波纹\(\delta_s\)(\(\alpha_s\)), 选择参数\(k\), 辨别参数\(k_1\)来衡量:
\[\begin{aligned}
& 1 - \delta_p \le |H_a(j\omega)| \le 1 + \delta_p, |\omega| \le \omega_p \\
& |H_a(j\omega)| \le \delta_s, |\omega_s| \le |\omega| \le \infty \\
& \alpha_p = -20 \log_{10}(1-\delta_p) \mathrm{dB} \\
& \alpha_s = -20 \log_{10}(\delta_s) \mathrm{dB} \\
& k = \frac{\omega_p}{\omega_s} \\
& k_1 = \frac{\epsilon}{\sqrt{A^2 - 1}}, A = \frac{1}{\delta_s}, \epsilon = \sqrt{\frac{2\delta_p}{1-\delta_p}}
\end{aligned}
\]
Butterworth低通滤波器的转移函数和幅频响应函数如下:
\[\begin{aligned}
& H_a(s) = \frac{z(s)}{p(s)} = \frac{\omega_c^N}{\prod_{l=1}^{N}(s-p_l)}, p_l=\omega_c e^{j[\pi(N+2l-1)/2N]} \\
& |H_a(j\omega)|^2 = \frac{1}{1+(\omega/\omega_c)^{2N}}
\end{aligned}
\]
Butterworth低通滤波器的阶数选择:
\[\begin{aligned}
& |H_a(j\omega_p)|^2 = \frac{1}{1+(\omega_p/\omega_c)^{2N}} = \frac{1}{\epsilon^2 + 1} \\
& |H_a(j\omega_s)|^2 = \frac{1}{1+(\omega_s/\omega_c)^{2N}} = \frac{1}{A^2} \\
& N = \frac{1}{2} \frac{\log_{10}[(A^2-1)/\epsilon^2]}{\log_{10}(\omega_s/\omega_p)} = \frac{\log_{10}(1/k_1)}{\log_{10}(1/k)}
\end{aligned}
\]
Chebyshev I型低通滤波器的幅频响应和转移函数如下:
\[\begin{aligned}
& |H_a(j\omega)|^2 = \frac{1}{1+\epsilon^2T_N^2(\omega/\omega_p)} \\
& T_N(\omega)=
\begin{cases}
\cos(N\cos^{-1}(\omega)), |\omega| \le 1 \\
\cosh(N\cosh^{-1}{\omega}), |\omega| \gt 1
\end{cases} \\
& \epsilon = \sqrt{\frac{2\delta_p}{1-\delta_p}}
\end{aligned}
\]
\[\begin{aligned}
& H_a(s) = \frac{z(s)}{p(s)} = \frac{\omega_c^N}{\prod_{l=1}^{N}(s-p_l)}, p_l=\delta_l+j\omega_l \\
& \delta_l = -\omega_p\xi\sin[\frac{(2l--1)\pi}{2N}], \omega_l=\omega_p\zeta\cos[\frac{(2l-1)\pi}{2N}] \\
& \xi=\frac{\gamma^2-1}{2\gamma}, \zeta=\frac{\gamma^2+1}{2\gamma}, \gamma=(\frac{1+\sqrt{1+\epsilon^2}}{\epsilon})^{1/N}
\end{aligned}
\]
Chebyshev I型低通滤波器的阶数选择:
\[N = \frac{\cosh^{-1}(1/k_1)}{\cosh^{-1}(1/k)}
\]
Chebyshev II型低通滤波器的幅频响应和转移函数如下:
\[\begin{aligned}
& |H_a(j\omega)|^2 = \frac{1}{1+\epsilon^2[\frac{T_N(\omega_s/\omega_p)}{T_N(\omega_s/\omega)}]^2}
\end{aligned}
\]
\[\begin{aligned}
& H_a(s) = \frac{z(s)}{p(s)} = \frac{z_l}{\prod_{l=1}^{N}(s-p_l)}, p_l=\delta_l+j\omega_l \\
\end{aligned}
\]
Chebyshev II型低通滤波器的阶数选择:
\[N = \frac{\cosh^{-1}(1/k_1)}{\cosh^{-1}(1/k)}
\]
椭圆近似低通滤波器的幅频响应和阶数选择如下:
\[\begin{aligned}
& |H_a(j\omega)|^2 = \frac{1}{1+\epsilon^2R_N^2(\omega/\omega_p)} \\
& N \cong \frac{2\log_{10}(4/k_1)}{\log_{10}(1/\rho)} \\
& \begin{cases}
k' = \sqrt{1-k^2} \\
\rho_0 = \frac{1-\sqrt{k'}}{2(1+\sqrt(k'))} \\
\rho = \rho_0 + 2(\rho_0)^5 + 15(\rho_0)^9 + 150(\rho_0)^{13}
\end{cases}
\end{aligned}
\]
Bessel低通滤波器转移函数如下:
\[\begin{aligned}
& H_a(s) = \frac{d_0}{B_N(s)}=\frac{d_0}{s^N+d_{N-1}s^{N-1}+\cdots+d_1s+d0} \\
& B_N(s) = (2N-1)B_{N-1}(s) + s^2B_{N-2}(s) \\
& B_1(s) = s+1 \\
& B_2(s) = s^2+3s+3 \\
& d_l = \frac{(2N-l)!}{2^{N-l}l!(N-l)!}
\end{aligned}
\]
记有低通滤波器\(H_{LP}(s)\).
从低通滤波器\(s\)映射到高通滤波器\(\hat{s}\), 有如下关系:
\[s=\frac{\omega_p \hat{\omega_p}}{\hat{s}}
\]
从低通滤波器\(s\)映射到带通滤波器\(\hat{s}\), \(\omega_o\)为带通中心点, \(\omega_{p2},\omega_{p1}\)为带通的边缘频率, 有如下关系:
\[s=\omega_p\frac{\hat{s^2}+\hat{\omega_o^2}}{\hat{s}(\hat{\omega_{p2} - \hat{\omega_{p1}}})}
\]
从低通滤波器\(s\)映射到带阻滤波器\(\hat{s}\), \(\omega_o\)为阻带中心点, \(\omega_{s2},\omega_{s1}\)为阻带的结束和起始频率, 有如下关系:
\[s = \omega_s\frac{\hat{s}(\hat{\omega_{s2}} - \hat{\omega_{s1}})}{\hat{\omega_s^2} + \hat{\omega_o^2}}
\]
见信号与线性系统.md, 有:
\[s=\frac{1}{T}\ln{z} \approx \frac{2}{T}(\frac{1-z^{-1}}{1+z^{-1}})
\]
故从s域的虚轴映射到z域的单位圆上(收敛的临界条件)有:
\[j\Omega = \frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{j\omega}} = j\frac{2}{T}\tan(\frac{\omega}{2}) \Rightarrow \Omega=\frac{2}{T}\tan(\frac{\omega}{2})
\]
因此, IIR数字滤波器设计步骤:
- 确定滤波器参数: 通带边缘频率, 阻带起始频率, 通带波纹, 阻带波纹, 选择参数和辨别参数;
- 根据滤波器参数选择模拟低通滤波器类型, 并确定滤波器的阶数, 得到滤波器的转移函数;
- 再将所得到的低通滤波器转义函数重新映射到所需要的高通/带通等转移函数;
- 将公式\(s=\frac{2}{T}(\frac{1-z^{-1}}{1+z^{-1}})\)带入到模拟滤波器的转移函数, 于是得到z域的转移函数;
有时会遇到这样的需求: 已经设计好了低通滤波器的结构, 但是需要改变该低通滤波器的截止频率. 此时, 一般没有必要重新设计滤波器结构, 可以再当前滤波器上再乘以一个全通滤波器即可完成截止频率的转换. 记需从截止频率为\(\omega_c\)的低通滤波器变换到\(\omega_c'\)的低通滤波器:
\[\begin{aligned}
& \because |H_{ap}(e^{j\omega})|^2 = 1 \\
& \therefore H_{ap}(e^{j\omega}) = \pm\frac{d_N+d_{N-1}e^{-j\omega}+\cdots+d_1e^{-j\omega(N-1)}+e^{-j\omega N}}{1+d_1e^{-j\omega}+\cdots+d_{N-1}e^{-j\omega(N-1)} +d_Ne^{-j\omega N}} \Rightarrow \\
& \quad H_{ap}(z) = \pm \frac{d_N+d_{N-1}z^{-1} + \cdots+d_1z^{-N+1}+z^{-N}}{1+d_1z^{-N}+\cdots+d_{N-1}z^{1-N}+d_Nz^{-N}} = \pm \prod_{i=1}^{N}\frac{-\lambda_i^* + z^{-1}}{1-\lambda_iz^{-1}} \\
& \quad H(z')=H(H_{ap}(z)) \Rightarrow H_{ap}(z) = \frac{1-\lambda z'}{z' - \lambda} \Rightarrow \\
& \quad e^{-j\omega}=\frac{e^{-j\omega'} - \lambda}{1-\lambda e^{-j\omega'}} \Rightarrow \tan(\omega_c/2)=\frac{1+\lambda}{1-\lambda}\tan(\omega_c'/2) \Rightarrow \\
& \quad \lambda = \frac{\sin(\frac{\omega_c-\omega_c'}{2})}{\sin(\frac{\omega_c-\omega_c'}{2})}
\end{aligned}
\]
FIR一般用来设计线性相位数字滤波器, 所谓线性相位需满足恒群延时\(\frac{\rm{d}\theta(\omega)}{\rm{d}\omega}=-\tau\), 或者恒群延时和恒相延时\(\frac{\theta(\omega)}{\omega}=-\tau\)都满足. 其中, \(\tau\)为常数.
当恒相延时和恒群延时同时成立时:
\[\begin{aligned}
& H(e^{j\omega}) = \sum_{n=0}^{N}h(n)\cdot(\cos(\omega n) - j\sin(\omega n)) \Rightarrow \\
& \tan(\theta(\omega))=\tan(-\tau \omega) \Rightarrow \frac{\sin(\tau \omega)}{\cos(\tau \omega)}=\frac{\sum_{n=0}^{N}h(n)\sin(\omega n)}{\sum_{n=0}^{N}h(n)\cos(\omega n)} \Rightarrow \\
& \sum_{n=0}^{N}h(n)sin(\omega (\tau - n)) = 0 \Rightarrow \\
& \begin{cases}
\tau = \frac{N}{2} \\
h(n) = h(N-n), 0 \le n \le N
\end{cases}
\end{aligned}
\]
当仅满足恒群延时时:
\[\begin{aligned}
& H(e^{j\omega}) = \sum_{n=0}^{N}h(n)\cdot(\cos(\omega n) - j\sin(\omega n)) \Rightarrow \\
& \tan(\theta(\omega))=\tan(\theta_0 -\tau \omega) \Rightarrow \frac{\sin(-\theta_0 + \tau \omega)}{\cos(-\theta_0 + \tau \omega)}=\frac{\sum_{n=0}^{N}h(n)\sin(\omega n)}{\sum_{n=0}^{N}h(n)\cos(\omega n)} \Rightarrow \\
& \sum_{n=0}^{N}h(n)sin(-\theta_0 + \omega (\tau - n)) = 0 \Rightarrow \\
& \begin{cases}
\tau = \frac{N}{2} \\
h(n) = -h(N-n), 0 \le n \le N
\end{cases}
\end{aligned}
\]
当转移函数满足线性相位时, 有:
\[\begin{aligned}
& H(z) = \pm\sum_{n=0}^{N}h(N-n)z^{-n} = \pm z^{1-N}H(z^{-1}) \\
& \exists\quad 共轭镜面对称的零点z_i, z_i^*, 1/z_i, 1/z_i^*
\end{aligned}
\]
综上, 根据采样点奇偶数不同, 及线性相位满足条件不同, 有4中线性相位FIR滤波器.
\[\begin{aligned}
& H(z) = \sum_{n=0}^{(N-1)/2}h(n)\cdot(z^{-n}+z^{-(N-n)}) \Rightarrow \\
& \begin{cases}
|H(\omega)| = 2\sum_{n=1}^{\frac{N+1}{2}}h(\frac{N+1}{2}-n)\cos((n-\frac{1}{2})\omega) \\
\theta(\omega) = -\omega (\frac{N}{2})
\end{cases}
\end{aligned}
\]
\(\cos((n-\frac{1}{2})\omega)\)在\(\omega=\pi\)时为0, 因此不适合用来设计高通/带阻滤波器(在\(\omega=\pi\)处不为0), 可用来设计低通/带通FIR数字滤波器.
\[\begin{aligned}
& H(z) = h(\frac{N}{2}) + \sum_{n=0}^{\frac{N}{2}-1}h(n)\cdot(z^{-n}+z^{-(N-n)}) \Rightarrow \\
& \begin{cases}
|H(\omega)| = h(\frac{N}{2}) + 2\sum_{n=1}^{\frac{N}{2}}h(\frac{N}{2}-n)\cos(n\omega) \\
\theta(\omega) = -\omega (\frac{N}{2})
\end{cases}
\end{aligned}
\]
\(\cos(n\omega)\)关于\(0,\pi,2\pi\)偶对称, 故可用来设计低通/高通/带通/带阻FIR数字滤波器.
\[\begin{aligned}
& H(z) = \sum_{n=0}^{(N-1)/2}h(n)\cdot(z^{-n}-z^{-(N-n)}) \\
& \begin{cases}
|H(\omega)| = 2\sum_{n=1}^{\frac{N+1}{2}}h(\frac{N+1}{2}-n)\sin((n-\frac{1}{2})\omega) \\
\theta(\omega) =\frac{\pi}{2} -\omega (\frac{N}{2})
\end{cases}
\end{aligned}
\]
\(\sin((n-\frac{1}{2})\omega)\)在\(\omega=0\)时为0, 故不可用来设计低通滤波器和带阻滤波器, 可用来设计高通/带通FIR数字滤波器.
\[\begin{aligned}
& H(z) = \sum_{n=0}^{\frac{N}{2}-1}h(n)\cdot(z^{-n}-z^{-(N-n)}) \\
& \begin{cases}
|H(\omega)| = 2\sum_{n=1}^{\frac{N}{2}}h(\frac{N}{2}-n)\sin(n\omega) \\
\theta(\omega) =\frac{\pi}{2} -\omega (\frac{N}{2})
\end{cases}
\end{aligned}
\]
\(\sin(n\omega)\)在\(\omega=0,\pi\)处为0, 故不可用来设计带阻/高通/低通滤波器, 可用来设计带通FIR数字滤波器.
理想的低通滤波器, 其时域冲击响应序列为\(h(n)=\frac{\sin(\omega_c n)}{\pi n}, -\infty\le n \le\infty\), 该序列并不绝对收敛, 但是能量主要集中在主瓣上. 因此实际应用中常常加窗函数截断主瓣或加上前几个旁瓣\(h(n)*W(n)\)作为系统的冲击响应序列.
假设所需要的抽样点数为\(N\), 则窗函数的窗口大小至少为\(2N+1\), 常用窗函数如下:
\[\begin{aligned}
& Rectangular(n) =\begin{cases}
1, -N \le n \le N \\
0, \quad others
\end{cases} \\
& Bartlett(n) = \begin{cases}
1 - \frac{|n|}{N+1}, -N \le n \le N \\
0, \quad others
\end{cases} \\
& Hann(n) = \begin{cases}
\frac{1}{2}[1-\cos(\frac{2\pi n}{2N+1})], -N \le n \le N \\
0, \quad others
\end{cases} \\
& Hanmming(n) = \begin{cases}
0.54 - 0.46\cos(\frac{2\pi n}{2N+1}), -N \le n \le N \\
0, \quad others
\end{cases} \\
& Blackman(n) = \begin{cases}
0.42 - 0.5\cos(\frac{2\pi n}{2N+1}) + 0.08\cos(\frac{4\pi n}{2N+1})
\end{cases}
\end{aligned}
\]
以上几种窗函数的属性如下:
窗类型 |
主瓣宽度 |
旁瓣峰值衰减(dB) |
最小阻带衰减(dB) |
过渡带宽度(\(\omega_p-\omega_s\)) |
Rectangular |
\(4\pi/(2N+1)\) |
13.3 |
20.9 |
0.92\(\pi/N\) |
Barlett |
\(4\pi/(2N+1)\) |
26.5 |
25 |
\(2.1\pi/N\) |
Hann |
\(8\pi/(2N+1)\) |
31.5 |
43.9 |
\(3.11\pi/N\) |
Hanmming |
\(8\pi/(2N+1)\) |
42.7 |
54.5 |
\(3.32\pi/N\) |
Blackman |
\(12\pi/(2N+1)\) |
58.1 |
75.3 |
\(5.56\pi/N\) |
综上, 设计策略: 通过过渡带宽度求出滤波器阶数\(N\), 根据滤波器类型/相位要求选择FIR类型, 求\(h(n)\cdot W(n)\)得到序列\(h_w(n)\), 将\(h_w(n)\)带入到所选择的FIR类型的转移函数中得到幅频响应.
参考资料
- Digtal Signal Processing, Sanjit K.Mitra.