用firwin设计滤波器

显然此频率响应和理想的低通滤波器相差甚远，并且即使增加FIR滤波器的系数也没有作用。因为我们舍弃了n<0的那一半系数，而这些系数有着相当大的影响，因此只截取n>=0的部分是不够的，如果我们将n<0的那一半系数也添加进滤波器的话，得到的频率响应将会有很大的改善。如下重新定义h_ideal函数，它返回中-n到n之间的系数：

def h_ideal(n, fc):
    return 2*fc*np.sinc(2*fc*np.arange(-n, n, 1.0))

下面是添加n<0系数之后的频率响应：

对称截取sinc函数的低通滤波器

这样做虽然改善了频率响应，但是给系统带来了许多延时，为了频率响应更好必须增加滤波器的点数，然而为了减少延时，必须减少点数，为了解决这个矛盾，我们给系数乘上一个窗函数，让它快速收敛。

SciPy提供了firwin用窗函数设计低通滤波器，firwin的调用形式如下：

firwin(N, cutoff, width=None, window='hamming')

其中N为滤波器的长度；cutoff为以正规化的频率；window为所使用的窗函数。

下面的程序用firwin设计低通滤波器，并且和上面的结果进行比较，注意由于firwin的cutoff频率是以取样频率/2正规化的，因此它是前面所介绍的的两倍。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

def h_ideal(n, fc):
    return 2*fc*np.sinc(2*fc*np.arange(-n, n, 1.0))

b = h_ideal(30, 0.25) # 以fs正规化的频率
b2 = signal.firwin(len(b), 0.5) # 以fs/2正规化的频率

w, h = signal.freqz(b)
w2, h2 = signal.freqz(b2)

pl.figure(figsize=(8,6))
pl.subplot(211)
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)), label=u"h_ideal")
pl.plot(w2/2/np.pi, 20*np.log10(np.abs(h2)), label=u"firwin")
pl.xlabel(u"正规化频率 周期/取样")
pl.ylabel(u"幅值(dB)")
pl.legend()
pl.subplot(212)
pl.plot(b, label=u"h_ideal")
pl.plot(b2, label=u"firwin")
pl.legend()
pl.show()

firwin使用窗函数设计的低通滤波器的频率响应和脉冲响应

使用firwin函数设计的滤波器并不是最优化的，为了实现同样效果频率响应，还存在长度更短的FIR滤波器。

用remez设计滤波器

remez函数能够帮助我们找到更优的滤波器系数。remez的调用形式如下：

remez(numtaps, bands, desired,
    weight=None, Hz=1, type='bandpass', maxiter=25, grid_density=16)

其中：

• numtaps : 所设计的FIR滤波器的长度
• bands ： 一个递增序列，它包括频率响应中的所有频带的边界，其值在0到Hz/2之间，如果参数Hz为缺省值1的话，那么可以把它当作是以取样频率正规化的频率
• desired : 长度为bands的一半的增益序列，它给出频率响应在bands中的每个频带的增益值
• weight : 长度和desired一样的权重序列，它给出desired中的每个增益所占的权重，即给出desired中的每个增益的重要性，值越大表示其越重要
• type : 'bandpass'或者'differentiator'，本书只介绍type为'bandpass'的情况

remez返回经过remez算法最优化之后的FIR滤波器的系数。此系数和用firwin所设计的结果一样是对称的。当numtaps为偶数时，所设计的滤波器对于取样频率/2的响应为0，因此无法设计出长度为偶数的高通滤波器。

下面的程序演示通过remez设计高通滤波器：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

for length in [11, 31, 51, 101, 201]:
    b = signal.remez(length, (0, 0.18,  0.2,  0.50), (0.01, 1))
    w, h = signal.freqz(b, 1)
    pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)), label=str(length))
pl.legend()
pl.xlabel(u"正规化频率 周期/取样")
pl.ylabel(u"幅值(dB)")
pl.title(u"remez设计高通滤波器 - 滤波器长度和频响的关系")
pl.show()

