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

Сам алгоритм практически  не изменился за исключением того, что все расстояния между элементами массива увеличились в M раз, согласно формуле (38). Приставка "M" перед  именами новых переменных означает умножение на M.

На вход этой функции будет передаваться массив со смещением. Тем самым происходит коррекция индекса на +h элементов, согласно формуле (38).

В основной функции  сначала анализируется параметр N. Если он нечетный, то мы вынуждены  отсечь один элемент (последний) и уменьшить N на единицу. Таким образом, на входе алгоритма требуется N ≥ 2.

Далее вычисляется L как максимальная степень двойки, кратная N, сама эта степень T и величина M.

Затем вычисляются  поворачивающие множители для "fft_step" и вызывается сам "fft_step" нужное число раз. Если оказывается, что N является степенью двойки, то сработает обычный "fft".

Следующим шагом  мы вычисляем поворачивающие коэффициенты уже для формулы (40). Этот фрагмент алгоритма очень похож на createWstore, поэтому в дополнительных комментариях не нуждается. Результат оказывается в массиве mult.

Далее вычисляется  формула (40). Смысл переменных: pX - указатель  на X_r+sL; rpsL = r + sL; mprM = m + rM; one - очередное слагаемое суммы (40).

Величина widx равна  остатку от деления m(r + sL) на N. Таким способом мы избегаем выхода за границы массива mult. Это можно делать, благодаря Теореме 1 - ведь поворачивающие множители как раз имеют период N.

Еще можно заметить, что в сумме (40) в первом множителе  поворачивающий коэффициент всегда равен 1, благодаря чему мы экономим несколько умножений.

И наконец, происходит корректировка результатов для  обратного ДПФ делением на N.

Функция БПФ  для четного N называется "fft2".

Пусть N = PQ. Этот представляет собой более общий  случай, чем рассмотренный в предыдущей главе. Теперь не требуется, чтобы P или Q были степенями двойки.

Напомним основную формулу:

    (1).

Рассмотрим формулу:

    (41).

Эта формула  эквивалентна формуле (1). В самом  деле:

Двойная сумма  по pQ + q - это все равно, что одинарная сумма по n, просто порядок слагаемых другой.

Вернемся к  формуле (41). В скобках там стоит  обычное прямое БПФ, только не для  всех x_n, а для последовательности из P штук, взятых с шагом Q. Количество последовательностей равно Q. Элементы последовательности выбираются формулой x_{pQ + q}, где p - номер элемента, q - номер последовательности.

Пусть мы умеем  вычислять БПФ для P элементов  за время Z(P). Нам понадобится вычислить Q таких БПФ, на что потратим время O(QZ(P)). Тем самым мы получим все множители в скобках из формулы (41). После этого надо будет вычислить N элементов X_k, каждый из которых потребует O(Q) действий, итого потребуется еще O(NQ) действий, а в сумме - O(QZ(P) + QN) действий.

В предыдущем алгоритме P - было степенью двойки, а для такого числа элементов мы умеем вычислять БПФ за O(Plog₂P) шагов, т.е. Z(P) = Plog₂P, и получится O(QZ(P) + QN) = O(QPlog₂P + QN) = O(Nlog₂P + N²/P) действий. Тем самым, мы несколько более простым путем получим ту же оценку сложности алгоритма.

Если P - не степень двойки, мы умеем вычислять ПФ за O(P²) шагов, т.е. Z(P) = P², и получится O(QZ(P) + QN) = O(QP² + QN) = O(N²/Q + N²/P) действий. Если P и Q достаточно велики, это лучше, чем O(N²).

В "Википедии" нашелся алгоритм сложности O(Nlog₂N)для произвольного N. Алгоритм основан на значительном увеличении числа элементов последовательности и дополнении ее нулями. Хотя такое расширение, естественно, будет компенсироваться высокой скоростью алгоритма.

Длина последовательности берется равной N', где N' - степень двойки ближайшая сверху к 2N. То есть, 4N > N' = 2^T ≥ 2N. Т.е. число элементов новой последовательности в худшем случае может превышать число элементов исходной последовательности почти в 4 раза. Например, если N = 9, то N' = 32.

Обозначим

Представим формулу (1) в виде:

Преобразуем степень:

kn = 2kn/2 = (−k² + k² + 2kn + n² − n²)/2 = (−k² + (k + n)² − n²)/2

Это используем для дальнейшего преобразования формулы (1):

Далее введем обозначения:

x'_n = x_nW ^−n2/2

X'_k = X_kW ^k2/2

Для n > N, полагаем x_n = 0.

Тогда

   (42)

И далее авторы алгоритма выполняют виртуозный "финт ушами". Пусть:

z_k = (2N − 2 − k)

w_n = W ^z_n^2/2

И превращают формулу (42) в:

   (43)

Проверим:

w_{zk − n} = W ^(z_{zk − n}^)2/2 = W ^{(2N − 2 − (z}_k ^{− n))2/2} = W ^{(2N − 2 − z}_k ^{+ n)2/2} = W ^{(2N − 2 − (2N − 2 − k) + n)2/2} = W ^{(2N − 2 − 2N + 2 + k + n)2/2} = W ^{(k + n)2/2}

Правая часть  формулы (43) является дискретной сверткой. В общем случае формула свертки  имеет вид:

   (44)

где f и g - функции  от целых переменных (т.е. дискретные функции). В данном случае роль таких функций играют x'_n и w_n, которые "зависят" от своих индексов как от целых переменных. Свертка обладает замечательным свойством:

F(f * g) = F(f) F(g)   (45)

где F(f) - это дискретная функция, которая получается из дискретной функции f с помощью быстрого преобразования Фурье.

Дальнейший порядок  вычисления таков. Сначала надо вычислить F(x'), то есть взять набор x'_n и вычислить БПФ для n = 0...N' − 1. Затем вычислить F(w), для чего взять набор w_n и вычислить БПФ для n = 0...N' − 1. Потом по формуле (45) простым перемножением получаем F(x' * w) = F(x') F(w). Просто берем i-й элемент, полученный после БПФ набора x'_n, и умножаем на i-й элемент, полученный после БПФ набора w_n. В результате получим БПФ свертки x' * w. Теперь у нас есть БПФ свертки, из него надо получить саму свертку. Для этого применяем обратное преобразование Фурье. Теперь нам надо найти X'_k. Для этого используем формулу (43). k-й элемент набора X_k будет соответствовать z_k-му элементу свертки. Наконец, останется перейти от X'_k к X_k.

В процессе вычислений мы применяли алгоритмы порядка  сложности N, N' и N'log₂N' (когда делали три БПФ). Поскольку N' не может превышать N более, чем в 4 раза, то весь алгоритм имеет сложность O(Nlog₂N)

Обратное преобразование требует величины W без минуса в показателе и в конце алгоритма - деления всех полученных элементов на N. Три алгоритма БПФ остаются в этом случае теми же: два прямых перобразования и одно обратное.

На следующей  страничке приведен листинг программы, выполняющей прямое и обратное преобразования по этому алгоритму. Оптимизация выполнена примерно теми же способами, что и в предыдущем случае с незначительными модификациями. Переменные x' в алгоритме обозначаются как x_, 2N как N2, N' как N_, z_k как N22, π/N (или −π/N при обратном преобразовании) как piN.