Алгоритм Кули — Тьюки (Glikjnmb Trln — M,Ztn)

Алгоритм Ку́ли — Тью́ки — наиболее часто используемый алгоритм вычисления быстрого преобразования Фурье. Алгоритм позволяет выразить дискретное преобразование Фурье длины, равной произвольному составному числу $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}$ и производится следующая замена индексов^[1]:

{\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$ . Тогда в терминах двумерных переменных формула преобразуется к виду^[2]

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}$ сводится к выполнению следующих действий:

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

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

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

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

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

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

{\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$ , то говорят об алгоритме Кули — Тьюки с прореживанием по частоте^[4]. В этом случае уравнения преобразуются к виду

{\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))$ ^[5]. Аналогичные формулы можно получить с использованием вещественных констант $1/(2\cos \left(2\pi k/N\right))$ ^[6].

История

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

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

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

Примечания

↑ Блейхут, 1989, с. 129.
↑ ¹ ² Блейхут, 1989, с. 130.
↑ ¹ ² ³ Блейхут, 1989, с. 133.
↑ Блейхут, 1989, с. 134.
↑ Блейхут, 1989, с. 138.
↑ Нуссбаумер, 1985, с. 92.
↑ Heideman, M. T., Johnson, D. H., Burrus, C. S.. Gauss and the history of the fast Fourier transform (англ.) // IEEE ASSP Magazine. — IEEE, 1984. — Vol. 1, iss. 4. — P. 14—21. — doi:10.1109/MASSP.1984.1162257. Архивировано 21 сентября 2023 года.
↑ Cooley, J. W., Tukey, J. W. An algorithm for the machine calculation of complex Fourier series (англ.) // Mathematics of Computation. — 1965. — Vol. 19. — P. 297—301. — doi:10.1090/S0025-5718-1965-0178586-1. Архивировано 7 мая 2021 года.
↑ Danielson, G. C., Lanczos, C. Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids (англ.) // Journal of the Franklin Institute. — 1942. — Vol. 233, iss. 4. — P. 365–380. — doi:10.1016/S0016-0032(42)90767-1.

Литература

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

См. также

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

[_dfbb5a8e6456d869-1] Блейхут, 1989, с. 129.

[_dfbb5b8e6456da33-2] ¹ ² Блейхут, 1989, с. 130.

[_dfbb5b8e6456da30-3] ¹ ² ³ Блейхут, 1989, с. 133.

[_dfbb5b8e6456da37-4] Блейхут, 1989, с. 134.

[_dfbb5b8e6456da3b-5] Блейхут, 1989, с. 138.

[_8f328b600a9b7875-6] Нуссбаумер, 1985, с. 92.

[7] Heideman, M. T., Johnson, D. H., Burrus, C. S.. Gauss and the history of the fast Fourier transform (англ.) // IEEE ASSP Magazine. — IEEE, 1984. — Vol. 1, iss. 4. — P. 14—21. — doi:10.1109/MASSP.1984.1162257. Архивировано 21 сентября 2023 года.

[8] Cooley, J. W., Tukey, J. W. An algorithm for the machine calculation of complex Fourier series (англ.) // Mathematics of Computation. — 1965. — Vol. 19. — P. 297—301. — doi:10.1090/S0025-5718-1965-0178586-1. Архивировано 7 мая 2021 года.

[9] Danielson, G. C., Lanczos, C. Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids (англ.) // Journal of the Franklin Institute. — 1942. — Vol. 233, iss. 4. — P. 365–380. — doi:10.1016/S0016-0032(42)90767-1.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]