Синтез алгоритмов обработки сигналов с ограничениями на минимальный параллелизм и объём памяти Салищев Сергей Игоревич

Содержание к диссертации

Введение

1 Факторы энергоэффективности 17

1.1 Общая модель энергопотребления тактируемых логических схем 24

1.2 Мощность КМОП устройств 25

1.3 Методы оптимизации энергопотребления 27

1.4 Методы логического синтеза 29

1.5 Типичная архитектура акселератора

1.5.1 Оценка рассеиваемой мощности акселератора при фиксированном размере задачи 32

1.5.2 Асимптотическая скорость роста мощности при росте размера задачи

1.6 Выбор оптимального типа памяти 36

1.7 Выбор оптимальной ширины и представления числовых данных 37

1.8 Рематериализация данных 42

1.9 Влияние перспективных технологий на энергоэффективность

1.9.1 Схемы с пороговым напряжением питания 43

1.9.2 Память на основе фазового перехода и магниторезистив-ная память 45

1.10 Влияние архитектуры программного обеспечения на энергоэффективность 45

1.10.1 Операционные системы с поддержкой страничной памяти 49

1.10.2 Операционные системы без поддержки страничной памяти 50

1.10.3 Cреды управляемого исполнения 51

1.11 Формальная верификация как средство повышения энергоэффективности

Аппаратное ускорение вычисления элементарных функций при помощи специализированных вычислительных блоков 59

2.1 Постановка задачи аппроксимации с заданной точностью 61

2.2 Задача уменьшения размера таблиц 62

2.3 Оценка точности аппроксимации на одном сегменте

2.3.1 Погрешность аппроксимации интерполяционным полиномом с ограничениями типа неравенства 64

2.3.2 Погрешность квадратичной и кубической аппроксимации интерполяционным полиномом 66

2.4 Расчёт таблиц при помощи целочисленного линейного програм мирования 68

2.4.1 Целочисленная аппроксимация на одном сегменте 69

2.4.2 Случай квазисплайна 70

2.4.3 Оптимизационная задача 71

2.5 Результаты прототипирования 75

3 Поточный алгоритм БПФ на многобанковой памяти 78

3.1 Общий подход к разработке поточных акселераторов БПФ 81

3.1.1 Синхронные графы потоков данных 81

3.1.2 Расщепляющее правило БПФ 84

3.1.3 Инверсия индексов 86

3.1.4 Формула БПФ произвольной размерности во временной области 88

3.1.5 Формула БПФ произвольной размерности в частотной области 91

3.1.6 Реализация круговой свёртки 95

3.2 Организация многобанковой памяти 96

3.2.1 Постановка задачи 97

3.2.2 Акселератор БПФ по смешанному основанию с 1r1w памятью 99

3.2.3 Акселератор БПФ по смешанному основанию с 1rw памятью 107

3.3 Самосортирующиеся БПФ 110

3.3.1 Акселератор самосортирующего БПФ с 1r1w памятью 111

3.4 Результаты прототипирования 117

4 Ускорение решения уравнений Юла–Уокера 118

4.1 Подавление дальнего эха с помощью линейного фильтра с длинной импульсной характеристикой 118

4.2 Обращение тёплицевой матрицы при помощи многочленов Сегё и Шура

4.2.1 Многочлены Сегё и факторизация обратной матрицы к тёплицевой 121

4.2.2 Проблема коэффициентов Шура 124

4.2.3 Спектральные плотности и функции Шура 126

4.2.4 Связь многочленов Сегё и Шура 128

4.3 Быстрый алгоритм Шура 131

4.3.1 Транзитивность многочленов Шура 131

4.3.2 Преобразование коэффициентов функций Шура 132

4.3.3 Структура бинарного дерева при расчёте параметров Шура133

4.3.4 Расчёт многочленов Шура по параметрам Шура 134

4.3.5 Расчёт остаточных членов 136

4.3.6 Формулировка быстрого алгоритма Шура 139

4.4 Сложность расчёта многочленов Шура 143

4.4.1 Общая оценка количества операций 143

4.4.2 Общая оценка количества адресуемой памяти 148

4.5 Оценка оптимального параллелизма и времени вычислений быстрого алгоритма Шура на устройстве с аппаратным ускорением БПФ 150

