SciPy的signal模块提供了信号处理方面的许多函数，包括卷积运算、B样条、滤波以及滤波器设计等方面的内容。

一、中值滤波

中值滤波能够比较有效地消除声音信号中的瞬间噪声或者图像中的斑点噪声。在signal模块中，medfilt( )对一维信号进行中值滤波，而medfilt2d( )对二维信号进行中值滤波。在scipy.ndimage模块中另有针对多维图像的中值滤波器，这里简单演示medfilt( )的效果。

t = np.arange(0, 20, 0.1)
x = np.sin(t)
x[np.random.randint(0, len(t), 20)] += np.random.standard_normal(20)*0.6 ?
x2 = signal.medfilt(x, 5) ?
x3 = signal.order_filter(x, np.ones(5), 2)
print np.all(x2 == x3)
pl.plot(t, x, label=u"带噪声的信号")
pl.plot(t, x2 + 0.5, alpha=0.6, label=u"中值滤波之后的信号")
pl.legend(loc="best")
True

?首先创建一个带有随机的瞬间噪声的正弦波，?然后调用medfilt( )进行中值滤波，第二个参数为计算中值的窗口大小，它必须是一个奇数。medfilt( )将信号中的每个元素都替换为其窗口内的中值。

最后绘制原始信号和滤波信号，为了便于比较，图中将滤波之后的信号统一向上偏移了0.5，结果如图1所示。中值滤波是排序滤波的一个特例。使用排序滤波可以将元素替换为其窗口内指定排序顺序的元素。其调用形式如下：

order_filter(a, domain, rank)

图1　使用中值滤波剔除瞬间噪声

其中a是一个多维数组，domain是维数和a相同的数组，它指定窗口的范围，rank是一个非负整数，用来选择窗口中元素排序后的值，0表示选择最小值，1表示选择第二小的值。中值滤波也可以用order_filter( )计算，注意domain参数是一个长度为5、值全为1的数组。

二、滤波器设计

signal模块提供了许多滤波器设计的函数。在下面的实例中，我们设计一个IIR带通滤波器，并查看其频率响应，最后使用它对频率扫描信号进行滤波计算。

sampling_rate = 8000.0

# 设计一个带通滤波器：
# 通带为0.2*4000 - 0.5*4000
# 阻带为<0.1*4000, >0.6*4000
# 通带增益的最大衰减值为2dB
# 阻带的最小衰减值为40dB
b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40) ?

# 使用freq计算滤波器的频率响应
w, h = signal.freqz(b, a) ?

# 计算增益
power = 20*np.log10(np.clip(np.abs(h), 1e-8, 1e100)) ?
freq = w / np.pi * sampling_rate / 2

?首先用iirdesign( )设计一个IIR带通滤波器。这个滤波器的通带为0.2f0到0.5f0，阻带为小于0.1f0和大于0.6f0，其中f0为信号取样频率的一半。如果取样频率为8kHz，那么这个带通滤波器的通带为800Hz到2kHz。通带的最大增益衰减为2dB，阻带的最小增益衰减为40dB，即通带的增益浮动在2dB之内，阻带至少有40dB的衰减。

iirdesgin( )返回两个数组b和a，它们分别是IIR滤波器的分子和分母部分的系数。其中a[0]恒等于1。?调用freqz( )计算所得到的滤波器的频率响应。freqz( )返回两个数组w和h，其中w是圆频率数组，通过ωf0/π可以计算出与其对应的实际频率。h是w中对应频率点的响应，它是一个复数数组，其幅值表示滤波器的增益特性，相角表示滤波器的相位特性。

?计算h的增益特性，并使用dB进行度量。由于h中存在幅值几乎为0的值，因此先用clip( )对其裁剪之后，再调用对数函数，避免计算出错。

在实际运用中为了测量未知系统的频率特性，经常将频率扫描波输入到系统中，观察系统的输出，从而计算其频率特性。下面让我们模拟这一过程：

# 产生两秒钟的取样频率为sampling_rate Hz的频率扫描信号
# 开始频率为0，结束频率为sampling_rate/2
t = np.arange(0, 2, 1/sampling_rate) ?
sweep = signal.chirp(t, f0=0, t1=2, f1=sampling_rate/2) ?
# 对频率扫描信号进行滤波
out = signal.lfilter(b, a, sweep) ?
# 将波形转换为能量
out = 20*np.log10(np.abs(out)) ?
# 找到所有局部最大值的下标
index = signal.argrelmax(out, order=3)  ?
# 绘制滤波之后的波形的增益
pl.figure(figsize=(8, 2.5))
pl.plot(freq, power, label=u"带通IIR滤波器的频率响应")
pl.plot(t[index]/2.0*4000, out[index], label=u"频率扫描波测量的频谱", alpha=0.6) ?
pl.legend(loc="best")

