Подтвердить что ты не робот

Алгоритм определения локальных максимумов

У меня есть данные, которые всегда выглядят примерно так:

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

Мне нужен алгоритм для поиска трех пиков.

Ось 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
4b9b3361

Ответ 1

Локальными максимумами была бы любая точка x, имеющая более высокое значение y, чем любой из ее левого и правого соседей. Чтобы устранить шум, вы можете установить некоторый порог допуска (например, точка x должна иметь более высокое значение y, чем n его соседей).

Чтобы избежать сканирования каждой точки, вы можете использовать один и тот же подход, но идти 5 или 10 точек за раз, чтобы получить грубое представление о том, где максимальная ложь. Затем вернитесь в эти области для более детального сканирования.

Ответ 2

Не могли бы вы двигаться по графику, регулярно вычисляя производную, и если он переключается с положительного на отрицательный, вы нашли пик?

Ответ 3

Вы можете попытаться поместить сплайн в данные, а затем найти экстремумы сплайна. Так как сплайн кусочно-полиномиальный, то точное расположение экстремумов можно найти с относительно простыми формулами.

Ответ 4

Я тренируюсь, я нашел то, что хорошо работает, это использовать операцию морфологии дилатации для создания расширенной версии вашей выборочной функции (точек данных), а затем для определения локального максимума сравните расширенную версию с оригиналом и где угодно, где расширенный вариант равен исходной версии должен быть локальным максимумом. Это работает очень хорошо, я нахожусь с данными 2D + (т.е. Изображениями), но поскольку у вас есть 1D-данные, может быть проще просто использовать различия между последовательными точками в качестве приближения к производной.

Обратите внимание, что если вы используете метод дилатации, структурирующий элемент (размер и форма), который вы используете в дилатации, значительно определяет типы пиков, которые вы ищете.

Также, если у вас есть шум в ваших данных, сгладьте его фильтрами нижних частот, как 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/dilate.htm

Ответ 5

прямой подход, что-то вроде этого:

#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) установить state = goup: идти вверх по кривой;
2) если значение previuos больше, чем текущее, тогда был максимум:
 распечатать и  установить состояние в будущее | 3), в то время как в состоянии godown дождитесь, пока предыдущий будет меньше текущего и переключится на (1)

для уменьшения шума используйте некоторую функцию сглаживания smooth()

Ответ 6

Знаете ли вы производную от данных? Если да, и вы можете символически решить систему, тогда вы можете найти точки, где производная равна нулю.

Если у вас нет формулы (которая, как представляется, имеет место из вашего OP), вы можете попробовать отфильтровать некоторый шум, а затем выполните следующие действия:

Если вы не можете решить символически, вы можете использовать что-то вроде метода Ньютона-Рафсона, чтобы добраться до локальных максимумов и случайно выбрать начальные точки из диапазона, чтобы попытаться захватить все максимумы.

Если у вас нет данных о производных, вы можете попробовать алгоритм альпинизма, который не требует производной и запускает его в нескольких разных случайно выбранных точках. Затем вы можете отслеживать точки, которые вы находите, когда итеративная альпинистская часть алгоритма заканчивается. Это будет только вероятностно найти все локальные максимумы, но это может быть достаточно хорошим для ваших целей.

EDIT: учитывая, что 3 пика будут примерно в тех же местах, которые вы должны попробовать, убедитесь, что начальная точка для этих алгоритмов находится рядом с этими точками, по крайней мере, в некоторые моменты, когда вы запускаете итеративный алгоритм.

Ответ 7

Вы можете попробовать полосовой фильтр, чтобы отклонить шум и облегчить надежно выбирать эти максимумы.

Точка полосы пропускания (а не низкочастотная) должна тянуть почти постоянную до нуля. Таким образом, самые высокие значения, которые вы найдете в отфильтрованном результате, вероятно, будут самыми ясными пиками.

Конечно, если вы можете определить узкий частотный диапазон для вашего сигнала и применить очень избирательный фильтр, он должен сделать довольно сложный алгоритм максимизации поиска - например, простое сканирование с образцом - более высокое, чем соседнее.

Возможно, вам не нужен сложный фильтр - вы можете уйти с помощью

Ответ 8

Как говорили другие, производные или по сравнению с местными соседями обычно работают для меня. Если вы беспокоитесь о шуме, я могу рекомендовать медианную фильтрацию как довольно быструю и надежную схему фильтрации. Я использую его и реверсирующую среднюю фильтрацию все время, чтобы шумоподавлять шумы в акустических датчиках, отлично работает.

Ответ 9

Другой метод - создать то, что я называю средним шагом ходьбы. Я не знаю, есть ли имя для него, но это похоже на то, что ваш набор данных, например, 1000 чисел, вы берете x [n] + x [n..1234567] номера скажем, 7 чисел ahaed, возьмите среднее число первых 3 и 3 последних используйте их, чтобы найти, будет ли линия, наложенная на них, идти вверх или вниз.

Когда он спускается, вы прошли номер горного загляда, после того, как один такой образец будет ждать, пока он снова не поднимется. Поэтому только после того, как вы поднимаете вверх, вы отслеживаете первый момент спуска. И повторите это.

Он обнаружит вершины и в зависимости от длины линии уклона (7).. 15..33 и т.д., вы также удалите шум.