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

Быстрое преобразование Фурье (БПФ, FFT) — алгоритм ускоренного вычисления дискретного преобразования Фурье, позволяющий получить результат за время, меньшее чем [math]\displaystyle{ O(N^2) }[/math] (требуемого для прямого, поформульного вычисления). Иногда под быстрым преобразованием Фурье понимается один из алгоритмов, называемый алгоритмом прореживания по частоте — времени, имеющий сложность [math]\displaystyle{ O(N\log(N)) }[/math].

Алгоритм применим к любым коммутативным ассоциативным кольцам с единицей, чаще применяют к полю комплексных чисел (c [math]\displaystyle{ \varepsilon=e^{2\pi i/n} }[/math]) и к кольцам вычетов (которым, в частности, является компьютерный целый тип).

Основной алгоритм

Основная статья: Алгоритм Кули — Тьюки

Амплитуды быстрого преобразования Фурье для разного количества компонент разложения [math]\displaystyle{ N }[/math]. Первый случай: [math]\displaystyle{ N = L-10 }[/math], где [math]\displaystyle{ L }[/math] — количество отсчётов сигнала; второй случай: [math]\displaystyle{ N = L }[/math]; третий: [math]\displaystyle{ N = L+10 }[/math]. В первом и последних случаях спектральные характеристики оцениваются менее точно. Данный эффект связан с явлением Гиббса.

При применении основного алгоритма дискретное преобразование Фурье может быть выполнено за [math]\displaystyle{ O(N(p_1+\cdots+p_n)) }[/math] действий при [math]\displaystyle{ N=p_1p_2\cdots p_n }[/math], в частности, при [math]\displaystyle{ N=2^n }[/math] понадобится [math]\displaystyle{ O(N\log(N)) }[/math] действий.

Дискретное преобразование Фурье преобразует набор чисел [math]\displaystyle{ a_0, \dots, a_{n-1} }[/math] в набор чисел [math]\displaystyle{ b_0, \dots, b_{n-1} }[/math], такой, что [math]\displaystyle{ b_i=\sum_{j=0}^{n-1}a_j\varepsilon^{ij} }[/math], где [math]\displaystyle{ \varepsilon }[/math] — первообразный корень из единицы, то есть [math]\displaystyle{ \varepsilon^n=1 }[/math] и [math]\displaystyle{ \varepsilon^k\neq 1 }[/math] при [math]\displaystyle{ 0\lt k\lt n }[/math]. Основной шаг алгоритма состоит в сведении задачи для [math]\displaystyle{ N }[/math] чисел к задаче с меньшим числом. Для [math]\displaystyle{ N=pq, p\gt 1, q\gt 1 }[/math] над полем комплексных чисел вводятся: [math]\displaystyle{ \varepsilon_\nu=e^{2\pi i/\nu} }[/math], [math]\displaystyle{ \varepsilon_\nu^{\nu}=1 }[/math], где [math]\displaystyle{ \nu }[/math] — любое число. Дискретное преобразование Фурье может быть представлено в виде [math]\displaystyle{ b_i=\sum_{k=0}^{p-1}\sum_{j=0}^{q-1}a_{kq+j}\varepsilon_N^{(kq+j)i} }[/math]. (Эти выражения могут быть легко получены, если исходную сумму разбить на меньшее число сумм с меньшим числом слагаемых, а после полученные суммы привести к одинаковому виду путём сдвига индексов и их последующего переобозначения).

Таким образом:

[math]\displaystyle{ b_i = \sum_{k=0}^{p-1}\sum_{j=0}^{q-1}a_{kq+j}\varepsilon_N^{(kq+j)i}=\sum_{j=0}^{q-1} \varepsilon_N^{ij} \left( \sum_{k=0}^{p-1}a_{kq+j}\varepsilon_N^{kiq} \right) }[/math].

С учётом того, что [math]\displaystyle{ \varepsilon_N^{kiq}=\varepsilon_{N/q}^{ki} }[/math] и [math]\displaystyle{ N/q=p }[/math], окончательная запись:

[math]\displaystyle{ b_i=\sum_{j=0}^{q-1} \varepsilon_N^{ij} \left( \sum_{k=0}^{p-1}a_{kq+j}\varepsilon_p^{ki} \right) }[/math].

Далее вычисляется каждое [math]\displaystyle{ b_i }[/math], где [math]\displaystyle{ i=\overline{0,p-1} }[/math], здесь по-прежнему требуется совершить [math]\displaystyle{ O(N) }[/math] действий, то есть на этом этапе производится [math]\displaystyle{ p\cdot O(N)=O(Np) }[/math] операций.

Далее считается [math]\displaystyle{ b_{i+mp} }[/math], где [math]\displaystyle{ i=\overline{0,p-1}, m=\overline{1,q-1} }[/math]. При замене [math]\displaystyle{ i\longmapsto i+mp }[/math] в последней формуле, выражения, стоящие в скобках, остались неизменными, а так как они уже были посчитаны на предыдущем шаге, то на вычисление каждого из них потребуется только [math]\displaystyle{ O(q) }[/math] действий. Всего [math]\displaystyle{ p(q-1)=N-p }[/math] чисел. Следовательно, операций на этом шаге [math]\displaystyle{ (N-p)\cdot O(q)=O((N-1)q)\cong O(Nq) }[/math]. Последнее с хорошей точностью верно при любых [math]\displaystyle{ N }[/math].