?为了调用chirp( )产生频率扫描波形的数据，首先需要产生一个表示取样时间的等差数组，这里产生两秒的取样频率为8kHz的取样时间数组。?然后调用chirp( )得到两秒的频率扫描波形的数据。频率扫描波的开始频率f0为0Hz，结束频率f1为4kHz，到达4kHz的时间为两秒，使用数组t作为取样时间点。?最后调用lfilter( )计算频率扫描波形经过带通滤波器之后的结果。

?为了和系统的增益特性图进行比较，需要获取输出波形的包络，因此先将输出波形数据转换为能量值。?为了计算包络，调用argrelmax( )找到out数组中所有局部最大值的下标，order参数指定局域最大值的范围，这里的值为3表示所有的局域最大值都是连续7个元素（前后各三个元素）中的最大值。?最后将时间转换为对应的频率，绘制所有局部最大点的能量值。

图2显示了freqz( )计算的频谱和频率扫描波得到的频率特性，可以看到结果是一致的。

图2　用频率扫描波测量的频率响应

三、连续时间线性系统

在上一节中，我们使用odeint( )对质量-弹簧-阻尼系统的微分方程组进行了数值积分，并且进行了PID控制模拟。该系统的微分方程为：。通过拉普拉斯变换可以将微分方程化为容易求解的代数方程：ms2X(s)+bsX(s)+kX(s)=F(s)。其中F(s)是f(t)的拉普拉斯变换，X(s)是x(t)的拉普拉数变换，而n次微分变成了sn。F(s)是输入信号，而X(s)是输出信号，将等式改写为输入除以输出的形式，就得到了系统的传递函数P(s)：

连续时间系统的传递函数是两个s的多项式的商。通过连续时间系统的传递函数，很容易计算某输入信号对应的输出信号。在下面的例子中，使用signal模块计算质量-弹簧-阻尼系统对阶跃信号以及正弦波信号的响应输出。?创建lti对象，可以使用控制理论中的多种形式表示连续时间线性系统，这里使用的是传递函数分子和分母多项式的系数。多项式的系数与numpy.poly1d的约定相同，即下标为0的元素是最高次项的系数。?调用lti.step( )方法计算系统的阶跃响应。T参数为计算响应的时间数组。?调用signal.lsim( )计算系统对正弦波信号的响应，它的第一个参数为lti对象，也可以直接传递(numerator, denominator)。U参数为保存输入信号的数组。step( )和lsim( )计算结果中的第二项为系统的输出信号，这里忽略其余的输出。

图3显示阶跃响应最终稳定在x=0.05处，这时的kx=1。

图3　系统的阶跃响应和正弦波响应

m, b, k = 1.0, 10, 20

numerator = [1]
denominator = [m, b, k]

plant = signal.lti(numerator, denominator) ?

t = np.arange(0, 2, 0.01)
_, x_step = plant.step(T=t) ?
_, x_sin, _ = signal.lsim(plant, U=np.sin(np.pi*t), T=t) ?

传递函数的代数运算可以表示由多个连续时间系统组成的系统，例如两个系统的级联的传递函数是各个系统的传递函数的乘积。而传递函数由分子和分母两个多项式构成，因此传递函数的四则运算可以使用NumPy的poly1d相关的函数实现。下面的SYS类通过定义__mul__、__add__、__sadd__、__div__等魔法方法，让它支持四则运算。

图4　反馈控制系统框图

?feedback( )方法计算与之对应的反馈系统的传递函数。在图4中，P是被控制的系统，C是控制器，C的输入信号是目标信号与实际输入的差xr-x。我们从xr到x的传递函数就是这个反馈系统的传递函数。根据图示可以列出如下拉普拉斯变换之后的代数方程：

整理可得：

如果将C(s)·P(s)看作系统Y(s)，那么可以得出反馈系统的传递函数为：

?为了让SYS对象能作为step( )、lsim( )等函数的第一个表示系统的参数，需要定义__iter__( )魔法方法返回传递函数的分子与分母的多项式系数。

from numbers import Real

