Быстрое преобразование Фурье

Привет всем!
Столкнулся я с необходимостью освоения быстрого преобразования Фурье. Освоение это идет с большим скрипом, мозги кипят - а понимание не приходит. Конечно же, существуют готовые библиотеки от разработчиков, но пользоваться ими нет желания, так как хочется вникнуть в суть вопроса, да и человек я такой, что люблю докопаться до истины. И ясного и понятного алгоритма нигде нет.
Реализовывать сначала буду на STM32F103 на имеющейся PB2, потом когда приедет отладочная плата - на STM32F429.
В любой литературе приведен вывод формул из ДПФ в БПФ. Нехитрыми математическими манипуляциями из формулы ДПФ получается вот это уравнение для первой половины коэффициентов:


<Изображение удалено>

И предполагается вычисление частотных составляющим по этим формулам. Но при расчете почему-то в качестве X(k) и Y(k) берутся исходные выборки. Но у нас же выборки это S(n), а не X(k) и Y(k) и согласно формуле Эйлера


<Изображение удалено>

Так что где правда? Может кто-нибудь четко объяснить суть алгоритма БПФ?

0

суть алгоритма бпф на кривых пальцах. берешь массив значений размером 2^n. обычное дискретное пф делается так: для первого коэффициента синус и косинус должны быть с периодом в этот 2^n. т.е. для k=1 будет 1 период на весь массив. для k=2 уложится 2 периода и т.д.. коэффициенты равны сумме попарных произведений этих двух массивов (сигнала и (ко)синуса).

если внимательно посмотреть на точечные графики синусов, мы увидим, что, во-первых, есть несколько точек во времени, когда синус равен одному и тому же числу. во-вторых, вторая половина синуса равна первой со знаком минус. вот подобие формулы дпф для k=1:

a1 = ... + X_n*Sin_n + ... + X_m*Sin_n + ... + X_i*(-Sin_n) + ... + X_j*(-Sin_n) + ...

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

0

C математикой БПФ для алгоритма по основанию 2 разобрался. Написал код на ассемблере, но результат получается неверный из-за переполнения 16 разрядного числа при накоплении результата в цикле (я использую младшие 16 разрядов в памяти для действительной части и старшие 16 разрядов для мнимой). Тут либо придется переписать под 32 разряда, либо как-то выкручиваться, чтобы не возникало переполнение.
Может кто поскажет как можно избежать переполнения?

0

посчитать по алгоритму действительно требующуюся разрядность.

0

Сообщение от Ymk посчитать по алгоритму действительно требующуюся разрядность.

Я не знаю, как по алгоритму это посчитать, но знаю, что результат д.б. в 128 раз больше исходного сигнала... т.е. на 8 разрядов. При входном сигнале АЦП в 12 разрядов получается 20 разрядов, что не помещается в пол слова (16 бит). Т.е. в очередном цикле, когда результат становится больше, происходит переполнение разрядов и результат искажается.
Вот я и подумал, можно ли как-то постепенно понижать разрядность, чтобы не было переполнения?
Не знаю правильно-ли, но может результат каждой бабочки предварительно стоит делить на 2 и только после этого сохранять в память? Но тут я вообще не уверен в правильности результата.
Самое простое использовать 32 бита и не парится особо, но все же для пользы дела думаю стоит и с этим разобраться.

0

Наконец-то осилил БПФ. Долгое время задача не решалась, поэтому даже пришлось забросить ее на пару месяцев. Зато после второго подхода нашлись незамеченные ошибки и все разрешилось влет! Конечно код получился не оптимальный, но все еще впереди... Главное результат правильный. Но было бы неплохо еще протестировать. Если есть желающие могу поделиться кодом.
Но теперь в моей задаче остается один вопросик.
После преобразования получается примерно такая картинка:
№ гармоники | Re | Im |
0..................Re1..Im1
1..................Re2..Im2
... .... ....
n..................Ren..Imn
........................................ ...................... ____________
Амплитуда k-й гармоники вычисляется как: A = V Re^2 + Im^2
Фаза k-й гармоники вычисляется как: ф = arctg(Re/Im)
Что такое фаза? Ведь для вычисления БПФ берется период выборок! Это угол k-й гармоники в начальный момент? Т.е. для 1й выборки?
Как вычислить arctg на ассемблере?
Как найти по вычисленным гармоникам амплитуду исходного сигнала в нужный момент времени?

0

Сообщение от wypuk Амплитуда k-й гармоники вычисляется как: A = V Re^2 + Im^2 Фаза k-й гармоники вычисляется как: ф = arctg(Re/Im) Что такое фаза? Ведь для вычисления БПФ берется период выборок! Это угол k-й гармоники в начальный момент? Т.е. для 1й выборки?

