Алгоритм БПФ Кули – Тьюки

Алгоритм быстрого преобразования Фурье

Алгоритм Кули – Тьюки, названный в честь Дж. W. Cooley и John Tukey, является наиболее распространенным алгоритмом быстрого преобразования Фурье (FFT). Он повторно выражает дискретное преобразование Фурье (ДПФ) произвольного составного размера N = N 1N2в терминах N 1 меньших ДПФ размером N 2, рекурсивно, чтобы сократить время вычисления до O (N log N) для сильно составных N (сглаженных чисел ). Из-за важности алгоритма конкретные варианты и стили реализации стали известны под своими именами, как описано ниже.

Поскольку алгоритм Кули – Тьюки разбивает ДПФ на более мелкие ДПФ, его можно произвольно комбинировать с любым другим алгоритмом ДПФ. Например, алгоритм Рейдера или Блюстайна может использоваться для обработки больших простых множителей, которые не могут быть разложены Кули-Тьюки, или алгоритм простых множителей может быть используется для большей эффективности при выделении относительно простых факторов.

Алгоритм и его рекурсивное приложение были изобретены Карлом Фридрихом Гауссом. Кули и Тьюки независимо друг от друга открыли и популяризировали его 160 лет спустя.

Содержание

• 1 История
• 2 Случай DIT с основанием 2
  • 2.1 Псевдокод
• 3 Идея
• 4 Варианты
• 5 Переупорядочение данных, изменение битов и алгоритмы на месте
• 6 Ссылки
• 7 Внешние ссылки

История

Этот алгоритм, включая его рекурсивное приложение, был изобретен около 1805 года Карлом Фридрихом Гауссом, который использовал его для интерполяции траекторий. из астероидов Паллада и Юноны, но его работа не получила широкого признания (была опубликована только посмертно и на неолатинском ). Однако Гаусс не анализировал асимптотическое время вычисления. Различные ограниченные формы также были повторно открыты несколько раз на протяжении 19 и начала 20 веков. БПФ стали популярными после того, как Джеймс Кули из IBM и Джон Тьюки из Принстон опубликовали в 1965 году статью, в которой заново изобретали алгоритм и описывали, как его выполнять. это удобно на компьютере.

Тьюки, как сообщается, придумал эту идею во время встречи Научного консультативного комитета президента Кеннеди, на котором обсуждались способы обнаружения испытаний ядерного оружия в Советском Союзе с помощью сейсмометров, расположенных за пределами страны. Эти датчики будут генерировать сейсмологические временные ряды. Однако для анализа этих данных потребуются быстрые алгоритмы вычисления ДПФ из-за количества датчиков и продолжительности времени. Эта задача имела решающее значение для ратификации предложенного запрета на ядерные испытания, чтобы любые нарушения могли быть обнаружены без необходимости посещения советских объектов. Другой участник той встречи, Ричард Гарвин из IBM, осознал потенциал этого метода и связал Тьюки с Кули, убедившись, что Кули не знал первоначальной цели. Вместо этого Кули сказали, что это необходимо для определения периодичности ориентации спинов в трехмерном кристалле гелия-3. Впоследствии Кули и Тьюки опубликовали свой совместный документ, который быстро получил широкое распространение благодаря одновременной разработке аналого-цифровых преобразователей, способных осуществлять дискретизацию с частотой до 300 кГц.

Тот факт, что Гаусс описал тот же алгоритм (хотя и без анализа его асимптотической стоимости), был реализован лишь через несколько лет после статьи Кули и Тьюки 1965 года. В своей статье они процитировали в качестве вдохновения только работу И. J. Good о том, что сейчас называется алгоритмом БПФ с простым коэффициентом (PFA); Хотя алгоритм Гуда изначально считался эквивалентным алгоритму Кули-Тьюки, быстро стало понятно, что PFA - это совсем другой алгоритм (работающий только для размеров, которые имеют относительно простые множители и полагаются на Китайская теорема об остатках, в отличие от поддержки любого составного размера в Кули-Тьюки).

Случай DIT с основанием 2

A счисление-2 прореживание во времени (DIT ) БПФ является самой простой и наиболее распространенной формой алгоритма Кули – Тьюки, хотя высокооптимизированные реализации Кули – Тьюки обычно используют другие формы алгоритма, как описано ниже. Radix-2 DIT делит DFT размера N на два перемежающихся DFT (отсюда и название "radix-2") размера N / 2 на каждом этапе рекурсии.

