Subversion Repositories Games.Chess Giants

/*-========================================================================-_

 |                                 - XDSP -                                 |

 |        Copyright (c) Microsoft Corporation.  All rights reserved.        |

 |~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~|

 |PROJECT: XDSP                         MODEL:   Unmanaged User-mode        |

 |VERSION: 1.1                          EXCEPT:  No Exceptions              |

 |CLASS:   N / A                        MINREQ:  WinXP, Xbox360             |

 |BASE:    N / A                        DIALECT: MSC++ 14.00                |

 |>------------------------------------------------------------------------<|

 | DUTY: DSP functions with CPU extension specific optimizations            |

 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^

  NOTES:

    1.  Definition of terms:

            DSP: Digital Signal Processing.

            FFT: Fast Fourier Transform.

    2.  All buffer parameters must be 16-byte aligned.

    3.  All FFT functions support only FLOAT32 mono audio.                  */

#pragma once

//--------------<D-E-F-I-N-I-T-I-O-N-S>-------------------------------------//

#include <windef.h> // general windows types

#include <math.h>   // trigonometric functions

#if defined(_XBOX)  // SIMD intrinsics

    #include <ppcintrinsics.h>

#else

    #include <emmintrin.h>

#endif

//--------------<M-A-C-R-O-S>-----------------------------------------------//

// assertion

#if !defined(DSPASSERT)

    #if DBG

        #define DSPASSERT(exp) if (!(exp)) { OutputDebugStringA("XDSP ASSERT: " #exp ", {" __FUNCTION__ "}\n"); __debugbreak(); }

    #else

        #define DSPASSERT(exp) __assume(exp)

    #endif

#endif

// true if n is a power of 2

#if !defined(ISPOWEROF2)

    #define ISPOWEROF2(n) ( ((n)&((n)-1)) == 0 && (n) != 0 )

#endif

//--------------<H-E-L-P-E-R-S>---------------------------------------------//

namespace XDSP {

#pragma warning(push)

#pragma warning(disable: 4328 4640) // disable "indirection alignment of formal parameter", "construction of local static object is not thread-safe" compile warnings

// Helper functions, used by the FFT functions.

// The application need not call them directly.

    // primitive types

    typedef __m128 XVECTOR;

    typedef XVECTOR& XVECTORREF;

    typedef const XVECTOR& XVECTORREFC;

    // Parallel multiplication of four complex numbers, assuming

    // real and imaginary values are stored in separate vectors.

    __forceinline void vmulComplex (__out XVECTORREF rResult, __out XVECTORREF iResult, __in XVECTORREFC r1, __in XVECTORREFC i1, __in XVECTORREFC r2, __in XVECTORREFC i2)

    {

        // (r1, i1) * (r2, i2) = (r1r2 - i1i2, r1i2 + r2i1)

        XVECTOR vi1i2 = _mm_mul_ps(i1, i2);

        XVECTOR vr1r2 = _mm_mul_ps(r1, r2);

        XVECTOR vr1i2 = _mm_mul_ps(r1, i2);

        XVECTOR vr2i1 = _mm_mul_ps(r2, i1);

        rResult = _mm_sub_ps(vr1r2, vi1i2); // real:      (r1*r2 - i1*i2)

        iResult = _mm_add_ps(vr1i2, vr2i1); // imaginary: (r1*i2 + r2*i1)

    }

    __forceinline void vmulComplex (__inout XVECTORREF r1, __inout XVECTORREF i1, __in XVECTORREFC r2, __in XVECTORREFC i2)

    {

        // (r1, i1) * (r2, i2) = (r1r2 - i1i2, r1i2 + r2i1)

        XVECTOR vi1i2 = _mm_mul_ps(i1, i2);

        XVECTOR vr1r2 = _mm_mul_ps(r1, r2);

        XVECTOR vr1i2 = _mm_mul_ps(r1, i2);

        XVECTOR vr2i1 = _mm_mul_ps(r2, i1);

        r1 = _mm_sub_ps(vr1r2, vi1i2); // real:      (r1*r2 - i1*i2)

        i1 = _mm_add_ps(vr1i2, vr2i1); // imaginary: (r1*i2 + r2*i1)

    }

    // Radix-4 decimation-in-time FFT butterfly.

    // This version assumes that all four elements of the butterfly are