4.5.1 Оценка количества комплексных чтений в вещественном алгоритме Шура 152

4.5.2 Оценка длины критического пути 155

4.5.3 Оценка времени вычислений быстрого алгоритма Шура для вещественных данных на 2 конвейерных процессорах157

4.5.4 Оценка оптимального параллелизма 159

4.6 Гибридный алгоритм фильтрации с малой задержкой 164

Заключение 166

Список иллюстраций 168

Список таблиц 169

Литература

• Оценка рассеиваемой мощности акселератора при фиксированном размере задачи
• Операционные системы с поддержкой страничной памяти
• Погрешность квадратичной и кубической аппроксимации интерполяционным полиномом
• Общий подход к разработке поточных акселераторов БПФ

Введение к работе

Актуальность темы. Энергоэффективность - одна из основных характеристик полупроводниковых беспроводных устройств, определяющая время работы устройства от батареи и его тепловой режим. Стандартные модели энергопотребления учитывают только активную мощность, в которой оценка затрат энергии пропорциональна сложности в элементарных операциях. Для схем с литографическими нормами 45 нм и менее на первый план выходят потери энергии в результате токов утечки, которые не зависят от вычислительной сложности. Они складываются из размера памяти и количества параллельных вычислительных элементов. Оптимальный параллелизм для минимизации энергопотребления рассматривался Ву и Ли [Woo, 2008] для многоядерных суперскалярных процессоров. В отличие от рассматриваемой нами задачи, процессоры работают на высокой частоте и при высоком напряжении питания, что позволяет не учитывать энергопотребление памяти. Одним из базовых блоков для алгоритмов цифровой обработки сигналов является вычисление элементарных функций In, exp, sin, cos, л/х, 1/х и т.п. Для аппаратных блоков вычислений с точностью 22 бита и выше обычно используется кусочно-полиномиальная аппроксимация. Проблемой энергоэффективной аппаратной реализации метода является оптимальный баланс между размером таблиц и степенью многочлена. Для элементарных функций таблицы являются избыточными и могут быть сокращены, используя гладкость функций. В работе [Strollo, 2011] предложен метод сокращения таблиц на 40% без существенного увеличения вычислений за счет использования двузвенного гладкого сплайна, однако вопрос точности решается эмпирически, при помощи тестирования на всех допустимых данных, что делает метод в описанном виде неприменимым для высоких точностей. Из более сложных алгоритмов наиболее часто используется Быстрое Преобразование Фурье. Существующие алгоритмы БПФ на специализированной аппаратуре можно разделить по асимптотике времени вычислений (9(1), (9(Inn), 0(п), (9(nlnn). При требовании переиспользования ресурсов и программно управляемой длины преобразования часто наиболее энергоэффективными оказываются аппаратные блоки БПФ с временем работы (9(nlnn) и памятью с произвольным доступом. К таким алгоритмам относится поточное БПФ [Johnson, 1992], которое может быть реализовано на реконфигурируемом вычислительном устройстве антимашине Хартенштейна [Hartenstein, 1991] Разработка поточных алгоритмов БПФ с ограничением на доступ к памяти может быть сведена к задаче составления расписания для синхронного графа потока данных [Parhi, Messerschmitt, 1991]. Алгоритм Джонсона применим только для

чистых оснований и для смешанного основания без параллелизма бабочек. Модификация алгоритма Джонсона, предложенная [Jo, Sunwoo, 2005], специализирована только для оснований 2/4. Во многих задачах более эффективными являются самосортирующие алгоритмы БПФ, такие, как алгоритм Джонсона-Буруса и Темплтона [Hegland, 1994], которые сформулированы только для однобанковой памяти. Также требуется модификация под наиболее эффективную архитектуру памяти, предоставляемую библиотекой компонентов, например, однопортовую память. Такие адаптации для алгоритмов БПФ без копирования пока не рассматривались в литературе. Другим часто используемым классом алгоритмов является факторизация тёплицевых и обратных к тёплицевым матриц. Задачи этого типа и большой размерности возникают в акустических задачах эхо- и шумоподавления и разделения источников сигнала. Для решения задачи факторизации обычно используются алгоритмы Левинсона и Шура, имеющие асимптотику сложности (²), где - длина вектора автокорреляций. Для задач большой размерности может использоваться быстрый алгоритм Шура, предложенный [Ammar, Gragg, 1987; Воеводин, Тыртышников, 1987] со сложностью (ln² ), который рассматривается как пример энергоэффективной реализации сложных алгоритмов с БПФ.

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

Для достижения поставленной цели необходимо было решить следующие задачи:

1. Сформулировать функционал сложности и другие требования, предъявляемые к алгоритмам для их эффективной реализации в аппаратуре;

2. Разработать эффективные целочисленные алгоритмы вычисления элементарных функций с заданной точностью и минимальной длиной таблицы;

3. Разработать расписание для реализации графа потока данных для БПФ на многобанковой памяти;

4. Исследовать оптимальный параллелелизм и размер памяти быстрого алгоритма Шура.

Основные положения, выносимые на защиту:

5. Метод качественной оценки мощности и выбора оптимального параллелизма для энергоэффективных специализированных КМОП вычислительных блоков (Лемма , Лемма ).

6. Метод вычисления элементарных функций при помощи почти гладкого четырехзвенного квазисплайна и оценка точности полиномиальной аппроксимации с коэффициентами с фиксированной точкой, ограниченной на равномерной сетке (Теорема ).

7. Теорема о размещении данных БПФ в многобанковой памяти при вычислении по произвольным смешанным основаниям (Теорема ).

8. Теорема о размещении данных и порядке вычисления самосортирующегося БПФ (Теорема ).

9. Теорема о размещении данных и порядке вычисления БПФ для одно-портовой памяти (Теорема ).

10. Анализ энергоэффективности алгоритма факторизации вещественных тёплицевых матриц на сверточном акселераторе для задачи эхоком-пенсации при помощи быстрого алгоритма Шура (Лемма , Лемма , Лемма , Лемма ).

Научная новизна: Все представленные результаты являются новыми.

11. Разработана модель энергопотребления для малопотребляющей цифровой схемы, выполняющей известный вычислительный алгоритм. Решена задача выбора оптимального параллелизма в данной модели.

12. Задача минимизации энергопотребления при расчёте значений стандартных функций сведена к минимизации длины таблиц. Разработаны новые методы аппроксимации квазисплайнами с неравномерным табулированием, удобные для аппаратной реализации.

13. Доказана теорема о размещении данных БПФ в многобанковой памяти при вычислении по произвольным смешанным основаниям, что обеспе-чивющая максимальную скорость вычислений при заданном параллелизме. Получены явные формулы БПФ в виде произведений Кронекера по стадиям произвольных порядков.

14. Доказана теорема о самосортирующейся модификации БПФ в многобанковой памяти по смешанным основаниям, а также аналогичная теорема для вычислительного устройства с однопортовой памятью.

5. Для быстрого алгоритма Шура найден минимальный объём памяти, вычислена длина критического пути и проведена оценка оптимального параллелизма.

Практическая значимость диссертационной работы состоит в сокращении площади и энергопотребления рассмотренных элементов и повышении их универсальности, что вносит существенный вклад в улучшение энергоэффективности автономных беспроводных устройств, реализуемых на базе специализированных полупроводниковых логических схем.

Достоверность изложенных в работе результатов обеспечивается практической реализацией предложенных схем и алгоритмов в виде полупроводниковых схем. Практическая реализация включала разработку моделей на языке SystemC, автоматическую верификацию с помощью системы “Aegis for SystemC”, логический синтез полупроводниковых схем уровня вентилей для акселераторов из моделей SystemC, логическое моделирование работы схем и синтез виртуальной топологии для малопотребляющего процесса производства полупроводниковых кристаллов с геометрическими нормами 22 нм.

Апробация работы. Основные результаты работы докладывались на: Международной конференции Общества Инженеров Акустиков (AES) (Россия, Санкт-Петербург, 2003), международной конференции по компьютерному анализу и моделированию (CDAM) (Беларусь, Минск, 2004), конференции молодых ученых “Гироскопия и Навигация” (Россия, Санкт-Петербург, 2004), семинаре кафедры теоретической кибернетики СПбГУ (Россия, Санкт-Петербург, 2015, 2016), семинарах лаборатории Intel Labs (2013 - 2015).

Личный вклад. Результаты, выносимые на защиту, получены автором самостоятельно.

Публикации. Основные результаты по теме диссертации изложены в 12 печатных публикациях [–], в том числе 4 [–] — в журналах, рекомендованных ВАК, 5 [5,,–11] — в тезисах докладов на международных конференциях на английском языке, из них 3 [–11] индексируются Scopus, [] заявка на патент США.

Работы [,5,–11] написаны в соавторстве. В работе [] автору принадлежит постановка задачи, формулировка всех теорем и их доказательство, кроме доказательства теоремы 4. В работе [, ] автору принадлежит раздел, посвященный практическому опыту реализации. В [] автору принадлежит постановка задачи, анализ существующих систем для выделения общих требований, раздел по использованию Java в системном программировании. В [, 11] автору принадлежит постановка задачи и разработка алгоритма обнаружения ошибок синхронизации с помощью анализа достижимости. В

работе [5] автор выполнял математическое моделирование.

Объем и структура работы. Полный объем диссертации — 209 страниц. Основной текст диссертации — 167 страниц. Диссертация состоит из введения, четырех глав и заключения с 9 рисунками, 14 таблицами и 5 приложениями. Список литературы содержит 85 наименований.

Оценка рассеиваемой мощности акселератора при фиксированном размере задачи

Рассмотрим факторы, влияющие на энергоэффективность вычислительных блоков при построении малопотребляющих полупроводниковых схем, специализированных для группы алгоритмов. Будем учитывать и другие два дополнительных критерия - скорость разработки и удобство композиции вычислительных блоков для реализации сложных алгоритмов.

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

Поскольку рассматриваемые схемы работают непрерывно, то для харак-теризации их энергопотребления используется мощность. Мощность полупроводниковой схемы складывается из активной и статической рассеиваемой мощности. р = Ра + Pg Активная мощность обусловлена переключением состояния полупроводниковой схемы и связана с интенсивностью вычислений. Статическая мощность не зависит от вычислений, а только от количества аппаратуры, подключенной к питанию.

Сложные вычислительные блоки разрабатываются в виде тактируемых полупроводниковых схем. Схема состоит из регистров, хранящих состояние, и комбинационных логических схем, описывающих переход между состояниями. Переход между состояниями происходит по стробу тактирующего сигнала. Такая схема может быть представлена как конечный автомат. F = (R, о", Е); ЕсйС {0,1}п; а : R ь- R Здесь R - множество допустимых состояний, Е - множество начальных состояний, а - функция перехода. Функция перехода в полупроводниковой схеме реализуется как направленный ациклический граф элементарных библиотечных функций или вентилей. Вентили могут реализовывать такие функции как AND, OR, NOR, XOR, мультиплексор, селектор и т.д. На более низком уровне вентили состоят из транзисторов. Состояние системы хранится в регистрах.

Такой уровень детализации описания логических схем называется уровнем регистровых передач (Register Transfer Level). 1.2 Мощность КМОП устройств Для оценки мощности будем опираться на базовую модель, описанную во второй главе [29]. Для одного вентиля мощность определяется как р = Ра + Pg. Здесь Ра - активная мощность, зависящая от данных, Ps - статическая мощность, не зависящая от данных. Ра = Pd + PSC) Ps = Pleak Здесь Pd - динамическая мощность, Р\еак - потери мощности от токов утечки, Psc - потери мощности в результате короткого замыкания при переключении транзисторов. Р aCV f (1.1) Здесь С - паразитная емкость, V - напряжение питания, / - тактовая частота, а - коэффициент частоты переключения, показывающий среднюю вероятность открытия вентилей на каждом такте. Для сигнала строба а = 1, для других в среднем 0.1. Максимальная тактовая частота работы схемы при данном напряжении питания связана с задержкой на критическом пути, то есть цепочке вентилей, имеющей самую большую задержку переключения. / 1/ СГЙ, dcrit = / di Задержка на одном вентиле также связана с напряжением по следующей формуле: СіУ С d(V) осттт — 77? т- (1.2) (у — VT) V Здесь CL - емкость нагрузки, которая состоит из перезаряжаемых входом емкостей затворов транзисторов и паразитных емкостей проводов. Ст = C0XWL Здесь Ст - емкость затвора транзистора, VT - пороговое напряжение, Сох -удельная емкость оксида затвора, W, L - ширина и длина затвора транзистора. Из-за дефектов изготовления цифровые схемы работают неустойчиво при напряжениях, близких к VT; существует минимальное напряжение Vo VT, обеспечивающее устойчивую работу схемы. При этом напряжении схема может работать на максимальной тактовой частоте /0 = 1/dcrtt(V0). Схема может работать и при меньшей тактовой частоте, если скорость ее работы достаточна для решения задачи.

Операционные системы с поддержкой страничной памяти

Предположим, что фиксирован следующий способ расчёта аппроксимаций заданной функции / на заданном промежутке определения А. Промежуток А разбивается на отрезки, и внутри каждого отрезка функция приближается полиномом. В памяти хранятся параметры, необходимые расчёта значений полиномов по всем отрезкам. Требуется минимизировать размер этой памяти.

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

В данном разделе представлены обозначения для одной из задач целочисленного линейного программирования, которые далее потребуются для основных формулировок. Предположим, что требуется аппроксимировать функцию / многочленом р(х) степени N на заданной сетке отсчётов: J є = maxІЄ[о..р] \f(xi) —p(xi)\,xo = 0,Xi Xi-i,xp = 1, = minpe,p = 2ik=oPkXk Задача оптимальной аппроксимации на сетке отсчётов можно следующим образом сформулировать в каноническом виде: є = min їх. Здесь V = {ж } =оА;=о - матрица Вандермонда, / = (/(жо),... ,f(xp)) -вектор значений в узлах, х = (ро, -PN, Є) - вектор переменных, їх - целевая функция, где / = (0,..., 0,1).

Для нахождения округленных коэффициентов полинома используется метод смешанного целочисленного линейного программирования. Вводятся дополнительные ограничения на переменные: х Є (Qim)N+1 х R+, где Qim — множество двоичных дробей с 1т цифрами после запятой. Использование целочисленного программирования обеспечивает существенно лучшую точность по сравнению с округлением коэффициентов к ближайшей представимой точке. Это позволяет сэкономить 1-2 бита на каждом коэффициенте.

В отличие от алгоритма Ремеза, метод на основе линейного программирования позволяет не только приближать оптимальные по точности интерполяционные полиномы, но и решать задачу в условиях ограничений на достаточную точность аппроксимации и на значения полиномов в узлах. Это дости гается ценой дополнительных вычислительных затрат на решение задачи на большом наборе узлов.

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

Идея уменьшения размера таблиц полиномиальной интерполяции связана с тем, что аппроксимирующие полиномы на соседних сегментах могут иметь избыточную информацию, которую можно учесть. Таким образом, каждая секция таблицы данных содержит параметры некоторого квазисплайна.

В конце данного раздела будет показано, что при помощи решения вспомогательной задачи линейного программирования и аналитических оценок по теореме 1 длины таблиц можно существенно сократить таким способом для ряда элементарных функций.

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

Рассмотрим аппроксимирующий квазисплайн степени -1 с 20 звеньями. Будем рассматривать двух- и четырехзвенные квазисплайны. Степени двойки используются, поскольку в этом случае старшие биты аргумента используются как адрес в таблице коэффициентов, а младшие — как аргумент интерполяции, без необходимости вычислений по модулю. Двузвенный квазисплайн позволяет переиспользовать все коэффициенты, не добавляя дополнительных вычислений. Квазисплайны с числом сегментов более 2 добавляют дополнительные сумматоры в вычислительный блок, таким образом, энергоэффективность переиспользования данных уменьшается за счет роста количества вычислений. Требование гладкости является чрезмерно жестким, и его выполнение может приводить к увеличению количества сегментов кусочно-полиномиальной интерполяции и удвоению размера таблиц, что не компенсируется их последующим сокращением. Для сохранения количества сегментов в узлах квазисплайна допускаются невязки, которые на практике равны 0 на большинстве сегментов. Поскольку ограничения не являются жесткими, такая оптимизация не приводит к уменьшению шага сетки по сравнению с задачей без ограничений.

Поскольку часть коэффициентов обнуляется, предложенный метод может считаться методом неравномерного табулирования, поскольку на сегментах, где функция имеет меньший модуль производной, для ее табулирования используется меньше данных. Для одного квазисплайна необходимо хранить коэффициентов полинома и (20 - 1) дополнительных значений невязок во внутренних узлах квазисплайна. Появление невязок приводит к необходимости дополнительных таблиц и дополнительных сумматоров для внесения поправок в коэффициенты.

Такой алгоритм аппроксимации является обобщением гладкого двузвен-ного сплайна, рассматривавшегося Стролло и др. [3]. Использование более двух звеньев квазисплайна и ослабление требований гладкости позволяет существенно сократить размер таблиц. Например, переход к четырехзвенному квазисплайну по сравнению с двузвенным дает сокращение таблиц на более чем 30%. Сформулируем оптимизационную задачу нахождения невязок для почти гладкого квазисплайна.

Погрешность квадратичной и кубической аппроксимации интерполяционным полиномом

Пусть п - натуральное число и ujn = е . Дискретное преобразование Фурье (ДПФ) порядка п является линейным оператором с матрицей J n = [bJn ]o k,l n Для произвольных целых N 0и0 j N введём обозначение е дг для j–го стандартного орта в RN. Далее g) обозначает произведение Кронекера: если С = А 8 В, размерности матриц А и В равны т х п и к х , соответственно, то Cik+j,Pe+q = Ai,jBPyQ при всех 0 і т, 0 j п, 0 р к, 0 q . Далее будем использовать следующие свойства произведения Кронекера: (А ) В) = А g) В , (A S В)(С S D) = (AC) S (BD), (A S B) S C = A g (B g)C). Второе равенство справедливо только при условии корректности произведений АС и BD. Из последнего равенства, в частности, следует, что (А (g) Ik) 8 1ц = А /н и Ik {h 8 А) = Iki g А. Пусть т,к — натуральные числа и п = кт. Набор векторов (е (g) dj,m) при 0 і к, 0 j т составляет стандартный базис в Rn. Введём обозначение для матрицы перестановки 1 и для диагональной матрицы W порядка п, определяемых формулами Щ{Єі,к ej,m) = j,m еА,к, 0 І &, 0 j Ш, Vm(ei,k 8 ej,m) = п{Єі,к ej,m), 0 І &, 0 j ГП. Перестановка Щ осуществляет следующее преобразование индексов: Ц. : гт + j — jk + і, 0 і к, 0 j т. Из определения L\ следует, что для любых векторов ж Є Rk и у є Rm: L (x (g) у) = у S х. Матрица может быть также явно представлена в виде

