定位局部最大值的算法

Question

我的数据总是看起来像这样：

alt text http://michaelfogleman.com/static/images/chart.png< /a>

我需要一种算法来定位三个峰值。

x 轴实际上是相机位置，y 轴是该位置处图像焦点/对比度的度量。三个不同距离处的特征可以聚焦，我需要确定这三个点的 x 值。

中间的驼峰总是比较难辨认，即使对于人类来说也是如此。

我有一个大部分有效的自制算法，但我想知道是否有一种标准方法可以从可能有一点噪音的函数中获取局部最大值。不过，峰值很容易克服噪音。

此外，作为相机数据，不需要扫描整个范围的算法可能会很有用。

编辑：发布我最终使用的Python代码。它使用我的原始代码来查找给定搜索阈值的最大值，并进行二分搜索以找到导致所需最大值数量的阈值。

编辑：下面的代码中包含示例数据。新代码的复杂度是 O(n) 而不是 O(n^2)。

def find_n_maxima(data, count):
    low = 0
    high = max(data) - min(data)
    for iteration in xrange(100): # max iterations
        mid = low + (high - low) / 2.0
        maxima = find_maxima(data, mid)
        if len(maxima) == count:
            return maxima
        elif len(maxima) < count: # threshold too high
            high = mid
        else: # threshold too low
            low = mid
    return None # failed

def find_maxima(data, threshold):
    def search(data, threshold, index, forward):
        max_index = index
        max_value = data[index]
        if forward:
            path = xrange(index + 1, len(data))
        else:
            path = xrange(index - 1, -1, -1)
        for i in path:
            if data[i] > max_value:
                max_index = i
                max_value = data[i]
            elif max_value - data[i] > threshold:
                break
        return max_index, i
    # forward pass
    forward = set()
    index = 0
    while index < len(data) - 1:
        maximum, index = search(data, threshold, index, True)
        forward.add(maximum)
        index += 1
    # reverse pass
    reverse = set()
    index = len(data) - 1
    while index > 0:
        maximum, index = search(data, threshold, index, False)
        reverse.add(maximum)
        index -= 1
    return sorted(forward & reverse)

