/*
 * BSD Licence:
 * Copyright (c) 2001, 2002 Ben Houston [ ben@exocortex.org ]
 * Exocortex Technologies [ www.exocortex.org ]
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, 
 * this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright 
 * notice, this list of conditions and the following disclaimer in the 
 * documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the <ORGANIZATION> nor the names of its contributors
 * may be used to endorse or promote products derived from this software
 * without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 */

using System;
using System.Diagnostics;
using System.IO;
using System.Text;

//using Exocortex.Imaging;

namespace Exocortex.DSP
{

    // Comments? Questions? Bugs? Tell Ben Houston at ben@exocortex.org
    // Version: May 4, 2002

    /// <summary>
    /// <p>Static functions for doing various Fourier Operations.</p>
    /// </summary>
    public class Fourier
    {

        //======================================================================================

        private Fourier()
        {
        }

        //======================================================================================

        static private void Swap(ref float a, ref float b)
        {
            float temp = a;
            a = b;
            b = temp;
        }
        static private void Swap(ref double a, ref double b)
        {
            double temp = a;
            a = b;
            b = temp;
        }
        static private void Swap(ref ComplexF a, ref ComplexF b)
        {
            ComplexF temp = a;
            a = b;
            b = temp;
        }
        static private void Swap(ref Complex a, ref Complex b)
        {
            Complex temp = a;
            a = b;
            b = temp;
        }

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

        private const int cMaxLength = 16777216;
        private const int cMinLength = 1;

        private const int cMaxBits = 24;
        private const int cMinBits = 0;


        static private bool IsPowerOf2(int x)
        {
            return (x & (x - 1)) == 0;
            //return	( x == Pow2( Log2( x ) ) );
        }
        static private int Pow2(int exponent)
        {
            if (exponent >= 0 && exponent < 31)
            {
                return 1 << exponent;
            }
            return 0;
        }
        static private int Log2(int x)
        {
            if (x <= 65536)
            {
                if (x <= 256)
                {
                    if (x <= 16)
                    {
                        if (x <= 4)
                        {
                            if (x <= 2)
                            {
                                if (x <= 1)
                                {
                                    return 0;
                                }
                                return 1;
                            }
                            return 2;
                        }
                        if (x <= 8)
                            return 3;
                        return 4;
                    }
                    if (x <= 64)
                    {
                        if (x <= 32)
                            return 5;
                        return 6;
                    }
                    if (x <= 128)
                        return 7;
                    return 8;
                }
                if (x <= 4096)
                {
                    if (x <= 1024)
                    {
                        if (x <= 512)
                            return 9;
                        return 10;
                    }
                    if (x <= 2048)
                        return 11;
                    return 12;
                }
                if (x <= 16384)
                {
                    if (x <= 8192)
                        return 13;
                    return 14;
                }
                if (x <= 32768)
                    return 15;
                return 16;
            }
            if (x <= 16777216)
            {
                if (x <= 1048576)
                {
                    if (x <= 262144)
                    {
                        if (x <= 131072)
                            return 17;
                        return 18;
                    }
                    if (x <= 524288)
                        return 19;
                    return 20;
                }
                if (x <= 4194304)
                {
                    if (x <= 2097152)
                        return 21;
                    return 22;
                }
                if (x <= 8388608)
                    return 23;
                return 24;
            }
            if (x <= 268435456)
            {
                if (x <= 67108864)
                {
                    if (x <= 33554432)
                        return 25;
                    return 26;
                }
                if (x <= 134217728)
                    return 27;
                return 28;
            }
            if (x <= 1073741824)
            {
                if (x <= 536870912)
                    return 29;
                return 30;
            }
            //	since int is unsigned it can never be higher than 2,147,483,647
            //	if( x <= 2147483648 )
            //		return	31;	
            //	return	32;	
            return 31;
        }

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

        static private int ReverseBits(int index, int numberOfBits)
        {
            Debug.Assert(numberOfBits >= cMinBits);
            Debug.Assert(numberOfBits <= cMaxBits);

            int reversedIndex = 0;
            for (int i = 0; i < numberOfBits; i++)
            {
                reversedIndex = (reversedIndex << 1) | (index & 1);
                index = (index >> 1);
            }
            return reversedIndex;
        }

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

