Поиск частоты

[direct]

Пришлось задачку решить - определить частоту цифрового сигнала (сигнал периодический одной частоты  от 0.0001 Гц до 1 Гц, пишется в файл с дискретом 10 раз в секунду). В принципе, если "тупо" найти спектр в указанном диапазоне, то частота определяется, но это много времени занимает. У кого есть соображения, как частоту быстрее можно найти (или сузить диапазон поиска частот)? И как быть, если не полный период сигнала записан? Спасибо за ответ!

[Alf]

Если под "цифровым сигналом" понимаются импульсы с двумя уровнями (0-1), то я бы не стал связываться с определением спектра (во-первых, спектральное разложение - не столь простая вещь, если нет готовых инструментов, во-вторых, спектр импульсного сигнала с крутыми фронтами содержит слишком много гармоник, которые нам к тому же совершенно не нужны). Достаточно определить период, измеряя интервал между переходами 0-1 и 1-0, а потом вычислить частоту как обратную периоду величину (тем более что именно так обычно и работают частотомеры при измерении очень низких частот).

Относительно второго вопроса (неполный период) уверен, что для цифрового сигнала методики не существует - невозможно, наблюдая постоянный уровень 0 или 1, предсказать время переключения заранее. Если же сигнал аналоговый, например, синусоида стабильной амплитуды, то не составит труда экстраполировать его до точки пересечения с осью Х, а там опять же определить частоту через период.

[direct]

Под "цифровым" сигналом понимается сигнал, полученный с АЦП (уровень кода от 0 до 65535). Я видимо, неверно выразился: это аналоговый оцифрованный сигнал! А про аналоговый сигнал известно, что он периодический (частота от 1 до 0.0001 Гц). Если искать спектр с таким шагом (0.0001 Гц), это вызывает большие накладные расходы. В интервал записи может попасть как несколько периодов, так и менее одного периода! Я более конкретно спрошу: что в такой ситуации выгоднее - оптимизировать расчет спектра или как-то сузить диапазон частот сигнала (как ) (а далее тупо считать спектр)  - вот в чем вопрос!

[Sla]

Ни фига себе диапазон частот у аналогового сигнала!
Это ж сколько точек?
С другой стороны, если данные записываются в файл, то я так понимаю, что о бытродействии преобразований речи и быть не должно. Все- таки не real-time.

[Alf]

direct, а вообще есть ли необходимость в расчете спектра? Если сигнал периодический, достаточно отследить переход через 0 (ну или через 32768, если он смещен на середину диапазона). Этим определится период первой гармоники, из него вычисляется частота. Остальные гармоники, я так понял, и не требуются?
Считать спектр сигнала на участке меньше одного периода в общем случае, по моему убеждению, затея безнадежная. Если амплитуда и форма сигнала известны, тогда имеет смысл воспользоваться интерполяцией, попытавшись совместить известный кусок с рассчитанным полным периодом, а потом определить период интерполированной функции тем же способом, что и выше.

[Sla]

Укажи частоту оцифровки сигнала, гладкость сигнала. Преобразование Фурье позволяет вычислить спектр имея по крайней мере 3 точки.

[Alf]

И сколько же гармоник можно определить по 3 точкам? Про количество периодов сигнала уже и спрашивать неловко...

[Sla]

неловко -не спрашивай
Если известна частота дискретизации - достатосчно 3! точек.
Количество гармоник - а сколько надо?
Надо только вспомнить преобразование Фурье.

[Alf]

...Надо только вспомнить преобразование Фурье.

Золотые слова.

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

Если же 3 точки приходится на несколько периодов сигнала, то старина Котельников должен в гробу перевернуться от таких методов обработки сигналов (да и Найквисту не поздоровится). Ты всерьез способен отличить синусоиду от меандра или там пилы по 3-м точкам?

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

[RXL]

Sla, см. первый пост:

пишется в файл с дискретом 10 раз в секунду

direct, расскажи детальнее: происхождение и форма сигнала, предпологаемый размер буфера, цель задачи.

[Sla]

пишется в файл 10 раз в секунду, - не есть частота оцифровки

Построить спеткр по 3 точкам - построим
с частотой сигнала - проблемка, согласен. и я говорил про спектр, а не про частоту.
Сигнал периодический!

[Alf]

...с частотой сигнала - проблемка, согласен. и я говорил про спектр, а не про частоту.

Так, это дело все интереснее становится... Как говорится, чем дальше в лес, тем толще партизаны.

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

То есть определить частоту проблематично, а получить спектр - запросто? Неувязочки...

Сигнал периодический!

И что, по трем отсчетам реально определить, когда они повторятся снова? Тут даже Нострадамус руки опустил бы.

[Sla]

послушай Зин! не трогай Нострадамуса.
частотные составляющие сигнала, но не частоту сигнала. Есть разница?

[Alf]

...частотные составляющие сигнала, но не частоту сигнала. Есть разница?

Нет разницы. Потому как в разложении периодического сигнала в ряд Фурье первая составляющая и есть частота повторения сигнала. Если ты действительно знаешь частотные составляющие сигнала, то автоматически получаешь и частоту самого сигнала.

Ладно, слишком много слов получается. Я тут, как и обещал, набросал не слишком сложный сигнал в Excel'е как сумму трех гармоник и взял три последовательных отсчета. Вот их значения:

Код:

0.264095659
0.255317054
0.299569688

Эксперимент достаточно чистый, поскольку частота выборки многократно превышает частоту третьей гармоники.

Принимая частоту выборки за единицу, жду данных о спектре сигнала.

[direct]

Sla, см. первый пост:

пишется в файл с дискретом 10 раз в секунду

