Мы используем карту сбора данных для считывания показаний с устройства, которое увеличивает свой сигнал до пика, а затем возвращается к исходному значению. Чтобы найти пиковое значение, мы в настоящее время ищем в массиве наибольшее значение и используем индекс для определения времени пикового значения, которое используется в наших расчетах.
Это хорошо работает, если наибольшее значение - это пик, который мы ищем, но если устройство работает неправильно, мы можем увидеть второй пик, который может быть выше, чем начальный пик. Мы принимаем 10 показаний в секунду с 16 устройств в течение 90 секунд.
Мои первые мысли заключаются в том, чтобы циклически проверять показания, чтобы увидеть, меньше ли предыдущая и следующая точки, чем текущая, чтобы найти пик и построить массив пиков. Возможно, мы должны смотреть на среднее число точек по обе стороны от текущей позиции, чтобы учесть шум в системе. Это лучший способ продолжить или есть лучшие методы?
Мы используем LabVIEW, и я проверил форумы LAVA, и есть много интересных примеров. Это часть нашего тестового программного обеспечения, и мы стараемся избегать использования слишком большого количества нестандартных библиотек VI, поэтому я надеялся получить отзывы о задействованных процессах / алгоритмах, а не о конкретном коде.
Вы можете попробовать усреднение сигнала, т.е. для каждой точки усредните значение с окружающими 3 или более точками. Если шумовые помехи огромны, то даже это может не помочь.
Я понимаю, что это не зависело от языка, но, предполагая, что вы используете LabView, есть много готовых VI для обработки сигналов, которые поставляются с LabView, которые вы можете использовать для сглаживания и уменьшения шума. На форумах NI являются отличным местом , чтобы получить более специализированную помощь на такого рода вещи.
Вы можете применить стандартное отклонение к вашей логике и заметить пики свыше x%.
Я думаю, что вы хотите взаимно коррелировать ваш сигнал с ожидаемым, образцовым сигналом. Но прошло много времени с тех пор, как я изучал обработку сигналов, и даже тогда я не обращал особого внимания.
Я не очень много знаю об инструментах, так что это может быть совершенно непрактично, но опять же, это может быть полезным в другом направлении. Если вы знаете, как показания могут потерпеть неудачу, и существует определенный интервал между пиками при таких сбоях, почему бы не сделать градиентный спуск в каждом интервале. Если спуск вернет вас в область, которую вы искали ранее, вы можете отказаться от нее. В зависимости от формы выбранной поверхности это также может помочь вам находить пики быстрее, чем поиск.
Эта проблема была изучена в некоторых деталях.
Существует множество очень современных реализаций в TSpectrum*
классах ROOT (инструмент анализа ядерной физики / физики элементарных частиц). Код работает в одно- и трехмерных данных.
Доступен исходный код ROOT, так что вы можете получить эту реализацию, если хотите.
Из документации класса TSpectrum :
Алгоритмы, используемые в этом классе, были опубликованы в следующих ссылках:
[1] М. Морхак и др.: Методы устранения фона для многомерного совпадения гамма-спектров. Ядерные приборы и методы в физических исследованиях A 401 (1997) 113-132.
[2] М. Морхак и др .: Эффективная одно- и двумерная деконволюция золота и ее применение к разложению гамма-спектров. Ядерные приборы и методы в физических исследованиях A 401 (1997) 385-408.
[3] М. Морхак и др .: Идентификация пиков в многомерном спектре гамма-излучения совпадений. Ядерные приборы и методы в физике исследований А 443 (2000), 108-125.
Документы связаны с документацией класса для тех из вас, у кого нет онлайн-подписки NIM.
Короткая версия того, что сделано, заключается в том, что гистограмма сглаживается для устранения шума, а затем локальные максимумы обнаруживаются с помощью грубой силы в сглаженной гистограмме.
Этот метод в основном из книги Дэвида Марра "Видение"
Gaussian Blur ваш сигнал с ожидаемой шириной ваших пиков. это избавляет от всплесков шума, и ваши данные фазы не повреждены.
Тогда обнаружение края (LOG сделает)
Тогда ваши края были краями элементов (например, вершины). ищите пики между краями, сортируйте пики по размеру, и все готово.
Я использовал варианты этого, и они работают очень хорошо.
Есть ли качественная разница между желаемым пиком и нежелательным вторым пиком? Если оба пика являются «острыми» - т.е. короткими по времени - при просмотре сигнала в частотной области (выполняя БПФ) вы получите энергию в большинстве полос. Но если в «хорошем» пике надежно присутствует энергия на частотах, которых нет в «плохом» пике, или наоборот, вы можете автоматически их дифференцировать таким образом.
Существует множество классических методов обнаружения пиков, любой из которых может сработать. Вы должны увидеть, что, в частности, ограничивает качество ваших данных. Вот основные описания:
Между любыми двумя точками в ваших данных,
(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 выше.
Я делал то, что вы делаете раньше: находя пики в данных последовательности ДНК, находя пики в производных, оцененных по измеренным кривым, и находя пики в гистограммах.
Я призываю вас внимательно следить за правильным базовым уровнем. Фильтрация Винера или другая фильтрация или простой анализ гистограммы часто являются простым способом определения исходного уровня при наличии шума.
Наконец, если ваши данные, как правило, зашумлены и вы получаете данные с карты как односторонний вывод без ссылок (или даже со ссылками, просто без разницы), и если вы усредняете множество наблюдений в каждой точке данных, попробуйте отсортировать их наблюдения и отбрасывание первого и последнего квартиля и усреднение того, что осталось. Есть множество таких тактик исключения, которые могут быть действительно полезными.
Я хотел бы внести в этот поток алгоритм, который я разработал сам :
Он основан на принципе дисперсии : если новая точка данных представляет собой заданное число стандартных отклонений от некоторого скользящего среднего, алгоритм подает сигнал (также называемый 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
. Например, a 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