    // adjacent in a single vector.

    //

    // Compute the product of the complex input vector and the

    // 4-element DFT matrix:

    //     | 1  1  1  1 |    | (r1X,i1X) |

    //     | 1 -j -1  j |    | (r1Y,i1Y) |

    //     | 1 -1  1 -1 |    | (r1Z,i1Z) |

    //     | 1  j -1 -j |    | (r1W,i1W) |

    //

    // This matrix can be decomposed into two simpler ones to reduce the

    // number of additions needed. The decomposed matrices look like this:

    //     | 1  0  1  0 |    | 1  0  1  0 |

    //     | 0  1  0 -j |    | 1  0 -1  0 |

    //     | 1  0 -1  0 |    | 0  1  0  1 |

    //     | 0  1  0  j |    | 0  1  0 -1 |

    //

    // Combine as follows:

    //          | 1  0  1  0 |   | (r1X,i1X) |         | (r1X + r1Z, i1X + i1Z) |

    // Temp   = | 1  0 -1  0 | * | (r1Y,i1Y) |       = | (r1X - r1Z, i1X - i1Z) |

    //          | 0  1  0  1 |   | (r1Z,i1Z) |         | (r1Y + r1W, i1Y + i1W) |

    //          | 0  1  0 -1 |   | (r1W,i1W) |         | (r1Y - r1W, i1Y - i1W) |

    //

    //          | 1  0  1  0 |   | (rTempX,iTempX) |   | (rTempX + rTempZ, iTempX + iTempZ) |

    // Result = | 0  1  0 -j | * | (rTempY,iTempY) | = | (rTempY + iTempW, iTempY - rTempW) |

    //          | 1  0 -1  0 |   | (rTempZ,iTempZ) |   | (rTempX - rTempZ, iTempX - iTempZ) |

    //          | 0  1  0  j |   | (rTempW,iTempW) |   | (rTempY - iTempW, iTempY + rTempW) |

    __forceinline void ButterflyDIT4_1 (__inout XVECTORREF r1, __inout XVECTORREF i1)

