Есть функция делающая БПФ (FFT) для комплексного float входа (для длин степени двойки).

У меня есть два вопроса. Искал ответы в интернете, но не нашёл. Нашёл только те же упоминания возможности. Особенно интересует второй вопрос.

1) Как можно сделать за один её вызов два БПФ для действительного входа той же длины?
Знаю что один вход надо записать в действительную часть, а второй в мнимую. Но что делать после вызова? Как разделить результаты?
(Кое что было написано тут, но непонятно: http://electronix.ru/forum/index.php?showtopic=71731)

2) Как можно сделать за один её вызов один БПФ для действительного входа двойной длины?
Этот вопрос важнее, чем первый.

Если можно напишите формулы.

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

Что именно здесь непонятно?

В книжке Р. Лайонс "Цифровая обработка сигналов" про это написано.

Это не на Си, и не математические формулы. Непонятны записи. Что например вот это: x(2:2:end) ?

Это матлаб.
x - вектор
x(2:2:end) - вектор из x(2) x(4) ... x(end)
x(1:2:end) - вектор из x(1) x(3) ... x(end-1)

Спасибо. Буду разбираться, пока до конца всё равно не понятно. Книжку я скачал только что. Тоже спасибо.

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

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

Рассмотрим самый примитивный пример (N=8), когда мы преобразуем комплексную функцию, разложенную на два массива:
double Re[8]; // это действительная часть исходной функции
double Im[8]; // это мнимая часть исходной функции
В том случае, когда преобразуется действительная функция, массив мнимых частей Im[] забит нулями. Но это не значит, что он не нужен, т.к. в нем мы получим амплитуды синусов.

После FFT-преобразования результат получится в той же паре массивов, а исходные значения испортятся. Получится вот что:
Re[0] - постоянная составляющая, Im[0] - обычно остается нулем.
Re[1] и Im[1] - амплитуды сos и sin 1-ой гармоники.
Re[2] и Im[2] - амплитуды сos и sin 2-ой гармоники.
Re[3] и Im[3] - амплитуды сos и sin 3-ей гармоники.
Re[4] - амплитуда сos 4-ой гармоники. Это частота Найквиста, т.к. более высокие таблично не представимы. Im[4] - равна нулю, т.к. у частоты Найквиста не может быть sin-составляющей.
Re[5] и Im[5] - мусор, копия Re[3] и -Im[3], соответственно.
Re[6] и Im[6] - мусор, копия Re[2] и -Im[2], соответственно.
Re[7] и Im[7] - мусор, копия Re[1] и -Im[1], соответственно.

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

Если у нас есть две действительные функции, то мы можем преобразовать их вдвоем сразу, за один прогон алгоритма FFT. Для этого надо первую действительную функцию запихнуть в Re[], как это и было раньше, а вторую записать в Im[], который раньше был обнулен.

После преобразования получим результат преобразования, соответствующий входной функции F1+i*F2 (множитель i из-за того, что F2 в мнимую часть была записана).
Однако, используя свойство зеркальности, можем выделить из суммы результаты, относящиеся в F1 и F2 по отдельности. Делается это так:

Re[0] – постоянная составляющая F1.
Im[0] – постоянная составляющая F2.

(Re[1]+Re[7])/2 - cos-амплитуда F1 для 1-ой гармоники, т.е. средняя арифметическая результата и зеркального мусора.
Почему так? А потому что cos-амплитуда F1 симметрична и от этой операции и останется прежней. Этой операцией она очистится от вклада sin-амплитуды F2, которая имеет противоположные по знаку амплитуды в Re[1] и Re[7].
(Re[1]-Re[7])/2 - sin-амплитуда F1 для 1-ой гармоники, тут уже симметричные косинусы взаимно уничтожились, а синусы остались живы.

Точно так же обрабатываем зеркальные пары:
(Re[2]+Re[6])/2 - cos-амплитуда F1 для 2-ой гармоники.
(Re[2]-Re[6])/2 - sin-амплитуда F1 для 2-ой гармоники.

(Re[3]+Re[5])/2 - cos-амплитуда F1 для 3-ой гармоники.
(Re[3]-Re[5])/2 - sin-амплитуда F1 для 4-ой гармоники.

Частоты Найквиста зеркального отражения не имеет, с ней проще:
Re[4] - cos-амплитуда F1 для частоты Найквиста.
Im[4] - cos-амплитуда F2 для частоты Найквиста.

В расчетах необходимо только уметь вычислять среднее арифметическое, и более никаких трудностей тут нет.

Про зеркальный спектр от действительных данных я знал. Благодаря этой особенности и линейности FFT от суммы двух сигналов я понял как вывести формулы для разделения спекторов. В общем с первым вопросом я разобался. Осталось разобраться со вторым. Всем, кто помогает - спасибо.

Появлился новый вопрос. Что такое частота Найквиста? Что она характеризует? И что будет, если в предыдущем примере Re[4] обнулить и вычислить обратное преобразование? Тоже получится чисто действительный сигнал? (С нулевой мнимой частью.)

Частота Найквиста выглядит как знакопеременный шум:
..., +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, +1, -1, ...
Оттого и формально считается косинусом, имеющим в первой половине периода +1, а во второй -1. Соответствующей ей синус был бы обязан начинаться с нуля, как и все синусы. И таковым он был бы и во второй части периода.

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

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

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

Поэтому частота Найквиста вполне реальна и пренебрегать ею не следует.

Какие хорошие формулы. Теперь всё понятно.

Странно, что люди изобрели БПФ, который делает вдвое больше работы, чем нужно.

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

Numerical Recipes in C
http://apps.nrbook.com/c/index.html

страница 511
"go to page: 511"

всё, что нужно

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

У меня оказывается эта книжка есть на винчестере.


Во время отладки можно заполнять только одну часть, удостовериться, что всё верно. А во время работы программы использовать обе части! Ведь так быстрее в 2 раза. Или тест написать, который будет с некоторым эпсилон сравнивать результаты вычислений (я так и сделал).

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

А я в Интернете нашел.
http://www2.units.it/ipl/students_area/imm...cal_Recipes.pdf

Это стандартный режим работы многих ЦСП (бит-реверсная адресация данных). Ускоряет работу алгоритма БПФ, делая ненужной сортировку. Просто реверс в нужные моменты включается и выключается.

Я имел в виду алгоритм программный.

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

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

Всё же почему в конце действительное число получается? Если X(k+N/2) является разностью вещественного и комплексного? Поясните пожалуйста.

Не понял вопроса.

X1(k) и X2(k) - комплексные числа? А X(k) в конце - зеркальный спектр? Тоже комплексный?

Другой вопрос. Можно ли так же быстро за вызов обратного FFT половинной длины посчитать обратную функцию к действительному преобразованию? (Кажется что ДА)