Бы́строе преобразова́ние Фурье́ (БПФ; англ. Fast Fourier Transform – FFT), класс алгоритмов, объединённых общим названием, позволяющих вычислить дискретное преобразование Фурье (ДПФ) любого размера $n,$ выполнив $O(n\,\log_2n)$ арифметических операций над комплексными числами. Идеи ускорения ДПФ восходят к работам К. Ф. Гаусса 1805 г., позднее многократно переоткрывались. Основополагающая современная работа по ускорению ДПФ появилась в 1965 г. (Cooley. 1965). Эта статья ссылается на работу (Good. 1958), описывающую алгоритм, который сегодня называют алгоритмом Гуда – Томаса (англ. Good – Thomas Prime Factor algorithm). Работа второго автора этого алгоритма, известного физика и прикладного математика Л. Томаса, опубликована в 1963 г. (Thomas. 1963).

Алгоритм Гуда – Томаса

Если $n=n_1n_2$ и наибольший общий делитель $n_1$ и $n_2$ равен $1,$ то алгоритм Гуда – Томаса с помощью китайской теоремы об остатках (см. Дедекиндово кольцо) сводит вычисление ДПФ размера $n$ к $n_1$ вычислениям ДПФ размера $n_2$ и $n_2$ вычислениям ДПФ размера $n_1.$

Алгоритм основан на широко известном алгебраическом факте о том, что если $n=n_1n_2$ и наибольший общий делитель $n_1$ и $n_2$ равен $1,$ то циклическая группа $\mathbb{Z}/n\mathbb{Z}$ изоморфна прямому произведению циклических групп $\mathbb{Z}/n_1\mathbb{Z}$ и $\mathbb{Z}/n_2\mathbb{Z}.$ И эта теорема, и алгоритм Гуда – Томаса обобщаются на случай, когда n разлагается в произведение нескольких попарно взаимно простых множителей.

Учитывая, что наивное вычисление ДПФ размера $k$ требует примерно $2k^2$ операций над комплексными числами, при выполнении, согласно алгоритму Гуда – Томаса, $n_1$ наивных ДПФ размера $n_2$ плюс $n_2$ наивных ДПФ размера $n_1,$ будет произведено примерно $n_1n_2{}^2 + n_2n_1{}^2 = (n_1n_2)(n_2 + n_1)$ операций над комплексными числами. При условии $n_1\approx n_2\approx \sqrt n$ это число операций равно $2n\sqrt{n},$ что существенно меньше, чем примерно $2n^2$ операций, необходимых для наивного вычисления ДПФ размера $n$ по определению.

Алгоритм Кули – Тьюки дискретного преобразования Фурье

В работе Дж. Кули и Дж. Тьюки (Cooley. 1965) фактически доказывается, что если $n=n_1n_2,$ то путём выполнения $(n_1-1)(n_2-1)$ умножений комплексных чисел можно свести вычисление $\text{ДПФ}_{n,w_n}$ размера $n$ к:

• $n_1$ вычислениям ДПФ размера $n_2$ плюс

• $(n_1-1)(n_2-1)$ умножениям на степени первообразного корня $w_n$ плюс

• $n_2$ вычислениям ДПФ размера $n_1.$

Это рассуждение не требует взаимной простоты $n_1$ и $n_2$ и его можно применять рекурсивно, что позволяет ускорить вычисление ДПФ размера $n$ в случаях, когда $n$ является произведением большого числа небольших множителей, например $n$ имеет вид $2^k$ или $3^k,$ $k$ – натуральное число.

Алгоритм Кули – Тьюки быстрого преобразования Фурье (n – степень числа 2 с натуральным показателем)

В этом случае ДПФ последовательно упрощается сведением ДПФ чётного размера $n$ к двум ДПФ половинного размера $n/2$ путём выполнения всего лишь $2n$ операций над комплексными числами. Конструкция такого сведения была опубликована ещё в 1942 г. и носит название леммы Ланцоша – Даниельсона. Если $n$ имеет вид $n=2^k,$ то, итерируя эту конструкцию, можно построить рекурсивный алгоритм, который выполняет ДПФ размера $n=2^k,$ затрачивая не более $2k 2^k = 2n\,\log_2n$ арифметических операций.

Рекурсивная реализация алгоритма Кули – Тьюки (n – степень числа 2 с натуральным показателем)

Следуя доказательству леммы Ланцоша – Даниельсона, можно реализовать такой алгоритм рекурсивной подпрограммой на псевдокоде, параметрами которой будут аргумент n типа цел (англ. integer) и модифицируемый аргумент A[0 . . n-1] – массив чисел типа компл (англ. complex) с индексами от $0$ до $n-1$:

