РЕФЕРАТ
по дисциплине «Анализ временных рядов»
на тему:
«Алгоритм фильтрации, пример на основе БПФ
»
СОДЕРЖАНИЕ
Введение
1. Основы алгоритмов БПФ
2. Алгоритм БПФ с прореживанием по времени
3. Программа и пример реализации алгоритма БПФ с прореживанием по времени
4. Алгоритм БПФ с прореживанием по частоте
5. Применение метода БПФ для вычисления обратного ДПФ (ОДПФ)
6. Применение БПФ для вычисления реакции ЦФ
7. Другие быстрые алгоритмы вычисления дискретного преобразования Фурье
7.1 Обобщенный алгоритм Кули-тьюки с произвольным основанием с множителями поворота
7.2 Алгоритм простых множителей
7.3 Алгоритм винограда
8. Анализ точности реализации алгоритмов БПФ
Заключение
Список использованных источников
ВВЕДЕНИЕ
В основе преобразования Фурье (ПФ) лежит чрезвычайно простая, но исключительно плодотворная идея – почти любую периодическую функцию можно представить суммой отдельных гармонических составляющих (синусоид и косинусоид с различными амплитудами A
, периодами Т
и, следовательно, частотами ω).
Неоспоримым достоинством ПФ является его гибкость – преобразование может использоваться как для непрерывных функций времени, так и для дискретных.
ПФ часто применяется при решении задач, возникающих в теории автоматического регулирования и управления, в теории фильтрации и т.д. Разберем один из примеров. Имеется некий линейный фильтр – изготовленный то ли в виде набора спаянных между собой резисторов, конденсаторов и катушек индуктивности, то ли в виде модульной конструкции интегральных микросхем. Известен также входной сигнал (на рис. 1 в качестве входного сигнала изображена дельта-функция, то есть импульс исчезающе короткой длительности и бесконечно большой амплитуды). Необходимо определить, какой сигнал появится на выходе нашего фильтра.
Рисунок 1 - Исследование линейного фильтра
Ход решения этой задачи зависит от того, какую позицию мы предпочтем. Выберем временной путь решения (верхняя половина рис. 2) – придется входной сигнал записать как функцию времени S
BX
(t
) и использовать импульсную характеристику фильтра h
(t
), то есть математическую запись его работы во времени. Отправимся по частотному пути (нижняя половина рис. 2) – нужно будет оперировать уже не с самим входным сигналом, а с его спектром g
bx
(ω). Δа и алгоритм работы нашего фильтра потребуется представить в частотной области – в виде частотной характеристики 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=log2
N,
а число требуемых арифметических операций для вычисления N-точечной ДПФ будет порядка Nv
, т.е. уменьшается примерно в N/log2
N раз. Так, при 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 = log2
N,
где 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).
Полный алгоритм БПФ с прореживанием по частоте и его программная реализация аналогичны рассмотренным выше для метода БПФ с прореживанием по времени.
Необходимо отметить, что в обоих алгоритмах БПФ - и с прореживанием по времени, и с прореживанием по частоте требуется примерно N
log2
N
операций (комплексных умножений) и оба алгоритма могут быть реализованы по способу с замещением, используя только один массив ячеек памяти. В обоих алгоритмах должна быть предусмотрена процедура двоичной инверсии - на входе (при прореживании по времени) или на выходе (при прореживании по частоте).
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=N1
N2
,
где N
- длина преобразования. Очевидно, приводимые ниже положения легко могут быть перенесены на более общий случай, когда
Обобщенный алгоритм Кули-Тьюки с произвольным основанием с множителями поворота. Итак, пусть N=N1
N2
,
гдеN1
и N2
- положительные целые. Покажем, что в этом случае вычисление исходного N-точечного ДПФ можно свести к вычислению N1
N
2
-точечных и N2
N1
-точечных ДПФ и N
умножениям на множители поворота . Для этого в выражения для ДПФ (1.1)
(7.1)
где , необходимо сделать подстановку:
k=k1
+k2
N2,
k1
=0, …., N2
-1; k2
=0, …., N1
-1;
(7.2)
n=n1
+n2
N2,
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; n1
k2
=0,
1,2; n2
k1
=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;
для выходной последовательности
(k1
N2
+k2
)=X((s1
k2
N1
+ s2
k1
N2
)mod N),
(7.6)
k1
=0,..., N1
-1; k2
=0,..., N2
-1,
где запись п
mod N
означает «остаток от деления п
на N»,
а s1
и s2
определяются из следующих уравнений в соответствии с китайской теоремой об остатках о восстановлении целого числа по его вычетам: s1
N1
= 1mod N2
,
s1
<N2
, s2
N2
== 1modN1
, s2
<N1
.
Тогда алгоритм N=N1
N2
- точечного
ДПФ представляется в виде
Таким образом, алгоритм простых множителей (АПМ) является способом представления одномерного ДПФ в виде многомерного, причем размерность зависит от числа взаимно простых сомножителей N.
Алгоритм простых множителей имеет ступенчатую форму объединения мало точечных преобразований. В данном случае на первой ступени производится N
1
N
2
-точечных ДПФ, а на второй ступени - 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(kv
2v-1
+... +k1
)=Xv
(kv
2v-1
+kv-1
2v-2
+.
.. +k1
).
Представим индекс элемента т-й
ступени в виде
(n1
,
..., nv-m
,km
,,.
.., k1
)=2v-1
n1
+2v-2
n2
+…
+2m
nv-m
+2m-1
km
+…+k1
=2m
l+q
(8.4)
Тогда число Xm(2ml+q)
можно рассматривать как q-й
элемент l
-го блока т-й
ступени.
Пример б. Рассмотрим описанную процедуру синтеза алгоритма БПФ с прореживанием по времени на примере 16-точечного ДПФ. В этом случае v=4. Индексы n
и k
представим следующим образом: n=n4
23
+n3
22
+n2
2+n1
, k=k4
23
+k3
22
+k2
2+k1.
Подставляя n
и k
в выражение для 16-точечного ДПФ, получаем
X(k4
k3
k2
k1
)=
(n4
nЗ
n2
n1
)
Теперь распишем алгоритм по ступеням:
т
= 0 - инверсия входной последовательности:
Х0
(n1
23
+n2
22
+n3
2+n4
)=x(n4
23
+n3'
22
+n2
2+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-2b
2
.
В случае 16-точечного БПФ «СК3 ошибки/СК3 выходного сигнала» = (2/3) 2-2b
2
.
В реферате рассмотрены вычислительные ошибки только одного из многочисленных алгоритмов БПФ для определенного вида представления чисел. С таким же успехом использованный подход может быть применен для анализа ошибок других алгоритмов. При этом, очевидно, СК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с., ил.