Добрый день.

Есть некий медленно меняющийся сигнал, отягощенный влиянием гармонической составляющей с переменным периодом (см картинку).
Нажмите для просмотра прикрепленного файла

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

В точке локального минимума (и максимума) производная непрерывной функции равна нулю.

Конечно! Спасибо!

в матлабе вот так можно

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

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

Если у Вас функция - суть сумма двух-трех синусов и косинусов, с заранее неизвестными фазами, амплитудами и периодами, и шума нет, или почти нет, можно сделать еще элегантнее:

загоняете оцифрованную функцию в вектор v_1, сдвигаете этот вектор на 1 вниз, сохраняя в v_2, и так далее столько раз, сколько у Вас синусов.
Подрезаете сверху и снизу вектора, чтобы одинаковое количество элементов было, и строите матрицу. Находите у нее правый сингулярный вектор, соответствующий минимальному сингулярному числу. Записываете элементы этого вектора в виде коэффициентов полинома и находите корни. Они все - Ваши периоды деленные на шаг дискретизации. Зная периоды, наименьшими квадратами находите амплитуды и фазы. Дальше - аналитически ищите все корни и решение у Вас на ладони.

Выглядит страшно, но зная лапак + мапл или математику, можно за несколько часов сделать и запрограммировать. В МК, на который лапак не компилится, реализовать будет ужасно сложно, думаю, что практически невозможно, в этом случае, остается только первый, не элегантный, метод.

Для начала я бы избавился от низкочастотной составляющей простейшим ФВЧ :-)

Дамы и господа,

виноват, не упомянул, что работаю с дискретным массивом данных. В самом начале там очень мало точек на период сигнала, т.е. прикладываться там Ньютоном не к чему. Пробовал дифференцировать данные еще до совета Татьяны, но опять-таки из-за малого кол-ва точек на период искать прохождение через 0 приходится "руками". И почему-то я думал, что есть какой-то красивый способ из области мат. физики чтоб получить все и сразу. Спасибо, что развеяли мои надежды.

iiv,
попробую, если будет возможность провести измерения с бОльшим разрешением.

serjj,
я понял, поиск по разнице в знаке. Спасибо.

Santik,
мне как раз и требуется получить красивую кривую этого НЧ сигнала, вычислив локальное среднее на базе амплитуды гармонического сигнала (да, именно так, всевозможные сглаживания применять нельзя).

Если я правильно понял, Вас интересует зависимость fвч(t). Рекомендую непрерывный вейвлет- анализ.
Смотреть здесь.
Если НЧ нужна составляющая - то ФНЧ можно применить.
Непонятно о каких сглаживаниях Вы говорите....

А в чём проблема применять фильтрацию? Медианное сглаживание с окном большим периода гармонического сигнала + скользящее среднее или просто узкополосный ФНЧ, который отрежет гармоническую составляющую? У вас есть априорные данные о соотношении частот НЧ и ВЧ составляющих в сигнале? Если да, то не виже никаких проблем с фильтрацией, т.к. сигналы разесены по частоте. Если же априорных данных нет, то можно посмотреть в сторону BSS (blind signal separation).

Что означает "переменный период"? Такого не бывает. Если Вы знаете (априорно), что мешающая Вам компонента описывается некоторой формулой с параметрами, то и нужно ее подогнать так, чтобы производная меняла знак один раз - на пике медленной составляющей функции.

Бывает. Пример - ЛЧМ (линейная частотная модуляция)
Модель:
do i=1,NN
t=(i-1)/Fd
FXSw(i)=sin(Pi2*(F1*t-F2*t*t/(2*T0)))+10*(1-cos(Pi2*1.*t))
end do
Картинка

Нет, не бывает... Математика - это такая конвенциональная игра. Конвенция это запрещает.

Запутал я всех капитально. Извините.


Все намного проще, приведенные данные - спектральный отклик фотосенсора, имеющего буферный слой SiO2 на поверхности. Из-за слоя SiO2 возникает паразитная интерференция, которая зависит от длины волны.

Решил задачу нахрапистым способом, как советовал serjj, сохраняя индекс экстремума если есть разница в знаке между соседними точками. Только я добавил в этот цикл условие, что ширина окна между экстремумами может варьироваться линейно в зав-ти от длины волны. Задача чисто инженерная. Вот результат:
исходные данные (красная линия), среднее значение (черная линия и точка) и массив экстремумов (синие звездочки).
Нажмите для просмотра прикрепленного файла

Еще раз прошу извинить что отнял ваше время. Спасибо за советы.

Ниччо не понял...
Фильтрация
Нажмите для просмотра прикрепленного файла
Модель + Фильтр Чебышева 2 рода 6-го порядка
А если всё-таки нужно смотреть как меняется частота ВЧ составляющей - это уже надо непрерывный вейвлет-анализ :-)