Дискретное преобразование Фурье (ДПФ) определяется по формуле:

$X\; k\; =\; \sum \; n\; =\; 0\; N\; -\; 1\; xne\; -\; 2\; \pi \; i\; N\; nk,\; \{\backslash \; displaystyle\; X\_\; \{k\}\; =\; \backslash \; sum\; \_\; \{n\; =\; 0\}\; ^\; \{N-1\}\; x\_\; \{n\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; nk\},\}$

где $k\; \{\backslash \; displaystyle\; k\}$- целое число от $0\; \{\backslash \; displaystyle\; 0\}$до $N\; -\; 1\; \{\backslash \; displaystyle\; N-1\}$.

Radix-2 DIT сначала вычисляет ДПФ входов с четным индексом $(x\; 2\; m\; =\; x\; 0,\; x\; 2,\dots ,\; x\; N\; -\; 2)\; \{\backslash \; displaystyle\; (x\_\; \{2m\}\; =\; x\_\; \{0\},\; x\_\; \{2\},\; \backslash \; ldots,\; x\_\; \{N-2\})\}$и входов с нечетным индексом $(x\; 2\; m\; +\; 1\; =\; x\; 1,\; x\; 3,\dots ,\; x\; N\; -\; 1)\; \{\backslash \; displaystyle\; (x\_\; \{\; 2m\; +\; 1\}\; =\; x\_\; \{1\},\; x\_\; \{3\},\; \backslash \; ldots,\; x\_\; \{N-1\})\}$, а затем объединяет эти два результата для получения ДПФ всей последовательности. Затем эту идею можно выполнить рекурсивно, чтобы сократить общее время выполнения до O (N log N). Эта упрощенная форма предполагает, что N является степенью двойки ; поскольку количество точек выборки N обычно может быть свободно выбрано приложением (например, путем изменения частоты дискретизации или окна, заполнения нулями и т. д.), это часто не является важным ограничением.

Алгоритм DIT radix-2 переупорядочивает ДПФ функции $xn\; \{\backslash \; displaystyle\; x\_\; \{n\}\}$на две части: сумму по четным индексам $n\; =\; 2\; m\; \{\backslash \; displaystyle\; n\; =\; \{2m\}\}$и сумма по нечетным индексам $n\; =\; 2\; m\; +\; 1\; \{\backslash \; displaystyle\; n\; =\; \{2m\; +\; 1\}\}$:

$X\; k\; =\; \sum \; m\; =\; 0\; N\; /\; 2\; -\; 1\; x\; 2\; me\; -\; 2\; \pi \; i\; N\; (2\; m)\; k\; +\; \sum \; m\; =\; 0\; N\; /\; 2\; -\; 1\; x\; 2\; m\; +\; 1\; e\; -\; 2\; \pi \; i\; N\; (2\; м\; +\; 1)\; К\; \{\backslash \; Displaystyle\; \{\backslash \; begin\; \{matrix\}\; X\_\; \{k\}\; =\; \backslash \; sum\; \backslash \; limits\; \_\; \{m\; =\; 0\}\; ^\; \{N\; /\; 2-1\}\; x\_\; \{2m\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; (2m)\; k\}\; +\; \backslash \; sum\; \backslash \; limits\; \_\; \{m\; =\; 0\}\; ^\; \{N\; /\; 2-1\}\; x\_\; \{2m\; +\; 1\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; (2m\; +\; 1)\; k\}\; \backslash \; end\; \{matrix\}\}\}$

Можно разложить на множитель общий множитель $e\; -\; 2\; \pi \; i\; N\; k\; \{\backslash \; displaystyle\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; k\}\}$из второй суммы, как показано в уравнении ниже. Тогда ясно, что две суммы - это ДПФ части с четным индексом $x\; 2\; m\; \{\backslash \; displaystyle\; x\_\; \{2m\}\}$и ДПФ части с нечетным индексом $x\; 2\; m\; +\; 1\; \{\backslash \; displaystyle\; x\_\; \{2m\; +\; 1\}\}$функции $xn\; \{\backslash \; displaystyle\; x\_\; \{n\}\}$. Обозначим ДПФ E Ven-индексированных входов $x\; 2\; m\; \{\backslash \; displaystyle\; x\_\; \{2m\}\}$как $E\; k\; \{\backslash \; displaystyle\; E\_\; \{k\}\}$и ДПФ O индексированных dd входов $x\; 2\; m\; +\; 1\; \{\backslash \; displaystyle\; x\_\; \{2m\; +\; 1\}\}$by $O\; k\; \{\backslash \; displaystyle\; O\_\; \{k\}\}$, и мы получаем:

