PythonLibs(7) - SciPy

基于numpy的一个高级模块，为数学，物理，工程等方面的科学计算提供无可替代的支持。
做重要的思想是：符号计算和函数向量化

1. 符号计算

demo: n次多项式的表示，以及计算

from scipy import poly1d

p = 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 + 5
   4      3      2
9 x + 24 x + 46 x + 40 x + 25
   3     2
1 x + 2 x + 5 x + 6
 
6 x + 4
[ 69 100]

2. 函数向量化

demo: 写一个支持向量操作函数

def addsubtract(a, b):
    # 结合scipy的函数实现向量操作
    if a > b:
        return a - b
    else:
        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
数字信号$W_n=\frac{\frac{1}{N_T}\ast 2}{1}$ 其中 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用法及代码示例