УТ БПФ - теперь быстрее двоичного и четверичного!

nabbla1 — 26.10.2013 Во введении я написал, что предлагаемое Уравновешенное Троичное БПФ требует на 10% больше арифметических действий, чем старое доброе двоичное, зато удобно в применении во многих задачах.

После многомесячного блуждания в трех соснах, торжественно заявляю:
можно так организовать вычисления атомарных операций УТ БПФ, что общее число арифметических операций составит (20/3)Nlog₃N≈4.2Nlog₂N. Для сравнения, двоичный алгоритм требует 5Nlog₂N операций, а четверичный: 4.25Nlog₂N. Алгоритмы более высоких степеней работают медленнее.

Итак, УТ БПФ - не только удобный, но и самый быстрый из всех алгоритмов с фиксированным основанием. БПФ со смешанным основанием 4/2 (split-radix FFT) работает еще быстрее, требуя 4Nlog₂N операций, но более сложен в реализации. В нем исходный массив делится на 3 неравные части, одна размером N/2 и две по N/4, и такая процедура вызывается рекурсивно. Но он не такой красивый)

Любопытно, что обычное троичное БПФ требует большего числа операций, около 5.89Nlog₂N, и лишь при его "уравновешивании" появляется новая симметрия (2 фазовых множителя становятся комплексно сопряженными друг к другу), которая и упрощает вычисления. Мало того, оказывается, похожий подход позволил в 2007 году Стивену Джонсону и Матью Фригго снизить на 6% количество операций в БПФ со смешанным основанием, побив мировой рекорд, державшийся более 20 лет! По сути, вместо деления массива X[k] на подмассивы X[2k],X[4k+1] и X[4k+3], они стали делить его на X[2k], X[4k-1] и X[4k+1], хотя из-за нумерации с нуля получился "заезд" в x[-1], который объявили тождественно равным значению X[N-1].

Мне казалось, что с математической точки зрения разница между обычным и уравновешенным БПФ тривиальна - циклический сдвиг - а вон оно как обернулось! Ладно, ближе к делу.

Раздел 3 - продолжение.

Асимптотические оценки вычислительной сложности БПФ по произвольному базису
Когда я в прошлом году прочел доклад на конференции МФТИ про УТ БПФ, мне присудили первое место в своей секции и пожелали изобрести еще для кучи пятиричное и семиричное БПФ! Понимаю, что это было сказано больше в шутку, однако когда я рассказал другу, он ответил – а что, может, они еще быстрее будут работать? Я на пальцах объяснил, что не будут, самый производительный алгоритм надо искать с основанием, близким к основанию натурального логарифма, 2.718… Не знаю, насколько стоит эти оценки приводить здесь, пока для интересу приведу. Тем более, что в самой первой статье по БПФ, с которой все началось, подобный расчет тоже был сделан с тем же результатом - троичное должно быть самым быстрым.

Атомарная операция для БПФ с фиксированным основанием R (когда на каждой итерации БПФ от R^q элементов сводится к линейной комбинации от R штук БПФ c R^q-1 элементов в каждом) состоит в умножении вектора из R элементов на матрицу RxR:

В общем случае для этого требуется R² умножений и (R-1)R сложений. При этом вычисляется сразу R значений, то есть на одно значение приходится R операций умножения. Всего итераций будет log_RN.
Полное число умножений равно

Найдя производную от R и приравняв ее к нулю, найдем, что минимум данного выражения достигается при R=e=2.71828…
Для ближайших целых чисел получаем M(2)=2.88NlnN и M(3)=2.73NlnN, то есть минимальное число умножений будет достигаться именно в троичном БПФ.
Проблема данной оценки в том, что не учтена возможность значительно убыстрить умножение на матрицы специального вида. Так, заведомо известно, что в данной матрице одна строка и один столбец будут единичными, поэтому требуется (R-1)(R-1) умножений, а не R². В этом случае аналитически найти минимум не удастся (он будет корнем нелинейного уравнения), однако прямая подстановка целых чисел покажет, что самым быстрым будет двоичное БПФ.
Но и это не может являться истиной в последней инстанции, как будет показано далее. Главное достоинство данной асимптотической оценки в том, что становится понятно - наиболее эффективны алгоритмы БПФ с малым базисом, порядка 2..4.

Подсчет числа операций в УТ БПФ
Приведем еще раз атомарную операцию УТ БПФ (2.1):