$X\; k\; =\; \sum \; m\; =\; 0\; N\; /\; 2\; -\; 1\; x\; 2\; me\; -\; 2\; \pi \; i\; N\; /\; 2\; mk\; ⏟\; DFT\; ofeven\; -\; indexedpartofxn\; +\; e\; -\; 2\; \pi \; i\; N\; k\; \sum \; m\; =\; 0\; N\; /\; 2\; -\; 1\; x\; 2\; m\; +\; 1\; e\; -\; 2\; \pi \; i\; N\; /\; 2\; mk\; ⏟\; DFT\; ofodd\; -\; indexedpartofxn\; =\; E\; k\; +\; e\; -\; 2\; \pi \; i\; N\; k\; O\; k.\; \{\backslash \; Displaystyle\; \{\backslash \; begin\; \{matrix\}\; X\_\; \{k\}\; =\; \backslash \; underbrace\; \{\backslash \; sum\; \backslash \; limits\; \_\; \{m\; =\; 0\}\; ^\; \{N\; /\; 2-1\}\; x\_\; \{2m\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\; /\; 2\}\}\; mk\}\}\; \_\; \{\backslash \; mathrm\; \{DFT\; \backslash ;\; of\; \backslash ;\; с\; четным\; индексом\; \backslash ;\; part\; \backslash ;\; of\; \backslash ;\}\; x\_\; \{n\}\}\; \{\}\; +\; e\; ^\; \{-\; \{\backslash \; frac\; \{\; 2\; \backslash \; pi\; i\}\; \{N\}\}\; k\}\; \backslash \; underbrace\; \{\backslash \; sum\; \backslash \; limits\; \_\; \{m\; =\; 0\}\; ^\; \{N\; /\; 2-1\}\; x\_\; \{2m\; +\; 1\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\; \}\; \{N\; /\; 2\}\}\; mk\}\}\; \_\; \{\backslash \; mathrm\; \{DFT\; \backslash ;\; of\; \backslash ;\; odd-indexed\; \backslash ;\; part\; \backslash ;\; of\; \backslash ;\}\; x\_\; \{n\}\}\; =\; E\_\; \{k\}\; +\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; k\}\; O\_\; \{k\}.\; \backslash \; end\; \{matrix\}\}\}$

Благодаря периодичности комплексной экспоненты, $X\; k\; +\; N\; 2\; \{\backslash \; displaystyle\; X\_\; \{k\; +\; \{\backslash \; frac\; \{N\}\; \{2\}\}\}\}$также получается из $E\; k\; \{\backslash \; displaystyle\; E\_\; \{k\}\}$и $О\; К\; \{\backslash \; Displaystyle\; O\_\; \{k\}\}$:

$Икс К + N 2 = ∑ м = 0 N / 2 - 1 x 2 меня - 2 π я N / 2 м (k + N 2) + е - 2 π i N (k + N 2) ∑ m = 0 N / 2 - 1 x 2 m + 1 e - 2 π i N / 2 m (k + N 2) = ∑ m = 0 N / 2 - 1 x 2 me - 2 π i N / 2 mke - 2 π mi + e - 2 π i N ke - π i ∑ m = 0 N / 2 - 1 x 2 m + 1 e - 2 π i N / 2 mke - 2 π mi = ∑ m = 0 N / 2 - 1 x 2 me - 2 π i N / 2 mk - e - 2 π i N k ∑ m = 0 N / 2 - 1 x 2 m + 1 e - 2 π i N / 2 мк знак равно Е К - е - 2 π я N К О К {\ Displaystyle {\ begin {выровнено} X_ {k + {\ frac {N} {2}}} = \ sum \ limits _ {m = 0} ^ {N / 2-1} x_ {2m} e ^ {- {\ frac {2 \ pi i} {N / 2}} m (k + {\ frac {N} {2}})} + e ^ { - {\ frac {2 \ pi i} {N}} (k + {\ frac {N} {2}})} \ sum \ limits _ {m = 0} ^ {N / 2-1} x_ {2m + 1} e ^ {- {\ frac {2 \ pi i} {N / 2}} m (k + {\ frac {N} {2}})} \\ = \ sum \ limits _ {m = 0} ^ {N / 2-1} x_ {2m} e ^ {- {\ frac {2 \ pi i} {N / 2}} mk} e ^ {- 2 \ pi mi} + e ^ {- {\ frac {2 \ pi i} {N}} k} e ^ {- \ pi i} \ sum \ limits _ {m = 0} ^ {N / 2-1} x_ {2m + 1} e ^ {- {\ frac {2 \ pi i} {N / 2}} mk} e ^ {- 2 \ pi mi} \\ = \ sum \ limits _ {m = 0} ^ {N / 2-1} x_ {2m} e ^ {- {\ frac {2 \ pi i} {N / 2}} mk} -e ^ {- {\ frac {2 \ pi i} {N}} k} \ sum \ limits _ {m = 0 } ^ {N / 2-1} x_ {2m + 1} e ^ {- {\ frac {2 \ pi i} {N / 2}} mk} \\ = E_ {k} -e ^ {- { \ frac {2 \ pi i} {N}} k} O_ {k} \ end {align}}}$

