Алгоритмы - построение и анализ (1021735), страница 193
Текст из файла (страница 193)
Из уравнения (30.4) можно записать ДПФ как произведение матриц у = К,а, где К, — матрица Вандермонда, содержащая соответствующие степени ь~„: 1 1 1 1 Мв 2 з Уо р-1 2(ч-1) о~„ З(ч — 1) шн а1 У1 ~ и 1'~а а2 аз Уз н 1 2(н — 1) З(и-1) (н-1)(н-1) 1 1од ~ть ьб . ь~и Ев-1 Уп — 1 Для выполнения обратной операции, которая записывается как а = РГТ„1 (у), необходимо умножить у на матрицу У'„1, обратную К,. Часть Ч11.
Избранные темы 942 Теорема 30.7. Для 3,6 = 0,1,...,п — 1, (т',к)-й элемент матрицы Ъ'„1 равен ~п "/п. Даказавыиьство. Покажем, что У„1У„= 1„, единичной матрице размера и х и. Рассмотрим (3',2')-й элемент матрицы Ъ;, ~у'„: и — 1 Ъ„], = '~у ~ы™3/и) (ы"1'/1 = ~у ю"0' 31/и . в=о ь=о Согласно лемме о суммировании (лемма 30.6), данная сумма равна 1, если у' = у, и 0 в противном случае. Заметим, что поскольку — (и — 1) < 3' — 3 < и — 1, у'-3 не кратно п, так что лемма применима. Если дана обратная матрица Ъ;, 1, обратное дискретное преобразование Фурье ПГТ,,1 (у) вычисляется по формуле 1 ь-1 а- = — ~ увы„ — ьу 1=О (30.11) для 1 = 0,1,...,и — 1.
Сравнивая уравнения (30.8) и (30.11), мы видим, что если модифицировать алгоритм БПФ вЂ” поменять ролями а и у, заменить ы„ на м,,1 и разделить каждый элемент результата на п — получится обратное ДПФ (см. упражнение 30.2-4). Таким образом, 13РТ,,1 также можно вычислить за время 6(п18п). Таким образом, с помощью прямого и обратного БПФ можно преобразовывать полинам степени не выше п из коэффициентной формы в форму значений в точках и обратно за время О (и 18 и).
Применительно к умножению полиномов, мы показали следующее. Теорема 30.8 (Теорема о свертке). Для любых двух векторов а и 6 длины п, где и является степенью 2, справедливо соотношение а З 6 = 11г'Т~,~ (13ГТз„(а) . ь1ГТз (6)), Упражнения 30.2-!. Докажите следствие 30.4. 30.2-2. Вычислите ДПФ вектора (О, 1, 2, 3).
где векторы а и 6 дополняются нулями до длины 2п, а " '* обозначает покомпо- нентное произведение двух 2п-элементных векторов. Глава 30. Полиномы и быстрое преобразование Фурье 943 30.2-3. Выполните упражнение 30.1-1, используя схему с временем выполнения 6(п1бт1). 30.2-4. Напишите псевдокод для вычисления 13ГТ„1 за время 6 (и 1к и). 30.2-5. Опишите обобщение процедуры БПФ для случая, когда п является степенью 3. Приведите рекуррентное соотношение для времени выполнения и решите его.
*30.2-б. Предположим, что вместо выполнения и-элементного БПФ над полем комплексных чисел (где п четно), мы используем кольцо Е„, целых чисел по модулю т = 2'и1з+1, где $ — произвольное положительное целое число. Используйте ь1 = 21 вместо ь~и в качестве главного значения корня и-й степени из единицы по модулю тп. Докажите, что прямое и обратное ДПФ в этой системе является вполне определенным. 30.2-7. Покажите, как для заданного списка значений го, з1,..., ги 1 (возможно, с повторениями) найти коэффициенты полинома Р (х) степени не выше и + 1, который имеет нули только в точках хо, г1,..., г„1 (возможно, с повторениями).
Процедура должна выполняться за время О (п 18~ п). (Указанисс полипом Р (х) имеет нулевое значение в точке х тогда и только тогда, когда Р (х) кратен (х — зу).) * 30.2-8. Чирп-преобразованием (сЫ1р 1гапз1опп) некоторого вектора а = (ао, а1, ..., аи 1) является вектор у = (уо, у1,..., у„1), где уя = 2 " а х"1, а з — неюторое комплексное число. Таким образом, ДПФ является частным случаем чирп-преобразования, полученного для х = ь1и. Докажите, что для любого комплексного числа з чирп-преобразование можно вычислить за время О (п1кп). (Указанигс воспользуйтесь уравнением п-1 у,ь'7з'У (,, '7 ) (,-1 -31'7 ) чтобы рассматривать чирп-преобразование как свертку.) 30.3 Эффективные реализации БПФ Поскольку практические приложения ДПФ, такие как обработка сигналов, требуют максимальной скорости, в данном разделе мы рассмотрим две наиболее эффективные реализации БПФ.
Сначала будет описана итеративная версия алгоритма БПФ, время выполнения которой составляет 6 (и 1к п), однаю при этом в 6 скрыта меньшая константа, чем в случае рекурсивной реализации, предложенной в разделе 30.2. После этого мы вернемся к интуитивным соображениям, приведшим нас к итеративной реализации, для разработки эффективной параллельной схемы БПФ. Часть Ч!1.
Избранные темы 944 Итеративная реализация БПФ Прежде всего заметим, что в цикле Рог в строках !0-13 процедуры Кис~)кь)чн РРТ значение ь~„у„вычисляется дважды. В терминологии компиляторов ь 1)) такое значение называется общим подвыражением (соттоп зпЬехргезз)оп).
Можно изменить цикл так, чтобы вычислять это значение только один раз, сохраняя его во временной переменной б: Рог/с -0 гоп/2 — 1 пот< — щу) !)) Уь ' — У), +г [о) )о) Уь+(и/з) ' Уь ы < — щмч Выполняемая в данном цикле операция (умножение м„уь, сохранение произвеь О! дения в переменной 1, прибавление и вычитание 1 из уь ), известна как лреобра1о) зование бабочки (Ьпцег))у орегайоп), схематически представленное на рис. 30.3.
Теперь покажем, как сделать структуру алгоритма БПФ итеративной, а не рекурсивной. На рис. 30.4 мы организовали входные векторы рекурсивных вызовов в обращении к процедуре Кис))кз)чн РРТ в древовидную структуру, в которой первый вызов производится для п = 8. Данное дерево содержит по одному узлу для каждого вызова процедуры, помеченному соответствующим входным вектором.
Каждое обращение к процедуре Кис))нз)чн РРТ производит два рекурсивных вызова, пока не получит 1-элементный вектор. Мы сделали первый вызов левым дочерним узлом, а второй вызов — правым. Посмотрев на дерево, можно заметить, что если удастся организовать элементы исходного вектора а в том порядке, в котором они появляются в листьях дерева, то выполнение процедуры Кнсикз)чн РРТ можно будет представить следующим образом. Берутся пары элементов, с помощью одного преобразования бабочки вычисляется ДПФ каждой пары, и пара заменяется своим ДПФ. В результате получается вектор, содержащий и/2 2-элементных ДПФ.
Затем эти ДПФ у,' ь .-,..... в-,....-.—.у„) .;~ ~*,'ч ' о2", Г,' а) б) Рис. 30.3. Преобразование бабочки. а) Слева поступают два вводимых значения, оз„" умножается на уь, соответствующие сумма и разность выводятся справа. б) Упрей) щенная схема преобразования бабочки, используемая в параллельной схеме БПФ Глава 30. Полиномы и быстрое преобразование Фурье 945 '. ~~ч.а~ а! аьа4 аа аа а7),2 Г (аа,а,,а„,аа) ~ ; (а„а~,а;,а.) ) ()аа,а~)~~ ~~(а.,а,) ) , )а,,а,) );)аьаг) ) .("'), .Ё (., ((аа)) ))а4)1 ~(аг)) !(а„)) Йа;)) Наа)~ Йаа)) (!а;) Рис. 30.4. Дерево входных векторов рекурсивных вызовов процедуры Квс))аз)чв РРТ.
Начальный вызов производится для и = 8 обьединяются в пары, с помощью двух преобразований бабочки вычисляются ДПФ для каждых четырех элементов вектора, объединенных в четверки, после чего два 2-элементных ДПФ заменяются одним 4-элементным ДПФ. Теперь вектор содержит п/4 4-элементных ДПФ. Процесс продолжается до тех пор, пока не получится два п/2-элементных ДПФ, которые с помощью п/2 преобразований бабочки можно объединить в юнечное п-элементное ДПФ. Чтобы превратить описанную последовательность действий в код, используем массив А [О..п — 1], который первоначально содержит элементы исходного вектора а в том порядке, в ютором они перечислены в листовых узлах дерева на рис. 30.4.
(Позже мы покажем, как определить этот порядок, который называется норазрядно обратной перестановкой (Ь)ьгечегза1 реппщайоп).) Поскольку объединение в пары необходимо выполнять на каждом уровне дерева, введем в качестве счетчика уровней переменную я, которая изменяется от 1 (на нижнем уровне, где составляются пары для вычисления 2-элементных ДПФ) до 1я п (в вершине, где объединяются два п/2-элементных ДПФ для получения окончательного результата).
Таким образом, алгоритм имеет следующую структуру: ! !ог в — 1 го )яп 2 г!о$ог)с -О!оп — 1Ьу2' 3 до Обьединяем два 2' '-элементных ДПФ из А [!а .. )а + 2' т — 1] и А[В + 2' т .. )а + 2' — 1] в одно 2'-элементное ДПФ в А[В .. )а + 2' — 1] Тело цикла (строка 3) можно более точно описать с помощью псевдокода. Копируем цикл !Ьг из процедуры Кассина)че РРТ, отождествив у[о) с А[)а..)а + 2' '— — 1] и у1)1 с А [)с+ 2' т..)а+ 2' — 1]. Поворачивающий множитель, используемый в каждом преобразовании бабочки, зависит от значения я; это степень ы где тп = 2'. (Мы ввели переменную т исключительно для того, чтобы формула была удобочитаемой.) Введем еще одну временную переменную и, которая позволит нам выполнять преобразование бабочки на месте, без привлечения до- 946 Часть ЧП.
Избранные темы полннтельной памяти. Заменив строку 3 в общей схеме телом цикла, получим следующий псевдокод, образующий базис параллельной реализации, которая будет представлена далее. В данном коде сначала вызывается вспомогательная процедура В1т КечекБе Согч(а, А), копирующая вектор а в массив А в том порядке, в котором нам нужны эти значения. 1теклт1че РРТ(а) 1 В1т Кечекзе СОРУ(а,А) 2 и — 1еиуГЛ[а] С> и является степенью 2. 3 аког в — 1 Го 1к и 4 доги -2' 5 2~Н/тл т б Гог 1с — 0 1о и — 1 Ьу т 7 бо ы — 1 8 тог2' -0 гога/2 — 1 9 бо 1 — ыА[Й+ 7'+ гп/2] 1О и — А[/с+ 2] 11 А[/с + 2] ~ — и + 1 12 А [к + 7 + гп/2] — и — ~ 13 со < — сом~д Каким образом процедура В1т Кечекзе Согч размещает элементы входного вектора а в массиве А в требуемом порядке? Порядок следования листов на рис.
30.4 является обратной перестановкой, или реверсом битов (Ь11-гечегза1 реппп1абоп), т.е. если геч (1с) — и-битовое целое число, образованное путем размещения битов исходного двоичного представления к в обратном порядке (" задом наперед"), то элемент аь следует поместить в массиве на место А [геч (/с)]. На рис. 30.4, например, листовые узлы следуют в таком порядке: О, 4, 2, 6, 1, 5, 3, 7; в двоичном представлении эта последовательность записывается как 000, 100, 010, 110, 001, 101, О! 1, 111; после записи битов каждого значения в обратном порядке получается обычная числовая последовательность: 000, 001, 010, 011, 100, 101, 110, 111.
Чтобы убедиться, что нам нужна именно обратная перестановка битов, заметим, что на верхнем уровне дерева индексы, младший бнт которых равен О, помещаются в левое поддерево, а индексы, младший бит которых равен 1, помещаются в правое поддерево. Исключая на каждом уровне младший бит, мы продолжаем процесс вниз по дереву, пока не получим в листовых узлах порядок, заданный обратной перестановкой битов.
Поскольку функция геч (/с) легко вычисляется, процедуру В ге Кечекзе Согч можно записать следующим образом: Глава 30. Полнномы н быстрое преобразование Фурье 947 В1т Кечейзе СОРУ(а, А) 1 п — 1епдса[а[ 2 Гог)с — Осоп — 1 3 с)о А[геч()с)) — аь Итеративная реализация БПФ выполняется за время О (и!8 п). Обращение к процедуре Вп' Кечекзе Сору(а,А) выполняется за время 0(п18п), поскольку проводится п итераций, а целое число от 0 до п — 1, содержащее 18 и битов, можно преобразовать к обратному порядку за время 0 (18 п)з. (На практике мы обычно заранее знаем начальное значение п, поэтому можно составить таблицу, отображающую )с в геч (lс), в результате процедура Взт Кечекзе Сору будет выполняться за время О (и) с малой скрытой константой. В качестве альтернативного подхода можно использовать схему амортизированного обратного бинарного битового счетчика, описанную в задаче 17-1.) Чтобы завершить доказательство того факта, что процедура 1тевлпче РРТ выполняется за время О (и !я п), покажем, что число Ь (п) выполнений тела самого внутреннего цикла (строки 8-13) составляет 6 (и!Еи).
Цикл Гог (строки 6-13) повторяется п/гп = и/2' раз для каждого значения в, а внутренний цикл (строки 8-13) повторяется ги/2 = 2' 1 раз. Отсюда !кп !яи Г.( )=')' — ", 2'-'= )'-"=О( 18 ) в=1 в=1 Параллельная схема БПФ Чтобы получить эффективный параллельный алгоритм БПФ, можно использовать многие свойства, которые позволили нам реализовать эффективный итеративный алгоритм БПФ. Мы будем представлять параллельный алгоритм БПФ в виде схемы, которая похожа на схемы сравнения, описанные в главе 27. Вместо компараторов в схеме БПФ используются преобразования бабочки, показанные на рис. 30.36. В данном случае также применимо понятие глубины, введенное нами в главе 27. Схема процедуры РдкАьЕЕь РРТ, вычисляющей БПФ для п введенных значений при и = 8, показана на рис.