    {

        // sign constants for radix-4 butterflies

        const static XVECTOR vDFT4SignBits1 = { 0.0f, -0.0f,  0.0f, -0.0f };

        const static XVECTOR vDFT4SignBits2 = { 0.0f,  0.0f, -0.0f, -0.0f };

        const static XVECTOR vDFT4SignBits3 = { 0.0f, -0.0f, -0.0f,  0.0f };

        // calculating Temp

        XVECTOR rTemp = _mm_add_ps( _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(1, 1, 0, 0)),                               // [r1X| r1X|r1Y| r1Y] +

                                    _mm_xor_ps(_mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 3, 2, 2)), vDFT4SignBits1) ); // [r1Z|-r1Z|r1W|-r1W]

        XVECTOR iTemp = _mm_add_ps( _mm_shuffle_ps(i1, i1, _MM_SHUFFLE(1, 1, 0, 0)),                               // [i1X| i1X|i1Y| i1Y] +

                                    _mm_xor_ps(_mm_shuffle_ps(i1, i1, _MM_SHUFFLE(3, 3, 2, 2)), vDFT4SignBits1) ); // [i1Z|-i1Z|i1W|-i1W]

        // calculating Result

        XVECTOR rZrWiZiW = _mm_shuffle_ps(rTemp, iTemp, _MM_SHUFFLE(3, 2, 3, 2));       // [rTempZ|rTempW|iTempZ|iTempW]

        XVECTOR rZiWrZiW = _mm_shuffle_ps(rZrWiZiW, rZrWiZiW, _MM_SHUFFLE(3, 0, 3, 0)); // [rTempZ|iTempW|rTempZ|iTempW]

        XVECTOR iZrWiZrW = _mm_shuffle_ps(rZrWiZiW, rZrWiZiW, _MM_SHUFFLE(1, 2, 1, 2)); // [rTempZ|iTempW|rTempZ|iTempW]

        r1 = _mm_add_ps( _mm_shuffle_ps(rTemp, rTemp, _MM_SHUFFLE(1, 0, 1, 0)), // [rTempX| rTempY| rTempX| rTempY] +

                         _mm_xor_ps(rZiWrZiW, vDFT4SignBits2) );                // [rTempZ| iTempW|-rTempZ|-iTempW]

        i1 = _mm_add_ps( _mm_shuffle_ps(iTemp, iTemp, _MM_SHUFFLE(1, 0, 1, 0)), // [iTempX| iTempY| iTempX| iTempY] +

                         _mm_xor_ps(iZrWiZrW, vDFT4SignBits3) );                // [iTempZ|-rTempW|-iTempZ| rTempW]

    }

    // Radix-4 decimation-in-time FFT butterfly.

    // This version assumes that elements of the butterfly are

    // in different vectors, so that each vector in the input

    // contains elements from four different butterflies.

    // The four separate butterflies are processed in parallel.

    //

    // The calculations here are the same as the ones in the single-vector

    // radix-4 DFT, but instead of being done on a single vector (X,Y,Z,W)

    // they are done in parallel on sixteen independent complex values.

    // There is no interdependence between the vector elements:

    // | 1  0  1  0 |    | (rIn0,iIn0) |               | (rIn0 + rIn2, iIn0 + iIn2) |

    // | 1  0 -1  0 | *  | (rIn1,iIn1) |  =   Temp   = | (rIn0 - rIn2, iIn0 - iIn2) |

    // | 0  1  0  1 |    | (rIn2,iIn2) |               | (rIn1 + rIn3, iIn1 + iIn3) |

    // | 0  1  0 -1 |    | (rIn3,iIn3) |               | (rIn1 - rIn3, iIn1 - iIn3) |

    //

    //          | 1  0  1  0 |   | (rTemp0,iTemp0) |   | (rTemp0 + rTemp2, iTemp0 + iTemp2) |

    // Result = | 0  1  0 -j | * | (rTemp1,iTemp1) | = | (rTemp1 + iTemp3, iTemp1 - rTemp3) |

    //          | 1  0 -1  0 |   | (rTemp2,iTemp2) |   | (rTemp0 - rTemp2, iTemp0 - iTemp2) |

    //          | 0  1  0  j |   | (rTemp3,iTemp3) |   | (rTemp1 - iTemp3, iTemp1 + rTemp3) |

    __forceinline void ButterflyDIT4_4 (__inout XVECTORREF r0,

                                        __inout XVECTORREF r1,

                                        __inout XVECTORREF r2,

                                        __inout XVECTORREF r3,

                                        __inout XVECTORREF i0,

                                        __inout XVECTORREF i1,

                                        __inout XVECTORREF i2,

                                        __inout XVECTORREF i3,

                                        __in_ecount(uStride*4) const XVECTOR* __restrict pUnityTableReal,

                                        __in_ecount(uStride*4) const XVECTOR* __restrict pUnityTableImaginary,

                                        const UINT32 uStride, const BOOL fLast)

    {

        DSPASSERT(pUnityTableReal != NULL);

        DSPASSERT(pUnityTableImaginary != NULL);

        DSPASSERT((UINT_PTR)pUnityTableReal % 16 == 0);

        DSPASSERT((UINT_PTR)pUnityTableImaginary % 16 == 0);

        DSPASSERT(ISPOWEROF2(uStride));

        XVECTOR rTemp0, rTemp1, rTemp2, rTemp3, rTemp4, rTemp5, rTemp6, rTemp7;

        XVECTOR iTemp0, iTemp1, iTemp2, iTemp3, iTemp4, iTemp5, iTemp6, iTemp7;

        // calculating Temp

        rTemp0 = _mm_add_ps(r0, r2);          iTemp0 = _mm_add_ps(i0, i2);

        rTemp2 = _mm_add_ps(r1, r3);          iTemp2 = _mm_add_ps(i1, i3);

        rTemp1 = _mm_sub_ps(r0, r2);          iTemp1 = _mm_sub_ps(i0, i2);

        rTemp3 = _mm_sub_ps(r1, r3);          iTemp3 = _mm_sub_ps(i1, i3);

        rTemp4 = _mm_add_ps(rTemp0, rTemp2);  iTemp4 = _mm_add_ps(iTemp0, iTemp2);

        rTemp5 = _mm_add_ps(rTemp1, iTemp3);  iTemp5 = _mm_sub_ps(iTemp1, rTemp3);

        rTemp6 = _mm_sub_ps(rTemp0, rTemp2);  iTemp6 = _mm_sub_ps(iTemp0, iTemp2);

        rTemp7 = _mm_sub_ps(rTemp1, iTemp3);  iTemp7 = _mm_add_ps(iTemp1, rTemp3);

        // calculating Result

        // vmulComplex(rTemp0, iTemp0, rTemp0, iTemp0, pUnityTableReal[0], pUnityTableImaginary[0]); // first one is always trivial

        vmulComplex(rTemp5, iTemp5, pUnityTableReal[uStride], pUnityTableImaginary[uStride]);

        vmulComplex(rTemp6, iTemp6, pUnityTableReal[uStride*2], pUnityTableImaginary[uStride*2]);

        vmulComplex(rTemp7, iTemp7, pUnityTableReal[uStride*3], pUnityTableImaginary[uStride*3]);

        if (fLast) {

            ButterflyDIT4_1(rTemp4, iTemp4);

            ButterflyDIT4_1(rTemp5, iTemp5);

            ButterflyDIT4_1(rTemp6, iTemp6);

            ButterflyDIT4_1(rTemp7, iTemp7);

        }

        r0 = rTemp4;    i0 = iTemp4;

        r1 = rTemp5;    i1 = iTemp5;

        r2 = rTemp6;    i2 = iTemp6;

        r3 = rTemp7;    i3 = iTemp7;

    }

