Automatic multiscale-based peak detection (AMPD)
本文代码的原始出处：推荐一个非常实用的峰值查找算法（peak detection） - 知乎 (zhihu.com)
import numpy as np
import matplotlib.pyplot as plt
# 自动峰检测算法，参考：https://zhuanlan.zhihu.com/p/549588865
def AMPD(data):
"""
实现AMPD算法
在AMPD算法中，将输入数据平均分成多个长度相等的窗口，并逐个对窗口进行峰值检测。首先选择一个窗口长度k（从1开始逐渐增大）。
然后算法会从第k个数据点开始遍历，对于每个数据点i，都会检查其左侧和右侧距离为k的数据点，如果i点的数值比这两个点都要大，那么就认为i点是一个峰值。如果在一个窗口内有多个峰值，则窗口内的峰值指数之和会越小，因为要同时满足窗口内所有峰值的条件。
这样就可以通过比较不同窗口长度的峰值指数之和，找到一个最合适的窗口大小，并用它来检测所有的波峰。
:param data: 1-D numpy.ndarray
:return: 波峰所在索引值的列表
"""
p_data = np.zeros_like(data, dtype=np.int32) # 创建一个和data具有相同形状的新数组，初始化全为0，数据类型为int32。这个数组用于在算法中记录每个数据点的峰值指数，即这个数据点处于多少个窗口的峰值位置。
count = data.shape[0] # 获取数据数组的长度。
arr_rowsum = [] # 用于存储每个窗口的峰值指数之和。因为算法会遍历不同长度的窗口，并计算每个窗口内的峰值的数量，最终累加到对应的arr_rowsum列表位置上，通过比较其中所有元素的大小，选择最小那个值对应的窗口长度作为最终判断波峰的阈值。
for k in range(1, count // 2 + 1): # 指定窗口长度的循环范围，从1到数据长度的1/2，使用循环变量k依次代表窗口长度。
row_sum = 0 # 当前窗口的峰值指数之和
for i in range(k, count - k): # 遍历当前窗口内的所有数据点，判断当前是否为峰值，如果当前点为峰值，则row_sum减去1，表示当前窗口内多了一个峰值。
if data[i] > data[i - k] and data[i] > data[i + k]:
row_sum -= 1
arr_rowsum.append(row_sum) # 将当前窗口的峰值指数之和添加到列表arr_rowsum中，以便后续比较窗口质量。
min_index = np.argmin(arr_rowsum) # 使用NumPy库的argmin函数找到arr_rowsum中最小值对应的索引，也就是最小的峰值指数之和所对应的窗口长度，赋值给max_window_length
max_window_length = min_index
for k in range(1, max_window_length + 1):
for i in range(k, count - k): # 再次判断当前点是否为峰值，并将p_data数组中当前点的值加1，表示当前点出现在一个峰值之中，标记了这个峰值的位置。
if data[i] > data[i - k] and data[i] > data[i + k]:
p_data[i] += 1
return np.where(p_data == max_window_length)[0] # 返回波峰所在索引值的列表
# 注意：第一次对点进行判断，是为了计算当前窗口的峰值个数，用来衡量当前窗口的质量，也就是计算峰值指数之和。第二次对点进行判断，是为了标记每个峰值的位置。
"""
用于模拟生成由3个信号叠加而成，其中包含了频率为100HZ和300HZ的信号以及随机噪声的长度为1000的数组y，来模拟真实环境下信号的特性
"""
def sim_data():
N = 1000
x = np.linspace(0, 200, N)
y = 2 * np.cos(2 * np.pi * 300 * x) \
+ 5 * np.sin(2 * np.pi * 100 * x) \
+ np.random.randn(N)
return y
def wav_data():
# 数据来源：https://github.com/LXP-Never/blog_data/tree/master/machine_learning_date
import scipy.io.wavfile as wf
sample_rate, noised_sigs = wf.read('../实验数据/music.wav')
print(sample_rate) # sample_rate：采样率44100
print(noised_sigs.shape) # noised_sigs:存储音频中每个采样点的采样位移(220500,)
times = np.arange(noised_sigs.size) / sample_rate # 该数组中每个元素表示对应采样点的时刻（单位为秒）。
return noised_sigs[:3000] # 只返回前3000个采样点的采样位移值，即只返回音频信号的前3000个采样点数据。
def vis():
y = sim_data()
# y = wav_data()
plt.plot(range(len(y)), y)
px = AMPD(y) # 获取信号中的峰值位置px。
plt.scatter(px, y[px], color="red") # 在原始信号图像上标注出峰值位置，峰值处的点具有红色。
a = plt.axes([.4, .2, .4, .4]) # 创建一个新的子图，显示信号的局部细节，并将峰值可视化。
a.scatter(px, y[px], color="red")
a.plot(range(len(y)), y, label=u"原信号", color='green')
a.set_xlim(100, 300) # 设置x轴的显示范围
plt.show()
vis()可以对该代码进一步优化，特别是在窗口长度的循环部分。一种可能的优化方法是使用numpy中的矩阵运算来替代循环操作，以提高计算效率。具体来说，可以使用np.maximum函数来找到每个窗口内的最大值，并将其与当前点的值进行比较，得到一个布尔型矩阵，其中True表示当前点处于峰值位置。然后使用np.sum函数对这个矩阵进行求和，即可得到当前窗口内的峰值数量：
def find_peaks(data):
p_data = np.zeros_like(data, dtype=np.int32)
count = data.shape[0]
arr_rowsum = []
for k in range(1, count // 2 + 1): # 使用numpy矩阵运算替代循环操作
max_data = np.maximum(data[k:count-k], data[k-1:count-k-1])
is_peak = (data[k:count-k] > max_data) & (data[k:count-k] > data[k-1:count-k-1])
row_sum = np.sum(is_peak)
arr_rowsum.append(row_sum)
min_index = np.argmin(arr_rowsum)
max_window_length = min_index
for k in range(1, max_window_length + 1): # 同样使用numpy矩阵运算替代循环操作
max_data = np.maximum(data[k:count-k], data[k-1:count-k-1])
is_peak = (data[k:count-k] > max_data) & (data[k:count-k] > data[k-1:count-k-1])
p_data[k:count-k][is_peak] += 1
return np.where(p_data == max_window_length)[0]