# 数字信号系统

## FIR和IIR滤波器

FIR滤波器根据如下公式进行计算：

IIR滤波器根据如下公式(直接1型)进行计算：

``````signal.lfilter(b, a, x, axis=-1, zi=None)
``````

``````y[m] = b[0]*x[m] + z[0,m-1]                         (1)
z[0,m] = b[1]*x[m] + z[1,m-1] - a[1]*y[m]           (2)
...
z[n-3,m] = b[n-2]*x[m] + z[n-2,m-1] - a[n-2]*y[m]
z[n-2,m] = b[n-1]*x[m] - a[n-1]*y[m]
``````

``````z[0,m-1] = b[1]*x[m-1] + z[1,m-2] - a[1]*y[m-1]     (3)
``````

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

# 某个均衡滤波器的参数
a = np.array([1.0, -1.947463016918843, 0.9555873701383931])
b = np.array([0.9833716591860479, -1.947463016918843, 0.9722157109523452])

# 44.1kHz， 1秒的频率扫描波
t = np.arange(0, 0.5, 1/44100.0)
x= signal.chirp(t, f0=10, t1 = 0.5, f1=1000.0)

# 直接一次计算滤波器的输出
y = signal.lfilter(b, a, x)

# 将输入信号分为50个数据一组
x2 = x.reshape((-1,50))

# 滤波器的初始状态为0， 长度是滤波器系数长度-1
z = np.zeros(max(len(a),len(b))-1, dtype=np.float)
y2 = [] # 保存输出的列表

for tx in x2:
# 对每段信号进行滤波，并更新滤波器的状态z
ty, z = signal.lfilter(b, a, tx, zi=z)
# 将输出添加到输出列表中
y2.append(ty)

# 将输出y2转换为一维数组
y2 = np.array(y2)
y2 = y2.reshape((-1,))

# 输出y和y2之间的误差
print np.sum((y-y2)**2)

# 绘图
pl.plot(t, y, t, y2)
pl.show()
``````

``````>>> impulse = np.zeros(1000, dtype=np.float)
>>> impulse[0] = 1
>>> h = signal.lfilter(b, a, impulse)
>>> h[-1]
-4.2666825205952273e-12
``````

``````>>> y3 = signal.lfilter(h, 1, x)
>>> np.sum((y-y3)**2)
3.7835244127856444e-17
``````

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

# 某个均衡滤波器的参数
a = np.array([1.0, -1.947463016918843, 0.9555873701383931])
b = np.array([0.9833716591860479, -1.947463016918843, 0.9722157109523452])

# 44.1kHz， 1秒的频率扫描波
t = np.arange(0, 0.5, 1/44100.0)
x= signal.chirp(t, f0=10, t1 = 0.5, f1=1000.0)
y = signal.lfilter(b, a, x)
ns = range(10, 1100, 100)
err = []

for n in ns:
# 计算脉冲响应
impulse = np.zeros(n, dtype=np.float)
impulse[0] = 1
h = signal.lfilter(b, a, impulse)

# 直接FIR滤波器的输出
y2 = signal.lfilter(h, 1, x)

# 输出y和y2之间的误差
err.append(np.sum((y-y2)**2))

# 绘图
pl.figure(figsize=(8,4))
pl.semilogy(ns , err, "-o")
pl.xlabel(u"脉冲响应长度")
pl.ylabel(u"FIR模拟IIR的误差")
pl.show()
``````

## FIR滤波器设计

``````# -*- 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)

w, h = signal.freqz(b)

pl.figure(figsize=(8,3))
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)))
pl.xlabel(u"正规化频率 周期/取样")
pl.ylabel(u"幅值(dB)")
pl.show()
``````

freqz用于计算数字滤波器的频率响应，它的调用方式如下：

``````freqz(b, a=1, worN=None, whole=0, plot=None)
``````

### 用firwin设计滤波器

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

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

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

``````# -*- 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使用窗函数设计的低通滤波器的频率响应和脉冲响应

### 用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是一种迭代算法，它能够找到一个n阶多项式，使得在指定的区间中此多项式和指定函数之间的最大误差最小化。由于FIR滤波器的频率响应实际上是一个多项式函数（请参考下节内容），因此可以用remez算法进行FIR滤波器系数设计。

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

``````# -*- 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设计的高通滤波器，长度越长频率响应越接近设计值

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

### 滤波器级联

``````# -*- 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()
``````

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

## IIR滤波器设计

### 巴特沃斯低通滤波器

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

``````# -*- 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()
``````

### 双线性变换

``````# -*- 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()
``````

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

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

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

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

