Нахождение положительных пиков функции

#python #алгоритм #обработка_сигналов

Требуется программно найти количество и местоположение положительных пиков на графике.
Известно, что количество пиков нечетное. Обычно 3, 5 или 7. В данном случае 3. Также
известно, что за каждым положительным пиком следует отрицательный.
Как лучше решить такую проблему на Python? Желательно без подбора коэффициентов
Дополнение: Данные всегда начинаются с положительного пика и заканчиваются отрицательным.
Данные представляют собой вертикальную проекцию изображения, обработанную функцией
np.gradient(projection)/projection. То есть отношение второй производной функции к
её высоте.
    

Ответы

Ответ 1

Поскольку просто локальный экстремум не устраивает, надо добавить проверку превышения
порога для отбраковки мелких пиков и артефактов. Порог можно вычислить как СКО выборки,
умноженное на некоторую константу, соответствующую доверительному интервалу для математического
ожидания.
На практике эту константу принимают равной 2.5.

UPD
С учётом дополнительных условий задачи пиком следует считать максимум на отрезке
от "верха" (начальной точки или первого значения выше q) до "низа" (первого значения
ниже -q).
При этом для массива размерности n возможен следующий алгоритм:


Найти СКО выборки и вычислить порог q.
Положить top=0.
Найти координату bottom первой точки отрезка [top, n-1] со значением меньше -q. Если
такой точки нет, перейти к п.7.
Найти координату максимума на отрезке [top, bottom-1]. Если значение максимума больше
q, записать его в массив пиков.
Найти координату top первой точки отрезка [bottom+1, n-1] со значением больше q.
Если такой точки нет, перейти к п.7.
Перейти к п.3.
Использовать массив пиков.


Ответ 2

Ничего лучшего, чем  поиск пиков в массиве в голову не приходит. Строгое чередование
положительных и отрицательных пиков облегчает задачу.

Т.е. алгоритм примерно такой.
Выбираем значения двух порогов: положительного и отрицательного. Ищем максимум пока
сигнал больше отрицательного порога. Найденный максимум и есть пик, а индекс в массиве
- его позиция. Ждем превышения положительного порога и ищем очередной положительный
пик. Пороги подбираются так, чтобы быть выше шума: колебания около нуля. Их значения
могут быть ясны исходя из природы сигнала (минимально возможная величина пика), либо
можно попытаться их вычислить статистически (найти среднеквадратическое отклонение
шума и установить порог в значение, например 3 СКО).
Боюсь, что готовой функции в библиотеках не будет. Однако реализация такого алгоритма
тривиальна. К сожалению, с Python мало знаком, поэтому приведу возможную реализацию
на C++ (думаю, что проблем с пониманием возникнуть не должно):

typedef double sample_type;
struct Peak {
    sample_type magnitude;
    std::size_t pos;
};
std::vector peaks(const std::vector &sig) {
    static sample_type threshold_lo = -0.05;
    static sample_type threshold_hi = 0.05;
    std::vector ret_peaks;
    std::vector::size_type isample = 0;
    while (1) {
        while ((isample < sig.size()) && (sig[isample] < threshold_hi)) {
            ++isample;
        }
        if (isample >= sig.size())
            break;
        ret_peaks.push_back(Peak(sig[isample], isample));
        while ((isample < sig.size()) && (sig[isample] > threshold_lo)) {
            if (sig[isample] > ret_peaks.back().magnitude) {
                ret_peaks.back().magnitude = sig[isample];
                ret_peaks.back().pos = isample;
            }
            ++isample;
        }
        if (isample >= sig.size())
            break;
    }
    return ret_peaks;
}


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