程序中，remez函数的bands参数给出两个频带(以取样频率正规化)：0到0.18和0.2到0.5，而desired给出两个频带的增益分别为0.01和1，因此它所设计的是一个通带频率为0.2、阻带增益为-40dB的高通滤波器。

此程序显示出滤波器长度和频率响应之间存在如下关系，可以看出滤波器越长，频率响应越接近设计值：

remez设计的高通滤波器，长度越长频率响应越接近设计值

下图显示权值和频率响应之间的关系，图中的滤波器长度为101。我们注意到，当权值为1, 0.01(红色曲线)时，两个频带的增益抖动量相同，这个权值正好和增益desired的设置相反。这时因为缺省情况下，增益越大的频带的频率响应要求越精确，而当权值和增益的乘积相等时，频率响应的误差也就相同了。

remez设计滤波器时权值影响频率响应

滤波器级联

假设有两个滤波器h1和h2，我们将h1的输出输入到h2，这样得到的滤波器称为h1和h2的级联。级联后的滤波器的脉冲响应为h1和h2的脉冲响应的卷积，而其频率响应为两个滤波器的频率响应的乘积。

下面的程序先用remez分别设计一个高通滤波器h1和一个低通滤波器h2，然后通过卷积计算出它们的级联滤波器h3的系数：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

h1 = signal.remez(201, (0, 0.18,  0.2,  0.50), (0.01, 1))
h2 = signal.remez(201, (0, 0.38,  0.4,  0.50), (1, 0.01))
h3 = np.convolve(h1, h2)

w, h = signal.freqz(h3, 1)
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)))

pl.legend()
pl.xlabel(u"正规化频率 周期/取样")
pl.ylabel(u"幅值(dB)")
pl.title(u"低通和高通级联为带通滤波器")
pl.show()

最后使用freqz函数计算h3的频率响应：

低通和高通滤波器级联之后是带通滤波器

可以看出，所得到的是一个带通滤波器。

我们也可以直接用remez设计带通滤波器：

>>> h4 = signal.remez(201, (0, 0.18, 0.2, 0.38, 0.4, 0.50), (0.01, 1, 0.01))

如果你观察此滤波器的频率响应的话，发现它和h3的基本一致，如果比较h3和h4的话，我们得到如下结果：

级联的滤波器和remz设计的带通滤波器的脉冲响应近似

可以看出虽然h3的长度几乎是h4的两倍，但是由于它的许多系数都接近于0，因此h3和h4的频率响应近似相同。

巴特沃斯低通滤波器

巴特沃斯低通滤波器的振幅的平方和频率之间的关系可以用如下公式表示：

其中n为滤波器的阶数，为振幅下降3dB时的截止频率。这个公式很容易理解：

• 当时，越小，振幅越接近于1
• 当时，越大，振幅越接近于0
• 随着n的增大，振幅接近于1或者0的速度将变快，即n越大，低通滤波器在阻频带的衰减速度将越快
• 当时，振幅的平方为1/2，即-3dB

下面我们推导出巴特沃斯低通滤波器的传递函数H(s)，其中，为复数平面上的点。

由于当时，，因此将带入到巴特沃斯低通滤波器的振幅平方公式中可以得到：

此公式有2n个极点，其中n个在左半平面，n个在右半平面，由于H(s)必须是稳定的，因此左半平面的n的极点属于H(s)。

最后得到的传递函数为：

其中为左半平面上的极点：

下面的程序绘制6、7阶巴特沃斯低通滤波器的S复平面上的极点：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# -*- coding: utf-8 -*-
from scipy import signal
import numpy as np
import matplotlib.pyplot as pl

pl.figure(figsize=(5,5))

b, a = signal.butter(6, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), '^', label=u"6阶巴特沃斯极点")

b, a = signal.butter(7, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), 's', label=u"7阶巴特沃斯极点")

