Алгоритм Винограда для расчета циклической свертки

Содержание

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

Распространяется под лицензией LGPL v3

Страница проекта на SourceForge libdspl-2.0 on SourceForge

Обнаружили ошибку? Выделите ее мышью и нажмите ctrl+enter
Введение

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

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

В данном разделе мы рассмотрим другой подход к сокращению вычислений при расчете циклических сверток: алгоритм Винограда, который обеспечивает минимально возможное число умножений при расчете циклической свертки. При этом наибольшее значение имеют алгоритмы коротких циклических сверток, которые, в свою очередь, позволяют построить максимально эффективные алгоритмы быстрого преобразования Фурье.

Вывод алгоритмов Винограда довольно громоздкий, но мы рассмотрим примеры данных алгоритмов для циклических сверток длины N = 2, N = 3 и N = 4, которые покажут методику и приемы проектирования алгоритмов Винограда циклических сверток произвольной длины.

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

В предыдущем разделе мы определили циклическую свертку двух последовательностей a(n) и b(n) равной длительности N отсчетов, n = 0\ldots N-1 как:

equation 1
(1)
которую можно представить в матричной форме:

equation 2
(2)
Квадратная матрица \mathbf{B} содержит только N различных элементов, причем каждая строка (столбец) матрицы получается путем циклического сдвига предыдущей строки (столбца) матрицы на один отсчет. Матрица такой структуры называется циркулянтом или циркуляционной матрицей.

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

Результат циклической свертки также может быть представлен через полиномиальное произведение вида [1]:

equation 3
(3)
где

equation 4
(4)
equation 5
(5)
equation 6
(6)
полиномы степени N-1, содержащие по N коэффициентов.

Таким образом, произведение полиномов A(z) и B(z) берется по модулю z^N-1, т.е. C(z) ничто иное как остаток от деления A(z)B(z) на z^N-1.

В результате C(z) в выражении (3) представляет собой полином степени N-1, чьи коэффициенты равны циклической свертке коэффициентов A(z) и B(z).

Само по себе выражение (3), на первый взгляд, не упрощает расчет циклической свертки, однако исследование модулярной алгебры, которое было проведено в 70-х годах XX века привело к разработке алгоритмов расчета циклической свертки с минимально возможным числом операций умножения.

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

В 1976 году Шмуэль Виноград доказал следующую теорему [2].

Теорема (Винограда) Количество операций умножения необходимое для расчета N-точечной циклической свертки согласно (3) равно 2N - k, где k — число неприводимых многочленов z^N-1.

Вычисление циклической свертки при N = 2 требует две операции умножения, так как 2N = 4 а многочлен z^2 - 1 = (z-1)(z+1) представляется произведением двух неприводимых многочленов.

Напомним, что неприводимым многочленом называют такой многочлен, который делится без остатка только на себя. Легко видеть, что многочлены z - 1 и z+1 являются неприводимыми.

Полиномы z^N-1 при N>1 не являются неприводимыми и могут быть разложены в произведение k неприводимых многочленов Q_i(z):

equation 7
(7)
Тогда, в соответствии с китайской теоремой об остатках, полином C(z) может быть представлен

equation 8
(8)
где

equation 9
(9)

equation 10
(10)

equation 11
(11)
а полиномы S_i(z) находятся из алгоритма Евклида как:

equation 12
(12)
а полином R_i(z) удовлетворяет условию:

equation 13
(13)
Заметим, что полиномы Q_i(z), T_i(z), R_i(z) и S_i(z) не зависят от входных полиномов A(z) и B(z), а только от величины N.

Правила деления полиномов

Для того чтобы получить алгоритм Винограда для циклической свертки необходимо научиться рассчитывать остаток от деления полиномов. Проще всего это сделать путем деления полиномов уголком, как это показано на рисунке 1.

Пример деления полиномов
Рисунок 1. Пример деления полиномов

Пример демонстрирует процесс деления полинома 4z^3 + 3z^2 - z + 2 на z^2 - z + 1. В результате частное равно 4z + 7, а остаток от деления равен 2z  -5. Тогда, можно записать, что

equation 14
(14)
а сам полином можно представить в виде:

equation 15
(15)
Рассмотрим еще один пример. Пусть требуется рассчитать Y(z) = (8z^3 + 2z -1) \Mod (z^2 - 2z). Это также можно посчитать через деление уголком, как это показно на рисунке 2 (отсутствующий квадратичный член 8z^3 + 2z -1 добавлен как 0 z^2).

Пример деления полиномов
Рисунок 2. Пример деления полиномов