алг БПФрекурс(арг цел n, аргрез компл A[0 . . n-1])

Подпрограмма получает в качестве аргумента целое число $n = 2^k,$ в качестве модифицируемого аргумента – массив A из $n$ комплексных чисел, содержащий исходные данные $(x_0, x_1, x_2, \ldots, x_{n-1}),$ находит числа $(y_0, y_1, y_2, \ldots, y_{n-1})$ по формулам $y_i =\sum_{j=0}^{n-1}{w_n}^{ij} x_j,$ $i = 0, 1, \dots, n-1,$ где $w_n=e^{–2\pi i /n}=\cos(–2\pi i /n) +i \sin(–2 \pi i/n),$ и записывает результат в массив A.

Для $n = 1$ программа БПФрекурс оставляет входной массив без изменений.

Выполнение программы БПФрекурс при $n > 1$ проводится в несколько этапов. Сначала создаётся массив W размера $n/2,$ который заполняется множителями бабочек Фурье, фигурирующими в доказательстве леммы Ланцоша – Даниельсона. Далее проводится т. н. прореживание входных данных по времени: в памяти создаются два новых массива размера $n/2$ каждый, которые заполняются элементами исходного массива с чётными и нечётными индексами. Затем к каждому из этих массивов размера $n/2$ применяется программа БПФрекурс, в результате чего каждый массив заменяется на своё дискретное преобразование Фурье. Наконец, с помощью выполнения $n/2$ бабочек Фурье с заранее вычисленными множителями результаты двух ДПФ размера $n/2$ объединяются в ДПФ исходного массива размера $n.$

алг БПФрекурс(арг цел n, аргрез компл A[0 . . n-1])

дано | n = 2**k и k >= 0

нач цел i

  если n = 1 то выход все

  компл w = (cos(-2*pi/n), sin(-2*pi/n))

  создать компл W[0 . . n/2-1]

  W[0] := 1

  нц для i от 1 до n/2-1 W[i] := W[i-1]*w кц

  создать компл Ачет [0 . . n/2-1]

  создать компл Анечет [0 . . n/2-1]

  нц для i от 0 до n/2-1

    Ачет[i] := A[2*i]

    Анечет[i] := A[2*i + 1]

  кц

  БПФрекурс(n/2, Aчет[0 . . n/2-1])

  БПФрекурс(n/2, Aнечет[0 . . n/2-1])

  нц для i от 0 до n/2-1

    | параллельное присваивание:

    A[i], A[i + n/2] := Aчет[i] + W[i]*Aнечет[i], Aчет[i] - W[i]*Aнечет[i]

  кц

кон

Индукцией по $k\ge0$ нетрудно проверить, что выполнение данной рекурсивной программы требует проведения не более $2k\,2^k = 2n\,\log_2n$ операций над комплекcными числами.

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

Быстрые алгоритмы для проведения ДПФ размера n существуют для любых значений $n,$ в том числе и для простых $n.$ Например, алгоритм Блюстейна выполняет ДПФ произвольного размера с помощью алгоритма быстрой свёртки (быстрого умножения многочленов) и для любого $n$ вычисляет ДПФ размера $n,$ выполняя $O(n\,\log_2n)$ арифметических операций над комплексными числами.

Нерекурсивная реализация алгоритма Кули – Тьюки (n – степень числа 2 с натуральным показателем)

Описанный выше рекурсивный алгоритм Кули – Тьюки требует выделения памяти для хранения промежуточных результатов. Легко заметить, что этот алгоритм проводит много пересылок элементов различных массивов до того, как над ними начнут выполняться какие-либо действия. Этой лишней работы позволяет избежать итеративная реализация БПФ Кули – Тьюки для случая, когда $n$ имеет вид $2^k,$ где $k$ – натуральное число.

Изобразим графически дерево входных массивов рекурсивного БПФ Кули – Тьюки для случая $n = 8 = 2^3,$ располагая каждую пару вновь создаваемых массивов ниже того массива, из которого они были получены, и размещая при этом в каждом случае подмассив чётных элементов левее массива нечётных элементов.