data = [
    1263.900, 1271.968, 1276.151, 1282.254, 1287.156, 1296.513,
    1298.799, 1304.725, 1309.996, 1314.484, 1321.759, 1323.988,
    1331.923, 1336.100, 1340.007, 1340.548, 1343.124, 1353.717,
    1359.175, 1364.638, 1364.548, 1357.525, 1362.012, 1367.190,
    1367.852, 1376.275, 1374.726, 1374.260, 1392.284, 1382.035,
    1399.418, 1401.785, 1400.353, 1418.418, 1420.401, 1423.711,
    1425.214, 1436.231, 1431.356, 1435.665, 1445.239, 1438.701,
    1441.988, 1448.930, 1455.066, 1455.047, 1456.652, 1456.771,
    1459.191, 1473.207, 1465.788, 1488.785, 1491.422, 1492.827,
    1498.112, 1498.855, 1505.426, 1514.587, 1512.174, 1525.244,
    1532.235, 1543.360, 1543.985, 1548.323, 1552.478, 1576.477,
    1589.333, 1610.769, 1623.852, 1634.618, 1662.585, 1704.127,
    1758.718, 1807.490, 1852.097, 1969.540, 2243.820, 2354.224,
    2881.420, 2818.216, 2552.177, 2355.270, 2033.465, 1965.328,
    1824.853, 1831.997, 1779.384, 1764.789, 1704.507, 1683.615,
    1652.712, 1646.422, 1620.593, 1620.235, 1613.024, 1607.675,
    1604.015, 1574.567, 1587.718, 1584.822, 1588.432, 1593.377,
    1590.533, 1601.445, 1667.327, 1739.034, 1915.442, 2128.835,
    2147.193, 1970.836, 1755.509, 1653.258, 1613.284, 1558.576,
    1552.720, 1541.606, 1516.091, 1503.747, 1488.797, 1492.021,
    1466.720, 1457.120, 1462.485, 1451.347, 1453.224, 1440.477,
    1438.634, 1444.571, 1428.962, 1431.486, 1421.721, 1421.367,
    1403.461, 1415.482, 1405.318, 1399.041, 1399.306, 1390.486,
    1396.746, 1386.178, 1376.941, 1369.880, 1359.294, 1358.123,
    1353.398, 1345.121, 1338.808, 1330.982, 1324.264, 1322.147,
    1321.098, 1313.729, 1310.168, 1304.218, 1293.445, 1285.296,
    1281.882, 1280.444, 1274.795, 1271.765, 1266.857, 1260.161,
    1254.380, 1247.886, 1250.585, 1246.901, 1245.061, 1238.658,
    1235.497, 1231.393, 1226.241, 1223.136, 1218.232, 1219.658,
    1222.149, 1216.385, 1214.313, 1211.167, 1208.203, 1206.178,
    1206.139, 1202.020, 1205.854, 1206.720, 1204.005, 1205.308,
    1199.405, 1198.023, 1196.419, 1194.532, 1194.543, 1193.482,
    1197.279, 1196.998, 1194.489, 1189.537, 1188.338, 1184.860,
    1184.633, 1184.930, 1182.631, 1187.617, 1179.873, 1171.960,
    1170.831, 1167.442, 1177.138, 1166.485, 1164.465, 1161.374,
    1167.185, 1174.334, 1186.339, 1202.136, 1234.999, 1283.328,
    1347.111, 1679.050, 1927.083, 1860.902, 1602.791, 1350.454,
    1274.236, 1207.727, 1169.078, 1138.025, 1117.319, 1109.169,
    1080.018, 1073.837, 1059.876, 1050.209, 1050.859, 1035.003,
    1029.214, 1024.602, 1017.932, 1006.911, 1010.722, 1005.582,
    1000.332, 998.0721, 992.7311, 992.6507, 981.0430, 969.9936,
    972.8696, 967.9463, 970.1519, 957.1309, 959.6917, 958.0536,
    954.6357, 954.9951, 947.8299, 953.3991, 949.2725, 948.9012,
    939.8549, 940.1641, 942.9881, 938.4526, 937.9550, 929.6279,
    935.5402, 921.5773, 933.6365, 918.7065, 922.5849, 939.6088,
    911.3251, 923.7205, 924.8227, 911.3192, 936.7066, 915.2046,
    919.0274, 915.0533, 910.9783, 913.6773, 916.6287, 907.9267,
    908.0421, 908.7398, 911.8401, 914.5696, 912.0115, 919.4418,
    917.0436, 920.5495, 917.6138, 907.5037, 908.5145, 919.5846,
    917.6047, 926.8447, 910.6347, 912.8305, 907.7085, 911.6889,
]

for n in xrange(1, 6):
    print 'Looking for %d maxima:' % n
    indexes = find_n_maxima(data, n)
    print indexes
    print ', '.join(str(data[i]) for i in indexes)
    print

输出：

Looking for 1 maxima:
[78]
2881.42

Looking for 2 maxima:
[78, 218]
2881.42, 1927.083

Looking for 3 maxima:
[78, 108, 218]
2881.42, 2147.193, 1927.083

Looking for 4 maxima:
[78, 108, 218, 274]
2881.42, 2147.193, 1927.083, 936.7066

Looking for 5 maxima:
[78, 108, 218, 269, 274]
2881.42, 2147.193, 1927.083, 939.6088, 936.7066

Answer 1

正如其他人提到的衍生品或与当地邻居进行比较通常对我有用。如果您担心噪音，我可以推荐中值过滤作为一种非常快速且可靠的过滤方案。我一直使用它和反向中值过滤来抑制声学传感器中的噪声，效果很好。

Answer 2

我想我在解决复利公式时也遇到过类似的问题。我认为如果点之间没有拐点，您可以尝试通过迭代获得局部最大值或最小值。

如果你不知道这个函数是怎么回事，我想除了求导之外没有解决办法。

我的解决方案将线段分成 4 个相等的部分（5 个点），并比较 3 个连续点的组以选择一个新的小线段。也许这是一些更好的实现。

