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

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

Сколько операций требуется  на проведение дискретного преобразования Фурье? Посчитав по определению (N раз суммировать N слагаемых), получаем величину порядка N ². Тем не менее, можно обойтись существенно меньшим числом операций.

Наиболее популярным из алгоритмов ускоренного вычисления ДПФ является т.н. метод Cooley-Tukey, позволяющий  вычислить ДПФ для числа отсчетов N = 2 ^k за время порядкаNlog₂ N (отсюда и название - быстрое преобразование Фурье, БПФ). Этот способ чем-то неуловимо напоминает быструю сортировку. В ходе работы алгоритма также проводится рекурсивное разбиение массива чисел на два подмассива и сведение вычисления ДПФ от целого массива к вычислению ДПФ от подмассивов в отдельности.

Изобретение БПФ  привело к потрясающему всплеску популярности преобразования Фурье. Целый  ряд важных задач раньше решался  за время порядка N ², но после проведения преобразования Фурье над исходными данными (за время порядка Nlog₂ N) решается практически мгновенно. Преобразование Фурье лежит в основе цифровых корелляторов и методов свертки, активно используется при спектральном анализе (практически в чистом виде), применяется при работе с длинными числами.

Широко распространено ошибочное мнение о том, что метод Cooley-Tukey - единственный существующий метод  выполнения БПФ, а само БПФ существует только для случая N = 2 ^k. На самом деле это не так - существуют алгоритмы БПФ для любого числа отсчетов. В одномерном случае, рассмотренном в этой статье, метод Винограда позволяет решить задачу для простого числа отсчетов N. Этот же алгоритм может быть легко обобщен на случай, когда N является степенью произвольного простого числа (а не только двойки), а также на случай, когда число N является произведением степеней простых чисел - т.е. N является произвольным числом, чье разложение на простые множители нам известно.

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

Как уже говорилось выше, существуют алгоритмы БПФ для  произвольного числа отсчетов, но наиболее широкое распространение  получил только алгоритм для случая N = 2 ^k, что является существенным ограничением. Почему же это произошло?

Причина этого в  том, что алгоритм, построенный по методу Cooley-Tukey, обладает рядом очень  хороших технологических свойств. Структура алгоритма и его  базовые операции не зависят от числа  отсчетов (меняется только число прогонов базовой операции "бабочка"). Алгоритм легко распараллеливается с использованием базовой операции и конвееризуется, а также легко каскадируется (коэфициенты  БПФ для 2N отсчетов могут быть легко  получены преобразованием коэфициентов двух БПФ по N отсчетов, полученных "прореживанием" через один исходных 2N отсчетов). Алгоритм прост и компактен, не требует  дополнительной оперативной памяти и допускает обработку данных "на месте". Существует целый ряд  оптимизированных именно для этого  алгоритма DSP-процессоров (это одновременно и причина, и следствие).

Всё это и обусловило популярность в инженерно/программистской  среде именно этого алгоритма, и, соответственно, выбора именно 2 ^k отсчетов при использовании БПФ. Правда, попутно это привело к незаслуженному забвению широкими массами альтернативных алгоритмов, некоторые из которых (что следует отметить) требуют меньше вещественных операций на один отсчет, чем алгоритм Cooley-Tukey.

 

Современноеa состояние вычислительной техники характеризуется

большим расширением сферы  цифровой обработки сигналов. Широкое

распространение получают исследования в области цифровых звукозапи-

сывающих и звуковоспроизводящих устройств, цифрового телевидения,

цифровых фильтров различного рода сигналов вплоть до радиосигналов.

Проектирование подобных устройств требует использования  и разработки

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

Одним из путей удовлетворения этих потребностей является выбор

разумно построенных алгоритмов. Вместо того чтобы повышать быстро-

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

миллионов умножений в  секунду, можно для некоторых  задач попытаться

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

ножений в секунду оказалось  достаточно. К настоящему моменту  имеется

хорошо разработанная  теория, позволяющая подойти к  решению задач с

этих позиций. Она хорошо известна теоретикам, специалистам в  данной

области, но на практике инженеры-конструкторы часто пренебрегают ею,

поскольку она пока недоступна в виде единого научного подхода. Инже-

нер-конструктор должен хорошо знать предмет, чтобы выбрать  алгоритм,

нужный в данном конкретном приложении, среди значительного  разнооб-

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

ния Фурье.

 

2.1. Свойства дискретного преобразования Фурье

Некоторые свойства ДПФ играют важную роль в практических во-

просах обработки сигналов. Ниже они будут в основном перечислены, де-

тали будут рассмотрены  только в случае необходимости.

Линейность

Если хр(n) и yp(n) - периодические последовательности (с периодом в

N отсчетов каждая), a Xр(k), Yp(k) - их ДПФ, то дискретное преобразование

Фурье последовательности хр(n) + yp(n) равно Xр(k) + Yp(k). Это положение

справедливо и для последовательностей  конечной длины.

 

Свойства симметрии

Если периодическая последовательность хр(n) с периодом в N отсчё-

