FFT变换 C源代码

FFT变换 C源代码

FFT C source code (Simple radix-2)

void fft_float (    unsigned  NumSamples,    int       InverseTransform,    float    *RealIn,    float    *ImagIn,    float    *RealOut,    float    *ImagOut ){    unsigned NumBits;    /* Number of bits needed to store indices */    unsigned i, j, k, n;    unsigned BlockSize, BlockEnd;    double angle_numerator = 2.0 * DDC_PI;    double tr, ti;     /* temp real, temp imaginary */    if ( !IsPowerOfTwo(NumSamples) )    {        fprintf (            stderr,            "Error in fft():  NumSamples=%u is not power of two\n",            NumSamples );        exit(1);    }    if ( InverseTransform )        angle_numerator = -angle_numerator;    CHECKPOINTER ( RealIn );    CHECKPOINTER ( RealOut );    CHECKPOINTER ( ImagOut );    NumBits = NumberOfBitsNeeded ( NumSamples );    /*    **   Do simultaneous data copy and bit-reversal ordering into outputs...    */    for ( i=0; i < NumSamples; i++ )    {        j = ReverseBits ( i, NumBits );        RealOut[j] = RealIn;        ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn;    }    /*    **   Do the FFT itself...    */    BlockEnd = 1;    for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )    {        double delta_angle = angle_numerator / (double)BlockSize;        double sm2 = sin ( -2 * delta_angle );        double sm1 = sin ( -delta_angle );        double cm2 = cos ( -2 * delta_angle );        double cm1 = cos ( -delta_angle );        double w = 2 * cm1;        double ar[3], ai[3];        double temp;        for ( i=0; i < NumSamples; i += BlockSize )        {            ar[2] = cm2;            ar[1] = cm1;            ai[2] = sm2;            ai[1] = sm1;            for ( j=i, n=0; n < BlockEnd; j++, n++ )            {                ar[0] = w*ar[1] - ar[2];                ar[2] = ar[1];                ar[1] = ar[0];                ai[0] = w*ai[1] - ai[2];                ai[2] = ai[1];                ai[1] = ai[0];                k = j + BlockEnd;                tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];                ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];                RealOut[k] = RealOut[j] - tr;                ImagOut[k] = ImagOut[j] - ti;                RealOut[j] += tr;                ImagOut[j] += ti;            }        }        BlockEnd = BlockSize;    }    /*    **   Need to normalize if inverse transform...    */    if ( InverseTransform )    {        double denom = (double)NumSamples;        for ( i=0; i < NumSamples; i++ )        {            RealOut /= denom;            ImagOut /= denom;        }    }}

int IsPowerOfTwo ( unsigned x ){    if ( x < 2 )        return FALSE;    if ( x & (x-1) )        // Thanks to 'byang' for this cute trick!        return FALSE;    return TRUE;}

unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo ){    unsigned i;    if ( PowerOfTwo < 2 )    {        fprintf (            stderr,            ">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too small.\n",            PowerOfTwo );        exit(1);    }    for ( i=0; ; i++ )    {        if ( PowerOfTwo & (1 << i) )            return i;    }}

unsigned ReverseBits ( unsigned index, unsigned NumBits ){    unsigned i, rev;    for ( i=rev=0; i < NumBits; i++ )    {        rev = (rev << 1) | (index & 1);        index >>= 1;    }    return rev;}

double Index_to_frequency ( unsigned NumSamples, unsigned Index ){    if ( Index >= NumSamples )        return 0.0;    else if ( Index <= NumSamples/2 )        return (double)Index / (double)NumSamples;    return -(double)(NumSamples-Index) / (double)NumSamples;}

转载于:https://www.cnblogs.com/xtrgm623/archive/2008/12/05/1348332.html

相关资源：图像处理 FFT快速傅里叶变换 C代码