Мы можем переписать $X\; k\; \{\backslash \; displaystyle\; X\_\; \{k\}\}$как:

$Икс\; К\; знак\; равно\; Е\; К\; +\; е\; -\; 2\; \pi \; я\; N\; К\; О\; К\; Икс\; К\; +\; N\; 2\; =\; Е\; К\; -\; е\; -\; 2\; \pi \; я\; N\; К\; О\; К\; \{\backslash \; Displaystyle\; \{\backslash \; begin\; \{matrix\}\; X\_\; \{k\; \}\; =\; E\_\; \{k\}\; +\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; k\}\; O\_\; \{k\}\; \backslash \backslash \; X\_\; \{k\; +\; \{\backslash \; frac\; \{N\}\; \{2\}\}\}\; =\; E\_\; \{k\}\; -e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\}\}\; \{k\}\}\; O\_\; \{k\}\; \backslash \; end\; \{matrix\}\}\}$

Этот результат, выражающий рекурсивное ДПФ длины N в терминах двух ДПФ размера N / 2 является ядром быстрого преобразования Фурье ДПФ с основанием 2. Алгоритм набирает скорость за счет повторного использования результатов промежуточных вычислений для вычисления нескольких выходных данных ДПФ. Обратите внимание, что окончательные результаты получаются +/- комбинацией $E\; k\; \{\backslash \; displaystyle\; E\_\; \{k\}\}$и $O\; k\; exp\; \; (-\; 2\; \pi \; ik\; /\; N)\; \{\backslash \; displaystyle\; O\_\; \{k\}\; \backslash \; exp\; (-2\; \backslash \; pi\; ik\; /\; N)\}$, который представляет собой просто ДПФ размера 2 (в данном контексте иногда называется бабочка ); когда это обобщается на более крупные координаты ниже, ДПФ размера 2 заменяется ДПФ большего размера (которое само по себе может быть оценено с помощью БПФ).

Схема потока данных для N = 8: БПФ с децимацией во времени по основанию 2 разбивает ДПФ длины N на два ДПФ длины N / 2, за которыми следует этап комбинирования, состоящий из множества ДПФ размера 2, называемый "бабочка" "операции (так называемые из-за формы диаграмм потоков данных).

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

Вышеупомянутое перевыражение ДПФ размера N в виде двух ДПФ размера N / 2 иногда называют Дэниелсоном - леммой Ланцоша , так как идентичность была отмечена этими двумя авторами в 1942 году (под влиянием работы Рунге 1903 года). Они применили свою лемму «обратным» рекурсивным способом, многократно удваивая размер ДПФ до тех пор, пока спектр преобразования не сойдется (хотя они, по-видимому, не осознавали линеарифмическую [т.е. порядок N log N] асимптотической сложности, которую они имели достигнуто). Работа Дэниэлсона – Ланцоша предшествовала широкому распространению механических или электронных компьютеров и требовала ручных вычислений (возможно, с помощью таких механических средств, как счетные машины ); они сообщили о времени вычисления 140 минут для ДПФ размером 64, работающего на реальных входных данных с точностью до 3-5 значащих цифр. В статье Кули и Тьюки 1965 года сообщается, что время обработки сложного ДПФ размером 2048 на IBM 7094 (вероятно, 36-битное одинарной точности, ~ 8 цифр) составляет 0,02 минуты. Если изменить масштаб времени на количество операций, это примерно соответствует коэффициенту ускорения около 800000. (Чтобы представить себе время ручного вычисления в перспективе, 140 минут для размера 64 соответствуют в среднем максимум 16 секундам на операцию с плавающей запятой, около 20% из которых являются умножениями.)