"НЧ-составляющая" на моей картинке - типичный спектральный отклик матрицы кремниевого фотосенсора. Исследуемый образец имеет довольно толстый слой SiO2 поверх Si. Границы воздух-SiO2 и SiO2-Si образует таким образом интерферометр Фабри-Перо, который в зависимости от длины волны либо усиливает либо ослабляет световой сигнал (интерференция), поглощаемый кремниевым фотодиодом.


звучит серьезно, если шеф одобрит, буду изучать вейвлет-анализ )

Да я не понял математических "танцев с бубном"! (Физику явления я не понял тем более :-)
А ЦОС надо изучать в независимости от одобрения Шефа :-)

Надо было использовать волновое число вместо длины волны... Было бы лучше... В предположении малой дисперсии...

если речь о сглаживании, то тоже подпишусь под Чебышевским фильтром, вейвлеты могут дать большую нестабильность. Если надо и низкую, и высокую частоту вытащить, возможно проще сделать Фурье, дальше для низких грубо так занулить все высокочастотные, а для высоких, - занулить низкочастотные гармоники, и потом сделать обратное Фурье.

Надо было:
1. Отфильтровать сигнал фильтром высокой частоты (ФВЧ), т.е. получить сигнал ВЧ-помехи.
2. Сделать взаимную корреляцию исходного сигнала и ВЧ- помехи (определение задержки)
3. С учётом задержки вычесть из исходного сигнала сигнал помехи.

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

да такую задачку можно решать не одним способом, просто перед реализацией надо задать себе вот такие вопросы:
1 - какие сигналы могут присутствовать в смеси? --имеются в виду и всевозможные шумы тоже
2 - при каком минимальном соотношении сигнал-шум и сигнал-помеха алгоритм должен корректно отрабатывать?
3 - какова допустимая ошибка эстиматора?
4 - на чём я буду этот алгоритм реализовывать? -- если это матлаб, то тут все средства хороши, но если это дохлый контроллер без fpu, то уже стоит задуматься
5 - нужен ли real time?
6 - если нужен real time, какой latency я могу себе позволить (это к вопросам о фильтрах)

Может быть я что-то забыл, но ИМХО это минимум, на который стоит найти ответы перед тем как сесть и запилить


Ну просто мы сигнальщики и привыкли работать с сигналами, не задаваясь вопросами откуда там взялась помеха или шум, просто берешь и обрабатываешь) Автор нарисовал картинку и сросил как из неё другую картинку получил, вот ему и советуют.

Автор снимает спектр действия фотодетектора в спектрометре... Медленно и печально... Снаружи приклеился паразитный эталон Фабри-Перо.

То есть фильтр может быть некаузальным. Это добавляет ещё с десяток возможных решений этой задачи :-)
Предложенный автором способ хорош только в случае если шкала наблюдения по длине волны неравномерная.

Шкала неравномерная? Нет, у него интерференционные максимумы и минимумы неравномерно лежат на его равномерной шкале. Они бы лежали равномерно на моей равномерной шкале волнового числа.

А почему не прокатил простой КИХ НЧ фильтр? по идее хоть прямоугольное хоть гаусовое окошки дадут вам нормальное среднее. ширина окна мне предполагается достаточна 2 периода ВЧ составляющей.

1 - какую полосу фильтра брать?
2 - какой порядок (и как следствие latency) из этого порядока последуют?
Имхо при высоком сигнал-шуме, это из пушки по воробьям.

Ну тогда я не вижу вообще никаких преимуществ в предлагаемом методе.
Насколько я понял метод заключается в проведении прямой от максимума до минимума и нахождения середины этого отрезка.
Это же сколько полезных отсчётов отбрасывается! А ведь максимум и минимум измеряется с определённой погрешностью, так что общая ошибка может быть существенной.
Мне не понятно, почему автор так упорно настроен против фильтрации. Да, действительно, при фильтрации в реальном времени получается сдвиг сигнала во времени - но это не наш случай - выяснилось, что фильтр может быть некаузальным! Так что создать фильтр без сдвига по шкале -это элементарная задача. (Ну может быть и не совсем элементарная, но гораздо проще, чем фильтрация в реальном времени).
Если максимумы/минимумы будут равномерно на шкале располагаться - это конечно лучше (узкополосную помеху проще вырезать), но не так уж и принципиально - полезный сигнал имеет частоту гораздо ниже, чем сигнал "помехи".
Преобразование Фурье - вот лучший способ для таких задач!

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

Ну я и занимаюсь ЦОС, просто вейвлеты - это уже для труъ специалистов

По поводу модели, вот картинка:
Нажмите для просмотра прикрепленного файла
слой SiO2 образует интерферометр (на самом деле, не только SiO2, но и подлегированная часть).

Хорошо что в теме есть Tanya

Ну не так и печально все. Очень даже весело