Алгоритм быстрого преобразования Фурье логично применять для [math]\displaystyle{ N\gg1 }[/math], потому как при малом числе отсчётов он даёт небольшой выигрыш в скорости по отношению к прямому расчёту дискретного преобразования Фурье. Таким образом, для того чтобы полностью перейти к набору чисел [math]\displaystyle{ b_0, \dots, b_{N-1} }[/math], необходимо [math]\displaystyle{ O(Np)+O(Nq) }[/math] действий. Следовательно, нет разницы, на какие два числа разбивать [math]\displaystyle{ N }[/math] — ответ от этого сильно не будет меняться. Уменьшено же число операций может быть только дальнейшим разбиением [math]\displaystyle{ N }[/math].

Скорость алгоритма [math]\displaystyle{ (N=pq) }[/math]:

1. [math]\displaystyle{ p\gg q \longrightarrow O(Np) }[/math]
2. [math]\displaystyle{ p\sim q \longrightarrow O(Nq) }[/math]
3. [math]\displaystyle{ p\ll q \longrightarrow O(Nq) }[/math]

То есть число операций при любом разбиении [math]\displaystyle{ N }[/math] на два числа, есть [math]\displaystyle{ O(Nc) }[/math], где [math]\displaystyle{ c=max(p,q) }[/math].

Обратное преобразование Фурье

Для обратного преобразования Фурье можно применять алгоритм прямого преобразования Фурье — нужно лишь использовать [math]\displaystyle{ \varepsilon^{-1} }[/math] вместо [math]\displaystyle{ \varepsilon }[/math] (или применить операцию комплексного сопряжения вначале к входным данным, а затем к результату, полученному после прямого преобразования Фурье), и окончательный результат поделить на [math]\displaystyle{ N }[/math].

Общий случай

Основная статья: Алгоритм Блуштайна

Общий случай может быть сведён к предыдущему. Для [math]\displaystyle{ 4N\gt 2^k\ge2N }[/math] имеет место:

[math]\displaystyle{ b_i=\varepsilon^{-i^2/2}\sum_{j=0}^{N-1}\varepsilon^{(i+j)^2/2}\varepsilon^{-j^2/2}a_j }[/math].

Обозначая [math]\displaystyle{ \bar{a}_i=\varepsilon^{-i^2/2}a_i, \bar{b}_i=\varepsilon^{i^2/2}b_i, c_i=\varepsilon^{(2N-2-i)^2/2} }[/math] получается:

[math]\displaystyle{ \bar{b}_i=\sum_{j=0}^{2N-2-i}\bar{a}_jc_{2N-2-i-j} }[/math],

если положить [math]\displaystyle{ \bar{a}_i=0 }[/math] при [math]\displaystyle{ i\ge N }[/math].

Таким образом задача сведена к вычислению свёртки, но это можно сделать с помощью трёх преобразований Фурье для [math]\displaystyle{ 2^k }[/math] элементов: выполняется прямое преобразование Фурье для [math]\displaystyle{ \{\bar{a}_i\}_{i=0}^{i=2^k-1} }[/math] и [math]\displaystyle{ \{c_i\}_{i=0}^{i=2^k-1} }[/math], поэлементно перемножаются результаты, и выполняется обратное преобразование Фурье.

Вычисления всех [math]\displaystyle{ \bar{a}_i }[/math] и [math]\displaystyle{ c_i }[/math] требуют [math]\displaystyle{ O(N) }[/math] действий, три преобразования Фурье требуют [math]\displaystyle{ O(N\log(N)) }[/math] действий, перемножение результатов преобразований Фурье требует [math]\displaystyle{ O(N) }[/math] действий, вычисление всех [math]\displaystyle{ b_i }[/math] зная значения свёртки требует [math]\displaystyle{ O(N) }[/math] действий. Итого для дискретного преобразования Фурье требуется [math]\displaystyle{ O(N\log(N)) }[/math] действий для любого [math]\displaystyle{ N }[/math].

Этот алгоритм быстрого преобразования Фурье может работать над кольцом только когда известны первообразные корни из единицы степеней [math]\displaystyle{ 2N }[/math] и [math]\displaystyle{ 2^k }[/math].

Вывод преобразования из дискретного

Наиболее распространённым алгоритмом быстрого преобразования Фурье является алгоритм Кули — Тьюки, при котором дискретное преобразование Фурье от [math]\displaystyle{ N=N_1 N_2 }[/math] выражается как сумма дискретных преобразований Фурье более малых размерностей [math]\displaystyle{ N_1 }[/math] и [math]\displaystyle{ N_2 }[/math] рекурсивно для того, чтобы достичь сложность [math]\displaystyle{ O(N\log(N)) }[/math]. Элементарный шаг сочленения меньших преобразований Фурье в этом алгоритме называется «бабочка». В вычислительной технике наиболее часто используется рекурсивное разложение преобразований надвое, то есть с основанием 2 (хотя может быть использовано любое основание), а количество входных отсчётов является степенью двойки. Для случаев, когда дискретное преобразование считается от количества отсчётов, являющихся простыми числами, могут быть использованы алгоритмы Блуштайна и Рейдера.

