Обнаружение пиков измеряемого сигнала
Мы используем карту сбора данных для считывания данных с устройства, которое увеличивает свой сигнал до пика, а затем возвращается к исходному значению. Чтобы найти пиковое значение, мы в настоящее время ищем в массиве наибольшее значение и используем индекс для определения времени пикового значения, которое используется в наших расчетах.
Это хорошо работает, если наибольшее значение - это пик, который мы ищем, но если устройство не работает правильно, мы можем увидеть второй пик, который может быть выше, чем начальный пик. Мы принимаем 10 показаний в секунду с 16 устройств в течение 90 секунд.
Мои первые мысли заключаются в том, чтобы циклически проверять показания, чтобы увидеть, меньше ли предыдущая и следующая точки, чем текущая, чтобы найти пик и построить массив пиков. Возможно, мы должны смотреть на среднее количество точек по обе стороны от текущей позиции, чтобы учесть шум в системе. Это лучший способ продолжить или есть лучшие методы?
Мы используем LabVIEW, и я проверил форумы LAVA, и есть много интересных примеров. Это часть нашего тестового программного обеспечения, и мы стараемся избегать использования слишком большого количества нестандартных библиотек VI, поэтому я надеялся получить отзывы о задействованных процессах / алгоритмах, а не о конкретном коде.
9 ответов
Вы можете попробовать усреднение сигнала, то есть для каждой точки усредните значение с окружающими 3 или более точками. Если шумовые помехи огромны, то даже это может не помочь.
Я понимаю, что это не зависело от языка, но, предположив, что вы используете LabView, есть много готовых VI для обработки сигналов, которые поставляются с LabView, которые вы можете использовать для сглаживания и уменьшения шума. Форумы NI - отличное место, чтобы получить более специализированную помощь по этому вопросу.
Существует множество классических методов обнаружения пиков, любой из которых может сработать. Вы должны увидеть, что, в частности, ограничивает качество ваших данных. Вот основные описания:
Между любыми двумя точками в ваших данных,
(x(0), y(0))а также(x(n), y(n)), сложитьy(i + 1) - y(i)за0 <= i < nи называть этоT("путешествие") и установитьR("подняться наy(n) - y(0) + kдля достаточно маленькихk,T/R > 1указывает на пик. Это работает нормально, если большое перемещение из-за шума маловероятно или если шум распределяется симметрично вокруг формы базовой кривой. Для вашего приложения примите самый ранний пик с оценкой выше заданного порогового значения или проанализируйте кривую перемещения по нарастанию для более интересных свойств.Используйте согласованные фильтры для оценки сходства со стандартной формой пика (по существу, используйте нормализованный точечный продукт для некоторой формы, чтобы получить косинометрическую метрику сходства)
Деконволюция по стандартной форме пика и проверка на высокие значения (хотя я часто нахожу 2, чтобы быть менее чувствительными к шуму для простого вывода инструментов).
Сгладить данные и проверить наличие тройки точек с одинаковым интервалом, если, если
x0 < x1 < x2, y1 > 0.5 * (y0 + y2)или проверьте евклидовы расстояния следующим образом:D((x0, y0), (x1, y1)) + D((x1, y1), (x2, y2)) > D((x0, y0),(x2, y2)), который опирается на неравенство треугольника. Использование простых соотношений снова предоставит вам механизм подсчета очков.Подберите к своим данным очень простую 2-гауссову модель смешения (например, в Numeric Recipes есть хороший готовый кусок кода). Возьмите более ранний пик. Это будет правильно работать с перекрывающимися пиками.
Найдите наилучшее совпадение в данных с простой кривой Гаусса, Коши, Пуассона или "что у вас есть". Оцените эту кривую в широком диапазоне и вычтите ее из копии данных, отметив ее пиковое местоположение. Повторение. Возьмем самый ранний пик, чьи параметры модели (вероятно, стандартное отклонение, но некоторые приложения могут заботиться о куртозе или других особенностях) соответствуют некоторому критерию. Остерегайтесь артефактов, оставшихся позади, когда пики вычитаются из данных. Наилучшее совпадение может быть определено видом подсчета совпадений, предложенным в #2 выше.
Я делал то, что вы делаете раньше: находя пики в данных последовательности ДНК, находя пики в производных, оцененных по измеренным кривым, и находя пики в гистограммах.
Я призываю вас внимательно следить за надлежащим базовым уровнем. Фильтрация Винера или другая фильтрация или простой анализ гистограммы часто являются простым способом определения исходного уровня при наличии шума.
Наконец, если ваши данные, как правило, зашумлены и вы получаете данные с карты как односторонний вывод без ссылок (или даже со ссылками, просто без разницы), и если вы усредняете множество наблюдений в каждой точке данных, попробуйте отсортировать их наблюдения и отбрасывание первого и последнего квартиля и усреднение того, что осталось. Есть множество таких тактик исключения, которые могут быть действительно полезными.
Я хотел бы внести в этот поток алгоритм, который я разработал сам:
Он основан на принципе дисперсии: если новая точка данных представляет собой заданное число x стандартных отклонений от некоторого скользящего среднего, алгоритм подает сигнал (также называемый z-счетом). Алгоритм очень надежен, потому что он строит отдельное скользящее среднее и отклонение, так что сигналы не нарушают порог. Поэтому будущие сигналы идентифицируются примерно с одинаковой точностью, независимо от количества предыдущих сигналов. Алгоритм занимает 3 входа: lag = the lag of the moving window, threshold = the z-score at which the algorithm signals а также influence = the influence (between 0 and 1) of new signals on the mean and standard deviation, Например, lag из 5 будет использовать последние 5 наблюдений для сглаживания данных. threshold 3,5 будет сигнализировать, если точка данных находится в 3,5 стандартных отклонениях от скользящего среднего. И influence 0,5 дает сигналы половину влияния, которое имеют нормальные точки данных. Аналогично influence "0" полностью игнорирует сигналы для пересчета нового порога: поэтому влияние "0" является наиболее надежным вариантом.
Это работает следующим образом:
ПСЕВДОКОД
# Let y be a vector of timeseries data of at least length lag+2
# Let mean() be a function that calculates the mean
# Let std() be a function that calculates the standard deviaton
# Let absolute() be the absolute value function
# Settings (the ones below are examples: choose what is best for your data)
set lag to 5; # lag 5 for the smoothing functions
set threshold to 3.5; # 3.5 standard deviations for signal
set influence to 0.5; # between 0 and 1, where 1 is normal influence, 0.5 is half
# Initialise variables
set signals to vector 0,...,0 of length of y; # Initialise signal results
set filteredY to y(1,...,lag) # Initialise filtered series
set avgFilter to null; # Initialise average filter
set stdFilter to null; # Initialise std. filter
set avgFilter(lag) to mean(y(1,...,lag)); # Initialise first value
set stdFilter(lag) to std(y(1,...,lag)); # Initialise first value
for i=lag+1,...,t do
if absolute(y(i) - avgFilter(i-1)) > threshold*stdFilter(i-1) then
if y(i) > avgFilter(i-1)
set signals(i) to +1; # Positive signal
else
set signals(i) to -1; # Negative signal
end
# Adjust the filters
set filteredY(i) to influence*y(i) + (1-influence)*filteredY(i-1);
set avgFilter(i) to mean(filteredY(i-lag,i),lag);
set stdFilter(i) to std(filteredY(i-lag,i),lag);
else
set signals(i) to 0; # No signal
# Adjust the filters
set filteredY(i) to y(i);
set avgFilter(i) to mean(filteredY(i-lag,i),lag);
set stdFilter(i) to std(filteredY(i-lag,i),lag);
end
end
демонстрация
> Оригинальный ответ
Эта проблема была изучена в некоторых деталях.
Есть множество очень современных реализаций в TSpectrum* классы ROOT (инструмент анализа ядерной физики / физики элементарных частиц). Код работает в одно- и трехмерных данных.
Исходный код ROOT доступен, так что вы можете получить эту реализацию, если хотите.
Из документации класса TSpectrum:
Алгоритмы, используемые в этом классе, были опубликованы в следующих ссылках:
[1] М. Морхак и др.: Методы устранения фона для многомерных спектров гамма-излучения совпадений. Ядерные приборы и методы в физических исследованиях A 401 (1997) 113-132.
[2] М. Морхак и др.: Эффективная одно- и двумерная деконволюция золота и ее применение к разложению гамма-спектров. Ядерные приборы и методы в физике исследований А 401 (1997) 385-408.
[3] М. Морхац и др.: Идентификация пиков в многомерных спектрах гамма-излучения совпадений. Ядерные приборы и методы в физике исследований А 443(2000), 108-125.
Документы связаны с документацией класса для тех из вас, у кого нет онлайн-подписки NIM.
Короткая версия того, что сделано, заключается в том, что гистограмма сглаживается для устранения шума, а затем локальные максимумы обнаруживаются с помощью грубой силы в сглаженной гистограмме.
Этот метод в основном из книги Дэвида Марра "Видение"
Gaussian Blur ваш сигнал с ожидаемой шириной ваших пиков. это избавляет от всплесков шума и ваши данные фазы не повреждены.
Тогда обнаружение края (LOG сделает)
Тогда ваши края были краями элементов (например, вершины). ищите пики между краями, сортируйте пики по размеру, и все готово.
Я использовал варианты этого, и они работают очень хорошо.
Я думаю, что вы хотите взаимно коррелировать ваш сигнал с ожидаемым, образцовым сигналом. Но прошло много времени с тех пор, как я изучал обработку сигналов, и даже тогда я не обращал особого внимания.
Есть ли качественная разница между желаемым пиком и нежелательным вторым пиком? Если оба пика являются "острыми" - т.е. короткими по времени - при просмотре сигнала в частотной области (выполняя БПФ) вы получите энергию в большинстве полос. Но если в "хорошем" пике надежно присутствует энергия на частотах, которых нет в "плохом" пике, или наоборот, вы можете автоматически их дифференцировать таким образом.
Вы можете применить стандартное отклонение к вашей логике и заметить пики свыше x%.
Я не очень много знаю об инструментах, так что это может быть совершенно непрактично, но опять же, это может быть полезным в другом направлении. Если вы знаете, как показания могут потерпеть неудачу, и существует определенный интервал между пиками при таких сбоях, почему бы не сделать градиентный спуск на каждом интервале. Если спуск вернет вас в область, которую вы искали ранее, вы можете отказаться от нее. В зависимости от формы выбранной поверхности это также может помочь вам находить пики быстрее, чем поиск.