$$\begin{matrix} \colorbox{#d7dbd8}{$x_0\qquad x_1\qquad x_2\qquad x_3\qquad x_4\qquad x_5\qquad x_6\qquad x_7$}\\ /\qquad\qquad\qquad\qquad\qquad\qquad\text{\textbackslash}\\ \colorbox{#d7dbd8}{$x_0\qquad x_2\qquad x_4\qquad x_6$}\qquad\colorbox{#d7dbd8}{$x_1\qquad x_3\qquad x_5\qquad x_7$} \\ /\qquad\qquad\text{\textbackslash}\qquad\qquad\qquad\qquad/\qquad\qquad\text{\textbackslash}\\ \colorbox{#d7dbd8}{$x_0\qquad x_4$}\qquad\colorbox{#d7dbd8}{$x_2\qquad x_6$}\qquad\colorbox{#d7dbd8}{$x_1\qquad x_5$}\qquad\colorbox{#d7dbd8}{$x_3\qquad x_7$} \\ /\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\\ \colorbox{#d7dbd8}{$x_0$}\qquad\colorbox{#d7dbd8}{$x_4$}\qquad\colorbox{#d7dbd8}{$x_2$}\qquad\colorbox{#d7dbd8}{$x_6$}\qquad\colorbox{#d7dbd8}{$x_1$}\qquad\colorbox{#d7dbd8}{$x_5$}\qquad\colorbox{#d7dbd8}{$x_3$}\qquad\colorbox{#d7dbd8}{$x_7$} \\ \end{matrix}$$Рис. 1. Дерево входных массивов рекурсивного БПФ для $n = 2^3.$

Зададим индексы элементов в первой и последней строках рисунка 1 в двоичном виде$$\begin{matrix}  x_{000} & x_{001} & x_{010} & x_{011} & x_{100} & x_{101} & x_{110} & x_{111} \\  x_{000} & x_{100} & x_{010} & x_{110} & x_{001} & x_{101} & x_{011} & x_{111} \end{matrix}$$Мы видим, что для $n = 2^3$ элементы в последней строке подвергнуты т. н. битреверсной перестановке и расположены в т. н. битреверсном порядке: двоичная запись индекса элемента в последней строке является перевернутой двоичной записью соответствующего элемента в исходной строке. Аналогичное утверждение верно для $n = 2^k$ при любом $k.$

Зададим теперь дерево выходных массивов алгоритма БПФрекурс. Вспомнив, что рекурсивный алгоритм БПФ записывает выходной массив на место входного массива, мы видим, что дерево выходных массивов можно получить, зеркально отразив дерево входных массивов относительно воображаемой горизонтальной линии (на рисунке 2 изображена пунктиром) и заменив каждый входной массив $a_0,  a_1,  a_2, \ldots, a_{m-1}$ на результат его дискретного преобразования Фурье, который мы будем обозначать $\overline{\overline{\colorbox{#42f5b6}{$a_0  a_1  a_2 \ldots a_{m-1}$}}}:$

$$\begin{matrix} \qquad\qquad\colorbox{#d7dbd8}{$x_0\qquad x_1\qquad x_2\qquad x_3\qquad x_4\qquad x_5\qquad x_6\qquad x_7$}\\ \qquad\qquad/\qquad\qquad\qquad\qquad\qquad\qquad\text{\textbackslash}\\ \qquad\qquad\colorbox{#d7dbd8}{$x_0\qquad x_2\qquad x_4\qquad x_6$}\qquad\colorbox{#d7dbd8}{$x_1\qquad x_3\qquad x_5\qquad x_7$} \\ \qquad\quad/\qquad\qquad\text{\textbackslash}\qquad\quad\qquad\qquad/\qquad\qquad\text{\textbackslash}\\ \qquad\qquad\colorbox{#d7dbd8}{$x_0\qquad x_4$}\qquad\colorbox{#d7dbd8}{$x_2\qquad x_6$}\qquad\colorbox{#d7dbd8}{$x_1\qquad x_5$}\qquad\colorbox{#d7dbd8}{$x_3\qquad x_7$} \\ \qquad\qquad/\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\qquad\qquad/\qquad\text{\textbackslash}\\ \qquad\qquad\colorbox{#d7dbd8}{$x_0$}\qquad\colorbox{#d7dbd8}{$x_4$}\qquad\colorbox{#d7dbd8}{$x_2$}\qquad\colorbox{#d7dbd8}{$x_6$}\qquad\colorbox{#d7dbd8}{$x_1$}\qquad\colorbox{#d7dbd8}{$x_5$}\qquad\colorbox{#d7dbd8}{$x_3$}\qquad\colorbox{#d7dbd8}{$x_7$} \\ \qquad\qquad-------------------------\\ \qquad\qquad\overline{\overline{\colorbox{#42f5b6}{$x_0$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_4$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_2$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_6$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_1$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_5$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_3$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_7$}}} \\ \text{Этап 0}\quad\quad\quad\text{\textbackslash}\qquad/\qquad\qquad\text{\textbackslash}\qquad/\qquad\qquad\text{\textbackslash}\qquad/\qquad\qquad\text{\textbackslash}\qquad/\quad\quad\\ \qquad\qquad\overline{\overline{\colorbox{#42f5b6}{$x_0\qquad x_4$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_2\qquad x_6$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_1\qquad x_5$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_3\qquad x_7$}}} \\ \text{Этап 1}\qquad\qquad\qquad\text{\textbackslash}\qquad\qquad/\qquad\quad\qquad\qquad\text{\textbackslash}\qquad\qquad/\qquad\qquad\qquad\\ \qquad\qquad\overline{\overline{\colorbox{#42f5b6}{$x_0\qquad x_2\qquad x_4\qquad x_6$}}}\qquad\overline{\overline{\colorbox{#42f5b6}{$x_1\qquad x_3\qquad x_5\qquad x_7$}}}\\ \text{Этап 2 }\qquad\qquad\qquad\qquad\text{\textbackslash}\qquad\qquad\qquad\qquad\qquad\qquad/\qquad\qquad\qquad\qquad\\ \qquad\qquad\overline{\overline{\colorbox{#42f5b6}{$x_0\qquad x_1\qquad x_2\qquad x_3\qquad x_4\qquad x_5\qquad x_6\qquad x_7$}}}\\ \end{matrix}$$Рис. 2. Деревья входных и выходных массивов рекурсивного БПФ.

Исходя из рисунка 2, можно построить эффективную нерекурсивную реализацию алгоритма Кули – Тьюки, которая проводит все вычисления без использования дополнительных ресурсов памяти (англ. in place), последовательно модифицируя единственный входно-выходной массив длины $n.$

алг БПФитер(арг цел n, аргрез компл A[0 . . n-1])

дано | n = 2**k и k >= 1

нач цел i, j, чет, нечет, k = log2(n); компл w

  создать компл W[0 . . n/2-1]

  W[0]:=1

  нц для i от 1 до n/2-1

    W[i] := W[i-1]*(cos(-2*pi/n), sin(-2*pi/n))

  кц

  цел s | номер этапа пробегает 0, 1,…, k-1

  цел l=1 | длина подмассива равна 2**s и пробегает 1,2,4,…,n/2

  цел n=n/2 | число пар подмассивов равно n/2**(s+1) и пробегает n/2,n/4,…,4,2,1

  нц для s от 0 до k-1

    нц для j от 0 до n-1

      чет := j*2*l; нечет := чет + l

      нц для i от 0 до l-1

        w := W[n*i]

        | параллельное присваивание:

        A[чет+i], A[нечет+i]:=A[чет+i]+w*A[нечет+i], A[чет+i]-w*A[нечет+i]

      кц

    кц

    n := n/2

    l := 2*l

  кц

кон

Предполагается, что на вход этой программы подаётся массив, содержащий вектор входных данных $x,$ уже переупорядоченный в битреверсном порядке: A[$\overline{b_{k-1}b_{k-2}\ldots b_1b_0}$] = $x_{\overline{b_0b_1\ldots b_{k-1}}},$ где $b_{k-1},b_{k-2},\ldots, b_1,b_0$ – двоичные цифры представления индекса, т. е. $\overline{b_{k-1}b_{k-2}\ldots b_1b_0}=\sum_{i=0}^{k-1}b_i2^i.$ В результате выполнения программы в массиве A размещается вектор $\text{ДПФ}_{n,w_n}(x).$

Алгоритм БПФитер проводит $k$ этапов модификации массива A с номерами от $0$ до $k-1.$ На этапе $s$ элементы массива А разбиваются на $n=2^{k-s-1}$ пар соседствующих подмассивов длины $l=2^s,$ и каждая пара подмассивов сливается in place в один массив удвоенной длины с помощью выполнения $l=2^s$ бабочек Фурье, согласно конструкции, описанной в лемме Ланцоша – Даниельсона. Множители бабочек выбираются из массива W с шагом $n.$ На каждом из $k$ этапов выполняется $n/2$ бабочек Фурье, и общее число арифметических операций в нерекурсивной реализации не превосходит $$\frac n2 + \frac {3n}2\cdot k = \frac n2 + \frac{3n}2\,\log_2n = O(n\,\log_2n).$$

Эффективность компьютерных программ, вычисляющих БПФ

Практиков интересует трудоёмкость и время выполнения ДПФ на конкретных компьютерах для конкретных значений $n.$ Это время зависит не только от используемого в программе алгоритма, но и от аппаратных особенностей компьютера, на котором выполняется программа. Среди свободно распространяемых программ, эффективно выполняющих БПФ, широкую известность получила программа FFTW (от англ. The Fastest Fourier Transform in The West) (Frigo. 1998). Эта программа в начале своего выполнения тестирует особенности компьютера, на котором она выполняется, и выбирает вариант алгоритма, учитывающий эти особенности.

До сих пор не известно, можно ли вычислить классическое ДПФ размера $n,$ выполнив (в обозначениях «о малое») $о(n\,\log_2n)$ арифметических операций.