        static private int[][] _reversedBits = new int[cMaxBits][];
        static private int[] GetReversedBits(int numberOfBits)
        {
            Debug.Assert(numberOfBits >= cMinBits);
            Debug.Assert(numberOfBits <= cMaxBits);
            if (_reversedBits[numberOfBits - 1] == null)
            {
                int maxBits = Fourier.Pow2(numberOfBits);
                int[] reversedBits = new int[maxBits];
                for (int i = 0; i < maxBits; i++)
                {
                    int oldBits = i;
                    int newBits = 0;
                    for (int j = 0; j < numberOfBits; j++)
                    {
                        newBits = (newBits << 1) | (oldBits & 1);
                        oldBits = (oldBits >> 1);
                    }
                    reversedBits[i] = newBits;
                }
                _reversedBits[numberOfBits - 1] = reversedBits;
            }
            return _reversedBits[numberOfBits - 1];
        }

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

        static private void ReorderArray(float[] data)
        {
            Debug.Assert(data != null);

            int length = data.Length / 2;

            Debug.Assert(Fourier.IsPowerOf2(length) == true);
            Debug.Assert(length >= cMinLength);
            Debug.Assert(length <= cMaxLength);

            int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    Fourier.Swap(ref data[(i << 1)], ref data[(swap << 1)]);
                    Fourier.Swap(ref data[(i << 1) + 1], ref data[(swap << 1) + 1]);
                }
            }
        }

        static private void ReorderArray(double[] data)
        {
            Debug.Assert(data != null);

            int length = data.Length / 2;

            Debug.Assert(Fourier.IsPowerOf2(length) == true);
            Debug.Assert(length >= cMinLength);
            Debug.Assert(length <= cMaxLength);

            int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    Fourier.Swap(ref data[i << 1], ref data[swap << 1]);
                    Fourier.Swap(ref data[i << 1 + 1], ref data[swap << 1 + 1]);
                }
            }
        }

        static private void ReorderArray(Complex[] data)
        {
            Debug.Assert(data != null);

            int length = data.Length;

            Debug.Assert(Fourier.IsPowerOf2(length) == true);
            Debug.Assert(length >= cMinLength);
            Debug.Assert(length <= cMaxLength);

            int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    Complex temp = data[i];
                    data[i] = data[swap];
                    data[swap] = temp;
                }
            }
        }

        static private void ReorderArray(ComplexF[] data)
        {
            Debug.Assert(data != null);

            int length = data.Length;

            Debug.Assert(Fourier.IsPowerOf2(length) == true);
            Debug.Assert(length >= cMinLength);
            Debug.Assert(length <= cMaxLength);

            int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    ComplexF temp = data[i];
                    data[i] = data[swap];
                    data[swap] = temp;
                }
            }
        }

        //======================================================================================

        private static int[][] _reverseBits = null;

        private static int _ReverseBits(int bits, int n)
        {
            int bitsReversed = 0;
            for (int i = 0; i < n; i++)
            {
                bitsReversed = (bitsReversed << 1) | (bits & 1);
                bits = (bits >> 1);
            }
            return bitsReversed;
        }

        private static void InitializeReverseBits(int levels)
        {
            _reverseBits = new int[levels + 1][];
            for (int j = 0; j < (levels + 1); j++)
            {
                int count = (int)Math.Pow(2, j);
                _reverseBits[j] = new int[count];
                for (int i = 0; i < count; i++)
                {
                    _reverseBits[j][i] = _ReverseBits(i, j);
                }
            }
        }

        private static int _lookupTabletLength = -1;
        private static double[,][] _uRLookup = null;
        private static double[,][] _uILookup = null;
        private static float[,][] _uRLookupF = null;
        private static float[,][] _uILookupF = null;

        private static void SyncLookupTableLength(int length)
        {
            //			Debug.Assert( length < 1024*10 );
            Debug.Assert(length >= 0);
            if (length > _lookupTabletLength)
            {
                int level = (int)Math.Ceiling(Math.Log(length, 2));
                Fourier.InitializeReverseBits(level);
                Fourier.InitializeComplexRotations(level);
                //_cFFTData	= new Complex[ Math2.CeilingBase( length, 2 ) ];
                //_cFFTDataF	= new ComplexF[ Math2.CeilingBase( length, 2 ) ];
                _lookupTabletLength = length;
            }
        }

        private static int GetLookupTableLength()
        {
            return _lookupTabletLength;
        }

        private static void ClearLookupTables()
        {
            _uRLookup = null;
            _uILookup = null;
            _uRLookupF = null;
            _uILookupF = null;
            _lookupTabletLength = -1;
        }

        private static void InitializeComplexRotations(int levels)
        {
            int ln = levels;
            //_wRLookup = new float[ levels + 1, 2 ];
            //_wILookup = new float[ levels + 1, 2 ];

            _uRLookup = new double[levels + 1, 2][];
            _uILookup = new double[levels + 1, 2][];

            _uRLookupF = new float[levels + 1, 2][];
            _uILookupF = new float[levels + 1, 2][];

            int N = 1;
            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                //float scale = (float)( 1 / Math.Sqrt( 1 << ln ) );

                // positive sign ( i.e. [M,0] )
                {
                    double uR = 1;
                    double uI = 0;
                    double angle = Math.PI / M * 1;
                    double wR = (double)Math.Cos(angle);
                    double wI = (double)Math.Sin(angle);

                    _uRLookup[level, 0] = new double[M];
                    _uILookup[level, 0] = new double[M];
                    _uRLookupF[level, 0] = new float[M];
                    _uILookupF[level, 0] = new float[M];

                    for (int j = 0; j < M; j++)
                    {
                        _uRLookupF[level, 0][j] = (float)(_uRLookup[level, 0][j] = uR);
                        _uILookupF[level, 0][j] = (float)(_uILookup[level, 0][j] = uI);
                        double uwI = uR * wI + uI * wR;
                        uR = uR * wR - uI * wI;
                        uI = uwI;
                    }
                }
                {


                    // negative sign ( i.e. [M,1] )
                    double uR = 1;
                    double uI = 0;
                    double angle = Math.PI / M * -1;
                    double wR = (double)Math.Cos(angle);
                    double wI = (double)Math.Sin(angle);

                    _uRLookup[level, 1] = new double[M];
                    _uILookup[level, 1] = new double[M];
                    _uRLookupF[level, 1] = new float[M];
                    _uILookupF[level, 1] = new float[M];

                    for (int j = 0; j < M; j++)
                    {
                        _uRLookupF[level, 1][j] = (float)(_uRLookup[level, 1][j] = uR);
                        _uILookupF[level, 1][j] = (float)(_uILookup[level, 1][j] = uI);
                        double uwI = uR * wI + uI * wR;
                        uR = uR * wR - uI * wI;
                        uI = uwI;
                    }
                }

            }
        }

        //======================================================================================
        //======================================================================================

        static private bool _bufferFLocked = false;
        static private float[] _bufferF = new float[0];

        static private void LockBufferF(int length, ref float[] buffer)
        {
            Debug.Assert(_bufferFLocked == false);
            _bufferFLocked = true;
            if (length >= _bufferF.Length)
            {
                _bufferF = new float[length];
            }
            buffer = _bufferF;
        }
        static private void UnlockBufferF(ref float[] buffer)
        {
            Debug.Assert(_bufferF == buffer);
            Debug.Assert(_bufferFLocked == true);
            _bufferFLocked = false;
            buffer = null;
        }

        private static void LinearFFT(float[] data, int start, int inc, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(start >= 0);
            Debug.Assert(inc >= 1);
            Debug.Assert(length >= 1);
            Debug.Assert((start + inc * (length - 1)) * 2 < data.Length);

            // copy to buffer
            float[] buffer = null;
            LockBufferF(length * 2, ref buffer);
            int j = start;
            for (int i = 0; i < length * 2; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferF(ref buffer);
        }

        private static void LinearFFT_Quick(float[] data, int start, int inc, int length, FourierDirection direction)
        {
            /*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) * 2 < data.Length );*/

            // copy to buffer
            float[] buffer = null;
            LockBufferF(length * 2, ref buffer);
            int j = start;
            for (int i = 0; i < length * 2; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT_Quick(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferF(ref buffer);
        }

        //======================================================================================
        //======================================================================================

        static private bool _bufferCFLocked = false;
        static private ComplexF[] _bufferCF = new ComplexF[0];

        static private void LockBufferCF(int length, ref ComplexF[] buffer)
        {
            Debug.Assert(length >= 0);
            Debug.Assert(_bufferCFLocked == false);

            _bufferCFLocked = true;
            if (length != _bufferCF.Length)
            {
                _bufferCF = new ComplexF[length];
            }
            buffer = _bufferCF;
        }
        static private void UnlockBufferCF(ref ComplexF[] buffer)
        {
            Debug.Assert(_bufferCF == buffer);
            Debug.Assert(_bufferCFLocked == true);

            _bufferCFLocked = false;
            buffer = null;
        }

        private static void LinearFFT(ComplexF[] data, int start, int inc, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(start >= 0);
            Debug.Assert(inc >= 1);
            Debug.Assert(length >= 1);
            Debug.Assert((start + inc * (length - 1)) < data.Length);

            // copy to buffer
            ComplexF[] buffer = null;
            LockBufferCF(length, ref buffer);
            int j = start;
            for (int i = 0; i < length; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferCF(ref buffer);
        }

        private static void LinearFFT_Quick(ComplexF[] data, int start, int inc, int length, FourierDirection direction)
        {
            /*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );	*/

            // copy to buffer
            ComplexF[] buffer = null;
            LockBufferCF(length, ref buffer);
            int j = start;
            for (int i = 0; i < length; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferCF(ref buffer);
        }

        //======================================================================================
        //======================================================================================

        static private bool _bufferCLocked = false;
        static private Complex[] _bufferC = new Complex[0];

        static private void LockBufferC(int length, ref Complex[] buffer)
        {
            Debug.Assert(length >= 0);
            Debug.Assert(_bufferCLocked == false);

            _bufferCLocked = true;
            if (length >= _bufferC.Length)
            {
                _bufferC = new Complex[length];
            }
            buffer = _bufferC;
        }
        static private void UnlockBufferC(ref Complex[] buffer)
        {
            Debug.Assert(_bufferC == buffer);
            Debug.Assert(_bufferCLocked == true);

            _bufferCLocked = false;
            buffer = null;
        }

        private static void LinearFFT(Complex[] data, int start, int inc, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(start >= 0);
            Debug.Assert(inc >= 1);
            Debug.Assert(length >= 1);
            Debug.Assert((start + inc * (length - 1)) < data.Length);

            // copy to buffer
            Complex[] buffer = null;
            LockBufferC(length, ref buffer);
            int j = start;
            for (int i = 0; i < length; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferC(ref buffer);
        }

        private static void LinearFFT_Quick(Complex[] data, int start, int inc, int length, FourierDirection direction)
        {
            /*Debug.Assert( data != null );
			Debug.Assert( start >= 0 );
			Debug.Assert( inc >= 1 );
			Debug.Assert( length >= 1 );
			Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );*/

            // copy to buffer
            Complex[] buffer = null;
            LockBufferC(length, ref buffer);
            int j = start;
            for (int i = 0; i < length; i++)
            {
                buffer[i] = data[j];
                j += inc;
            }

            FFT_Quick(buffer, length, direction);

            // copy from buffer
            j = start;
            for (int i = 0; i < length; i++)
            {
                data[j] = buffer[i];
                j += inc;
            }
            UnlockBufferC(ref buffer);
        }

        //======================================================================================
        //======================================================================================

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT(float[] data, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(data.Length >= length * 2);
            Debug.Assert(Fourier.IsPowerOf2(length) == true);

            Fourier.SyncLookupTableLength(length);

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                float[] uRLookup = _uRLookupF[level, signIndex];
                float[] uILookup = _uILookupF[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    float uR = uRLookup[j];
                    float uI = uILookup[j];

                    for (int evenT = j; evenT < length; evenT += N)
                    {
                        int even = evenT << 1;
                        int odd = (evenT + M) << 1;

                        float r = data[odd];
                        float i = data[odd + 1];

                        float odduR = r * uR - i * uI;
                        float odduI = r * uI + i * uR;

                        r = data[even];
                        i = data[even + 1];

                        data[even] = r + odduR;
                        data[even + 1] = i + odduI;

                        data[odd] = r - odduR;
                        data[odd + 1] = i - odduI;
                    }
                }
            }
        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT_Quick(float[] data, int length, FourierDirection direction)
        {
            Debug.Assert(data != null);
            Debug.Assert(data.Length >= length * 2);
            Debug.Assert(Fourier.IsPowerOf2(length) == true);

            Fourier.SyncLookupTableLength(length);

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                float[] uRLookup = _uRLookupF[level, signIndex];
                float[] uILookup = _uILookupF[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    float uR = uRLookup[j];
                    float uI = uILookup[j];

                    for (int evenT = j; evenT < length; evenT += N)
                    {
                        int even = evenT << 1;
                        int odd = (evenT + M) << 1;

                        float r = data[odd];
                        float i = data[odd + 1];

                        float odduR = r * uR - i * uI;
                        float odduI = r * uI + i * uR;

                        r = data[even];
                        i = data[even + 1];

                        data[even] = r + odduR;
                        data[even + 1] = i + odduI;

                        data[odd] = r - odduR;
                        data[odd + 1] = i - odduI;
                    }
                }
            }
        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT(ComplexF[] data, int length, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < length)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
            }
            if (Fourier.IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

            Fourier.SyncLookupTableLength(length);

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;

            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                float[] uRLookup = _uRLookupF[level, signIndex];
                float[] uILookup = _uILookupF[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    float uR = uRLookup[j];
                    float uI = uILookup[j];

                    for (int even = j; even < length; even += N)
                    {
                        int odd = even + M;

                        float r = data[odd].Re;
                        float i = data[odd].Im;

                        float odduR = r * uR - i * uI;
                        float odduI = r * uI + i * uR;

                        r = data[even].Re;
                        i = data[even].Im;

                        data[even].Re = r + odduR;
                        data[even].Im = i + odduI;

                        data[odd].Re = r - odduR;
                        data[odd].Im = i - odduI;
                    }
                }
            }

        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT_Quick(ComplexF[] data, int length, FourierDirection direction)
        {
            /*if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < length ) {
				throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );
			}
			if( Fourier.IsPowerOf2( length ) == false ) {
				throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );
			}

			Fourier.SyncLookupTableLength( length );*/

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;

            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                float[] uRLookup = _uRLookupF[level, signIndex];
                float[] uILookup = _uILookupF[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    float uR = uRLookup[j];
                    float uI = uILookup[j];

                    for (int even = j; even < length; even += N)
                    {
                        int odd = even + M;

                        float r = data[odd].Re;
                        float i = data[odd].Im;

                        float odduR = r * uR - i * uI;
                        float odduI = r * uI + i * uR;

                        r = data[even].Re;
                        i = data[even].Im;

                        data[even].Re = r + odduR;
                        data[even].Im = i + odduI;

                        data[odd].Re = r - odduR;
                        data[odd].Im = i - odduI;
                    }
                }
            }

        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="direction"></param>
        public static void FFT(ComplexF[] data, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            Fourier.FFT(data, data.Length, direction);
        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT(Complex[] data, int length, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < length)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
            }
            if (Fourier.IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

            Fourier.SyncLookupTableLength(length);

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;

            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                double[] uRLookup = _uRLookup[level, signIndex];
                double[] uILookup = _uILookup[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    double uR = uRLookup[j];
                    double uI = uILookup[j];

                    for (int even = j; even < length; even += N)
                    {
                        int odd = even + M;

                        double r = data[odd].Re;
                        double i = data[odd].Im;

                        double odduR = r * uR - i * uI;
                        double odduI = r * uI + i * uR;

                        r = data[even].Re;
                        i = data[even].Im;

                        data[even].Re = r + odduR;
                        data[even].Im = i + odduI;

                        data[odd].Re = r - odduR;
                        data[odd].Im = i - odduI;
                    }
                }
            }

        }

        /// <summary>
        /// Compute a 1D fast Fourier transform of a dataset of complex numbers.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void FFT_Quick(Complex[] data, int length, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < length)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
            }
            if (Fourier.IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

            Fourier.SyncLookupTableLength(length);

            int ln = Fourier.Log2(length);

            // reorder array
            Fourier.ReorderArray(data);

            // successive doubling
            int N = 1;
            int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;

            for (int level = 1; level <= ln; level++)
            {
                int M = N;
                N <<= 1;

                double[] uRLookup = _uRLookup[level, signIndex];
                double[] uILookup = _uILookup[level, signIndex];

                for (int j = 0; j < M; j++)
                {
                    double uR = uRLookup[j];
                    double uI = uILookup[j];

                    for (int even = j; even < length; even += N)
                    {
                        int odd = even + M;

                        double r = data[odd].Re;
                        double i = data[odd].Im;

                        double odduR = r * uR - i * uI;
                        double odduI = r * uI + i * uR;

                        r = data[even].Re;
                        i = data[even].Im;

                        data[even].Re = r + odduR;
                        data[even].Im = i + odduI;

                        data[odd].Re = r - odduR;
                        data[odd].Im = i - odduI;
                    }
                }
            }

        }

        /// <summary>
        /// Compute a 1D real-symmetric fast fourier transform.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="direction"></param>
        public static void RFFT(float[] data, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            Fourier.RFFT(data, data.Length, direction);
        }

        /// <summary>
        /// Compute a 1D real-symmetric fast fourier transform.
        /// </summary>
        /// <param name="data"></param>
        /// <param name="length"></param>
        /// <param name="direction"></param>
        public static void RFFT(float[] data, int length, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < length)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
            }
            if (Fourier.IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

            float c1 = 0.5f, c2;
            float theta = (float)Math.PI / (length / 2);

            if (direction == FourierDirection.Forward)
            {
                c2 = -0.5f;
                FFT(data, length / 2, direction);
            }
            else
            {
                c2 = 0.5f;
                theta = -theta;
            }

            float wtemp = (float)Math.Sin(0.5 * theta);
            float wpr = -2 * wtemp * wtemp;
            float wpi = (float)Math.Sin(theta);
            float wr = 1 + wpr;
            float wi = wpi;

            // do / undo packing
            for (int i = 1; i < length / 4; i++)
            {
                int a = 2 * i;
                int b = length - 2 * i;
                float h1r = c1 * (data[a] + data[b]);
                float h1i = c1 * (data[a + 1] - data[b + 1]);
                float h2r = -c2 * (data[a + 1] + data[b + 1]);
                float h2i = c2 * (data[a] - data[b]);
                data[a] = h1r + wr * h2r - wi * h2i;
                data[a + 1] = h1i + wr * h2i + wi * h2r;
                data[b] = h1r - wr * h2r + wi * h2i;
                data[b + 1] = -h1i + wr * h2i + wi * h2r;
                wr = (wtemp = wr) * wpr - wi * wpi + wr;
                wi = wi * wpr + wtemp * wpi + wi;
            }

            if (direction == FourierDirection.Forward)
            {
                float hir = data[0];
                data[0] = hir + data[1];
                data[1] = hir - data[1];
            }
            else
            {
                float hir = data[0];
                data[0] = c1 * (hir + data[1]);
                data[1] = c1 * (hir - data[1]);
                Fourier.FFT(data, length / 2, direction);
            }
        }

        /// <summary>
        /// Compute a 2D fast fourier transform on a data set of complex numbers (represented as pairs of floats)
        /// </summary>
        /// <param name="data"></param>
        /// <param name="xLength"></param>
        /// <param name="yLength"></param>
        /// <param name="direction"></param>
        public static void FFT2(float[] data, int xLength, int yLength, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < xLength * yLength * 2)
            {
                throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * 2' parameter");
            }
            if (Fourier.IsPowerOf2(xLength) == false)
            {
                throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(yLength) == false)
            {
                throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
            }

            int xInc = 1;
            int yInc = xLength;

            if (xLength > 1)
            {
                Fourier.SyncLookupTableLength(xLength);
                for (int y = 0; y < yLength; y++)
                {
                    int xStart = y * yInc;
                    Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
                }
            }

            if (yLength > 1)
            {
                Fourier.SyncLookupTableLength(yLength);
                for (int x = 0; x < xLength; x++)
                {
                    int yStart = x * xInc;
                    Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
                }
            }
        }

        /// <summary>
        /// Compute a 2D fast fourier transform on a data set of complex numbers
        /// </summary>
        /// <param name="data"></param>
        /// <param name="xLength"></param>
        /// <param name="yLength"></param>
        /// <param name="direction"></param>
        public static void FFT2(ComplexF[] data, int xLength, int yLength, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < xLength * yLength)
            {
                throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter");
            }
            if (Fourier.IsPowerOf2(xLength) == false)
            {
                throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(yLength) == false)
            {
                throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
            }

            int xInc = 1;
            int yInc = xLength;

            if (xLength > 1)
            {
                Fourier.SyncLookupTableLength(xLength);
                for (int y = 0; y < yLength; y++)
                {
                    int xStart = y * yInc;
                    Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
                }
            }

            if (yLength > 1)
            {
                Fourier.SyncLookupTableLength(yLength);
                for (int x = 0; x < xLength; x++)
                {
                    int yStart = x * xInc;
                    Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
                }
            }
        }

        /// <summary>
        /// Compute a 2D fast fourier transform on a data set of complex numbers
        /// </summary>
        /// <param name="data"></param>
        /// <param name="xLength"></param>
        /// <param name="yLength"></param>
        /// <param name="direction"></param>
        public static void FFT2(Complex[] data, int xLength, int yLength, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < xLength * yLength)
            {
                throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter");
            }
            if (Fourier.IsPowerOf2(xLength) == false)
            {
                throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(yLength) == false)
            {
                throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
            }

            int xInc = 1;
            int yInc = xLength;

            if (xLength > 1)
            {
                Fourier.SyncLookupTableLength(xLength);
                for (int y = 0; y < yLength; y++)
                {
                    int xStart = y * yInc;
                    Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
                }
            }

            if (yLength > 1)
            {
                Fourier.SyncLookupTableLength(yLength);
                for (int x = 0; x < xLength; x++)
                {
                    int yStart = x * xInc;
                    Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
                }
            }
        }

        /// <summary>
        /// Compute a 3D fast fourier transform on a data set of complex numbers
        /// </summary>
        /// <param name="data"></param>
        /// <param name="xLength"></param>
        /// <param name="yLength"></param>
        /// <param name="zLength"></param>
        /// <param name="direction"></param>
        public static void FFT3(ComplexF[] data, int xLength, int yLength, int zLength, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < xLength * yLength * zLength)
            {
                throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter");
            }
            if (Fourier.IsPowerOf2(xLength) == false)
            {
                throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(yLength) == false)
            {
                throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(zLength) == false)
            {
                throw new ArgumentOutOfRangeException("zLength", zLength, "must be a power of 2");
            }

            int xInc = 1;
            int yInc = xLength;
            int zInc = xLength * yLength;

            if (xLength > 1)
            {
                Fourier.SyncLookupTableLength(xLength);
                for (int z = 0; z < zLength; z++)
                {
                    for (int y = 0; y < yLength; y++)
                    {
                        int xStart = y * yInc + z * zInc;
                        Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
                    }
                }
            }

            if (yLength > 1)
            {
                Fourier.SyncLookupTableLength(yLength);
                for (int z = 0; z < zLength; z++)
                {
                    for (int x = 0; x < xLength; x++)
                    {
                        int yStart = z * zInc + x * xInc;
                        Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
                    }
                }
            }

            if (zLength > 1)
            {
                Fourier.SyncLookupTableLength(zLength);
                for (int y = 0; y < yLength; y++)
                {
                    for (int x = 0; x < xLength; x++)
                    {
                        int zStart = y * yInc + x * xInc;
                        Fourier.LinearFFT_Quick(data, zStart, zInc, zLength, direction);
                    }
                }
            }
        }

        /// <summary>
        /// Compute a 3D fast fourier transform on a data set of complex numbers
        /// </summary>
        /// <param name="data"></param>
        /// <param name="xLength"></param>
        /// <param name="yLength"></param>
        /// <param name="zLength"></param>
        /// <param name="direction"></param>
        public static void FFT3(Complex[] data, int xLength, int yLength, int zLength, FourierDirection direction)
        {
            if (data == null)
            {
                throw new ArgumentNullException("data");
            }
            if (data.Length < xLength * yLength * zLength)
            {
                throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter");
            }
            if (Fourier.IsPowerOf2(xLength) == false)
            {
                throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(yLength) == false)
            {
                throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
            }
            if (Fourier.IsPowerOf2(zLength) == false)
            {
                throw new ArgumentOutOfRangeException("zLength", zLength, "must be a power of 2");
            }

            int xInc = 1;
            int yInc = xLength;
            int zInc = xLength * yLength;

            if (xLength > 1)
            {
                Fourier.SyncLookupTableLength(xLength);
                for (int z = 0; z < zLength; z++)
                {
                    for (int y = 0; y < yLength; y++)
                    {
                        int xStart = y * yInc + z * zInc;
                        Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
                    }
                }
            }

            if (yLength > 1)
            {
                Fourier.SyncLookupTableLength(yLength);
                for (int z = 0; z < zLength; z++)
                {
                    for (int x = 0; x < xLength; x++)
                    {
                        int yStart = z * zInc + x * xInc;
                        Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
                    }
                }
            }

            if (zLength > 1)
            {
                Fourier.SyncLookupTableLength(zLength);
                for (int y = 0; y < yLength; y++)
                {
                    for (int x = 0; x < xLength; x++)
                    {
                        int zStart = y * yInc + x * xInc;
                        Fourier.LinearFFT_Quick(data, zStart, zInc, zLength, direction);
                    }
                }
            }
        }

    }
}