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

Содержание

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

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

Страница библиотеки на GitHub

Обнаружили ошибку? Выделите ее мышью и нажмите 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]


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

equation 88
(88)

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

equation 89
(89)
equation 90
(90)
equation 91
(91)

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

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

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

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

Скрипт 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] Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. Москва, Мир, 1989, 448 c.

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

Последнее изменение страницы: 14.11.2020 (11:14:06)
Страница создана Latex to HTML translator ver. 5.20.11.14