Таким образом

equation 16
(16)
Заметим, что если порядок полинома X(z) меньше порядка полинома Y(z), то

equation 17
(17)
Для расчета полиномов A_i(z), B_i(z) и других полезно вспомнить теорему Безу [3, стр. 45], которая говорит о том, что

equation 18
(18)
т.е. остаток от деления полинома X(z) на z - g, равен значению полинома X(z) при z = g.

Циклическая свертка N=2

Применение вышеизложенных формул рассмотрим на простейшем примере для N=2. Тогда

equation 19
(19)
представляется произведением двух неприводимых многочленов Q_1(z) = z-1 и Q_2(z) = z+1. Входными данными для расчета циклической свертки являются последовательности [a_0, \,\,\,a_1] и [b_0, \,\,\, b_1]. Тогда полином A(z) = a_0 + a_1 z, соответственно B(z) = b_0 + b_1 z и в соответствии с (3) и (8) можно записать:

equation 20
(20)
Шаг 1. Рассчитываем полиномы A_i(z) и B_i(z), использую теорему Безу (18) и выражение (10):

equation 21
(21)
equation 22
(22)
Аналогично используя (11):

equation 23
(23)
equation 24
(24)
Шаг 2. Формируем полиномы C_i(z) согласно (9):

equation 25
(25)
equation 26
(26)
Здесь s_1 и s_2 — промежуточные суммы, а m_1 и m_2 — промежуточные умножения.

Расчет полиномов C_1(z) и C_2(z) потребовал две операции умножения.

Шаг 3. Рассчитываем полиномы T_i(z) используя (12):

equation 27
(27)
equation 28
(28)
Шаг 4. Расчет полиномов R_i(z) из уравнения (13).

equation 29
(29)
equation 30
(30)
Шаг 5. Формируем полиномы S_i(z) согласно (12) на основе рассчитанных R_i(z) и T_i(z):

equation 31
(31)
Шаг 6. Восстанавливаем полином C(z) согласно 8:
equation 32
(32)
где [c_0, \,\,\, c_1] — результат двухточечной циклической свертки последовательностей [a_0, \,\,\,a_1] и [b_0, \,\,\, b_1].

Заметим, что коэффициенты \frac{1}{2} можно учесть при вычислении двух произведений m_1 и m_2 на шаге 2.

Тогда окончательно двухточечную циклическую свертку

equation 33
(33)
можно получить через алгоритм Винограда вида:

equation 34
(34)
Таким образом, получили два умножения, 6 сложений и масштабирование с коэффициентом \frac{1}{2}.

Скрипт GNU Octave cconv2.m , реализующий расчет двухточечной циклической свертки приведен в следующем листинге:


clear all; close all; clc;

pkg load signal; % Для запуска в Matlab эту строку надо удалить!

% ------------------------------------------------------------------------------
% Исходные данные - a и b случайные комплексные векторы размера [3 x 1]
% ------------------------------------------------------------------------------
a = randn(2, 1) +  1i* randn(2, 1);
b = randn(2, 1) +  1i* randn(2, 1);


% ------------------------------------------------------------------------------
% Алгоритм Винограда расчета циклической свертки N = 4
% ------------------------------------------------------------------------------
  
s1 = b(2) + b(1);   s2 = b(2) - b(1); 

m1 = (a(2) + a(1)) * s1 / 2;
m2 = (a(2) - a(1)) * s2 / 2;

c1 = m1 - m2;
c0 = m1 + m2;


% ------------------------------------------------------------------------------
% ПРОВЕРКА РЕЗУЛЬТАТА
% ------------------------------------------------------------------------------

% Встроенная функция циклической свертки
z0 = cconv(a, b, 2)

% Циклическая свертка через FFT
z1 = ifft(fft(a).*fft(b))  

% Циклическая свертка через матричное умножение
ind = [1 2;
       2 1;];
z2 = b(ind)*a

% Циклическая свертка через алгоритм Винограда
z3 = [c0; c1]


Алгоритм для двухточечной циклической свертки можно получить напрямую из выражения (33) без использования долгих полиномиальных преобразований, раскрыв матричное произведение и вынеся общие множители за скобки. Однако для более длинных сверток простым вынесением за скобки уже не обойтись, и необходимо использовать полиномиальные преобразования алгоритма Винограда.

Циклическая свертка N=3

Рассмотрим теперь алгоритм Винограда для N=3. Тогда

