mcs/class/System.Drawing/Test/DrawingTest/Exocortex.DSP/src/Fourier.cs

   1 /*\r
   2  * BSD Licence:\r
   3  * Copyright (c) 2001, 2002 Ben Houston [ ben@exocortex.org ]\r
   4  * Exocortex Technologies [ www.exocortex.org ]\r
   5  * All rights reserved.\r
   6  *\r
   7  * Redistribution and use in source and binary forms, with or without \r
   8  * modification, are permitted provided that the following conditions are met:\r
   9  *\r
  10  * 1. Redistributions of source code must retain the above copyright notice, \r
  11  * this list of conditions and the following disclaimer.\r
  12  * 2. Redistributions in binary form must reproduce the above copyright \r
  13  * notice, this list of conditions and the following disclaimer in the \r
  14  * documentation and/or other materials provided with the distribution.\r
  15  * 3. Neither the name of the <ORGANIZATION> nor the names of its contributors\r
  16  * may be used to endorse or promote products derived from this software\r
  17  * without specific prior written permission.\r
  18  *\r
  19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"\r
  20  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\r
  21  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\r
  22  * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR\r
  23  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL\r
  24  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR\r
  25  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER\r
  26  * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT \r
  27  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY\r
  28  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH\r
  29  * DAMAGE.\r
  30  */\r
  31 \r
  32 using System;\r
  33 using System.Diagnostics;\r
  34 using System.IO;\r
  35 using System.Text;\r
  36 \r
  37 //using Exocortex.Imaging;\r
  38 \r
  39 namespace Exocortex.DSP {\r
  40 \r
  41         // Comments? Questions? Bugs? Tell Ben Houston at ben@exocortex.org\r
  42         // Version: May 4, 2002\r
  43 \r
  44         /// <summary>\r
  45         /// <p>Static functions for doing various Fourier Operations.</p>\r
  46         /// </summary>\r
  47         public class Fourier {\r
  48                 \r
  49                 //======================================================================================\r
  50 \r
  51                 private Fourier() {\r
  52                 }\r
  53 \r
  54                 //======================================================================================\r
  55 \r
  56                 static private void Swap( ref float a, ref float b ) {\r
  57                         float temp = a;\r
  58                         a = b;\r
  59                         b = temp;\r
  60                 }\r
  61                 static private void Swap( ref double a, ref double b ) {\r
  62                         double temp = a;\r
  63                         a = b;\r
  64                         b = temp;\r
  65                 }\r
  66                 static private void Swap( ref ComplexF a, ref ComplexF b ) {\r
  67                         ComplexF temp = a;\r
  68                         a = b;\r
  69                         b = temp;\r
  70                 }\r
  71                 static private void Swap( ref Complex a, ref Complex b ) {\r
  72                         Complex temp = a;\r
  73                         a = b;\r
  74                         b = temp;\r
  75                 }\r
  76                 \r
  77                 //-------------------------------------------------------------------------------------\r
  78                 \r
  79                 private const int       cMaxLength      = 4096;\r
  80                 private const int       cMinLength      = 1;\r
  81 \r
  82                 private const int       cMaxBits        = 12;\r
  83                 private const int       cMinBits        = 0;\r
  84         \r
  85 \r
  86                 static private bool     IsPowerOf2( int x ) {\r
  87                         return  (x & (x - 1)) == 0;\r
  88                         //return        ( x == Pow2( Log2( x ) ) );\r
  89                 }\r
  90                 static private int      Pow2( int exponent ) {\r
  91                         if( exponent >= 0 && exponent < 31 ) {\r
  92                                 return  1 << exponent;\r
  93                         }\r
  94                         return  0;\r
  95                 }\r
  96                 static private int      Log2( int x ) {\r
  97                         if( x <= 65536 ) {\r
  98                                 if( x <= 256 ) {\r
  99                                         if( x <= 16 ) {\r
 100                                                 if( x <= 4 ) {  \r
 101                                                         if( x <= 2 ) {\r
 102                                                                 if( x <= 1 ) {\r
 103                                                                         return  0;\r
 104                                                                 }\r
 105                                                                 return  1;      \r
 106                                                         }\r
 107                                                         return  2;                              \r
 108                                                 }\r
 109                                                 if( x <= 8 )\r
 110                                                         return  3;                      \r
 111                                                 return  4;                              \r
 112                                         }\r
 113                                         if( x <= 64 ) {\r
 114                                                 if( x <= 32 )\r
 115                                                         return  5;      \r
 116                                                 return  6;                              \r
 117                                         }\r
 118                                         if( x <= 128 )\r
 119                                                 return  7;              \r
 120                                         return  8;                              \r
 121                                 }\r
 122                                 if( x <= 4096 ) {       \r
 123                                         if( x <= 1024 ) {       \r
 124                                                 if( x <= 512 )\r
 125                                                         return  9;                      \r
 126                                                 return  10;                             \r
 127                                         }\r
 128                                         if( x <= 2048 )\r
 129                                                 return  11;                     \r
 130                                         return  12;                             \r
 131                                 }\r
 132                                 if( x <= 16384 ) {\r
 133                                         if( x <= 8192 )\r
 134                                                 return  13;                     \r
 135                                         return  14;                             \r
 136                                 }\r
 137                                 if( x <= 32768 )\r
 138                                         return  15;     \r
 139                                 return  16;     \r
 140                         }\r
 141                         if( x <= 16777216 ) {\r
 142                                 if( x <= 1048576 ) {\r
 143                                         if( x <= 262144 ) {     \r
 144                                                 if( x <= 131072 )\r
 145                                                         return  17;                     \r
 146                                                 return  18;                             \r
 147                                         }\r
 148                                         if( x <= 524288 )\r
 149                                                 return  19;                     \r
 150                                         return  20;                             \r
 151                                 }\r
 152                                 if( x <= 4194304 ) {\r
 153                                         if( x <= 2097152 )\r
 154                                                 return  21;     \r
 155                                         return  22;                             \r
 156                                 }\r
 157                                 if( x <= 8388608 )\r
 158                                         return  23;             \r
 159                                 return  24;                             \r
 160                         }\r
 161                         if( x <= 268435456 ) {  \r
 162                                 if( x <= 67108864 ) {   \r
 163                                         if( x <= 33554432 )\r
 164                                                 return  25;                     \r
 165                                         return  26;                             \r
 166                                 }\r
 167                                 if( x <= 134217728 )\r
 168                                         return  27;                     \r
 169                                 return  28;                             \r
 170                         }\r
 171                         if( x <= 1073741824 ) {\r
 172                                 if( x <= 536870912 )\r
 173                                         return  29;                     \r
 174                                 return  30;                             \r
 175                         }\r
 176                         //      since int is unsigned it can never be higher than 2,147,483,647\r
 177                         //      if( x <= 2147483648 )\r
 178                         //              return  31;     \r
 179                         //      return  32;     \r
 180                         return  31;\r
 181                 }\r
 182 \r
 183                 //-------------------------------------------------------------------------------------\r
 184                 //-------------------------------------------------------------------------------------\r
 185 \r
 186                 static private int      ReverseBits( int index, int numberOfBits ) {\r
 187                         Debug.Assert( numberOfBits >= cMinBits );\r
 188                         Debug.Assert( numberOfBits <= cMaxBits );\r
 189 \r
 190                         int reversedIndex = 0;\r
 191                         for( int i = 0; i < numberOfBits; i ++ ) {\r
 192                                 reversedIndex = ( reversedIndex << 1 ) | ( index & 1 );\r
 193                                 index = ( index >> 1 );\r
 194                         }\r
 195                         return reversedIndex;\r
 196                 }\r
 197 \r
 198                 //-------------------------------------------------------------------------------------\r
 199                 \r
 200                 static private int[][]  _reversedBits   = new int[ cMaxBits ][];\r
 201                 static private int[]            GetReversedBits( int numberOfBits ) {\r
 202                         Debug.Assert( numberOfBits >= cMinBits );\r
 203                         Debug.Assert( numberOfBits <= cMaxBits );\r
 204                         if( _reversedBits[ numberOfBits - 1 ] == null ) {\r
 205                                 int             maxBits = Fourier.Pow2( numberOfBits );\r
 206                                 int[]   reversedBits = new int[ maxBits ];\r
 207                                 for( int i = 0; i < maxBits; i ++ ) {\r
 208                                         int oldBits = i;\r
 209                                         int newBits = 0;\r
 210                                         for( int j = 0; j < numberOfBits; j ++ ) {\r
 211                                                 newBits = ( newBits << 1 ) | ( oldBits & 1 );\r
 212                                                 oldBits = ( oldBits >> 1 );\r
 213                                         }\r
 214                                         reversedBits[ i ] = newBits;\r
 215                                 }\r
 216                                 _reversedBits[ numberOfBits - 1 ] = reversedBits;\r
 217                         }\r
 218                         return  _reversedBits[ numberOfBits - 1 ];\r
 219                 }\r
 220 \r
 221                 //-------------------------------------------------------------------------------------\r
 222 \r
 223                 static private void ReorderArray( float[] data ) {\r
 224                         Debug.Assert( data != null );\r
 225 \r
 226                         int length = data.Length / 2;\r
 227                         \r
 228                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 229                         Debug.Assert( length >= cMinLength );\r
 230                         Debug.Assert( length <= cMaxLength );\r
 231 \r
 232                         int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );\r
 233                         for( int i = 0; i < length; i ++ ) {\r
 234                                 int swap = reversedBits[ i ];\r
 235                                 if( swap > i ) {\r
 236                                         Fourier.Swap( ref data[ (i<<1) ], ref data[ (swap<<1) ] );\r
 237                                         Fourier.Swap( ref data[ (i<<1) + 1 ], ref data[ (swap<<1) + 1 ] );\r
 238                                 }\r
 239                         }\r
 240                 }\r
 241 \r
 242                 static private void ReorderArray( double[] data ) {\r
 243                         Debug.Assert( data != null );\r
 244 \r
 245                         int length = data.Length / 2;\r
 246                         \r
 247                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 248                         Debug.Assert( length >= cMinLength );\r
 249                         Debug.Assert( length <= cMaxLength );\r
 250 \r
 251                         int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );\r
 252                         for( int i = 0; i < length; i ++ ) {\r
 253                                 int swap = reversedBits[ i ];\r
 254                                 if( swap > i ) {\r
 255                                         Fourier.Swap( ref data[ i<<1 ], ref data[ swap<<1 ] );\r
 256                                         Fourier.Swap( ref data[ i<<1 + 1 ], ref data[ swap<<1 + 1 ] );\r
 257                                 }\r
 258                         }\r
 259                 }\r
 260 \r
 261                 static private void ReorderArray( Complex[] data ) {\r
 262                         Debug.Assert( data != null );\r
 263         \r
 264                         int length = data.Length;\r
 265                         \r
 266                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 267                         Debug.Assert( length >= cMinLength );\r
 268                         Debug.Assert( length <= cMaxLength );\r
 269                         \r
 270                         int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );\r
 271                         for( int i = 0; i < length; i ++ ) {\r
 272                                 int swap = reversedBits[ i ];\r
 273                                 if( swap > i ) {\r
 274                                         Complex temp = data[ i ];\r
 275                                         data[ i ] = data[ swap ];\r
 276                                         data[ swap ] = temp;\r
 277                                 }\r
 278                         }\r
 279                 }\r
 280 \r
 281                 static private void ReorderArray( ComplexF[] data ) {\r
 282                         Debug.Assert( data != null );\r
 283 \r
 284                         int length = data.Length;\r
 285                         \r
 286                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 287                         Debug.Assert( length >= cMinLength );\r
 288                         Debug.Assert( length <= cMaxLength );\r
 289 \r
 290                         int[] reversedBits = Fourier.GetReversedBits( Fourier.Log2( length ) );\r
 291                         for( int i = 0; i < length; i ++ ) {\r
 292                                 int swap = reversedBits[ i ];\r
 293                                 if( swap > i ) {\r
 294                                         ComplexF temp = data[ i ];\r
 295                                         data[ i ] = data[ swap ];\r
 296                                         data[ swap ] = temp;\r
 297                                 }\r
 298                         }\r
 299                 }\r
 300 \r
 301                 //======================================================================================\r
 302 \r
 303                 private static int[][]  _reverseBits = null;\r
 304 \r
 305                 private static int      _ReverseBits( int bits, int n ) {\r
 306                         int bitsReversed = 0;\r
 307                         for( int i = 0; i < n; i ++ ) {\r
 308                                 bitsReversed = ( bitsReversed << 1 ) | ( bits & 1 );\r
 309                                 bits = ( bits >> 1 );\r
 310                         }\r
 311                         return bitsReversed;\r
 312                 }\r
 313 \r
 314                 private static void     InitializeReverseBits( int levels ) {\r
 315                         _reverseBits = new int[levels + 1][];\r
 316                         for( int j = 0; j < ( levels + 1 ); j ++ ) {\r
 317                                 int count = (int) Math.Pow( 2, j );\r
 318                                 _reverseBits[j] = new int[ count ];\r
 319                                 for( int i = 0; i < count; i ++ ) {\r
 320                                         _reverseBits[j][i] = _ReverseBits( i, j );\r
 321                                 }\r
 322                         }\r
 323                 }\r
 324 \r
 325                 private static int _lookupTabletLength = -1;\r
 326                 private static double[,][]      _uRLookup       = null;\r
 327                 private static double[,][]      _uILookup       = null;\r
 328                 private static  float[,][]      _uRLookupF      = null;\r
 329                 private static  float[,][]      _uILookupF      = null;\r
 330 \r
 331                 private static void     SyncLookupTableLength( int length ) {\r
 332                         Debug.Assert( length < 1024*10 );\r
 333                         Debug.Assert( length >= 0 );\r
 334                         if( length > _lookupTabletLength ) {\r
 335                                 int level = (int) Math.Ceiling( Math.Log( length, 2 ) );\r
 336                                 Fourier.InitializeReverseBits( level );\r
 337                                 Fourier.InitializeComplexRotations( level );\r
 338                                 //_cFFTData     = new Complex[ Math2.CeilingBase( length, 2 ) ];\r
 339                                 //_cFFTDataF    = new ComplexF[ Math2.CeilingBase( length, 2 ) ];\r
 340                                 _lookupTabletLength = length;\r
 341                         }\r
 342                 }\r
 343 \r
 344                 private static int      GetLookupTableLength() {\r
 345                         return  _lookupTabletLength;\r
 346                 }\r
 347 \r
 348                 private static void     ClearLookupTables() {\r
 349                         _uRLookup       = null;\r
 350                         _uILookup       = null;\r
 351                         _uRLookupF      = null;\r
 352                         _uILookupF      = null;\r
 353                         _lookupTabletLength     = -1;\r
 354                 }\r
 355                 \r
 356                 private static void InitializeComplexRotations( int levels ) {\r
 357                         int ln = levels;\r
 358                         //_wRLookup = new float[ levels + 1, 2 ];\r
 359                         //_wILookup = new float[ levels + 1, 2 ];\r
 360                         \r
 361                         _uRLookup = new double[ levels + 1, 2 ][];\r
 362                         _uILookup = new double[ levels + 1, 2 ][];\r
 363 \r
 364                         _uRLookupF = new float[ levels + 1, 2 ][];\r
 365                         _uILookupF = new float[ levels + 1, 2 ][];\r
 366 \r
 367                         int N = 1;\r
 368                         for( int level = 1; level <= ln; level ++ ) {\r
 369                                 int M = N;\r
 370                                 N <<= 1;\r
 371 \r
 372                                 //float scale = (float)( 1 / Math.Sqrt( 1 << ln ) );\r
 373 \r
 374                                 // positive sign ( i.e. [M,0] )\r
 375                                 {\r
 376                                         double  uR = 1;\r
 377                                         double  uI = 0;\r
 378                                         double  angle = (double) Math.PI / M * 1;\r
 379                                         double  wR = (double) Math.Cos( angle );\r
 380                                         double  wI = (double) Math.Sin( angle );\r
 381 \r
 382                                         _uRLookup[level,0] = new double[ M ];\r
 383                                         _uILookup[level,0] = new double[ M ];\r
 384                                         _uRLookupF[level,0] = new float[ M ];\r
 385                                         _uILookupF[level,0] = new float[ M ];\r
 386 \r
 387                                         for( int j = 0; j < M; j ++ ) {\r
 388                                                 _uRLookupF[level,0][j] = (float)( _uRLookup[level,0][j] = uR );\r
 389                                                 _uILookupF[level,0][j] = (float)( _uILookup[level,0][j] = uI );\r
 390                                                 double  uwI = uR*wI + uI*wR;\r
 391                                                 uR = uR*wR - uI*wI;\r
 392                                                 uI = uwI;\r
 393                                         }\r
 394                                 }\r
 395                                 {\r
 396 \r
 397 \r
 398                                 // negative sign ( i.e. [M,1] )\r
 399                                         double  uR = 1;\r
 400                     double      uI = 0;\r
 401                                         double  angle = (double) Math.PI / M * -1;\r
 402                                         double  wR = (double) Math.Cos( angle );\r
 403                                         double  wI = (double) Math.Sin( angle );\r
 404 \r
 405                                         _uRLookup[level,1] = new double[ M ];\r
 406                                         _uILookup[level,1] = new double[ M ];\r
 407                                         _uRLookupF[level,1] = new float[ M ];\r
 408                                         _uILookupF[level,1] = new float[ M ];\r
 409 \r
 410                                         for( int j = 0; j < M; j ++ ) {\r
 411                                                 _uRLookupF[level,1][j] = (float)( _uRLookup[level,1][j] = uR );\r
 412                                                 _uILookupF[level,1][j] = (float)( _uILookup[level,1][j] = uI );\r
 413                                                 double  uwI = uR*wI + uI*wR;\r
 414                                                 uR = uR*wR - uI*wI;\r
 415                                                 uI = uwI;\r
 416                                         }\r
 417                                 }\r
 418 \r
 419                         }\r
 420                 }\r
 421                 \r
 422                 //======================================================================================\r
 423                 //======================================================================================\r
 424 \r
 425                 static private bool             _bufferFLocked  = false;\r
 426                 static private float[]  _bufferF                = new float[ 0 ];\r
 427 \r
 428                 static private void             LockBufferF( int length, ref float[] buffer ) {\r
 429                         Debug.Assert( _bufferFLocked == false );\r
 430                         _bufferFLocked = true;\r
 431                         if( length >= _bufferF.Length ) {\r
 432                                 _bufferF        = new float[ length ];\r
 433                         }\r
 434                         buffer =        _bufferF;\r
 435                 }\r
 436                 static private void             UnlockBufferF( ref float[] buffer ) {\r
 437                         Debug.Assert( _bufferF == buffer );\r
 438                         Debug.Assert( _bufferFLocked == true );\r
 439                         _bufferFLocked = false;\r
 440                         buffer = null;\r
 441                 }\r
 442 \r
 443                 private static void     LinearFFT( float[] data, int start, int inc, int length, FourierDirection direction ) {\r
 444                         Debug.Assert( data != null );\r
 445                         Debug.Assert( start >= 0 );\r
 446                         Debug.Assert( inc >= 1 );\r
 447                         Debug.Assert( length >= 1 );\r
 448                         Debug.Assert( ( start + inc * ( length - 1 ) ) * 2 < data.Length );\r
 449                         \r
 450                         // copy to buffer\r
 451                         float[] buffer = null;\r
 452                         LockBufferF( length * 2, ref buffer );\r
 453                         int j = start;\r
 454                         for( int i = 0; i < length * 2; i ++ ) {\r
 455                                 buffer[ i ] = data[ j ];\r
 456                                 j += inc;\r
 457                         }\r
 458 \r
 459                         FFT( buffer, length, direction );\r
 460 \r
 461                         // copy from buffer\r
 462                         j = start;\r
 463                         for( int i = 0; i < length; i ++ ) {\r
 464                                 data[ j ] = buffer[ i ];\r
 465                                 j += inc;\r
 466                         }\r
 467                         UnlockBufferF( ref buffer );\r
 468                 }\r
 469 \r
 470                 private static void     LinearFFT_Quick( float[] data, int start, int inc, int length, FourierDirection direction ) {\r
 471                         /*Debug.Assert( data != null );\r
 472                         Debug.Assert( start >= 0 );\r
 473                         Debug.Assert( inc >= 1 );\r
 474                         Debug.Assert( length >= 1 );\r
 475                         Debug.Assert( ( start + inc * ( length - 1 ) ) * 2 < data.Length );*/\r
 476                         \r
 477                         // copy to buffer\r
 478                         float[] buffer = null;\r
 479                         LockBufferF( length * 2, ref buffer );\r
 480                         int j = start;\r
 481                         for( int i = 0; i < length * 2; i ++ ) {\r
 482                                 buffer[ i ] = data[ j ];\r
 483                                 j += inc;\r
 484                         }\r
 485 \r
 486                         FFT_Quick( buffer, length, direction );\r
 487 \r
 488                         // copy from buffer\r
 489                         j = start;\r
 490                         for( int i = 0; i < length; i ++ ) {\r
 491                                 data[ j ] = buffer[ i ];\r
 492                                 j += inc;\r
 493                         }\r
 494                         UnlockBufferF( ref buffer );\r
 495                 }\r
 496 \r
 497                 //======================================================================================\r
 498                 //======================================================================================\r
 499                 \r
 500                 static private bool                     _bufferCFLocked = false;\r
 501                 static private ComplexF[]       _bufferCF               = new ComplexF[ 0 ];\r
 502 \r
 503                 static private void             LockBufferCF( int length, ref ComplexF[] buffer ) {\r
 504                         Debug.Assert( length >= 0 );\r
 505                         Debug.Assert( _bufferCFLocked == false );\r
 506                         \r
 507                         _bufferCFLocked = true;\r
 508                         if( length != _bufferCF.Length ) {\r
 509                                 _bufferCF       = new ComplexF[ length ];\r
 510                         }\r
 511                         buffer =        _bufferCF;\r
 512                 }\r
 513                 static private void             UnlockBufferCF( ref ComplexF[] buffer ) {\r
 514                         Debug.Assert( _bufferCF == buffer );\r
 515                         Debug.Assert( _bufferCFLocked == true );\r
 516                         \r
 517                         _bufferCFLocked = false;\r
 518                         buffer = null;\r
 519                 }\r
 520 \r
 521                 private static void     LinearFFT( ComplexF[] data, int start, int inc, int length, FourierDirection direction ) {\r
 522                         Debug.Assert( data != null );\r
 523                         Debug.Assert( start >= 0 );\r
 524                         Debug.Assert( inc >= 1 );\r
 525                         Debug.Assert( length >= 1 );\r
 526                         Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );\r
 527                         \r
 528                         // copy to buffer\r
 529                         ComplexF[]      buffer = null;\r
 530                         LockBufferCF( length, ref buffer );\r
 531                         int j = start;\r
 532                         for( int i = 0; i < length; i ++ ) {\r
 533                                 buffer[ i ] = data[ j ];\r
 534                                 j += inc;\r
 535                         }\r
 536 \r
 537                         FFT( buffer, length, direction );\r
 538 \r
 539                         // copy from buffer\r
 540                         j = start;\r
 541                         for( int i = 0; i < length; i ++ ) {\r
 542                                 data[ j ] = buffer[ i ];\r
 543                                 j += inc;\r
 544                         }\r
 545                         UnlockBufferCF( ref buffer );\r
 546                 }\r
 547 \r
 548                 private static void     LinearFFT_Quick( ComplexF[] data, int start, int inc, int length, FourierDirection direction ) {\r
 549                         /*Debug.Assert( data != null );\r
 550                         Debug.Assert( start >= 0 );\r
 551                         Debug.Assert( inc >= 1 );\r
 552                         Debug.Assert( length >= 1 );\r
 553                         Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length ); */\r
 554                         \r
 555                         // copy to buffer\r
 556                         ComplexF[]      buffer = null;\r
 557                         LockBufferCF( length, ref buffer );\r
 558                         int j = start;\r
 559                         for( int i = 0; i < length; i ++ ) {\r
 560                                 buffer[ i ] = data[ j ];\r
 561                                 j += inc;\r
 562                         }\r
 563 \r
 564                         FFT( buffer, length, direction );\r
 565 \r
 566                         // copy from buffer\r
 567                         j = start;\r
 568                         for( int i = 0; i < length; i ++ ) {\r
 569                                 data[ j ] = buffer[ i ];\r
 570                                 j += inc;\r
 571                         }\r
 572                         UnlockBufferCF( ref buffer );\r
 573                 }\r
 574 \r
 575                 //======================================================================================\r
 576                 //======================================================================================\r
 577                 \r
 578                 static private bool                     _bufferCLocked  = false;\r
 579                 static private Complex[]        _bufferC                = new Complex[ 0 ];\r
 580 \r
 581                 static private void             LockBufferC( int length, ref Complex[] buffer ) {\r
 582                         Debug.Assert( length >= 0 );\r
 583                         Debug.Assert( _bufferCLocked == false );\r
 584                         \r
 585                         _bufferCLocked = true;\r
 586                         if( length >= _bufferC.Length ) {\r
 587                                 _bufferC        = new Complex[ length ];\r
 588                         }\r
 589                         buffer =        _bufferC;\r
 590                 }\r
 591                 static private void             UnlockBufferC( ref Complex[] buffer ) {\r
 592                         Debug.Assert( _bufferC == buffer );\r
 593                         Debug.Assert( _bufferCLocked == true );\r
 594                         \r
 595                         _bufferCLocked = false;\r
 596                         buffer = null;\r
 597                 }\r
 598 \r
 599                 private static void     LinearFFT( Complex[] data, int start, int inc, int length, FourierDirection direction ) {\r
 600                         Debug.Assert( data != null );\r
 601                         Debug.Assert( start >= 0 );\r
 602                         Debug.Assert( inc >= 1 );\r
 603                         Debug.Assert( length >= 1 );\r
 604                         Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );\r
 605                         \r
 606                         // copy to buffer\r
 607                         Complex[]       buffer = null;\r
 608                         LockBufferC( length, ref buffer );\r
 609                         int j = start;\r
 610                         for( int i = 0; i < length; i ++ ) {\r
 611                                 buffer[ i ] = data[ j ];\r
 612                                 j += inc;\r
 613                         }\r
 614 \r
 615                         FFT( buffer, length, direction );\r
 616 \r
 617                         // copy from buffer\r
 618                         j = start;\r
 619                         for( int i = 0; i < length; i ++ ) {\r
 620                                 data[ j ] = buffer[ i ];\r
 621                                 j += inc;\r
 622                         }\r
 623                         UnlockBufferC( ref buffer );\r
 624                 }\r
 625 \r
 626                 private static void     LinearFFT_Quick( Complex[] data, int start, int inc, int length, FourierDirection direction ) {\r
 627                         /*Debug.Assert( data != null );\r
 628                         Debug.Assert( start >= 0 );\r
 629                         Debug.Assert( inc >= 1 );\r
 630                         Debug.Assert( length >= 1 );\r
 631                         Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );*/\r
 632                         \r
 633                         // copy to buffer\r
 634                         Complex[]       buffer = null;\r
 635                         LockBufferC( length, ref buffer );\r
 636                         int j = start;\r
 637                         for( int i = 0; i < length; i ++ ) {\r
 638                                 buffer[ i ] = data[ j ];\r
 639                                 j += inc;\r
 640                         }\r
 641 \r
 642                         FFT_Quick( buffer, length, direction );\r
 643 \r
 644                         // copy from buffer\r
 645                         j = start;\r
 646                         for( int i = 0; i < length; i ++ ) {\r
 647                                 data[ j ] = buffer[ i ];\r
 648                                 j += inc;\r
 649                         }\r
 650                         UnlockBufferC( ref buffer );\r
 651                 }\r
 652 \r
 653                 //======================================================================================\r
 654                 //======================================================================================\r
 655 \r
 656                 /// <summary>\r
 657                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).\r
 658                 /// </summary>\r
 659                 /// <param name="data"></param>\r
 660                 /// <param name="length"></param>\r
 661                 /// <param name="direction"></param>\r
 662                 public static void      FFT( float[] data, int length, FourierDirection direction ) {\r
 663                         Debug.Assert( data != null );\r
 664                         Debug.Assert( data.Length >= length*2 );\r
 665                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 666 \r
 667                         Fourier.SyncLookupTableLength( length );\r
 668                         \r
 669                         int ln = Fourier.Log2( length );\r
 670 \r
 671                         // reorder array\r
 672                         Fourier.ReorderArray( data );\r
 673 \r
 674                         // successive doubling\r
 675                         int N = 1;\r
 676                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 677                         for( int level = 1; level <= ln; level ++ ) {\r
 678                                 int M = N;\r
 679                                 N <<= 1;\r
 680 \r
 681                                 float[] uRLookup = _uRLookupF[ level, signIndex ];\r
 682                                 float[] uILookup = _uILookupF[ level, signIndex ];\r
 683 \r
 684                                 for( int j = 0; j < M; j ++ ) {\r
 685                                         float uR = uRLookup[j];\r
 686                                         float uI = uILookup[j];\r
 687                                 \r
 688                                         for( int evenT = j; evenT < length; evenT += N ) {\r
 689                                                 int even = evenT << 1;\r
 690                                                 int odd = ( evenT + M ) << 1;\r
 691                                                 \r
 692                                                 float r = data[ odd ];\r
 693                                                 float i = data[ odd+1 ];\r
 694 \r
 695                                                 float odduR = r * uR - i * uI;\r
 696                                                 float odduI = r * uI + i * uR;\r
 697 \r
 698                                                 r = data[ even ];\r
 699                                                 i = data[ even+1 ];\r
 700 \r
 701                                                 data[ even ]    = r + odduR;\r
 702                                                 data[ even+1 ]  = i + odduI;\r
 703 \r
 704                                                 data[ odd ]             = r - odduR;\r
 705                                                 data[ odd+1 ]   = i - odduI;\r
 706                                         }\r
 707                                 }\r
 708                         }\r
 709                 }\r
 710                 \r
 711                 /// <summary>\r
 712                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).\r
 713                 /// </summary>\r
 714                 /// <param name="data"></param>\r
 715                 /// <param name="length"></param>\r
 716                 /// <param name="direction"></param>\r
 717                 public static void      FFT_Quick( float[] data, int length, FourierDirection direction ) {\r
 718                         /*Debug.Assert( data != null );\r
 719                         Debug.Assert( data.Length >= length*2 );\r
 720                         Debug.Assert( Fourier.IsPowerOf2( length ) == true );\r
 721 \r
 722                         Fourier.SyncLookupTableLength( length );*/\r
 723                         \r
 724                         int ln = Fourier.Log2( length );\r
 725 \r
 726                         // reorder array\r
 727                         Fourier.ReorderArray( data );\r
 728 \r
 729                         // successive doubling\r
 730                         int N = 1;\r
 731                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 732                         for( int level = 1; level <= ln; level ++ ) {\r
 733                                 int M = N;\r
 734                                 N <<= 1;\r
 735 \r
 736                                 float[] uRLookup = _uRLookupF[ level, signIndex ];\r
 737                                 float[] uILookup = _uILookupF[ level, signIndex ];\r
 738 \r
 739                                 for( int j = 0; j < M; j ++ ) {\r
 740                                         float uR = uRLookup[j];\r
 741                                         float uI = uILookup[j];\r
 742                                 \r
 743                                         for( int evenT = j; evenT < length; evenT += N ) {\r
 744                                                 int even = evenT << 1;\r
 745                                                 int odd = ( evenT + M ) << 1;\r
 746                                                 \r
 747                                                 float r = data[ odd ];\r
 748                                                 float i = data[ odd+1 ];\r
 749 \r
 750                                                 float odduR = r * uR - i * uI;\r
 751                                                 float odduI = r * uI + i * uR;\r
 752 \r
 753                                                 r = data[ even ];\r
 754                                                 i = data[ even+1 ];\r
 755 \r
 756                                                 data[ even ]    = r + odduR;\r
 757                                                 data[ even+1 ]  = i + odduI;\r
 758 \r
 759                                                 data[ odd ]             = r - odduR;\r
 760                                                 data[ odd+1 ]   = i - odduI;\r
 761                                         }\r
 762                                 }\r
 763                         }\r
 764                 }\r
 765 \r
 766                 /// <summary>\r
 767                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.\r
 768                 /// </summary>\r
 769                 /// <param name="data"></param>\r
 770                 /// <param name="length"></param>\r
 771                 /// <param name="direction"></param>\r
 772                 public static void      FFT( ComplexF[] data, int length, FourierDirection direction ) {\r
 773                         if( data == null ) {\r
 774                                 throw new ArgumentNullException( "data" );\r
 775                         }\r
 776                         if( data.Length < length ) {\r
 777                                 throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );\r
 778                         }\r
 779                         if( Fourier.IsPowerOf2( length ) == false ) {\r
 780                                 throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );\r
 781                         }\r
 782 \r
 783                         Fourier.SyncLookupTableLength( length );\r
 784 \r
 785                         int ln = Fourier.Log2( length );\r
 786                         \r
 787                         // reorder array\r
 788                         Fourier.ReorderArray( data );\r
 789                         \r
 790                         // successive doubling\r
 791                         int N = 1;\r
 792                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 793                         \r
 794                         for( int level = 1; level <= ln; level ++ ) {\r
 795                                 int M = N;\r
 796                                 N <<= 1;\r
 797 \r
 798                                 float[] uRLookup = _uRLookupF[ level, signIndex ];\r
 799                                 float[] uILookup = _uILookupF[ level, signIndex ];\r
 800 \r
 801                                 for( int j = 0; j < M; j ++ ) {\r
 802                                         float uR = uRLookup[j];\r
 803                                         float uI = uILookup[j];\r
 804 \r
 805                                         for( int even = j; even < length; even += N ) {\r
 806                                                 int odd  = even + M;\r
 807                                                 \r
 808                                                 float   r = data[ odd ].Re;\r
 809                                                 float   i = data[ odd ].Im;\r
 810 \r
 811                                                 float   odduR = r * uR - i * uI;\r
 812                                                 float   odduI = r * uI + i * uR;\r
 813 \r
 814                                                 r = data[ even ].Re;\r
 815                                                 i = data[ even ].Im;\r
 816                                                 \r
 817                                                 data[ even ].Re = r + odduR;\r
 818                                                 data[ even ].Im = i + odduI;\r
 819                                                 \r
 820                                                 data[ odd ].Re  = r - odduR;\r
 821                                                 data[ odd ].Im  = i - odduI;\r
 822                                         }\r
 823                                 }\r
 824                         }\r
 825 \r
 826                 }\r
 827 \r
 828                 /// <summary>\r
 829                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.\r
 830                 /// </summary>\r
 831                 /// <param name="data"></param>\r
 832                 /// <param name="length"></param>\r
 833                 /// <param name="direction"></param>\r
 834                 public static void      FFT_Quick( ComplexF[] data, int length, FourierDirection direction ) {\r
 835                         /*if( data == null ) {\r
 836                                 throw new ArgumentNullException( "data" );\r
 837                         }\r
 838                         if( data.Length < length ) {\r
 839                                 throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );\r
 840                         }\r
 841                         if( Fourier.IsPowerOf2( length ) == false ) {\r
 842                                 throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );\r
 843                         }\r
 844 \r
 845                         Fourier.SyncLookupTableLength( length );*/\r
 846 \r
 847                         int ln = Fourier.Log2( length );\r
 848                         \r
 849                         // reorder array\r
 850                         Fourier.ReorderArray( data );\r
 851                         \r
 852                         // successive doubling\r
 853                         int N = 1;\r
 854                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 855                         \r
 856                         for( int level = 1; level <= ln; level ++ ) {\r
 857                                 int M = N;\r
 858                                 N <<= 1;\r
 859 \r
 860                                 float[] uRLookup = _uRLookupF[ level, signIndex ];\r
 861                                 float[] uILookup = _uILookupF[ level, signIndex ];\r
 862 \r
 863                                 for( int j = 0; j < M; j ++ ) {\r
 864                                         float uR = uRLookup[j];\r
 865                                         float uI = uILookup[j];\r
 866 \r
 867                                         for( int even = j; even < length; even += N ) {\r
 868                                                 int odd  = even + M;\r
 869                                                 \r
 870                                                 float   r = data[ odd ].Re;\r
 871                                                 float   i = data[ odd ].Im;\r
 872 \r
 873                                                 float   odduR = r * uR - i * uI;\r
 874                                                 float   odduI = r * uI + i * uR;\r
 875 \r
 876                                                 r = data[ even ].Re;\r
 877                                                 i = data[ even ].Im;\r
 878                                                 \r
 879                                                 data[ even ].Re = r + odduR;\r
 880                                                 data[ even ].Im = i + odduI;\r
 881                                                 \r
 882                                                 data[ odd ].Re  = r - odduR;\r
 883                                                 data[ odd ].Im  = i - odduI;\r
 884                                         }\r
 885                                 }\r
 886                         }\r
 887 \r
 888                 }\r
 889 \r
 890                 /// <summary>\r
 891                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.\r
 892                 /// </summary>\r
 893                 /// <param name="data"></param>\r
 894                 /// <param name="direction"></param>\r
 895                 public static void      FFT( ComplexF[] data, FourierDirection direction ) {\r
 896                         if( data == null ) {\r
 897                                 throw new ArgumentNullException( "data" );\r
 898                         }\r
 899                         Fourier.FFT( data, data.Length, direction );\r
 900                 }\r
 901                 \r
 902                 /// <summary>\r
 903                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.\r
 904                 /// </summary>\r
 905                 /// <param name="data"></param>\r
 906                 /// <param name="length"></param>\r
 907                 /// <param name="direction"></param>\r
 908                 public static void      FFT( Complex[] data, int length, FourierDirection direction ) {\r
 909                         if( data == null ) {\r
 910                                 throw new ArgumentNullException( "data" );\r
 911                         }\r
 912                         if( data.Length < length ) {\r
 913                                 throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );\r
 914                         }\r
 915                         if( Fourier.IsPowerOf2( length ) == false ) {\r
 916                                 throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );\r
 917                         }\r
 918 \r
 919                         Fourier.SyncLookupTableLength( length );\r
 920 \r
 921                         int ln = Fourier.Log2( length );\r
 922                         \r
 923                         // reorder array\r
 924                         Fourier.ReorderArray( data );\r
 925                         \r
 926                         // successive doubling\r
 927                         int N = 1;\r
 928                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 929                         \r
 930                         for( int level = 1; level <= ln; level ++ ) {\r
 931                                 int M = N;\r
 932                                 N <<= 1;\r
 933 \r
 934                                 double[] uRLookup = _uRLookup[ level, signIndex ];\r
 935                                 double[] uILookup = _uILookup[ level, signIndex ];\r
 936 \r
 937                                 for( int j = 0; j < M; j ++ ) {\r
 938                                         double uR = uRLookup[j];\r
 939                                         double uI = uILookup[j];\r
 940 \r
 941                                         for( int even = j; even < length; even += N ) {\r
 942                                                 int odd  = even + M;\r
 943                                                 \r
 944                                                 double  r = data[ odd ].Re;\r
 945                                                 double  i = data[ odd ].Im;\r
 946 \r
 947                                                 double  odduR = r * uR - i * uI;\r
 948                                                 double  odduI = r * uI + i * uR;\r
 949 \r
 950                                                 r = data[ even ].Re;\r
 951                                                 i = data[ even ].Im;\r
 952                                                 \r
 953                                                 data[ even ].Re = r + odduR;\r
 954                                                 data[ even ].Im = i + odduI;\r
 955                                                 \r
 956                                                 data[ odd ].Re  = r - odduR;\r
 957                                                 data[ odd ].Im  = i - odduI;\r
 958                                         }\r
 959                                 }\r
 960                         }\r
 961 \r
 962                 }\r
 963                 \r
 964                 /// <summary>\r
 965                 /// Compute a 1D fast Fourier transform of a dataset of complex numbers.\r
 966                 /// </summary>\r
 967                 /// <param name="data"></param>\r
 968                 /// <param name="length"></param>\r
 969                 /// <param name="direction"></param>\r
 970                 public static void      FFT_Quick( Complex[] data, int length, FourierDirection direction ) {\r
 971                         /*if( data == null ) {\r
 972                                 throw new ArgumentNullException( "data" );\r
 973                         }\r
 974                         if( data.Length < length ) {\r
 975                                 throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );\r
 976                         }\r
 977                         if( Fourier.IsPowerOf2( length ) == false ) {\r
 978                                 throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );\r
 979                         }\r
 980 \r
 981                         Fourier.SyncLookupTableLength( length );   */\r
 982 \r
 983                         int ln = Fourier.Log2( length );\r
 984                         \r
 985                         // reorder array\r
 986                         Fourier.ReorderArray( data );\r
 987                         \r
 988                         // successive doubling\r
 989                         int N = 1;\r
 990                         int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;\r
 991                         \r
 992                         for( int level = 1; level <= ln; level ++ ) {\r
 993                                 int M = N;\r
 994                                 N <<= 1;\r
 995 \r
 996                                 double[] uRLookup = _uRLookup[ level, signIndex ];\r
 997                                 double[] uILookup = _uILookup[ level, signIndex ];\r
 998 \r
 999                                 for( int j = 0; j < M; j ++ ) {\r
1000                                         double uR = uRLookup[j];\r
1001                                         double uI = uILookup[j];\r
1002 \r
1003                                         for( int even = j; even < length; even += N ) {\r
1004                                                 int odd  = even + M;\r
1005                                                 \r
1006                                                 double  r = data[ odd ].Re;\r
1007                                                 double  i = data[ odd ].Im;\r
1008 \r
1009                                                 double  odduR = r * uR - i * uI;\r
1010                                                 double  odduI = r * uI + i * uR;\r
1011 \r
1012                                                 r = data[ even ].Re;\r
1013                                                 i = data[ even ].Im;\r
1014                                                 \r
1015                                                 data[ even ].Re = r + odduR;\r
1016                                                 data[ even ].Im = i + odduI;\r
1017                                                 \r
1018                                                 data[ odd ].Re  = r - odduR;\r
1019                                                 data[ odd ].Im  = i - odduI;\r
1020                                         }\r
1021                                 }\r
1022                         }\r
1023 \r
1024                 }\r
1025 \r
1026                 /// <summary>\r
1027                 /// Compute a 1D real-symmetric fast fourier transform.\r
1028                 /// </summary>\r
1029                 /// <param name="data"></param>\r
1030                 /// <param name="direction"></param>\r
1031                 public static void      RFFT( float[] data, FourierDirection direction ) {\r
1032                         if( data == null ) {\r
1033                                 throw new ArgumentNullException( "data" );\r
1034                         }\r
1035                         Fourier.RFFT( data, data.Length, direction );\r
1036                 }\r
1037                 \r
1038                 /// <summary>\r
1039                 /// Compute a 1D real-symmetric fast fourier transform.\r
1040                 /// </summary>\r
1041                 /// <param name="data"></param>\r
1042                 /// <param name="length"></param>\r
1043                 /// <param name="direction"></param>\r
1044                 public static void      RFFT( float[] data, int length, FourierDirection direction ) {\r
1045                         if( data == null ) {\r
1046                                 throw new ArgumentNullException( "data" );\r
1047                         }\r
1048                         if( data.Length < length ) {\r
1049                                 throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );\r
1050                         }\r
1051                         if( Fourier.IsPowerOf2( length ) == false ) {\r
1052                                 throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );\r
1053                         }\r
1054 \r
1055                         float   c1 = 0.5f, c2;\r
1056                         float   theta   = (float) Math.PI / (length/2);\r
1057 \r
1058                         if( direction == FourierDirection.Forward ) {\r
1059                                 c2 = -0.5f;\r
1060                                 FFT( data, length/2, direction );\r
1061                         }\r
1062                         else {\r
1063                                 c2 = 0.5f;\r
1064                                 theta = - theta;\r
1065                         }\r
1066 \r
1067                         float   wtemp = (float) Math.Sin( 0.5*theta );\r
1068                         float   wpr = -2 * wtemp*wtemp;\r
1069                         float   wpi     =(float)  Math.Sin( theta );\r
1070                         float   wr      = 1 + wpr;\r
1071                         float   wi      = wpi;\r
1072 \r
1073                         // do / undo packing\r
1074                         for( int i = 1; i < length/4; i ++ ) {\r
1075                                 int a = 2*i;\r
1076                                 int b = length - 2*i;\r
1077                                 float   h1r     = c1 * ( data[ a ] + data[ b ] );\r
1078                                 float   h1i     = c1 * ( data[ a+1 ] - data[ b+1 ] );\r
1079                                 float   h2r     = -c2 * ( data[ a+1 ] + data[ b+1 ] );\r
1080                                 float   h2i     = c2 * ( data[ a ] - data[ b ] );\r
1081                                 data[ a ]       = h1r + wr*h2r - wi*h2i;\r
1082                                 data[ a+1 ]     = h1i + wr*h2i + wi*h2r;\r
1083                                 data[ b ]       = h1r - wr*h2r + wi*h2i;\r
1084                                 data[ b+1 ]     = -h1i + wr*h2i + wi*h2r;\r
1085                                 wr = (wtemp = wr) * wpr - wi * wpi + wr;\r
1086                                 wi = wi * wpr + wtemp * wpi + wi;\r
1087                         }\r
1088                         \r
1089                         if( direction == FourierDirection.Forward ) {\r
1090                                 float  hir = data[0];\r
1091                                 data[0] = hir + data[1];\r
1092                                 data[1] = hir - data[1];\r
1093                         }\r
1094                         else {\r
1095                                 float  hir = data[0];\r
1096                                 data[0] = c1 * ( hir + data[1] );\r
1097                                 data[1] = c1 * ( hir - data[1] );\r
1098                                 Fourier.FFT( data, length/2, direction );\r
1099                         }\r
1100                 }\r
1101                 \r
1102                 /// <summary>\r
1103                 /// Compute a 2D fast fourier transform on a data set of complex numbers (represented as pairs of floats)\r
1104                 /// </summary>\r
1105                 /// <param name="data"></param>\r
1106                 /// <param name="xLength"></param>\r
1107                 /// <param name="yLength"></param>\r
1108                 /// <param name="direction"></param>\r
1109                 public static void      FFT2( float[] data, int xLength, int yLength, FourierDirection direction ) {\r
1110                         if( data == null ) {\r
1111                                 throw new ArgumentNullException( "data" );\r
1112                         }\r
1113                         if( data.Length < xLength*yLength*2 ) {\r
1114                                 throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * 2' parameter" );\r
1115                         }\r
1116                         if( Fourier.IsPowerOf2( xLength ) == false ) {\r
1117                                 throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );\r
1118                         }\r
1119                         if( Fourier.IsPowerOf2( yLength ) == false ) {\r
1120                                 throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );\r
1121                         }\r
1122 \r
1123                         int xInc = 1;\r
1124                         int yInc = xLength;\r
1125 \r
1126                         if( xLength > 1 ) {\r
1127                                 Fourier.SyncLookupTableLength( xLength );\r
1128                                 for( int y = 0; y < yLength; y ++ ) {\r
1129                                         int xStart = y * yInc;\r
1130                                         Fourier.LinearFFT_Quick( data, xStart, xInc, xLength, direction );\r
1131                                 }\r
1132                         }\r
1133                         \r
1134                         if( yLength > 1 ) {\r
1135                                 Fourier.SyncLookupTableLength( yLength );\r
1136                                 for( int x = 0; x < xLength; x ++ ) {\r
1137                                         int yStart = x * xInc;\r
1138                                         Fourier.LinearFFT_Quick( data, yStart, yInc, yLength, direction );\r
1139                                 }\r
1140                         }\r
1141                 }\r
1142 \r
1143                 /// <summary>\r
1144                 /// Compute a 2D fast fourier transform on a data set of complex numbers\r
1145                 /// </summary>\r
1146                 /// <param name="data"></param>\r
1147                 /// <param name="xLength"></param>\r
1148                 /// <param name="yLength"></param>\r
1149                 /// <param name="direction"></param>\r
1150                 public static void      FFT2( ComplexF[] data, int xLength, int yLength, FourierDirection direction ) {\r
1151                         if( data == null ) {\r
1152                                 throw new ArgumentNullException( "data" );\r
1153                         }\r
1154                         if( data.Length < xLength*yLength ) {\r
1155                                 throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter" );\r
1156                         }\r
1157                         if( Fourier.IsPowerOf2( xLength ) == false ) {\r
1158                                 throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );\r
1159                         }\r
1160                         if( Fourier.IsPowerOf2( yLength ) == false ) {\r
1161                                 throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );\r
1162                         }\r
1163 \r
1164                         int xInc = 1;\r
1165                         int yInc = xLength;\r
1166 \r
1167                         if( xLength > 1 ) {\r
1168                                 Fourier.SyncLookupTableLength( xLength );\r
1169                                 for( int y = 0; y < yLength; y ++ ) {\r
1170                                         int xStart = y * yInc;\r
1171                                         Fourier.LinearFFT_Quick( data, xStart, xInc, xLength, direction );\r
1172                                 }\r
1173                         }\r
1174                         \r
1175                         if( yLength > 1 ) {\r
1176                                 Fourier.SyncLookupTableLength( yLength );\r
1177                                 for( int x = 0; x < xLength; x ++ ) {\r
1178                                         int yStart = x * xInc;\r
1179                                         Fourier.LinearFFT_Quick( data, yStart, yInc, yLength, direction );\r
1180                                 }\r
1181                         }\r
1182                 }\r
1183 \r
1184                 /// <summary>\r
1185                 /// Compute a 2D fast fourier transform on a data set of complex numbers\r
1186                 /// </summary>\r
1187                 /// <param name="data"></param>\r
1188                 /// <param name="xLength"></param>\r
1189                 /// <param name="yLength"></param>\r
1190                 /// <param name="direction"></param>\r
1191                 public static void      FFT2( Complex[] data, int xLength, int yLength, FourierDirection direction ) {\r
1192                         if( data == null ) {\r
1193                                 throw new ArgumentNullException( "data" );\r
1194                         }\r
1195                         if( data.Length < xLength*yLength ) {\r
1196                                 throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter" );\r
1197                         }\r
1198                         if( Fourier.IsPowerOf2( xLength ) == false ) {\r
1199                                 throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );\r
1200                         }\r
1201                         if( Fourier.IsPowerOf2( yLength ) == false ) {\r
1202                                 throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );\r
1203                         }\r
1204 \r
1205                         int xInc = 1;\r
1206                         int yInc = xLength;\r
1207 \r
1208                         if( xLength > 1 ) {\r
1209                                 Fourier.SyncLookupTableLength( xLength );\r
1210                                 for( int y = 0; y < yLength; y ++ ) {\r
1211                                         int xStart = y * yInc;\r
1212                                         Fourier.LinearFFT_Quick( data, xStart, xInc, xLength, direction );\r
1213                                 }\r
1214                         }\r
1215                         \r
1216                         if( yLength > 1 ) {\r
1217                                 Fourier.SyncLookupTableLength( yLength );\r
1218                                 for( int x = 0; x < xLength; x ++ ) {\r
1219                                         int yStart = x * xInc;\r
1220                                         Fourier.LinearFFT_Quick( data, yStart, yInc, yLength, direction );\r
1221                                 }\r
1222                         }\r
1223                 }\r
1224 \r
1225                 /// <summary>\r
1226                 /// Compute a 3D fast fourier transform on a data set of complex numbers\r
1227                 /// </summary>\r
1228                 /// <param name="data"></param>\r
1229                 /// <param name="xLength"></param>\r
1230                 /// <param name="yLength"></param>\r
1231                 /// <param name="zLength"></param>\r
1232                 /// <param name="direction"></param>\r
1233                 public static void      FFT3( ComplexF[] data, int xLength, int yLength, int zLength, FourierDirection direction ) {\r
1234                         if( data == null ) {\r
1235                                 throw new ArgumentNullException( "data" );\r
1236                         }\r
1237                         if( data.Length < xLength*yLength*zLength ) {\r
1238                                 throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter" );\r
1239                         }\r
1240                         if( Fourier.IsPowerOf2( xLength ) == false ) {\r
1241                                 throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );\r
1242                         }\r
1243                         if( Fourier.IsPowerOf2( yLength ) == false ) {\r
1244                                 throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );\r
1245                         }\r
1246                         if( Fourier.IsPowerOf2( zLength ) == false ) {\r
1247                                 throw new ArgumentOutOfRangeException( "zLength", zLength, "must be a power of 2" );\r
1248                         }\r
1249 \r
1250                         int xInc = 1;\r
1251                         int yInc = xLength;\r
1252                         int zInc = xLength * yLength;\r
1253 \r
1254                         if( xLength > 1 ) {\r
1255                                 Fourier.SyncLookupTableLength( xLength );\r
1256                                 for( int z = 0; z < zLength; z ++ ) {\r
1257                                         for( int y = 0; y < yLength; y ++ ) {\r
1258                                                 int xStart = y * yInc + z * zInc;\r
1259                                                 Fourier.LinearFFT_Quick( data, xStart, xInc, xLength, direction );\r
1260                                         }\r
1261                                 }\r
1262                         }\r
1263                         \r
1264                         if( yLength > 1 ) {\r
1265                                 Fourier.SyncLookupTableLength( yLength );\r
1266                                 for( int z = 0; z < zLength; z ++ ) {\r
1267                                         for( int x = 0; x < xLength; x ++ ) {\r
1268                                                 int yStart = z * zInc + x * xInc;\r
1269                                                 Fourier.LinearFFT_Quick( data, yStart, yInc, yLength, direction );\r
1270                                         }\r
1271                                 }\r
1272                         }\r
1273                         \r
1274                         if( zLength > 1 ) {\r
1275                                 Fourier.SyncLookupTableLength( zLength );\r
1276                                 for( int y = 0; y < yLength; y ++ ) {\r
1277                                         for( int x = 0; x < xLength; x ++ ) {\r
1278                                                 int zStart = y * yInc + x * xInc;\r
1279                                                 Fourier.LinearFFT_Quick( data, zStart, zInc, zLength, direction );\r
1280                                         }\r
1281                                 }\r
1282                         }\r
1283                 }\r
1284 \r
1285                 /// <summary>\r
1286                 /// Compute a 3D fast fourier transform on a data set of complex numbers\r
1287                 /// </summary>\r
1288                 /// <param name="data"></param>\r
1289                 /// <param name="xLength"></param>\r
1290                 /// <param name="yLength"></param>\r
1291                 /// <param name="zLength"></param>\r
1292                 /// <param name="direction"></param>\r
1293                 public static void      FFT3( Complex[] data, int xLength, int yLength, int zLength, FourierDirection direction ) {\r
1294                         if( data == null ) {\r
1295                                 throw new ArgumentNullException( "data" );\r
1296                         }\r
1297                         if( data.Length < xLength*yLength*zLength ) {\r
1298                                 throw new ArgumentOutOfRangeException( "data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter" );\r
1299                         }\r
1300                         if( Fourier.IsPowerOf2( xLength ) == false ) {\r
1301                                 throw new ArgumentOutOfRangeException( "xLength", xLength, "must be a power of 2" );\r
1302                         }\r
1303                         if( Fourier.IsPowerOf2( yLength ) == false ) {\r
1304                                 throw new ArgumentOutOfRangeException( "yLength", yLength, "must be a power of 2" );\r
1305                         }\r
1306                         if( Fourier.IsPowerOf2( zLength ) == false ) {\r
1307                                 throw new ArgumentOutOfRangeException( "zLength", zLength, "must be a power of 2" );\r
1308                         }\r
1309 \r
1310                         int xInc = 1;\r
1311                         int yInc = xLength;\r
1312                         int zInc = xLength * yLength;\r
1313 \r
1314                         if( xLength > 1 ) {\r
1315                                 Fourier.SyncLookupTableLength( xLength );\r
1316                                 for( int z = 0; z < zLength; z ++ ) {\r
1317                                         for( int y = 0; y < yLength; y ++ ) {\r
1318                                                 int xStart = y * yInc + z * zInc;\r
1319                                                 Fourier.LinearFFT_Quick( data, xStart, xInc, xLength, direction );\r
1320                                         }\r
1321                                 }\r
1322                         }\r
1323                         \r
1324                         if( yLength > 1 ) {\r
1325                                 Fourier.SyncLookupTableLength( yLength );\r
1326                                 for( int z = 0; z < zLength; z ++ ) {\r
1327                                         for( int x = 0; x < xLength; x ++ ) {\r
1328                                                 int yStart = z * zInc + x * xInc;\r
1329                                                 Fourier.LinearFFT_Quick( data, yStart, yInc, yLength, direction );\r
1330                                         }\r
1331                                 }\r
1332                         }\r
1333                         \r
1334                         if( zLength > 1 ) {\r
1335                                 Fourier.SyncLookupTableLength( zLength );\r
1336                                 for( int y = 0; y < yLength; y ++ ) {\r
1337                                         for( int x = 0; x < xLength; x ++ ) {\r
1338                                                 int zStart = y * yInc + x * xInc;\r
1339                                                 Fourier.LinearFFT_Quick( data, zStart, zInc, zLength, direction );\r
1340                                         }\r
1341                                 }\r
1342                         }\r
1343                 }\r
1344                 \r
1345         }\r
1346 }\r