f3434s 0 23 апреля, 2021 Опубликовано 23 апреля, 2021 · Жалоба В архиве исходный код и описание метода. Код: Скрытый текст #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; } Furie.rar Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
kovigor 5 23 апреля, 2021 Опубликовано 23 апреля, 2021 · Жалоба Вы бы хоть указали, из какой книги взята статья во вложении ... Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
f3434s 0 23 апреля, 2021 Опубликовано 23 апреля, 2021 · Жалоба Алгоритм БПФ составной длины Там есть и список литературы. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться