几种常见窗函数的特性

解决一下上一篇的遗留问题

上一篇谈了FIR滤波器的加窗是什么,但是自觉得还是有些地方没有说明白,就好比模拟加窗过程都没有讲明白。

不过今天再看书时,又有所得,忽然明白了为什么模拟加窗过程出了问题,这是由于我之前对于数字滤波器的设计过程还很模糊。

窗函数设计法的设计思路:

  • 给定要求的理想频响\(H_d(e^{j\omega})\),一般为分段常数
  • 转为时域设计,所以需要求出\(h_d(n)\)

\[h_d(n) = IDTFT[H_d(e^{j\omega})]\]

  • 由于\(h_d(n)\)是无限时长的,所以需要加窗截断为\(h(n)\),窗的长度为N

因为窗函数是在时域内截断,所以需要将理想滤波器转换到时域来处理。

\[h(n) = h_d(n)\cdotp w(n), \quad 0\le n\le N-1\]

  • 求出加窗后的实际频响\(H(e^{j\omega})\)

\[H(e^{j\omega}) = DTFT[h(n)]\]

  • 检查\(H(e^{j\omega})\)是否满足\(H_d(e^{j\omega})\),不满足就重复3、4步骤

我使用的子程序是书上提供的,可以产生一个理想滤波器。

今天突然对这个程序有了新的理解。因为这个是子程序产生的时域的滤波器,为了能够在电脑中存储,本来就已经加了一个矩形窗,因此它的频域波形存在波纹;而我又加一个矩形窗,对这个已经加了矩形窗的滤波器当然是一点效果都没有咯。

几种常见窗函数

回到今天的主题,今天探讨一下几种常见滤波器的特性和使用场景。

翻了很多遍书,发现对于滤波器的设计,主要关心的是过渡带宽(Transition bandwidth)通带边沿衰减(Passband ripple)阻带最小衰减(Minimum stopband attenuation),而且大部分的参数都是用dB作单位。

使用dB做单位的好处有:

  • 数值变小。由于分贝是取对数值,所以能很方便的表示大的数量的变化
  • 运算方便。放大器级联时,总的放大倍数是各级相乘。用分贝做单位时,总增益就是相加。
  • 方便感知。人对强度(光照、声音)的感知,接近于于强度的对数的正比。

至于为什么要这么多种窗呢?那是因为不同的窗特性不一样,比如最简单的矩形窗,虽然完成了截断工作,但是通带衰减大、阻带衰减小,导致能量的浪费;而之后的多种窗则或多或少的弥补了这些缺点。

矩形窗

矩形窗的定义为

\[ w(n) = \begin{cases} 1, & 0\le n \le M-1 \\\ 0, & otherwise \end{cases} \]

频率响应函数为

\[ W(e^{j\omega}) = W(\omega)e^{j\theta(\omega)} = \frac{\sin(\omega M /2)}{\sin(\omega /2)}e^{-j\frac{N-1}{2}\omega} \]

因此

\[ W(\omega) = \frac{\sin(\omega M /2)}{\sin(\omega /2)} \]

下面分析窗函数的主要参数:

  • 幅度响应\(W_r(\omega)\)第一个零点在\(\omega=\omega_1\)

\[ \frac{\omega_1 M}{2} = \pi ~~ or ~~ \omega_1 = \frac{2\pi}{M} \]

因此,主瓣宽度为为\(2\omega_1 = 4\pi /M\),因此传输带宽近似于\(4\pi /M\)

  • 第一个旁瓣大概在\(\omega=3\pi/M\)的位置,因此它的幅值为

\[ \left| W_r\left(\omega = \frac{3\pi}{M}\right)\right| = \left|\frac{\sin(3\pi /2)}{\sin(3\pi /2M)}\right| \simeq \frac{2M}{3\pi} \]

对比主瓣的幅值,旁瓣幅值峰值为

\[ \frac{2}{3\pi} \approx 21.22\% \equiv 13 dB \]

三角形窗

由于吉布斯现象,矩形窗存在一个0到1的越变;而三角形窗则提供了一个比较缓慢的变化,它的定义式为:

\[ w(n) = \begin{cases} \frac{2n}{M-1}, & 0\le n \le\frac{M-1}{2} \\\ 2-\frac{2n}{M-1}, & \frac{M-1}{2}\le n\le M-1 \\\ 0, & otherwise \end{cases} \]

谱密度函数表达式如下,’≈’仅当\(M\gg 1\)时成立

\[ W(\omega) \approx \frac{2}{M} \frac{\sin(M\omega/4)}{\sin(\omega/4)} \]

主瓣宽度为\(8\pi /M\),旁瓣峰值衰减为25dB

汉宁窗(Hanning)

这是一个升余弦窗,被定义为

\[ w(n) = \begin{cases} 0.5\left[1-\cos\left(\frac{2\pi n}{M-1}\right)\right] & 0\le n \le M-1 \\\ 0 & otherwise \end{cases} \]

主瓣宽度为\(8\pi /M\),旁瓣峰值衰减为31dB

海明窗(Hamming)

海明窗和汉宁窗很像,不同的是它有一部分是不连续的,被定义为

\[ w(n)= \begin{cases} 0.54 - 0.46\cos\left(\frac{2\pi n}{M-1}\right) & 0\le n \le M-1 \\\ 0 & otherwise \end{cases} \]

主瓣宽度为\(8\pi /M\),旁瓣峰值衰减为41d

布莱克曼窗(Blackman)

这个窗函数和前两个窗函数很像,不过增加了升余弦的二次谐波分量,被定义为

\[ w(n) = \begin{cases} 0.42 - 0.5\cos\left(\frac{2\pi n}{M-1}\right) + 0.08\cos\left(\frac{4\pi n}{M-1}\right) & 0\le n \le M-1 \\\ 0 & otherwise \end{cases} \]

主瓣宽度为\(12\pi /M\),旁瓣峰值衰减为57dB

凯泽窗(Kaiser)

这是一个非常有用的窗函数,它可以同时调整主瓣宽度与旁瓣宽度,这是其他窗函数不具备的,被定义为

\[ w(n) = \frac{I_0\left[\beta \sqrt{1-\left(1-\frac{2n}{M-1}\right)^2}\right]}{I_0[\beta]} \]

\(I_0\)是第一类零阶贝塞尔函数,\(\beta\)是用来调整窗函数性能的参数

本人使用的\(\beta = 8.5\)

如何选择窗函数

选择窗函数可以参考前文中的窗函数设计法

先确定自己的需求,然后根据窗函数的极限性能,做出选择,最后再验证这个窗函数是否符合需求


本文中的代码已上传本人的github


参考书籍:

  • _DigitalSignalProcessingUsingMatlabv4.0_JohnG.Proakis
  • 《数字信号处理教程》程佩青