pl.axis("equal")
pl.legend(loc="center right")
pl.show()

程序中，使用butter函数设计巴特沃斯滤波器，缺省情况下它设计的是数字滤波器，为了设计模拟滤波器，需要传递关键字参数analog=1。获得传递函数的b和a的系数之后，通过tf2zpk函数将它们转换为零点和极点：

巴特沃斯低通滤波器在S复平面上的极点分布

双线性变换

有了连续时间的传递函数H(s)，下一步就是如何将它转换为离散时间的传递函数H(z)。转换的方法有几种，其中最常用的是双线性变换，其变换公式为：

其中T为离散时间的取样周期。双线性变换公式的推导过程请参考下面的链接：

双线性变换公式推导： http://en.wikipedia.org/wiki/Bilinear_transform

双线性变换实际上是s复平面和z复平面上的点的映射变换，他将s复平面上的竖线变换成z复平面上的圆，而s复平面上的Y轴对应于z复平面上的单位圆。下面的程序演示了这一对应关系：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl

def stoz(s):
    """
    将s复平面映射到z复平面
    为了方便起见，假设取样周期T=1
    """
    return (2+s)/(2-s)
    
def make_vline(x):
    return x + 1j*np.linspace(-100.0,100.0,20000)
    
fig = pl.figure(figsize=(7,3))    
axs = pl.subplot(121)
axz = pl.subplot(122)
for x in np.arange(-3, 4, 1):
    s = make_vline(x)
    z = stoz(s)
    axs.plot(np.real(s), np.imag(s))
    axz.plot(np.real(z), np.imag(z))

axs.set_xlim(-4,4)
axz.axis("equal")
axz.set_ylim(-3,3)

pl.show()

程序中的stoz函数是将s变换为z的变换公式，只需要对上述的双线性变换公式稍作变形即可得到，这里为了方便起见，假设取样周期T=1。

双线性变换将s平面(左图)上的竖线变换为z平面上的圆(右图)

通过双线性变换之后，滤波器的频率响应会发生变化。在下一节中我们会介绍，离散时间的滤波器的频率响应是将其传递函数H(z)用进行替换。将其带入到双线性变换公式得到：

对于s平面上的点来说，过原点的竖线就是连续时间传递函数H(s)的频率响应，因此双线性变换将离散时间的频率通过如下的公式转换为连续时间的频率:

而其反函数为：

下面让我们用程序来验证这个频率转换公式。首先我们载入所需要的库，并且定义离散时间的取样频率fs为8kHz，设计的巴特沃斯低通滤波器通带截至频率为1kHz：

>>> from scipy import signal
>>> from numpy import *
>>> fs = 8000.0
>>> f = 1000.0

下面使用butter函数设计一个3阶的巴特沃斯滤波器，注意关键字参数analog=1，表示设计连续时间传递函数H(s)的系数，由于通带频率参数为圆频率，因此需要乘以2*pi：

>>> b, a = signal.butter(3, 2*pi*f, analog=1)

然后调用双线性变换函数bilinear，将系数转换为离散时间的传递函数系数，通过关键字参数fs指定取样频率：

>>> b2, a2 = signal.bilinear(b,a,fs=fs)

接下来调用freqz函数得到此数字滤波器的频率响应，为了得到尽可能精确的值，我们通过worN关键字参数让它计算10000点的频率响应：

>>> w2, h2 = signal.freqz(b2,a2,worN=1000)

接下来将h2转换为增益，并且找到增益为-3dB(精确值为10*log10(0.5))时所对应的正规化圆频率w的下标idx，w/(2*pi)*fs就是其对应的实际频率值：

>>> p2 = 20*log10(abs(h2))
>>> idx = argmin(abs(p2-10*log10(0.5)))
>>> w2[idx]/2/pi*8000
952.8

通过频率转换公式得到的频率为：

>>> 2*fs*arctan(2*pi*f/2/fs) /2/pi
952.8840223