Псевдокод

В псевдокоде можно записать следующую процедуру:

X0,..., N − 1 ← ditfft2 (x, N, s) : ДПФ (x 0, x s, x 2s,..., x (N-1) s): если N = 1, то X0← x 0 базовый случай ДПФ размером 1 иначе X 0,..., N / 2−1 ← ditfft2 (x, N / 2, 2s) ДПФ (x 0, x 2s, x 4s,...) X N/2,...,N−1 ← ditfft2 (x + s, N / 2, 2s) ДПФ (x s, x s+ 2s, x s+ 4s,...) для k = 0 до N / 2−1 do объединить ДПФ двух половин в полное ДПФ: t ← X kXk← t + exp (−2πi k / N) X k+ N / 2 Xk+ N / 2 ← t - exp (−2πi k / N) X k+ N / 2 конец для end if

Здесь ditfft2(x, N, 1), вычисляет X = DFT (x) не на своем месте по основанию 2 D IT FFT, где N - целая степень двойки, а s = 1 - шаг входного массива x . x + s обозначает массив, начинающийся с x s.

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

Высокопроизводительные реализации БПФ вносят много изменений в реализацию такого алгоритма по сравнению с этим простым псевдокодом. Например, можно использовать больший базовый случай, чем N = 1, чтобы амортизировать накладные расходы на рекурсию, коэффициенты вращения $exp\; \; [-\; 2\; \pi \; ik\; /\; N]\; \{\; \backslash \; displaystyle\; \backslash \; exp\; [-2\; \backslash \; pi\; ik\; /\; N]\}$может быть вычислено заранее, и для cache часто используются большие радиусы; вместе эти и другие оптимизации могут улучшить производительность на порядок или больше. (Во многих реализациях из учебников рекурсия в глубину полностью устранена в пользу нерекурсивного подхода в ширину, хотя утверждалось, что рекурсия в глубину имеет лучшую память для locality.) Некоторые из этих идей более подробно описаны ниже.

Идея

Базовый шаг БПФ Кули – Тьюки для общих факторизаций можно рассматривать как переинтерпретацию 1d DFT как нечто вроде 2d DFT. Входной массив 1d длины N = N 1N2переинтерпретируется как матрица 2d N 1×N2, хранящаяся в порядке по столбцам. Один выполняет меньшие 1d DFT в направлении N 2 (несмежное направление), затем умножается на фазовые коэффициенты (коэффициенты вращения) и, наконец, выполняет 1d DFT вдоль направления N 1. Шаг транспонирования может быть выполнен в середине, как показано здесь, или в начале или в конце. Это делается рекурсивно для меньших преобразований. Иллюстрация порядка строк и столбцов

В более общем смысле, алгоритмы Кули-Тьюки рекурсивно повторно выражают ДПФ составного размера N = N 1N2как:

1. Выполнить N 1 ДПФ размера N 2.
2. Умножить на комплексные корни из единицы (часто называемые множителями ).
3. Выполнить N 2 ДПФ размера N 1.

Обычно либо N 1, либо N 2 является малым множителем (не обязательно простым), называемым основанием (которое может различаются между этапами рекурсии). Если N 1 является основанием системы счисления, это называется алгоритмом децимации во времени (DIT), тогда как если N 2 - это по основанию системы счисления, это прореживание по частоте (DIF, также называемый алгоритмом Сэнде – Тьюки). Представленная выше версия была алгоритмом DIT с основанием 2; в окончательном выражении фаза умножения нечетного преобразования равна коэффициент поворота и +/- комбинация (бабочка) четных и нечетных преобразований - это ДПФ размера 2 (малое ДПФ системы счисления является некоторым время, известное как бабочка, так называемое из-за формы диаграммы потока данных для случая с основанием 2.)

Варианты