Здесь W_Nⁿ и W_N^-n– фазовые множители (англ. tweedle factors), на которые необходимо домножить F₊[n] и F_-[n] перед выполнением дальнейших действий. Они всегда возникают в быстрых преобразованиях с фиксированным основанием. Существуют алгоритмы БПФ от N=p₁p₂…p_q, где p_n – взаимно простые числа, в которых фазовые множители не возникают, например алгоритм Винограда [см. P. Duhamel, M. Vetterly]. Однако соответствующие алгоритмы гораздо сложнее, и хотя в них достигается теоретический минимум по числу умножений, они требуют существенно большего числа сложений, в результате чего они до сих пор не нашли широкого практического применения.

Множители

соответствуют поворотам точки на комплексной плоскости на углы ±π/3. В двоичном БПФ появляются повороты на ±π, то есть просто смена знака, а в четверичном – еще дополнительно на ±π/2 (умножение на ±i), что тоже выполняется тривиально.
Поворот на π/3 можно сделать тривиальной операцией, если перейти к базису Эйзенштейна, в котором комплексные числа представлены в виде x+wy, w=exp(2πi/3). Но и в стандартном декартовом базисе можно прийти к экономичному способу вычисления (2.1), использующему симметрию задачи.
В первую очередь, домножаем F₊[n] и F_-[n] на фазовые множители:
X₊=W_NⁿF₊[n],
X_-=W_N^-nF_-[n].
Для единообразия также напишем
X₀=F₀[n].
Пока мы использовали 2 комплексных умножения, или 8 действительных умножений и 4 сложения. Для краткости будем далее писать 8m+4a (8 multiplications + 4 additions) чтобы обозначить количество действительных операций.
Теперь мы можем найти F[n]:
F[n]=X₀+X₊+X_-.
На это ушло 2 комплексных сложения, или 4a.
Запишем формулы для F[n±N/3]:

Величина X₊+X_- была посчитана выше, кроме нее необходимо посчитать X₊-X_- (требует 2a), затем помножить первое на ½ (требует 2m), а второе на (i√3)/2 (требует 2m), после чего сложить (требует 6a).
Итого, для выполнения атомарной операции (2.1) требуется 12 действительных умножений и 16 сложений. Полное число умножений УТБПФ составляет:

Для сравнения, двоичный БПФ требует 2Nlog₂N действительных умножений.
Полное число сложений составляет:

Двоичный БПФ требует 3Nlog₂N действительных сложений.
Число операций в рассмотренной выше реализации УТ БПФ лишь на 12% выше, чем в двоичном БПФ.

Это еще одна причина, почему троичные БПФ не рассматривались в качестве серьезной альтернативы для двоичных, а тем более, для схем со смешанным основанием – split-radix – которые требуют лишь 4Nlog₂N операций. Рассмотрим, как можно снизить вычислительные затраты УТ БПФ и получить рекордно низкое число операций.
Работа в базисе Эйзенштейна
Любое комплексное число можно записать в виде
– кубический корень из единицы (см. рисунок)



Перевод из обычной записи и назад производится следующим образом:
(переход к базису Эйзенштейна),
(переход к прямоугольному, Гауссову базису)
Сложение комплексных чисел в таком базисе выполняется путем сложения отдельных компонент:
(a+bw)+(c+dw)=(a+c)+(b+d)w
Приведем свойства числа w:

w³=1
Умножение произвольного комплексного числа в базисе Эйзенштейна на w или w-1 тривиально:
(a+bw)w=aw-b-bw=-b+(a-b)w,
(a+bw) w^-1=a(-1-w)+b=(b-a)-aw
Именно благодаря этому удается значительно снизить количество действительных умножений в УТ БПФ.
Запишем правило умножения двух комплексных чисел в базисе Эйзенштейна:
(a+bw)(c+dw)=ac+(ad+bc)w-(1+w)bd=(ac-bd)+(ad+bc-bd)w
В такой записи требуется 4 действительных умножения и 3 сложения – на одно больше, чем в прямоугольном базисе. Есть возможность одно умножение заменить сложением. Вычисляем сначала вспомогательные величины:
α=b(c-d),
β=c(a+b),
γ=ad.
Теперь искомое выражение примет вид:
(a+bw)(c+dw)=(α+β)+(α+γ)w
Как видно, достаточно 3 умножений и 4 сложений, причем одно из них может быть выполнено заранее для фазовых множителей, на подготовительном этапе работы.
Похожий метод существует и для декартового базиса:
α=a(c+d),
β=d(a+b),
γ=c(b-a),
(a+bi)(c+di)=(α-β)+(α+γ)i
Он позволяет заменить 4u+2a на 3u+5a, причем два сложения можно выполнить заранее. Когда выполнение умножения было делом существенно более длительным, чем сложение, подобные методы часто применялись. При работе в базисе Эйзенштейна, как можно видеть, такая процедура предпочтительна всегда, поскольку снижает количество арифметических операций внутри цикла.
Наконец, найдем количество арифметических операций УТ БПФ при работе в базисе Эйзенштейна. Умножение на фазовые множители: 6m+6a. Умножение на w и w^-1: 4a. Сложение всех операндов: 12a. Полное число операций для вычисления (2.1): 6m+22a.
Получаем число умножений для УТ БПФ: M=2Nlog₃N≈1.26Nlog₂N.
Число сложений: A=(22/3)Nlog₃N≈4.62Nlog₂N
Как видим, полное количество операций одинаково в обеих реализациях, но при работе в базисе Эйзенштейна требуется в 2 раза меньше умножений. Учитывая, что в современных компьютерах умножение и сложение выполняются за сравнимое время, а переход к базису Эйзенштейна и обратный переход требуют дополнительных накладных расходов, его использование не рекомендуется.
Использование симметрии преобразования в алгоритме прореживания по частоте
Поскольку прямое и обратное преобразование Фурье отличаются лишь знаком в экспоненте и, возможно, постоянным множителем, то граф операций можно выполнять и в обратном порядке. Данный метод называют прореживанием по частоте. Именно в нем наглядно видно, что комплексно-сопряженные фазовые множители позволяют упростить вычисления.

Запишем базовую операцию для алгоритма прореживания по частоте (пожалуй, стоит в разделе 2 привести подробный вывод и этого алгоритма):

Здесь звездочка означает комплексное сопряжение, а индексы опущены, дабы не загромождать выкладки.
Перепишем выражения для F₊, F_-, учитывая, что 1+w+w*=0:

Величину Z=wW можно вычислить заранее и на этапе вычисления УТ БПФ брать его значение из таблицы.
Сначала вычисляем значения в скобках:
α=X₀-X_-=a+bi,
β=X₊-X_-=c+di.
На это требуется 4 действительных сложения (4a). Далее, обозначим
W=W_r+iW_i,
Z=Z_r+iZ_i
Тогда,
F₊=(aW_r+cZ_r)-(bW_i+dZ_i)+[(bW_r+dZ_r)+(aW_i+cZ_i)]i
F_-=(aW_r+cZ_r)+(bW_i+dZ_i)+[(bW_r+dZ_r)-(aW_i+cZ_i)]i
На вычисление каждой из 4 скобок требуется 2m+a операций, затем еще 4a для их попарного сложения, всего 8m+8a.
Итого, в одной атомарной операции выполняется 8 действительных умножений и 12 действительных сложений. Полное число умножений равно

(для двоичного, напомним, требуется 2Nlog₂N умножений, а для четверичного – 1.5Nlog₂N)
Число сложений:

Общее число операций: 4.2Nlog₂N.

Библиография по данному разделу
Использовать базис Эйзенштейна для троичного БПФ было предложено в [P. Dubois and A.N. Venetsanopoulos, "A new algorithm for the radix-3 FFT", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-26, June 1978, pp. 222-225.] я не читал этой статьи - денег за нее просят, выводил сам по аннотации. Сравнив с тем, что пишут в других местах, ссылаясь на нее, понял - то же самое)

Быстрый алгоритм прореживания по частоте предложен по крайней мере в [Todd Mateer - FAST FOURIER TRANSFORM ALGORITHMS WITH APPLICATIONS - 2008]. Это диссертация на PhD, использующая необычный формализм для описания быстрых преобразований. Из двух алгоритмов (прореживания по времени и по частоте) один, к сожалению, ошибочен. Зато второй - просто зашибись! Приведен комментарий - чтобы фазовые множители были комплексно сопряженными, надо изменить способ инверсии. Как это сделать - читателям предлагается придумать самостоятельно.

Нигде не встречал упоминания об уравновешенных преобразованиях и связи УТ БПФ с уравновешенной троичной системой счисления и алгоритмов ур. троичной инверсии. Похоже, это сугубо моё!


UPD. Похоже, ложная тревога, забыл еще 4 сложения посчитать. С ними получается 5.04Nlog₂N - на 0.9% дольше, чем двоичный. Тоже неплохо, но никакого рекорда.

Оставить комментарий

Популярные посты:

• Архив