Реализация БИХ-фильтров с помощью набора инструкций Intel® AVX для комплексных чисел с плавающей точкой

Введение

В данной статье описывается реализация сложных фильтров с бесконечной импульсной характеристикой (БИХ, в английской литературе - IIR Infinite Impulse Response) с помощью набора SIMD инструкций Intel® AVX (Single Instruction Multiple Data – технология одновременной обработки однотипных данных одной инструкцией). Основная трудность при реализации таких алгоритмов заключается в обратной связи, присущей БИХ-фильтрам. Каждое последующее выходное значение фильтра зависит от предыдущего выходного значения, т.е. несколько выходных значений не могут быть посчитаны одновременно (параллельно). В настоящем документе изложен общий подход к векторизации БИХ и специфика ее реализации на основе архитектуры набора команд Intel AVX. Функциональные аспекты такой реализации рассматриваются ниже.

Реализация БИХ-фильтров с помощью набора инструкций Intel® AVX для комплексных чисел с плавающей точкой

Рассмотрим широко известную формулу БИХ-фильтра:

Уравнение 1

Рисунок 1. Структура БИХ-фильтра произвольного порядка.

Где “a” и “b”являются коэффициентами БИХ, “x” – вектор входных данных, “y” – выходной вектор, “order” – порядок фильтра. Из иллюстрации видна непосредственная зависимость  значениея функции от предыдущих выходных значений. Цель последующих преобразований - устранить обратную связь, т.е. зависимость выходных данных от предыдущих значений. Мы будем оперировать числами в формате с плавающей точкой и разрядностью, соответствующей инструкциям Intel AVX. Таким образом, мы заинтересованы в векторизации цикла с уменьшением количества итераций в 8 раз (или 4 при обработке комплексных чисел). В силу универсальности используемого подхода, данный метод также может быть применен к любому другому типу данных, например, к вещественным (real) или данным с двойной точностью (double) с учетом фактора векториции (8 или 4 соответственно). Последующее изложение относится к вычислениям в формате комплексных чисел с плавающей запятой, поэтому все операции по сложению, вычитанию и умножению в тексте ниже представляют собой арифметические операции с комплексными числами. Вследствие этого, наша цель – пересчет коэффициентов “a" для векторизации вычисления “y" (для устранения задержки обратной связи). Уравнение 1) может быть представлено следующим образом:

Уравнение 2

Где Fx(n) – функция, которая зависит только от “x” и представляет собой простой фильтр с конечной импульсной характеристикой (КИХ, или FIR - Finite Impulse Response). Выполним двухступенчатую последовательную фильтрацию: сначала по коэффициенту "x" (КИХ), а затем по коэффициенту "y". Часть КИХ для каждой из четырех комплексных точек может быть выражена следующими формулами (предположим для простоты, что “order”=k+1):     

Уравнение 3

КИХ-часть БИХ фильтра имеет относительно простую реализацию с использованием набора инструкций Intel AVX и векторизацией на 8 точек одинарной (float) точности (4 комплексных числа) по длине вектора. Восемь копий реальной части “b.re” и восемь копий мнимой части “b.im”, выровненные на границу 32 байта, сохраняются в структуре состояния БИХ, что является общим подходом для алгоритма Intel AVX (см. функция инициализации).

Рисунок 2. Конфигурация коэффициентов КИХ-части

b0.re, b0.re, b0.re, b0.re, b0.re, b0.re, b0.re, b0.re

-b0.im, b0.im,-b0.im, b0.im,-b0.im, b0.im,-b0.im, b0.im

b1.re, b1.re, b1.re, b1.re, b1.re, b1.re, b1.re, b1.re

-b1.im, b1.im,-b1.im, b1.im,-b1.im, b1.im,-b1.im, b1.im

…………………..

bk.re, bk.re, bk.re, bk.re, bk.re, bk.re, bk.re, bk.re

-bk.im, bk.im,-bk.im, bk.im,-bk.im, bk.im,-bk.im, bk.im

Ниже изложена схема КИХ-части БИХ-алгоритма, реализованная на базе расширенных 256-битных регистров (полный код ассемблера приводится в разделе Функция ассемблера для Коэффициентов B):

Рисунок 3. КИХ-часть алгоритма

Load VMMreg = {x3.im, x3.re, x2.im, x2.re, x1.im, x1.re, x0.im, x0.re}

vmulps {bk.re, bk.re, bk.re, bk.re, bk.re, bk.re, bk.re, bk.re} +

vshufps VMMreg = {x3.re, x3.im, x2.re, x2.im, x1.re, x1.im, x0.re, x0.im}

vmulps {-bk.im, bk.im,-bk.im, bk.im,-bk.im, bk.im,-bk.im, bk.im} +

…………………………………………………………………………………………… …………………………………………………………………………………………………… +