Есть много других вариантов алгоритма Кули – Тьюки. Реализации со смешанным основанием обрабатывают составные размеры с различными (обычно небольшими) факторами в дополнение к двум, обычно (но не всегда) используя алгоритм O (N) для простых базовых случаев рекурсии (это также возможно использовать алгоритм N log N для простых базовых случаев, такой как алгоритм Rader или Bluestein ). Разделение системы счисления объединяет основания 2 и 4, используя тот факт, что первое преобразование системы счисления 2 не требует коэффициента поворота, чтобы достичь того, что долгое время было наименьшим известным числом арифметических операций для размеров степени двойки, хотя последние вариации достигают еще меньшего числа. (На современных компьютерах производительность определяется больше соображениями кеша и конвейера ЦП, чем строгим подсчетом операций; хорошо оптимизированные реализации БПФ часто используют большие радиусы и / или жестко запрограммированные базовые преобразования значительного размера.).

Другой способ взглянуть на алгоритм Кули-Тьюки состоит в том, что он повторно выражает одномерное ДПФ размера N как N 1 на N 2 двух- размерное ДПФ (плюс твидллы), где выходная матрица транспонирована. Чистый результат всех этих транспозиций для алгоритма radix-2 соответствует инверсии битов входных (DIF) или выходных (DIT) индексов. Если вместо использования малой системы счисления используется основание системы счисления примерно √N и явные транспозиции матрицы ввода / вывода, это называется четырехшаговым алгоритмом (или шестиэтапным, в зависимости от количества транспозиции), первоначально предложенные для улучшения локальности памяти, например для оптимизации кэша или внеядерной операции, и позже было показано, что это оптимальный алгоритм без учета кеша.

Общая факторизация Кули-Тьюки переписывает индексы k и n как $К\; =\; N\; 2\; К\; 1\; +\; К\; 2\; \{\backslash \; Displaystyle\; к\; =\; N\_\; \{2\}\; k\_\; \{1\}\; +\; k\_\; \{2\}\}$и $n\; =\; N\; 1\; n\; 2\; +\; n\; 1\; \{\backslash \; displaystyle\; n\; =\; N\_\; \{1\}\; n\_\; \{2\}\; +\; n\_\; \{1\}\}$, соответственно, где индексы k a и n a начинаются с 0..N a -1 (для 1 или 2). То есть он повторно индексирует вход (n) и выход (k) как N 1 на N 2 двумерных массивов в основных столбцах и порядок старших строк соответственно; разница между этими индексами - это транспозиция, как упоминалось выше. Когда это повторное индексирование подставляется в формулу ДПФ для nk, $N\; 1\; n\; 2\; N\; 2\; k\; 1\; \{\backslash \; displaystyle\; N\_\; \{1\}\; n\_\; \{2\}\; N\_\; \{2\}\; k\_\; \{1\}\}$перекрестный член равен нулю (его экспонента равна единице), а оставшиеся члены дают

$XN\; 2\; k\; 1\; +\; k\; 2\; =\; \sum \; n\; 1\; =\; 0\; N\; 1\; -\; 1\; \sum \; n\; 2\; =\; 0\; N\; 2\; -\; 1\; x\; N\; 1\; n\; 2\; +\; N\; 1\; е\; -\; 2\; \pi \; я\; N\; 1\; N\; 2\; \cdot \; (N\; 1\; N\; 2\; +\; N\; 1)\; \cdot \; (N\; 2\; К\; 1\; +\; К\; 2)\; \{\backslash \; Displaystyle\; X\_\; \{N\_\; \{2\}\; k\_\; \{1\}\; +\; k\_\; \{2\; \}\}\; =\; \backslash \; sum\; \_\; \{n\_\; \{1\}\; =\; 0\}\; ^\; \{N\_\; \{1\}\; -1\}\; \backslash \; sum\; \_\; \{n\_\; \{2\}\; =\; 0\}\; ^\; \{N\_\; \{2\}\; -1\}\; x\_\; \{N\_\; \{1\}\; n\_\; \{2\}\; +\; n\_\; \{1\}\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{1\}\; N\_\; \{2\}\}\}\; \backslash \; cdot\; (N\_\; \{1\}\; n\_\; \{2\}\; +\; n\_\; \{1\})\; \backslash \; cdot\; (N\_\; \{2\}\; k\_\; \{1\}\; +\; k\_\; \{2\})\}\}$