тов является действительной, то ее ДПФ Xр(k) удовлетворяет следующим

условиям симметрии:

Re[Xp(k)] = Re[Xp(N - k)];

Im[Xp(k)] = -Im[Xp(N - k)];

½Xp(k)½ = ½Xp(N - k)½;

arg[Xp(k)] = - arg[Xp(N - k)].

Аналогичные равенства справедливы  и для конечной действитель-

ной последовательности x(n), имеющей N-точечное ДПФ X(k). Если ввести

дополнительное условие  симметрии последовательности xp(n), т.е. считать,

что xp(n) = xp(N - n), тo окажется, что Xр(k) может быть только действитель-

ной.

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

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

двух последовательностей, используя свойства симметрии. Рассмотрим

действительные периодические  последовательности хр(n) и ур(n) с перио-

дом в N отсчётов и N -точечными ДПФ Xр(k) и Yр(k) соответственно. Вве-

дем комплексную последовательность zp(n) вида zp(n)= xp(n)+j×yp(n).

Ee ДПФ

( ) j N nk

N

n

Z p k xp n j yp n e (2 / )

1

0

[ ( ) ( )] - p

-

=

= å + × × ;

Zp(k) = Xp(k) + j×Yp(k).

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

Re[Zp(k)] = Re[Xp(k)] - Im[Yp(k)];

Im[Zp(k)] = Im[Xp(k)] + Re[Yp(k)].

Действительные части Xp(k) и Yp(k) cиммeтpичны, а мнимые анти-

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

и вычитания:

[ ( )]

2

Re[ ( )] Re[ ( )]

Re

Z k Z N k

X k p p

p

+ -

= ;

[ ( )]

2

Re[ ( )] Re[ ( )]

Im

Z N k Z k

Y k p p

p

- -

= ;

[ ( )]

2

Im[ ( )] Im[ ( )]

Re

Z k Z N k

Y k p p

p

+ -

= ;

[ ( )]

2

Im[ ( )] Im[ ( )]

Im

Z k Z N k

X k p p

p

- -

= .

– 28 –

Итак, вычисляя одно N-точечное ДПФ, удается преобразовать сразу

две действительные последовательности длиной по N/2 отсчетов. Если эти

последовательности являются еще и симметричными, то число  операций,

необходимых для получения  их ДПФ, можно сократить еще больше.

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

Из выражения (2.2) видно, что  ДПФ является трудоемкой

пpoцeдуpoй, включaющeй (N - 1)2 комплексных умножений и N×(N - 1)

комплексных сложений. Это  обстоятельство сильно ограничивает приме-

нение прямого ДПФ при  цифровой обработке сигналов.

В настоящее время появилось  множество методов ускоренного  вы-

числения дискретного  преобразования Фурье , отличающихся от методов

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

но следует отметить алгоритм быстрого преобразований Фурье (БПФ)[1, 3,

20, 23, 25] , алгоритм Винограда  преобразования Фурье (АВПФ) [1, 20] и

метод двойного преобразования [22].

Основная идея быстрого преобразования Фурье (БПФ) состоит в

том, чтобы разбить исходную N -точечную последовательность на две бо-

лее короткие последовательности, ДПФ которых могут быть скомбиниро-

ваны таким образом, чтобы  получилось ДПФ исходной N-точечной после-

довательности. Далее полученные короткие последовательности аналогич-

но разбивают на еще  более короткие и т.д. В результате прямое ДПФ не-

обходимо выполнять только над r-точечными последовательностями, где

r - основание алгоритма БПФ. Следует отметить, что прямое ДПФ для r = 2

и r = 4 не содержит нетривиальных умножений. Однако, БПФ для N > 4

содержат комплексные  умножения, поcкольку либо перед r-точечным

ДПФ (для алгоритма с  прореживанием по времени), либо после  него (для

алгоритма с прореживанием  по частоте) отсчеты БПФ следует  умножать

на фазовые множители [23]. Таким образом, N-точечное БПФ  будет вы-

полняться за m = logr N этапов (на каждом этапе необходимо выполнить

N/r прямых ДПФ и умножить входные или выходные отсчеты на фазовые

множители).

АВПФ основан на применении теоретико-числовых методов

(полиномиальной алгебры,  китайской теоремы об остатках) при

вычислении ДПФ.

Метод двойного преобразования заключается в том, что сначала  ис-

ходный сигнал подвергается преобразованию Уолша, Хаара или  Адамара

(преобразованиям в базисе  прямоугольных функций, не требующим  не-

тривиальных умножений), а  затем осуществляется переход в  базис Фурье

умножением на определенную матрицу перехода. Этот метод эффективен

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

Несмотря на обилие методов  преобразования Фурье [1, 23], в на-

– 29 –

стоящее время не существует универсального алгоритма ДПФ, отличаю-

щегося во всех случаях  минимальной трудоемкостью. Далее  будет рас-

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

синтезированному в работе [18]. Этот метод позволяет решать задачу ДПФ

с приемлемым быстродействием.

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