Load VMMreg = {x3+k.im, x3+k.re, x2+k.im, x2+k.re, x1+k.im, x1+k.re, xk.im, xk.re}

vmulps {b0.re, b0.re, b0.re, b0.re, b0.re, b0.re, b0.re, b0.re}+

Vshufps VMMreg = {x3+k.re, x3+k.im, x2+k.re, x2+k.im, x1+k.re, x1+k.im, xk.re, xk.im}

vmulps {-b0.im, b0.im,-b0.im, b0.im,-b0.im, b0.im,-b0.im, b0.im} =

Fx(n+0).re,im, Fx(n+1). re,im, Fx(n+2). re,im, Fx(n+3). re,im

Temporally store to destination vector.

После первого шага все значения Fx(n) сохраняем в результирующем векторе, и получаем:

Уравнение 4

Чтобы устранить задержку обратной связи, подставим  yn из первой строки в три последующие строки, затем yn+1 из второй строки в две последующие строки, а yn+2 из третьей строки в четвертую строку. Например, выражение для точки yn+2 будет выглядеть следующим образом:

Уравнение 5

В левой части первых четырех строк мы получим новые выходные значения - yn , yn+1 , yn+2 , yn+3, а в правой – функции, заданные от «старых» координат (yn-1 ...yn-k) и от Fx(n+0)...Fx(n+3) с новыми коэффициентами, которые рассчитываются при инициализации структуры pState (см. функция инициализации): 

Уравнение 6

Уравнение 7

Уравнение 8

Уравнение 9

Для комплексных значений одинарной точности и с учетом ширины AVX регистров для векторизации следует использовать фактор 4. Таким образом, имеем 4 новых набора комплексных коэффициентов “a” (для “y”), связанных между собой следующим образом:

  • Первый набор a0n определяется Уравнением 6;
  • Второй набор a1n определяется Уравнением 7;
  • Третий набор a2n определяется Уравнением 8;
  • Четвертый набор a3n определяется Уравнением 9;

Также необходимо выполнить ряд дополнительных вычислений для функций Fx(n+0).re,im, Fx(n+1).re,im, Fx(n+2).re,im, Fx(n+3).re,im – при этом коэффициенты для каждого набора таковы:

Рисунок 4. Дополнительный треугольник коэффициентов для БИХ-части

1.

0,

0,

0,

1;

2.

a01,

0,

0,

1;

3.

a11,

a01,

0,

1;

4.

a21,

a11,

a01,

1;

Аналогичным образом, как и на первой стадии КИХ, все новые коэффициенты рассчитываются на стадии инициализации (см. функцию инициализации) и сохраняются в естественном для Intel AVX порядке с надлежащим выравниванием (часть, состоящая из одних единиц, не требуется и опускается):

Рисунок 5. Конфигурация коэффициентов БИХ-части

a01.re,im,  a11.re,im, a21.re,im, a31.re,im

a02.re,im,  a12.re,im, a22.re,im, a32.re,im

…....................

a0k.re,im, a1k.re,im, a2k.re,im, a3k.re,im

  0.re,im, a01.re,im, a11.re,im, a21.re,im

  0.re,im,   0.re,im, a01.re,im, a11.re,im

  0.re,im,   0.re,im,    0.re,im, a01.re,im

