[Из песочницы] Достижение максимальной производительности Быстрого Преобразования Фурье на основе управления данными
Сейчас проблема производительности приложений с цифровой обработкой сигналов решается использованием специализированного оборудования. При этом такому оборудованию требуются еще и оптимизированные под него программы.
Есть и несколько подходов для подгонки алгоритмов под машинную архитектуру, один из которых — управление данными (data driven). Каждый, кто сталкивался с ручным программированием данных, знает, что это дело это не простое. Однако в большинстве случаев можно спроектировать пре-компилятор, который существенно упростит задачу. В статье описана методика построения двух управляемых данными алгоритмов БПФ и способы достижения максимальной производительности, превосходящей теоретическую.
Немного математики
N временных отсчетов любого сигнала могут быть преобразованы в N фазо-частот при помощи Прямого Дискретного Преобразования Фурье.
Xk = n=0N-1∑ xn* ωNkn, k= 0…N-1, ωN= e-2πi/N, i = 2√-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
#include
#include
#include
#include
#define M_PI 3.14159265358979323846
using namespace std;
typedef complex w_type;
static vector fft(const vector &In)
{
int i = 0, wi = 0;
int n = In.size();
vector A(n / 2), B(n / 2), Out(n);
if (n == 1) {
return vector(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 At = fft(A);
vector 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 In(1 << ln);
std::transform(argv + 1, argv + argc, In.begin(),[&](const char* arg) {
return w_type(atof(arg), 0);
});
vector Out = fft(In);
for (vector::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; |
Логарифм по основанию 2 от размера сигнала определяет число рекурсий. Несколько стилей n2ln функции приведены ниже.
| Прототип | int n2ln( int n ) |
| С++ | return (int)floor( log(n) / log(2.0) ); |
| Классический С | int ln = 0; |
| Встроенный С | |
В последней реализации наихудший случай исполнения равен среднему случаю, который С-компилятор может реализовать при помощи 3-х сдвигов и 4-х условных переходов. Кроме того эта реализация может быть представлена в виде макроса для расчета значений констант во время компиляции.
БПФ на С
Классическая реализация рекурсивного БПФ на С похожа на С++ вариант
#include
#include
#include
#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
#include
#include
#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
#include
#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
#include
#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