``````>>> 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
``````

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

### 滤波器的频带转换

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

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

• 为0，则替代之后的频率为无穷大，而低通滤波器无穷大处的频率响应为0，即转换之后的滤波器在0处的频率响应为0；
• 为无穷大，则替代之后的频率为0，因此转换之后的滤波器在无穷大处频率响应为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\.        ])
``````

``````# -*- 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)")

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

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

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

``````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.semilogx(w, h)
``````

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

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

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

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

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

## 滤波器的频率响应

``````def freqz(b, a=1, worN=None, whole=0, plot=None):
b, a = map(atleast_1d, (b,a))
if whole:
lastpoint = 2*pi
else:
lastpoint = pi
if worN is None:
N = 512
w = numpy.arange(0,lastpoint,lastpoint/N)
elif isinstance(worN, types.IntType):
N = worN
w = numpy.arange(0,lastpoint,lastpoint/N)
else:
w = worN
w = atleast_1d(w)
zm1 = exp(-1j*w)
h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
if not plot is None:
plot(w, h)
return w, h
``````

``````w = numpy.arange(0,pi,pi/N)
zm1 = exp(-1j*w)
h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
``````

``````w = numpy.arange(0,pi,pi/N)
``````

``````zm1 = exp(-1j*w)
``````

``````h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
``````

polyval(p, x)函数对于数组x中的每个元素计算多项式p的值，其计算公式如下：

``````p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
``````

freqz中在Z平面单位圆上所取的点是等距线性的，然而我们经常需要在绘制频率响应图表时要求频率坐标为对数坐标，对于对数坐标，等距的频率点会造成低频过疏，高频过密的问题，因此我们可以如下改造freqz函数，使其更适合计算对数频率坐标的频率响应。

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

def logfreqz(b, a, f0, f1, fs, N):
"""
以对数频率坐标计算滤波器b,a的频率响应
f0, f1: 计算频率响应的开始频率和结束频率
fs: 取样频率
"""
w0, w1 = np.log10(f0/fs*2*np.pi), np.log10(f1/fs*2*np.pi)
# 不包括结束频率
w = np.logspace(w0, w1, N, endpoint=False)
zm1 = np.exp(-1j*w)
h = np.polyval(b[::-1], zm1) / np.polyval(a[::-1], zm1)
return w/2/np.pi*fs, h

for n in range(1, 6):
# 设计n阶的通频为0.1*4000 = 400Hz的高通滤波器
b, a = signal.iirfilter(n, [0.1, 1])
f, h = logfreqz(b, a, 10.0, 4000.0, 8000.0, 400)
gain = 20*np.log10(np.abs(h))
pl.semilogx(f, gain, label="N=%s" % n)
slope = (gain[100]-gain[10]) / (np.log2(f[100]) - np.log2(f[10]))
print "N=%s, slope=%s dB" % (n, slope)
pl.ylim(-100, 20)
pl.xlabel(u"频率(Hz)")
pl.ylabel(u"增益(dB)")
pl.legend()
pl.show()
``````

iirfilter设计5个不同阶的IIR高通滤波器

``````# -*- coding: utf-8 -*-

import scipy.signal as signal
import pylab as pl
import math
import numpy as np

def design_equalizer(freq, Q, gain, Fs):
'''设计二次均衡滤波器的系数'''
A = 10**(gain/40.0)
w0 = 2*math.pi*freq/Fs
alpha = math.sin(w0) / 2 / Q

b0 = 1 + alpha * A
b1 = -2*math.cos(w0)
b2 = 1 - alpha * A
a0 = 1 + alpha / A
a1 = -2*math.cos(w0)
a2 = 1 - alpha / A
return [b0/a0,b1/a0,b2/a0], [1.0, a1/a0, a2/a0]

pl.figure(figsize=(8,4))
for freq in [1000, 2000, 4000]:
for q in [0.5, 1.0]:
for p in [5, -5, -10]:
b,a = design_equalizer(freq, q, p, 44100)
w, h = signal.freqz(b, a)
pl.semilogx(w/np.pi*44100, 20*np.log10(np.abs(h)))
pl.xlim(100, 44100)
pl.xlabel(u"频率(Hz)")
pl.ylabel(u"振幅(dB)")
pl.show()
``````

ScrubberEditor的BUG

``````# Establish the slider increment:
increment = self.factory.increment
if increment <= 0.0:
if (low is None) or (high is None) or isinstance( low, int ):
increment = 1.0
else:
increment = pow( 10, round( log10( (high - low) / 100.0 ) ) )

self.increment = increment # ** 将此行移出if作用域
``````