Алгоритм для определения местоположения локальных максимумов
У меня есть данные, которые всегда выглядят примерно так:
сопроводительный текст 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
8 ответов
Локальным максимумом будет любая точка x, которая имеет большее значение y, чем любой из ее соседей слева и справа. Чтобы устранить шум, можно ввести некоторый порог допустимости (например, точка x должна иметь большее значение y, чем n ее соседей).
Чтобы не сканировать каждую точку, можно использовать тот же подход, но просмотреть 5-10 точек за раз, чтобы получить приблизительное представление о том, где находится максимум. Затем вернуться к этим областям для более детального сканирования.
Не могли бы вы двигаться вдоль графика, регулярно вычисляя производную, и если она переходит от положительного значения к отрицательному, вы нашли пик?
Знаете ли вы производную данных? Если да и вы можете символически решить систему, то вы можете найти точки, в которых производная равна нулю.
Если у вас нет формулы (что, по-видимому, так и есть, судя по вашему OP), то вы можете попробовать отфильтровать шум, а затем сделать следующее:
Если вы не можете решить символически, то вы можете использовать что-то вроде метода Ньютона-Рафсона, чтобы добраться до локальных максимумов и выбрать начальные точки случайным образом из диапазона, чтобы попытаться захватить все максимумы.
Если у вас нет данных о производной, вы можете попробовать алгоритм подъема по холму, который не требует производной, и начать его в нескольких различных случайно выбранных точках. Затем вы можете отслеживать точки, которые вы найдете, когда итеративная часть алгоритма завершится. Это только вероятностно найдет все локальные максимумы, но может быть достаточно хорошо для ваших целей.
EDIT: учитывая, что 3 вершины будут находиться примерно в одних и тех же местах, вы должны попытаться гарантировать, что начальная точка для этих алгоритмов находится рядом с этими точками, по крайней мере, в некоторых случаях, когда вы запускаете итерационный алгоритм.
Можно попробовать подогнать сплайн к данным, а затем найти экстремумы сплайна. Поскольку сплайн является кусочно-полиномиальным, точное расположение экстремумов можно найти с помощью относительно простых формул.
На практике я обнаружил, что хорошо работает использование операции морфологии расширения для создания расширенной версии вашей выборочной функции (точки данных), затем для определения локальных максимумов сравните расширенную версию с оригиналом и в любом месте, где расширенная версия равна оригинальной, должен быть локальный максимум. Это хорошо работает с 2D+ данными (т.е. изображениями), но поскольку у вас одномерные данные, может быть проще просто использовать разницу между последовательными точками как приближение к производной.
Обратите внимание, если вы используете метод расширения, структурирующий элемент (как размер, так и форма), который вы используете в расширении, в значительной степени определяет типы пиков, которые вы ищете.
Также, если в ваших данных есть шум, перед поиском сгладьте его с помощью фильтра низких частот, например, одномерного гауссова фильтра.
Больше информации о расширении можно найти здесь:
вот идея, реализованная в 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
Вы можете попробовать полосовой фильтр, чтобы отсеять шум и облегчить надежное выделение максимумов.
Смысл полосового фильтра (а не низкочастотного) в том, чтобы вытянуть почти постоянную до нуля. Таким образом, самые высокие значения, которые вы найдете в отфильтрованном результате, вероятно, будут самыми четкими пиками.
Конечно, если вы можете определить узкий частотный диапазон для вашего сигнала и применить очень избирательный фильтр, это должно сделать практичным довольно несложный алгоритм поиска максимумов - например, простое сканирование по принципу "образец-больше-чем-любой-сосед".
Сложный фильтр может и не понадобиться - можно обойтись мексиканским вейвлетом на одном хорошо выбранном масштабе. Один масштаб, вероятно, означает, что это уже не совсем вейвлет - просто полосовой FIR-фильтр.
EDIT
Существует асимметричный вейвлет (забыл название), который, если мексиканская шляпа является аналогом косинуса, берет на себя роль синуса. Я упоминаю его, поскольку он сочетает полосовую фильтрацию с производной - нулевые пересечения в результате являются стационарными точками полосовой фильтрованной версии исходного сигнала.
Затем "отлаженное" сканирование может идентифицировать надежные максимумы путем поиска точек пересечения в этом "производном" сигнале.
прямой подход, примерно так:
#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;
}
output:
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) если предыдущее значение больше текущего, то было максимальное:
распечатайте это и
установить состояние "Год назад"
3) в состоянии Godown дождитесь, пока предыдущее не станет меньше текущего, и переключитесь на (1)
, чтобы уменьшить шум, используйте некоторую функцию сглаживания smooth ()
Как уже упоминали другие, у меня обычно работают производные или сравнение с местными соседями. Если вас беспокоит шум, я могу порекомендовать медианную фильтрацию как довольно быструю и надежную схему фильтрации. Я все время использую его и обратную медианную фильтрацию, чтобы подавить шум в акустических датчиках, отлично работает.