Эти матрицы обладают следующими свойствами. Лемма 5 Пусть = . Тогда ІJ J І J_J І ІJ І = I J_J I m к к VV —— j-j VV j-j Доказательство. Утверждения непосредственно следует из определений. Произведение Кронекера некоммутативно. Перестановка 1киТтв произведении Кронекера осуществляется при помощи матриц перестановок L. Лемма 6 Пусть п = km. Тогда Щ{1к 3 т)1 ат = J m Ik Доказательство. Пусть 0 i kиO j m. Тогда Щ{1к 7m) m(ej,m Є і,к) = Щ{ і,к {J m j,m)) = {J m j,m) ei,k = (J m Ik){e-j,m Єі,к), что доказывает утверждение леммы. Алгоритмы быстрого преобразования Фурье (БПФ) основаны на факторизации матрицы ДПФ при помощи расщепляющего правила, сформулированного в следующем утверждении. Лемма 7 Пусть п = кт, где к, т — натуральные числа. Тогда J-n = Щ(1к Fra m k 1т) Доказательство. Пусть 0 і к, 0 j т. Тогда Поскольку (Lk)T = L, то при 0 p к, 0 q m аналогично получаем, что (Т Т (Ik J m) т(ЄЧ,т ер,к) ) = ( (ер,к тЄЧ,т) ) к q,m Отсюда Из определения следует, что = 1. Поэтому ipm-\-pj-\-kjq (p-\-Qk)(im-\-j) — iqn П UJ = иО %т+3 ЩП = (eg) Є k)J-n(eiyk 8 e,j,m), что доказывает утверждение леммы. п Следствие 2 Пусть п = кт, где к, т — натуральные числа. Тогда J-n = Щ(1} g) J7TO)H/ L (/TO g) Тк)І к Доказательство. Формула получается из леммы 7 при помощи подстановки утверждения леммы 6. 3.1.3 Инверсия индексов Общая формула алгоритма БПФ получается итерированием расщепляющего правила из следствия 2. Эта формула включает отображение, порождённое инверсией мультииндексов, которое определяется ниже.