//--------------<F-U-N-C-T-I-O-N-S>-----------------------------------------//

      ////

      // DESCRIPTION:

      //  4-sample FFT.

      //

      // PARAMETERS:

      //  pReal      - [inout] real components, must have at least uCount elements

      //  pImaginary - [inout] imaginary components, must have at least uCount elements

      //  uCount     - [in]    number of FFT iterations

      //

      // RETURN VALUE:

      //  void

      ////

    __forceinline void FFT4 (__inout_ecount(uCount) XVECTOR* __restrict pReal, __inout_ecount(uCount) XVECTOR* __restrict pImaginary, const UINT32 uCount=1)

    {

        DSPASSERT(pReal != NULL);

        DSPASSERT(pImaginary != NULL);

        DSPASSERT((UINT_PTR)pReal % 16 == 0);

        DSPASSERT((UINT_PTR)pImaginary % 16 == 0);

        DSPASSERT(ISPOWEROF2(uCount));

        for (UINT32 uIndex=0; uIndex<uCount; ++uIndex) {

            ButterflyDIT4_1(pReal[uIndex], pImaginary[uIndex]);

        }

    }

      ////

      // DESCRIPTION:

      //  8-sample FFT.

      //

      // PARAMETERS:

      //  pReal      - [inout] real components, must have at least uCount*2 elements

      //  pImaginary - [inout] imaginary components, must have at least uCount*2 elements

      //  uCount     - [in]    number of FFT iterations

      //

      // RETURN VALUE:

      //  void

      ////

    __forceinline void FFT8 (__inout_ecount(uCount*2) XVECTOR* __restrict pReal, __inout_ecount(uCount*2) XVECTOR* __restrict pImaginary, const UINT32 uCount=1)

    {

        DSPASSERT(pReal != NULL);

        DSPASSERT(pImaginary != NULL);

        DSPASSERT((UINT_PTR)pReal % 16 == 0);

        DSPASSERT((UINT_PTR)pImaginary % 16 == 0);

        DSPASSERT(ISPOWEROF2(uCount));

        static XVECTOR wr1 = {  1.0f,  0.70710677f,  0.0f, -0.70710677f };

        static XVECTOR wi1 = {  0.0f, -0.70710677f, -1.0f, -0.70710677f };

        static XVECTOR wr2 = { -1.0f, -0.70710677f,  0.0f,  0.70710677f };

        static XVECTOR wi2 = {  0.0f,  0.70710677f,  1.0f,  0.70710677f };

        for (UINT32 uIndex=0; uIndex<uCount; ++uIndex) {

            XVECTOR* __restrict pR = pReal      + uIndex*2;

            XVECTOR* __restrict pI = pImaginary + uIndex*2;

            XVECTOR oddsR  = _mm_shuffle_ps(pR[0], pR[1], _MM_SHUFFLE(3, 1, 3, 1));

            XVECTOR evensR = _mm_shuffle_ps(pR[0], pR[1], _MM_SHUFFLE(2, 0, 2, 0));

            XVECTOR oddsI  = _mm_shuffle_ps(pI[0], pI[1], _MM_SHUFFLE(3, 1, 3, 1));

            XVECTOR evensI = _mm_shuffle_ps(pI[0], pI[1], _MM_SHUFFLE(2, 0, 2, 0));

            ButterflyDIT4_1(oddsR, oddsI);

            ButterflyDIT4_1(evensR, evensI);

            XVECTOR r, i;

            vmulComplex(r, i, oddsR, oddsI, wr1, wi1);

            pR[0] = _mm_add_ps(evensR, r);

            pI[0] = _mm_add_ps(evensI, i);

            vmulComplex(r, i, oddsR, oddsI, wr2, wi2);

            pR[1] = _mm_add_ps(evensR, r);

            pI[1] = _mm_add_ps(evensI, i);

        }

    }

      ////

      // DESCRIPTION:

      //  16-sample FFT.

      //

      // PARAMETERS:

      //  pReal      - [inout] real components, must have at least uCount*4 elements

      //  pImaginary - [inout] imaginary components, must have at least uCount*4 elements

      //  uCount     - [in]    number of FFT iterations

      //

      // RETURN VALUE:

      //  void

      ////

    __forceinline void FFT16 (__inout_ecount(uCount*4) XVECTOR* __restrict pReal, __inout_ecount(uCount*4) XVECTOR* __restrict pImaginary, const UINT32 uCount=1)

    {

        DSPASSERT(pReal != NULL);

        DSPASSERT(pImaginary != NULL);

        DSPASSERT((UINT_PTR)pReal % 16 == 0);

        DSPASSERT((UINT_PTR)pImaginary % 16 == 0);

        DSPASSERT(ISPOWEROF2(uCount));

        XVECTOR aUnityTableReal[4]      = { 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 0.92387950f, 0.70710677f, 0.38268343f, 1.0f, 0.70710677f, -4.3711388e-008f, -0.70710677f, 1.0f, 0.38268343f, -0.70710677f, -0.92387950f };

        XVECTOR aUnityTableImaginary[4] = { -0.0f, -0.0f, -0.0f, -0.0f, -0.0f, -0.38268343f, -0.70710677f, -0.92387950f, -0.0f, -0.70710677f, -1.0f, -0.70710677f, -0.0f, -0.92387950f, -0.70710677f, 0.38268343f };

        for (UINT32 uIndex=0; uIndex<uCount; ++uIndex) {

            ButterflyDIT4_4(pReal[uIndex*4],

                            pReal[uIndex*4 + 1],

                            pReal[uIndex*4 + 2],

                            pReal[uIndex*4 + 3],

                            pImaginary[uIndex*4],

                            pImaginary[uIndex*4 + 1],

                            pImaginary[uIndex*4 + 2],

                            pImaginary[uIndex*4 + 3],

                            aUnityTableReal,

                            aUnityTableImaginary,

                            1, TRUE);

        }

    }

      ////

      // DESCRIPTION:

      //  2^N-sample FFT.

      //

      // REMARKS:

      //  For FFTs length 16 and below, call FFT16(), FFT8(), or FFT4().

      //

      // PARAMETERS:

      //  pReal       - [inout] real components, must have at least (uLength*uCount)/4 elements

      //  pImaginary  - [inout] imaginary components, must have at least (uLength*uCount)/4 elements

      //  pUnityTable - [in]    unity table, must have at least uLength*uCount elements, see FFTInitializeUnityTable()

      //  uLength     - [in]    FFT length in samples, must be a power of 2 > 16

      //  uCount      - [in]    number of FFT iterations

      //

      // RETURN VALUE:

      //  void

      ////

    inline void FFT (__inout_ecount((uLength*uCount)/4) XVECTOR* __restrict pReal, __inout_ecount((uLength*uCount)/4) XVECTOR* __restrict pImaginary, __in_ecount(uLength*uCount) const XVECTOR* __restrict pUnityTable, const UINT32 uLength, const UINT32 uCount=1)

    {

        DSPASSERT(pReal != NULL);

        DSPASSERT(pImaginary != NULL);

        DSPASSERT(pUnityTable != NULL);

        DSPASSERT((UINT_PTR)pReal % 16 == 0);

        DSPASSERT((UINT_PTR)pImaginary % 16 == 0);

        DSPASSERT((UINT_PTR)pUnityTable % 16 == 0);

        DSPASSERT(uLength > 16);

        DSPASSERT(ISPOWEROF2(uLength));

        DSPASSERT(ISPOWEROF2(uCount));

        const XVECTOR* __restrict pUnityTableReal      = pUnityTable;

        const XVECTOR* __restrict pUnityTableImaginary = pUnityTable + (uLength>>2);

        const UINT32 uTotal         = uCount * uLength;

        const UINT32 uTotal_vectors = uTotal >> 2;

        const UINT32 uStage_vectors = uLength >> 2;

        const UINT32 uStride        = uStage_vectors >> 2; // stride between butterfly elements

        const UINT32 uSkip          = uStage_vectors - uStride;

        for (UINT32 uIndex=0; uIndex<(uTotal_vectors>>2); ++uIndex) {

            UINT32 n = (uIndex/uStride) * (uStride + uSkip) + (uIndex % uStride);

            ButterflyDIT4_4(pReal[n],

                            pReal[n + uStride],

                            pReal[n + uStride*2],

                            pReal[n + uStride*3],

                            pImaginary[n ],

                            pImaginary[n + uStride],

                            pImaginary[n + uStride*2],

                            pImaginary[n + uStride*3],

                            pUnityTableReal      + n % uStage_vectors,

                            pUnityTableImaginary + n % uStage_vectors,

                            uStride, FALSE);

        }

        if (uLength > 16*4) {

            FFT(pReal, pImaginary, pUnityTable+(uLength>>1), uLength>>2, uCount*4);

        } else if (uLength == 16*4) {

            FFT16(pReal, pImaginary, uCount*4);

        } else if (uLength == 8*4) {

            FFT8(pReal, pImaginary, uCount*4);

        } else if (uLength == 4*4) {

            FFT4(pReal, pImaginary, uCount*4);

        }

    }

