Алгоритм Кули — Тьюки

Алгоритм Ку́ли — Тью́ки — наиболее часто используемый алгоритм вычисления быстрого преобразования Фурье. Алгоритм позволяет выразить дискретное преобразование Фурье длины, равной произвольному составному числу $N$, через определённое количество преобразований меньшей длины с помощью рекурсии, понижая таким образом сложность вычислений до $O(N\log N)$ для гладких $N$. Назван в честь Дж. Кули и Дж. Тьюки.

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

Схема общего алгоритма Кули — Тьюки

Для произвольного натурального числа $N$ дискретным преобразованием Фурье комплексного вектора $x(j)$, где $j=0,\ldots ,N-1$, называется комплексный вектор $X(k)$, где $k=0,\ldots ,N-1$, компоненты которого задаются формулой

$X(k)=\sum \limits _{j=0}^{N-1}x(j)\omega ^{kj},\;k=0,\ldots ,N-1,$

где ${\textstyle \omega =e^{-{\frac {2\pi i}{N}}}}$.

Для построения алгоритма делается предположение, что $N=N_{1}N_{2}$ для некоторых натуральных $N_{1}$ и $N_{2}$ и производится следующая замена индексов:

${\begin{array}{lll}&j=j_{1}+N_{1}j_{2},\ &j_{1}=0,\ldots ,N_{1}-1,\ &j_{2}=0,\ldots ,N_{2}-1,\\&k=k_{2}+N_{2}k_{1},\ &k_{1}=0,\ldots ,N_{1}-1,\ &k_{2}=0,\ldots ,N_{2}-1,\end{array}}$

в результате которой получается

$X(k_{2}+N_{2}k_{1})=\sum \limits _{j_{1}=0}^{N_{1}-1}\sum \limits _{j_{2}=0}^{N_{2}-1}x(j_{1}+N_{1}j_{2})\omega ^{(j_{1}+N_{1}j_{2})(k_{2}+N_{2}k_{1})}.$

Далее векторы входных и выходных данных преобразуются в двумерные массивы $x'$ и $X'$, задающиеся равенствами

${\begin{array}{lll}&x'(j_{1},j_{2})=x(j_{1}+N_{1}j_{2}),\ &j_{1}=0,\ldots ,N_{1}-1,\ &j_{2}=0,\ldots ,N_{2}-1,\\&X'(k_{1},k_{2})=X(k_{2}+N_{2}k_{1}),\ &k_{1}=0,\ldots ,N_{1}-1,\ &k_{2}=0,\ldots ,N_{2}-1.\end{array}}$

Стоит заметить, что компоненты $X$ упорядочены не так, как компоненты $x$.

Далее пусть $\omega ^{N_{1}}=\gamma$ и $\omega ^{N_{2}}=\beta$. Очевидно, что $\omega ^{N_{1}N_{2}j_{2}k_{1}}=1$. Тогда в терминах двумерных переменных формула преобразуется к виду

$X'(k_{1},k_{2})=\sum \limits _{j_{1}=0}^{N_{1}-1}\beta ^{j_{1}k_{1}}\left(\omega ^{j_{1}k_{2}}\sum \limits _{j_{2}=0}^{N_{2}-1}x'(j_{1},j_{2})\gamma ^{j_{2}k_{2}}\right).$

Таким образом, вычисление преобразования длины $N=N_{1}N_{2}$ сводится к выполнению следующих действий:

1. Вычисление $N_{1}$ преобразований длины $N_{2}$.
2. Умножение на комплексные «поворачивающие» множители.
3. Вычисление $N_{2}$ преобразований длины $N_{1}$.

При этом вместо $N^{2}$ комплексных сложений и $N^{2}$ комплексных умножений исходной формулы итоговая схема содержит не более $N(N_{1}+N_{2}-2)$ комплексных сложений и $N(N_{1}+N_{2}+1)$ комплексных умножений.

Каждое из преобразований длины $N_{1}$ или $N_{2}$ можно вычислять с помощью различных быстрых алгоритмов, в частности, снова по вышеприведённой схеме. В этом случае преобразование длины ${\textstyle N=\prod \limits _{i=1}^{n}N_{i}}$ может быть представлено в форме, требующей выполнения ${\textstyle N\sum \limits _{i=1}^{n}N_{i}}$ комплексных умножений.

Алгоритм по основанию два

Во многих приложениях длина преобразования равна степени двойки: $N=2^{m}$. Тогда в вышеприведённой схеме возможны варианты $N_{1}=2$ или $N_{2}=2$. В этом случае говорят об алгоритме Кули — Тьюки по основанию два (англ. radix-2).

Если $N_{1}=2$, то говорят об алгоритме Кули — Тьюки с прореживанием по времени. В этом случае уравнения преобразуются к виду