Совершенно верно! Можно было перейти к волновым числам, определение экстремумов в этом случае легче.

Господа (и дамы), ну вы тут разошлись. Задача решалась на ПК скриптом на питоне. Никаких МК для вычислений, слава аллаху, не предполагается.
Всем спасибо за дискуссию, много интересного узнал
Пожалуй, надо закрывать тему, а то скоро предложат облачные вычисления, IoT и Raspberry PI )

А Фурье здесь не причём!
Всегда "для существенного разрешения по частоте требуется довольно большое количество отсчётов"!
Это, если так можно сказать, закон природы!
Методов определения частоты сигнала (с большой точностью) длительностью 1 период - не существует, сколько отсчётов ни возьми!

Пожалуй да, но нужно кое-что прояснить, раз такое дело

Вы неправы, вот вам пример:
Нажмите для просмотра прикрепленного файла
и
Нажмите для просмотра прикрепленного файла
Господа Марпл, Хайкин и пр. в своё время много сил положили на спектральный анализ коротких выборок и/или разделения близко расположенных по частоте сигналов, а вы говорите, что не существует. Разрешающая способность преобразования Фурье по частоте не является непреодолимым и фундаментальным, она определяет разрешение только для этого метода и производных от него (периодограммы, коррелограммы), а не всего спектрального анализа. Другое дело, что у каждого метода есть свои ограничения и подводные камни, и Фурье самый универсальный и "неприхотливый" по условиям применения и ресурсам. Но когда есть задача по малому количеству отсчётов определить частоты гармоник, входящих в смесь, удтверждать, что задача физически нереализуемая только на основании теории Фурье - неправильно.

вам потребуется ких с симметричным окном. и latency такого фильтра однозначна D = W/2 - где D- задержка[отсчеты], W - ширина окна - т.е. результат находится в центре окна. по сравнению с Вашим методом - полагаю что ширины окна в 2 периода ВЧ составляющей хватит. форма окна конечно важна - от нее зависит амплитуда пролезающих на выход ВЧ компонент. так или иначе - даже прямоугольное окно, имхо, будет не намного хуже того что Вы сделали.
в отличие от Вашего фильтра КИХ примитивы хорошо отлажены, оптимизированы, и за счет отсутсвия ветвлений имеют довольно высокую скорость обработки. а в некоторых ДСП встречаются даже отдельные ядра именно для обработки КИХ.
Еще один минус Вашего фильтра - чувствительность к ВЧ помехам - так как вы семплируете их лишь в точках экстремумов, то у Вас нет возможности их убрать, они прямо лезут в ваш результат. и наконец, реальный экстремум может находиться между отсчетами, и если ваши отсчеты от него далеко то вы получаете дополнительную погрешность (апертурная ее чтоли назвать?)
вобщем если период ВЧ сигнала у вас укладывается в 20-30 семплов, то КИХ - стоит рассматривать.

Я эти вопросы привёл не потому что хочу получить на них ответы, а чтобы проиллюстрировать, что нужно держать в голове приступая к задаче выделения сигнала из смеси в случае одномерной задачи.
То, что вы написали - прописная истина, я же имел в виду цепь рассуждений - есть смесь из двух или более сигналов, если мы хотим разделить их фильтром, то как выбрать полосу фильтра (знаем ли мы априори соотношение частот или нет) - далее, если с полосой определились, мы можем оценить требуемый порядок для выбранной структуры и подавления за полосой - порядок в свою очередь даст нам latency, по формуле, которую вы представили.
Так вот, может оказаться, что а) мы ничего не знаем о соотношении полос, б) полоса слишком узкая, что даст большую latency, в) задача не real time (как у ТС), г) шумами можно пренебречь, поэтому возможно использовать простые численные методы и т.д. В этих случаях КИХ неприменим или избыточен.
Только после проработки данных вопросов можно приступить к самому КИХ фильтру (делов на 5 минут). Нужно задачу с головы решать
А у ТС вообще всё просто оказалось.

я делал подобный фильтр как раз для реалтайм приложения. его скорость меня огорчила. правда у меня стояла более жесткая задача - интерполировать позиции экстремумов и нулей в сетке Fs*8. отдельных траблов добавили все же вч шумы, которые у меня в сигнале таки присутствовали. даже 1дискрет АЦП вызывал плохо-осознаваемые артефакты, практически неотслеживаемые.
насколько я понял, Вы собрались сделать некий фильтр-убийцу способный работать на любой паре смешиваемых частот. полагаю для случая 2х частот - Вы правы. Ваш метод сработает и будет быстр. но уже на 3х частотах - он неработоспособен. Ваш фильтр сносен для очень узкой задачи.

В цифровом мире это звучит иначе:
При прохождении через точку экстремума производная меняет знак.