//--------------------------------------------------------------------------//

  ////

  // DESCRIPTION:

  //  Initializes unity roots lookup table used by FFT functions.

  //  Once initialized, the table need not be initialized again unless a

  //  different FFT length is desired.

  //

  // REMARKS:

  //  The unity tables of FFT length 16 and below are hard coded into the

  //  respective FFT functions and so need not be initialized.

  //

  // PARAMETERS:

  //  pUnityTable - [out] unity table, receives unity roots lookup table, must have at least uLength elements

  //  uLength     - [in]  FFT length in samples, must be a power of 2 > 16

  //

  // RETURN VALUE:

  //  void

  ////

inline void FFTInitializeUnityTable (__out_ecount(uLength) XVECTOR* __restrict pUnityTable, UINT32 uLength)

{

    DSPASSERT(pUnityTable != NULL);

    DSPASSERT(uLength > 16);

    DSPASSERT(ISPOWEROF2(uLength));

    FLOAT32* __restrict pfUnityTable = (FLOAT32* __restrict)pUnityTable;

    // initialize unity table for recursive FFT lengths: uLength, uLength/4, uLength/16... > 16

    do {

        FLOAT32 flStep = 6.283185307f / uLength; // 2PI / FFT length

        uLength >>= 2;

        // pUnityTable[0 to uLength*4-1] contains real components for current FFT length

        // pUnityTable[uLength*4 to uLength*8-1] contains imaginary components for current FFT length

        for (UINT32 i=0; i<4; ++i) {

            for (UINT32 j=0; j<uLength; ++j) {

                UINT32 uIndex = (i*uLength) + j;

                pfUnityTable[uIndex]             = cosf(FLOAT32(i)*FLOAT32(j)*flStep);  // real component

                pfUnityTable[uIndex + uLength*4] = -sinf(FLOAT32(i)*FLOAT32(j)*flStep); // imaginary component

            }

        }

        pfUnityTable += uLength*8;

    } while (uLength > 16);

}

  ////

  // DESCRIPTION:

  //  The FFT functions generate output in bit reversed order.

  //  Use this function to re-arrange them into order of increasing frequency.

  //

  // PARAMETERS:

  //  pOutput     - [out] output buffer, receives samples in order of increasing frequency, must have at least (1<<uLog2Length) elements

  //  pInput      - [in]  input buffer, samples in bit reversed order as generated by FFT functions, must have at least (1<<uLog2Length) elements

  //  uLog2Length - [in]  LOG (base 2) of FFT length in samples, must be > 0

  //

  // RETURN VALUE:

  //  void

  ////