equation 35
(35)
представляется произведением k=2 неприводимых многочленов Q_1(z) = z-1 и Q_2(z) = z^2+z +1. Входными данными для расчета циклической свертки являются последовательности [a_0, \,\,\,a_1, \,\,\, a_2] и [b_0, \,\,\, b_1, \,\,\, b_2]. Тогда полиномы A(z) и B(z) равны

equation 36
(36)
В соответствии с (3) и (8) можно записать:

equation 37
(37)
Шаг 1. Рассчитываем полиномы A_i(z) и B_i(z), использую теорему Безу (18), процедуру деления полиномов уголком и выражение (10):

equation 38
(38)
equation 39
(39)
Аналогично используя (11):

equation 40
(40)
equation 41
(41)
Шаг 2. Формируем полиномы C_i(z) согласно (9):

equation 42
(42)
equation 43
(43)
Рассмотрим произведение:
equation 44
(44)
Нетрудно показать, что

equation 45
(45)
где

equation 46
(46)
Тогда после деления уголком получим
equation 47
(47)
Расчет полиномов C_1(z) и C_2(z) потребовал четыре операции умножения.

Шаг 3. Рассчитываем полиномы T_i(z) используя (12):

equation 48
(48)
equation 49
(49)
Шаг 4. Расчет полиномов R_i(z) из уравнения (13).

equation 50
(50)
Полином R_2(z) найдем из уравнения

equation 51
(51)
equation 52
(52)
Чтобы выполнить это условие, полином R_2(z) должен иметь степень не ниже первой, т.е. можно записать R_2(z) = \alpha z + \beta, где \alpha и \beta можно найти из уравнения (52). Учтем, что

equation 53
(53)
тогда

equation 54
(54)
Тогда можно составить систему уравнений:

equation 55
(55)
решением которой будет \alpha = -\frac{1}{3} и \beta = -\frac{2}{3}.

Окончательно R_2(z) = -\frac{1}{3}(z + 2).

Шаг 5. Формируем полиномы S_i(z) согласно (12) на основе рассчитанных R_i(z) и T_i(z):

equation 56
(56)
equation 57
(57)

Шаг 6. Восстанавливаем полином C(z) согласно 8:

equation 58
(58)
Рассчитываем

equation 59
(59)
equation 60
(60)
Тогда
equation 61
(61)
Коэффициент \frac{1}{3} учтем в умножениях m_1 \ldots m_4, выполняемых на шаге 2. Тогда после деления уголком получим окончательно:

equation 62
(62)
где [c_0, \,\,\, c_1, \,\,\, c_2] — результат трехточечной циклической свертки последовательностей [a_0, \,\,\,a_1, \,\,\,a_2] и [b_0, \,\,\, b_1, \,\,\, b_2].

Заметим, что коэффициент \frac{1}{3} также можно учесть при вычислении умножений m_1\ldots m_4. Тогда окончательно трехточечную циклическую свертку

equation 63
(63)
можно получить через алгоритм Винограда вида:

equation 64
(64)
Полученный алгоритм требует 19 сложений и 4 умножения для расчета трехточечной циклической свертки.

Скрипт GNU Octave cconv3.m , реализующий расчет трехточечной циклической свертки приведен в следующем листинге:


clear all; close all; clc;

pkg load signal; % Для запуска в Matlab эту строку надо удалить!

% ------------------------------------------------------------------------------
% Исходные данные - a и b случайные комплексные векторы размера [3 x 1]
% ------------------------------------------------------------------------------
a = randn(3, 1) +  1i* randn(3, 1);
b = randn(3, 1) +  1i* randn(3, 1);


% ------------------------------------------------------------------------------
% Алгоритм Винограда расчета циклической свертки N = 4
% ------------------------------------------------------------------------------     
s1 = a(1) + a(2);   s2 = b(1) + b(2); 

m1 = (  s1 + a(3)) * (  s2 + b(3)) / 3;
m2 = (a(2) - a(3)) * (b(2) - b(3)) / 3;
m3 = (a(1) - a(2)) * (b(2) - b(1)) / 3;
m4 = (a(1) - a(3)) * (b(1) - b(3)) / 3;

q0 = m4 - m2;      q1 = m4 + m3;
s3 = m1 - q0;      s4 = m1 + q0;
s5 = m1 + q1;      s6 = q1 - q0;

c2 = s3 - q1; 
c1 = s5 + s6;
c0 = s4 - s6;


% ------------------------------------------------------------------------------
% ПРОВЕРКА РЕЗУЛЬТАТА
% ------------------------------------------------------------------------------

% Встроенная функция циклической свертки
z0 = cconv(a, b, 3)

% Циклическая свертка через FFT
z1 = ifft(fft(a).*fft(b))  