Фаза - это угол гармоники в сигнале, по отношению к Вашим осям координат. Для каждой гармоники фаза своя.

Сообщение от wypuk Как вычислить arctg на ассемблере?

Проще всего - таблично, но сначала желательно выделить нужную четверть.

Сообщение от wypuk Как найти по вычисленным гармоникам амплитуду исходного сигнала в нужный момент времени?

"В нужный момент времени" это значит внутри периода выборок? Если да, то никак. Фурье дает среднюю амплитуду.

0

Сообщение от ZIvS Фаза - это угол гармоники в сигнале, по отношению к Вашим осям координат. Для каждой гармоники фаза своя.

Я понимаю, что это угол. Но в какой момент времени? Например возьмем определенный отрезок сигнала, такой, что в нем помещается к примеру 1 период 1-й гармоники 2 периода 2-й и т.д. У каждой из них имеется какой-то начальный сдвиг фазы относительно "начала координат" (первая выборка). Так вот фаза - это этот самый сдвиг фазы?
На рисунке дФ?


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

0

Сообщение от wypuk Так вот фаза - это этот самый сдвиг фазы? На рисунке дФ?

Да

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

- собственно, восстановление сигнала по спектру - это ОПФ (ПДПФ, или ОБПФ).
Вы же, грубо говоря, получили спектр периодического сигнала с периодом, равным длине выборки. применяя обратно преобразование - из спектра получите этот сигнал.

0

Сообщение от wypuk ... У каждой из них имеется какой-то начальный сдвиг фазы относительно "начала координат" (первая выборка). Так вот фаза - это этот самый сдвиг фазы?

Да. Получаемый угол есть угол разницы между опорной косинусойдой/синусойдой и гармоникой в сигнале.

Фактически - фаза в сигнале всегда вращается, но так-как мы умножаем сигнал на косинусоиду/синусоиду мы переносим гармонику на нулевую частоту (как это происходит в приемниках с нулевой ПЧ и квадратурных демодуляторах), далее выделяем постоянную составляющую - амплитуду этой нулевой частоты.

Представьте ситуацию - у Вас есть одна косинусоида в сигнале, и Вы ее оцифровали и провели Фурье. Если бы Вам так повезло, что косинусоида в сигнале совпала-бы по фазе с опорной косинусоидой преобразования, то на оси Im Вы получите ноль и соответственно вся амплитуда сигнала есть Re, а фаза будет равна нулю.
Вот этот сдвиг фаз и есть получаемый угол.

Да, первая выборка здесь играет роль, Вы правы. Но нужно убедиться, что для всех частот в первой выборке опорные косинус/синус имеют равную свою фазу.
Это зависит от топологии преобразования.

0

Сообщение от Mykysoft - собственно, восстановление сигнала по спектру - это ОПФ (ПДПФ, или ОБПФ). Вы же, грубо говоря, получили спектр периодического сигнала с периодом, равным длине выборки. применяя обратно преобразование - из спектра получите этот сигнал.

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


<Изображение удалено>

0

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

Сообщение от virt Этот ряд сходится крайне медленно. Поищите, есть очень быстрые ряды для arctg - буквально 4-5 слагаемых с ошибкой в пятом знаке.

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

0

Вернусь послезавтра домой, гляну в справочнике по рядам, если сами раньше не найдёте.
Неужели в интернете не одного талмуда не выложено?

Вот более калькуляберный ряд


<Изображение удалено>

Сообщение от virt Вот более калькуляберный ряд

что-то этот ряд при t > 1 не сходится совершенно. Получается что я могу вычислять углы только до 45 градусов с его помощью! Не катит.

0

Сообщение от wypuk Сообщение от virt Вот более калькуляберный ряд

что-то этот ряд при t > 1 не сходится совершенно. Получается что я могу вычислять углы только до 45 градусов с его помощью! Не катит.

Ну, можете при t>1 поменять местами синус и косинус и вычесть результат из 90 градусов.

То есть arctg(t) = 90-arctg(1/t)

0

Вообще-то для лучшей сходимости ряда аргумент уменьшают. Например так: http://algorithm.narod.ru/el/arctan.html

0

Сообщение от wypuk ... углы ...

У Вас аргумент сколько разрядов?

0

Сообщение от ZIvS У Вас аргумент сколько разрядов?

Незнаю, я пока не разобрался как его посчитать. А какое это имеет значение?

Я тут подумал, может мне вообще не стоит вычислять аргумент, мне он не особо-то и нужен... для моей задачи нужно знать углы векторов синусоидального напряжения, чтобы можно было поворачивать эти вектора на определенные углы. Re и Im я вычислил, длину вектора тоже научился вычислять. sin и cos рассчитать несложно. Может стоить вспомнить тригонометрию и: sin(a+b) = sin(a)cos(b) + sin(b)cos(a). Так наверное даже проще будет.

0