//I know always f(-1)<f(1) and no inflection point between -1 and 1.


  function maximo(capital_inicial,periodos, aportacion){
  var min=-1;var max=1;var c=[];
  c[0]=capitalfinal(capital_inicial,periodos,min,aportacion);
  c[2]=capitalfinal(capital_inicial,periodos,(min+max)/2,aportacion);
  c[4]=capitalfinal(capital_inicial,periodos,max,aportacion);
  while(max-min>0.0001){
    c[1]=capitalfinal(capital_inicial,periodos,(3*min+max)/4,aportacion);
    if(segmento(c[0],c[1],c[2])==-1){
      c[3]=capitalfinal(capital_inicial,periodos,(min+max*3)/4,aportacion);
      if(segmento(c[1],c[2],c[3])==-1){
        min=(min+max)/2;
        c[0]=c[2];
        c[2]=c[3];
      }else{
        min=(3*min+max)/4;
        max=(min+max*3)/4
        c[0]=c[1];
        c[4]=c[3];
      }
    }else{
      max=(min+max)/2;
      c[4]=c[2];
      c[2]=c[1];
    }
  }
  return (max+min)/2;
}
function segmento(a,b,c){
  return a>b && b>c ? -1 : a<b && b<c ? 1:0;
}

https://docs.google.com/spreadsheets/ d/1crPdZfsOSbkfqya9C7Ezg-BZcVZzf46_avX-NVPsaqw/edit#gid=2021465223

Answer 3

您可以按两个具有不同 alpha 的指数平均值进行过滤。然后在两个交叉点之间获取数据，其中数据的值高于交叉点。然后你可以使用点的导数。

Answer 4

局部最大值是任何具有比其左右邻居更高 y 值的 x 点。为了消除噪声，您可以输入某种容差阈值（例如，x 点的 y 值必须高于其相邻的 n 个点）。

为了避免扫描每个点，您可以使用相同的方法，但一次扫描 5 或 10 个点，以大致了解最大值在哪里。然后返回这些区域进行更详细的扫描。

Answer 5

难道你不能沿着图表移动，定期计算导数，如果它从正值变为负值，你就找到了峰值吗？

Answer 6

另一种方法是创建我所说的步行坡度平均值。
我不知道它是否有一个名称，但它是这样的，您的数据集例如是 1000 个数字，您采用 x[n] +x[n..1234567] 数字
说出 7 个数字，取前 3 个和后 3 个的平均值
用它们来确定放在它们上面的线是向上还是向下。

当它下降时，你会经过一个山峰数字，在一个这样的样本之后，继续等待，直到它再次上升。因此，只有在上升之后，您才能追踪下降的第一个时刻。并重复一遍。
可能会像下面这样编码（python），但其他语言也是可能的，只需坚持这个想法。

import numpy as np

def slope_detection(x, window_size):
  """
  This function detects the slope of a given signal and identifies peaks.

  Args:
    x: A numpy array of the signal.
    window_size: The number of samples to use for averaging.

  Returns:
    A list of tuples, where each tuple contains the index of a peak and its slope.
  """

  peaks = []

  # Iterate through the signal, starting from the window_size-th sample
  for i in range(window_size, len(x)):
    # Calculate the average of the previous window_size samples
    avg_prev = np.mean(x[i-window_size:i])

    # Calculate the average of the next window_size samples
    avg_next = np.mean(x[i:i+window_size])

    # Calculate the slope
    slope = avg_next - avg_prev

    # Identify slope direction and peak
    if slope > 0:
      # Slope up
      if x[i] > avg_prev:
        peaks.append((i, "up"))
    elif slope < 0:
      # Slope down
      if x[i] < avg_next:
        peaks.append((i, "down"))
    else:
      # Slope is 0, consider it a peak if it's the highest point in the vicinity
      if i > window_size and i < len(x) - window_size and x[i] >= x[i-window_size:i+window_size].max():
        peaks.append((i, "peak"))

  return peaks

# Example usage
x = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1])
window_size = 3
peaks = slope_detection(x, window_size)

print("Peaks:")
for peak in peaks:
  print(f"  - Index: {peak[0]}, Slope: {peak[1]}")


slope detection with an avg function, can be ideal for more noisy data (like stock trade data, and senorry data)

It will detect tops, and depending on slope line length (7) .. 15 ..33 etc, you also remove noise. 

Answer 7

您可以尝试将样条线拟合到数据，然后找到样条线的极值。由于样条是分段多项式，因此可以使用相对简单的公式找到极值的确切位置。

Answer 8

我在实践中发现，效果好的方法是使用扩张形态运算来生成采样函数（数据点）的扩张版本，然后识别局部最大值，将扩张版本与原始版本以及扩张版本等于的任何地方进行比较原始版本应该是局部最大值。我发现这对于 2D+ 数据（即图像）非常有效，但由于您有 1D 数据，因此使用连续点之间的差异作为导数的近似值可能会更容易。