% Циклическая свертка через матричное умножение
ind = [1 3 2;
       2 1 3;
       3 2 1];
z2 = b(ind)*a

% Циклическая свертка через алгоритм Винограда
z3 = [c0; c1; c2]


Заметим, что для многих практических приложений вектор \mathbf{b} считается постоянным, поэтому операции сложения, использующие только b_k, где k = 0, \,\,1, \,\, 2 могут быть предварительно просчитаны и могут не учитываться в оценке вычислительной сложности алгоритма. Таким образом алгоритм может быть переписан в следующем виде:

equation 65
(65)
Операции на пересчет коэффициентов r_1 \ldots r_4 можно не учитывать. Тогда полученный алгоритм требует 4 умножения и 14 сложений.

Количество сложений при этом можно оптимизировать, если использовать следующий эквивалентный алгоритм [4, стр. 66]:

equation 66
(66)
который также требует 4 умножения, но всего 11 сложений (операции на расчет r_1 \ldots r_4 не учитываются).

Циклическая свертка N=4

\label{subsection:winograd_cconv4} Ну и наконец рассмотрим еще один пример алгоритма Винограда для расчета циклической свертки размера N=4.

Циклическую свертку последовательностей [a_0, \,\, a_1, \,\, a_2,\,\,a_3] и [b_0, \,\, b_1, \,\, b_2,\,\,b_3] можно представить в матричной форме как можно записать как:

equation 67
(67)
Тогда в модулярной арифметике можно представить:

equation 68
(68)
equation 69
(69)
Полином z^4 - 1 можно разложить на произведение трех полиномов:

equation 70
(70)
Тогда k=3 согласно теореме Винограда, 4-точечная циклическая свертка может быть реализована с использованием 2N - k = 5 умножителей.

Получим алгоритм Винограда для 4-точечной циклической свертки.

Шаг 1. Рассчитываем полиномы A_i(z) и B_i(z), использую теорему Безу (18), процедуру деления полиномов уголком и выражение (10):

equation 71
(71)
где s_1 = a_0 + a_2, s_2 = a_1 + a_3. Тогда
equation 72
(72)
Третий полином можно получить в результате деления уголком:
equation 73
(73)
где s_3 = a_0 - a_2, s_4 = a_1 - a_3.

Аналогично используя (11):

equation 74
(74)
где q_1 = b_0 + b_2, q_2 = b_1 + b_3. Тогда
equation 75
(75)
Третий полином можно получить в результате деления уголком:
equation 76
(76)
где q_3 = b_0 - b_2, q_4 = b_1 - b_3.

Шаг 2. Формируем полиномы C_i(z) согласно (9):

Первые два полинома вырождаются в умножения m_1 и m_2:

equation 77
(77)
equation 78
(78)
Третий полином C_3(z) равен
equation 79
(79)
В результате деления уголком получим остаток:

equation 80
(80)
Заметим, что коэффициенты C_3(z) можно получить из матричного уравнения:

equation 81
(81)
в котором произведение матрицы на вектор описывает двухточечную циклическую свертку.

Применим алгоритм Винограда для расчета двухточечной циклической свертки.

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

equation 82
(82)
equation 83
(83)
Согласно алгоритму Винограда для двухточечной свертки и выражению (81)

equation 84
(84)
Таким образом, расчет полиномов C_1(z), C_2(z) и C_3(z) требует 5 операций умножения.

Шаг 3. Рассчитываем полиномы T_i(z) используя (12):

equation 85
(85)
equation 86
(86)
equation 87
(87)
Шаг 4. Расчет полиномов R_i(z) из уравнения (13).

equation 88
(88)
откуда R_1(z) = \frac{1}{4}. Аналогично

equation 89
(89)
откуда R_2(z) = -\frac{1}{4}, а также

equation 90
(90)

Шаг 5. Формируем полиномы S_i(z) согласно (12) на основе рассчитанных R_i(z) и T_i(z):

equation 91
(91)
equation 92
(92)
equation 93
(93)

Шаг 6. Восстанавливаем полином C(z) согласно 8:

equation 94
(94)
Введем обозначения
equation 95
(95)
а также учтем масштабирующие коэффициенты \frac{1}{4} при расчете множителей m_1 и m_2, а коэффициенты \frac{1}{2} при p_1 и p_0 учтем в множителях m_3, m_4 и m_5.

Тогда окончательно алгоритм Винограда циклической свертки размера N = 4 имеет вид:

equation 96
(96)
Алгоритм требует 25 операций сложения и 5 умножений.