$=\; \sum \; n\; 1\; =\; 0\; N\; 1\; -\; 1\; [e\; -\; 2\; \pi \; i\; N\; 1\; N\; 2\; n\; 1\; k\; 2]\; (\sum \; n\; 2\; знак\; равно\; 0\; N\; 2\; -\; 1\; Икс\; N\; 1\; N\; 2\; +\; N\; 1\; е\; -\; 2\; \pi \; я\; N\; 2\; N\; 2\; К\; 2)\; е\; -\; 2\; \pi \; я\; N\; 1\; N\; 1\; К\; 1\; \{\backslash \; Displaystyle\; =\; \backslash \; sum\; \_\; \{n\_\; \{1\}\; =\; 0\}\; ^\; \{N\_\; \{1\}\; -1\}\; \backslash \; left\; [e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{1\}\; N\_\; \{2\}\}\}\; n\_\; \{1\}\; k\_\; \{2\}\}\; \backslash \; right\; ]\; \backslash \; left\; (\backslash \; sum\; \_\; \{n\_\; \{2\}\; =\; 0\}\; ^\; \{N\_\; \{2\}\; -1\}\; x\_\; \{N\_\; \{1\}\; n\_\; \{2\}\; +\; n\_\; \{1\}\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{2\}\}\}\; n\_\; \{2\}\; k\_\; \{2\}\}\; \backslash \; right)\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{1\}\}\}\; n\_\; \{1\}\; k\_\; \{1\}\; \}\}$

$=\; \sum \; n\; 1\; =\; 0\; N\; 1\; -\; 1\; (\sum \; n\; 2\; =\; 0\; N\; 2\; -\; 1\; x\; N\; 1\; n\; 2\; +\; n\; 1\; e\; -\; 2\; \pi \; i\; N\; 2\; n\; 2\; k\; 2)\; e\; -\; 2\; \pi \; iN\; 1\; N\; 2\; N\; 1\; (N\; 2\; К\; 1\; +\; К\; 2)\; \{\backslash \; displaystyle\; =\; \backslash \; sum\; \_\; \{n\_\; \{1\}\; =\; 0\}\; ^\; \{N\_\; \{1\}\; -1\}\; \backslash \; left\; (\backslash \; sum\; \_\; \{n\_\; \{2\}\; =\; 0\}\; ^\; \{N\_\; \{2\}\; -1\}\; x\_\; \{N\_\; \{1\}\; n\_\; \{2\}\; +\; n\_\; \{1\}\}\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{2\}\}\}\; n\_\; \{\; 2\}\; k\_\; \{2\}\}\; \backslash \; right)\; e\; ^\; \{-\; \{\backslash \; frac\; \{2\; \backslash \; pi\; i\}\; \{N\_\; \{1\}\; N\_\; \{2\}\}\}\; n\_\; \{1\}\; (N\_\; \{2\}\; k\_\; \{1\}\; +\; k\_\; \{\; 2\})\}\}$.

где каждая внутренняя сумма представляет собой ДПФ размера N 2, каждая внешняя сумма представляет собой ДПФ размера N 1, а [...] термин в квадратных скобках - это фактор вращения.

Может использоваться произвольное основание системы счисления r (а также смешанное основание), как было показано Кули и Тьюки, а также Гауссом (который привел примеры шагов с основанием счисления-3 и основание-6). Кули и Тьюки первоначально предположили, что основание системы счисления требует O (r) работы, и, следовательно, сочли сложность для системы счисления r равной O (r N / r log r N) = O (N log 2 (N) r / log 2 r); из расчета значений r / log 2 r для целых значений r от 2 до 12 найдено оптимальное основание системы счисления 3 (ближайшее целое число к e, которое минимизирует r / log 2 r). Однако этот анализ был ошибочным: основание-бабочка также является ДПФ и может быть выполнено с помощью алгоритма БПФ за O (r log r) операций, следовательно, основание r фактически сокращается в сложности O (r log (r) N / r log r N), а оптимальное r определяется из более сложных соображений. На практике очень большое значение r (32 или 64) важно для эффективного использования, например большое количество регистров процессора на современных процессорах, и даже неограниченное основание системы счисления r = √N также достигает сложности O (N log N) и имеет теоретические и практические преимущества для больших N, как упоминалось выше.

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

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