Например, для вычисления быстрого преобразования Фурье по алгоритму Кули — Тьюки с основанием 2 для вектора [math]\displaystyle{ \vec x }[/math], состоящего из [math]\displaystyle{ N }[/math] элементов:

[math]\displaystyle{ \vec X = \hat A \vec x }[/math],

с элементами [math]\displaystyle{ \hat A }[/math] вида:

[math]\displaystyle{ a_{N}^{mn} = e^ { -\frac{2\pi i}{N} mn} }[/math]

дискретное преобразование можно выразить как сумму двух частей: сумму чётных индексов [math]\displaystyle{ m={2n} }[/math] и сумму нечётных индексов [math]\displaystyle{ m={2n+1} }[/math]:

[math]\displaystyle{ X_m=\sum_{n=0}^{N-1} x_n a_{N}^{mn} = \sum_{n=0}^{N/2-1}x_{2n} a_{N}^{2nm} + \sum_{n=0}^{N/2-1}x_{2n+1} a_{N}^{(2n+1)m} }[/math].

Коэффициенты [math]\displaystyle{ a_{N}^{2nm} }[/math] и [math]\displaystyle{ a_{N}^{(2n+1)m} }[/math] можно переписать следующим образом:

[math]\displaystyle{ a_{N}^{2nm} = e^\left( -2\pi i \frac{2mn}{N} \right) = e^ \left( -2\pi i \frac{mn}{N/2} \right) = a_{N/2}^{nm} }[/math]

[math]\displaystyle{ a_{N}^{(2n+1)m} = e^ \left( -2\pi i \frac{m}{N} \right) a_{N/2}^{nm} }[/math]

В результате:

[math]\displaystyle{ X_m=\sum_{n=0}^{N/2-1} x_{2n} a_{N/2}^{nm} + e^ { -\frac{2\pi i}{N} m} \sum_{n=0}^{N/2-1} x_{2n+1} a_{N/2}^{nm} }[/math]

Вычисление данного выражения можно упростить, используя:

• свойство периодичности ДПФ:

  [math]\displaystyle{ a_{N/2}^{(m+\frac{N}{2})n} = a_{N/2}^{nm} }[/math],

• коэффициент поворота БПФ удовлетворяет следующему равенству:

  [math]\displaystyle{ \begin{matrix} e^{\frac{-2\pi i}{N} (m + N/2)} & = & e^{\frac{-2\pi i m}{N} - {\pi i}} \\ & = & e^{-\pi i} e^{\frac{-2\pi i m}{N}} \\ & = & -e^{\frac{-2\pi i m}{N}} \end{matrix} }[/math]

В результате упрощений, обозначив дискретное преобразование Фурье чётных индексов [math]\displaystyle{ x_{2m} }[/math] через [math]\displaystyle{ E_m }[/math] и преобразование нечётных индексов [math]\displaystyle{ x_{2m + 1} }[/math] через [math]\displaystyle{ O_m }[/math], для [math]\displaystyle{ 0 \leqslant m \lt \frac{N}{2} }[/math] получается:

[math]\displaystyle{ \begin{matrix} X_m & = & E_m + e^{-\frac{2\pi i}{N}m} O_m \\ X_{m+\frac{N}{2}} & = & E_m - e^{-\frac{2\pi i}{N}m} O_m \end{matrix} }[/math]

Данная запись является базой алгоритма Кули — Тьюки с основанием 2 для вычисления быстрого преобразования Фурье, то есть дискретное преобразование от вектора, состоящего из [math]\displaystyle{ N }[/math] отсчётов, приведено к линейной композиции двух дискретных преобразований Фурье от [math]\displaystyle{ \frac N2 }[/math] отсчётов, и, если для первоначальной задачи требовалось [math]\displaystyle{ N^2 }[/math] операций, то для полученной композиции — [math]\displaystyle{ \frac{N^2}{2} }[/math] (за счёт повторного использования промежуточных результатов вычислений [math]\displaystyle{ E_m }[/math] и [math]\displaystyle{ O_m }[/math]). Если [math]\displaystyle{ N }[/math] является степенью двух, то это разделение можно продолжать рекурсивно до тех пор, пока не доходит до двухточечного преобразования Фурье, которое вычисляется по следующим формулам:

[math]\displaystyle{ \begin{cases} X_0=x_0+x_1\\ X_1=x_0-x_1 \end{cases} }[/math]

При рекурсивном делении дискретного преобразования Фурье от [math]\displaystyle{ N }[/math] входных значений на сумму 2 дискретных преобразований по [math]\displaystyle{ N/2 }[/math] входных значений сложность алгоритма становится равной [math]\displaystyle{ O(N\log(N)) }[/math].

Ссылки

• Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления свёрток. — М.: Радио и связь, 1985.