6.scipy Signal

1import pylab as pl
2import numpy as np
3from scipy import signal

1import matplotlib as mpl
2mpl.rcParams['font.sans-serif'] = ['SimHei']

 1#%fig=使用中值滤波剔除瞬间噪声
 2t = np.arange(0, 20, 0.1)
 3x = np.sin(t)
 4x[np.random.randint(0, len(t), 20)] += np.random.standard_normal(20)*0.6 #❶
 5x2 = signal.medfilt(x, 5) #❷
 6x3 = signal.order_filter(x, np.ones(5), 2)
 7print (np.all(x2 == x3))
 8pl.plot(t, x, label=u"带噪声的信号")
 9pl.plot(t, x2 + 0.5, alpha=0.6, label=u"中值滤波之后的信号")
10pl.legend(loc="best");

 1sampling_rate = 8000.0
 2
 3# 设计一个带通滤波器：
 4# 通带为0.2*4000 - 0.5*4000
 5# 阻带为<0.1*4000, >0.6*4000
 6# 通带增益的最大衰减值为2dB
 7# 阻带的最小衰减值为40dB
 8b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40) #❶
 9
10# 使用freq计算滤波器的频率响应
11w, h = signal.freqz(b, a) #❷
12
13# 计算增益
14power = 20*np.log10(np.clip(np.abs(h), 1e-8, 1e100)) #❸
15freq = w / np.pi * sampling_rate / 2

 1#%fig=用频率扫描波测量的频率响应
 2# 产生2秒钟的取样频率为sampling_rate Hz的频率扫描信号
 3# 开始频率为0， 结束频率为sampling_rate/2
 4t = np.arange(0, 2, 1/sampling_rate) #❶
 5sweep = signal.chirp(t, f0=0, t1=2, f1=sampling_rate/2) #❷
 6# 对频率扫描信号进行滤波
 7out = signal.lfilter(b, a, sweep) #❸
 8# 将波形转换为能量
 9out = 20*np.log10(np.abs(out)) #❹
10# 找到所有局部最大值的下标
11index = signal.argrelmax(out, order=3)  #❺
12# 绘制滤波之后的波形的增益
13pl.figure(figsize=(8, 2.5))
14pl.plot(freq, power, label=u"带通IIR滤波器的频率响应") 
15pl.plot(t[index]/2.0*4000, out[index], label=u"频率扫描波测量的频谱", alpha=0.6) #❻
16pl.legend(loc="best")
17#%hide
18pl.title(u"频率扫描波测量的滤波器频谱")
19pl.ylim(-100,20)
20pl.ylabel(u"增益(dB)")
21pl.xlabel(u"频率(Hz)");

 1#%fig=系统的阶跃响应和正弦波响应
 2m, b, k = 1.0, 10, 20
 3
 4numerator = [1]
 5denominator = [m, b, k]
 6
 7plant = signal.lti(numerator, denominator)  #❶
 8
 9t = np.arange(0, 2, 0.01)
10_, x_step = plant.step(T=t)  #❷
11_, x_sin, _ = signal.lsim(plant, U=np.sin(np.pi * t), T=t)  #❸
12#%hide
13pl.plot(t, x_step, label=u"阶跃响应")
14pl.plot(t, x_sin, label=u"正弦波响应")
15pl.legend(loc="best")
16pl.xlabel(u"时间（秒）")
17pl.ylabel(u"位移（米）")