Таким образом, на каждые четыре выходных значения мы затрачиваем дополнительно 10 операций сложения и 6 операций умножения с комплексными числами («треугольник» дополнитеных  коэффициентов (рис 4) умножается на вектор {Fx(n+0).re,im, Fx(n+1).re,im, Fx(n+2).re,im, Fx(n+3).re,im} . Таковы вычислительные затраты алгоритма, позволяющего устранить обратную. связь:

Рисунок 6. БИХ-часть алгоритма

vbroadcast to VMM = {yn-1.re, yn-1.re, yn-1.re, yn-1.re,yn-1.re, yn-1.re, yn-1.re, yn-1.re}

vmulps {a01.re, a01.im, a11.re, a11.im, a21.re, a21.im, a31.re, a31.im} +

vbroadcast to VMM = {yn-1.im, yn-1.im, yn-1.im, yn-1.im,yn-1.im, yn-1.im, yn-1.im, yn-1.im}

vmulps {a01.im, a01.re, a11.im, a11.re, a21.im, a21.re, a31.im, a31.re} +

……………… ……………… +

vbroadcast to VMM = {yn-k.re, yn-k.re, yn-k.re, yn-k.re,yn-k.re, yn-k.re, yn-k.re, yn-k.re}

vmulps {a0k.re, a0k.im, a1k.re, a1k.im, a2k.re, a2k.im, a3k.re, a3k.im} +

vbroadcast to VMM = {yn-k.im, yn-k.im, yn-k.im, yn-k.im,yn-k.im, yn-k.im, yn-k.im, yn-k.im}

vmulps {a0k.im, a0k.re, a1k.im, a1k.re, a2k.im, a2k.re, a3k.im, a3k.re} +

vbroadcast to VMM = { Fx(n+0), Fx(n+0), Fx(n+0), Fx(n+0)}re.im

vmulps { a0k.re,im, a1k.re,im, a2k.re,im, a3k.re,im } + (extra)

vbroadcast to VMM = { Fx(n+1), Fx(n+1), Fx(n+1), Fx(n+1)}re.im

vmulps { 0.re,im, a01.re,im, a11.re,im, a21.re,im } + (extra)

vbroadcast to VMM = { Fx(n+2), Fx(n+2), Fx(n+2), Fx(n+2)}re.im

vmulps { 0.re,im, 0.re,im, a01.re,im, a11.re,im } + (extra)

vbroadcast to VMM = { Fx(n+3), Fx(n+3), Fx(n+3), Fx(n+3)}re.im

vmulps { 0.re,im, 0.re,im, 0.re,im, a01.re,im } + (extra)

vload to VMM = { Fx(n+3), Fx(n+2), Fx(n+1), Fx(n+0)}re.im = vstore VMM = {Yn+0, Yn+1, Yn+2, Yn+3}re.im

Temporally store to destination vector

Данный подход позволяет устранить обратную связь по выходным значениям фильтра и повысить производительность БИХ-фильтров для всех типов SIMD-архитектуры более чем в два раза.

Ниже приводится полный листинг вышеописанного алгоритма БИХ, используемого в библиотеках Intel® IPP (Intel® Integrated Performance Primitives).

Рисунок 7. Блок-схема алгоритма БИХ

IppStatus ippsIIRInitAlloc_32fc( IppsIIRState_32fc** ppState, const Ipp32fc* pTaps, int order, const Ipp32fc* pDlyLine);

Order of taps = B0,B1,...,Border,A0,A1,...,Aorder

  1. This function allocates memory and initializes an arbitrary IIR filter state.

Taps are normed: a1=A1/A0 , a2=A2/A0 , ..., border=Border/A0 .

Taps are recalculated according to (6) – (9) and stored by eight coefficients according to Figure 1 and Figure 4 order with appropriate alignment (on 32 byte boundary) in pState structure.

IppStatus ippsIIR_32fc( const Ipp32fc *pSrc, Ipp32fc *pDst, int len, IppsIIRState_32fc *pState )

{

The first stage: filtering “by x” as simple FIR filter does with bk coefficients (the code written in Intel AVX). Temporal output is written to pDst vector (Figure 2 and Assembler Function for B Coefficients).

The second stage: filtering “by y” by bundles of four (the code written in Intel AVX) with recalculated ak coefficients (Figure 5 and Assembler Function for A Coefficients)

}

Заключение

В силу наличия обратной связи предыдущие реализации БИХ-фильтров для архитектуры IA32 имели скалярную реализацию. Это означает, что производительность векторной функции на выходную точку практически соответствовало производительности для скалярной (точечной) функции. SSE векторизация осуществлялась по «порядку фильтра», а не по длине вектора. Заметное повышение производительности при таком подходе в сравнение с простейшей C-реализацией наблюдалось лишь для фильтров выше 32-го порядка. БИХ-фильтры столь больших порядков практически не используются при обработке сигналов. Для порядков ниже 32 традиционная реализация БИХ-фильтра посредством команд Intel SSE не связана с каким-либо увеличением производительности из-за того, что латентности инструкций суммируется из-за наличия зависимости по данным. Представленный подход к устранению обратной связи, реализованный посредством набора инструкций Intel AVX, позволяет фильтру 6-го порядка  достичь ускорения в ~1,6 раз по сравнению с использованием аналогичного кода команд Intel SSE. Необходимо принять во внимание, что данный показатель был получен с помощью программного симулятора. Фактический  результат на реальных процессорах может варьироваться в зависимости от различных обстоятельств.

Приложение A: Код функции

Функция инициализации

Код функции инициализации (скачать)

Функция фильтрации

Код функции фильтрации (скачать)  

Функция ассемблера для коэффициентов B

Код функции ассемблера для коэффициентов B (скачать)

Функция ассемблера для коэффициентов A

Код функции ассемблера для коэффициентов A (скачать)

Об авторе

Игорь Астахов – руководитель отдела разработки ПО в Группе программных продуктов и решений (подразделения ПО визуальных вычислений, CIP и библиотек Intel IPP). Его специализацией является разработка библиотек Intel IPP для всех существующих и будущих типов архитектуры Intel. Игорь работает в компании Intel в течение пяти лет.

Categories:
For more complete information about compiler optimizations, see our Optimization Notice.