Что значит руками? Если производная поменяла знак -- значит экстремум между последних двух точек. Вам и скрипт написали. Если недосаточно точек для этого, то никакая математика их не заменит.

"руками" и есть сканировать скриптом на предмет прохождения через 0 анализируя знак пары точек производной. Это был ответ к сообщению Tanya.


производная в любом мире ведет себя одинаково. То что Вы описываете касается случая, когда экстремум в наборе данных имеет недостаточно точек (случай коротких длин волн в моих спектрах).

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

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

Это точно улучшит результаты.

Марпл , Хайкин... А всё равно кошерных методов "и/или разделения близко расположенных по частоте сигналов" - не существует!
Предлагаю рассмотреть элементарную задачу - протонный магнитометр. За 0.5 - 1 сек предлагается определить частоту прецессии (2-3 кГц) с точностью 0.001 Гц

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

Ну есть в оптике критерий Рэлея.
Разрешение преобразования Фурье определяется длиной исследуемой временной последовательности.
Длиннее последовательность -- больше точек (они чаще стоят, а значит рассояние между ними меньше) на выходе. Ведь длина спектра не изменится.



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

А я не согласна, что Вы дополнили.

Как вам будет угодно, сударыня.

Вот главное слово здесь - " длиной исследуемой временной последовательности".
А не количеством точек ... Я АЦП могу взять 1 МГц для оцифровки 3кГц сигнала - точек много - эффекта 0.

А может, Вы просто не умеете добиваться положительного (> 0) эффекта?

Меж тем, при тех условиях, что были озвучены ранее, добиться "положительного эффекта" можно уже просто измерив "2-3 кГц сигнал" всего в трех точках..

Пусть, к примеру, нам известно, что значение "2-3 кГц сигнала" в точках t1, t2 и t3 равно S1, S2 и, S3, соответственно..

Причем, мы выбираем моменты измерения сигнала таким образом, чтобы:

t3 - t2 = t2 - t1,

и при этом:

f*(t3 - t1) < 1, и

S2 ≠ 0.

Пусть сигнал представляет собой синусоиду с неизвестной частотой "ω", амплитудой "A" и фазой "φ":

S(t) = A*sin[ω*t + φ].

Тогда находим:

S1 = A*sin[ω*t1 + φ],
S2 = A*sin[ω*t2 + φ],
S3 = A*sin[ω*t3 + φ].

Складывая первое уравнение с третьим, находим:

S3 + S1 = A*{sin[ω*t3 + φ] + sin[ω*t1 + φ]}.

Тогда:

S3 + S1 = 2*A*sin[ω*(t3 + t1)/2 + φ]*cos[ω*(t3 - t1)/2] = 2*A*sin[ω*t2 + φ]*cos[ω*(t3 - t1)/2] = 2*S2*cos[ω*(t3 - t1)/2].

Или:

cos[ω*(t3 - t1)/2] = (S3 + S1)/(2*S2).

После чего находим частоту:

f = [1/[pi*(t3 - t1)]]*arccos[(S3 + S1)/(2*S2)].

Как-то так..

Опять же, что из себя представляет сигнал этого протонного магнитометра? Если такой сигнал можно рассматривать как комплексную экспоненту с шумом, то все ваши фантазии будут строго ограничены пределом Крамера-Рао, который зависит от длины выборки и соотношения С/Ш и, используя различные методы интерполяции можно достаточно близко приблизиться к этой границе (тема поднималась неоднократно). Если же сигнал имеет более сложную структуру, то однозначного ответа нет, надо смотреть по ситуации.

Решение красивое, но что будет, если а) добавить шум, б) добавить интерференцию
Можно привести пример. Из пространственного спектроанализа, но сути это не меняет. В АФАР можно определить направление на цель всего по 2-м сенсорам, оценив задержку сигнала, полученного с одного сенсора относительно другого. Однако эта схема работает с относительно сильными и узкополосными сигналами, не содержащими интерференцию. Если же это условие не выполняется (что происходит довольно часто), то придётся а) наращивать количество сенсоров (т.е. количество сэмплов для анализа) и б) усложнять алгоритм (в отличие от временного спектроанализа, где мы в основном можем увеличить длину выборки, в пространственном размер решётки ограничен, т.е. Фурье с его малым разрешением неприемлим).


Давайте рассмотрим аналогичную задачу (но сделаем её ещё сложнее):
Нажмите для просмотра прикрепленного файла
символьная частота = 100 кГц, ОСШ = 60 дБ, частота сигнала = 80,348 Гц, количество символов для анализа = 1024 + 128 (как видно из картинки - менее одного периода), погрешность ~ 0,001 Гц.

Ну мои сексуальные фантазии так далеко не заходят... :-)

печалька , не богатая значит у Вас фантазия, трудно подобрать более простой сигнал, чем exp(jwt)

Ну тогда порадуй старика новым оригиальным решением :-)

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