使用sciPy计算信号的波峰和波谷

使用sciPy计算信号的波峰和波谷

使用scipy.signal.find_peaks计算信号的波峰

使用语法：

scipy.signal.find_peaks(x, height=None, threshold=None, distance=None, prominence=None, width=None, wlen=None, rel_height=0.5, plateau_size=None)

参数介绍：

• x(sequence):需要计算峰值的信号序列
• height(number or ndarray or sequence, optional):设置峰值的最小高度，如果是2维元素序列,第一个值解释为最小高度，第二个值为最大高度
• threshold(number or ndarray or sequence, optional):峰值的阈值，与相邻峰值的垂直距离，如果是2维元素序列,第一个值解释为最小值，第二个值为最大值
• distance(number, optional):在峰值间的最小水平距离（>=1）,较小的峰值首先移除,直到满足所有剩余峰值的条件
• prominence(number or ndarray or sequence, optional):达到高峰的高度，如果是2维元素序列,第一个值解释为最小高度，第二个值为最大高度
• width(number or ndarray or sequence, optional):样本中所需波峰的宽度，如果是2维元素序列,第一个值解释为最小高度，第二个值为最大高度
• wlen(int, optional):用于计算峰值突起,因此,只有在参数突出或宽度的情况下才使用
• rel_height(float, optional):用于计算峰宽,因此只在给定宽度时使用。
• plateau_size:number or ndarray or sequence, optional:在样本中平顶波峰所需要的尺寸大小

代码实例：

from scipy.signal import find_peaks
from scipy.misc import electrocardiogram
import matplotlib.pyplot as plt
%matplotlib inline
value = electrocardiogram()[1000:5000]
peak, _ = find_peaks(value, height=0)
plt.plot(value)
plt.plot(peak, value[peak], "*") # peak为横坐标，value[peak]为对应纵坐标
plt.plot(np.zeros_like(value), "--", color="green")
plt.show()

效果图：

使用scipy.signal.peak_prominence()计算信号的波峰高度

主要作用是计算每个波峰的峰值高度。
使用语法：
scipy.signal.peak_prominences(x, peaks, wlen=None)
参数介绍：

• x(sequence):信号数据
• peaks(sequence):峰值索引
• wlen(int, optional):窗口长度，可以可选择地限定每个波峰的评估区域到信号子集
  Returns
prominences(ndarray):计算每个波峰的峰值高度

代码实例：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
%matplotlib inline
array_data = np.linspace(0, 5 * np.pi, 1100)
sig = np.sin(array_data) + 0.5 * np.sin(1.7 * array_data)
sig_peaks, _ = signal.find_peaks(sig)
promin = signal.peak_prominences(sig, sig_peaks)[0]
promin

结果输出：array([1.3780961 , 0.04026852, 2.11918643, 0.3302705 ])

contur_height = sig[sig_peaks] - promin
plt.plot(sig)
plt.plot(sig_peaks, sig[sig_peaks], "*")
plt.vlines(x=sig_peaks, ymin=contur_height, ymax=sig[sig_peaks])
plt.show()

效果图：

使用scipy.signal.find_peaks_cwt()计算1-D信号序列的波峰

使用语法：
scipy.signal.find_peaks_cwt(vector, widths, wavelet=None, max_distances=None, gap_thresh=None, min_length=None, min_snr=1, noise_perc=10, window_size=None)
参数介绍：

• vector(ndarray): 1-D数组
• widths(float or sequence):用于计算CWT矩阵的宽度，通常，这个范围应该覆盖预期的峰值宽度
• wavelet(callable, optional):应该取两个参数并返回一个1-D数组来卷积向量,第一个参数决定了返回的小波阵列的点数,第二个参数是小波的尺度(宽度)。应该归一化和对称。默认是ricker小波。
• max_distances(ndarray, optional):在每一行，只有当在row[n]的相对最大值在max_distances[n]，则从row[n+1]的相对最大值连接脊线。默认值是widths/4
• gap_thresh(float, optional):如果在max_distances里没有找到相对最大值，就会出现间隙，如果超过gap_thresh点而没有连接一个新的相对最大值，则脊线中断。默认是宽度数组的第一个值，也就是widths[0]
• min_length(int, optional):脊线的最小长度必须是可接受的,默认是cwt.shape[0]/4,也就是widths数量的1/4
• min_snr(float, optional):最小信噪比，默认1
• noise_perc(float, optional):在计算噪音下界时,低于认为噪声的要检查的数值点的百分比。使用stats.scoreatpercentile计算。默认是10。
• window_size(int, optional):用于计算噪声下界的窗口大小，默认 cwt.shape[1] / 20
Returns:peaks_indices(ndarray)发现峰值的位置索引。列表被排序。

代码实例：

from scipy.signal import find_peaks_cwt
import numpy as np
x_data = np.arange(0, np.pi, 0.06)
sin_data = np.sin(x_data)
peak_indices = signal.find_peaks_cwt(sin_data, np.arange(2,11))
peak_indices, x_data[peak_indices], sin_data[peak_indices]

