Сейчас проблема производительности приложений с цифровой обработкой сигналов решается использованием специализированного оборудования. При этом такому оборудованию требуются еще и оптимизированные под него программы.
Есть и несколько подходов для подгонки алгоритмов под машинную архитектуру, один из которых – управление данными(data driven). Каждый, кто сталкивался с ручным программированием данных, знает, что это дело это не простое. Однако в большинстве случаев можно спроектировать пре-компилятор, который существенно упростит задачу. В статье описана методика построения двух управляемых данными алгоритмов БПФ и способы достижения максимальной производительности, превосходящей теоретическую.
N временных отсчетов любого сигнала могут быть преобразованы в N фазо-частот при помощи Прямого Дискретного Преобразования Фурье.
Xk = n=0N-1? xn* ?Nkn, k= 0..N-1, ?N= e-2?i/N, i = 2v-1 ;
Для того чтобы получить комплексный фазо-частотный спектр {X} нужно и временные отсчеты {x} перевести в комплексный вид с нулевой мнимой частью. Поворачивающийся множитель ? представляет собой “вращающийся вектор” на единичной окружности в комплексной плоскости.
Для упрощения вычислений используется Быстрое Преобразование Фурье, где N является степенью от 2. В таком случае комплексный сигнал X=(x0, x1 … xn-1) можно поделить на два A=(x0, x2 … xn-2) и B=(x1, x3 … xn-1) уже из N/2 отсчетов.
Ak = n=0N/2-1? x2 n* ?Nk(2n) ;
Bk= n=0N/2-1? x2 n+1 * ?Nk(2n+1) = n=0N/2-1? x2 n+1 * ?Nk* ?Nk(2n) = ?Nk* n=0N/2-1? x2 n+1* ?Nk(2n);
Теперь можно вывести нужную нам формулу: Xk = Ak + ?Nk * Bk;
Здесь N фазо-частот нашего сигнала можно получить, применяя БПФ по очереди для N/2 четных отсчетов и для N/2 нечетных отсчетов. На следующем шаге A и B опять разбиваются пополам. Стоит заметить, что этот подход применим к Обратному Преобразованию Фурье, которое здесь не рассматривается.
Рекурсивная реализация[1] БПФ по приведенным выше формулам проста и приведена ниже:
Действительные части комплексных значений входного вектора состоят из входных параметров, количество которых усекается до ближайшей степени 2-ки. Первый уровень рекурсии разбивает сигнал In из n-значений на 2 сигнала A и B, каждый из которых состоит из n/2 значений. Следующий уровень разбивает данные на 4 сигнала по n/4 значения. Рекурсия прекращается, когда остается только сигнал из одного значения. Сигналы A и B преобразуются в At и Bt, а затем в возвращаемый вектор сигналов Out. В конце программа выводит реальные и мнимые части фаза-частотного сигнала из выходного вектора.
Чем выше уровень абстракции языка приложения, тем меньше усилий требуется для программирования. Но тем ниже производительность. Можно писать более эффективный код, снижая уровень абстракции, но тратя больше усилий. Разберем это на основе элементарной операции БПФ, которая может быть представлена как C = a + w * b в комплексном виде. Эта формула в общем случае реализуется через умножения с накоплением, поэтому имеет аббревиатуру MAC(multiply–accumulate operation):
Логарифм по основанию 2 от размера сигнала определяет число рекурсий. Несколько стилей n2ln функции приведены ниже.
В последней реализации наихудший случай исполнения равен среднему случаю, который С-компилятор может реализовать при помощи 3-х сдвигов и 4-х условных переходов. Кроме того эта реализация может быть представлена в виде макроса для расчета значений констант во время компиляции.
Классическая реализация рекурсивного БПФ на С похожа на С++ вариант
Относительно С++ реализации здесь сделано несколько алгоритмических изменений. Во-первых, комплексные расчеты производятся вручную. Во-вторых, для экономии памяти входной и выходной вектор из реализации C++ объединены в один буфер. Буфера выделяются явно, потому размер преобразования нужно передавать в рекурсивную функцию.
Пример рекурсивного БПФ на основе управления данными, сделанный из предыдущей реализации на С.
Здесь буфера для вычислений содержит не сами комплексные значения, а указатели на них. В Mac массив фактически последовательно записываются отложенные элементарные БНФ операции, для того чтобы их сделать потом. Другими словами, это байт-код элементарных БНФ операций.
Двумерный массив DataTrace используется для поддержки этих операций. После вызова рекурсивной процедуры пользователь должен вызвать calculate_macs для исполнения сгенерированного байт-кода. Эта процедура имеет всего один цикл, но может применяться многократно. Это максимальная производительность для теоретической сложности БПФ n*log2(n). Но проблема тут в памяти — Mac и DataTrace так же имеют n*log2(n) элементов. И это слишком много для недорогих встроенных решений.
Теперь пришло время разделить предыдущую реализацию БПФ на две программы. Первая будет генерировать С структуры байт-кода, а вторая их исполнять. При таком подходе генератор С структур фактически является пре-компилятором, в котором можно, не жалея ресурсов, реализовывать различные оптимизационные стратегии, например, для оптимизации памяти RAM. Ранее и DataTrace массив N*log2(N) элементов, в программе ниже его аналог массив XY имеет всего N+2 элементов.
Это пример БПФ для 16 отсчетов. Один элемент tbl массива содержит комплексное значение поворачивающего множителя и три смещения для организации вычислений на ячейках массива XY. При этом сам код имеет всего один “for” цикл.
Следующий пример основан на группировке элементарных БПФ операций относительно поворачивающего множителя.
Это пример БПФ для 32 отсчетов. Количество элементарных БПФ операций — N. Массив XY организован по схеме свопа и его размер — 2*N. Это дает возможность условной инструкции XMAC оперировать с 3-мя банками памяти одновременно, хотя на практике компиляторами это игнорируется. Тем не менее, XMAC может быть теоретически реализован даже одной машинной инструкцией.
Самый известный БПФ алгоритм использует сложную логику перестановки адресов, которая в графическом виде напоминает крылья бабочки. А вот в данном примере новые адреса получаются простым инкрементом от старых, потому с названием все просто — это метод гусеницы. Стоит заметить, что длинный switch оператор так же напоминает гусеницу.
Элементарная комплексная операция БПФ (c = a + ? * b) состоит из 4-х обычных операций:
В литературе по БПФ [2] для этих вычислений можно найти несколько стратегий оптимизаций, которые классифицированы ниже:
Почти все эти оптимизации БПФ применимы к методу гусеницы.
Возможны также архитектурно-машинные оптимизации. Например, элементарная БПФ операция в коде выше реализована как макрос XMAC. Она может быть распараллелена при помощи SIMD инструкций для x86-процессоров Intel AVX:
Приведенный макрос поддерживает две операции с плавающей точкой одновременно при помощи 128-битных регистров. Но стоит взглянуть на гусеницу еще раз – из-за особенностей адресации ближайшие XMAC-и могут быть объединены вместе, например, для реализации при помощи 512-регистров(AVX3).
В завершении стоит сказать, что направлений развития пре-компилятора куда больше чем возможностей. Поэтому целью данной статьи является сбор дополнительных требований и выявления областей, где этот подход может быть полезен.
[1] Notes on Recursive FFT (Fast Fourier Transform) algorithm, Fall 2010, COSC 511, Susan Haynes
[2] Notes on the FFT, C. S. Burrus, Department of Electrical and Computer Engineering
Rice University, Houston, TX 77251-1892
[3] Caterpillar Implementation Based on Generated Code
[4] Some fft_caterpillar examples, https://github.com/r35382/fft_caterpillar
Есть и несколько подходов для подгонки алгоритмов под машинную архитектуру, один из которых – управление данными(data driven). Каждый, кто сталкивался с ручным программированием данных, знает, что это дело это не простое. Однако в большинстве случаев можно спроектировать пре-компилятор, который существенно упростит задачу. В статье описана методика построения двух управляемых данными алгоритмов БПФ и способы достижения максимальной производительности, превосходящей теоретическую.
Немного математики
N временных отсчетов любого сигнала могут быть преобразованы в N фазо-частот при помощи Прямого Дискретного Преобразования Фурье.
Xk = n=0N-1? xn* ?Nkn, k= 0..N-1, ?N= e-2?i/N, i = 2v-1 ;
Для того чтобы получить комплексный фазо-частотный спектр {X} нужно и временные отсчеты {x} перевести в комплексный вид с нулевой мнимой частью. Поворачивающийся множитель ? представляет собой “вращающийся вектор” на единичной окружности в комплексной плоскости.
Для упрощения вычислений используется Быстрое Преобразование Фурье, где N является степенью от 2. В таком случае комплексный сигнал X=(x0, x1 … xn-1) можно поделить на два A=(x0, x2 … xn-2) и B=(x1, x3 … xn-1) уже из N/2 отсчетов.
Ak = n=0N/2-1? x2 n* ?Nk(2n) ;
Bk= n=0N/2-1? x2 n+1 * ?Nk(2n+1) = n=0N/2-1? x2 n+1 * ?Nk* ?Nk(2n) = ?Nk* n=0N/2-1? x2 n+1* ?Nk(2n);
Теперь можно вывести нужную нам формулу: Xk = Ak + ?Nk * Bk;
Здесь N фазо-частот нашего сигнала можно получить, применяя БПФ по очереди для N/2 четных отсчетов и для N/2 нечетных отсчетов. На следующем шаге A и B опять разбиваются пополам. Стоит заметить, что этот подход применим к Обратному Преобразованию Фурье, которое здесь не рассматривается.
БПФ на C++
Рекурсивная реализация[1] БПФ по приведенным выше формулам проста и приведена ниже:
#include <complex>
#include <vector>
#include <algorithm>
#include <iostream>
#include <math.h>
#define M_PI 3.14159265358979323846
using namespace std;
typedef complex<double> w_type;
static vector<w_type> fft(const vector<w_type> &In)
{
int i = 0, wi = 0;
int n = In.size();
vector<w_type> A(n / 2), B(n / 2), Out(n);
if (n == 1) {
return vector<w_type>(1, In[0]);
}
i = 0;
copy_if( In.begin(), In.end(), A.begin(), [&i] (w_type e) {
return !(i++ % 2);
} );
copy_if( In.begin(), In.end(), B.begin(), [&i] (w_type e) {
return (i++ % 2);
} );
vector<w_type> At = fft(A);
vector<w_type> Bt = fft(B);
transform(At.begin(), At.end(), Bt.begin(), Out.begin(), [&wi, &n]
(w_type& a, w_type& b) {
return a + b * exp(w_type(0, 2 * M_PI * wi++ / n));
});
transform(At.begin(), At.end(), Bt.begin(), Out.begin() + n / 2, [&wi, &n]
(w_type& a, w_type& b) {
return a + b * exp(w_type(0, 2 * M_PI * wi++ / n));
});
return Out;
}
void main(int argc, char* argv[])
{
int ln = (int)floor( log(argc - 1.0) / log(2.0) );
vector<w_type> In(1 << ln);
std::transform(argv + 1, argv + argc, In.begin(),[&](const char* arg) {
return w_type(atof(arg), 0);
});
vector<w_type> Out = fft(In);
for (vector<w_type>::iterator itr = Out.begin(); itr != Out.end(); itr++) {
cout << *itr << endl;
}
}
Действительные части комплексных значений входного вектора состоят из входных параметров, количество которых усекается до ближайшей степени 2-ки. Первый уровень рекурсии разбивает сигнал In из n-значений на 2 сигнала A и B, каждый из которых состоит из n/2 значений. Следующий уровень разбивает данные на 4 сигнала по n/4 значения. Рекурсия прекращается, когда остается только сигнал из одного значения. Сигналы A и B преобразуются в At и Bt, а затем в возвращаемый вектор сигналов Out. В конце программа выводит реальные и мнимые части фаза-частотного сигнала из выходного вектора.
Уровни абстракции
Чем выше уровень абстракции языка приложения, тем меньше усилий требуется для программирования. Но тем ниже производительность. Можно писать более эффективный код, снижая уровень абстракции, но тратя больше усилий. Разберем это на основе элементарной операции БПФ, которая может быть представлена как C = a + w * b в комплексном виде. Эта формула в общем случае реализуется через умножения с накоплением, поэтому имеет аббревиатуру MAC(multiply–accumulate operation):
| Прототип | void w_mac( w_type* cc, w_type a, w_type w, w_type b) |
| С++ | *cc = a + b * exp(w_type(0, 2 * M_PI * w / n)) |
| Классический С | cc->r = a.r + w.r * b.r - w.i * b.i; |
| Оптимизированный С | register double reg; |
| Прототип | int n2ln( int n ) |
| С++ | return (int)floor( log(n) / log(2.0) ); |
| Классический С | int ln = 0; |
| Встроенный С | |
БПФ на С
Классическая реализация рекурсивного БПФ на С похожа на С++ вариант
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define M_PI 3.14159265358979323846
typedef struct { double r; double i; } w_type;
int n2ln( int n );
void w_mac( w_type* cc, w_type a, w_type w, w_type b );
static void fft0( w_type InOut[], int n )
{
int i;
w_type w, *A, *B;
if (n == 1) return;
A = malloc( sizeof(w_type) * n / 2 );
B = malloc( sizeof(w_type) * n / 2 );
for (i = 0; i < n / 2; i++) {
A[i] = InOut[i * 2];
B[i] = InOut[i * 2 + 1];
}
fft0( A, n / 2 );
fft0( B, n / 2 );
for (i = 0; i < n; i++) {
w.r = cos(2 * M_PI * i / n);
w.i = sin(2 * M_PI * i / n);
w_mac( &InOut[i], A[i % (n / 2)], w, B[i % (n / 2)] );
}
free( A );
free( B );
}
void main( int argc, char * argv[] )
{
int i;
int ln = n2ln(argc - 1);
w_type* InOut = malloc( (1 << ln) * sizeof(w_type) );
for (i = 0; i < (1 << ln); i++) {
InOut[i].r = atof(argv[i+1]);
InOut[i].i = 0;
}
fft0( InOut, 1 << ln );
for(i = 0; i < (1 << ln); i++) {
printf("%.4f %.4f\n", InOut[i].r, InOut[i].i);
}
}Относительно С++ реализации здесь сделано несколько алгоритмических изменений. Во-первых, комплексные расчеты производятся вручную. Во-вторых, для экономии памяти входной и выходной вектор из реализации C++ объединены в один буфер. Буфера выделяются явно, потому размер преобразования нужно передавать в рекурсивную функцию.
С реализация БПФ на основе управления данными
Пример рекурсивного БПФ на основе управления данными, сделанный из предыдущей реализации на С.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define M_PI 3.14159265358979323846
#define LN_FFT 4 /* number of samples is 1 << LN_FFT */
typedef struct { double r; double i; } w_type;
int n2ln( int n );
void w_mac( w_type* cc, w_type a, w_type w, w_type b );
static struct tMac {
w_type *c, *a, *b, w;
} Mac[LN_FFT * (1 << LN_FFT)];
static int nMac;
static w_type DataTrace[LN_FFT + 1][1 << LN_FFT];
static int BusyDataTrace[LN_FFT + 1];
static void calculate_macs()
{
int i;
for (i = 0; i < nMac; i++) {
w_mac(Mac[i].c, *Mac[i].a, Mac[i].w, *Mac[i].b);
}
}
static void record_mac( w_type** cc, w_type* a, w_type w, w_type *b, int n )
{
int ln = n2ln(n);
int i = BusyDataTrace[ln]++;
*cc = &DataTrace[ln][i];
Mac[nMac].c = &DataTrace[ln][i];
Mac[nMac].w = w;
Mac[nMac].a = a;
Mac[nMac].b = b;
nMac++;
}
static void fft0( w_type* InOut[], int n )
{
int i;
w_type w, **A, **B;
if (n == 1) return;
A = malloc( sizeof(w_type*) * n / 2 );
B = malloc( sizeof(w_type*) * n / 2 );
for (i = 0; i < n / 2; i++) {
A[i] = InOut[i * 2];
B[i] = InOut[i * 2 + 1];
}
fft0( &A[0], n / 2 );
fft0( &B[0], n / 2 );
for (i = 0; i < n; i++) {
w.r = cos(2 * M_PI * i / n);
w.i = sin(2 * M_PI * i / n);
record_mac( &InOut[i], A[i % (n / 2)], w, B[i % (n / 2)], n );
}
free(A);
free(B);
}
void main( int argc, char* argv[] )
{
int i;
w_type** InOut = malloc( sizeof(w_type*) * (1 << LN_FFT) );
for (i = 0; i < (1 << LN_FFT); i++) {
InOut[i] = &DataTrace[0][i];
DataTrace[0][i].r = atof( argv[i % (argc - 1) + 1] );
DataTrace[0][i].i = 0;
}
fft0( InOut, 1 << LN_FFT );
calculate_macs();
for(i = 0; i < (1 << LN_FFT); i++) {
printf("%.4f %.4f\n", DataTrace[LN_FFT][i].r, DataTrace[LN_FFT][i].i);
}
free(InOut);
} Здесь буфера для вычислений содержит не сами комплексные значения, а указатели на них. В Mac массив фактически последовательно записываются отложенные элементарные БНФ операции, для того чтобы их сделать потом. Другими словами, это байт-код элементарных БНФ операций.
Двумерный массив DataTrace используется для поддержки этих операций. После вызова рекурсивной процедуры пользователь должен вызвать calculate_macs для исполнения сгенерированного байт-кода. Эта процедура имеет всего один цикл, но может применяться многократно. Это максимальная производительность для теоретической сложности БПФ n*log2(n). Но проблема тут в памяти — Mac и DataTrace так же имеют n*log2(n) элементов. И это слишком много для недорогих встроенных решений.
Табличная реализация
Теперь пришло время разделить предыдущую реализацию БПФ на две программы. Первая будет генерировать С структуры байт-кода, а вторая их исполнять. При таком подходе генератор С структур фактически является пре-компилятором, в котором можно, не жалея ресурсов, реализовывать различные оптимизационные стратегии, например, для оптимизации памяти RAM. Ранее и DataTrace массив N*log2(N) элементов, в программе ниже его аналог массив XY имеет всего N+2 элементов.
#include <stdlib.h>
#include <stdio.h>
#define LN_FFT 4 /* power of 2 is number of fft points */
/* */
#define W_0_02 1.000000000000000 /* angle 0.00 dg */
#define W_1_04 0.000000000000000 /* angle 90.00 dg */
#define W_1_08 0.707106781186548 /* angle 45.00 dg */
#define W_1_16 0.923879532511287 /* angle 22.50 dg */
#define W_3_16 0.382683432365090 /* angle 67.50 dg */
typedef struct { double r; double i; } w_type;
static const struct fft_tbl {
double wr, wi;
unsigned char c, a, b;
} tbl[] = {
{ W_0_02,+W_1_04,16, 0, 8}, {-W_0_02,+W_1_04,17, 0, 8},
{ W_0_02,+W_1_04, 0, 4,12}, {-W_0_02,+W_1_04, 8, 4,12},
{ W_0_02,+W_1_04, 4, 2,10}, {-W_0_02,+W_1_04,12, 2,10},
{ W_0_02,+W_1_04, 2, 6,14}, {-W_0_02,+W_1_04,10, 6,14},
{ W_0_02,+W_1_04, 6, 1, 9}, {-W_0_02,+W_1_04,14, 1, 9},
{ W_0_02,+W_1_04, 1, 5,13}, {-W_0_02,+W_1_04, 9, 5,13},
{ W_0_02,+W_1_04, 5, 3,11}, {-W_0_02,+W_1_04,13, 3,11},
{ W_0_02,+W_1_04, 3, 7,15}, {-W_0_02,+W_1_04,11, 7,15},
{ W_0_02,+W_1_04, 7,16, 0}, {-W_0_02,+W_1_04,15,16, 0},
{ W_1_04,+W_0_02,16,17, 8}, {-W_1_04,-W_0_02, 0,17, 8},
{ W_0_02,+W_1_04,17, 4, 2}, {-W_0_02,+W_1_04, 8, 4, 2},
{ W_1_04,+W_0_02, 4,12,10}, {-W_1_04,-W_0_02, 2,12,10},
{ W_0_02,+W_1_04,12, 6, 1}, {-W_0_02,+W_1_04,10, 6, 1},
{ W_1_04,+W_0_02, 6,14, 9}, {-W_1_04,-W_0_02, 1,14, 9},
{ W_0_02,+W_1_04,14, 5, 3}, {-W_0_02,+W_1_04, 9, 5, 3},
{ W_1_04,+W_0_02, 5,13,11}, {-W_1_04,-W_0_02, 3,13,11},
{ W_0_02,+W_1_04,13, 7,17}, {-W_0_02,+W_1_04,11, 7,17},
{ W_1_08,+W_1_08, 7,16, 4}, {-W_1_08,-W_1_08,17,16, 4},
{ W_1_04,+W_0_02,16,15, 8}, {-W_1_04,-W_0_02, 4,15, 8},
{-W_1_08,+W_1_08,15, 0, 2}, { W_1_08,-W_1_08, 8, 0, 2},
{ W_0_02,+W_1_04, 0,12,14}, {-W_0_02,+W_1_04, 2,12,14},
{ W_1_08,+W_1_08,12, 6, 5}, {-W_1_08,-W_1_08,14, 6, 5},
{ W_1_04,+W_0_02, 6,10, 9}, {-W_1_04,-W_0_02, 5,10, 9},
{-W_1_08,+W_1_08,10, 1, 3}, { W_1_08,-W_1_08, 9, 1, 3},
{ W_0_02,+W_1_04, 1,13, 0}, {-W_0_02,+W_1_04, 3,13, 0},
{ W_1_16,+W_3_16,13, 7,12}, {-W_1_16,-W_3_16, 0, 7,12},
{ W_1_08,+W_1_08, 7,16, 6}, {-W_1_08,-W_1_08,12,16, 6},
{ W_3_16,+W_1_16,16,15,10}, {-W_3_16,-W_1_16, 6,15,10},
{ W_1_04,+W_0_02,15,11, 2}, {-W_1_04,-W_0_02,10,11, 2},
{-W_3_16,+W_1_16,11,17,14}, { W_3_16,-W_1_16, 2,17,14},
{-W_1_08,+W_1_08,17, 4, 5}, { W_1_08,-W_1_08,14, 4, 5},
{-W_1_16,+W_3_16, 4, 8, 9}, { W_1_16,-W_3_16, 5, 8, 9},
};
static const unsigned char OutOrder[]={
1,13,7,16,15,11,17,4,3,0,12,6,10,2,14,5,};
static struct { double r; double i; } XY[(1 << LN_FFT) + 2];
void fft_table()
{
int i;
register const struct fft_tbl* t;
for (i = 0, t = tbl; i < sizeof(tbl) / sizeof(tbl[0]); i++, t++) {
XY[t->c].r = XY[t->a].r + t->wr * XY[t->b].r - t->wi * XY[t->b].i;
XY[t->c].i = XY[t->a].i + t->wr * XY[t->b].i + t->wi * XY[t->b].r;
}
}
void main(int argc, char* argv[])
{
int i;
for (i = 0; i < (1 << LN_FFT); i++) {
XY[i].r = atof( argv[i % (argc - 1) + 1] );
XY[i].i = 0;
}
fft_table();
for(i = 0; i < (1 << LN_FFT); i++) {
printf("%.4f %.4f\n", XY[OutOrder[i]].r, XY[OutOrder[i]].i);
}
}Это пример БПФ для 16 отсчетов. Один элемент tbl массива содержит комплексное значение поворачивающего множителя и три смещения для организации вычислений на ячейках массива XY. При этом сам код имеет всего один “for” цикл.
Метод Гусеницы
Следующий пример основан на группировке элементарных БПФ операций относительно поворачивающего множителя.
#include <stdlib.h>
#include <stdio.h>
#define LN_FFT 5 /* power of 2 is number of fft points */
/* */
#define W_0_02 1.000000000000000 /* angle 0.00 dg */
#define W_0_04 0.000000000000000 /* angle 90.00 dg */
#define W_0_08 0.707106781186547 /* angle 45.00 dg */
#define W_0_16 0.923879532511287 /* angle 22.50 dg */
#define W_1_16 0.382683432365090 /* angle 67.50 dg */
#define W_0_32 0.980785280403230 /* angle 11.25 dg */
#define W_1_32 0.831469612302545 /* angle 33.75 dg */
#define W_2_32 0.555570233019602 /* angle 56.25 dg */
#define W_3_32 0.195090322016128 /* angle 78.75 dg */
typedef struct { double r; double i; } w_type;
#define X2Y(a) (a + (1 << (LN_FFT - 1)))
#define XMAC(c, a, wr, wi) c->r = a->r + wr * X2Y(a)->r - wi * X2Y(a)->i; c->i = a->i + wr * X2Y(a)->i + wi * X2Y(a)->r;
static w_type XY[2][(1 << LN_FFT)];
static const w_type* pOut = LN_FFT % 2 ? &XY[1][0] : &XY[0][0];
static const unsigned char OutOrder[]={31,15,23,14,27,13,22,12,29,11,21,
10,26,9,20,8,30,7,19,6,25,5,18,4,28,3,17,2,24,1,16,0,};
void fft_caterpillar()
{
int i, j, lim;
register w_type *pc, *pa; /* pb = a + (1 << (LN_FFT - 1)) */
for (i = 1; i <= LN_FFT; i++) {
pc = i % 2 ? &XY[1][0] : &XY[0][0];
pa = i % 2 ? &XY[0][0] : &XY[1][0];
lim = 1 << (LN_FFT - i);
for (j = 0; j < lim; j++) {
switch (i) {
case 5:
XMAC(pc, pa, W_0_32, -W_3_32); pc++; pa += 1;
XMAC(pc, pa, W_1_32, -W_2_32); pc++; pa += 1;
XMAC(pc, pa, W_2_32, -W_1_32); pc++; pa += 1;
XMAC(pc, pa, W_3_32, -W_0_32); pc++; pa += 1;
XMAC(pc, pa, -W_3_32, -W_0_32); pc++; pa += 1;
XMAC(pc, pa, -W_2_32, -W_1_32); pc++; pa += 1;
XMAC(pc, pa, -W_1_32, -W_2_32); pc++; pa += 1;
XMAC(pc, pa, -W_0_32, -W_3_32); pc++; pa += 1;
pa -= 8;
XMAC(pc, pa, -W_0_32, +W_3_32); pc++; pa += 1;
XMAC(pc, pa, -W_1_32, +W_2_32); pc++; pa += 1;
XMAC(pc, pa, -W_2_32, +W_1_32); pc++; pa += 1;
XMAC(pc, pa, -W_3_32, +W_0_32); pc++; pa += 1;
XMAC(pc, pa, W_3_32, +W_0_32); pc++; pa += 1;
XMAC(pc, pa, W_2_32, +W_1_32); pc++; pa += 1;
XMAC(pc, pa, W_1_32, +W_2_32); pc++; pa += 1;
XMAC(pc, pa, W_0_32, +W_3_32); pc++; pa += 1;
case 4:
XMAC(pc, pa, W_0_16, -W_1_16); pc++; pa += 1;
XMAC(pc, pa, W_1_16, -W_0_16); pc++; pa += 1;
XMAC(pc, pa, -W_1_16, -W_0_16); pc++; pa += 1;
XMAC(pc, pa, -W_0_16, -W_1_16); pc++; pa += 1;
pa -= 4;
XMAC(pc, pa, -W_0_16, +W_1_16); pc++; pa += 1;
XMAC(pc, pa, -W_1_16, +W_0_16); pc++; pa += 1;
XMAC(pc, pa, W_1_16, +W_0_16); pc++; pa += 1;
XMAC(pc, pa, W_0_16, +W_1_16); pc++; pa += 1;
case 3:
XMAC(pc, pa, W_0_08, -W_0_08); pc++; pa += 1;
XMAC(pc, pa, -W_0_08, -W_0_08); pc++; pa += 1;
pa -= 2;
XMAC(pc, pa, -W_0_08, +W_0_08); pc++; pa += 1;
XMAC(pc, pa, W_0_08, +W_0_08); pc++; pa += 1;
case 2:
XMAC(pc, pa, -W_0_04, -W_0_02); pc++; pa += 1;
pa -= 1;
XMAC(pc, pa, W_0_04, +W_0_02); pc++; pa += 1;
case 1:
XMAC(pc, pa, -W_0_02, +W_0_04); pc++; pa += 1;
pa -= 1;
case 0:
XMAC(pc, pa, W_0_02, +W_0_04); pc++; pa += 1;
}
}
}
}
void main(int argc, char* argv[])
{
int i;
for (i = 0; i < (1 << LN_FFT); i++) {
XY[0][i].r = atof( argv[i % (argc - 1) + 1] );
XY[0][i].i = 0;
}
fft_caterpillar();
for(i = 0; i < (1 << LN_FFT); i++) {
printf("%.4f %.4f\n", pOut[OutOrder[i]].r, pOut[OutOrder[i]].i);
}
}Это пример БПФ для 32 отсчетов. Количество элементарных БПФ операций — N. Массив XY организован по схеме свопа и его размер — 2*N. Это дает возможность условной инструкции XMAC оперировать с 3-мя банками памяти одновременно, хотя на практике компиляторами это игнорируется. Тем не менее, XMAC может быть теоретически реализован даже одной машинной инструкцией.
Самый известный БПФ алгоритм использует сложную логику перестановки адресов, которая в графическом виде напоминает крылья бабочки. А вот в данном примере новые адреса получаются простым инкрементом от старых, потому с названием все просто — это метод гусеницы. Стоит заметить, что длинный switch оператор так же напоминает гусеницу.
Дальнейшие оптимизации
Элементарная комплексная операция БПФ (c = a + ? * b) состоит из 4-х обычных операций:
reg = a.real + w.real * b.real;
c.real = reg - w.imag * b.imag;
reg = a.imag + w.real * b.imag;
c.imag = reg + w.imag * b.real;В литературе по БПФ [2] для этих вычислений можно найти несколько стратегий оптимизаций, которые классифицированы ниже:
- Оптимизации, основанные на особенностях входа/выхода БПФ
- Мнимые части входного сигнала нулевые
- Выходной фаза-частотный сигнал дублируется, по факту из него нужно N/2+1 значений
- Оптимизации, основанные на вырожденном поворачивающем множителе
- Для ? со значениями 0 и 1 достаточно операции сложения
- Для ? со значением 0.7071 (угол 45 градусов) одна операция умножения может быть сэкономлена
- Оптимизации, основанные на особенностях потребителя БПФ
- БПФ может брать и выдавать данные в требуемом порядке, задаваемом в настройках пре-компилятора
- В БПФ может быть встроена нормализация или подгонка по окну входных/выходных данных
Почти все эти оптимизации БПФ применимы к методу гусеницы.
Возможны также архитектурно-машинные оптимизации. Например, элементарная БПФ операция в коде выше реализована как макрос XMAC. Она может быть распараллелена при помощи SIMD инструкций для x86-процессоров Intel AVX:
#define XMAC(c, a, wr, wi) vec1 = _mm_setr_pd(wr, wi); vec2 = *( __m128d*)&X2Y(a)->r; vec1 = _mm_hadd_pd(_mm_mul_pd(vec1, _mm_mul_pd(vec2, neg)), _mm_mul_pd(vec1, _mm_shuffle_pd(vec2, vec2, 1))); *( __m128d*)c = _mm_add_pd(*( __m128d*)a, vec1);Приведенный макрос поддерживает две операции с плавающей точкой одновременно при помощи 128-битных регистров. Но стоит взглянуть на гусеницу еще раз – из-за особенностей адресации ближайшие XMAC-и могут быть объединены вместе, например, для реализации при помощи 512-регистров(AVX3).
В завершении стоит сказать, что направлений развития пре-компилятора куда больше чем возможностей. Поэтому целью данной статьи является сбор дополнительных требований и выявления областей, где этот подход может быть полезен.
Источники
[1] Notes on Recursive FFT (Fast Fourier Transform) algorithm, Fall 2010, COSC 511, Susan Haynes
[2] Notes on the FFT, C. S. Burrus, Department of Electrical and Computer Engineering
Rice University, Houston, TX 77251-1892
[3] Caterpillar Implementation Based on Generated Code
[4] Some fft_caterpillar examples, https://github.com/r35382/fft_caterpillar
Комментарии (15)
leshabirukov
27.12.2017 17:32Так времянки-то где?
fronda Автор
27.12.2017 18:00-1Какие? C чем? И для чего нужны? Статья, собственно, выложена для сбора требований
RomanArzumanyan
27.12.2017 18:08+2Заявлено о достижении максимальной производительности, значит должно быть сравнение с prior art с подробнейшим разбором по косточкам.
Вы даже не сказали, что понимаете под максимальной производительностью: минимизация числа арифметических операций, или числа тактов процессора (какого, кстати?), минимизация числа гейтов для asic и т. п.
leshabirukov
27.12.2017 18:25Какие? C чем?
На любой на ваш выбор машине запустите ДПФ на 256 точек тысячу раз подряд, и замерьте время работы для каждого приведённого варианта.fronda Автор
27.12.2017 18:47Это мне поможет найти заказчика?
Может лучше научиться делать ассемблерную версию или с фиксированной точкой?
RomanArzumanyan
27.12.2017 18:03+1Самый известный БПФ алгоритм использует сложную логику перестановки адресов, которая в графическом виде напоминает крылья бабочки
Там проще всё. Плотную матрицу оператора преобразования заменяют произведением нескольких разреженных матриц, в число которых входит матрица перестановок. Вот и получаются бабочки. Сравнение с отцами-основателями (метод Чена или метод Кули-Тьюки) есть?
способы достижения максимальной производительности, превосходящей теоретическую
Хотелось бы видеть развёрнутое сравнение на эту тему, если заявлено в шапке статьи.fronda Автор
27.12.2017 18:38Не заявлено. Заявлена попытка.
Для сравнения с отцами основателями нужно решить что делать с w (поворотные множители вычисляются в run-time), т.е. преобразовывать конкретный код или получит откуда-нибудь референс-код. А грубо и так понятно — достаточно сравнить ассемблерный цикл там и ХМАС здесь.
В планах бенчмарки с fftw, но…
aamonster
27.12.2017 20:35Ага, и ещё заявлены "способы достижения максимальной производительности, превосходящей теоретическую" (посмеялся с человека, не знающего, что такое O(n), но пишущего статьи по оптимизации).
Зря вы так, люди зашли в надежде увидеть что-то интересное (я поначалу подумал, что вы просто ошиблись в процитированной фразе и имелось в виду "быстрее реализации в лоб", но вы даже это не проверили, увы) и обломались, сейчас вам минусов насуют...fronda Автор
27.12.2017 20:47Способы перечислены в главе «Дальнейшие оптимизации». — Тем кто проглядывает, а не читает.
Тем кто читает, скажу, что там нет новизны — все{ пардон, почти все} способы описаны в той или иной литературе. По отдельности. Часто, с очень сложной реализацией.
У меня их реализация ну очень проста. Хотя это конечно imho.aamonster
27.12.2017 21:57Я к тому, что "производительность, превосходящая теоретическую" в данном случае — это быстрее, чем за O(n*log(n)).
Причём O-нотация имеет смысл, когда n стремится к бесконечности, так что замеры для 16 или 64 отсчётов тут не годятся (хотя они интересны на практике — но для их оптимизации есть свои подходы… И тут может быть интересно исследование вида "взять алгоритм и доказать, что он реализует преобразование Фурье за минимально возможное количество операций", а потом построить ещё один алгоритм, который на реальном CPU работает быстрее за счёт минимизации cache misses. Но в любом случае в исследовании должны быть или замеры, или точный расчёт числа операций).RomanArzumanyan
28.12.2017 00:38За счёт суперскалярности и инструкций типа FMA, процессор может выполнить алгоритм за меньшее число инструкций / циклов, чем шагов в алгоритме. Фактическая производительность в смысле числа элементарных шагов будет выше теоретической.
Для БПФ O-нотация не имеет большого смысла, количество инструкций конечно и известно заранее, ветвлений нет, регистров нужно много, инструкции нужны широкие. В общем, всё это отлично ложится на архитектуру VLIW, которая и популярна в DSP. Время исполнения кода в таком случае будет известно ещё до запуска программы.
Промахов по кэшу тоже не особо должно быть — прочли несколько кэш-линий из массива, и муссируем эти числа пока не выполним преобразование. Если регистров много и / или преобразование маленького размера (например, преобразования 4х4, 8х8 в видеокодеках), то БПФ делается целиком в регистрах. Полкило прочли, долго считаем в регистрах, полкило записали — идеальный алгоритм для тщательной оптимизации. Очень мало в нём непредсказуемых факторов.aamonster
28.12.2017 02:26В главном (предсказуемость бпф, возможности его оптимизаций под железо)я с вами согласен, но:
- O-нотация имеет смысл. Но только для больших размеров окна (как, собственно, всегда и было), для малых количество операций считается в лоб (ну и там куча оптимизаций типа того, чему равны синусы/косинусы 90 и 45 градусов).
- Суперскалярность и т.п. не влияют на O-оценки — только на константу перед n log(n). Но, надеюсь, вы помните, что O(100 n)==O(n)?
- Раз уж мы умеем оценивать число операций (неважно, в виде O(...) или точно) — о какой "производительности, превосходящей теоретическую" идёт речь? В идеале производительность будет в точности равна теоретической. Можно лишь сравнить с другими реализациями по бенчмаркам — но именно этого автор не сделал, хотя это главное.
DustCn
Вообще то MAC это в AVX512 будет FMA (Fused Multiply-Add). Так что еще проще выйдет