КРАТКИЙ ВЫВОД АЛГОРИТМА КУЛИ–ТЬЮКИ Текст научной статьи по специальности «Математика»

УДК: 517.58 MSC2010: 65T50

КРАТКИЙ ВЫВОД АЛГОРИТМА КУЛИ-ТЬЮКИ © М. С. Беспалов

ВЛАДИМИРСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕрСИТЕТ ИМЕНИ

Александра Григорьевича и Николая Григорьевича Столетовых ИНСТИТУТ прикладной МАТЕМАТИКИ, ФИЗИКИ и информатики ул. Горького, 87, Владимир, 600000, Российская Федерация e-mail: [email protected]

Brief derivation of the Cooley-Tukey algorithm.

Bespalov M. S.

Abstract. In the development of digital information processing are distinguished the emergence of a fast algorithm for the implementation of the discrete Fourier transform (DFT), known as the Fast Fourier Transform (FFT).

The description of this famous algorithm appeared in a paper by Cooley and Tukey, and is a flow of concepts from Good's ideas for the discrete Walsh transform to DFT.

The paper gives the final form of the matrix notation of the FFT algorithm for an arbitrary composite order. It is proposed that the algorithm should start with a reverse permutation rather than include a perfect reshuffle in each step of the algorithm.

The inverse permutation matrix is presented as a b-product of unit matrices (the b-product is a new type of tensor product of matrices introduced earlier by the author). The complete derivation of the FFT algorithm is given.

The methodical aspect of the article is that it can be used in the educational process, the practical significance is in that the results can be used for software development.

The new in the matrix form of FFT is the representation of the L permutation matrix as the b-product of the identity matrices. This allows the algorithm to start with a reverse rearrangement, without being distracted by subsequent permutations. In the algorithms which are presented in the Internet FFT permutation in the form of the perfect shuffle offer to do at every step.

The program and algorithm of the composite order FFT is based on theorem 1 or its consequences, similar to theorem 2. The flowchart for N =16 was represented earlier.

Keywords: Fast Fourier Transform,, tensor product of matrices.

Введение

В качестве революционного скачка в развитии цифровой обработки информации выделяют появление быстрого алгоритма реализации дискретного преобразования Фурье (ДПФ), известного под названием быстрое преобразование Фурье (БПФ).

Описание этого знаменитого алгоритма появилось в статье [1] Кули и Тьюки и является переносом идей, предложенных Гудом [2] для дискретного преобразования Уолша, на ДПФ.

При практическом применении чаще рассматривают ДПФ порядка N = 2П. В статье предлагается короткий (за счет формализации) вывод этого знаменитого алгоритма для общего случая составного порядка N = щп2 ... пк. Матричная форма ДПФ из обзора [3], изложенная в [4], дополнена новым вариантом тензорного произведения матриц, введенным в [5], и новым подходом к доказательству.

В качестве основного аппарата используются два вида тензорного произведения матриц, сравнительная характеристика которых приведена в [6]. Это кронекерово поизведение А < В и 5-произведение А 0 В. Напомним, что кронекерово произведение матриц А = ) и В представляется [7] в виде блочной матрицы с блоками вида а^- ■ В. Также в виде блочной матрицы с блоками в виде произведения столбца первой на строку второй матриц определяется [5] и 5-произведение. Поясним определение с помощью примера с матрицами

Н 0 Е =

V

Н =

1 1 1 1

1 1 1 -1

(10) (01)

)

и

Е=

- i

1

(10) (01)

\

/

1 0

0 1

/ 1 0 1 0 \

1 0 -1 0

0 1 0 1

V 0 1 0 -1 у

(1)

(2)

То есть, в блок-строке меняется столбец первого сомножителя, а в блок-столбце меняется строка второго сомножителя.

В [6] доказано, что каждое из этих двух видов тензорного произведения ассоциативно и для них справедлива

Лемма. Если размеры матрицы таковы, что обычные умножения возможны, то (А ■ В) <8) (С ■ В) = (А <8) С) ■ (В <8) В), (3)

(А ■ В) 0 (С ■ В) = (А 0 С) ■ (В < В) = (С < А) ■ (В 0 В). (4)

Эта лемма применяется для факторизации, рассматриваемой в данной статье, матрицы Фурье, являющейся матрицей обратного ДПФ

Вм =

и Г,

V ) ],к=о

где = ехр

2пг

N.

(5)

1

1

Прямое ДПФ задается матрицей FN или формулой (5) с заменой wN на wN = ш-1. ДПФ вида (5) переводит входной сигнал x, представленный в виде N-мерного столбца, в сигнал на выходе по формуле

y = Fn • x,

что в координатной форме записывается так

N-1

