Быстрое преобразование Фурье

цы размерности nk (n'k). Векторы произведения размерности n'k (nk) запи-

сываются на место исходных в m-мерную таблицу.

Отсюда следует, что минимальный  размер m-мерной таблицы дан-

ных N' должен быть равен произведению чисел строк матриц предсложе-

ний (чисел столбцов матриц постсложений) базовых алгоритмов

' .

1

0

' Õ-

=

=

m

k

N nk

При обработке m-мерной таблицы данных m-мерной таблицей Bb

происходит поэлементное умножение таблиц с записью результатов  в

m-мерной таблице данных.

При синтезе обобщенного  гнездового алгоритма целесообразно

m-мерные таблицы представлять одномерными, поскольку m – априорно

неизвестный параметр алгоритма. Соответствие между одномерным l и

многомерными lk, k = 0, ..., m - 1 индексами таблицы F можно записать

(2.47):

[] [ ] å Õ

-

=

-

=

º - = + ×

1

1

1

0

'

0 , 1,..., 1 , 0

m

i

i

j

F l F l l lm l l li nj , (2.47)

где n'k – размерность k-го индекса.

Если размерности индексов m-мерной таблицы выбрать с запасом,

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

таблицы (2.47) будут только тривиальные умножения.

Введем обозначения:

m – количество вложений в схеме алгоритма;

– 55 –

F[] – m-мерная таблица для хранения комплексных данных. В начале

алгоритма в нее записываются (с переиндексацией на входе) исходные

данные;

r1[], r2[] – рабочие массивы для хранения векторов вещественных

или комплексных данных;

P[k], k = 0, ..., m - 1 – массив, элементы которого определяются

[ ] Õ-

=

=

1

0

'

k

j

P k nj для k > 0, P[0] = 1;

L – массив целых чисел размерности m для хранения многомерных

индексов m-мерной таблицы, L[] = (l0, ..., lm-1);

nk, k = 0, ..., m - 1 – вложения алгоритма или число столбцов в матри-

цах предсложений (или число  строк в матрицах постсложений) базовых

алгоритмов;

n'k, k = 0, ..., m - 1 – число строк в матрицах предсложений (число

столбцов в матрицах постсложений) базовых алгоритмов;

nak, k = 0, ..., m - 1 – рабочий массив целых чисел;

m1, m2 – рабочие переменные целого типа для индексации m-мерных

таблиц как одномерных;

f (L) – функция для вычисления одномерного индекса m-мерной таб-

лицы по массиву многомерных  индексов L (см. (2.47));

Q(L, k) – процедура изменения массива многомерных индексов L,

обеспечивающих перебор  всех векторов m-мерной таблицы, параллельных

k-му измерению.

Переиндексация на входе  записывается в виде

F[l0, ..., lm-1] := x[l], (2.48)

lk = l (mod nk), k = 0,...,m - 1; l = 0, ..., N - 1,

где x[l], l = 0, ..., N - 1 – входная последовательность ДПФ.

Переиндексацию на выходе можно представить в виде

X[l] := F[l0 ,..., lm-1]; (2.49)

( N)

n

l l N

m

k k

k mod

1

0

å-

=

×

= ;

lk = 0, ..., nk -1; k = 0, ..., m -1,

где X[l], l = 0,...,N - 1 – коэффициенты дискретного спектра Фурье.

В принятых обозначениях приведем алгоритм формирования

m-мерной таблицы Bb из диагональных матриц Bk, k=0, ..., m - 1 базовых

алгоритмов ДПФ Винограда

L[i]:= 0,i = 0,...,m -1;

( )

[ ] [ ]

( )

0,..., 1;

,0 ;

1 : , , 1: 1 1; 0,..., 1;

1: ;

'0

'

'0

0 = -

ïþ

ïý

ü

= = + = -

=

n

j N

Q L

Bb m B l l m m l n

m f L

– 56 –

[ ]

( )

[ ] [ ] [ ]

[ ]

( )

0,..., 1.

0,..., ,

, ;

1: 1 ; 0,..., 1;

1 : 1 , ,

1: ;

: 0; 0,..., 1;

'

'

'

= -

ï ï ï

þ

ï ï ï

ý

ü

