#define mag_sqrd(re,im) (re*re+im*im)
#define Decibels(re,im) (20.0 * log10(double(mag_sqrd(re,im))))


inline double GetFrequencyIntensity(double re, double im)
{
    return sqrt((re*re)+(im*im));
}


void fft_double (unsigned int p_nSamples, bool p_bInverseTransform, 
    double *p_lpRealIn, double *p_lpImagIn, 
    double *p_lpRealOut, double *p_lpImagOut)
{

    if(!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;


    unsigned int NumBits;
    unsigned int i, j, k, n;
    unsigned int BlockSize, BlockEnd;

    double angle_numerator = 2.0 * PI;
    double tr, ti;

    if( !IsPowerOfTwo(p_nSamples) )
    {
        return;
    }

    if( p_bInverseTransform ) angle_numerator = -angle_numerator;

    NumBits = NumberOfBitsNeeded ( p_nSamples );


    for( i=0; i < p_nSamples; i++ )
    {
        j = ReverseBits ( i, NumBits );
        p_lpRealOut[j] = p_lpRealIn[i];
        p_lpImagOut[j] = (p_lpImagIn == NULL) ? 0.0 : p_lpImagIn[i];
    }


    BlockEnd = 1;
    for( BlockSize = 2; BlockSize <= p_nSamples; 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];

        for( i=0; i < p_nSamples; 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]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
                ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];

                p_lpRealOut[k] = p_lpRealOut[j] - tr;
                p_lpImagOut[k] = p_lpImagOut[j] - ti;

                p_lpRealOut[j] += tr;
                p_lpImagOut[j] += ti;

            }
        }

        BlockEnd = BlockSize;

    }


    if( p_bInverseTransform )
    {
        double denom = (double)p_nSamples;

        for ( i=0; i < p_nSamples; i++ )
        {
            p_lpRealOut[i] /= denom;
            p_lpImagOut[i] /= denom;
        }
    }

}



///////////////////////////////////////////////////////////
// check is a number is a power of 2
///////////////////////////////////////////////////////////

bool IsPowerOfTwo( unsigned int p_nX )
{

    if( p_nX < 2 ) return false;

    if( p_nX & (p_nX-1) ) return false;

    return true;

}


///////////////////////////////////////////////////////////
// return needed bits for fft
///////////////////////////////////////////////////////////

unsigned int NumberOfBitsNeeded( unsigned int p_nSamples )
{

    int i;

    if( p_nSamples < 2 )
    {
        return 0;
    }

    for ( i=0; ; i++ )
    {
        if( p_nSamples & (1 << i) ) return i;
    }

}



///////////////////////////////////////////////////////////
// ?
///////////////////////////////////////////////////////////

unsigned int ReverseBits(unsigned int p_nIndex, unsigned int p_nBits)
{

    unsigned int i, rev;

    for(i=rev=0; i < p_nBits; i++)
    {
        rev = (rev << 1) | (p_nIndex & 1);
        p_nIndex >>= 1;
    }

    return rev;

}



///////////////////////////////////////////////////////////
// return a frequency from the basefreq and num of samples
///////////////////////////////////////////////////////////

double Index_to_frequency(unsigned int p_nBaseFreq, 
    unsigned int p_nSamples, unsigned int p_nIndex)
{

    if(p_nIndex >= p_nSamples)
    {
        return 0.0;
    }
    else if(p_nIndex <= p_nSamples/2)
    {
        return ( (double)p_nIndex / (double)p_nSamples * p_nBaseFreq );
    }
    else
    {
        return ( -(double)(p_nSamples-p_nIndex) / (double)p_nSamples * p_nBaseFreq );
    }

}




main {{{{{{{{{{{{{
    DWORD nCount = 0;
    for (DWORD dw = 0; dw < FFT_LEN; dw++)
    {
        {
            //copy audio signal to fft real component for left channel
            finleft[nCount] = (double)((short*)pwh->lpData)[dw++];
            //copy audio signal to fft real component for right channel
            finright[nCount++] = (double)((short*)pwh->lpData)[dw];
        }
    }
    //Perform FFT on left channel
    fft_double(FFT_LEN/2,0,finleft,NULL,fout,foutimg);
    float re,im;
    for(int i=0;i < FFT_LEN/2;i++)
    {
        re = fout[i];
        im = foutimg[i];
        //get frequency intensity and scale to 0..256 range
        fdraw[i]=(GetFrequencyIntensity(re,im))/256;
    }

    // Perform FFT on right channel
    fft_double(FFT_LEN/2,0,finright,NULL,fout,foutimg);
    for(i=0;i < FFT_LEN/2;i++)
    {
        re = fout[i];
        im = foutimg[i];
        //get frequency intensity and scale to 0..256 range
        fdraw[i]=(GetFrequencyIntensity(re,im))/256;
}}}}}}}}}}}}}}}}}}}}}