inline void FFTUnswizzle (__out_ecount(1<<uLog2Length) FLOAT32* __restrict pOutput, __in_ecount(1<<uLog2Length) const FLOAT32* __restrict pInput, const UINT32 uLog2Length)

{

    DSPASSERT(pOutput != NULL);

    DSPASSERT(pInput != NULL);

    DSPASSERT(uLog2Length > 0);

    const UINT32 uLength = UINT32(1 << uLog2Length);

    if ((uLog2Length & 0x1) == 0) {

        // even powers of two

        for (UINT32 uIndex=0; uIndex<uLength; ++uIndex) {

            UINT32 n = uIndex;

            n = ( (n & 0xcccccccc) >> 2 )  | ( (n & 0x33333333) << 2 );

            n = ( (n & 0xf0f0f0f0) >> 4 )  | ( (n & 0x0f0f0f0f) << 4 );

            n = ( (n & 0xff00ff00) >> 8 )  | ( (n & 0x00ff00ff) << 8 );

            n = ( (n & 0xffff0000) >> 16 ) | ( (n & 0x0000ffff) << 16 );

            n >>= (32 - uLog2Length);

            pOutput[n] = pInput[uIndex];

        }

    } else {

        // odd powers of two

        for (UINT32 uIndex=0; uIndex<uLength; ++uIndex) {

            UINT32 n = (uIndex>>3);

            n = ( (n & 0xcccccccc) >> 2 )  | ( (n & 0x33333333) << 2 );

            n = ( (n & 0xf0f0f0f0) >> 4 )  | ( (n & 0x0f0f0f0f) << 4 );

            n = ( (n & 0xff00ff00) >> 8 )  | ( (n & 0x00ff00ff) << 8 );

            n = ( (n & 0xffff0000) >> 16 ) | ( (n & 0x0000ffff) << 16 );

            n >>= (32 - (uLog2Length-3));

            n |= ((uIndex & 0x7) << (uLog2Length - 3));

            pOutput[n] = pInput[uIndex];

        }

    }

}

  ////

  // DESCRIPTION:

  //  Convert complex components to polar form.

  //

  // PARAMETERS:

  //  pOutput         - [out] output buffer, receives samples in polar form, must have at least uLength/4 elements

  //  pInputReal      - [in]  input buffer (real components), must have at least uLength/4 elements

  //  pInputImaginary - [in]  input buffer (imaginary components), must have at least uLength/4 elements

  //  uLength         - [in]  FFT length in samples, must be a power of 2 >= 4

  //

  // RETURN VALUE:

  //  void

  ////