yk = ^ XjwN , k = 0,1,... ,N — 1. (6)

j=0

При анализе ДПФ и в других сферах часто встречается перестановочная матрица идеальной тасовки (perfect shuffles) [8], которая N = nm упорядоченных чисел от 0 до N — 1 переводит из одного лексикографического порядка (относительно двухразрядной по основаниям n и m системы счисления) в другой. А именно, число k + sm переходит в s + kn, где k пробегает от 0 до m — 1, а s пробегает от 0 до n — 1. В докладе [4] эти матрицы обозначены L^n. Легко заметить, что матрица идеальной тасовки выражается через b-произведение единичных матриц

L(n) = 1 (7)1

mn т ^ п'

При умножении этой матрицы на массив x выбираются из x подмассивы объема m с шагом n и помещаются подряд.

1. ПРОСТЕЙШИЙ СЛУЧАЙ БПФ СОСТАВНОГО ПОРЯДКА

Пусть N = пт. Применяя соотношение Ч' = I е N 1 = шт к формуле (6) для

)п

= шт

начальных m координат выходного сигнала, получим

т—1 п—1 т—1 / п—1 \

V, = V V т- ■ ^ V х ■(№] < 3

ук = /_^ТПп+П = I /_^тПп+П I шт •

32=0 л=0 ¿2=0 \^1=0 )

Общий вид номера отсчета возьмем зт + к, где 0 < к < т, 0 < в < п, и, используя сокращение а3- = х3'Ч3 и соотношение ш' = 1, получим:

N-1 N-1 т—1п-1

Ук+вт = ^ Х3ШМ+8т3 = )= "}2а32П+Пш№2П+Л).

j=0 3=0 32=0 31=0

' - 1 т- 1 ' - 1 т- 1

Ук+вт = ^ а32'+31 Ч3 Ч^2 = ^ Ч'31 Чу1 х32'+31 шт2 • (7)

31=0 32 =0 31=0 32=0

Запись данной формулы (7) в матричном виде и составляет основное

Предложение. Матрица ДПФ порядка N = пш факторизуется (раскладывается на множители) в виде

1„т = (Щ 0 1т) ■ ТПт ■ (4 0 1т) ■ (1т 0 1П), (8)

где Ттп — диагональная матрица порядка N = пш с элементами на диагонали ¿¿,.+/т = ш1^, 1П — единичная матрица порядка п.

Поясним этот переход. Внутренняя сумма в (7) равная ^™=о Хгп+л^т2 в матричном виде вычисляется по формуле

(1П 0 1т) ■ (1т 0 1п) ■ X.

Сначала осуществляется перестановка элементов в виде идеальной тасовки, а потом ДПФ порядка ш для каждой из п пачек.

Умножение в (7) на поправочные множители ш^/1 в матричном виде записыва-

Т-г(п)

ется с помощью умножения на диагональную матрицу ТПт с этими элементами на диагонали.

Внешняя сумма в (7) есть ш раз повторенное ДПФ порядка п с выбором элементов с шагом ш и размещением результата вычислений в те же ячейки. Это параллельное вычисление в матричном виде оформляется с использованием операции кронекерова произведения в виде оператора с матрицей 1п 01т. При этом результат вычисления получается в нужном порядке: сначала при в = 0 для всех к от 0 до ш — 1, потом при в = 1 при том же наборе к и т. д.

Формулу (8) можно немного сократить применением свойства (4), до

^пт (1п 0 1т) ' Т(гга ' (1т 0 1п).

Приведенный вывод БПФ отличен от предложенного в [4]. Формула (7) удобнее (8) при составлении алгоритма, а описанный переход полезен при построении блок-схемы в виде стандартных многокрылых «бабочек».

2. Овщий случай БПФ составного порядка

Пусть N = п.... п2п0.

Введем обозначения Ш/ = поп2 ... П/, М/ = П/+1П/+2 ... п^, где ш0 = 1, М^ = 1, N = Ш/М/. Для предыдущего уровня обозначим N = п^_о... п2по и М/ = П/+1П/+2 ... П._1, где ММ^-о = 1.

Теорема 1. Матрица ДПФ порядка N = п^.Пд._о... п2по факторизуется в виде

^ = ■ ■ ... ■ 51 ■ 4 (9)

