数字信号处理·滤波器之章
你说得对，但是《数字信号处理》是电子信息工程学院独立开设的一门核心专业课，课程发生在一个被称作D221的幻想世界，在这里，被神选中的序列将被授予“傅里叶变换”，导引频域分解之力。玩家将扮演一位名为“DSP”的神秘角色，在自由的系统中邂逅性格各异、能力独特的信号们，和他们一起分析频谱，找回失散的频段——同时，逐步发掘“线性系统”的真相。
数字滤波器设计
数字滤波器设计是通过幅频响应\(|H(e^{j\omega})|\)设计系统函数\(H(z)\)，在数字系统中编程实现。利用模拟滤波器设计数字滤波器，就是要寻求某种映射，把\(s\)平面映射到\(z\)平面，为了因果稳定性，这个映射应该满足以下两个条件：
- \(s\)平面的虚轴\(j\Omega\)映射为\(z\)平面的单位圆
- \(s\)平面的左半平面映射为\(z\)平面的单位圆内
冲激响应不变法设计IIR滤波器
冲激响应不变法的根本思想是从系统的冲激响应出发，保证数字系统的冲激响应是模拟系统的冲激响应的采样，即： \[ h[n]=h_c(t)|_{t=nT_s}=h_c(nT_s) \] 其中\(T_s\)是采样周期。
对模拟系统冲激响应进行理想采样，有： \[ h_s(t)=\sum_{n=-\infty}^{\infty} h_c\left(n T_s\right) \delta\left(t-n T_{\mathrm{s}}\right) \] 其拉普拉斯变换为： \[ H_s(s)=\sum_{n=-\infty}^\infty h[n]e^{-snT_s} \] 数字系统冲激响应的\(z\)变换为： \[ H(z)=\sum_{n=-\infty}^\infty h[n]z^{-n} \] 所以，有： \[ H_s(s)=H(z)|_{z=e^{sT_s}} \] 把模拟信号的拉普拉斯变换在\(s\)平面上沿着虚轴以\(\Omega_s=2\pi/T_s\)为周期进行周期延拓后，按照\(z=e^{sT_s}\)，即\(s=\ln z/T_s\)的规则替换变元，即可得到\(H(z)\)。需要注意，和\(H(z)\)有关的并不是一开始的系统函数\(H(s)\)，而是它周期延拓后的结果\(H_s(s)\)，如果只知道系统函数，那么要么从冲激响应的路径解决问题，要么先进行周期延拓。
这个映射可以保证把\(s\)左半平面映射到\(z\)平面单位圆的内部，保证稳定的模拟滤波器能变换出稳定的数字滤波器。
因为上述过程中存在周期延拓，因此需要保证模拟滤波器是带限的，但是实际上这并不可能做到，因此用冲激响应不变法必定产生频域混叠。
【例】现有系统： \[ H_c(s)=\sum_{k=1}^N \frac{A_k}{s-s_k} \] 用冲激响应不变法设计数字滤波器，求其\(H(z)\)
【解】对系统函数作反变换，得冲激响应： \[ h_c(t)=\sum_{k=1}^N A_ke^{s_kt}u(t) \] 采样，得： \[ h_s[n]=\sum_{k=1}^N T_n A_k\left(\mathrm{e}^{s_k T}\right)^n \mathrm{u}[n] \] 作\(z\)变换，得： \[ H(z)=\sum_{k=1}^N \frac{T_s A_k}{1-\mathrm{e}^{s_k T_s} z^{-1}} \]
观察上例，\(s=s_k\)处的极点变换成Z平面中\(e^{s_kT_s}\)处的极点，所以左半S平面的极点可以映射到Z平面的单位圆内；但离散时间系统函数中的零点是部分分式展开式中的极点和系数的函数，通常并不按照与极点相同的方式进行映射，因此冲激响应不变法不能保证最小相位系统映射后的仍然是最小相位系统。
双线性变换法设计IIR滤波器
为了解决冲激响应不变法的频域混叠现象，需要研究双线性变换法。双线性变换法有两步：
第一步，把\(s\)平面的\(j\Omega\)轴压缩到\(s_1\)平面的\(j\Omega_1\)轴\([-\pi/T_s,\pi/T_s]\)的一段。可以通过正切变换，即： \[ \Omega=c\tan \frac{\Omega_1T_s}{2} \] 把这个关系拓展到整个\(s,s_1\)平面，有： \[ s=\frac{2}{T_s}\frac{1-\mathrm{e}^{-s_1 T_s}}{1+\mathrm{e}^{-s_1 T_s}}=\frac{2}{T_{\mathrm{s}}} \tanh \frac{s_1 T_{s}}{2} \] 再通过\(z=e^{s_1T_s}\)变换到\(z\)平面，有： \[ s=\frac 2{T_s}\frac{1-z^{-1}}{1+z^{-1}} \] 或者说， \[ z=\frac{1+(T_s/2)s}{1-(T_s/2)s} \] 如果\(z=re^{j\omega},s=\sigma+j\Omega\)，有： \[ \begin{gathered} r=\left[\frac{\left(2 / T_{\mathrm{s}}+\sigma\right)^2+\Omega^2}{\left(2 / T_{\mathrm{s}}-\sigma\right)^2+\Omega^2}\right] \\ \omega=\arctan \frac{\Omega}{2 / T_{\mathrm{s}}+\sigma}+\arctan \frac{\Omega}{2 / T_{s}-\sigma} \end{gathered} \] 因为变换的过程中频率产生了畸变，所以设计时可以进行“预畸变”，根据 \[ \Omega=\frac{2}{T_s}\tan\frac{\omega}{2} \] 双线性变换的零极点是按照同种方式映射的，因此可以保证最小相位、全通性质在映射前后不变。而且它完全规避了频域混叠，生成的离散滤波器和连续滤波器有一一对应的关系，这一点也是和冲激响应不变法不同的。
窗函数法设计FIR滤波器
基本原理
FIR滤波器的设计任务是，选择有限长度的单位脉冲响应\(h[n]\)，使得系统的频响\(H(e^{j\omega})\)满足技术要求。
窗函数法的基本思想是，首先对理想零相位低通滤波器，画出其单位脉冲响应\(h[n]\)（应该是一个偶序列），如下图所示：
现在把这个\(h[n]\)向右移动\(M/2\)位，变成\(h[n-M/2]\)，然后加一个窗函数\(w[n]=R_M[n]\)，形成如下图所示的序列\(h_d[n]\)：
由此可见，因果可实现的FIR系统的幅度函数\(|H(e^{j\omega})|\)是理想滤波器的幅度函数和窗函数\(w[n]\)的幅度函数\(W_g(e^{j\omega})\)的卷积，而其相频响应来自于位移所产生的\(e^{-j\frac{\omega M}{2}}\).有以下结论：
- 窗函数旁瓣影响通带内波纹
- 窗函数旁瓣影响阻带内波纹
- 主瓣宽度影响过渡带宽度
- 增加窗的长度可以减小主瓣宽度，从而减小过渡带，但是不能改变肩峰的相对值
常用窗函数
矩形窗 \[ w[n]=R_M[n] \]
主瓣宽度 旁瓣峰值 过渡带宽度 通带阻带波纹最大起伏 \(\frac{4\pi}{M+1}\) -13dB \(\frac{1.8\pi}{M+1}\) -21dB 
Barlett窗 \[ w[n]= \begin{cases}2 n / M, & 0 \leqslant n \leqslant M / 2 \\ 2-2 n / M, & M / 2<n \leqslant M \\ 0, & \text {others}\end{cases} \]
主瓣宽度 旁瓣峰值 过渡带宽度 通带阻带波纹最大起伏 \(\frac{8\pi}{M+1}\) -25dB \(\frac{4.1\pi}{M+1}\) -25dB 
Hanning窗 \[ w[n]= \begin{cases}0.5-0.5 \cos (2 \pi n / M), & 0 \leqslant n \leqslant M \\ 0, & \text {others}\end{cases} \] | 主瓣宽度 | 旁瓣峰值 | 过渡带宽度 | 通带阻带波纹最大起伏 | | :----------------: | :------: | :------------------: | :------------------: | | \(\frac{8\pi}{M+1}\) | -41dB | \(\frac{6.2\pi}{M+1}\) | -44dB |
它能使能量更有效地集中在主瓣，但主瓣更长
Kaiser窗 \[ w[n]= \begin{cases}\frac{I_0\left[\beta\left(1-[(n-\alpha) / \alpha]^2\right)^{\frac{1}{2}}\right]}{I_0(\beta)}, & 0 \leqslant n \leqslant M \\ 0, & \text {others}\end{cases} \] 其中 \[ \beta= \begin{cases}0.1102(A-8.7), & A>50 \\ 0.5842(A-21)^{0.4}+0.07886(A-21), & 21 \leqslant A \leqslant 50 \\ 0, & A<21\end{cases} \]
\[ \begin{align} \alpha&=\frac M2\\ M&=\frac{A-8}{2.285(\omega_s-\omega_p)}\\ A&=-20\lg \min\{\delta_1,\delta_2\} \end{align} \]
\(I_0\)是第一类零阶修正贝塞尔函数。于是，得到的因果滤波器为： \[ h_{\text {lpd }}[n]=\frac{\sin \left[\omega_c(n-M / 2)\right]}{\pi(n-M / 2)} w[n] \] 其中 \[ \omega_c=\frac{\omega_s+\omega_p}{2} \]
频率采样法设计FIR滤波器
这个方法的思想是，根据给定的性能确定系统的频响，然后频响进行采样得到采样值\(H[k]\)，对其进行IDFT得到\(h[n]\)。具体来说，步骤如下：
根据阻带最小衰减选择过渡带采样点数\(m\)
有\(A=-20\lg\delta \mathrm{ dB}\)，根据表：
过渡带采样点数 1 2 3 阻带衰减/dB 44-54 54-75 85-95 选择\(m\)
确定过渡带宽度\(\Delta \omega=\omega_s-\omega_p\)，根据 \[ \Delta \omega \geq \frac{2\pi(m+1)}{M+1} \] 确定滤波器长度\(M+1\)
根据待设计滤波器的性能指标，确定是第几类广义线性相位系统，有以下规则：
- \(M\)是偶数那么在I、III中选，否则在II、IV中选
- 第I类可做各种滤波器，第II类可做低通带通、第III类可做带通、第IV类可做高通带阻
确定理想滤波器的幅频\(A(e^{j\omega})\)，相频\(e^{j(\beta-\alpha \omega)}\)
对频响进行采样，对于不同的广义线性相位系统，有以下公式成立：
第I类，\(M\)是偶数，\(\beta=0\ or\ \pi\) \[ H[k]= \begin{cases}A[k] \mathrm{e}^{-j \frac{M\pi k}{M+1}}, & k \in\left[0, \frac{M}{2}\right] \\ H^*[M+1-k], & k \in\left[\frac{M}{2}+1, M\right]\end{cases} \]
第II类，M是奇数，\(\beta = 0\ or\ \pi\) \[ H[k]= \begin{cases}A[k] \mathrm{e}^{-j \frac{M\pi k}{M+1}}, & k \in\left[0, \frac{M-1}{2}\right] \\ 0, & k=\frac{M+1}{2} \\ H^*[M+1-k], & k \in\left[\frac{M+3}{2}, M\right]\end{cases} \]
第III类，M是偶数。\(\beta = \pi/2 \ or \ 3\pi/2\) \[ H[k]= \begin{cases}0, & k=0 \\ {jA}[k] \mathrm{e}^{-j\frac{M\pi k}{M+1}},& k \in\left[1, \frac{M}{2}\right] \\ H^*[M+1-k], & k \in\left[\frac{M}{2}+1, M\right]\end{cases} \]
第IV类，M是奇数，\(\beta = \pi/2 \ or \ 3\pi/2\) \[ H[k]= \begin{cases}0, & k=0 \\ {jA}[k] \mathrm{e}^{-j\frac{M\pi k}{M+1}}, & k \in\left[1, \frac{M+1}{2}\right] \\ H^*[M+1-k], & k \in\left[\frac{M+3}{2}, M\right]\end{cases} \]
求\(H[k]\)的\(M+1\)点IDFT，得到\(h[n]\)