=

ï ï

þ

ï ï

ý

ü

= + = -

= ×

=

= = -

k m

n

j N

Q L k

m m P k l n

Bb m Bb m B l l

m f L

L i i m

k k

k

Сам обобщенный гнездовой  алгоритм БПФ по схеме Винограда  име-

ет вид

ПЕРЕИНДЕКСАЦИЯ НА ВХОДЕ (2.48)

УМНОЖЕНИЕ МАТРИЦ Ak НА ВЕКТОРЫ m-МЕРНОЙ ТАБЛИЦЫ

1) формирование m-мерной таблицы Bb;

2) переиндексация на входе  последовательности x[l], l = 0,K, N -1 в

соответствии с выражением (2.48);

3) обработка m-мерной таблицы F матрицами предсложений

[ ]

( )

[] [ ] [ ]

[ ] [ ] [ ]

( )

0,..., 1;

0,..., 1;

: ;

, ;

2 : 2 , 2 : 2 ; 0,..., ;

2 : 1;

1 : 1 , 1: 1 ; 0,..., 1;

1: ; 2 : 1;

: 0, 0,..., 1;

: , 0,..., 1;

'

'

= -

ï ï ï ï ï

þ

ï ï ï ï ï

ý

ü

= -

ï ï ï ï þ ï ï  ï ï

ý

ü

= ×

= = + =

= ×

= = + = -

= =

= = -

= = -

k m

n

j N

n

n

N N

Q L k

F m r l m m P k l n

r A r

r l F m m m P k l n

m f L m m

L i i m

na n k m

k

k

k

k

k

k

k k

4) поэлементное умножение m-мерных таблиц

F[l]:= F[l]× Bb[l],l = 0,..., N' -1;

5) обработка m-мерной таблицы F матрицами постсложений

– 57 –

[ ]

( )

[] [ ] [ ]

[ ] [ ] [ ]

( )

0,..., 1;

0,..., 1;

: ;

, ;

2 : 2 , 2 : 2 ; 0,..., 1;

2 : 1;

1 : 1 , 1: 1 ; 0,..., 1;

1: ; 2 : 1;

: 0, 0,..., 1;

: ; 0,..., 1;

'

'

'

'

= -

ï ï ï ï ï

þ

ï ï ï ï ï

ý

ü

= -

ï ï ï ï

þ

ï ï ï ï

ý

ü

= ×

= = + = -

= ×

= = + = -

= =

= = -

= = -

k m

n

j N

n

n

N N

Q L k

F m r l m m P k l n

r C r

r l F m m m P k l n

m f L m m

L i i m

na n k m

k

k

k

k

k

k

k k

6) переиндексация на выходе (2.49)

Примеры использования рассмотренного алгоритма для длины блока

данных N = 6 приведены в прилож. 2 (пример реализации гнездового алго-

ритма ДПФ Винограда) и  прилож. 3 (пример реализации гнездового алго-

ритма ДПФ Винограда в  терминах кронекеровских произведений матриц).

Трудоемкость обобщенного  гнездового алгоритма БПФ в вещест-

венных умножениях для  вещественных исходных данных M и сложениях A

можно представить

Õ -

=

=

1

0

m

k

M Mk , (2.50)

,

... ... ...

... ...

1

2

0

2

1

1

1

1

0

1

1

0

0 1 2 3 1 0 1 2 1

0 1 2 1 0 0 2 1

-

-

=

-

=

-

= +

-

=

-

=

- - -

- -

× ÷

÷

ø

ö

ç ç

è

æ

+

ú úû

ù

ê êë

é

÷ ÷

ø

ö

ç ç

è

æ

× × ÷

÷

ø

ö

ç ç

è

æ

= × +

+ × × × × × + + ×  × × =

= × × × × + × ×  × × +

Õ å Õ Õ Õ m

m

j

j

m

i

m

j i

i j

i

j

j

m

j

j

m m m

m m

A n M A n M A

M M A n n M M M A

A A n n n M A n n

(2.51)

где Ak , k = 0,...,m -1 – число сложений в k-м базовом алгоритме Винограда

приведено в [1, 25]; Mk , k = 0,...,m -1 – число умножений в k-м базовом ал-

