Перейти к содержанию
    

Найти экстремумы функции с гармонической составляющей

Пожалуй, надо закрывать тему, а то скоро предложат облачные вычисления, IoT и Raspberry PI

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

А Фурье здесь не причём!

Всегда "для существенного разрешения по частоте требуется довольно большое количество отсчётов"!

Это, если так можно сказать, закон природы!

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

Вы неправы, вот вам пример:

post-81866-1438241292_thumb.png

и

post-81866-1438241300_thumb.png

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

1 - какую полосу фильтра брать?

2 - какой порядок (и как следствие latency) из этого порядока последуют?

Имхо при высоком сигнал-шуме, это из пушки по воробьям.

вам потребуется ких с симметричным окном. и latency такого фильтра однозначна D = W/2 - где D- задержка[отсчеты], W - ширина окна - т.е. результат находится в центре окна. по сравнению с Вашим методом - полагаю что ширины окна в 2 периода ВЧ составляющей хватит. форма окна конечно важна - от нее зависит амплитуда пролезающих на выход ВЧ компонент. так или иначе - даже прямоугольное окно, имхо, будет не намного хуже того что Вы сделали.

в отличие от Вашего фильтра КИХ примитивы хорошо отлажены, оптимизированы, и за счет отсутсвия ветвлений имеют довольно высокую скорость обработки. а в некоторых ДСП встречаются даже отдельные ядра именно для обработки КИХ.

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

вобщем если период ВЧ сигнала у вас укладывается в 20-30 семплов, то КИХ - стоит рассматривать.

Изменено пользователем AlexRayne

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

вам потребуется ких с симметричным окном. и latency такого фильтра однозначна D = W/2 - где D- задержка[отсчеты], W - ширина окна - т.е. результат находится в центре окна. по сравнению с Вашим методом - полагаю что ширины окна в 2 периода ВЧ составляющей хватит. форма окна конечно важна - от нее зависит амплитуда пролезающих на выход ВЧ компонент. так или иначе - даже прямоугольное окно, имхо, будет не намного хуже того что Вы сделали.

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

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

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

Только после проработки данных вопросов можно приступить к самому КИХ фильтру (делов на 5 минут). Нужно задачу с головы решать :rolleyes:

А у ТС вообще всё просто оказалось.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

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

Только после проработки данных вопросов можно приступить к самому КИХ фильтру (делов на 5 минут). Нужно задачу с головы решать :rolleyes:

А у ТС вообще всё просто оказалось.

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

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

Изменено пользователем AlexRayne

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

В цифровом мире это звучит иначе:

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

 

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

 

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

 

iiv,

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

 

serjj,

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

 

Santik,

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

 

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

 

В цифровом мире это звучит иначе:

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

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

 

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

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

 

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

Изменено пользователем Tarbal

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Марпл , Хайкин... А всё равно кошерных методов "и/или разделения близко расположенных по частоте сигналов" - не существует!

Предлагаю рассмотреть элементарную задачу - протонный магнитометр. За 0.5 - 1 сек предлагается определить частоту прецессии (2-3 кГц) с точностью 0.001 Гц

Изменено пользователем Santik

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Я вот не зря написала про непрерывную функцию. Вижу, что не всем понятно...

Поясню. Предполагалось, что в данном случае нужно выбрать несколько точек (5 -7...)

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

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Марпл , Хайкин... А всё равно кошерных методов "и/или разделения близко расположенных по частоте сигналов" - не существует!

 

Ну есть в оптике критерий Рэлея.

Разрешение преобразования Фурье определяется длиной исследуемой временной последовательности.

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

 

Я вот не зря написала про непрерывную функцию. Вижу, что не всем понятно...

Поясню. Предполагалось, что в данном случае нужно выбрать несколько точек (5 -7...)

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

Это в первом приближении...

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

 

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

Изменено пользователем Tarbal

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Я вас не оспаривал, а просто дополнил. Полностью согласен с тем, что вы написали.

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Разрешение преобразования Фурье определяется длиной исследуемой временной последовательности.

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

Вот главное слово здесь - " длиной исследуемой временной последовательности".

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

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Марпл , Хайкин... А всё равно кошерных методов "и/или разделения близко расположенных по частоте сигналов" - не существует!

Предлагаю рассмотреть элементарную задачу - протонный магнитометр. За 0.5 - 1 сек предлагается определить частоту прецессии (2-3 кГц) с точностью 0.001 Гц

Вот главное слово здесь - " длиной исследуемой временной последовательности".

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

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

 

Меж тем, при тех условиях, что были озвучены ранее, добиться "положительного эффекта" можно уже просто измерив "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)].

 

Как-то так..

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Предлагаю рассмотреть элементарную задачу - протонный магнитометр. За 0.5 - 1 сек предлагается определить частоту прецессии (2-3 кГц) с точностью 0.001 Гц

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...