请注意，如果您确实使用膨胀技术，则在膨胀中使用的结构元素（大小和形状）很大程度上决定了您要查找的峰的类型。

此外，如果数据中存在噪声，请在搜索之前使用低通滤波器（例如一维高斯滤波器）对其进行平滑处理。

有关膨胀的更多信息可以在这里找到：

这是在 matlab 中实现的想法： http://www.mathworks.com/matlabcentral/fileexchange/14498-local-maxima-minima

如果您不知道什么是膨胀：
http://en.wikipedia.org/wiki/Dilation_%28morphology%29

（一旦你理解了，它就非常简单了，这是一个非常简单的解释）
http://homepages.inf.ed.ac.uk/rbf/ HIPR2/扩张.htm

Answer 9

直接方法，如下所示：

#include <stdio.h>
#include <stdlib.h>

#define MAXN 100

double smooth(double arr[], int n, int i)
{
        double l,r,smoo;
        l = (i - 1 < 0)?arr[0]:arr[i-1];
        r = (i + 1 >= n)?arr[n-1]:arr[i+1];
        smoo = (l + 2*arr[i] + r)/4;
        return smoo;
}

void findmax(double arr[], int n)
{
        double prev = arr[0];
        int i;
        int goup = 1;

        for(i = 0; i < n; i++)
        {
                double cur = smooth(arr,n,i);
                if(goup) {
                        if(prev > cur && i > 0) {
                                printf("max at %d = %lf\n", i-1, arr[i-1]);
                                goup = 0;
                        }
                } else {
                        if(prev < cur)
                                goup = 1;
                }
                prev = cur;
        }
}

int main()
{
        double arr[MAXN] = {0,0,1,2,3,4,4,3,2,2,3,4,6,8,7,8,6,3,1,0};
        int n = 20, i;

        for(i = 0; i < n; i++)
                printf("%.1lf ",arr[i]);
        printf("\n");

        findmax(arr,n);
        return 0;
}

输出：

0.0 0.0 1.0 2.0 3.0 4.0 4.0 3.0 2.0 2.0 3.0 4.0 6.0 8.0 7.0 8.0 6.0 3.0 1.0 0.0
max at 6 = 4.000000
max at 14 = 7.000000

1）set state = goup：沿着曲线向上；
2) 如果先前的值大于当前值，则存在最大值：
打印它并
将状态设置为 godown
3）在godown状态下等待，直到前一个小于当前并切换到（1）

以减少噪音，使用一些平滑函数smooth()

Answer 10

你知道数据的导数吗？
如果是，并且您可以象征性地求解该系统，那么您可以找到导数等于零的点。

如果你没有公式（这似乎是你的OP中的情况），那么你可能想尝试过滤掉一些噪音，然后执行以下操作：

如果你无法符号解决，那么你可以使用牛顿-拉夫森之类的东西方法获得局部最大值并从范围中随机选择起点以尝试捕获所有最大值。

如果您没有导数数据，那么您可能想尝试不需要导数的爬山算法，并从多个不同的随机选择点开始。然后，您可以跟踪算法的迭代爬山部分终止时找到的点。这只会在概率上找到所有局部最大值，但它可能足以满足您的目的。

编辑：考虑到 3 个峰值将大致位于相同的位置，您应该尝试保证这些算法的起点至少在运行迭代算法的某些时候接近这些点。

Answer 11

您可以尝试使用带通滤波器来抑制噪声，并更轻松地可靠地选择这些最大值。

带通（而不是低通）的目的是将近常数拉至零。这样，您在过滤结果中找到的最高值可能是最清晰的峰值。

当然，如果您可以为信号定义一个较窄的频率范围并应用一个非常有选择性的滤波器，那么它应该使一个相当简单的最大值查找算法变得实用 - 例如，一个简单的样本扫描高于任何一个邻居扫描。

复杂的过滤器可能不是必需的 - 您可以一次使用 墨西哥帽子小波精心挑选的规模。一个尺度可能意味着它不再是真正的小波——只是一个带通 FIR 滤波器。

编辑

有一个不对称小波（我忘了名字），如果墨西哥帽类似于余弦，它就扮演正弦的角色。我提到它是因为它结合了带通滤波和一种导数 - 结果中的零交叉是原始信号的带通滤波版本的驻点。

然后，“去抖”扫描可以通过寻找该“导数”信号中的交叉点来识别可靠的最大值。