实际使用scipy.signal库设计IIR滤波器没有这么麻烦，因为它所提供的滤波器设计函数缺省都是直接设计数字滤波器。这些函数设计数字滤波器时采用的取样频率为2，即以香农频率fs/2为1进行正规化。因此要设计取样频率为fs、通带频率为f的滤波器需要将通带频率正规化为f/(fs/2)，下面调用butter函数设计数字低通滤波器，这里使用上述计算所得的通带频率：

>>>  b3,a3 = signal.butter(3, 952.8840223/(fs/2))
>>> sum(abs(b3-b2))
1.3226225670237568e-13
>>> sum(abs(a3-a2))
7.0876637892069994e-13

所得的滤波器的系数b3和a3与手工通过bilinear函数计算的系数b2和a2是一致的。在signal库设计数字滤波器时，其内部会先通过频率转换公式对频率进行转换，然后设计连续时间的传递函数系数，最后通过bilinear函数进行系数转换。有兴趣的读者可以查看signal.iirfilter函数的源代码。

滤波器的频带转换

只要知道了低通滤波器的传递函数H(s)，就很容易利用变量替换设计出同样阶数的高通、带通或者其它通带频率的低通滤波器。让我们来看看低通滤波器的变换。

假设我们使用巴特沃斯低通滤波器设计公式，设计出通带频率为1弧度/秒的标准低通滤波器：

>>> b, a = signal.butter(2, 1.0, analog=1)
>>> np.real(b)
array([ 1.])
>>> np.real(a)
array([ 1.        ,  1.41421356,  1.        ])

为了让它变为通带频率为的低通滤波器，只需要进行如下替换：

由于当时，H(s)就是滤波器的频率响应。因此所设计的标准低通滤波器H(s)在时振幅下降3dB的，而则在处下降3dB。下面是通带频率为2弧度/秒的2阶低通滤波器的系数：

>>> b2, a2 = signal.butter(2, 2.0, analog=1)
>>> np.real(b2)
array([ 4.])
>>> np.real(a2)
array([ 1.        ,  2.82842712,  4.        ])

可以看出将代入到前面的H(s)中即可得到这些系数。

低通滤波器转高通滤波器的替代公式为：

此替代公式很容易理解：

• 若为0，则替代之后的频率为无穷大，而低通滤波器无穷大处的频率响应为0，即转换之后的滤波器在0处的频率响应为0；
• 若为无穷大，则替代之后的频率为0，因此转换之后的滤波器在无穷大处频率响应为1。

下面设计通带频率为1弧度/秒的高通滤波器：

>>> b3,a3 = signal.butter(2,1.0,btype="high",analog=1)
>>> np.real(b3)
array([ 1.,  0.,  0.])
>>> np.real(a3)
array([ 1.        ,  1.41421356,  1.        ])

可以看出这些系数是将带入到H(s)中之后，上下分母乘以之后得到的。

低通滤波器还可以转换为带通滤波器，这可能有点难以理解，让我们先来看看替代公式，假设带通滤波器的高低通带频率为和：

其中，。 被称为通带带宽，而则是通带的中心频率。

让我们通过下面的程序来研究一下为何这种替代能够将低通映射为带通滤波器：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

# -*- coding: utf-8 -*-
import numpy as np
from scipy import signal
import pylab as pl

b, a = signal.butter(2, 1.0, analog=1)

# 低通->带通的频率变换函数
w1 = 1.0 # 低通带频率
w2 = 2.0 # 高通带频率
dw = w2 - w1 # 通带宽度
w0 = np.sqrt(w1*w2) # 通带中心频率

# 产生10**-2到10**2的频率点
w = np.logspace(-2, 2, 1000)

# 使用频率变换公式计算出转换之后的频率
nw = np.imag(w0/dw*(1j*w/w0 + w0/(1j*w)))

_, h = signal.freqs(b, a, worN=nw)
h = 20*np.log10(np.abs(h))