Самый известный метод переупорядочения включает явное перестановку битов для алгоритмов с основанием 2 на месте. Реверс битов - это перестановка, где данные с индексом n записываются в двоичном с цифрами b 4b3b2b1b0(например, 5 цифр для N = 32 входов), передается в индекс с перевернутыми цифрами b 0b1b2b3b4. Рассмотрим последний этап алгоритма DIT с основанием 2, подобный представленному выше, где выходные данные записываются на месте поверх входных: when $E\; k\; \{\backslash \; displaystyle\; E\_\; \{k\}\}$и $O\; k\; \{\backslash \; displaystyle\; O\_\; \{k\}\}$комбинируются с ДПФ размера 2, эти два значения перезаписываются выходными данными. Однако два выходных значения должны находиться в первой и второй половине выходного массива, что соответствует старшему разряду b 4 (для N = 32); тогда как два входа $E\; k\; \{\backslash \; displaystyle\; E\_\; \{k\}\}$и $O\; k\; \{\backslash \; displaystyle\; O\_\; \{k\}\}$чередуются в четном и нечетном элементы, соответствующие младшему значащему биту b 0. Таким образом, чтобы получить вывод в правильном месте, b 0 должен занять место b 4, и индекс станет b 0b4b3b2b1. И для следующего рекурсивного этапа эти 4 младших бита станут b 1b4b3b2. Если вы включите все рекурсивные этапы алгоритма Radix-2 DIT, все биты должны быть перевернуты, и, следовательно, необходимо предварительно обработать ввод ( или постобработка вывода) с небольшим изменением направления, чтобы получить вывод в упорядоченном виде. (Если каждое субпреобразование размера N / 2 должно работать с непрерывными данными, вход DIT предварительно обрабатывается путем инвертирования битов.) Соответственно, если вы выполняете все шаги в обратном порядке, вы получаете алгоритм DIF с основанием 2. с инверсией битов при постобработке (или предварительной обработке соответственно).

Логарифм (журнал), используемый в этом алгоритме, представляет собой логарифм по основанию 2.

Ниже приведен псевдокод для итеративного алгоритма БПФ с основанием 2, реализованного с использованием перестановки с обращением битов.

алгоритм iterative-fft isinput: Массив a из n комплексных значений, где n является степенью 2. вывод: Массив A ДПФ a. bit-reverse-copy (a, A) n ← a.length для s = 1 до log (n) do m ← 2 ω m ← exp (−2πi / m) для k = 0 до n-1 bym doω ← 1 для j = 0 до m / 2 - 1 do t ← ω A [k + j + m / 2] u ← A [k + j] A [k + j] ← u + t A [ k + j + m / 2] ← u - t ω ← ω ω mreturn A

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

алгоритм бит-обратная копия (a, A) isвход: Массив a из n комплексных значений, где n - степень 2. выход: Массив A из размер n. n ← a.length для k = от 0 до n - 1 do A [rev (k)]: = a [k]

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

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

Проблема значительно упрощается, если он неуместен : выходной массив отличается от входного массива или, что эквивалентно, доступен вспомогательный массив равного размера. Алгоритм автосортировки Стокхэма выполняет каждый этап БПФ вне места, обычно выполняя обратную и прямую запись между двумя массивами, транспонируя одну «цифру» индексов на каждом этапе, и был особенно популярен. на архитектурах SIMD. Еще большие потенциальные преимущества SIMD (большее количество последовательных обращений) были предложены для алгоритма Пиза, который также меняет порядок на каждом этапе, но этот метод требует отдельного обращения бит / цифра и O (N log N) хранение. Также можно напрямую применить определение факторизации Кули-Тьюки с явной (глубиной ) рекурсией и малыми корнями, что дает естественный выход вне места без отдельного шага перестановки (как в псевдокоде выше), и можно утверждать, что он имеет преимущества локальности без учета кеша в системах с иерархической памятью.

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

Ссылки

Внешние ссылки

• «Быстрое преобразование Фурье - БПФ». Техника Кули-Тьюки. Статья. 10. Простой педагогический алгоритм radix-2 на C ++
• «KISSFFT». GitHub. Простая реализация Кули – Тьюки со смешанным основанием в C
• Dsplib на GitHub
• «Алгоритм прореживания Radix-2 во времени БПФ». Архивировано из оригинала 31 октября 2017 года. "Алгоритм БПФ по основанию два с прореживанием по времени" (на русском языке).
• "Прореживание Radix-2 в алгоритме частотного БПФ «. Архивировано из оригинала 14 ноября 2017 года. «Алгоритм БПФ по основанию два с прореживанием по частоте» (на русском языке).