结果输出：(array([27], dtype=int64), array([1.62]), array([0.99878974]))

使用scipy.signal.peak_widths计算信号序列的波峰宽度

使用语法：
scipy.signal.peak_widths(x, peaks, rel_height=0.5, prominence_data=None, wlen=None)
参数介绍：

• x(sequence):信号序列或数组
• peaks(sequence):信号序列的波峰索引
• rel_height(float, optional):选择峰值宽度的相对高度,作为其峰高的百分比。1.0计算峰值的宽度在其最低等高线上,0.5值在显著高度的一半。大于等于0。
• prominence_data(tuple, optional):当调用相同参数x和峰值时,三种数组的元组与peak_height的输出相匹配。如果没有提供,这些数据是在内部计算。
• wlen(int, optional):样本中窗口长度传给peak_prominences作为prominence_data内部计算的可选参数。如果prominence_data参数给出，则可忽略
  Returns:
• widths(ndarray):样本中每个峰值的宽度
• width_heights(ndarray):在等高线宽度被评估的高度
• left_ips, right_ips(ndarray):在各自评价高度的水平线的左和右交点的插值位置。

代码实例：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
%matplotlib inline
array_data = np.linspace(0, 5 * np.pi, 1100)
sig = np.sin(array_data) + 0.5 * np.sin(1.7 * array_data)
sig_peaks, _ = signal.find_peaks(sig)
promin = signal.peak_prominences(sig, sig_peaks)[0]
promin
half_peak_res = signal.peak_widths(sig,sig_peaks,rel_height=0.5)
full_peak_res = signal.peak_widths(sig,sig_peaks,rel_height=1)
plt.plot(sig)
plt.plot(sig_peaks, sig[sig_peaks], "*")
plt.hlines(*half_peak_res[1:], color="C4")
plt.hlines(*full_peak_res[1:], color="C5")
plt.show()

效果图：

使用scipy.signal.argrelextrema() 计算信号序列的波峰和波谷

主要作用计算数据的极值
使用语法：
scipy.signal.argrelextrema(data, comparator, axis=0, order=1, mode='clip')
参数介绍：

• data(ndarray): 需要找到极值的数组
• comparator(callable):使用函数来比较两个数据点。应该以两个数组作为参数。
• axis(int, optional):从数据中选择的轴。默认值为0
• order(int, optional):每边有多少个点用来考虑comparator(n,n + x)是真的。
• mode(str, optional):如何处理向量的边缘。wrap(wrap around)或clip(将溢出与最后(或第一个)元素相同)。默认是clip。
  Returns:
  extrema(tuple of ndarrays):在整数数组中最大值的索引。极值[k]是数据轴k的数组。注意,返回值是一个元组,即使数据是1-D
  代码实例：

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

data_x = np.arange(start = 0, stop = 40, step = 1, dtype='int')
data_y = np.array([98,96,97,100,95,105,75,50,45,42,
                    51,85,90,92,91,89,101,62,65,52,
                    47,58,55,75,89,92,94,91,89,79,
                    85,65,42,55,48,50,85,88,95,100])

# Find peaks
# order:两侧使用多少点进行比较
peak_indexes = signal.argrelextrema(data_y, np.greater, order=1)
peak_indexes = peak_indexes[0]

# Find valleys
# order:两侧使用多少点进行比较
valley_indexes = signal.argrelextrema(data_y, np.less, order=1)
valley_indexes = valley_indexes[0]

(fig, ax) = plt.subplots()

# Plot all data
ax.plot(data_x, data_y)

# Plot peaks
peak_x = peak_indexes
peak_y = data_y[peak_indexes]
ax.scatter(peak_x, peak_y, marker='o', color='red', label="Peaks_Pointers")

# Plot valleys
valley_x = valley_indexes
valley_y = data_y[valley_indexes]
ax.scatter(valley_x, valley_y, marker='o', color='green', label="Valleys_Pointers")
ax.plot(peak_x, peak_y, marker='*', linestyle='dashed', color='orange', label="peaks_connetted_line")
ax.plot(valley_x, valley_y, marker='*', linestyle='dashed', color='blue', label="valleys_connetted_line")
# 添加标题
plt.title('Find peaks and valleys using argrelextrema()')
# 添加图例
plt.legend(loc='best')
# 保存图像
# plt.savefig('peaks-valleys.png')
# 显示图像
plt.show()

效果显示：

当order变大时，则图像上的点变得很稀疏，（order=3时）

使用scipy.signal.argrelextrema() 计算信号序列的波谷(极小值)

代码实例：

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('Agg')
from scipy.signal import argrelextrema
array_data = np.arange(start = 0, stop = 40, step = 1, dtype='int')
y_array_data = np.random.random(40)*5
minima_ind = argrelextrema(y_array_data, np.less)
minima_ind = minima_ind[0]
minima_ind
(fig, ax) = plt.subplots()
ax.plot(array_data, y_array_data)

x_minima = minima_ind
y_minima = y_array_data[minima_ind]
ax.plot(x_minima, y_minima, marker='*', linestyle='dashed', color='green', label="Minima")

效果显示：