SciPy 简易使用教程
- 1. 符号计算
- 2. 函数向量化
- 3. 波形处理scipy.signal
- 3.1 滤波器
- 3.2 波峰定位
基于numpy的一个高级模块,为数学,物理,工程等方面的科学计算提供无可替代的支持。
做重要的思想是:符号计算和函数向量化
1. 符号计算
demo: n次多项式的表示,以及计算
from scipy import poly1dp = poly1d([3, 4, 5]) # p(x) = 3x^2 + 4x + 5
print(p)
print(p * p) # p*p = 9x^4 + 24x^3 + 46x^2 + 40x + 25
print(p.integ(k=6)) # 求p(x)的不定积分,指定常数项为6
print(p.deriv()) # 求p(x)的一阶导数: 6x + 4
print(p([4, 5])) # 求p(4),p(5)的结果
输出
2
3 x + 4 x + 54 3 2
9 x + 24 x + 46 x + 40 x + 253 2
1 x + 2 x + 5 x + 66 x + 4
[ 69 100]
2. 函数向量化
demo: 写一个支持向量操作函数
def addsubtract(a, b):# 结合scipy的函数实现向量操作if a > b:return a - belse:return a + b
vec_addsubtract = np.vectorize(addsubtract)
print(vec_addsubtract([0, 3, 6, 9], [1, 3, 5, 7]))
vec_poly1d = np.vectorize(p)
print(vec_poly1d([4, 5])) # 输出结果和p([4, 5]) 是一致的
3. 波形处理scipy.signal
3.1 滤波器
基本低通、高通、带通、带阻滤波器,关键点是计算Wn。假设采样频率为1000hz,信号本身最大的频率为500hz,要滤除400hz以上频率成分,即截止频率为400hz,则wn=2*400/1000=0.8。Wn=0.8
数字信号Wn=1NT∗21W_n =\frac{\frac{1}{N_T} * 2}{1}Wn=1NT1∗2 其中 N_T 为一个周期内的采样点数
滤波器的阶数: 物理上理解,有几个滤波器串联。
from scipy import signal
b, a = signal.butter(8, 0.8, 'lowpass') #8 为滤波器的阶数
filtedData = signal.filtfilt(b, a, data) #data为要过滤的信号b, a = signal.butter(8, 0.2, 'highpass') #高通
b, a = signal.butter(8, [0.2,0.8], 'bandpass') #带通
b, a = signal.butter(8, [0.2,0.8], 'bandstop') #带阻
疑问:滤波器的数字实现、频率滤波器和均值滤波器的关系。
参考文档:Python实现信号滤波(基于scipy)
3.2 波峰定位
from scipy.signal import find_peaks
x = [......] # 波形信号
peaks, _ = find_peaks(x, height=0) # 高度大于0的峰
peaks, _ = find_peaks(x, height=(-border, border)) # 高度范围
peaks, _ = find_peaks(x, distance=150) # 要求至少距离150个样本的峰
peaks, properties = find_peaks(x, prominence=(None, 0.6)) # 嘈杂信号的识别,突出基线信号0.6 高度的峰
peaks, properties = find_peaks(x, prominence=1, width=20) # 宽度要求# 配套显示工作
import matplotlib.pyplot as plt
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
突出程度:峰的突出度衡量一个峰从信号的周围基线中突出多少,并定义为峰与其最低轮廓线之间的垂直距离,
参考文档:python scipy signal.find_peaks用法及代码示例