желательно оформить с пояснениями ...
C++:
// Быстрое преобразование Фурье. Метод разложения на множители.
// Описание: https://ru.dsplib.org/content/fft_composite/fft_composite.html
#include "math.h"
#include <map>
#include <vector>
#define M_PI 3.14159265358979323846
class ShortComplex2
{
public:
double re, im;
ShortComplex2() { re = 0; im = 0; }
ShortComplex2( const ShortComplex2& co ) { re = co.re; im = co.im; }
ShortComplex2( double re_, double im_ ) { re = re_; im = im_; }
ShortComplex2& operator=( const ShortComplex2& co ) { re = co.re; im = co.im; return *this; }
ShortComplex2& operator+=( const ShortComplex2& co ) { re += co.re; im += co.im; return *this; }
ShortComplex2 operator+( const ShortComplex2& co ) { ShortComplex2 co2( re + co.re, im + co.im ); return co2; }
ShortComplex2 operator-( const ShortComplex2& co ) { ShortComplex2 co2( re - co.re, im - co.im ); return co2; }
ShortComplex2 operator*( const ShortComplex2& co ) { ShortComplex2 co2( re * co.re - im * co.im, re * co.im + im * co.re ); return co2; }
};
/*
* Быстрое преобразование Фурье. Метод разложения на простые множители.
*
* pCoIn - Массив (указатель) исходной последовательности
* pCoOut - Массив (указатель) для результата
* pBufSize - Массив (указатель) для работы метода длины Size.
* Size - Размер последовательности
* bDir - true: для прямого преобразования, false: для обратного преобразования
*/
void SmpNm_FFT( ShortComplex2 *pCoIn, ShortComplex2 *pCoOut, ShortComplex2 *pBufSize, int Size, bool bDir );
/*
* Обычное дискретное преобразование Фурье
*/
void DPF( ShortComplex2 *pCoIn, ShortComplex2 *pCoOut, int Size, bool bDir );
/*
* Пример программы с тестом и сравнением с обычным ДПФ
*/
int _tmain(int argc, _TCHAR* argv[]);
/*
* строит вектор простых чисел до данного числа
*/
void CheckSmplNumbs( std::vector<size_t> &vSmplNum, size_t UpLevel );
/*
* Раскладывает число на простые множители
*/
size_t DecomposeNum( size_t uV, std::vector<size_t> &vSmplM );
/*
* Рабочая функция для БПФ
*/
void GNext( int N, int ppPred, int pLast, ShortComplex2 *pF1, ShortComplex2 *pF2, std::vector<ShortComplex2> &vW );
void DPF( ShortComplex2 *pCoIn, ShortComplex2 *pCoOut, int Size, bool bDir )
{
ShortComplex2 c1;
double d2 = 2. * M_PI / Size, d3;
int i, k;
if( Size < 1 )
return;
if( bDir ) d2 = - d2;
for( k = 0; k < Size; k++ )
{
ShortComplex2 & pCoRez = pCoOut[ k ];
pCoRez = pCoIn[ 0 ];
for( i = 1; i < Size; i++ )
{
d3 = k * i * d2;
c1.re = cos( d3 );
c1.im = sin( d3 );
pCoRez.re += ( c1.re * pCoIn[ i ].re - c1.im * pCoIn[ i ].im );
pCoRez.im += ( c1.re * pCoIn[ i ].im + c1.im * pCoIn[ i ].re );
}
}
if( bDir == 0 )
{
for( k = 0; k < Size; k++ )
{
pCoOut[ k ].re /= Size;
pCoOut[ k ].im /= Size;
}
}
}
void CheckSmplNumbs( std::vector<size_t> &vSmplNum, size_t UpLevel )
{
size_t i, k;
if( vSmplNum.empty() )
{
vSmplNum.reserve( 64 );
vSmplNum.resize( 3 );
vSmplNum[ 0 ] = 2;
vSmplNum[ 1 ] = 3;
vSmplNum[ 2 ] = 5;
}
for( k = vSmplNum.back() + 2; ; k += 2 )
{
for( i = 0; ; i++ )
{
if( i >= vSmplNum.size() )
{
vSmplNum.push_back( k );
if( k > UpLevel )
return;
break;
}
if( ( k % vSmplNum[ i ] ) == 0 )
break;
}
}
}
size_t DecomposeNum( size_t uV, std::vector<size_t> &vSmplM )
{
static std::vector<size_t> vSmplNum;
static std::map< size_t, std::vector<size_t> > mpSmplNum;
std::vector<size_t> & vSmplNum2 = mpSmplNum[ uV ];
std::map<size_t,size_t> mp;
size_t i, j = 0;
if( i = vSmplNum2.size() )
{
vSmplM = vSmplNum2;
return i;
}
CheckSmplNumbs( vSmplNum, uV );
m_begin:;
for( i = 0; ; i++ )
{
if( uV == vSmplNum[ i ] )
{
mp[ vSmplNum[ i ] ] ++;
j ++;
goto m_end;
}
if( ( uV % vSmplNum[ i ] ) == 0 )
{
mp[ vSmplNum[ i ] ] ++;
uV /= vSmplNum[ i ];
j ++;
goto m_begin;
}
}
m_end:;
vSmplM.resize( j );
i = 0;
for( std::map<size_t,size_t>::reverse_iterator P = mp.rbegin(); P != mp.rend(); P++ )
{
for( j = 0; j < P->second; j++ )
vSmplM[ i++ ] = P->first;
}
vSmplNum2 = vSmplM;
return i;
}
void GNext( int N, int ppPred, int pLast, ShortComplex2 *pF1, ShortComplex2 *pF2, std::vector<ShortComplex2> &vW )
{
int N2, n, p2, l1, N3, k, k2, l2;
p2 = ppPred * pLast;
N2 = N / ppPred;
N3 = N / p2;
for( n = 0; n < N2; n++ )
{
k = ( n % N3 ) * p2;
for( l1 = 0; l1 < ppPred; l1++ )
{
ShortComplex2 & pCoRez = pF1[ ppPred * n + l1 ];
pCoRez = pF2[ k2 = k + l1 ];
for( l2 = 1, k2 += ppPred; l2 < pLast; l2++, k2 += ppPred )
pCoRez += pF2[ k2 ] * vW[ ( n * ppPred * l2 ) % N ];
}
}
}
void SmpNm_FFT( ShortComplex2 *pCoIn, ShortComplex2 *pCoOut, ShortComplex2 *pBufSize, int Size, bool bDir )
{
std::map<int,std::vector<ShortComplex2>> mpW;
std::vector<size_t> vSmplM;
ShortComplex2 *pBuf[ 2 ];
int ppPred, p2, i;
if( Size < 2 )
return;
std::vector<ShortComplex2> &vW = mpW[ Size ];
if( vW.size() < (size_t)Size )
{
double d;
vW.resize( Size );
for( i = 0; i < Size; i++ )
{
d = - M_PI * ( 2 * i ) / Size;
vW[ i ].re = cos( d );
vW[ i ].im = sin( d );
}
}
DecomposeNum( Size, vSmplM );
i = vSmplM.size() - 1;
p2 = vSmplM[ i ];
ppPred = Size / p2;
if( bDir )
{
pBuf[ 0 ] = pCoOut;
pBuf[ 1 ] = pBufSize;
GNext( Size, ppPred, p2, pBuf[ i & 1 ], pCoIn, vW );
}
else
{
pBuf[ 0 ] = pBufSize;
pBuf[ 1 ] = pCoOut;
ShortComplex2 *pBuf2 = pBuf[ ( i + 1 ) & 1 ];
for( int j = 0; j < Size; j++ ) { pBuf2[ j ].re = pCoIn[ j ].re; pBuf2[ j ].im = - pCoIn[ j ].im; }
GNext( Size, ppPred, p2, pBuf[ i & 1 ], pBuf2, vW );
}
while( i-- > 0 )
{
p2 = vSmplM[ i ];
ppPred /= p2;
GNext( Size, ppPred, p2, pBuf[ i & 1 ], pBuf[ ( i + 1 ) & 1 ], vW );
}
if( ! bDir )
{
for( int j = 0; j < Size; j++ ) { pCoOut[ j ].re = pBufSize[ j ].re / Size; pCoOut[ j ].im = - pBufSize[ j ].im / Size; }
}
}
int _tmain(int argc, _TCHAR* argv[])
{
ShortComplex2 dSpec[2048], dSpec2[2048], dSpec3[2048], dSpecNew2[2048], dSpecNew2Buf[2048], dSpecNew3[2048];
DWORD dw1 = 0, dw2 = 0;
double dF1, DT, b;
int i, sz;
ZeroMemory(dSpec,sizeof(dSpec));
ZeroMemory(dSpec2,sizeof(dSpec2));
ZeroMemory(dSpec3,sizeof(dSpec3));
ZeroMemory(dSpecNew2,sizeof(dSpecNew2));
ZeroMemory(dSpecNew2Buf,sizeof(dSpecNew2Buf));
ZeroMemory(dSpecNew3,sizeof(dSpecNew2Buf));
#if 0
sz = 15; // 7 * 11;
DT = 1. / 4000.;
dF1 = 1. / ( sz * DT );
b = 1. / DT;
double f1 = 500, f2 = 1500, A1 = 1, A2 = 2;
double w1 = 2*M_PI*f1*DT; // Угловые частоты
double w2 = 2*M_PI*f2*DT;
for( i = 0; i < sz; i++ )
{
dSpec[ i ].re = A1*cos(w1*i+M_PI/4.)+A2*cos(w2*i+M_PI/8.);
}
#endif
sz = 1722; // 1025; // 7 * 11;
DT = 1. / 200;
dF1 = 1. / ( sz * DT );
b = 1. / DT;
double f1 = 2, f2 = 4, A1 = 1, A2 = 2;
double w1 = 2*M_PI*f1*DT; // Угловые частоты
double w2 = 2*M_PI*f2*DT;
for( i = 0; i < sz; i++ )
{
dSpec[ i ].re = A1*cos(w1*i+M_PI/4.)+A2*cos(w2*i+M_PI/8.);
}
dw1 = GetTickCount();
SmpNm_FFT( dSpec, dSpecNew2, dSpecNew2Buf, sz, true );
for( i = 0; i < sz; i++ )
{
dSpecNew2[ i ].re *= DT;
dSpecNew2[ i ].im *= DT;
}
SmpNm_FFT( dSpecNew2, dSpecNew3, dSpecNew2Buf, sz, false );
for( i = 0; i < sz; i++ )
{
dSpecNew3[ i ].re *= b;
dSpecNew3[ i ].im *= b;
}
dw1 = GetTickCount() - dw1;
printf( "\n\n\n******************************************************\n\n\n");
dw2 = GetTickCount();
DPF( dSpec, dSpec2, sz, true );
for( i = 0; i < sz; i++ )
{
dSpec2[ i ].re *= DT;
dSpec2[ i ].im *= DT;
}
DPF( dSpec2, dSpec3, sz, false );
for( i = 0; i < sz; i++ )
{
dSpec3[ i ].re *= b;
dSpec3[ i ].im *= b;
}
dw2 = GetTickCount() - dw2;
for( i = 0; i < sz; i++ )
printf( "\n %3i. %8.10lf -> : %8.5lf : [%8.5lf,%8.5lf] -> %8.10lf, %8.5lf -|- [%8.5lf,%8.5lf] -> %8.10lf, %8.5lf",
i + 1, dSpec[ i ].re, dF1 * i, dSpec2[ i ].re, dSpec2[ i ].im, dSpec3[ i ].re, dSpec3[ i ].im, dSpecNew2[ i ].re, dSpecNew2[ i ].im, dSpecNew3[ i ].re, dSpecNew3[ i ].im );
printf("\n\n Time: %u - %u", dw1, dw2 );
scanf("%i",&i);
return 0;
}