Jump to content

    
Sign in to follow this  
f3434s

Быстрое преобразование Фурье. Метод разложения на множители.

Recommended Posts

В архиве исходный код и описание метода.

 

Код:

Скрытый текст



#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

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Sign in to follow this