Нахождение приближенных локальных максимумов с зашумленными данными в Matlab

Question

Часто задаваемые вопросы по Matlab описывают однострочный метод для нахождения локальных максимумов:

index = find( diff( sign( diff([0; x(:); 0]) ) ) < 0 );

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

Было бы неплохо также однострочное решение.

Редактировать: Я работаю с шумными биологическими изображениями, которые я пытаюсь разделить на отдельные разделы.

Answer 1

Я не уверен, с каким типом данных вы имеете дело, но вот метод, который я использовал для обработки речевых данных, который может помочь вам найти локальные максимумы. Он использует три функции из панели инструментов обработки сигналов: ГИЛЬБЕРТ , БАБОЧКА и FILTFILT .

data = (...the waveform of noisy data...);
Fs = (...the sampling rate of the data...);
[b,a] = butter(5,20/(Fs/2),'low');  % Create a low-pass butterworth filter;
                                    %   adjust the values as needed.
smoothData = filtfilt(b,a,abs(hilbert(data)));  % Apply a hilbert transform
                                                %   and filter the data.

Затем вы выполните поиск максимума для smoothData . Использование HILBERT сначала создает положительный конверт в данных, затем FILTFILT использует коэффициенты фильтра от BUTTER для низкочастотной фильтрации конверта данных.

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

ЧТО-ТО СЛУЧАЙНО ПОПРОБОВАТЬ:

Сначала это была случайная идея ... вы можете попытаться повторить процесс, найдя максимумы максим:

index = find(diff(sign(diff([0; x(:); 0]))) < 0);
maxIndex = index(find(diff(sign(diff([0; x(index); 0]))) < 0));

Однако, в зависимости от отношения сигнал / шум, было бы неясно, сколько раз это нужно будет повторить, чтобы получить интересующие вас локальные максимумы. Это просто случайная опция без фильтрации. =) * * Один тысяча двадцать-шесть

НАЙТИ МАКСИМА:

На всякий случай вам может понравиться еще один алгоритм нахождения максимума в одну строку, который я видел (в дополнение к указанному вами):

index = find((x > [x(1) x(1:(end-1))]) & (x >= [x(2:end) x(end)]));

Answer 2

В зависимости от того, что вы хотите, часто полезно фильтровать зашумленные данные. Взгляните на MEDFILT1 или используйте CONV вместе с FSPECIAL . В последнем подходе вы, вероятно, захотите использовать «тот же» аргумент для CONV и «гауссов» фильтр, созданный FSPECIAL.

После того, как вы выполните фильтрацию, пропустите его через поиск максимума.

РЕДАКТИРОВАТЬ: сложность во время выполнения

Допустим, входной вектор имеет длину X, а длина ядра фильтра - K.

Медианный фильтр может работать, выполняя сортировку вставки, поэтому он должен быть O (X K + K log K). Я не смотрел на исходный код и возможны другие реализации, но в основном это должно быть O (X K).

Когда K мало, conv использует прямой алгоритм O (X * K). Когда X и K почти одинаковы, тогда быстрее использовать быстрое преобразование Фурье. Эта реализация O (X log X + K log K). Matlab достаточно умен, чтобы автоматически выбирать правильный алгоритм в зависимости от размеров ввода.

Answer 3

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

Другой подход к такой проблеме - тот, который я предпочитаю. Здесь используется сплайн-модель наименьших квадратов для сглаживания данных и получения аппроксимации вместо интерполяции. Такой сплайн по методу наименьших квадратов имеет то преимущество, что позволяет пользователю с большой степенью контроля вводить свои знания о проблеме в модель. Например, часто ученый или инженер имеет информацию, такую ​​как монотонность, об изучаемом процессе. Это может быть встроено в модель сплайна наименьших квадратов. Другим связанным вариантом является использование сглаживающего сплайна. Они также могут быть построены с помощью регуляризирующих ограничений. Если у вас есть набор инструментов сплайнов, то spap2 будет полезен для сплайн-модели. Тогда fnmin найдет минимизатор. (Максимизатор легко получается из кода минимизации.)

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

Answer 4

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

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