${\begin{array}{rcl}X(k)&=&\sum \limits _{j=0}^{N/2-1}x(2j)\omega ^{2kj}+\omega ^{k}\sum \limits _{j=0}^{N/2-1}x(2j+1)\omega ^{2kj}\\X(k+N/2)&=&\sum \limits _{j=0}^{N/2-1}x(2j)\omega ^{2kj}-\omega ^{k}\sum \limits _{j=0}^{N/2-1}x(2j+1)\omega ^{2kj},\ k=0,\ldots ,N/2-1.\end{array}}$

Схема реализации операции «бабочки»

Если ввести обозначения $E(k)=\sum \limits _{j=0}^{N/2-1}x(2j)\omega ^{2kj}$ и $O(k)=\sum \limits _{j=0}^{N/2-1}x(2j+1)\omega ^{2kj}$, то уравнения можно переписать как

${\begin{array}{rcl}X(k)&=&E(k)+\omega ^{k}O(k)\\X(k+N/2)&=&E(k)-\omega ^{k}O(k),\ k=0,\ldots ,N/2-1.\end{array}}$

Такая операция обычно называется «бабочкой».

Данная процедура может быть применена к входному вектору рекурсивно. На каждом шаге $N$-точечное преобразование разбивается на два $(N/2)$-точечных преобразования, которые, в свою очередь, разбиваются таким же образом до тех пор, пока длина преобразования не станет равна единице. Затем происходит обратный ход, на каждом шаге длины результатов преобразований удваиваются с помощью «бабочек». При такой реализации выполняется $(N/2)\log _{2}N$ комплексных умножений и $N\log _{2}N$ комплексных сложений.

Этот процесс является примером применения методики «разделяй и властвуй». При этом во многих реализациях прямой рекурсии избегают и вместо неё дерево вычислений проходится в порядке поиска в ширину.

Если $N_{2}=2$, то говорят об алгоритме Кули — Тьюки с прореживанием по частоте. В этом случае уравнения преобразуются к виду

${\begin{array}{rcl}X(2k)&=&\sum \limits _{j=0}^{N/2-1}{\bigg (}x(j)+x(j+N/2){\bigg )}\omega ^{2kj},\\X(2k+1)&=&\sum \limits _{j=0}^{N/2-1}{\bigg (}x(j)-x(j+N/2){\bigg )}\omega ^{j(2k+1)},\ k=0,\ldots ,N/2-1.\end{array}}$

Алгоритм Рейдера — Бреннера

Пусть

$a(k)=\left\lbrace {\begin{array}{rl}0,&k=0,\\{\dfrac {x(k)-x(k+N/2)}{2i\sin \left(2\pi k/N\right)}},&k=1,\ldots ,N/2-1\end{array}}\right.$

и пусть

$A(k)=\sum \limits _{j=0}^{N/2-1}a(j)\omega ^{2kj},\ k=0,\ldots ,N/2-1.$

С использованием формул алгоритма с прореживанием по частоте нетрудно убедиться, что выполняется следующее соотношение:

$X(2k+1)=A(k)-A(k+1)+(x(0)-x(N/2)),\ k=0,\ldots ,N/2-1.$

Такая модификация алгоритма по основанию два называется алгоритмом Рейдера — Бреннера. Она позволяет уменьшить вычислительную сложность за счёт более простых умножений на чисто мнимые константы $1/(2i\sin \left(2\pi k/N\right))$. Аналогичные формулы можно получить с использованием вещественных констант $1/(2\cos \left(2\pi k/N\right))$.

История

Алгоритм и его рекурсивная реализация были изобретены около 1805 года К. Гауссом при интерполировании траекторий астероидов Паллада и Юнона. Тогда открытие не получило широкого распространения и было опубликовано лишь после смерти учёного на новой латыни. Результат Гаусса несколько раз переоткрывался в различных формах в течение последующих 150 лет и стал популярным после публикации в 1965 году статьи Дж. Кули из IBM и Дж. Тьюки из Принстона, в которой алгоритм был в очередной раз переоткрыт, а также описывалась удобная реализация для ЭВМ.

Тот факт, что первооткрывателем алгоритма является Гаусс, был обнаружен лишь через несколько лет после публикации Кули и Тьюки. В своей статье они ссылались только на работу И. Дж. Гуда, в которой был описан алгоритм Гуда — Томаса.

Выражение преобразования Фурье длины $N$ через два преобразования длины $N/2$ иногда называют леммой Дэниэльсона — Ланцоша, так как оно было получено этими двумя авторами в 1942 году.

Литература

• Блейхут, Р.. Быстрые алгоритмы цифровой обработки сигналов (рус.). — М.: Мир, 1989. — 448 с.
• Нуссбаумер, Г.. Быстрое преобразование Фурье и алгоритмы вычисления свёрток. — М.: Радио и связь, 1985.
• Рабинер, Л., Гоулд, Б.. Теория и применение цифровой обработки сигналов. — М.: Мир, 1978.

См. также

• Алгоритм Блуштайна