и прореживанием по времени

Если N - четное число, а исходная N-точечная последовательность

разбита на две (N/2) -точечные последовательности, то для вычисления ис-

комого N-точечного ДПФ потребуется около (N×/2)2×2 = N 2 /2 комплексных

умножений, т.е. в два раза меньше по сравнению с прямым вычислением,

здесь множитель (N/2)2 дает число умножений, необходимое для прямого

вычисления (N/2)-точечного ДПФ, а множитель 2 соответствует двум

ДПФ, которые должны быть вычислены. Эту операцию можно повторить,

определяя вместо (N/2)-точечного ДПФ два (N/4)-точечных ДПФ (предпо-

лагая, что ((N/2) - четное) и сокращая тем самым объем вычислений еще в

два раза. Выигрыш в два  раза является приближенным, поскольку  не учи-

тывается, каким образом  из ДПФ меньшего размера образуется искомое

N-точечное ДПФ.

Проиллюстрируем описанную  методику для N-точечной последова-

тельности {x(n)} считая, что N равно степени 2. Введем две (N/2)-точечные

последовательности {x1(n)} и {x2(n)} из четных и нечетных членов x(n) со-

ответственно, т.е. x1(n) = x(2n), n = 0, 1, …, (N /2) – 1; x2(n) = x(2n + 1),

n = 0, 1, …, (N /2) – 1.

При этом N-точечное ДПФ последовательности {x(n)} можно запи-

сать как

( ) å å å å

-

=

+

-

=

-

=

-

=

= + = + +

/ 2 1

0

(2 1)

/ 2 1

0

2

1

0

2

1

0

1( ) ( ) (2 ) (2 1)

N

n

n k

N

N

n

nk

N

N

n

nk

N

N

n

nk

X k х n WN х n W х n W х n W

С учетом того, что [ ( ) ] ( ( )) / 2 ,

2 2 / 2 2 / / 2

N

j N j N

WN = e = e =W - p - p

получим

() ( ) ( )

() ( ) ( ) ïþ

ïý

ü

= + ×

= å × + å ×

-

=

×

-

=

×

,

;

1 2

/ 2 1

0

2 / 2

/ 2 1

0

1 / 2

X k X k W X k

X k x n W W x n W

k

N

N

n

n k

N

N

n

k

N

n k

N , (2.6)

где X1(k) и X2(k) равны (N /2) -точечным ДПФ последовательностей x1(n) и

x2(n).

Из формул (2.6) следует, что N-точечное ДПФ X(k) может быть раз-

ложено на два (N/2)-точечных ДПФ, результаты которых объединяются.

Если бы (N/2)-точечные ДПФ вычислялись обычным способом, то для вы-

числения N-точечного ДПФ потребовалось бы, очевидно, (N 2/2 + N) ком-

плексных умножений. При  больших N (когда N 2/2 >> N) это позволяет со-

– 30 –

кратить время вычисления на 50 %.

Поскольку X(k) рассчитано при 0 £ k < N, a X1(k) и X2(k) рассчитаны

при 0 £ k £ N /2 - 1 , то необходимо доопределить формулу (2.6) для

k ³ N/2. Это определение очевидно и может быть записано следующим

образом:

( )

( ) ( )

ï ïî

ï ïí

ì

- £ £ ÷ø

ö

çè

+ × æ - ÷ø

ö

çè

æ -

+ × £ £ -

=

k N 1.

2

, N

2 2

1;

2

, 0 k N

1 2

1 2

X k N W X k N

X k W X k

X k

k

N

k

N

(2.7)

На рис. 2.1 с помощью  направленного графа представлена последо-

вательность операций при  выполнении восьмиточечного ДПФ  с использо-

ванием двух четырехточечных  преобразований. Входную последователь-

ность х(n) сначала разбивают на две последовательности х1(n) и x2(n) из

четных и нечетных членов x(n), после чего paccчитывaют их преобразова-

ния X1(k) и X2(k). Затем в соответствии с формулой (2.7) получают X(k).

3

W8

2

W8

1

W8

0

W8

x1(0) = x(0)

x1(1) = x(2)

x1(2) = x(4)

x1(3) = x(6)

x2(0) = x(1)

x2(1) = x(3)

x2(2) = x(5)

x2(3) = x(7)

Четырех-

точечное

ДПФ

Четырех-

точечное

ДПФ

X1(0)

X1(1)

X1(2)

X1(3)

X2(0)

X2(1)

X2(2)

X2(3)

X(0)

X(1)

X(2)

X(3)

X(4)

X(5)

X(6)

X(7)

a Wk 8 ×

a

b

a + b

a - b

a

Рис.2.1. Вычисление восьмиточечного  дискретного преобразования Фурье

через два четырехточечных  дискретных преобразования Фурье

Рассмотренная схема вычислений может быть использована для рас-

чета N/2 -точечных ДПФ в соответствии с формулами (2.6) и (2.7). Каждая

из последовательностей х1(n) и х2(n) разбивается на две поcледовательно-