def as_sys(s):
if isinstance(s, Real):
return SYS([s], [1])
return s

class SYS(object):
def __init__(self, num, den):
self.num = num
self.den = den

def feedback(self): ?
return self / (self + 1)

def __mul__(self, s):
s = as_sys(s)
num = np.polymul(self.num, s.num)
den = np.polymul(self.den, s.den)
return SYS(num, den)

def __add__(self, s):
s = as_sys(s)
den = np.polymul(self.den, s.den)
num = np.polyadd(np.polymul(self.num, s.den),
np.polymul(s.num, self.den))
return SYS(num, den)

def __sadd__(self, s):
return self + s

def __div__(self, s):
s = as_sys(s)
return self * SYS(s.den, s.num)

def __iter__(self): ?
return iter((self.num, self.den))

下面我们用SYS类计算使用PI控制器控制质量-弹簧-阻尼系统时的阶跃响应。PI控制器的传递函数为：

注意上节中介绍的PI控制器是离散时间的，使用累加器近似计算积分器的输出，而本节采用连续时间系统的系统响应模拟控制系统。

?质量-弹簧-阻尼系统的传递函数为plant，?PI控制器的传递函数为pi_ctrl，为了step( )不抛出LinAlgError异常，这里将PI控制器的传递函数的分母常数项设置为一个非常小的值。?计算反馈系统的传递函数feedback。由图5可以看出Ki为0时，系统的输出位移与目标位移之间存在一定的差距，Kp越大差距越小，但是会出现过冲现象。适当调节Kp与Ki可以减弱过冲现象，但是仍然会有超过目标位移的时刻。

M, b, k = 1.0, 10, 20
plant = SYS([1], [M, b, k]) ?

pi_settings = [(10, 1e-10), (200, 1e-10),
(200, 100),  (50, 100)]

fig, ax = pl.subplots(figsize=(8, 3))

for pi_setting in pi_settings:
pi_ctrl = SYS(pi_setting, [1, 1e-6]) ?
feedback = (pi_ctrl * plant).feedback( ) ?
_, x = signal.step(feedback, T=t)    
label = "$K_p={:d}, K_i={:3.0f}#34;.format(*pi_setting)
ax.plot(t, x, label=label)

ax.legend(loc="best", ncol=2)
ax.set_xlabel(u"时间(秒)")
ax.set_ylabel(u"位移(米)")

图5　使用PI控制器的控制系统的阶跃响应

为了计算施加于质量的控制力，可以将误差信号传递给lsim( )计算控制器的输出：

_, f, _ = signal.lsim(pi_ctrl, U=1-x, T=t)

为了彻底消除过冲现象，需要使用PID控制，PID控制器的传递函数为：

下面计算PID控制器构成的反馈系统的阶跃响应。由于PID控制器需要对输入信号进行微分，而阶跃输入信号会导致PID的输出中包含脉冲输出，即时间无限短、值无限大的信号。

kd, kp, ki = 30, 200, 400
pid_ctrl = SYS([kd, kp, ki], [1, 1e-6])
feedback = (pid_ctrl * plant).feedback( )
_, x2 = signal.step(feedback, T=t)

为了让PID控制器的输出在限定的范围之内，可以在反馈系统之前添加一个低通滤波器，一阶低通滤波器的传递函数为：

添加低通滤波器之后，PID控制器的输入就是连续信号了，如图6所示。

图6　带低通滤波器的反馈控制系统框图

lp = SYS([1], [0.2, 1])
lp_feedback = lp * (pid_ctrl * plant).feedback( )
_, x3 = signal.step(lp_feedback, T=t)

由于PID控制器的传递函数的分子阶数高于分母阶数，因此无法使用lsim( )计算。我们可以把x_r当作系统的输入，把f当作输出，通过下面的方程计算从x_r到f的传递函数：

F(s)=(Xr(s)·LP(s)-F(s)·P(s))·C(s)

得到的传递函数为：

下面根据上面的公式计算带低通滤波器的控制系统中控制器的输出：

pid_out = (pid_ctrl * lp) / (pid_ctrl * plant + 1)
_, f3 = signal.step(pid_out, T=t)

图7显示了上述PI控制、PID控制以及带低通滤波的PID控制等系统中滑块的位移以及控制力。由于PID控制的控制力存在脉冲信号，因此无法在图中正确显示。由位移曲线可以看出低通+PID控制可以有效抑制过冲现象。

图7　滑块的位移以及控制力