где = (1М. 0 0 1т,._1) ■ (1М. 0 Т^')), Ь = 1щ 0 1П2 0 ... 0 1Пк.

Доказательство. Формула (8) при п1 = ш, п2 = п составляет базу индукции. Индуктивное предположение состоит в верности формулы

= ■ «s'fc-2 ■ • • • ■ ■ L j (10)

■(Im

Так как N = nkN, то по (8)

где Sj = (Im. 0 Fn. 0 Im.-!) ■ (Ijm. 0 T^')), L = /щ 0 42 0 • • • 0 Infc_i•

fn = (Fnfc 0 iN) ■ TNnk) ■ (4fc 0 fJ) ■ (iN 0 znfc). (11)

Первые два сомножителя (11) составляют Sk. В третий сомножитель подставим (9) и k раз применим формулу (3)

/Пк 0 fJ = (/nfc)k 0 (Sk-1 ■ Sfc-2 ■... ■ Si ■ L) = = (4fc 0 Sk-i) ■ (/nfc 0 Sfc-2) ■... ■ (/nfc 0 Si) ■ (4fc 0 L).

Каждый, кроме последнего, множитель при j от 0 до k — 1 по формуле (3) приводит к

Ink 0 Sj = (4fc)2 0 ((/MM. 0 Fnj 0 /((-) ■ (/д, 0 Tmjj= = (/nk 0 /Mj 0 Fnj 0 /m,,!) ■ (/nk 0 /Mj 0 )) = Sj,

так как 4k 0 /д = /м,.

А последний объединим с последним из (11) и применим (4)

(/nk 0 L) ■ (/n 0 /nk) = (L ■ /j) 0 (/nk)2 = L 0 /nk = L. □

Следствие 1. Основную формулу БПФ (9) можно слегка сократить за счет совмещения предварительной перестановки L с S1 — первым шагом вычислений (формула (2) есть простейший пример этого совмещения):

Si ■ L = Fm 0 /n2 0 ... 0 /nk.

Доказательство. Так как T((n1l) = тПП1) = /ni, то S1 = /Ml 0 Fni. Обозначим Li = /n2 0 ... 0 /nk. По формуле (4)

Si ■ L = (/MI 0 Fni) ■ (/ni 0 Li) = Fni 0 Li.

Следствие 2. При N = nknk-i... n2ni возможен следующий метод факторизации в обратном порядке

fN = LT ■ Si ■ S2 ■... ■ Sk, (12)

где Sj = (/м, 0 )) ■ (/м, 0 Fn, 0 /(j-i), LT = /nk 0 ... 0 /n2 0 /щ.

Доказательство. Относительно операции транспонирования различные виды умножения ведут себя по разному, что для тензорных произведений доказано в [6]:

(Л ■ Б)т = Бт ■ Лт, (Л 0 Б)т = Лт 0 Бт, (Л 0 Б)т = Бт 0 Лт.

Изменим обозначение на Д, = (/м,0/т_1 )■(/м,0ТТт'7 в (9) и протранспонируем

^ = )т = ■ ■... ■ Д1 ■ Ь)т = ьт ■ №)т ■... ■ (5к)т.

Так как матрицы /п и симметричны, то

(д,)т = (/м, 0 ттп))т ■ (/м, 0 ^ 0 /т,_1 )т = .

3. Наиболее распространенный случай БПФ

На практике наибольшую популярность получил случай БПФ порядка N = 2П. Если будем использовать обозначения Н и Е из (1) и сокращения Е, = Е® = /2,, то формулы (9) и (12) перепишутся в виде.

Следствие 3. При N = 2П верны формулы факторизации

^ = 5п ■ 5п-1 ■ ... ■ 51 ■ Ь

где = Еп-,-1 0 ((Н 0 Е,-1) ■ Т®), Ь = Еп0;

^ = Ь ■ 51 ■ ■... ■ 5П,

где = Еп-,-1 0 (Т2(2) ■ (Н 0 Е,-1)), Ь = Еп0.

В [6] доказано, что преобразование с матрицей Ь есть перестановка реверс, подробное описание которой приведено в [9]. Если ввести обозначение х + 1 для операции добавления ко всем элементам массива х единицы, то программа построения перестановки в реверсном порядке множества чисел {0,1,..., N — 1}, где N = 2П, выглядит так:

( 0\ -1 1 / 2х \

х := ; для ? от 1 до п — 1 выполняем х : =

^ ^2х + 1 у

Линейное преобразование с матрицей Н = служит не тольло ДПФ порядка 2, но и примитивным дискретным преобразованием Хаара, и представляет собой пере/ а \ (а +

ход от входного вектора х = к вектору на выходе у = по формуле

\Ч \а — ъ)

умножения матриц у = Н ■ х, что обозначим у = Н(х). Условно это преобразование можно записать в виде перехода

(а, Ъ) —' (а + Ъ, а — Ъ). (13)

Кроме (13) в дальнейшем будем использовать поправленное дискретное преобразование Хаара у = Нд(х) с поправкой Л с условным обозначением

(а, Ъ)д —У (а + ЛЪ, а — ЛЪ).

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

То есть исходный сигнал х при каждом заданном j будем разбивать на 2П-1 непересекающихся подсигнала (дизъюнктивное разбиение).

Если j = 0, то прореживание отсутствует. Тогда определим

хт[0] := (х2т х2т+1) для т € 0 : 2П-1 — 1.

Координаты в подсигнал брали с шагом 20 = 1.

В общем случае прореживания уровня j координаты в подсигнал берем с шагом 2-? и организуем сквозную нумерацию подсигналов так, что сигнал с меньшими координатами получает меньший номер:

хт(,+к [?] := (хт!У+1+к хт(,+1+^+27) для т € 0 : 2П-5-1 — 1, к € 0 : 23 — 1.

Теорема 2. Алгоритмический вид быстрого алгоритма дискретного преобразования Фурье порядка N = 2П состоит из предварительного шага в виде реверсной перестановки и следующих основных шагов быстрого алгоритма:

для j от 0 до п — 1, для т от 0 до 2п--?-1 — 1 и для к от 0 до 2-? — 1 выполняем

ат27+к-[?] := ^^ (ат27 +к-[j]).

Заключение

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

Новым в матричной записи БПФ служит представление перестановочной матрицы Ь в виде Ъ-произведений единичных матриц. Это позволяет начать алгоритм с реверсной перестановки, не отвлекаясь в дальнейшем на перестановки. В представленных в интернете вариантах БПФ перестановку в виде идеальной тасовки предлагают делать на каждом шаге.

Программа и алгоритм БПФ составного порядка строится на основе теоремы 1 или ее следствий по аналогии с теоремой 2. Блок-схема для N =16 представлена в [10, с. 65].

СПИСОК ЛИТЕРАТУРЫ

1. COOLEY, J. W. & TUKEY, J. W. (1965) An algorithm for the machine calculation of complex Fourier series. Math. Comput.. 19 (90). Pp. 297-301.

2. GOOD, I. J. (1958) The interaction algorithm and practical Fourier analysis. J. Royal Stat. Soc.. Ser.B. (20). Pp. 361-372.

3. JOHNSON, J., JOHNSON, R. W., RODRIGUEZ, D. & TOLIMIERI, R. (1990) A methodology for designing, modifying and implementing Fourier transform algorithms on various architectures. Circuits, Systems and Signal Processing. 9 (4). Pp. 449-500.

4. Малоземов, В. Н., Просеков, О. В. Факторизация Кули-Тьюки матрицы Фурье / Избранные главы дискретного гармонического анализа и геометрического моделирования. Часть первая. Изд. 2-е. Под ред. проф. В. Н. Малоземова // СПб.: Изд-во ВВМ. - 2014. - C. 20-29.

MALOZEMOV, V. N. & PROSEKOV, J. V. (2014) Cooley-Tukey factorization of the Fourier matrix. Selected chapters of discrete harmonic analysis and geometric modeling. Part one. Ed. prof. V. N. Malozemov. St.Petersburg: VVM. Pp. 20-29.

5. Беспалов, М. С. Математические методы в информатике и вычислительной технике. В 2-х ч. Ч. 2. Введение в прикладной гармонический анализ / М. С. Беспалов. - Владимир: ВлГУ, 2007. - 244 c.

BESPALOV, M. S. (2007) Mathematical methods in informatics and computer science. Part two. Introduction to applied harmonic ayalysis. Vladimir: VlSU..

6. Беспалов, М. С. О свойствах тензорного произведения матриц // Журнал вычислительной математики и математической физики. — 2014. — 54(4). — C. 547-561.

BESPALOV, M. S. (2014) On the Properties of a New Tensor Product of Matrices. Computation Mathematics and Mathematical Physics. 54 (4). Pp. 531-574.

7. Беллман, Р. Введение в теорию матриц / Р. Беллман. — М.: Наука, . — c.1969368

BELLMAN, R. (1960) Introduction to Matrix Analysis. New York, Toronto, London: MCGRAY-HILL BOOK COMPANY.

8. DIACONIS, P., GRAHAM, R. L. & KANTOR, W. M. (1983) The mathematics of perfect shuffles. Adv. in Appl. Math.. 4 (2). Pp. 175-196.

9. Малоземов, В. Н. Основы дискретного гармонического анализа / В. Н. Малоземов, С. М. Машарский. - СПб.: Лань, 2012. - 304 с. MALOZEMOV, V. N. & MACHARSKI, S. M. (2012) Discrete harmonic analysis basics. St.Petersburg: Lan..

10. Беспалов, М. С. Дискретные и вероятностные модели / М. С. Беспалов. — Владимир: ВлГУ, 2017. - 84 с.

BESPALOV, M. S. (2017) Discrete and probabilistic models. Vladimir: VlSU.