direct, расскажи детальнее: происхождение и форма сигнала, предпологаемый размер буфера, цель задачи.

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

[direct]

Я вот что хотел сказать: это сигнал обладает шумами, поэтому идею поиска перехода через ноль я пока "похоронил" (не знаю, как ее реализовать). Идею "по 3 точкам" хороню сразу (обратили внимание на диапазон частот?), если частота 1 Гц будет, то по трем точкам не определишь частоты! Что касается поиска спектра, то я по максимуму спектра нахожу основную частоту, но, блин. медленно работает (слишком широкий диапазон). У кого есть идеи по сужению диапазона частот этого сигнала?

[Lex]

direct, есл у тебя буфер 256 точек, то неужели дискретное преобразование для 256 точек медлено работает?

Кстати несколько непоняток. Как с таким буфером ты хочешь ловить частоту в 0.0001Гц ? период займет 100000 точек.

[Dimka]

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

[direct]

direct, есл у тебя буфер 256 точек, то неужели дискретное преобразование для 256 точек медлено работает?

Кстати несколько непоняток. Как с таким буфером ты хочешь ловить частоту в 0.0001Гц ? период займет 100000 точек.

Только одним способом - увеличивать интервал между точками! Будет не 10 точек в секунду, а, скажем, 1 точка за 10 сек:). Если период не получится, то придется по неполному периоду частоту искать!
Но 256 точек  - это для ЦПС, где будет работать алгоритм в реальном времени. А в файло я пока могу и побольше точек записать (главное, чтобы в принципе частота определялась!)

[direct]

Расчет спектра можно ускорить, применив ШПФ (шустрое преобразование фурье )! А если в буфере данных меньше периода, то расчет спектра дает лишь оценку основной частоты (искажения большие!). Но ведь я видел, как цифровые осциллографы выводят частоту и амплитуду сигналов (даже с неполным периодом), какие алгоритмы там применяются?

[Lex]

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

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

И точность этих измерений +- лапоть.

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

[direct]

По ходу дела, так и придется делать. Сначала фильтром урезать полосу частот, далее поиск спектра. Благодарю за обсуждение вопроса!

[Alf]

Я вот что хотел сказать: это сигнал обладает шумами, поэтому идею поиска перехода через ноль я пока "похоронил" (не знаю, как ее реализовать).

А не проще подавить шумы фильтром и искать-таки переход через 0 сглаженной кривой? Зачем возиться со спектром, если все равно информация о гармониках, кроме основной, не нужна?

Идею "по 3 точкам" хороню сразу (обратили внимание на диапазон частот?), если частота 1 Гц будет, то по трем точкам не определишь частоты!

Что похоронил - это правильно, уж больно порочная идея. Только при чем тут диапазон частот? А если частота 1 ГГц будет, то определишь? Тут в самом принципе дело.

[RXL]

Хлчу еще напомнить неплохую книгу: http://www.dspguide.com

[direct]

А не проще подавить шумы фильтром и искать-таки переход через 0 сглаженной кривой? Зачем возиться со спектром, если все равно информация о гармониках, кроме основной, не нужна?

Что похоронил - это правильно, уж больно порочная идея. Только при чем тут диапазон частот? А если частота 1 ГГц будет, то определишь? Тут в самом принципе дело.

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

[RXL]

direct, можно считать спектр ступенчато:
1) выборка 10Гц на 256 точек: 5-0.039Гц
2) из предыдущей выборки выбрать с некоторым шагом несколько семплов (напр. 2 точки) и делать это каждый раз при выполнении 1-го шага. Когда наберется полный буфер (скажем, те же 256 точек), то выполнить fft и для него (получается 10Гц/128 : 0.039-0.0003Гц).

Т.е., каждые 25.6 секунды у тебя вычисляется первый спектр, а каждые 3276.8 - второй. Кстати, перед выборкой семплов для второго шага полезно будет прогнать его через фильтр, да бы обрезать частоты выше половины частоты семплирования.

[Alf]

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

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

Обработка будет выглядеть примерно так.

1. Набить буфер размером n байт входными данными.

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

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

При таком подходе даже не понадобится буфер на 256 точек, он здесь явно избыточен. Вполне достаточно будет буфера размером 2*k+1, где, как я уже говорил, k<10.

Буфер в данном случае организован в виде очереди FIFO, т.е. при считывании очередного слова из входного файла все данные сдвигаются, выталкивая самое старое значение.

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

[Lex]

Согласен с методом Alf'а. Тут вся тонкость будет найти правильное пересечени ноля. Скорее всего придется немного поиграться с разными методами фильтрации и усреднения сигнала. А так наиболее правильный метод из предложенных для детектирования такого диапазона частот.

[direct]

Очень привлекательный метод. Согласен, с фильтрами придется повозиться (ведь надо монотонный сигнал получить). Ну а если в буфер попадет, скажем, только одно пересечение с нулем (меньше, чем полупериод), что тогда делать?

[Alf]

direct, я, наверное, не до конца прояснил свою идею.

Представь себе буфер небольшого размера (порядка десятка отсчетов) в виде очереди FIFO. Он используется как скользящее окно для просмотра измеренных значений. Содержимое буфера используется для фильтрации шумов, после чего по нему определяется единственное значение функции.

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

Когда при очередном отсчете найден переход функции через 0:

1. Отмечаешь текущее время отсчета.
2. Находишь разницу между временем текущего перехода через 0 и предыдущего. Как догадываешься, это и есть период сигнала.
3. Находишь частоту как величину, обратную периоду.
4. Сохраняешь найденное время переходя через 0 как предыдущее и ищешь следующий переход.

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

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

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