diff options
Diffstat (limited to 'mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs')
| -rw-r--r-- | mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs | 1346 |
1 files changed, 1346 insertions, 0 deletions
diff --git a/mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs b/mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs new file mode 100644 index 00000000000..87bbd971ec7 --- /dev/null +++ b/mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs @@ -0,0 +1,1346 @@ +/*
+ * 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 = 4096;
+ private const int cMinLength = 1;
+
+ private const int cMaxBits = 12;
+ 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 = (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;
+ }
+ }
+
+ }
+ }
+
+ //======================================================================================
+ //======================================================================================
+
+ 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 );
+ }
+ }
+ }
+ }
+
+ }
+}
|