Преобразование Фурье

Сообщение №9181 от 13 ноября 2003 г. 11:25
Тема: Преобразование Фурье


Отклики на это сообщение:

Доброе время суток ВСЕМ!
Вопрос следующий, имеется некое устройство АЦП, которое собирает данные с частотой дискретизации 260.870 кгц, так вот потом со снятым сигналом требуется сделать бпф и построить спектр короче. Ближайшее подходящее значение степени двойки 262144, и следовательно частотный шаг будет не 1 гц, а мне требуется 1 гц, как тут быть?


> Доброе время суток ВСЕМ!
> Вопрос следующий, имеется некое устройство АЦП, которое собирает данные с частотой дискретизации 260.870 кгц, так вот потом со снятым сигналом требуется сделать бпф и построить спектр короче. Ближайшее подходящее значение степени двойки 262144, и следовательно частотный шаг будет не 1 гц, а мне требуется 1 гц, как тут быть?

Дополнить секундный отрезок нулями.
Или посчитать, как есть, а потом интерполировать спектр на нужные точки. Принципиально результат не отличается, хотя в реализации возможны различия (из-за окон и т.п.)


> Доброе время суток ВСЕМ!
> Вопрос следующий, имеется некое устройство АЦП, которое собирает данные с частотой дискретизации 260.870 кгц, так вот потом со снятым сигналом требуется сделать бпф и построить спектр короче. Ближайшее подходящее значение степени двойки 262144, и следовательно частотный шаг будет не 1 гц, а мне требуется 1 гц, как тут быть?

Дополнить секундный отрезок нулями.
Или посчитать, как есть, а потом интерполировать спектр на нужные точки. Принципиально результат не отличается, хотя в реализации возможны различия (из-за окон и т.п.)


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


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

Есть. Например, алгоритм Винограда.


Большое спасибо! А не посоветуете ли вы, где в Интернете можно найти удобоваримое описание этого алгоритма? А то попытка поюзать Google/Rambler толковых результатов не дала... Кстати, если кому понадобится реализовать обычный БПФ, то вот тут http://psi-logic.narod.ru/fft/fft.htm имеется вариант на C.


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

Есть.
1. Хорошо работает БПФ для степеней тройки при "нетрадиционном" представлении комплексных чисел
2. Общая идеология такова. Если число отсчетов - составное, то есть метод (декомпозиция Гуда-Томаса) сведения к БПФ длин, равным различным взаимно простым сомножителям.
3. Слухи о Винограде - сильно раздутая реклама частных результатов. Кстати, эти алгоритмы ориентированы на минимизацию только числа умножений, возможно и за счет непропорционального увеличения числа сложений.
4. Есть ОЧЕНЬ плохие простые длины, для которых попытка В ПРИНЦИПЕ сделать что-то наукообразное приводит к алгоритмам со сложностью, выше тривиальной оценки (квадрат числа отсчетов).
5. Несколько лет назад мы исследовали задачу об оптимальном пополнении нулями до ближайшей "хорошей" длины.
Если есть конкретные значения параметров - напишите, определим оптимальную структуру алгоритма (только до 30.11. я буду в отъезде)
6. Из Интернетовского ширпотреба есть неплохая по подбору алгоритмов инсталлируемая библиотека FFT

http://fftw.org


Конкретная задача состоит в анализе дискретного звукового сигнала. Требуется взять звуковой фрагмент

F(t), t= [x1, x2]

(к сожалению произвольной длины) и разложить его на сумму гармоник:

Fj(A*sin(ω t + φ)

Притом крайне желательно свести число гармоник к минимуму и максимально точно определить ω.

Для этого я сейчас использую "медленный" алгоритм преобразования Фурье. В результате получается спектр. В его первой половине обнаруживается максимум (главная гармоника). Берутся также соседние два отсчета и делается аппроксимация многочленом по трем точкам. Обнаруживается предполагаемый максимум - его частота и амплитуда. Аналогично находится фаза. Затем из спектра ПФ вычитается спектр для найденного обертона и ищется следующий максимум. Не знаю, насколько этот алгоритм оптимален, советы и замечания очень приветствуются. Ясно, что узкое место алгоритма - ПФ. БПФ я запрограммировал, но он работает только на числе отсчетов, равном степени двойки. А если дополнять нулями, то в спектре появляются "лишние" частоты, которые ухо вовсе и не различает. Можно было бы дополнить сигнал так, как будто бы дополнительная часть имеет точно такой же спектр, что и исходная часть. Но для этого прежде надо найти спектр исходной части, а для этого его надо дополнить. Замкнутый круг :) Ускорить обычный ПФ непонятно как. Я даже не могу использовать умножения вместо возведения экспонент в степень, потому что при многократных умножениях накапливается ошибка округления.

В общем, программист просит любых советов от математиков, которые могут им показаться полезными :) Возможно, я отстал от жизни и звук вообще анализируется нынче не спектральным анализом по Фурье?


> Конкретная задача состоит в анализе дискретного звукового сигнала. Требуется взять звуковой фрагмент

> F(t), t= [x1, x2]

> (к сожалению произвольной длины) и разложить его на сумму гармоник:

> Fj(A*sin(ω t + φ)

> Притом крайне желательно свести число гармоник к минимуму и максимально точно определить ω.

> Для этого я сейчас использую "медленный" алгоритм преобразования Фурье. В результате получается спектр. В его первой половине обнаруживается максимум (главная гармоника). Берутся также соседние два отсчета и делается аппроксимация многочленом по трем точкам. Обнаруживается предполагаемый максимум - его частота и амплитуда. Аналогично находится фаза. Затем из спектра ПФ вычитается спектр для найденного обертона и ищется следующий максимум. Не знаю, насколько этот алгоритм оптимален, советы и замечания очень приветствуются. Ясно, что узкое место алгоритма - ПФ. БПФ я запрограммировал, но он работает только на числе отсчетов, равном степени двойки. А если дополнять нулями, то в спектре появляются "лишние" частоты, которые ухо вовсе и не различает. Можно было бы дополнить сигнал так, как будто бы дополнительная часть имеет точно такой же спектр, что и исходная часть. Но для этого прежде надо найти спектр исходной части, а для этого его надо дополнить. Замкнутый круг :) Ускорить обычный ПФ непонятно как. Я даже не могу использовать умножения вместо возведения экспонент в степень, потому что при многократных умножениях накапливается ошибка округления.

> В общем, программист просит любых советов от математиков, которые могут им показаться полезными :) Возможно, я отстал от жизни и звук вообще анализируется нынче не спектральным анализом по Фурье?

Понял задачу и проблемы паразитных частот при дополнении.
Все же сигнал у Вас длины не совсем произвольной.
Чуть пошевелить можно. Или обрезать немного , если пополнять не хочется.
Вы все же посмотрите библиотеку на
fftw.org

Вполне приличная для одномерного ДПФ

А минимизация числа гармоник - проблема не длины, а сигнала. Его корр.свойств.
Тут уж как получится. Где энергия локализуется...

Извините, я сейчас на чемоданах.
Если необходимо - после приезда продолжим.

Можно и подробно о проблемах мне намылить.
Приеду - прочту.


Не могу найти подходящую литературу, а самостоятельно решить не получается. Может кто знает метод решения или посоветует литературу в Инете?

Итак. Был исходный гармонический сигнал на отрезке [0, T]:

f(x)= A cos(2πν + φ)

Этот сигнал подвергли дискретизации, разбив отрезок на N частей. Теперь задача: восстановить исходные амплитуду A, частоту ν=m/T и фазу φ. Если частота случайно угодит на точку дискретизации (m - целое) - то все клево. Делаем дискретное преобразование Фурье, получаем в нем единственный не ноль на месте искомой частоты и приехали.

Иначе произойдет так наз. "растекание" и получится такая феня:

Xk = A/2[exp(jφ) P + exp(-jφ) R]
P = (exp(p) - 1) / (exp(p/N) - 1)
R = (exp(r) - 1) / (exp(r/N) - 1)
p = 2πj(m-k)
r = -2πj(m+k)

Здесь Xk, k= 0..N-1 - результаты дискретного преобразования Фурье. Максимальное Xk будет соответствовать частоте k/T и окажется слева или справа от настоящей частоты m/T.

Я могу найти A, ν и φ тупыми численными методами. А можно ли получить результат аналитически или хотя бы как-то упростить задачу?

26 ноября 2003




> f(N)= A cosfN+BsinfN
Я переформулировал несущественно. Видно, что 1/обычно достаточно знать 3 точки, 2/решений для частоты много.
Алгебраизовать можно, думаю, переходом к тангенсу половинного угла - тут и масштаб (по А,В) уберется. Впрочем, этот путь не пробовал.
Практически, с учетом ошибок, можно минимизировать ошибку приближения - после простого выделения диапазонов. Возможно, целесообразно рассмотреть функцию (от частоты)"ошибка метода наименьших квадратов".
nikvic@newtech/ru


> > Доброе время суток ВСЕМ!
> > Вопрос следующий, имеется некое устройство АЦП, которое собирает данные с частотой дискретизации 260.870 кгц, так вот потом со снятым сигналом требуется сделать бпф и построить спектр короче. Ближайшее подходящее значение степени двойки 262144, и следовательно частотный шаг будет не 1 гц, а мне требуется 1 гц, как тут быть?

> Дополнить секундный отрезок нулями.
> Или посчитать, как есть, а потом интерполировать спектр на нужные точки. Принципиально результат не отличается, хотя в реализации возможны различия (из-за окон и т.п.)

Я делал БПФ где нужно число элементов кратное 8 интересует обращайся MixaD@Xaker.ru


Здравствуйте! Я генерирую отсчёты косинуса, делаю БПФ, нахожу максимальную по модулю гармонику - она достаточно точно совпадает с тем, что я задал, проблема в том, что очень сильно искажается фаза.
Пример: частота дискретизации=1кГц, сигнал: 100cos(2PI*100Гц*t + 0 рад)
моя программа находит в ПФ на 1048576 отсчётах: отсчёт №104858, частота 100.000381Гц, модуль 75.683 (вместо 100), фаза -72 градуса (-1.257 рад) (вместо 0).
Вопрос: почему так?
24 июля 2004 г. 14:40:


Физика в анимациях - Купить диск - Тесты по физике - Графики on-line

Реклама:
Rambler's Top100