Если вектор \mathbf{b} считать постоянным, то коэффициенты q_1 \ldots q_8 могут быть предварительно рассчитанными и тогда алгоритм требует 5 умножений и 17 сложений.

Оптимизация алгоритма для минимизации сложений при постоянных значениях b_0 \ldots b_3 позволяет рассчитывать 4-точечную свертку при использовании 15 умножений как [4, стр. 67]:

equation 97
(97)
Здесь коэффициенты r_1 \ldots r_7 рассчитываются однократно, и операции для расчета данных коэффициентов можно не учитывать.

При построении малоточечных ДПФ Винограда целесообразнее использовать оптимизированные алгоритмы циклических сверток с минимальным числом умножителей и сумматоров. В этом случае коэффициенты r_1 \ldots r_7 будут постоянными поворотными коэффициентами, которые будут входить в выражение как предварительно просчитанные константы.

Скрипт GNU Octave cconv4.m , реализующий расчет четырехточечной циклической свертки приведен в следующем листинге:


clear all; close all; clc;

pkg load signal; % Для запуска в Matlab эту строку надо удалить!

% ------------------------------------------------------------------------------
% Исходные данные - a и b случайные комплексные векторы размера [4 x 1]
% ------------------------------------------------------------------------------
a = randn(4, 1) +  1i* randn(4, 1);
b = randn(4, 1) +  1i* randn(4, 1);


% ------------------------------------------------------------------------------
% Алгоритм Винограда расчета циклической свертки N = 4
% ------------------------------------------------------------------------------
% Расчет коэффициентов С1(z) = m1 и C2(z) = m2      
s1 = a(1) + a(3);   s2 = a(2) + a(4);   q1 = b(1) + b(3);   q2 = b(2) + b(4); 

s3 = a(1) - a(3);   s4 = a(2) - a(4);   q3 = b(1) - b(3);   q4 = b(2) - b(4); 

s5 = s1 + s2;       s6 = s1 - s2;       q5 = q1 + q2;       q6 = q1 - q2;

m1 = s5 * q5 / 4;   m2 = s6 * q6 / 4;


% Расчет коэффициентов С3(z) = p1 * z + p0 через двухточ. циклическую свертку  
s7 = s4 + s3;       s8 = s4 - s3;       q7 = q4 + q3;       q8 = q4 - q3;

m3 = s7 * q7 / 4;   m4 = s8 * q8 / 4;   m5 = s4*q4;
 
q9 = m3 + m4;       p0 = q9 - m5;       p1 = m3 - m4;

% восстановление C(z) =  (C1(z)*S1(z) + C2(z)*S2(z) + C3(z)*S3(z)) mod (z^4-1).
% C(z) = c0 + c1*z + c2*z^2 + c3*z^4.
% Коэффициенты 1/4 учтены в m1 и m2, а 1/2 в m3, m4 и m5.
s9  = m1 - m2;      s10 = m1 + m2; 

c0 = s10 + p0;      c1 = s9  + p1;      c2 = s10 - p0;      c3 = s9  - p1;


% ------------------------------------------------------------------------------
% ПРОВЕРКА РЕЗУЛЬТАТА
% ------------------------------------------------------------------------------

% Встроенная функция циклической свертки
z0 = cconv(a, b, 4)

% Циклическая свертка через FFT
z1 = ifft(fft(a).*fft(b))  

% Циклическая свертка через матричное умножение
ind = [1 4 3 2;
       2 1 4 3;
       3 2 1 4;
       4 3 2 1];
z2 = b(ind)*a

% Циклическая свертка через алгоритм Винограда
z3 = [c0; c1; c2; c3]


странице обсуждения статьи

Выводы

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

Мы рассмотрели методику разработки алгоритмов для произвольной длины свертки, а также рассмотрели примеры вывода алгоритмов Винограда для сверток длины N = 2, N = 3 и N = 4.

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

Список литературы
[1] Winograd S. On Computing the Discrete Fourier Transform. MATHEMATICS OF COMPUTATION, 1978, Vol. 32, Num. 141, pp. 175–189

[2] Winograd S. Some bilinear forms whose multiplicative complexity depends on the field of constants. Mathematical systems theory, 1976, Vol. 10, pp. 169–180. doi 10.1007/BF01683270

[3] Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. Москва, Наука, 1970, 720 с.

[4] Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток. Москва, Радио и связь, 1985.

[5] Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. Москва, Мир, 1989, 448 c.

Последнее изменение страницы: 12.05.2022 (19:41:21)
Страница создана Latex to HTML translator ver. 5.20.11.14