Быстрое преобразование Фурье (стр. 7 из 7)
Функция complex_mul выполняет умножение комплексных чисел z1 и z2 и записывает результат в ячейку z.
Мы собираемся вызывать алгоритм "fft" M раз для одного и того же числа элементов L. Для оптимизации мы должны вынести за пределы алгоритма "fft" те действия, которые не зависят от выбора конкретных элементов в массиве.
К таким действиям можно отнести выделение памяти для хранения поворачивающих множителей и освобождение ее. Деление результата на N также можно выполнить позже для всего массива сразу. А главное, мы можем всего лишь однажды вычислить все поворачивающие множители. Для решения этой задачи мы напишем отдельную функцию createWstore, которая и вычисляет множители, и выделает массив ячеек памяти для них. Позднее этот массив будет передаваться как параметр в новый вариант алгоритма "fft".
Функция createWstore представляет собой ту часть основного цикла "fft", которая отвечала за рассчет поворачивающих множителей. Однако, она оптимизирована. Число итераций меньше, чем в исходном алгоритме, поскольку не на всех итерациях выполнялось вычисление нового множителя, иногда происходило только копирование. Поэтому шаг приращений индекса в массиве вдвое больше: Skew2 вместо Skew. Переменная k для определения четности/нечетности нам больше не нужна, так что внутренний цикл останавливается при достижении границы массива WstoreEnd.
Порядок вычисления множителей 
будет понятнее на примере. Если L = 32, то порядок изменения k такой (точкой с запятой разделены итерации внешнего цикла, а запятой - внутреннего):
k = 0; 16; 8, 24; 4, 12, 20, 28; 2, 6, 10, 14, 18, 22, 26; 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23...31
Далее, вместо старого алгоритма "fft" мы напишем новый - "fft_step" с учетом того, что поворачивающие множители считать не надо, и не надо делить результат на N.
Сам алгоритм практически не изменился за исключением того, что все расстояния между элементами массива увеличились в 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 - указатель на Xr+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). В скобках там стоит обычное прямое БПФ, только не для всех xn, а для последовательности из P штук, взятых с шагом Q. Количество последовательностей равно Q. Элементы последовательности выбираются формулой xpQ + q, где p - номер элемента, q - номер последовательности.
Пусть мы умеем вычислять БПФ для P элементов за время Z(P). Нам понадобится вычислить Q таких БПФ, на что потратим время O(QZ(P)). Тем самым мы получим все множители в скобках из формулы (41). После этого надо будет вычислить N элементов Xk, каждый из которых потребует O(Q) действий, итого потребуется еще O(NQ) действий, а в сумме - O(QZ(P) + QN) действий.
В предыдущем алгоритме P - было степенью двойки, а для такого числа элементов мы умеем вычислять БПФ за O(Plog2P) шагов, т.е. Z(P) = Plog2P, и получится O(QZ(P) + QN) = O(QPlog2P + QN) = O(Nlog2P + N2/P) действий. Тем самым, мы несколько более простым путем получим ту же оценку сложности алгоритма.
Если P - не степень двойки, мы умеем вычислять ПФ за O(P2) шагов, т.е. Z(P) = P2, и получится O(QZ(P) + QN) = O(QP2 + QN) = O(N2/Q + N2/P) действий. Если P и Q достаточно велики, это лучше, чем O(N2).
В "Википедии" нашелся алгоритм сложности O(Nlog2N)для произвольного N. Алгоритм основан на значительном увеличении числа элементов последовательности и дополнении ее нулями. Хотя такое расширение, естественно, будет компенсироваться высокой скоростью алгоритма.
Длина последовательности берется равной N', где N' - степень двойки ближайшая сверху к 2N. То есть, 4N > N' = 2T ≥ 2N. Т.е. число элементов новой последовательности в худшем случае может превышать число элементов исходной последовательности почти в 4 раза. Например, если N = 9, то N' = 32.
Обозначим
Представим формулу (1) в виде:
Преобразуем степень:
kn = 2kn/2 = (−k2 + k2 + 2kn + n2 − n2)/2 = (−k2 + (k + n)2 − n2)/2
Это используем для дальнейшего преобразования формулы (1):
Далее введем обозначения:
x'n = xnW −n2/2
X'k = XkW k2/2
Для n > N, полагаем xn = 0.
Тогда
(42)И далее авторы алгоритма выполняют виртуозный "финт ушами". Пусть:
zk = (2N − 2 − k)
wn = W zn2/2
И превращают формулу (42) в:
(43)Проверим:
wzk − n = W (zzk − n)2/2 = W (2N − 2 − (zk − n))2/2 = W (2N − 2 − zk + 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 и wn, которые "зависят" от своих индексов как от целых переменных. Свертка обладает замечательным свойством:
F(f * g) = F(f) F(g) (45)
где F(f) - это дискретная функция, которая получается из дискретной функции f с помощью быстрого преобразования Фурье.
Дальнейший порядок вычисления таков. Сначала надо вычислить F(x'), то есть взять набор x'n и вычислить БПФ для n = 0...N' − 1. Затем вычислить F(w), для чего взять набор wn и вычислить БПФ для n = 0...N' − 1. Потом по формуле (45) простым перемножением получаем F(x' * w) = F(x') F(w). Просто берем i-й элемент, полученный после БПФ набора x'n, и умножаем на i-й элемент, полученный после БПФ набора wn. В результате получим БПФ свертки x' * w. Теперь у нас есть БПФ свертки, из него надо получить саму свертку. Для этого применяем обратное преобразование Фурье. Теперь нам надо найти X'k. Для этого используем формулу (43). k-й элемент набора Xk будет соответствовать zk-му элементу свертки. Наконец, останется перейти от X'k к Xk.
В процессе вычислений мы применяли алгоритмы порядка сложности N, N' и N'log2N' (когда делали три БПФ). Поскольку N' не может превышать N более, чем в 4 раза, то весь алгоритм имеет сложность O(Nlog2N)
Обратное преобразование требует величины W без минуса в показателе и в конце алгоритма - деления всех полученных элементов на N. Три алгоритма БПФ остаются в этом случае теми же: два прямых перобразования и одно обратное.
На следующей страничке приведен листинг программы, выполняющей прямое и обратное преобразования по этому алгоритму. Оптимизация выполнена примерно теми же способами, что и в предыдущем случае с незначительными модификациями. Переменные x' в алгоритме обозначаются как x_, 2N как N2, N' как N_, zk как N22, π/N (или −π/N при обратном преобразовании) как piN.