горитме Винограда приведено  в работах [1, 25].

Если исходные данные комплексные, то для реализации обобщенно-

го гнездового алгоритма  БПФ потребуется всего в два  раза больше веще-

ственных умножений и  сложений.

Из приведенной записи гнездового алгоритма следует, что  его можно

выполнять конвейерно, разбив на следующие этапы:

1) переиндексация на входе;

2) умножение m-мерной таблицы  на матрицы предсложений базовых

алгоритмов;

3) поэлементное умножение m-мерной таблицы на m-мерную табли-

цу Bb;

– 58 –

4) умножение m-мерной таблицы на матрицы постсложений базовых

алгоритмов;

5) переиндексация на выходе.

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

реализации каждого из 5 этапов было приблизительно одинаковым. Оче-

видно, трудоемкость этапов конвейеризации различна. Поэтому каждый

этап необходимо реализовывать  различным количеством процессоров.

Рассмотрим, каким образом  можно распараллелить выполнение эта-

пов 2 – 4 на L процессоров. Наиболее эффективно (в смысле быстродейст-

вия) поставленную задачу реализовать  на основе системы с разделяемой L

входовой оперативной  памятью для хранения m-мерной таблицы данных.

Программное обеспечение  и матрицы базовых алгоритмов можно хранить

как в разделяемой, так  и локальной памяти каждого процессора. Предпо-

ложим, нам следует вычислить N-точечное ДПФ

.

1

0 Õ-

=

=

m

k

N nk

Этап 2 выполняется следующим  образом. При реализации k-го вло-

жения каждый из процессоров  умножает матрицу предсложений Аk на оп-

ределенные N/nk /L векторов m-мерной таблицы F, а результирующие век-

торы записывает на место  исходных в разделяемую память. Наиболее эф-

фективной реализация алгоритма  получается, когда векторы m-мерной

таблицы можно разделить  между L процессорами поровну.

При выполнении этапа 3 каждый из процессоров поэлементно умно-

жает N'/L определенных элементов m-мерных таблиц F и Bb и результат

записывает в таблицу F, хранящуюся в разделяемой памяти.

Этап 4 выполняется аналогично этапу 2. Только вместо матриц пред-

сложений Ak в алгоритме используются матрицы постсложений Ck и цикл

по k проходит от m - 1 до 0 с декриминированием k.

При организации таким  образом мультипроцессорной системы  осо-

бенные требования предъявляются  к L входовой разделяемой оперативной

памяти системы. Для эффективного функционирования системы в смысле

быстродействия необходимо, чтобы быстродействие L входовой памяти в

L раз превосходило быстродействие каждого из L процессоров.

Возможны и другие организации  мультипроцессорных систем для

реализации гнездовых  алгоритмов. Однако предлагаемая система  обеспе-

чивает в ряде случаев  увеличение быстродействия гнездового алгоритма в

L раз по сравнению с однопроцессорной системой.

– 59 –

2.9. Алгоритм быстрого преобразования Фурье двумерных сигналов

Рассмотрим двумерное N1 ´ N2-точечное ДПФ

= 0, ..., 1; = 0, ..., 1;

,

2 1

1

0

1

0

, ,

1 2

1 2

- -

= å å × ×

-

=

-

=

× ×

l N k N

X x W W

N

n

N

m

m lN

k n

k l n m N (2.52)

где W e j 2 / N , j = -1

N

= - × p .

Традиционный построчно-столбцевой метод [23 и др.] состоит в вы-

числении Xk,l как N2 одномерных ДПФ по индексу l и N1 одномерных ДПФ

по индексу k. В этом методе Xk,l отображается в N2 N1-точечных и в

N2 N1-точечных ДПФ. Введем обозначения:

My1(Ni), i = 1, 2 – количество умножений в Ni-точечном одномерном

алгоритме БПФ; Nl1(Ni), i = 1, 2 – количество сложений в Ni-точечном од-

номерном алгоритме БПФ; My2(N1 ´ N2) – количество умножений в

(N1 ´ N2)-точечном двумерном алгоритме БПФ; Nl2(N1 ´ N2) – количество

сложений в (N1 ´ N2)-точечном двумерном алгоритме БПФ.

