Быстрое преобразование Фурье (код на C++)

Файл test.cpp
/*
* Long multiplication with the FFT.
*
* Daniel Shved, MIPT, 2010.
* danshved [at] gmail.com
*/
#include "dft.h"
#include
#include
using namespace std;

const int base = 10;

/*
* Turns a string with a number into an array of complex numbers (digits)

*/
comp *make_polynomial(const string &str, int length)
{
comp *result = new comp[length];
for(int i = 0; i < str.length(); i++)
result[i] = comp((double)(str[str.length() - 1 - i] - '0'), 0.0);
for(int i = str.length(); i < length; i++)
result[i] = comp(0.0, 0.0);
return result;
}

/*
* This is used to turn complex numbers which are known to be almost real

* integers into integers.
*/
int round(comp x)
{
return (int)(x.real()+0.5);
}

/*
* Given a polynomial ``a'', returns the string with number a(10).

*/
string make_number(comp *a, int n)
{
// Turn coefficients into digits (carry)
for(int i = 0; i < n - 1; i++) {
a[i+1] += round(a[i]) / base;
a[i] = round(a[i]) % base;
}

// Determine length

int realLength = 0;
for(int i = 0; i < n; i++)
if(round(a[i]) != 0)
realLength = i+1;
if(realLength == 0)
return string("0");
// Make a string
string result(realLength, '0');
for(int i = 0; i < realLength; i++)
result[realLength - 1 - i] = '0' + (char)round(a[i]);
return result;
}

/*
* Long multiplication.

*/
string multiply(string number1, string number2)
{
// Determine the size of polynomials
int deg = number1.length() + number2.length();
int n = 1;
while(n<=deg) n*= 2;

// Turn numbers into complex polynomials

comp *a = make_polynomial(number1, n);
comp *b = make_polynomial(number2, n);

// Multiply polynomials with the FFT

comp *product = new comp[n];
convolution(a, b, product, n);

// Turn the polynomial back into a number

string result = make_number(product, n);

// Clean up and exit

delete a;
delete b;
delete product;
return result;
}

/*
* A test. Reads two non-negative integers from the input and prints their

* product.
*/
int main()
{
string a, b, res;
cin >> a >> b;
res = multiply(a, b);
cout << res << endl;
return 0;
}
БПФ на 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; cc->i = a.i + w.r * b.i + w.i * b.r;
Оптимизированный С	register double reg; reg = mac(a.r, w.r, b.r); cc->r = mac(reg, -w.i, b.i); reg = mac(a.i, w.r, b.i ); cc->i = mac(reg, w.i * b.r );

Логарифм по основанию 2 от размера сигнала определяет число рекурсий. Несколько стилей n2ln функции приведены ниже.

Прототип	int n2ln( int n )
С++	return (int)floor( log(n) / log(2.0) );
Классический С	int ln = 0; while (n >>= 1) ln++;
Встроенный С	return n/256 ?n/4048 ?n/16348 ?n/32768?15:14 :n/8096?13:12 :n/1024 ?n/2048?11:10 :n/512?9:8 :n/16 ?n/64 ?n/128?7:6 :n/32?5:4 :n/4 ?n/8?3:2 :n/2?1: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). 

В завершении стоит сказать, что направлений развития пре-компилятора куда больше чем возможностей. Поэтому целью данной статьи является сбор дополнительных требований и выявления областей, где этот подход может быть полезен.
Дискретное преобразование Фурье на С++. Округляет - в чем дело?


#include "stdafx.h"
#include
#include
#include
#include
#include


const double pi2 = 6.283185307179586;


using namespace std;


int k, i, n, K, N, kol;
double S[1000][2], ugol, Nd, kold, kd, nd;
double a[1000], b[1000], c[1000], fi[1000], m, d, x;
double coord[1000][2];


void main(void)
{


setlocale(LC_CTYPE, "rus");//используем русские символы


int kol;
cout << "Количество точек=";//ввели количество точек, на которое разбиваем 2п
cin >> kol;


kold = kol;


d = 6.283185307179586 / (kold - 1.0);//разбили интервал на части
x = 0.0;
for (i = 0; i <(kol+1); i++)
{
coord[i][0] = i;
coord[i][1] = sin(x);
x = x + d;
}


fstream f;//обозначили поток f


f.open("C:\\coord\\coord.txt", fstream::out);//открыли файл
for (i = 0; i < kol; i++)
f << coord[i][0] << " " << coord[i][1]<<"\n";
f.close();


cout << "Число гармоник=";
cin >> K;//число гармоник, на которые разобьем интервал
N = 0;//переменная для кол-ва символов. пока =0
fstream i_f;//обозначили поток f
i_f.open("C:\\coord\\coord.txt", fstream::in);//открыли файл


while (!i_f.eof())//пока не конец файла
{
i_f >> S[N][0] >> S[N][1];//записываем в файл массив из данных координат
N++;//число симв в файле
}
i_f.close();


N=N-1;
Nd=N;//привели int в double
fstream o_f;//обозначили поток f


//считаем а и b
for (k = 0; k {kd=k;
for (n = 0; n {nd=n;
ugol = (pi2*kd*nd)/Nd;


a[k] = a[k] + (2.0/Nd)*S[n][1]*cos(ugol);//*
b[k] = b[k] + (2.0/Nd)*S[n][1]*sin(ugol);
}
c[k] = sqrt(pow(a[k],2) + pow(b[k],2));
}
o_f.close();


fstream o_f;//обозначили поток f
o_f.open("C:\\coord\\ачх.txt", fstream::out);//открыли файл
for (i = 0; i {
o_f< }
o_f.close();


system("pause");
}