inline void FFTPolar (__out_ecount(uLength/4) XVECTOR* __restrict pOutput, __in_ecount(uLength/4) const XVECTOR* __restrict pInputReal, __in_ecount(uLength/4) const XVECTOR* __restrict pInputImaginary, const UINT32 uLength)

{

    DSPASSERT(pOutput != NULL);

    DSPASSERT(pInputReal != NULL);

    DSPASSERT(pInputImaginary != NULL);

    DSPASSERT(uLength >= 4);

    DSPASSERT(ISPOWEROF2(uLength));

    FLOAT32 flOneOverLength = 1.0f / uLength;

    // result = sqrtf((real/uLength)^2 + (imaginary/uLength)^2) * 2

        XVECTOR vOneOverLength = _mm_set_ps1(flOneOverLength);

        for (UINT32 uIndex=0; uIndex<(uLength>>2); ++uIndex) {

            XVECTOR vReal      = _mm_mul_ps(pInputReal[uIndex], vOneOverLength);

            XVECTOR vImaginary = _mm_mul_ps(pInputImaginary[uIndex], vOneOverLength);

            XVECTOR vRR        = _mm_mul_ps(vReal, vReal);

            XVECTOR vII        = _mm_mul_ps(vImaginary, vImaginary);

            XVECTOR vRRplusII  = _mm_add_ps(vRR, vII);

            XVECTOR vTotal  = _mm_sqrt_ps(vRRplusII);

            pOutput[uIndex] = _mm_add_ps(vTotal, vTotal);

        }

}

#pragma warning(pop)

}; // namespace XDSP

//---------------------------------<-EOF->----------------------------------//