在 R 中查找具有最小峰宽的峰 - 类似于 MATLAB 函数

Posted 2023-03-25

【中文标题】在 R 中查找具有最小峰宽的峰 - 类似于 MATLAB 函数

【英文标题】：Finding peaks with minimum peak width in R - similar to MATLAB function

【发布时间】：2021-08-14 08:08:48

【问题描述】：

我需要在时间序列数据中找到峰值，但结果需要等于 MATLAB 中 findpeaks 函数的结果，参数“MinPeakWidth”设置为 10。我已经尝试了很多函数为了实现这一点：pracma::findpeaks、fluoR::find_peaks、splus2R::peaks、IDPmisc::peaks（这个有一个关于峰宽的论点，但结果不一样）。我已经查看了其他函数以及，包括用于生物导体中的色谱和光谱分析的软件包。除此之外，我还尝试了 *** 中另一个问题的功能（以及很少的改动）：Finding local maxima and minima

MATLAB 中的 findpeaks 函数用于查找局部最大值，具有以下特点：

找到局部最大值。峰值按出现的顺序输出。尽管是最大值，但不包括第一个样本。对于平坦的峰值，该函数仅返回索引最低的点。

MATLAB 网站中“MinPeakWidth”参数的解释是

最小峰宽，指定为逗号分隔的对，由 'MinPeakWidth' 和一个正实标量组成。使用此参数仅选择宽度至少为“MinPeakWidth”的峰。 如果您指定位置向量 x，则 'MinPeakWidth' 必须以 x 表示。如果 x 是日期时间数组，则将 'MinPeakWidth' 指定为持续时间标量或以天数表示的数值标量。 如果您指定采样率 Fs，则“MinPeakWidth”必须以时间单位表示。 如果您既不指定 x 也不指定 Fs，则“MinPeakWidth”必须以样本为单位表示。

数据类型：双精度 |单身 |持续时间

这是数据：

valores <- tibble::tibble(V1 = c(
0.04386573, 0.06169861, 0.03743560, 0.04512523, 0.04517977, 0.02927114, 0.04224937, 0.06596527, 2.15621006, 0.02547804, 0.03134409, 0.02867694,
0.08251871, 0.03252856, 0.06901365, 0.03201109, 0.04214851, 0.04679828, 0.04076178, 0.03922274, 1.65163662, 0.03630282, 0.04146608, 0.02618668,
0.04845364, 0.03202031, 0.03699149, 0.02811389, 0.03354410, 0.02975296, 0.03378896, 0.04440788, 0.46503730, 0.06128226, 0.01934736, 0.02055138,
0.04233819, 0.03398005, 0.02528630, 0.03694652, 0.02888223, 0.03463824, 0.04380172, 0.03297124, 0.04850558, 0.04579087, 1.48031231, 0.03735059,
0.04192204, 0.05789367, 0.03819694, 0.03344671, 0.05867103, 0.02590745, 0.05405133, 0.04941912, 0.63658824, 0.03134409, 0.04151859, 0.03502503,
0.02182294, 0.15397702, 0.02455722, 0.02775277, 0.04596132, 0.03900906, 0.03383408, 0.03517160, 0.02927114, 0.03888822, 0.03077891, 0.04236406,
0.05663730, 0.03619537, 0.04294887, 0.03497815, 0.03995837, 0.04374904, 0.03922274, 0.03596561, 0.03157820, 0.26390591, 0.06596527, 0.04050374,
0.02888223, 0.03824380, 0.05459656, 0.02969611, 0.86277224, 0.02385613, 0.03888451, 0.06496997, 0.03930725, 0.02931837, 0.06021005, 0.03330982,
0.02649659, 0.06600261, 0.02854480, 0.03691669, 0.06584168, 0.02076757, 0.02624355, 0.03679596, 0.03377049, 0.03590172, 0.03694652, 0.03575540,
0.02532416, 0.02818711, 0.04565318, 0.03252856, 0.04121822, 0.03147210, 0.05002047, 0.03809792, 0.02802299, 0.03399243, 0.03466543, 0.02829443,
0.03339476, 0.02129232, 0.03103367, 0.05071605, 0.03590172, 0.04386435, 0.03297124, 0.04323263, 0.03506247, 0.06225121, 0.02862442, 0.02862442,
0.06032925, 0.04400082, 0.03765090, 0.03477973, 0.02024540, 0.03564245, 0.05199116, 0.03699149, 0.03506247, 0.02129232, 0.02389752, 0.04996414,
0.04281258, 0.02587514, 0.03079668, 0.03895791, 0.02639014, 0.07333564, 0.02639014, 0.04074970, 0.04346211, 0.06032925, 0.03506247, 0.04950545,
0.04133673, 0.03835127, 0.02616212, 0.03399243, 0.02962473, 0.04800780, 0.03517160, 0.04105323, 0.03649472, 0.03000509, 0.05367187, 0.03858981,
0.03684529, 0.02941408, 0.04733265, 0.02590745, 0.02389752, 0.02385495, 0.03649472, 0.02508245, 0.02649659, 0.03152265, 0.02906310, 0.04950545,
0.03497815, 0.04374904, 0.03610649, 0.03799523, 0.02912771, 0.03694652, 0.05105353, 0.03000509, 0.02902378, 0.06425520, 0.05660319, 0.03065341,
0.04449069, 0.03638436, 0.02582273, 0.03753463, 0.02756006, 0.07215131, 0.02418869, 0.03431030, 0.04474425, 0.42589279, 0.02879489, 0.02872819,
0.02512494, 0.02450022, 0.03416346, 0.04560013, 1.40417366, 0.04784363, 0.04950545, 0.04685682, 0.03346052, 0.03255004, 0.07296053, 0.04491526,
0.02910482, 0.05448995, 0.01934736, 0.02195528, 0.03506247, 0.03157064, 0.03504810, 0.03754736, 0.03301058, 0.06886929, 0.03994190, 0.05130644,
0.21007323, 0.05630628, 0.02893721, 0.03683226, 0.03825290, 0.02494987, 0.02633410, 0.02721408, 0.03798986, 0.33473991, 0.04236406, 0.02389752,
0.03562747, 0.04662421, 0.02373767, 0.04918125, 0.04478894, 0.02418869, 0.03511514, 0.02871556, 0.05586166, 0.49014922, 0.03406339, 0.84823093,
0.03416346, 0.08729506, 0.03147210, 0.02889640, 0.06181828, 0.04940672, 0.03666858, 0.03019139, 0.03919279, 0.04864613, 0.03720420, 0.04726722,
0.04141298, 0.02862442, 0.29112744, 0.03964319, 0.05657445, 0.03930888, 0.04400082, 0.02722065, 0.03451685, 0.02911419, 0.02831578, 0.04001334,
0.05130644, 0.03134409, 0.03408579, 0.03232126, 0.03624218, 0.04708792, 0.06291741, 0.05663730, 0.03813209, 0.70582932, 0.04149421, 0.03607614,
0.03201109, 0.02055138, 0.03727305, 0.03182562, 0.02987404, 0.04142461, 0.03433624, 0.04264550, 0.02875086, 0.05797661, 0.04248705, 0.04476514))