pl.figure(figsize=(8,5))

pl.subplot(221)
pl.semilogx(w, nw) # X轴使用log坐标绘图
pl.xlabel(u"变换前圆频率(弧度/秒)")
pl.ylabel(u"变换后圆频率(弧度/秒)")

pl.subplot(222)
pl.plot(h, nw)
pl.xlabel(u"低通滤波器的频率响应(dB)")

pl.subplot(212)
pl.semilogx(w, h)
pl.xlabel(u"变换前圆频率(弧度/秒)")
pl.ylabel(u"带通滤波器的频率响应(dB)")

pl.subplots_adjust(wspace=0.3, hspace=0.3, top=0.95, bottom=0.14)

print "center:", w[np.argmin(np.abs(nw))]
pl.show()

程序中先使用butter函数设计一个模拟的二阶标准低通滤波器：

b, a = signal.butter(2, 1.0, analog=1)

我们要将其转换为通带频率为w1=1到w2=2的带通滤波器：

# 低通->带通的频率变换函数
w1 = 1.0 # 低通带频率
w2 = 2.0 # 高通带频率
dw = w2 - w1 # 通带带宽
w0 = np.sqrt(w1*w2) # 通带中心频率

假设我们关心的频率响应的频率段为0.01到100，使用logspace函数产生一个这个区间的等比数列w：

w = np.logspace(-2, 2, 1000)

通过带通频率转换公式将其转换为新的频率序列nw：

nw = np.imag(w0/dw*(1j*w/w0 + w0/(1j*w)))

使用此新的频率序列nw计算出每个频率点对应的低通滤波器的频率响应，注意我们通过worN关键字传输让freqs函数计算指定频率的频率响应：

_, h = signal.freqs(b, a, worN=nw)
h = 20*np.log10(np.abs(h))

下面的语句就绘制出w和h的关系，也就是带通滤波器的频率响应：

pl.semilogx(w, h)

下图是程序的输出。其中左上图绘制的是频率转换函数，右上图绘制的是低通滤波器的频率响应(X轴为响应，Y轴为频率)，最下面绘制的是最终的带通滤波器的频率响应。

使用频率转换公式将低通变为带通滤波器

由于带通的频率转换公式将0到无穷大映射到正无穷到负无穷，而低通滤波器在正负无穷处的频率响应都为0，因此可以想象转换后的滤波器是一个带通滤波器。而转换之后的频率nw为0的点所对应的原始频率就是带通滤波器的中心频率，此处的频率响应为1，下面的程序找到nw绝对值最小的下标，并输出其对应的转换前的频率，我们看到它和w0是一致的：

>>> print w[np.argmin(np.abs(nw))]
1.4130259906

事实上，scipy.signal库已经为我们提供了频带转换的函数：

• lp2lp : 低通转低通
• lp2hp : 低通转高通
• lp2bp : 低通转带通
• lp2bs : 低通转带阻，转带阻的公式留给读者思考

下面以lp2bp为例简要说明一下函数的用法，假设b,a为二阶标准低通滤波器，下面的语句将转换为通带为1到2弧度的带通滤波器，前两个参数为滤波器的系数，后两个参数分别为中心频率和通带带宽：

>>> b, a = signal.butter(2, 1.0, analog=1)
>>> b3, a3 = signal.lp2bp(b,a,np.sqrt(2), 1)

我们也可以直接调用butter设计一个低通滤波器：

>>> b4, a4 = signal.butter(2, [1,2], btype='bandpass', analog=1)

两个结果是完全一致的：

>>> np.all(b3==b4)
True
>>> np.all(a3==a4)
True

>>> import scipy.signal
>>> signal.freqz??

N=1, slope=5.9955865774 dB
N=2, slope=12.0417201051 dB
N=3, slope=18.0630802032 dB
N=4, slope=24.0841135443 dB
N=5, slope=30.1051375912 dB