степени меньше га в корнях степени га из единицы, то есть в точках ийп,и\,и2п, (Напомним, что для умножения двух многочленов степени меньше т мы использовали их значения в 2т точках, так что нынешнее значение га удвоено по сравнению с разделом 32.1). Мы предполагаем, что га является степенью 2 (этого для наших целей достаточно, так как к многочлену всегда можно добавить нулевые старшие коэффициенты). Итак, нам задан вектор коэффициентов а = (ао, а\,... , ara-i) и надо вычислить

п-1

у* = А(а;*) = J>rf.(32.8)

з=о

для к = 0,1,... , га - 1.

Вектор у = (уо, у\,... , yn-i) называется дискретным преобразованием Фурье (Discrete Fourier Transform, DFT) вектора a = (ао, a\,... , an i). Обозначение: у = DFTra(a).

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

(Fast Fourier Transform, FFT) представляет собой метод быстрого вычисления дискретного преобразования Фурье, использующий свойства комплексных корней из единицы и требующий времени ©(ralgra) (а не 0(га2), как получится при использовании формул (32.8)).

Быстрое преобразование Фурье использует метод «разделяй и властвуй». Выделим в многочлене А отдельно члены чётных и нечётных степеней, записав

А(х) = А(х2) + хА(х2),(32.9)

А И (ж) = ао + а2х + а4х2 + ... + ап 2хП2~1,

где

А(ж) = а\ + а3х + а5Ж2 + ... + ап \хп12 1. Тем самым задача вычисления А(ж) в точках и°а, ип, ..., и.

п-1

сводится к

1.вычислению значений многочленов А И и А степени меньше га/2 в точках

(с°)2,()2,...,К-1)2;(32.10)

2.комбинации результатов по формуле (32.9).

Лемма о делении пополам гарантирует, что список (32.10) содержит всего га/2 различных чисел, а именно, комплексных корней степени га/2 из единицы (каждый входит дважды). Таким образом, нам нужно вычислить значения многочленов А И и А (степени меньше га/2) в га/2 комплексных корнях степени га/2 из единицы. Эти подзадачи имеют тот же вид, что исходная, но вдвое меньший

размер. Мы приходим к такому рекурсивному алгоритму: вычисления преобразования Фурье вектора а = (ао, а\,... , ara-i) (гДе п - степень 2):

\textsc{Recursive-FFT}$(a)$У, Рекурсивное БПФ

1$n\leftarrow length[а]$ \qquad $n$ --- степень $2$

2if $n=l$

3\quad then return $a$

4$\omega n\leftarrow e~{2\pi i/n}$

5$\omega\leftarrow 1$

6$a~{[0]}\leftarrow (a 0,a 2,\ldots,a {n-2>)$

7$a"{[l]}\leftarrow (a l,a 3,\ldots,a {n-l>)$

8$y~{[0]}\leftarrow$ \textsc{Recursive-FFT}$(a~{[0]})$

9$y~{[l]}\leftarrow$ \textsc{Recursive-FFT}$(a~{[1]>)$

10for $k\leftarrowO$ to $n/2-l$

11\quad do $y k\leftarrow y k~{[0]} + \omega y k~{[l]}$

12\qquad $y {k+(n/2)>\leftarrow y k~{[0]} - \omega y k~{[l]}$

13\qquad $\omega\leftarrow\omega\omega n$

14return $y$

Процедура Recursive-FFT работает следующим образом. Строки 2-3 образуют «базис рекурсии»: дискретным преобразованием Фурье для вектора длины 1 является сам этот вектор, так как

уо = а0° = а01 = а0.

В строках 6-7 формируются вектора коэффициентов многочленов А И иСтроки 4, 5 и 13 гарантируют, что и = ojk в момент

выполнения строк 11-12 (мы экономим время, не вычисляя значение и>кп каждый раз заново) В строках 8-9 рекурсивно вычисляются значения

yf = An/2), или (поскольку и>п/2 = ипк по лемме о сокращении)

уГ] = аМ(<0, у[ц = аЩ).

В строках 11-12 собираются вместе результаты рекурсивных вычислений DFTra/2- В строке 11 для уо, у\, , уп/2-1 получается

ук = yf + Jj? = A[0]Uk) += AU);

последнее равенство следует из (32.9). Для уп/2, Уп/2+i, , Уп-1 строка 12 даёт (при к = О,1,... , га/2 - 1)

Ук+(п/2)

[0]

У к

= уГ]++(п/2Ы1]

= A\ulk)+oJkn+(nl2)A[l\ojlk) = А(+("/2)).

Второе равенство верно, поскольку С0п+п2 = -ш. Четвёртое равенство верно, поскольку из cj™ = 1 следует, что со2к = со2к+п. Последнее равенство следует из уравнения (32.9). Таким образом, вектор у, возвращаемый процедурой Recursive-FFT, действительно есть дискретное преобразование Фурье входного вектора а.

Сколько времени занимает процедура Recursive-FFT? Если не учитывать рекурсивные вызовы, она требует времени в (га), где га - длина входного вектора. С учётом рекурсивных вызовов получаем такое соотношение для времени Т(п) работы процедуры

Г(га) = 2Г(га/2) + 6(га) = 6(ralgra).

Таким образом, приведённый алгоритм (рекурсивный вариант быстрого преобразования Фурье) позволяет вычислить значения многочлена степени меньше га в комплексных корнях степени га из единицы за время в (га lgra).

Интерполяция по значениям в корнях из единицы

Нам осталось показать, как перейти обратно от значений многочлена в комплексных корнях из единицы к его коэффициентам. Для этого мы представим преобразование Фурье как умножение на матрицу, и найдём обратную матрицу.

Согласно уравнению (32.4), дискретное преобразование Фурье можно записать как матричное умножение у = Vna, где Vn - это матрица Вандермонда, составленная из степеней ujn:

(

Уо У\

У2

Уз

\

\ Уп-1 /

LO,

п-1 2(гг-1) 3(гг-1)

LO,

LO.„

LO,

a0 a3

V an i /

Элемент матрицы Vn с индексами (k,j) равен LokJ (при j,k = 0,1,... , га - 1); показатели степеней в матрице Vn образуют «таблицу умножения».