根据上面的数据，我使用 pracma::findpeaks 函数获得了 22 个峰值，代码如下：

picos_r <- pracma::findpeaks(-valores$V1, minpeakdistance = 10)

使用 MATLAB 函数

picos_matlab = findpeaks(-dado_r, 'MinPeakWidth', 10);

我得到了11个峰，如下：

picos_matlab <- c(-0.02547804, -0.02618668, -0.01934736, -0.02182294, -0.0245572200000000, -0.0202454, -0.02385495, -0.01934736, -0.02373767, -0.02862442, -0.02722065)

我使用 pracma::findpeaks 是因为它已经在我正在编写的函数的另一部分给出了相同的结果。我已经尝试更改 pracma::findpeaks 的代码，但收效甚微。

【参考方案1】：

cardidates 包包含一个启发式峰值搜索算法，可以使用参数xmax、minpeak 和mincut 进行微调。它是为一个特殊的问题而设计的，但也可以用于其他事情。举个例子：

library("cardidates")

p <- peakwindow(valores$V1)
plot(p) # detects 14 peaks

p <- peakwindow(valores$V1, minpeak=0.18)
plot(p) # detects 11 peaks

详细信息在package vignette 和https://doi.org/10.1007/s00442-007-0783-2 中进行了描述

另一种选择是在峰值检测之前运行更平滑。

【讨论】：

谢谢！我要试试这个！此外，该值已经平滑。它们来自另一个函数。

如前所述，该算法与您对 Matlab 包的描述不同，并且具有其他调整启发式方法。 Cardidates 算法对我来说效果很好，但除此之外，我还对您在其他 R 包中找到的峰值检测算法感兴趣。但是，如果您需要精确复制 Matlab 算法，可能需要查看其源代码（如果可用）。

@tptzoldt，我需要 MATLAB 函数的精确复制。我查看了 MATLAB 的源代码，甚至复制了一个函数，但我仍然理解它。如果您有 MATLAB，它位于：计算机中的 MATLAB 文件夹 > 版本文件夹 > 工具箱 > 信号 > 信号 > findpeaks。

@tptzoldt，对我来说，处理复制值的函数是 pracma::findpeaks。 IDPmisc::peaks 也很有趣，因为它有一个 stepf 参数，可以更改以查找峰值。他们中的大多数使用 diff() 来检测峰值，但在我上面发布的线程中，有一些以不同的方式计算。

还有 'gsignal' 包，它是在 R 中重新实现的 'signal' Octave 包，看起来很有前途。它包含一个findpeaks() 函数。我没有尝试过，因为我不清楚您的主要测试用例是什么。

【参考方案2】：

我不确定您的测试用例是什么：-valores$V1、valores$V1 或 -dado_r（那是什么）？ 如果你这样做，我认为pracma::findpeaks() 做得很好：

x <- valores$V1
P <- pracma::findpeaks(x,
        minpeakdistance = 10, minpeakheight = sd(x))

plot(x, type = 'l', col = 4)
grid()
points(P[,2], P[, 1], pch=20, col = 2)

它会发现 11 个突出的峰，而其他四五个因太近而无法计数。所有较小的（标准差）都被忽略了。

【讨论】：

我的测试用例是-valores$V1。 -dado_r 是 MATLAB 中的变量名称，仅供参考。我将对此进行测试并提供反馈。感谢您的帮助！