Пусть а = (пк,пк-і,...,щ) - некоторый упорядоченный набор натуральных чисел, называемый далее мультииндексом, и N = Ylk=onk. Каждое целое число от 0 до N — 1 единственным образом представляется в системе счисления, порождённой мультииндексом а, в виде мультииндекса Р = (рк,Ро) из условия п = Po + o(pi + - +пк-2(рк-і+пк-іРк) ), 0 Pj nj, 0 j К. В этом случае связь между мультииндексом р и числом п обозначается следующим образом: р = па, п = ра. Мультииндекс а будем называть порождающим систему счисления, а мультииндекс р — согласованным с этой системой. Количество кодируемых чисел обозначим N =\а\.

Из определения сразу следует, что для любых порождающих мультиин-дексов а, (3 и для любых согласованных с ними мультииндексов р, q, соответственно, ера,\а\ eql3,\fS\ = e(p,q)(a P),\a\\/3\ Инверсией мультииндекса р = (рк, Ро) назовём тот же набор компонент, но в обратном порядке: р = (ро,Рк). Мультииндекс а = (по,... ,пк) также порождает нумерацию чисел от 0 до N — 1. Перестановкой инверсии мультииндекса а назовём перестановку Ра чисел от 0 до N — 1, определяемую правилом

Это же правило можно записать в виде Рап = ((па ))а. Таким образом, перестановка Ра отображает числа от 0 до \а\ — 1, записанные в системе счисления, порождённой мультииндексом а , в числа с инвертированным цифровым представлением в системе счисления, порождённой мультииндексом а.

Перестановка Ра порождает линейное отображение в RN, которое будем обозначать Sa. Оно полностью характеризуется значениями Saen = ераЩм при 0 п N.

Из определения матрицы перестановки Ьп следует, что если а = (т,п) имеет только две компоненты, то Sa = LI Лемма 8 Пусть а - мультииндекс и N = \а\. Пусть М 0 и /3 = (М, а) расширенный мультииндекс. Тогда S(M,a) = (їм Sa)LN , S(aM) = LM (їм Sa). Доказательство. Пусть 0 га Ми0 п ІУ. Для доказательства первого равенства введём обозначение для мультииндекса (3 = (М, а). Выполним вспомогательные преобразования: rnN+Рап = (m, (na Y) = [( a 5 ) ] = [{{пМ-\-т)я ) у = Ря(т-\-пМ). Подставим это равенство в следующем преобразовании: (1м Sa)L%M{en N g eTOjM) = {їм 5 а)(ет,м 8 еПіЛг) = eTOjM 8 {Saen,N) &т,М У &Pan,N &mN-\-Pan,M N &Pp(m-\-nM) (H m+nM,MN (H\ n,N У m,M)i что равносильно первому утверждению леммы 8. Аналогично доказывается второе утверждение. Обозначим 7 = («, М). Тогда т + МРап = ((па У тр1 = [(m,na ) ]7 = [((m7V + n)7 ) ]7 = P7(m7V + n). Подставим это равенство в следующем преобразовании: М V М У Ja)\&m,M У &n,N) М v т,М У a&n,N) М \&т,М У &Pan,N) &Pan,N У &т,М &NPan-\-m,MN &Pp(mN-\-n) - [H&mN-\-n,MN Э[!1\&т,М У &n,N)i что равносильно второму утверждению леммы 8.

Используем расщепляющее правило для приведения алгоритма БПФ к виду, пригодному для реализации на архитектурном шаблоне рис. 3.2.

Пусть размерность N ДПФ можно разложить на множители: N = [ [к=0 пк. Тогда алгоритм расчёта ДПФ может содержать только бабочки длиной п0,... пК. На этом свойстве основан алгоритм БПФ, который обычно применяется для размерностей N, равных степени числа 2.

Общий подход к разработке поточных акселераторов БПФ

Доказательство. Пусть і — номер стадии по основанию R. За один такт обрабатывается ровно одна бабочка. Старший разряд номеров всех отсчётов одной бабочки одинаковый. Чётность этого разряда определяет группу из R банков памяти, из которых производится считывание на данном такте. Старший разряд вместе с группой из R банков памяти меняются на каждом такте в соответствии с определением Ti(k) и разложения номера к в системе счисления, порождённой мультииндексом а . Поскольку длина конвейера нечётна, то на каждом такте записать результатов ”на месте” осуществляется не в ту группу банков, из которой производится считывание. Поэтому не происходит конфликта памяти между операциями считывания и записи.

Для того чтобы номера банков двух отсчётов совпадали, rri2(k) = m t) требуется, чтобы совпадали номера т(к) = т() из теоремы 9 и, кроме того, чтобы совпадала чётность младшей компоненты в системе счисления, порождённой мультииндексом а .

На стадии по основанию R порядок обхода бабочек тот же, что и в теореме 9. Следовательно, конфликта памяти при считывании не возникает по теореме 9. На каждой стадии один элемент памяти участвует только в одной бабочке.

Младшим разрядом номера бабочки рг является число ко, определяющее чётность номера банка памяти. Чётность номера такта совпадает с чётностью номера банка памяти в операции считывания. Если запаздывание по записи нечётное, то операция записи выполняется в банки памяти с противоположной чётностью. Отсюда не возникает конфликта операций считывания и записи.

Инверсия индексов Sa в алгоритме БПФ требует дополнительного прохода по памяти. Поскольку точки с инвертированными индексами могут лежать в одном банке памяти, то из-за конфликтов обращения к памяти время этого прохода равно времени выполнения двух стадий БПФ. При этом во многих алгоритмах, использующих БПФ, задействовать на этом проходе вычислительные блоки невозможно. Такой простой не является энергоэффективным, и его желательно избегать.

Для решения этой проблемы существует класс само сортирующих алгоритмов БПФ, не требующих явного шага инверсии индексов. Для этого была разработана модификация самоcортирующего алгоритма Джонсона-Буруса [8]. В алгоритме БПФ Джонсона-Буруса стадия по малому основанию выполняется посредине. С точки зрения потоковой архитектуры этот алгоритм может быть представлен следующим образом:

Половина стадий не вычисляется на месте. Однако если рассматривать группы бабочек по R элементов, то они выполняются на месте. Для отсутствия конфликтов предшествования между бабочками в одной группе достаточно, чтобы внутренний конвейер вычислительного блока удовлетворял условию В данном разделе излагается обобщение алгоритма Джонсона-Буруса на БПФ с выполнением малой стадии в начале.

Для реализации этого алгоритма на архитектурном шаблоне рис. 3.2 достаточно изменить функции адресации и увеличить длину внутреннего конвейера.

Лемма 12 Пусть а и Ь — квадратные матрицы, и матрицы А = In (g а, В = Ь (g 1т имеют одинаковый размер N. алгоритма Тогда если тп делится на N, то АВ = В А. Доказательство. Размер матрицы а обозначим к, а размер матрицы В — . Из условия следует, что N = кп = т. По условию, q = mn/N — целое число. Отсюда так что малая стадия по основанию г = R/q выполняется в начале. В алгоритмах БПФ как во временной, так и в частотной области, индексы компонент входного массива записываются в системе счисления, порождённой мультиин-дексом а = (г, Л,..., R). По теореме 6 алгоритм БПФ в частотной области описывается формулой TN = SaEn_i7aEn_2,a E0jCn где Sa — перестановка компонент вектора в соответствии с инверсией мульти-индекса а. Цель самосортирующегося состоит в замене начальной перестановки на постепенные частные перестановки на каждой стадии, осуществляемые при минимальной длине конвейера.

Введём оператор перестановки начальной и конечной цифр в индексе произвольного вектора. Пусть вектор имеет длину N и число N делится на к2, где к 1. Тогда матрица перестановки Y определяется уравнениями