CC BY
43
Оптимизация преобразования Фурье под архитектуру Эльбрус
П.А. Ишин ЗАО МЦСТ, старший научный сотрудник, ishin@mcst.ru Введение
Дискретные преобразования Фурье (ДПФ) имеют очень широкое распространение. Они применяется в алгоритмах цифровой обработки сигналов, при сжатии звука и изображения. ДПФ помогает решать частные дифференциальные уравнения и выполнять такие операции, как свертки. ДПФ также активно используются в статистике, при анализе временных рядов. В этой статье рассказано о способе реализации ДПФ для микропроцессора Эльбрус, позволяющего эффективно использовать возможности этого микропроцессора.
Особенности архитектуры микропроцессора Эльбрус
Микропроцессор Эльбрус характеризуется высокой архитектурной скоростью. Это свойство поддерживается технологией оптимизирующей языковой компиляции, для которой, в свою очередь, предусмотрена специальная аппаратная поддержка [1].
Рассмотрим особенности, поддерживающие высокую архитектурную скорость:
1) Большой набор исполняющих устройств с возможностью явно управлять их запуском с использованием широкой команды;
2) Большой регистровый файл, организованный в виде процедурных окон произвольной длины с автоматической откачкой/подкачкой окна при переполнении/исчерпании регистрового файла;
3) Наличие специальных операций и режимов исполнения, позволяющих преодолевать ограничивающие параллелизм зависимости:
a) режим спекулятивных (упреждающих) вычислений, которые нужны, чтобы уменьшить зависимости, связанные с условной передачей управления;
b) исполнение операций под управлением предиката позволяет избавиться во многих случаях от операций передачи управления, что ведет к увеличению параллелизма и, как следствие, к росту архитектурной скорости;
c) для преодоления статически неопределяемых зависимостей по данным предусмотрены операции спекулятивного (упреждающего) считывания данных в обход предполагаемой зависимости с последующим контролем корректности перестановки операций.
4) Специальные операции подготовки выполнения переходов, позволяющие распараллелить передачу управления по нескольким направлениям без потери тактов на выбор соответствующего пути;
5) Операции, поддерживающие лучшую адаптацию программы к иерархической памяти, которые представляют собой аппаратные средства синхронной и асинхронной предподкачки данных при регулярной обработке, а также средства управления размещением данных в кэшах различного уровня;
6) Аппаратная поддержка программной конвейеризации циклов, включающая в себя средства переименования регистров с использованием специальной вращающейся базы, а также счетчик цикла и специальные предикаты, управляющие циклом.
Микропроцессор Эльбрус позволяет выполнять за один такт до 10 арифметических операций, из них до 4 операций с плавающей арифметикой, до 7 команд считывания из памяти (или 6 считываний и одну запись в память, или 4 считывания и две записи). Общее количество вращающихся регистров в конвейеризованном цикле 128.
Дискретное преобразование Фурье и его свойства
ДПФ конечного вектора ^^^ , называется выра-
жение:
Х(к)= £х(п)е Ы ,к = 0,1,...., N-1, (1)
N-1 --кП
или в другой форме
N-1
Х(к)=X ФЖ„к, (2)
п=0
2пп
где = е
= „ N
N
Коэффициент называется поворачивающим множителем.
В случае, когда Х(к) является комплексным, при прямом вычислении ^точечного ДПФ нужно выполнить для каждого значения к (N-1) умножений и (N-1) сложений комплексных чисел или 4(И-1) умножений и 2(И-1) сложений действительных чисел. Для всех N значений к=0,1,...,И-1 требуется (Ы-1)2 умножений и N(N-1) сложений комплексных чисел. Для больших значений N (порядка нескольких сотен или
0
п=
тысяч) прямое вычисление ДПФ по выражению (1) требует выполнения весьма большого числа арифметических операций умножения и сложения, что затрудняет реализацию вычислений в реальном масштабе времени.
Быстрое преобразование Фурье (БПФ)
Основная идея БПФ состоит в том, чтобы разбить исходный N отсчетный вектор X(k) на два более коротких вектора, ДПФ которых могут быть скомбинированы таким образом, чтобы получить ДПФ исходного ^отсчетного вектора.
Разобьем исходный вектор x(n) на два Ж2-отсчетных вектора xl(n) и X2(n), составленных соответственно из четных и нечетных отсчетов исходного вектора x(n)
N
хх(п) = х( 2п ),п = 0,1,...., — -1,
х2(п) = х( 2п +1 ),п = 0,1,...,у -1.
Тогда выражение (2) можно записать следующим образом:
N2-1 N2-1
Х(к)= £х(п)ШтПк ^ £х2(п)Шы
или
Х(к) = Х1(к)+ЖкХ2(К), (3)
где Xl(k) и X2(k) - Ж2-отчетные ДПФ векторов xl(n) и X2(n) соответственно.
Получаем, что вычисление X(k) сводится к вычислению X1(k) и X2(k). Далее для вычисления X1(k) и X2(k) можем применить этот же метод, разбив эти вектора пополам. В результате получаем общее число
ЛГ
операций — N .
Базовая операция алгоритма (так называемая «бабочка») состоит в том, что два входных числа A и B объединяются для получения двух выходных чисел X и Y по правилу
п=
п=
При комплексном векторе, на вычисление одной «бабочки» требуется 10 арифметических операций - 4 умножения и 6 сложений. Более подробно особенности алгоритмов ДПФ и БПФ описаны в [2] (раздел "Быстрое Преобразование Фурье").
Изменение алгоритма БПФ под архитектуру Эльбрус
Рассмотрим стандартный алгоритм БПФ с 16 точками (Краткое описание алгоритма описано выше, пример реализации приведен в [3], разделы "Реализация" и "Улучшенная реализация"). Вот так выглядит схема преобразований (см. Рис. 1)
Рис. 1. Стандартная структура БПФ
Как видим, в данном алгоритме при переходе к обработке следующего слоя удваивается шаг бабочки. В архитектуре Эльбрус реализована асинхронная подкачка данных(АРВ), которая позволяет заранее подкачивать значения как из Ь2 кеша, так и из общей памяти. Наиболее эффективно подкачка работает, когда подкачиваемые элементы располагаются в памяти подряд, поэтому стандартный алгоритм БПФ будет не столь эффективен, из-за того, что там данные располагаются достаточно далеко в памяти.
Чтобы избежать этих проблем, данные записываются в память в другой последовательности. Результат вычисления бабочки располагаются так, чтобы при обработке следующего слоя считывания из памяти происходили подряд. На Рис.2 схематично изображен новый порядок сохранения данных в память.
Рис. 2. Измененная структура БПФ
При такой структуре, происходит считывание и запись в память подряд (более подробно о данной структуре и её преимуществах в [4]).
Правда, для реализации этого алгоритма требуется дополнительная память, так как, в отличие от стандартной реализации БПФ, результат вычислений сохраняется в другом порядке.
Обработка первого слоя
При реализации алгоритма необходима перестановка отсчетов входного вектора, чтобы выходная последовательность имела естественный порядок расположения отсчетов, т.е. к=0,1,...,М-1. В приведенном примере для 16-точечного БПФ для этого требовался следующий порядок размещения отсчетов входного вектора: х(0), х(8), х(4), х(12), х(2), х(10), х(6), х(14), х(1), х(9), х(5), х(13), х(3), х(11), х(7), х(15). Закономерность перестановки заключается в том, что отсчеты входного вектора должны
быть размещены в памяти в двоично-инверсном порядке. Это означает, что требуемый номер ячейки памяти для размещения очередного отсчета входного вектора определяется обратной перестановкой двоичных разрядов в двоичном представлении номера отсчета.
Для этого задается специальная таблица BitReversTab, с помощью которой загружается нужный элемент входного вектора. При помощи APB считывание данных таблицы не замедляет работу программы.
Вычисление поворачивающих множителей.
Есть, по крайней мере, три способа вычисления поворачивающих множителей.
Первый способ - это вычисление поворачивающих множителей напрямую по формуле
2k nk
Wk=e N = cos[(2nN)k] -/sin [(2nN)k].
Этот способ вычисления связан с большими временными затратами, так как используются "тяжелые" операции синуса и косинуса.
Второй это заполнение специальной таблицы множителей. Тогда при вычислениях достаточно просто считывать нужные значения. Но чем выше порядок преобразования, тем больше требуется таблица, а значит, работа с памятью становится все проблематичнее.
Третий способ - это вычисление с помощью рекуррентной формулы
wN+! = WNWN . (4)
Этот способ менее затрачен по времени, чем первый и не требует обращение к памяти. Правда, у него более низкая точность результата.
При реализации алгоритма была использована комбинация второго и третьего способа. Для вычисления Фурье малых порядков (меньше 14) используется заранее вычисленная таблица. Для вычисления больших порядков, используется рекуррентная формула.
При этом если множитель W^ есть в таблице, то происходит считывание его оттуда, а последующие множители , W¡k+2 ,W¡k+3 ... вычисляются по рекуррентной формуле (4) до тех пор, пока не дойдем до следующего множителя из таблицы. Такой способ позволяет увеличить точность вычисления множителей и при этом не хранить большую таблицу поворачивающих множителей.
Комбинированные операции
Как было описано выше, для вычисления одной «бабочки» требуется 10 арифметических операций - 4 умножения и 6 сложений. Архитектура Эльбрус позволяет объединять последовательные операции MUL и ADD(SUB) в одну операцию MUL_ADD(SUB).Т.е. последовательность команд:
a3 = MUL (a1, a2)
a5 = ADD (a3, a4)
заменяется одной командой:
a5 = MUL_ADD (a1,a2,a4).
Это позволяет сократить количество арифметических операций до 8(4 команды ADD и 2 MUL и 2 MUL_ADD). Широкая команда позволяет нам выполнять 4 плавающие арифметические операции за такт. Это означает, что для вычисления одной «бабочки» нужно два такта в формате 64 битных чисел с плавающей точкой (тип double).
Объединение слоев
Для ускорения вычислений решено обрабатывать сразу несколько слоев одновременно. Как было описано выше, для вычисления одной «бабочки» требуется 8 арифметических операций. Для вычисления сразу 2 слоев нам потребуется вычислить сразу 4 бабочки, или 32 арифметические операции, для 3 слоев - 12 бабочек и 96 операции, для 4 слоев 32 бабочки и 256 операций. Всего в конвейеризированном цикле используется не более 128 вращающихся регистров (всего количество регистров 192). Исходя из этого, получается, что максимально можно объединить 3 слоя, дабы цикл смог конвейеризоваться. Такое объединение позволяет сэкономить на операциях считывания и записи в память. Так же объединение сразу 3 слоев позволяет сэкономить 4 арифметические операции за счет раскрытия скобок, при умножении на поворачивающие множители.
Вычисление БПФ для векторов больших размеров
Описанный выше алгоритм хорошо работает, для векторов, которые помещаются в L2 кеш, так как темп работы с L2 кешем позволяет подкачивать оттуда данные без задержек. Когда же вектор не помещается в L2 кеш, то APB подкачивает данные из памяти, темп работы с которой гораздо ниже. Чтобы их минимизировать, применяется блочный алгоритм БПФ.
Из свойства преобразования Фурье следует, что можно сделать преобразование над четной и нечетной частями вектора. Затем по формуле (3) посчитать БПФ для исходного вектора. Аналогично можно разде-
лить вектор на 4, 8 и.т.д. частей. На основе этого свойства, реализован следующий алгоритм:
Исходный вектор делится на несколько частей, которые позволяют вести вычисления, не выходя за границы L2 кеша. Далее, для каждого из таких векторов, производится БПФ нужного порядка. Когда происходит обработка последнего слоя, для каждого из этих векторов происходит подкачка данных для выполнения БПФ для следующего вектора. Это позволяет уменьшить затраты при работе с памятью. Затем вектора постепенно объединяются, и происходит обработка оставшихся слоев.
Заключение
Описанный метод реализации БПФ позволяет эффективно использовать возможности микропроцессора Эльбрус.
1) Комбинированные операции позволяют уменьшить количество арифметических операций.
2) С помощью асинхронной подкачки данных можно работать с L2 кешем без задержек.
3) Наличие конвейеризированных циклов и большое количество вращающихся регистров позволяют обрабатывать сразу по 3 слоя, благодаря чему экономится арифметические операции и обращения в память.
4) При обработке больших векторов, применяется блочный метод и подкачка данных, что позволяет уменьшить потери при обращении в память.
Как было написано выше, на одну бабочку приходится 8 арифметических операций. Получается, что теоретический максимум для микропроцессора Эльбрус это 2 такта на бабочку для комплексных векторов типа double.
Данный алгоритм позволяет достичь на размере 1024(10 порядок) показателя 2.15 такта на бабочку (93% от пиковой). На размере 16284(14 порядок) 4.13 тактов на бабочку (48% от пиковой).
Такая большая разница в производительности между 10 и 14 порядками, объясняется выходом за пределы L2 кеша и снижением скорости подкачки данных из памяти по сравнению с L2 кешем.
Для сравнения, стандартный алгоритм БПФ показывает на размере 1024(10 порядок) 3.78 такта на бабочку (53% от пиковой). На размере 16284(14 порядок) 10.37 тактов на бабочку (19% от пиковой). Лучший известный алгоритм FFTW [5] показывает следующие результаты 1024(10 порядок) 2.93 такта на бабочку (68% от пиковой). На размере 16284(14 порядок) 9.31 тактов на бабочку (21% от пиковой).
Благодаря дополнительной возможности векторизации, для переменных типа float получаются результаты примерно в два раза лучше,
чем для типа double(для размера 1024 1.15 тактов на бабочку), но подробное описание этой реализации выходят за рамки данной статьи.
Литература
1. В.Ю. Волконский. Оптимизирующие компиляторы для архитектуры я явным параллелизмом команд и аппаратной поддержкой двоичной совместимости. 2004
2. Цифровая обработка сигнала. Лекция. К.М. Денисов. . http://www.ets.ifmo.ru/denisov/dsp/index.htm
3. Быстрое преобразование Фурье за O (N log N). Применение к умножению двух полиномов или длинных чисел. http://e-maxx. ru/algo/fft_multiply
4. Boris Lerner. Writing Efficient Floating-Point FFTs for ADSP-TS201 Ti-gerSHARC® Processors. 2004
5. "Fastest Fourier Transform in the West.". http://www.fftw.org/
