using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace PQS.Server.MathProcessor.Utils
{

    internal enum FourierDirection
    {
        /// <summary>
        /// Forward direction.  Usually in reference to moving from temporal
        /// representation to frequency representation
        /// </summary>
        Forward = 1,
        /// <summary>
        /// Backward direction. Usually in reference to moving from frequency
        /// representation to temporal representation
        /// </summary>
        Backward = -1,
    }

    internal sealed class FourierUtils
    {

        private const int MAXLENGTH = 8192; // was 4096
        private const int MINLENGTH = 1;
        private const int MAXBITS = 13; //was 12
        private const int MINBITS = 0;

        private int[][] _reversedBits = new int[MAXBITS][];
        private int[][] _reverseBits = null;
        private int _lookupTabletLength = -1;
        private double[,][] _uRLookup = null;
        private double[,][] _uILookup = null;
        private float[,][] _uRLookupF = null;
        private float[,][] _uILookupF = null;
        //		private bool _bufferFLocked	    = false;
        private float[] _bufferF = new float[0];

        public FourierUtils(int length)
        {
            this.SyncLookupTableLength(length);
        }

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

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

        static public bool IsPowerOf2(int x)
        {
            return (x & (x - 1)) == 0;
        }

        private static int Pow2(int exponent)
        {
            if (exponent >= 0 && exponent < 31)
            {
                return 1 << exponent;
            }
            return 0;
        }

        public static 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;
        }

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

        private int[] GetReversedBits(int numberOfBits)
        {
            if (_reversedBits[numberOfBits - 1] == null)
            {
                int maxBits = FourierUtils.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];
        }

        private void ReorderArray(float[] data)
        {
            int length = data.Length / 2;
            int[] reversedBits = this.GetReversedBits(FourierUtils.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    FourierUtils.Swap(ref data[(i << 1)], ref data[(swap << 1)]);
                    FourierUtils.Swap(ref data[(i << 1) + 1], ref data[(swap << 1) + 1]);
                }
            }
        }

        private void ReorderArray(double[] data)
        {
            int length = data.Length / 2;
            int[] reversedBits = this.GetReversedBits(FourierUtils.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    FourierUtils.Swap(ref data[i << 1], ref data[swap << 1]);
                    FourierUtils.Swap(ref data[i << 1 + 1], ref data[swap << 1 + 1]);
                }
            }
        }

        private void InitializeReverseBits(int levels)
        {
            _reverseBits = new int[levels + 1][];
            int LevelsPlusOne = levels + 1;
            for (int j = 0; j < LevelsPlusOne; 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 void SyncLookupTableLength(int length)
        {
            if (length > _lookupTabletLength)
            {
                int level = (int)Math.Ceiling(Math.Log(length, 2));
                InitializeReverseBits(level);
                InitializeComplexRotations(level);
                //_cFFTData	= new Complex[ Math2.CeilingBase( length, 2 ) ];
                //_cFFTDataF	= new ComplexF[ Math2.CeilingBase( length, 2 ) ];
                _lookupTabletLength = length;
            }
        }

        //private int	GetLookupTableLength() {
        //    return	_lookupTabletLength;
        //}

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

        private 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 = (double)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 = (double)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;
                    }
                }

            }
        }

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

        public void FFT(float[] data, int length, FourierDirection direction)
        {

            // this.SyncLookupTableLength( length );

            // int ln = this.Log2( length );

            this.ReorderArray(data);

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

            int FourierLogResult = FourierUtils.Log2(length);
            for (int level = 1; level <= FourierLogResult; 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;
                    }
                }
            }
        }

        public void RFFT(float[] data, int length, FourierDirection direction)
        {

#if DEBUG

            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 (IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

#endif

            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]);
                FFT(data, length / 2, direction);
            }
        }

    }

    internal unsafe sealed class uFourier
    {

        private const int MAXLENGTH = 8192; // was 4096
        private const int MINLENGTH = 1;
        private const int MAXBITS = 13; //was 12
        private const int MINBITS = 0;

        private int[][] _reversedBits = new int[MAXBITS][];
        private int[][] _reverseBits = null;
        private int _lookupTabletLength = -1;
        private double[,][] _uRLookup = null;
        private double[,][] _uILookup = null;
        private float[,][] _uRLookupF = null;
        private float[,][] _uILookupF = null;
        //		private bool _bufferFLocked	    = false;
        private float[] _bufferF = new float[0];

        public uFourier(int length)
        {
            this.SyncLookupTableLength(length);
        }

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

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

        static public bool IsPowerOf2(int x)
        {
            return (x & (x - 1)) == 0;
        }

        private static int Pow2(int exponent)
        {
            if (exponent >= 0 && exponent < 31)
            {
                return 1 << exponent;
            }
            return 0;
        }

        public static 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;
        }

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

        private int[] GetReversedBits(int numberOfBits)
        {
            if (_reversedBits[numberOfBits - 1] == null)
            {
                int maxBits = uFourier.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];
        }

        private void ReorderArray(float* data, int Length)
        {
            fixed (int* ptrReversedBits = this.GetReversedBits(FourierUtils.Log2(Length)))
            {
                for (int i = 0; i < Length; i++)
                {
                    int swap = ptrReversedBits[i];
                    if (swap > i)
                    {
                        uFourier.Swap(ref data[(i << 1)], ref data[(swap << 1)]);
                        uFourier.Swap(ref data[(i << 1) + 1], ref data[(swap << 1) + 1]);
                    }
                }
            }
        }

        private void ReorderArray(double* data, int Length)
        {
            int length = Length / 2;
            int[] reversedBits = this.GetReversedBits(FourierUtils.Log2(length));
            for (int i = 0; i < length; i++)
            {
                int swap = reversedBits[i];
                if (swap > i)
                {
                    uFourier.Swap(ref data[i << 1], ref data[swap << 1]);
                    uFourier.Swap(ref data[i << 1 + 1], ref data[swap << 1 + 1]);
                }
            }
        }

        private void InitializeReverseBits(int levels)
        {
            _reverseBits = new int[levels + 1][];
            int LevelsPlusOne = levels + 1;
            for (int j = 0; j < LevelsPlusOne; 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 void SyncLookupTableLength(int length)
        {
            if (length > _lookupTabletLength)
            {
                int level = (int)Math.Ceiling(Math.Log(length, 2));
                InitializeReverseBits(level);
                InitializeComplexRotations(level);
                //_cFFTData	= new Complex[ Math2.CeilingBase( length, 2 ) ];
                //_cFFTDataF	= new ComplexF[ Math2.CeilingBase( length, 2 ) ];
                _lookupTabletLength = length;
            }
        }

        //private int	GetLookupTableLength() {
        //    return	_lookupTabletLength;
        //}

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

        private 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 = (double)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 = (double)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;
                    }
                }

            }
        }

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

        public void FFT(float* data, int length, FourierDirection direction)
        {

            this.ReorderArray(data, length);

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

            int FourierLogResult = FourierUtils.Log2(length);
            for (int level = 1; level <= FourierLogResult; level++)
            {
                int M = N;
                N <<= 1;

                fixed (float* uRLookup = _uRLookupF[level, signIndex], 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;
                        }
                    }
                }
            }
        }

        public void RFFT(float* data, int length, FourierDirection direction)
        {

#if DEBUG

            if (data == null)
            {
                throw new ArgumentNullException("data");
            }

            if (IsPowerOf2(length) == false)
            {
                throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
            }

#endif

            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]);
                FFT(data, length / 2, direction);
            }
        }

    }

    public static class FourierWrapper
    { // sealed for optimization

        #region Members

        private static FourierUtils _Fourier16384;
        private static FourierUtils _Fourier8192;
        private static FourierUtils _Fourier4096;
        private static FourierUtils _Fourier2048;
        private static FourierUtils _Fourier1024;
        private static FourierUtils _Fourier512;
        private static FourierUtils _Fourier256;
        private static FourierUtils _Fourier128;
        private static FourierUtils _Fourier64;
        private static FourierUtils _Fourier32;
        private static FourierUtils _Fourier16;
        private static FourierUtils _Fourier8;
        private static FourierUtils _Fourier4;
        private static FourierUtils _Fourier2;

        private volatile static object _SyncLock = new object();

        #endregion

        ///<summary>
        ///Preforms IFFT on a flat 1D signal. The length of the flat signal should be mod2 and the values
        ///should look like |Re|Im|Re|Im|Re|Im|Re|Im|Re|Im....
        ///</summary>
        ///<param name="Sig">The Harmonic Vector.</param>
        ///<returns>A real signal of the wave form.</returns>	
        //public void IFFTReal(ref float[] Sig) {

        //    int FFTNumber = Sig.Length >> 1; // Both Image and real

        //    switch (FFTNumber) {                
        //        case (32):
        //            if (_Fourier32 == null) _Fourier32 = new Fourier(32);
        //            _Fourier32.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (64):
        //            if (_Fourier64 == null) _Fourier64 = new Fourier(64);
        //            _Fourier64.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (128):
        //            if (_Fourier128 == null) _Fourier128 = new Fourier(128);
        //            _Fourier128.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (256):
        //            if (_Fourier256 == null) _Fourier256 = new Fourier(256);
        //            _Fourier256.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (512):
        //            if (_Fourier512 == null) _Fourier512 = new Fourier(512);
        //            _Fourier512.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (1024):
        //            if (_Fourier1024 == null) _Fourier1024 = new Fourier(1024);
        //            _Fourier1024.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (2048):
        //            if (_Fourier2048 == null) _Fourier2048 = new Fourier(2048);
        //            _Fourier2048.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (4096):
        //            if (_Fourier4096 == null) _Fourier4096 = new Fourier(4096);
        //            _Fourier4096.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;
        //        case (8192):
        //            if (_Fourier8192 == null) _Fourier8192 = new Fourier(8192);
        //            _Fourier8192.FFT(Sig, FFTNumber, FourierDirection.Backward);
        //            break;                    
        //        default:
        //            throw new ArgumentException("FFT " + FFTNumber + " not supported");
        //    }

        //    float[] FFTSigReal = new float[FFTNumber];
        //    // Take Real part only.
        //    int FFTSigRealLength = FFTSigReal.Length;
        //    for (int i = 0; i < FFTSigRealLength; i++) FFTSigReal[i] = Sig[i << 1];            

        //    Sig = FFTSigReal;
        //    return;
        //}

        /// <summary>
        /// Preforms an IFFT operation on a real signal (only real values are returned).
        /// </summary>
        /// <param name="Sig">Value array</param>
        /// <remarks>Note that the lock on the creation of the FFT objects does not prevent doing the same creation multiple times but it is much
        /// better than doing the lock on the 'if'</remarks>
        public static void IFFTReal(ref float[] Sig)
        {

            int FFTNumber = Sig.Length;

            switch (FFTNumber)
            {
                case (2):
                    if (_Fourier2 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2 == null)
                            {
                                _Fourier2 = new FourierUtils(2);
                            }
                        }
                    }
                    _Fourier2.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (4):
                    if (_Fourier4 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4 == null)
                            {
                                _Fourier4 = new FourierUtils(4);
                            }
                        }
                    }
                    _Fourier4.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (8):
                    if (_Fourier8 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8 == null)
                            {
                                _Fourier8 = new FourierUtils(8);
                            }
                        }
                    }
                    _Fourier8.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (16):
                    if (_Fourier16 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16 == null)
                            {
                                _Fourier16 = new FourierUtils(16);
                            }
                        }
                    }
                    _Fourier16.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (32):
                    if (_Fourier32 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier32 == null)
                            {
                                _Fourier32 = new FourierUtils(32);
                            }
                        }
                    }
                    _Fourier32.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (64):
                    if (_Fourier64 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier64 == null)
                            {
                                _Fourier64 = new FourierUtils(64);
                            }
                        }
                    }
                    _Fourier64.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (128):
                    if (_Fourier128 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier128 == null)
                            {
                                _Fourier128 = new FourierUtils(128);
                            }
                        }
                    }
                    _Fourier128.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (256):
                    if (_Fourier256 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier256 == null)
                            {
                                _Fourier256 = new FourierUtils(256);
                            }
                        }
                    }
                    _Fourier256.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (512):
                    if (_Fourier512 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier512 == null)
                            {
                                _Fourier512 = new FourierUtils(512);
                            }
                        }
                    }
                    _Fourier512.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (1024):
                    if (_Fourier1024 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier1024 == null)
                            {
                                _Fourier1024 = new FourierUtils(1024);
                            }
                        }
                    }
                    _Fourier1024.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (2048):
                    if (_Fourier2048 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2048 == null)
                            {
                                _Fourier2048 = new FourierUtils(2048);
                            }
                        }
                    }
                    _Fourier2048.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (4096):
                    if (_Fourier4096 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4096 == null)
                            {
                                _Fourier4096 = new FourierUtils(4096);
                            }
                        }
                    }
                    _Fourier4096.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (8192):
                    if (_Fourier8192 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8192 == null)
                            {
                                _Fourier8192 = new FourierUtils(8192);
                            }
                        }
                    }
                    _Fourier8192.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (16384):
                    if (_Fourier16384 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16384 == null)
                            {
                                _Fourier16384 = new FourierUtils(16384);
                            }
                        }
                    }
                    _Fourier16384.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                default:
                    throw new ArgumentException("FFT " + FFTNumber + " not supported");
            }

            #region Fix response

            float FixFactor = 2.0f / (float)Sig.Length;
            for (int i = 0; i < Sig.Length; i++)
            {
                Sig[i] *= FixFactor;
            }

            #endregion

            Array.Reverse(Sig, 1, Sig.Length - 1);

            return;
        }

        /// <summary>
        /// Preforms an FFT operation on a real signal (only real values are returned).
        /// </summary>
        /// <param name="Sig">Value array</param>
        /// <remarks>Note that the lock on the creation of the FFT objects does not prevent doing the same creation multiple times but it is much
        /// better than doing the lock on the 'if'</remarks>
        public static void FFTReal(ref float[] Sig)
        {

            int FFTNumber = Sig.Length;

            switch (FFTNumber)
            {
                case (2):
                    if (_Fourier2 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2 == null)
                            {
                                _Fourier2 = new FourierUtils(2);
                            }
                        }
                    }
                    _Fourier2.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (4):
                    if (_Fourier4 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4 == null)
                            {
                                _Fourier4 = new FourierUtils(4);
                            }
                        }
                    }
                    _Fourier4.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (8):
                    if (_Fourier8 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8 == null)
                            {
                                _Fourier8 = new FourierUtils(8);
                            }
                        }
                    }
                    _Fourier8.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (16):
                    if (_Fourier16 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16 == null)
                            {
                                _Fourier16 = new FourierUtils(16);
                            }
                        }
                    }
                    _Fourier16.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (32):
                    if (_Fourier32 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier32 == null)
                            {
                                _Fourier32 = new FourierUtils(32);
                            }
                        }
                    }
                    _Fourier32.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (64):
                    if (_Fourier64 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier64 == null)
                            {
                                _Fourier64 = new FourierUtils(64);
                            }
                        }
                    }
                    _Fourier64.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (128):
                    if (_Fourier128 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier128 == null)
                            {
                                _Fourier128 = new FourierUtils(128);
                            }
                        }
                    }
                    _Fourier128.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (256):
                    if (_Fourier256 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier256 == null)
                            {
                                _Fourier256 = new FourierUtils(256);
                            }
                        }
                    }
                    _Fourier256.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (512):
                    if (_Fourier512 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier512 == null)
                            {
                                _Fourier512 = new FourierUtils(512);
                            }
                        }
                    }
                    _Fourier512.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (1024):
                    if (_Fourier1024 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier1024 == null)
                            {
                                _Fourier1024 = new FourierUtils(1024);
                            }
                        }
                    }
                    _Fourier1024.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (2048):
                    if (_Fourier2048 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2048 == null)
                            {
                                _Fourier2048 = new FourierUtils(2048);
                            }
                        }
                    }
                    _Fourier2048.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (4096):
                    if (_Fourier4096 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4096 == null)
                            {
                                _Fourier4096 = new FourierUtils(4096);
                            }
                        }
                    }
                    _Fourier4096.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (8192):
                    if (_Fourier8192 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8192 == null)
                            {
                                _Fourier8192 = new FourierUtils(8192);
                            }
                        }
                    }
                    _Fourier8192.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (16384):
                    if (_Fourier16384 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16384 == null)
                            {
                                _Fourier16384 = new FourierUtils(16384);
                            }
                        }
                    }
                    _Fourier16384.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                default:
                    throw new ArgumentException("FFT " + FFTNumber + " not supported");
            }
            return;
        }

        public static void Clear()
        {
            lock (_SyncLock)
            {
                _Fourier16384 = null;
                _Fourier8192 = null;
                _Fourier4096 = null;
                _Fourier2048 = null;
                _Fourier1024 = null;
                _Fourier512 = null;
                _Fourier256 = null;
                _Fourier128 = null;
                _Fourier64 = null;
                _Fourier32 = null;
            }
        }
    }

    public unsafe static class uFourierWrapper
    { // sealed for optimization

        #region Members

        private static uFourier _Fourier16384;
        private static uFourier _Fourier8192;
        private static uFourier _Fourier4096;
        private static uFourier _Fourier2048;
        private static uFourier _Fourier1024;
        private static uFourier _Fourier512;
        private static uFourier _Fourier256;
        private static uFourier _Fourier128;
        private static uFourier _Fourier64;
        private static uFourier _Fourier32;
        private static uFourier _Fourier16;
        private static uFourier _Fourier8;
        private static uFourier _Fourier4;
        private static uFourier _Fourier2;

        private volatile static object _SyncLock = new object();

        #endregion

        /// <summary>
        /// Preforms an IFFT operation on a real signal (only real values are returned).
        /// </summary>
        /// <param name="Sig">Value array</param>
        /// <remarks>Note that the lock on the creation of the FFT objects does not prevent doing the same creation multiple times but it is much
        /// better than doing the lock on the 'if'</remarks>
        public static void IFFTReal(float* Sig, int Size)
        {

            int FFTNumber = Size;

            switch (FFTNumber)
            {
                case (2):
                    if (_Fourier2 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2 == null)
                            {
                                _Fourier2 = new uFourier(2);
                            }
                        }
                    }
                    _Fourier2.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (4):
                    if (_Fourier4 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4 == null)
                            {
                                _Fourier4 = new uFourier(4);
                            }
                        }
                    }
                    _Fourier4.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (8):
                    if (_Fourier8 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8 == null)
                            {
                                _Fourier8 = new uFourier(8);
                            }
                        }
                    }
                    _Fourier8.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (16):
                    if (_Fourier16 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16 == null)
                            {
                                _Fourier16 = new uFourier(16);
                            }
                        }
                    }
                    _Fourier16.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (32):
                    if (_Fourier32 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier32 == null)
                            {
                                _Fourier32 = new uFourier(32);
                            }
                        }
                    }
                    _Fourier32.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (64):
                    if (_Fourier64 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier64 == null)
                            {
                                _Fourier64 = new uFourier(64);
                            }
                        }
                    }
                    _Fourier64.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (128):
                    if (_Fourier128 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier128 == null)
                            {
                                _Fourier128 = new uFourier(128);
                            }
                        }
                    }
                    _Fourier128.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (256):
                    if (_Fourier256 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier256 == null)
                            {
                                _Fourier256 = new uFourier(256);
                            }
                        }
                    }
                    _Fourier256.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (512):
                    if (_Fourier512 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier512 == null)
                            {
                                _Fourier512 = new uFourier(512);
                            }
                        }
                    }
                    _Fourier512.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (1024):
                    if (_Fourier1024 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier1024 == null)
                            {
                                _Fourier1024 = new uFourier(1024);
                            }
                        }
                    }
                    _Fourier1024.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (2048):
                    if (_Fourier2048 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2048 == null)
                            {
                                _Fourier2048 = new uFourier(2048);
                            }
                        }
                    }
                    _Fourier2048.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (4096):
                    if (_Fourier4096 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4096 == null)
                            {
                                _Fourier4096 = new uFourier(4096);
                            }
                        }
                    }
                    _Fourier4096.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (8192):
                    if (_Fourier8192 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8192 == null)
                            {
                                _Fourier8192 = new uFourier(8192);
                            }
                        }
                    }
                    _Fourier8192.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                case (16384):
                    if (_Fourier16384 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16384 == null)
                            {
                                _Fourier16384 = new uFourier(16384);
                            }
                        }
                    }
                    _Fourier16384.RFFT(Sig, FFTNumber, FourierDirection.Backward);
                    break;
                default:
                    throw new ArgumentException("FFT " + FFTNumber + " not supported");
            }

            #region Fix response

            float FixFactor = 2.0f / (float)Size;
            for (int i = 0; i < Size; i++)
            {
                Sig[i] *= FixFactor;
            }

            #endregion

            return;
        }

        /// <summary>
        /// Preforms an FFT operation on a real signal (only real values are returned).
        /// </summary>
        /// <param name="Sig">Value array</param>
        /// <remarks>Note that the lock on the creation of the FFT objects does not prevent doing the same creation multiple times but it is much
        /// better than doing the lock on the 'if'</remarks>
        public static void FFTReal(float* Sig, int Size)
        {

            int FFTNumber = Size;

            switch (FFTNumber)
            {
                case (2):
                    if (_Fourier2 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2 == null)
                            {
                                _Fourier2 = new uFourier(2);
                            }
                        }
                    }
                    _Fourier2.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (4):
                    if (_Fourier4 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4 == null)
                            {
                                _Fourier4 = new uFourier(4);
                            }
                        }
                    }
                    _Fourier4.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (8):
                    if (_Fourier8 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8 == null)
                            {
                                _Fourier8 = new uFourier(8);
                            }
                        }
                    }
                    _Fourier8.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (16):
                    if (_Fourier16 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16 == null)
                            {
                                _Fourier16 = new uFourier(16);
                            }
                        }
                    }
                    _Fourier16.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (32):
                    if (_Fourier32 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier32 == null)
                            {
                                _Fourier32 = new uFourier(32);
                            }
                        }
                    }
                    _Fourier32.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (64):
                    if (_Fourier64 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier64 == null)
                            {
                                _Fourier64 = new uFourier(64);
                            }
                        }
                    }
                    _Fourier64.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (128):
                    if (_Fourier128 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier128 == null)
                            {
                                _Fourier128 = new uFourier(128);
                            }
                        }
                    }
                    _Fourier128.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (256):
                    if (_Fourier256 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier256 == null)
                            {
                                _Fourier256 = new uFourier(256);
                            }
                        }
                    }
                    _Fourier256.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (512):
                    if (_Fourier512 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier512 == null)
                            {
                                _Fourier512 = new uFourier(512);
                            }
                        }
                    }
                    _Fourier512.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (1024):
                    if (_Fourier1024 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier1024 == null)
                            {
                                _Fourier1024 = new uFourier(1024);
                            }
                        }
                    }
                    _Fourier1024.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (2048):
                    if (_Fourier2048 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier2048 == null)
                            {
                                _Fourier2048 = new uFourier(2048);
                            }
                        }
                    }
                    _Fourier2048.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (4096):
                    if (_Fourier4096 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier4096 == null)
                            {
                                _Fourier4096 = new uFourier(4096);
                            }
                        }
                    }
                    _Fourier4096.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (8192):
                    if (_Fourier8192 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier8192 == null)
                            {
                                _Fourier8192 = new uFourier(8192);
                            }
                        }
                    }
                    _Fourier8192.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                case (16384):
                    if (_Fourier16384 == null)
                    {
                        lock (_SyncLock)
                        {
                            if (_Fourier16384 == null)
                            {
                                _Fourier16384 = new uFourier(16384);
                            }
                        }
                    }
                    _Fourier16384.RFFT(Sig, FFTNumber, FourierDirection.Forward);
                    break;
                default:
                    throw new ArgumentException("FFT " + FFTNumber + " not supported");
            }
            return;
        }

        public static void Clear()
        {
            lock (_SyncLock)
            {
                _Fourier16384 = null;
                _Fourier8192 = null;
                _Fourier4096 = null;
                _Fourier2048 = null;
                _Fourier1024 = null;
                _Fourier512 = null;
                _Fourier256 = null;
                _Fourier128 = null;
                _Fourier64 = null;
                _Fourier32 = null;
            }
        }
    }

}
