Быстрое преобразование Фурье по составному основанию

#include<iostream> 
#include<cmath> 
#include<stdlib.h> 
#include<complex> 
 
//template< class T > class complex; 
//template<> class complex<double>; 
 
using namespace std; 
 
#define Ngarm 30 
#define PI 3.14159265 
#define Freq 1000 
#define Amp 1 
#define Fd (2*Freq) 
#define p1 2 
#define p2 3 
#define p3 5 
 
 
 
const complex<long double> J (0, 1); 
 
complex<long double> Wnk(int M, int p, int q) 
{ 
return(exp(-J * (long double)(2 * PI / Ngarm * M * p * q))); 
} 
 
int main(void) 
{ 
int a = 0; 
 
long double Sig[p1][p2*p3]; //Выделение памяти под массив значений сигнала 
//complex<double> FFT1[]; 
for (int y = 0; y < p1; y++) //фиксируем строку 
{ 
for (int x = 0; x < p2*p3; x++) //заполнение по столбцам 
{ 
Sig[y][x] = Amp * sin(2 * a * PI*Freq / Fd); 
a++; 
//cout « Sig[i][k] « '/t'; 
} 
//cout « endl; 
} 
 
complex<long double> FFT2[p1][p2*p3]; //Выделяем память под массив двухточечных дпф всех стольцов 
 
for (int x = 0; x < p2*p3; x++) //фиксируем столбец 
{ 
for (int s = 0; s < p1; s++) //заносим значения в строчки FFT5 строка
{ 
for (int y = 0; y < p1; y++) //повторяем сумму 
{ 
FFT2[s][x] = FFT2[s][x] + Sig[x][y] * Wnk(p2*p3, s, y); 
} 
FFT2[s][x] = FFT2[s][x] * Wnk(1, s, x); //домножаем на поворачивающие множители 
} 
} 
 
complex<long double> FFTFull[p1][p2][p3]; //выделяем память под 3-массив,чтобы далее разложить ДПФ строки в матрицы по базису p2 и p3 
 
for (int y = 0; y < p1; y++) //поехали по строчкам 
{ 
for (int x = 0; x < p2; x++) //столбцы 
{ 
for (int z = 0; z < p3; z++) //в глубину 
{ 
int n = p3 * x + z; //номер элемента в матрице (ширина*глубина) по порядку строки FFT5 
FFTFull[y][x][z] = FFT2[y][n]; //заполнение 
} 
} 
} 
 
//complex<long double> FFT3[p1][p2][p3]; //выделяем память под ДПФ трёхточечных стобцов 
 
for (int y = 0; y < p1; y++) //по строкам срезали пласт рассматриваемой матрицы p2xp3 
{ 
for (int x = 0; x < p3; x++) //по столбцам p2xp3 
{ 
for (int s = 0; s < p2; s++) //Что будем заносить в FFT3 
{ 
for (int z = 0; z < p2; z++) //суммируем по строчкам p2xp3 
{ 
FFTFull[y][x][s] = FFTFull[y][x][s] + FFTFull[y][x][z] * Wnk(p3, s, z); 
} 
FFTFull[y][x][s] = FFTFull[y][x][s] * Wnk(1, x, s); //домножаем на поворачивающие множители 
} 
} 
} 
 
//complex<long double> FFT5[p1][p2][p3]; //выделяем память под пятиточечнуб ДПФ 
 
for (int y = 0; y < p1; y++) 
{ 
for (int z = 0; z < p2; z++) 
{ 
for (int s = 0; s < p3; s++) 
{ 
for (int x = 0; x < p3; x++) 
{ 
FFTFull[y][s][z] = FFTFull[y][s][z] + FFTFull[y][x][z] * Wnk(p2, s, x); 
} 
} 
} 
} 
 
return 0; 
}