РЕФЕРАТ
по дисциплине «Анализ временныхрядов»
на тему:
«Алгоритмфильтрации, пример на основе БПФ»
СОДЕРЖАНИЕ
Введение
1. Основы алгоритмов БПФ
2. АлгоритмБПФ с прореживанием по времени
3. Программаи пример реализации алгоритма БПФ с прореживанием по времени
4. АлгоритмБПФ с прореживанием по частоте
5. Применениеметода БПФ для вычисления обратного ДПФ (ОДПФ)
6. ПрименениеБПФ для вычисления реакции ЦФ
7. Другиебыстрые алгоритмы вычисления дискретного преобразования Фурье
7.1Обобщенный алгоритм Кули-тьюки с произвольным основанием с множителями поворота
7.2 Алгоритмпростых множителей
7.3 Алгоритмвинограда
8. Анализточности реализации алгоритмов БПФ
Заключение
Списокиспользованных источников
ВВЕДЕНИЕ
В основепреобразования Фурье(ПФ) лежит чрезвычайно простая, но исключительно плодотворная идея – почтилюбую периодическую функцию можно представить суммой отдельных гармоническихсоставляющих (синусоид и косинусоид с различными амплитудами A,периодами Т и, следовательно, частотами ω).
Неоспоримымдостоинством ПФ является его гибкость – преобразование может использоваться какдля непрерывных функций времени, так и для дискретных.
ПФ частоприменяется при решении задач, возникающих в теории автоматическогорегулирования и управления, в теории фильтрации и т.д. Разберем один изпримеров. Имеется некий линейный фильтр – изготовленный то ли в виде набораспаянных между собой резисторов, конденсаторов и катушек индуктивности, то ли ввиде модульной конструкции интегральных микросхем. Известен также входнойсигнал (на рис. 1 в качестве входного сигнала изображена дельта-функция, тоесть импульс исчезающе короткой длительности и бесконечно большой амплитуды).Необходимо определить, какой сигнал появится на выходе нашего фильтра.
/>
Рисунок 1 — Исследование линейногофильтра
Ход решенияэтой задачи зависит от того, какую позицию мы предпочтем. Выберем временнойпуть решения (верхняя половина рис. 2) – придется входной сигнал записать какфункцию времени SBX(t) и использовать импульснуюхарактеристику фильтра h(t), то есть математическую запись егоработы во времени. Отправимся по частотному пути (нижняя половина рис. 2) –нужно будет оперировать уже не с самим входным сигналом, а с его спектром gbx(ω).Δа и алгоритм работы нашего фильтра потребуется представить в частотнойобласти – в виде частотной характеристики K(ω). Δля этоговоспользуемся помощью «магического стекла» ПФ.
/>
Рисунок 2 — Быстрое преобразованиеФурье
Итак, двапути – какой из них избрать? По-видимому, тот, который проще. Во всяком случае,в большинстве практических задач предпочтение отдается частотному направлению.
Есливыполнять ДПФ входной последовательности, впрямую – строго по исходной формуле,то потребуется много времени (особенно если количество входных отсчетоввелико). Конструктивнее использовать принцип «разделяй и властвуй», лежащий воснове алгоритма БПФ. Согласно ему входная последовательность делится на группы(например, четные и нечетные отсчеты), и для каждой из них выполняется ДПФ, азатем полученные результаты объединяются. В итоге получается ДПФ входнойпоследовательности – и существенная экономия времени. Поэтому описанныйалгоритм так и назвали – быстрое преобразование Фурье.
В данномреферате рассмотрим более подробно быстрое преобразование Фурье.
1. ОСНОВЫАЛГОРИТМОВ БПФ
Дискретное преобразование Фурье (ДПФ) X(k) конечнойпоследователъности х(пТ), п=О, 1,..., N-1 определяется согласно(1.1), (1.2):
/> (1.1)
/> (1.2)/> />
причем является /> периодическойпоследовательностью с периодом N, так как /> т=О, />1, />2… Непосредственное вычисление ДПФ (5.1) прикомплексных значениях х(пТ) требует для каждого значения k (N — 1)умноженийи (N — 1) сложений комплексных чисел или 4 (N -1) умножений и (2N — 2) сложений действительных чисел, а для всех N значений k=0, 1,…,N -1 требуется примерно N2умножений и N2сложений комплексных чисел. Таким образом, для больших значений N (порядканескольких сотен или тысяч) прямое вычисление ДПФ (1.1) требует выполнениявесьма большого числа арифметических операций умножения и сложения, чтозатрудняет реализацию вычисления в реальном масштабе времени процессов испектров.
Быстрым преобразованием Фурье называют набор алгоритмов, реализациякоторых приводит к существенному уменьшению вычислительной сложности ДПФ (1.1).Исходная идея этих алгоритмов состоит в том, что N-точечная последовательностьразбивается, на две более короткие, например на две (N/2)точечныепоследовательности, вычисляются ДПФ для этих более коротких последовательностейи из этих ДПФ конструируется ДПФ исходной последовательности. Для двух(N/2)-точечных последовательностей требуется примерно /> умножений комплексных чисел, т. е. число умножений(а также сложений) уменьшается примерно в 2 раза. Аналогично вместо вычисления ДПФ(Н/2)-точечной последовательности можно вычислить ДПФ для двух (Н/4)-точечныхпоследовательностей и таким образом вновь уменьшить требуемое число умножений исложений. Если N=2v, v>O и целое, то процесс уменьшенияразмера ДПФ может быть продолжен до тех пор, пока не останутся только2-точечные ДПФ. При этом общее число этапов вычисления ДПФ будет равно v=log2N, а число требуемых арифметических операций для вычисления N-точечнойДПФ будет порядка Nv, т.е. уменьшается примерно в N/log2Nраз. Так, при N=1000 для прямого вычисления ДПФ согласно (1.1) требуетсяпримерно N2= 106 операций комплексных умножений исложений, а при использовании алгоритмов БПФ таких операций требуется всего порядка104, т. е. объем вычислений сокращается примерно на два порядка.
Рассмотрим два алгоритма БПФ: с прореживанием по времени (вкоторых требуется перестановка отсчетов входной последовательности х(пТ)) ис прореживанием по частоте (в которых требуется перестановка отсчетов выходнойпоследовательности Х(k)).
2. АЛГОРИТМ БПФ С ПРОРЕЖИВАНИЕМ ПО ВРЕМЕНИ
Пусть требуется вычислить ДПФ (1.1) при N=2v,гдеv>O. где v>O целое (если N />2V, то можно последовательность х(пТ) дополнитьв конце нулевыми элементами так, чтобы длина результирующей последовательностибыла степенью 2).
Разобьем исходную N-точечную последовательность х(пТ)=хv(п), где v =log2 N. п=О...., N -1, на две(N/2)-точечные последовательности Хv-1,0(п) и Хv-1,1(п), состоящие соответственно из четных и нечетных членов х (пТ), т.е.
/> (2.1)
При этом N-точечное ДПФ (1.1) можно записать в виде
/>(2.2)
/>Учитывая, что получаем
/> (2.3)
Где /> и /> — (N/2)-точечные ДПФсоответственно последовательностей /> и />:/> />
/> ; .
Так как Xv (k) должно быть определено для N точек(k=0, 1,..., N-1), а Хv-1,0(k)и Хv-1,1(k) определяются толькo для N/2 точек (k=0,1,..., N/2-1), доопределим (2.3) для значений k=N/2, N/2+1,...,N-1; учитывая,что Хv-1,0(k) и Хv-1,1(k)периодическиефункции с периодом Н/2, можно записать
Xv (k+N/2)= Хv-1,0 (k+N/2) + /> Хv-1,1 (k+N/2)= Хv-1,0 (k)- /> Хv-1,1 (k),
/> (2.4)
Так как />= =-1.
Формулы (2.3) и (3.4) дают алгоритм вычисления N-точечной ДПФчерез (N/2)-точечных ДПФ. Этот алгоритм можно представить направленным графом,имеющим вид «бабочки»
/>
Рисунок 3 – граф имеющий вид «бабочки»
(рис.3, а), в котором выходные числа с и d получаютсяиз входных чисел а и b по правилам
/> (2.5)
В качестве примера граф на рис.3, б представляет операции(2.3) и (2.4). Аналогично можно теперь выразить (N/2)-точечные ДПФХv-1,0(k) и Хv-1,1(k) через (N/4)-точечные ДПФ:
/>(2.6а)
И
/> (2.7б)
где Хv-2,0(k) и Хv-2,1(k) — соответственно (N/4)-точечные ДПФ четных Хv-2,0(n) и нечетныхХv-2,1(п)
членов последовательности Хv-1,0 (п), а Хv-2,2(k)и Хv-2,3(k)-соответственно (N/4)-точечные ДПФ четных Хv-2,2(п) и нечетных Хv-2,3 (п) членов последовательности Хv-1,1(п).
Процесс уменьшения размера ДПФ от М до M/2, где Мравна степени 2, продолжается до тех пор, пока на v-м шаге (v = log2N, где N-исходный размер ДПФ) не окажутся только 2-точечные ДПФ Ф(k),k=0,1, для двухточечных последовательностей />(п), п = 0,1, определяемые из соотношений/> />
(2.8)
Последние вычисляются без операции умножения.
Пример 1. Построим алгоритм БПФ с прореживанием по времени для N=8=23,v=3, т. е. для последовательности х(пТ), п=0,1,2, З,4.5,6,7. Разобьемсогласно (2.1) исходную последовательность х (пТ) = х3 (п)на две последовательности: x2,0 (п) и х2,1(n),-состоящиесоответственно из четных и нечетных членов х3 (п):
/>
\ (2.9)
/>
Рисунок 4 — Алгоритм 8-точечного БПФ
Теперь вновь разобьем последовательности (2.9) напоследовательности из нечетных и четных членов последовательностей (2.9):
/>(2.10)
Последовательности (2.10) являются уже двухточечными.
Теперь, используя алгоритм, представленный графом «бабочка» (см.рис.3, а), строим алгоритм 8-точечного БПФ (рис.4). Вначале построимисходный массив. Как видно из (2.10), он из элементов последовательности х(n)=х(nТ),n=0,1… 7, причем на входах первого графа «бабочка» первой ступени помещаютсячисла х(0) и х(4). На входах второго графа «бабочка» -числа х(2)и х(6), на входах третьей «бабочки» — х(l) и х(5) и на входахчетвертой «бабочки» — x(3) и x (7).
Таким образом, если предположить, что последовательность Х(п) записываетсяв массив ячеек памяти, то удобно осуществить размещение х(п) в следующемпорядке (рис.2): х(0), х(4), Х(2), х(6), х(1), х(5), х(3), х(7). Легкозаметить, что элементы этой последовательности получаются из исходной х(п) всоотnетствии с двоичной инверсией номеров, т. е. число х(п) с номеромв двоичном представлении п=(пv-1,..., n0)запоминается в ячейке памяти с номером />=(n0,..., пv-1). Так, число х(4) сномером в двоичном представлении 4(10) = 100(2)запоминается в ячейке с номером 001(2)=1(10), а число х(3),где 3(10) =011(2), запоминается в ячейке с номером110(з) =6(10) и т. д. Итак, можно считать, что начальнаяступень преобразования Х0(k), k=0,I… 7, получатся просто врезультате прореживания (в указанном смысле) исходной временнойпоследовательности х(пТ), п=0, I...7, т. е Х0 (k) =х (/>T), где k = /> — двоично-инверсное представление номера п.
На выходах N/2 = 4 «бабочек» т=l-й ступениобразовываются значения Х2(k), являющиеся входными числами«бабочек» т=2-й ступени. На выходах последней значения выходнойпоследовательности ХЗ (k)=X(k), k=0.… 7. Выходнаяпоследовательность X(k), k=0,1...7, получается в естественномпорядке следования.
Как показано в рассмотренном примере, все входные числа «бабочек» Х0(k) на начальной ступени являются элементами заданной последовательности х(п),п=0.… N-1, причем получаются из х(п) в соответствии сдвоичной инверсией номеров т. е. число х(пТ)=х(п) с двоичнымпредставлением номера п является входным числом Х0(k) «бабочки»с номером k, равным инверсному двоичному представлению номера п.
Заметим, что в рассмотренном алгоритме БПФ можно выполнитьвычисления по способу с замещением. Если разместить входную последовательность Х0(k) «бабочек» в массиве из 2v ячеек памяти,то после вычисления выходов «бабочек» входные элементы становятся ненужными и вуказанные ячейки памяти могут быть записаны вычисленные выходные числа. Наследующей ступени вновь вычисленные значения выходов «бабочек» записываются вячейки массива вместо использованных входных чисел, и в конце вычислений вовходном массиве окажутся записанными значения X(k) в естественномпорядке, т. е. значения ДПФ при k=0, I, 2… N-1.
3. ПРОГРАММА И ПРИМЕР РЕАЛИЗАЦИИ АЛГОРИТМА БПФ С ПРОРЕЖИВАНИЕМ ПОВРЕМЕНИ
Ниже при водится программа вычисления БПФ с прореживанием повремени по способу с замещением и рассматриваются при меры реализации этойпрограммы.
Программа 1 — быстрое преобразование Фурье с основанием два ипрореживанием по времени. Программа осуществляет алгоритм БПФ с основанием дваи прореживанием по времени комплексной или вещественной последовательности х(п)длиной N отсчетов. Вещественные составляющие отсчетов исходной последовательностизаписываются в массив А1(N), а мнимые — массив А2(N). В программедля ознакомления с ее работой предусмотрено формирование входнойпоследовательности, соответствующей отсчетам полигармонического сигнала/> />
(3.1)
строки (80-240). При использовании программы для выполнения БПФпроизвольной последовательности необходимо заменить строки 80-240, организовавввод исходной последовательности.
/>
/>
Основными этапами обработки являются: ввод исходных данных (строки50-240), двоично-инверсная перестановка исходной последовательности (строки250-350), собственно алгоритм БПФ (строки 360-510), расчет амплитуд и фазанализируемого сигнала по результатам БПФ (строки 520-590) и вывод результатов(строки 600-690). Пользователю выводятся в виде таблицы значения номеракомпоненты (гармоники) БПФ, вещественная и мнимая ее составляющие [Аl (1) и А2 (1)],амплитуда и фаза соответствующей гармоники [R (1) и Fl (1)].
Пример 2. Реализация БПФ вещественного сигнала /> содержащего три составляющие при значенияхпараметров: А0=2, w0=/>0=0, А1=I, w1=0,125,/>1=0,7854, А2=3,w2=0,3125, />2=1,57.
/>
В качестве исходных данных последовательно вводятся значения:
N=16; J=3; А(0)=2; w(0)=0; w1(0)=0; A(1)= 1; w(1)=0,125;w1 (1)=0,7854; А (2)=З; w(2)=0,3125; w1(2)= 1,57; I 9=1;
Пример 3. Реализация БПФ комплексного сигнала (3.1), содержащего трисоставляющие (J=3), при значениях параметров Ak, wkи />kтаких же, как
/>
/>
в примере 2. Ввод исходных данных аналогичен примеру 2, заисключением того, что значение I 9=0.
4. АЛГОРИТМ БПФ С ПРОРЕЖИВАНИЕМ ПО ЧАСТОТЕ
Рассматриваемый ниже алгоритм вычисления ДПФ (1.1) отличается тем,что входная последовательность х(пТ), п=0,..., N -1, разбиваетсяна две последовательности посередине (т. е. одна последовательность для n=0...N/2-1,а другая — для п=N/2...N-1) и эта процедура продолжается для каждойновой последовательности до тех пор, пока не получается искомая выходная одноэлементнаяпоследовательность Х (k); при этом величины Х (k) уже оказываютсяв выходном массиве в «прореженном» порядке и их приведение к естественномупорядку связано с инверсией двоичного представления индексов k ввычисленных значениях Х (k).
Итак, запишем ДПФ (1.1) в виде:
/> (4.1)
Учитывая, что />=/>получаем
/>. (4.2)
Подставив вместо k в (4.2) значение 2k или (2k+ 1),получим выражения для четных и нечетных отсчетов ДПФ: .
/>; (4.3)
/>, (4.4)
где теперь для значений/>:
Х0(п)=х(п) +x(n+N/2); (4.5)
Х1(п)=х(п) -х(n+N/2). (4.6)
/>
Рисунок 5 — Выполнениебазовой операции «бабочка»
Следовательно, вычисление N-точечного ДПФ X(k) сводитсяк вычислению двух N/2-точечных ДПФ причетных и нечетных значениях k дляфункций х0(п) и x1(п) и выполнениюбазовойоперации «бабочка» (рис.5)отличие операции «бабочка» здесьзаключаетсяв том, что комплексноеумножение выполняется после операциисложения-вычитания.
Ту же процедуру можно теперь применить к x0(п) их1(п) и перейти от N/2-точечных ДПФ кN/4-точечным ДПФ и, такимобразом, свести вычисление Х (2k) и Х(2k + 1) через Х (4k), X(4k+2), X(4k+ 1), X(4k+3). Продолживэтот процесс, перейдем в конечном итоге к 2-точечным ДП Ф с последующим прямымвычислением всех выходных отсчетов Х (k). Полный алгоритм БПФ спрореживанием по частоте и его программная реализация аналогичны рассмотреннымвыше для метода БПФ с прореживанием по времени.
Необходимо отметить, что в обоих алгоритмах БПФ — и спрореживанием по времени, и с прореживанием по частоте требуется примерно Nlog2N операций (комплексных умножений) и оба алгоритма могут бытьреализованы по способу с замещением, используя только один массив ячеек памяти.В обоих алгоритмах должна быть предусмотрена процедура двоичной инверсии — навходе (при прореживании по времени) или на выходе (при прореживании почастоте).
5. ПРИМЕНЕНИЕ МЕТОДА БПФ ДЛЯ ВЫЧИСЛЕНИЯ ОБРАТНОГО ДПФ (ОДПФ)
По определению (1.2) ОДПФ х(пТ) N-точечнойпоследовательности X(k), k=0, 1,..., N-1, выражается соотношением
/>(5.1)
причем в общем случае и х(пТ), и Х(k)-комплексные. Пусть х(пТ)и Х*(k)-последовательности, комплексно сопряженныесоответственно с х(пТ) и X(k). Согласно (5.1) можно записать
/> (5.2)
Но выражение суммы в правой части (5.2) есть прямое ДПФ последовательностиХ*(k), k=0,..., N-1, и, следовательно, эта сумма может бытьвычислена при помощи рассмотренных алгоритмов и программ БПФ.
Таким образом обеспечивается вычисление последовательности Nx*(пТ) и для определения х(пТ) остается взять комплексно сопряженное сNx*(пТ) выражение и разделить его на N:
/> />
(5.3)
6. ПРИМЕНЕНИЕ БПФ ДЛЯ ВЫЧИСЛЕНИЯ РЕАКЦИИ ЦФ
Вычисление реакции у(пТ) ЦФ с импульсной характеристикой h(пТ),п=0, 1,..., N-1, на входное воздействие х(пТ), п=О, 1,...M-1, может быть выполнено на основе алгоритма свертки
/> (6.1)
при п=0, 1,..., N+M-2.
Применение алгоритмов БПФ позволяет выполнить эффективноевычисление выходной последовательности у(пТ) ЦФ. С этой целью следуетопределить ДПФ H(k) и X(k) в N+M-1 точках дляпоследовательностей h(пТ) и х(пТ), затем определить ДПФ Y(k)=H(k)X(k)выходной последовательности у(пТ). Вычисление у (пТ) по ОДПФ Y(k)выполняется, например, по алгоритму (5.3). Для вычисления ДПФ и ОДПФиспользуются алгоритмы БПФ. Отметим, что если длина М последовательностих(пТ) велика, то реализация упомянутого выше алгоритма вычисления у(пТ)связана со значительной временной задержкой (для накопления всех М выборокх(пТ)). С целью уменьшения этой задержки можно входнуюпоследовательность х(пТ) разбить на отрезки xi(пТ) каждыйдлиной L и обрабатывать каждый из них независимо от других. Представим
/> (6.2)
Тогда можно (6.1) записать в виде/> />
(6.3)
Где частная свертка/> />
(6.4)
Таким образом можно начинать расчет методами БПФ частных сверток иформировать у(пТ) путем соответствующего суммирования элементов частныхсверток [2].
7. ДРУГИЕ БЫСТРЫЕ АЛГОРИТМЫ ВЫЧИСЛЕНИЯ ДИСКРЕТНОГО ПРЕОБРАЗОВАНИЯФУРЬЕ
7.1 Обобщенный алгоритм Кули-Тьюки с произвольным основанием смножителями поворота
Кроме рассмотренных выше классических алгоритмов БПФ известных какалгоритмы Кули-Тьюки по основанию 2, известно множество других. Некоторые изних позволяют существенно повысить эффективность вычисления дискретногопреобразовании Фурье. Так, алгоритмы Винограда при равном числе сложении требуютпримерно в 5 раз меньше умножений, чем алгоритмы Кули-Тьюки.
В основе всех известных алгоритмов лежит принцип разбиенииисходного ДПФ на совокупность мало точечных. Различие заключается в способахвычисления мало точечных алгоритмов и последующего объединения частичныхрезультатов. При этом размер преобразования не обязательно равен степени двух,т. е. становится возможным БПФ произвольной длины, что очень важно для рядапрактических задач. Так, в технике связи при цифровом преобразованиимногоканальных сигналов размер БПФ определяется числом объединяемых каналов.
Кратко рассмотрим только некоторые, наиболее важные алгоритмы, наоснове которых впоследствии возникло множество различных эффективныхмодификаций. Это: 1) обобщенный алгоритм Кули-Тьюки с произвольным основанием смножителями поворота; 2) алгоритм простых множителей Гуда-Винограда; 3)алгоритм Винограда.
Для простоты изложения везде будем полагать N=N1N2,где N — длина преобразования. Очевидно, приводимые ниже положениялегко могут быть перенесены на более общий случай, когда />
Обобщенный алгоритм Кули-Тьюки с произвольным основанием смножителями поворота. Итак, пусть N=N1N2, гдеN1и N2 — положительные целые. Покажем, что в этом случаевычисление исходного N-точечного ДПФ можно свести к вычислению N1N2-точечных и N2N1-точечныхДПФ и N умножениям на множители поворота />. Для этого в выражения для ДПФ (1.1)/> />
(7.1)
где />,необходимо сделать подстановку:
k=k1 +k2N2, k1=0, ….,N2-1; k2=0, …., N1-1; (7.2)
n=n1 +n2N2, n1=0, ….,N2-1; n2=0, …., N2-1. (7.3)
/>
Рисунок 6 — Сигнальный граф алгоритма
Тогда ДПФ (7.1) преобразуется к виду/> />
(7.4)
Таким образом, полученный алгоритм включает в себя две основныеступени: на первой ступени переставленные в соответствии с (7.3) входныевыборки подвергаются N2-точечному преобразованию Фурье. На второйступени производится вычисление N1-точечных ДПФ. Между первой ивторой ступенями осуществляется операция поворота путем умножения наповорачивающие множители />. Полученнаяпоследовательность на выходе ДПФ должна быть переставлена в соответствии с (7.2).
Пример 4. Пусть N=6, N1=3, N2=2. Положимk=k1+k2*2; n=n1+ +n2*3; n1k2=0,1,2;n2k1=0, 1.
Тогда/> />
Соответствующий сигнальный граф алгоритма изображен на рис.6.
Рассмотренный подход может быть положен в основу синтеза алгоритмовБПФ Кули-Тьюки с произвольным постоянным основанием. Наибольшую популярностьполучили алгоритмы с основаниями 4 и 8, позволяющие повысить эффективность вычисленияДПФ по сравнению с классическими алгоритмами по основанию 2. Отметим, чтоалгоритмы с основанием 2 также могут быть получены с использованиемрассмотренного подхода. Таким образом, рассмотренный метод синтеза являетсяобщим и позволяет синтезировать различные модификации алгоритма Кули — Тьюки спроизвольными постоянным и смешанным основаниями.
7.2 Алгоритм простых множителей
В случае, когда N представимо произведением взаимно простыхмножителей, имеется возможность избавиться от поворачивающих множителей вразложении (7.4). Тем самым можно достигнуть еще большей экономии числаопераций.
Чтобы избавиться от множителей поворота, нужно произвестиперестановку входной и выходной последовательностей, отличную от (7.2) и (7.3).Такой перестановкой может быть следующая:
для входной последовательности
/> (7.5)
n1.=0,...., N1-1; п2=0,..., N2-1;
для выходной последовательности
/>(k1N2+k2)=X((s1k2N1 + s2k1N2)modN), (7.6)
k1=0,..., N1-1; k2=0,...,N2 -1,
где записьп mod N означает «остаток от деленияпна N», а s1и s2 определяются изследующих уравнений в соответствии с китайской теоремой об остатках овосстановлении целого числа по его вычетам: s1N1=1mod N2, s1== 1modN1,s2 Тогда алгоритмN=N1N2 — точечного ДПФ представляется в виде/> />
Таким образом, алгоритм простых множителей (АПМ) является способомпредставления одномерного ДПФ в виде многомерного, причем размерность зависитот числа взаимно простых сомножителей N. Алгоритм простых множителейимеет ступенчатую форму объединения мало точечных преобразований. В данномслучае на первой ступени производится N1N2-точечныхДПФ, а на второй ступени — N2 N1-точечных ДПФ.
Впервые АПМ был предложен Гудом [3]. В том случае, когдаиспользуемые мало точечные алгоритмы синтезированы оптимальным образом пометоду Винограда, получается алгоритм Гуда — Винограда [3]. Оптимальные мало точечныеалгоритмы БПФ синтезируются путем сведения мало точечного ДПФ к совокупностициклических сверток. Для последних Виноградом [2] доказана теорема осуществовании алгоритма вычисления с минимальным количеством умножений и былпредложен метод синтеза, основанный на последовательном вычисленииполиномиальных вычетов по неприводимым полиномам в поле рациональных чисел всоответствии с полиномиальным вариантом китайской теоремы об остатках [3].
7.3 Алгоритм Винограда
Дальнейшей экономии вычислений в случае разложения N навзаимно простые множители можно достичь, если ступенчатый характер объединениячастичных мало точечных преобразований заменить вложенным. В этом и заключаетсяидея алгоритма Винограда. Идею вложения мало точечных алгоритмов легче всегопонять на примере.
Пример 5. Рассмотрим случай 6-точечного ДПФ, т. е. N=6. ПустьN1=2, N2=3. Приведем сначала алгоритмы мало точечных(2- и 3-точечных) ДПФ, синтезированные оптимальным образом по методу Винограда.
Алгоритм 2-точечного ДПФ
/>
имеет вид
/> (7.8)
где si-сложения, mi-умножения; Аiи ai — выходные и входные числа.
Алгоритм 3-точечного ДПФ
/>
Имеет вид
/>, />, />
/>, />, /> (7.9)
/>, />, />
Преобразуем исходную 6-точечную последовательность в двумернуюточечную в соответствии с (7.5) и (7.6). Тогда матрицу 6-точечного ДПФ можнопредставить в виде прямого произведения 2- и 3-точечных ДПФ и преобразованиеможно записать в виде:
/>, (7.10)
/>, />, />, />,
Где
/>
Применим к матричному преобразованию (7.10) алгоритм 2-точечногоБПФ (7.8). В результате найдем векторы/>и />:
/>
Для вычисления векторов />и /> используем алгоритм3-точечного БПФ (7.9).
Таким образом, мы как бы «вложили» алгоритм 3-точечного БПФ вструктуру 2-точечного, который оперирует 3-точечными векторами. Характернойособенностью «вложенных» алгоритмов является то, что требуемое число умноженийдля N-точечного алгоритма равно произведению числа умножений, требуемыхдля каждого из частичных алгоритмов. Этот факт легко проверяется на приведенномвыше алгоритме.
В заключение отметим, что принцип «вложения» малоточечныхалгоритмов применяется также для вычисления N-точечных циклических сверток,когда N разлагается на взаимно простые множители, если имеютсядостаточно эффективные алгоритмы малоточечных сверток. Вложенные алгоритмыциклических сверток получили название по имени авторов алгоритмов Агарвала-Кули[3]. По сравнению с традиционными алгоритмами вычисления свертки с использованиемБПФ алгоритмы Агарвала-Кули позволяют сэкономить число умножений почти напорядок.
Другими классами еще более экономных алгоритмов ДПФ и сверткиявляются алгоритмы, основанные на теоретико-числовых и полиномиальныхпреобразованиях, с которыми можно познакомиться в [3].
8. АНАЛИЗ ТОЧНОСТИ РЕАЛИЗАЦИИ АЛГОРИТМОВ БПФ
При реализации алгоритма БПФ, как и, при реализации другихалгоритмов цифровой обработки сигналов, возникают вычислительные ошибки, округлением(усечением) произведений, квантованием коэффициентов и. возможно процедуроймасштабирования чисел (с целью устранения переполнения сумматоров). При анализеошибок будем принимать допущения об их характере, т. е. будем рассчитыватьвыходную ошибку БПФ как суперпозицию ошибок, вызванных каждым независимымисточником ошибок.
Методику оценки вычислительных ошибок БПФ рассмотрим на примеререализации БПФ по основанию 2 и с прореживанием по времени. Рассматриваемаяметодика может быть применена и для анализа ошибок других алгоритмов БПФ. Будемпредполагать, что: обрабатываемые числа представляются с помощью b1+1разрядов, а коэффициенты — с помощью b2+1 разрядов с учетом знака;для аппроксимации произведений используется операция округления;масштабирование промежуточных результатов производится на входе каждой операции«бабочка» путем сдвига чисел на один разряд вправо (деление на два); входныеданные пронормированы таким образом, что/>, и подчиняются равномерному закону распределения т. е. имеютматематическое ожидание равное нулю, и дисперсию /> равную 1/3. Следовательно, среднеквадратическое значение (СК3)входной последовательности равно также 1/3. В соответствии с теоремой Парсеваля
/>
выходная последовательность X(k) БПФ будет иметь СК3 N/3,гдеN-размер преобразования. С целью уяснения методики анализаошибок получим алгоритм БПФ в аналитическом виде. Для этого в выражение для N=2v-точечного ДПФ (5.1)
/> (8.1)
необходимо подставить двоичное разложение коэффициентов п иk:
/> (8.2)
в результате алгоритм БПФ можно представить, как ранее убедились,в виде v+1-ступенчатого процесса.
На ступени т=0 производится двоично-инверсная перестановкавходной последовательности
/>
На каждой из остальных v ступеней (т=1, 2,….,v) производитсяпреобразование типа «бабочка» выходной последовательности предыдущей ступени.
Так, на ступени т = 1 производится преобразованиепоследовательности
X0(n1,…,пv):
/>
На ступени m=2 -преобразование последовательности Х1 (п1,...,nv-1, k1)
/>
На m-й ступени
/> (8.3)
Так постепенно в двоичном представлении индекса последовательностиХmс увеличением т происходит замена коэффициентовni на kj
Наконец при т = v
/>
Выходная последовательность последней ступени является искомой:
X(k)=X(kv2v-1 +… +k1)=Xv(kv2v-1 +kv-12v-2+.… +k1).
Представим индекс элемента т-й ступени в виде
(n1,..., nv-m,km,,..., k1)=2v-1n1+2v-2n2+…+2mnv-m+2m-1km+…+k1=2ml+q (8.4)
Тогда число Xm(2ml+q) можно рассматривать как q-й элементl-го блока т-й ступени.
Пример б. Рассмотрим описанную процедуру синтеза алгоритма БПФ спрореживанием по времени на примере 16-точечного ДПФ. В этом случае v=4.Индексы n и k представим следующим образом: n=n423+n322+n22+n1,k=k423 +k322 +k22+k1.
Подставляя n и k в выражение для 16-точечного ДПФ,получаем
X(k4k3k2k1)= /> /> /> />(n4nЗn2n1) />
/>
/>
Теперь распишем алгоритм по ступеням:
т = 0 — инверсия входной последовательности:
Х0(n1 23+n2 22+n32+n4)=x(n4 23+n3'22+n22+n1);
т= l — преобразование «бабочка» последовательности Х0:
/>
т = 2 — преобразование «бабочка» последовательности X1:
/>
m=3 — преобразование «бабочка» последовательности Х2:
/>
m=4 — преобразование «бабочка» последовательности ХЗ:
/>
Искомая выходная последовательность
/>
/>
Рисунок 7
Направленный граф полученного в примере 16-точечного БПФ суказанием источников элементарных ошибок округления произведений имасштабирования, а также путей их распространения изображен на рис.7. Ошибкиокругления появляются в местах умножения комплексных чисел на нетривиальныеповорачивающие множители.
Ошибки масштабирования появляются на обоих входах каждой операции«бабочка» из-за сдвига входных чисел на один разряд вправо (деление на два).
Элементарные ошибки округления и масштабирования, возникающие наразличных ступенях алгоритма, приводят к тому, что отсчеты ДПФ X(k) навыходе вычисляются неточно.
Обозначим ошибку вычисления k-ro отсчета ДПФ через e(k)=X'(k)-X(k), где Х'(k)-вычисленное значение отсчета.
Анализ траекторий распространения элементарных ошибок, возникающихна различных ступенях алгоритма, позволяет сделать следующие выводы:
1. Элементарная ошибка с дисперсией />, возникающая на т-й ступени алгоритма. вызывает ошибку сдисперсией /> в 2v-m выходныхточках БПФ.
2. Дисперсия ошибки каждого выходного отсчета алгоритма,обусловленной округлением произведений, равна сумме ошибок, 2v-m изкоторых вызывается ошибками т-й ступени.
Перейдем к анализу арифметических ошибок. В силу того, чтоматематическое ожидание всех элементарных ошибок равно нулю, математическоеожидание результирующей ошибки E(e(k))=0 для всех k. Тогда СКЗошибок ДПФ будет определяться только дисперсией элементарных ошибок.
Дисперсия ошибок округления произведения двух комплексных чисел нат-й ступени равна 4/>. Множитель /> появляется из-замасштабирования на т предыдущих ступенях путем деления пополам.
Вычислим СКЗ ошибок e(k), k=0,...,N-1. Из анализа алгоритмаБПФ (8.3) и соответствующего направленного графа следует, что нетривиальныеумножения, связанные с вычислением отсчета Х (k), появляются, начиная соступени
/> (8.5)
где. /> Это объясняется тем, чтопервое появление
единицы в двоичном представлении />, />, соответствует умножениюна коэффициент -1 на ступени т, тогда на следующей т+ l-й ступениналичие коэффициента /> = 1 приведет к умножениюна j, а уже на т+2 и далее ступенях все умножения будутнетривиальными. Это хорошо проследить, анализируя выражение (8.3), описывающее т-юступень алгоритма. Можно показать путем вычислений s (k), что дляопределения отсчетов с номерами k=0, N/4; N/2; 3N/4 не требуетсявыполнение нетривиальных умножений, все умножения здесь на />; />. Ошибка округления этих отсчетов равна нулю. А длявычисления СКЗ ошибок округления отсчетов с другими k воспользуемсявторым выводом. В результате получим
/> (8.6)
Пример 7. Рассмотрим вычисление ошибок округления отсчетов16-точечного БПФ. Для этого выпишем для каждого номера отсчета его двоичноепредставление />. Затем пользуясьвыражением (8.6), вычисляем ошибки округления. Результаты расчетов приводятся втабл.1. Заметим, что отсчеты с номерами k=4, 8, 12 вычисляются точнотак, как s(k)>4, где 4-максимальное число возможных ступенейалгоритма, т. е. в формировании отсчетов с номерами 0, 4, 8, 12 участвуютумножения только на тривиальные множители /> (см. также граф на рис.5).
Вычислим СКЗ ошибки, усредненное по всем k. На первых двухступенях алгоритма (т=1,2) все умножения являются точными, так какмножители тривиальны />; на выходах т-й ступени(т = 3,…, у) появляется N произведений в /> блоках, причем четыре произведения в каждомблоке являются точными в результате умножения также на тривиальные множители.
Таблица 1
/>
Тогда, пользуясь выводом 1, получаем
/>(8.7)
В случае 16-точечного БПФ
/>
Определим СКЗ выходной ошибки алгоритма, обусловленноймасштабированием промежуточных результатов. Ошибка масштабирования комплексногочисла на m-й ступени алгоритма имеет нулевое среднее и дисперсию />. Множитель/> появляется из-за масштабирования на m предыдущих ступенях.
В соответствии с выводом 2 СКЗ индивидуальной ошибки равно
/> (8.8)
в случае 16-точечного БПФ
/>
Усредненное по всем выходам k СК3 также равно
/> (8.9)
Среднеквадратическое значение суммарной ошибки, обусловленнойокруглением и масштабированием, вычисляется как сумма отдельных СК3:
/> (8.10)
СК3 суммарной ошибки, усредненное по всем выходам алгоритма, равно
/> (8.11)
В случае 1б-точечного БПФ
/>
Результаты вычисления СК3 суммарных ошибок также приведены в табл.1.Точность алгоритма принято оценивать по отношению «СК3 ошибки / СК3 выходного сигнала».В данном случае оно составляет
____«СКЗ ошибки»______ = /> (8.12)
«СКЗ выходного сигнала»
В случае 1б-точечного БПФ
___«СК3 ошибки»_______ = />
«СК3 выходного сигнала»
Как видно из полученных выражений, точность алгоритма зависит отдвух параметров: длины преобразования N и разрядности представлениячисел b1.Если известны требования по точности вычисления иразмер преобразования, разрядность процессора БПФ можно определить изследующего выражения, полученного логарифмированием выражения (8.12):
/> (8.13)
Теперь перейдем к анализу ошибок БПФ, вызванных неточнымпредставлением значений поворачивающих множителей /> Пусть
/> -точные значениякоэффициентов Фурье, а />'(k) — неточные,полученные при условии представления коэффициентов конечным числом разрядов
/>.
В рассматриваемом алгоритме каждый элемент /> представляет собой произведение v квантованныхкоэффициентов. Действительно, можно показать, что каждый отсчет входнойпоследовательности, продвигаясь к k-мувыходу (см. рис.5), накаждой ступени алгоритма претерпевает умножение только на один поворачивающиймножитель, т. е.
/> (8.14)
где
/> (8.15)
Дисперсия элементарных ошибок округления /> комплексных коэффициентов равна />
Индивидуальная ошибка БПФ, обусловленная округлениемкоэффициентов, равна
/> (8.16)
Вычитая (8.15) из (8.14), получаем
/> членывысшего порядка.
Пренебрегая членами высшего порядка и подставляя последнеевыражение в (8.16), получаем следующее СК3 ошибки, вызванной квантованиемкоэффициентов:
/> (8.17)
По теореме Парсеваля СК3 выходного сигнала равно
/> (8.18)
Отсюда «СК3 ошибки / СК3 выходного сигнала» = (v/б) 2-2b2.В случае 16-точечного БПФ «СК3 ошибки/СК3 выходного сигнала» = (2/3) 2-2b2.
В реферате рассмотрены вычислительные ошибки только одного из многочисленныхалгоритмов БПФ для определенного вида представления чисел. С таким же успехомиспользованный подход может быть применен для анализа ошибок других алгоритмов.При этом, очевидно, СК3 ошибок будет существенно зависеть от конфигурациинаправленного графа алгоритма, способа представления чисел и методамасштабирования.
ЗАКЛЮЧЕНИЕ
Досередины 1960-х для представления спектрального разложения использовалисьточные формулы для нахождения параметров синусов и косинусов. Соответствующиевычисления требовали как минимум N**2 (комплексных) умножений. Таким образом,даже сегодня высокоскоростному компьютеру потребовалось бы очень много временидля анализа даже небольшого временного ряда (для 8,000 наблюдений потребовалосьбы, по меньшей мере 64 миллиона умножений).
Ситуациякардинально изменилась с открытием так называемого алгоритма быстрого преобразования Фурье, илиБПФ для краткости. Достаточно сказать, что при применении алгоритма БПФ времявыполнения спектрального анализа ряда длины N стало пропорционально N*log2(N)что конечно является огромным прогрессом.
Однаконедостаток стандартного алгоритма БПФ состоит в том, что число данных рядадолжно быть равным степени 2 (т.е. 16, 64, 128, 256,...). Обычно это приводит кнеобходимости добавлять нули во временной ряд, который, как описано выше, вбольшинстве случаев не меняет характерные пики периодограммы или оценкиспектральной плотности. Тем не менее, в некоторых случаях, когда единицавремени значительна, добавление констант во временной ряд может сделатьрезультаты более громоздкими.
СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ
1. Цифроваяобработка сигналов: Учебн. Пособие для вузов/Л. М. Гольденберг, Б.Д. Матюшкин,М. Н. Поляк. – 2изд., перераб. и доп. – М.: Радио и связь, 1990. – 256 с.: ил.
2. РабинерЛ., Гоулд Б. Теория и применения цифровой обработки сигналов. – М.: Мир,1978.-848с.
3. МакклепланД.Х., Рейдер Ч. М. Применение теории чисел в цифровой обработке сигналов: Пер.с англ. – М.: Радио и связь, 1983.-264с.
4. ГольденбергЛ.М., Матюшкин Б.Д., Поляк М. Н. Цифровая обработка сигналов: Справочник.- М.:Радио и связь, 1985. – 312с., ил.