Тогда сложность двумерного ДПФ можно представить выражениями

(2.53)

My2(N1 ´ N2) = N1 · My1(N2) + N2 · My1(N1),

Nl2(N1 ´ N2) = N1 · Nl1(N2) + N2 · Nl1(N1). (2.53)

При этом важно, чтобы одномерные БПФ реализовывались эффек-

тивными в смысле трудоемкости алгоритмами. Если N1 = N2 = N, то слож-

ность двумерного ДПФ записывается (2.54)

My2(N ´ N) = 2 ×N×My1(N),

Nl2(N ´ N) = 2 ×N ×Nl1(N). (2.54)

Когда My1(N) < 2 × N, более эффективно вычислять Xk,l с помощью

гнездового алгоритма  Винограда (см. подразд. 2.3). В этом случае Xk,l вы-

числяется как N-точечное ДПФ, в котором каждый скаляр заменяется век-

тором из N точек и каждое умножение заменяется N-точечным ДПФ

[1, 25]. В этом случае  для N1 = N2 = N сложность N ´ N точечного ДПФ за-

писывается выражениями (2.55)

My2(N ´ N) = [My1(N)]2,

Nl2(N ´ N) = [N + My1(N)] × Nl1(N), (2.55)

а для N1 ¹ N2 выражениями (2.56)

My2(N1 ´ N2) = My1(N1) ×My1(N2),

Nl2(N1 ´ N2) = N1 ×Nl1(N2) + My1(N2) ×Nl1(N1). (2.56)

Реализация построчно-столбцевого  метода не вызывает затруднений,

если известна реализация эффективного одномерного алгоритма  БПФ.

– 60 –

Реализация гнездового алгоритма  для двумерных сигналов освещена

в литературе [1 и др.] только в терминах кронекеровского произведения

матриц. Недостаток метода заключается  в том, что кронекеровские произ-

ведения матриц базовых алгоритмов БПФ для большого класса задач, об-

ладающих практической ценностью, имеют большую размерность и  соот-

ветственно требует большой  памяти для их хранения.

Ниже представлен гнездовой  алгоритм БПФ двумерных сигналов,

свободный от указанного недостатка. Основная идея алгоритма заключает-

ся в следующем. Двумерная  таблица данных обрабатывается матрицами

предсложений и постсложений базовых алгоритмов БПФ, отличающихся

малой размерностью. Из диагональных матриц базовых алгоритмов БПФ

формируется двумерная таблица  размерности My1(N1) × My1(N2) способом,

аналогичным формированию многомерной  таблицы (m = 2) Bb для одно-

мерных сигналов. Только вместо базовых алгоритмов Bk, k = 1, 2 исполь-

зуются одномерные представления m1 и m2-мерных таблиц Bb одномерных

гнездовых алгоритмов БПФ. Алгоритм записывается следующим образом:

1) переиндексация столбцов  исходной двумерной таблицы данных

(по алгоритму переиндексации  на входе одномерного алгоритма);

2) обработка столбцов  двумерной таблицы данных матрицами  пред-

сложений для размерности N1 с записью результирующих векторов на ме-

сто исходных столбцов (как  в одномерном алгоритме);

3) переиндексация строк  двумерной таблицы данных (по  алгоритму

переиндексации на входе  одномерного алгоритма);

4) обработка строк двумерной  таблицы данных матрицами предсло-

жений для размерности N2 с записью результирующих векторов на место

исходных строк (как в  одномерном алгоритме);

5) поэлементное умножение  таблицы данных матрицами на  на дву-

мерную таблицу, сформированную из диагональных матриц базовых алго-

ритмов, с записью произведений в двумерную таблицу данных;

6) обработка строк двумерной  таблицы данных матрицами постсло-

жений для размерности N2 с записью результирующих векторов на место

исходных строк (как в  одномерном алгоритме);

7) переиндексация строк  двумерной таблицы данных (по  алгоритму

переиндексации на выходе одномерного алгоритма);

8) обработка столбцов  двумерной таблицы данных матрицами  пост-

сложений для размерности N1 с записью результирующих векторов на ме-

сто исходных столбцов (как  в одномерном алгоритме);

9) переиндексация столбцов  исходной двумерной таблицы данных