1536 lines
52 KiB
C#
1536 lines
52 KiB
C#
/*
|
|
* BSD Licence:
|
|
* Copyright (c) 2001, 2002 Ben Houston [ ben@exocortex.org ]
|
|
* Exocortex Technologies [ www.exocortex.org ]
|
|
* All rights reserved.
|
|
*
|
|
* Redistribution and use in source and binary forms, with or without
|
|
* modification, are permitted provided that the following conditions are met:
|
|
*
|
|
* 1. Redistributions of source code must retain the above copyright notice,
|
|
* this list of conditions and the following disclaimer.
|
|
* 2. Redistributions in binary form must reproduce the above copyright
|
|
* notice, this list of conditions and the following disclaimer in the
|
|
* documentation and/or other materials provided with the distribution.
|
|
* 3. Neither the name of the <ORGANIZATION> nor the names of its contributors
|
|
* may be used to endorse or promote products derived from this software
|
|
* without specific prior written permission.
|
|
*
|
|
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
|
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
|
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
|
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
|
|
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
|
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
|
|
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
|
|
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
|
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
|
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
|
|
* DAMAGE.
|
|
*/
|
|
|
|
using System;
|
|
using System.Diagnostics;
|
|
using System.IO;
|
|
using System.Text;
|
|
|
|
//using Exocortex.Imaging;
|
|
|
|
namespace Exocortex.DSP
|
|
{
|
|
|
|
// Comments? Questions? Bugs? Tell Ben Houston at ben@exocortex.org
|
|
// Version: May 4, 2002
|
|
|
|
/// <summary>
|
|
/// <p>Static functions for doing various Fourier Operations.</p>
|
|
/// </summary>
|
|
public class Fourier
|
|
{
|
|
|
|
//======================================================================================
|
|
|
|
private Fourier()
|
|
{
|
|
}
|
|
|
|
//======================================================================================
|
|
|
|
static private void Swap(ref float a, ref float b)
|
|
{
|
|
float temp = a;
|
|
a = b;
|
|
b = temp;
|
|
}
|
|
static private void Swap(ref double a, ref double b)
|
|
{
|
|
double temp = a;
|
|
a = b;
|
|
b = temp;
|
|
}
|
|
static private void Swap(ref ComplexF a, ref ComplexF b)
|
|
{
|
|
ComplexF temp = a;
|
|
a = b;
|
|
b = temp;
|
|
}
|
|
static private void Swap(ref Complex a, ref Complex b)
|
|
{
|
|
Complex temp = a;
|
|
a = b;
|
|
b = temp;
|
|
}
|
|
|
|
//-------------------------------------------------------------------------------------
|
|
|
|
private const int cMaxLength = 16777216;
|
|
private const int cMinLength = 1;
|
|
|
|
private const int cMaxBits = 24;
|
|
private const int cMinBits = 0;
|
|
|
|
|
|
static private bool IsPowerOf2(int x)
|
|
{
|
|
return (x & (x - 1)) == 0;
|
|
//return ( x == Pow2( Log2( x ) ) );
|
|
}
|
|
static private int Pow2(int exponent)
|
|
{
|
|
if (exponent >= 0 && exponent < 31)
|
|
{
|
|
return 1 << exponent;
|
|
}
|
|
return 0;
|
|
}
|
|
static private int Log2(int x)
|
|
{
|
|
if (x <= 65536)
|
|
{
|
|
if (x <= 256)
|
|
{
|
|
if (x <= 16)
|
|
{
|
|
if (x <= 4)
|
|
{
|
|
if (x <= 2)
|
|
{
|
|
if (x <= 1)
|
|
{
|
|
return 0;
|
|
}
|
|
return 1;
|
|
}
|
|
return 2;
|
|
}
|
|
if (x <= 8)
|
|
return 3;
|
|
return 4;
|
|
}
|
|
if (x <= 64)
|
|
{
|
|
if (x <= 32)
|
|
return 5;
|
|
return 6;
|
|
}
|
|
if (x <= 128)
|
|
return 7;
|
|
return 8;
|
|
}
|
|
if (x <= 4096)
|
|
{
|
|
if (x <= 1024)
|
|
{
|
|
if (x <= 512)
|
|
return 9;
|
|
return 10;
|
|
}
|
|
if (x <= 2048)
|
|
return 11;
|
|
return 12;
|
|
}
|
|
if (x <= 16384)
|
|
{
|
|
if (x <= 8192)
|
|
return 13;
|
|
return 14;
|
|
}
|
|
if (x <= 32768)
|
|
return 15;
|
|
return 16;
|
|
}
|
|
if (x <= 16777216)
|
|
{
|
|
if (x <= 1048576)
|
|
{
|
|
if (x <= 262144)
|
|
{
|
|
if (x <= 131072)
|
|
return 17;
|
|
return 18;
|
|
}
|
|
if (x <= 524288)
|
|
return 19;
|
|
return 20;
|
|
}
|
|
if (x <= 4194304)
|
|
{
|
|
if (x <= 2097152)
|
|
return 21;
|
|
return 22;
|
|
}
|
|
if (x <= 8388608)
|
|
return 23;
|
|
return 24;
|
|
}
|
|
if (x <= 268435456)
|
|
{
|
|
if (x <= 67108864)
|
|
{
|
|
if (x <= 33554432)
|
|
return 25;
|
|
return 26;
|
|
}
|
|
if (x <= 134217728)
|
|
return 27;
|
|
return 28;
|
|
}
|
|
if (x <= 1073741824)
|
|
{
|
|
if (x <= 536870912)
|
|
return 29;
|
|
return 30;
|
|
}
|
|
// since int is unsigned it can never be higher than 2,147,483,647
|
|
// if( x <= 2147483648 )
|
|
// return 31;
|
|
// return 32;
|
|
return 31;
|
|
}
|
|
|
|
//-------------------------------------------------------------------------------------
|
|
//-------------------------------------------------------------------------------------
|
|
|
|
static private int ReverseBits(int index, int numberOfBits)
|
|
{
|
|
Debug.Assert(numberOfBits >= cMinBits);
|
|
Debug.Assert(numberOfBits <= cMaxBits);
|
|
|
|
int reversedIndex = 0;
|
|
for (int i = 0; i < numberOfBits; i++)
|
|
{
|
|
reversedIndex = (reversedIndex << 1) | (index & 1);
|
|
index = (index >> 1);
|
|
}
|
|
return reversedIndex;
|
|
}
|
|
|
|
//-------------------------------------------------------------------------------------
|
|
|
|
static private int[][] _reversedBits = new int[cMaxBits][];
|
|
static private int[] GetReversedBits(int numberOfBits)
|
|
{
|
|
Debug.Assert(numberOfBits >= cMinBits);
|
|
Debug.Assert(numberOfBits <= cMaxBits);
|
|
if (_reversedBits[numberOfBits - 1] == null)
|
|
{
|
|
int maxBits = Fourier.Pow2(numberOfBits);
|
|
int[] reversedBits = new int[maxBits];
|
|
for (int i = 0; i < maxBits; i++)
|
|
{
|
|
int oldBits = i;
|
|
int newBits = 0;
|
|
for (int j = 0; j < numberOfBits; j++)
|
|
{
|
|
newBits = (newBits << 1) | (oldBits & 1);
|
|
oldBits = (oldBits >> 1);
|
|
}
|
|
reversedBits[i] = newBits;
|
|
}
|
|
_reversedBits[numberOfBits - 1] = reversedBits;
|
|
}
|
|
return _reversedBits[numberOfBits - 1];
|
|
}
|
|
|
|
//-------------------------------------------------------------------------------------
|
|
|
|
static private void ReorderArray(float[] data)
|
|
{
|
|
Debug.Assert(data != null);
|
|
|
|
int length = data.Length / 2;
|
|
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
Debug.Assert(length >= cMinLength);
|
|
Debug.Assert(length <= cMaxLength);
|
|
|
|
int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
int swap = reversedBits[i];
|
|
if (swap > i)
|
|
{
|
|
Fourier.Swap(ref data[(i << 1)], ref data[(swap << 1)]);
|
|
Fourier.Swap(ref data[(i << 1) + 1], ref data[(swap << 1) + 1]);
|
|
}
|
|
}
|
|
}
|
|
|
|
static private void ReorderArray(double[] data)
|
|
{
|
|
Debug.Assert(data != null);
|
|
|
|
int length = data.Length / 2;
|
|
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
Debug.Assert(length >= cMinLength);
|
|
Debug.Assert(length <= cMaxLength);
|
|
|
|
int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
int swap = reversedBits[i];
|
|
if (swap > i)
|
|
{
|
|
Fourier.Swap(ref data[i << 1], ref data[swap << 1]);
|
|
Fourier.Swap(ref data[i << 1 + 1], ref data[swap << 1 + 1]);
|
|
}
|
|
}
|
|
}
|
|
|
|
static private void ReorderArray(Complex[] data)
|
|
{
|
|
Debug.Assert(data != null);
|
|
|
|
int length = data.Length;
|
|
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
Debug.Assert(length >= cMinLength);
|
|
Debug.Assert(length <= cMaxLength);
|
|
|
|
int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
int swap = reversedBits[i];
|
|
if (swap > i)
|
|
{
|
|
Complex temp = data[i];
|
|
data[i] = data[swap];
|
|
data[swap] = temp;
|
|
}
|
|
}
|
|
}
|
|
|
|
static private void ReorderArray(ComplexF[] data)
|
|
{
|
|
Debug.Assert(data != null);
|
|
|
|
int length = data.Length;
|
|
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
Debug.Assert(length >= cMinLength);
|
|
Debug.Assert(length <= cMaxLength);
|
|
|
|
int[] reversedBits = Fourier.GetReversedBits(Fourier.Log2(length));
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
int swap = reversedBits[i];
|
|
if (swap > i)
|
|
{
|
|
ComplexF temp = data[i];
|
|
data[i] = data[swap];
|
|
data[swap] = temp;
|
|
}
|
|
}
|
|
}
|
|
|
|
//======================================================================================
|
|
|
|
private static int[][] _reverseBits = null;
|
|
|
|
private static int _ReverseBits(int bits, int n)
|
|
{
|
|
int bitsReversed = 0;
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
bitsReversed = (bitsReversed << 1) | (bits & 1);
|
|
bits = (bits >> 1);
|
|
}
|
|
return bitsReversed;
|
|
}
|
|
|
|
private static void InitializeReverseBits(int levels)
|
|
{
|
|
_reverseBits = new int[levels + 1][];
|
|
for (int j = 0; j < (levels + 1); j++)
|
|
{
|
|
int count = (int)Math.Pow(2, j);
|
|
_reverseBits[j] = new int[count];
|
|
for (int i = 0; i < count; i++)
|
|
{
|
|
_reverseBits[j][i] = _ReverseBits(i, j);
|
|
}
|
|
}
|
|
}
|
|
|
|
private static int _lookupTabletLength = -1;
|
|
private static double[,][] _uRLookup = null;
|
|
private static double[,][] _uILookup = null;
|
|
private static float[,][] _uRLookupF = null;
|
|
private static float[,][] _uILookupF = null;
|
|
|
|
private static void SyncLookupTableLength(int length)
|
|
{
|
|
// Debug.Assert( length < 1024*10 );
|
|
Debug.Assert(length >= 0);
|
|
if (length > _lookupTabletLength)
|
|
{
|
|
int level = (int)Math.Ceiling(Math.Log(length, 2));
|
|
Fourier.InitializeReverseBits(level);
|
|
Fourier.InitializeComplexRotations(level);
|
|
//_cFFTData = new Complex[ Math2.CeilingBase( length, 2 ) ];
|
|
//_cFFTDataF = new ComplexF[ Math2.CeilingBase( length, 2 ) ];
|
|
_lookupTabletLength = length;
|
|
}
|
|
}
|
|
|
|
private static int GetLookupTableLength()
|
|
{
|
|
return _lookupTabletLength;
|
|
}
|
|
|
|
private static void ClearLookupTables()
|
|
{
|
|
_uRLookup = null;
|
|
_uILookup = null;
|
|
_uRLookupF = null;
|
|
_uILookupF = null;
|
|
_lookupTabletLength = -1;
|
|
}
|
|
|
|
private static void InitializeComplexRotations(int levels)
|
|
{
|
|
int ln = levels;
|
|
//_wRLookup = new float[ levels + 1, 2 ];
|
|
//_wILookup = new float[ levels + 1, 2 ];
|
|
|
|
_uRLookup = new double[levels + 1, 2][];
|
|
_uILookup = new double[levels + 1, 2][];
|
|
|
|
_uRLookupF = new float[levels + 1, 2][];
|
|
_uILookupF = new float[levels + 1, 2][];
|
|
|
|
int N = 1;
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
//float scale = (float)( 1 / Math.Sqrt( 1 << ln ) );
|
|
|
|
// positive sign ( i.e. [M,0] )
|
|
{
|
|
double uR = 1;
|
|
double uI = 0;
|
|
double angle = Math.PI / M * 1;
|
|
double wR = (double)Math.Cos(angle);
|
|
double wI = (double)Math.Sin(angle);
|
|
|
|
_uRLookup[level, 0] = new double[M];
|
|
_uILookup[level, 0] = new double[M];
|
|
_uRLookupF[level, 0] = new float[M];
|
|
_uILookupF[level, 0] = new float[M];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
_uRLookupF[level, 0][j] = (float)(_uRLookup[level, 0][j] = uR);
|
|
_uILookupF[level, 0][j] = (float)(_uILookup[level, 0][j] = uI);
|
|
double uwI = uR * wI + uI * wR;
|
|
uR = uR * wR - uI * wI;
|
|
uI = uwI;
|
|
}
|
|
}
|
|
{
|
|
|
|
|
|
// negative sign ( i.e. [M,1] )
|
|
double uR = 1;
|
|
double uI = 0;
|
|
double angle = Math.PI / M * -1;
|
|
double wR = (double)Math.Cos(angle);
|
|
double wI = (double)Math.Sin(angle);
|
|
|
|
_uRLookup[level, 1] = new double[M];
|
|
_uILookup[level, 1] = new double[M];
|
|
_uRLookupF[level, 1] = new float[M];
|
|
_uILookupF[level, 1] = new float[M];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
_uRLookupF[level, 1][j] = (float)(_uRLookup[level, 1][j] = uR);
|
|
_uILookupF[level, 1][j] = (float)(_uILookup[level, 1][j] = uI);
|
|
double uwI = uR * wI + uI * wR;
|
|
uR = uR * wR - uI * wI;
|
|
uI = uwI;
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
//======================================================================================
|
|
//======================================================================================
|
|
|
|
static private bool _bufferFLocked = false;
|
|
static private float[] _bufferF = new float[0];
|
|
|
|
static private void LockBufferF(int length, ref float[] buffer)
|
|
{
|
|
Debug.Assert(_bufferFLocked == false);
|
|
_bufferFLocked = true;
|
|
if (length >= _bufferF.Length)
|
|
{
|
|
_bufferF = new float[length];
|
|
}
|
|
buffer = _bufferF;
|
|
}
|
|
static private void UnlockBufferF(ref float[] buffer)
|
|
{
|
|
Debug.Assert(_bufferF == buffer);
|
|
Debug.Assert(_bufferFLocked == true);
|
|
_bufferFLocked = false;
|
|
buffer = null;
|
|
}
|
|
|
|
private static void LinearFFT(float[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
Debug.Assert(data != null);
|
|
Debug.Assert(start >= 0);
|
|
Debug.Assert(inc >= 1);
|
|
Debug.Assert(length >= 1);
|
|
Debug.Assert((start + inc * (length - 1)) * 2 < data.Length);
|
|
|
|
// copy to buffer
|
|
float[] buffer = null;
|
|
LockBufferF(length * 2, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length * 2; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferF(ref buffer);
|
|
}
|
|
|
|
private static void LinearFFT_Quick(float[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
/*Debug.Assert( data != null );
|
|
Debug.Assert( start >= 0 );
|
|
Debug.Assert( inc >= 1 );
|
|
Debug.Assert( length >= 1 );
|
|
Debug.Assert( ( start + inc * ( length - 1 ) ) * 2 < data.Length );*/
|
|
|
|
// copy to buffer
|
|
float[] buffer = null;
|
|
LockBufferF(length * 2, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length * 2; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT_Quick(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferF(ref buffer);
|
|
}
|
|
|
|
//======================================================================================
|
|
//======================================================================================
|
|
|
|
static private bool _bufferCFLocked = false;
|
|
static private ComplexF[] _bufferCF = new ComplexF[0];
|
|
|
|
static private void LockBufferCF(int length, ref ComplexF[] buffer)
|
|
{
|
|
Debug.Assert(length >= 0);
|
|
Debug.Assert(_bufferCFLocked == false);
|
|
|
|
_bufferCFLocked = true;
|
|
if (length != _bufferCF.Length)
|
|
{
|
|
_bufferCF = new ComplexF[length];
|
|
}
|
|
buffer = _bufferCF;
|
|
}
|
|
static private void UnlockBufferCF(ref ComplexF[] buffer)
|
|
{
|
|
Debug.Assert(_bufferCF == buffer);
|
|
Debug.Assert(_bufferCFLocked == true);
|
|
|
|
_bufferCFLocked = false;
|
|
buffer = null;
|
|
}
|
|
|
|
private static void LinearFFT(ComplexF[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
Debug.Assert(data != null);
|
|
Debug.Assert(start >= 0);
|
|
Debug.Assert(inc >= 1);
|
|
Debug.Assert(length >= 1);
|
|
Debug.Assert((start + inc * (length - 1)) < data.Length);
|
|
|
|
// copy to buffer
|
|
ComplexF[] buffer = null;
|
|
LockBufferCF(length, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferCF(ref buffer);
|
|
}
|
|
|
|
private static void LinearFFT_Quick(ComplexF[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
/*Debug.Assert( data != null );
|
|
Debug.Assert( start >= 0 );
|
|
Debug.Assert( inc >= 1 );
|
|
Debug.Assert( length >= 1 );
|
|
Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length ); */
|
|
|
|
// copy to buffer
|
|
ComplexF[] buffer = null;
|
|
LockBufferCF(length, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferCF(ref buffer);
|
|
}
|
|
|
|
//======================================================================================
|
|
//======================================================================================
|
|
|
|
static private bool _bufferCLocked = false;
|
|
static private Complex[] _bufferC = new Complex[0];
|
|
|
|
static private void LockBufferC(int length, ref Complex[] buffer)
|
|
{
|
|
Debug.Assert(length >= 0);
|
|
Debug.Assert(_bufferCLocked == false);
|
|
|
|
_bufferCLocked = true;
|
|
if (length >= _bufferC.Length)
|
|
{
|
|
_bufferC = new Complex[length];
|
|
}
|
|
buffer = _bufferC;
|
|
}
|
|
static private void UnlockBufferC(ref Complex[] buffer)
|
|
{
|
|
Debug.Assert(_bufferC == buffer);
|
|
Debug.Assert(_bufferCLocked == true);
|
|
|
|
_bufferCLocked = false;
|
|
buffer = null;
|
|
}
|
|
|
|
private static void LinearFFT(Complex[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
Debug.Assert(data != null);
|
|
Debug.Assert(start >= 0);
|
|
Debug.Assert(inc >= 1);
|
|
Debug.Assert(length >= 1);
|
|
Debug.Assert((start + inc * (length - 1)) < data.Length);
|
|
|
|
// copy to buffer
|
|
Complex[] buffer = null;
|
|
LockBufferC(length, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferC(ref buffer);
|
|
}
|
|
|
|
private static void LinearFFT_Quick(Complex[] data, int start, int inc, int length, FourierDirection direction)
|
|
{
|
|
/*Debug.Assert( data != null );
|
|
Debug.Assert( start >= 0 );
|
|
Debug.Assert( inc >= 1 );
|
|
Debug.Assert( length >= 1 );
|
|
Debug.Assert( ( start + inc * ( length - 1 ) ) < data.Length );*/
|
|
|
|
// copy to buffer
|
|
Complex[] buffer = null;
|
|
LockBufferC(length, ref buffer);
|
|
int j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
buffer[i] = data[j];
|
|
j += inc;
|
|
}
|
|
|
|
FFT_Quick(buffer, length, direction);
|
|
|
|
// copy from buffer
|
|
j = start;
|
|
for (int i = 0; i < length; i++)
|
|
{
|
|
data[j] = buffer[i];
|
|
j += inc;
|
|
}
|
|
UnlockBufferC(ref buffer);
|
|
}
|
|
|
|
//======================================================================================
|
|
//======================================================================================
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT(float[] data, int length, FourierDirection direction)
|
|
{
|
|
Debug.Assert(data != null);
|
|
Debug.Assert(data.Length >= length * 2);
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
|
|
Fourier.SyncLookupTableLength(length);
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
float[] uRLookup = _uRLookupF[level, signIndex];
|
|
float[] uILookup = _uILookupF[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
float uR = uRLookup[j];
|
|
float uI = uILookup[j];
|
|
|
|
for (int evenT = j; evenT < length; evenT += N)
|
|
{
|
|
int even = evenT << 1;
|
|
int odd = (evenT + M) << 1;
|
|
|
|
float r = data[odd];
|
|
float i = data[odd + 1];
|
|
|
|
float odduR = r * uR - i * uI;
|
|
float odduI = r * uI + i * uR;
|
|
|
|
r = data[even];
|
|
i = data[even + 1];
|
|
|
|
data[even] = r + odduR;
|
|
data[even + 1] = i + odduI;
|
|
|
|
data[odd] = r - odduR;
|
|
data[odd + 1] = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers (as pairs of float's).
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT_Quick(float[] data, int length, FourierDirection direction)
|
|
{
|
|
Debug.Assert(data != null);
|
|
Debug.Assert(data.Length >= length * 2);
|
|
Debug.Assert(Fourier.IsPowerOf2(length) == true);
|
|
|
|
Fourier.SyncLookupTableLength(length);
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
float[] uRLookup = _uRLookupF[level, signIndex];
|
|
float[] uILookup = _uILookupF[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
float uR = uRLookup[j];
|
|
float uI = uILookup[j];
|
|
|
|
for (int evenT = j; evenT < length; evenT += N)
|
|
{
|
|
int even = evenT << 1;
|
|
int odd = (evenT + M) << 1;
|
|
|
|
float r = data[odd];
|
|
float i = data[odd + 1];
|
|
|
|
float odduR = r * uR - i * uI;
|
|
float odduI = r * uI + i * uR;
|
|
|
|
r = data[even];
|
|
i = data[even + 1];
|
|
|
|
data[even] = r + odduR;
|
|
data[even + 1] = i + odduI;
|
|
|
|
data[odd] = r - odduR;
|
|
data[odd + 1] = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT(ComplexF[] data, int length, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < length)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(length) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
|
|
}
|
|
|
|
Fourier.SyncLookupTableLength(length);
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
float[] uRLookup = _uRLookupF[level, signIndex];
|
|
float[] uILookup = _uILookupF[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
float uR = uRLookup[j];
|
|
float uI = uILookup[j];
|
|
|
|
for (int even = j; even < length; even += N)
|
|
{
|
|
int odd = even + M;
|
|
|
|
float r = data[odd].Re;
|
|
float i = data[odd].Im;
|
|
|
|
float odduR = r * uR - i * uI;
|
|
float odduI = r * uI + i * uR;
|
|
|
|
r = data[even].Re;
|
|
i = data[even].Im;
|
|
|
|
data[even].Re = r + odduR;
|
|
data[even].Im = i + odduI;
|
|
|
|
data[odd].Re = r - odduR;
|
|
data[odd].Im = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT_Quick(ComplexF[] data, int length, FourierDirection direction)
|
|
{
|
|
/*if( data == null ) {
|
|
throw new ArgumentNullException( "data" );
|
|
}
|
|
if( data.Length < length ) {
|
|
throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );
|
|
}
|
|
if( Fourier.IsPowerOf2( length ) == false ) {
|
|
throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );
|
|
}
|
|
|
|
Fourier.SyncLookupTableLength( length );*/
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
float[] uRLookup = _uRLookupF[level, signIndex];
|
|
float[] uILookup = _uILookupF[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
float uR = uRLookup[j];
|
|
float uI = uILookup[j];
|
|
|
|
for (int even = j; even < length; even += N)
|
|
{
|
|
int odd = even + M;
|
|
|
|
float r = data[odd].Re;
|
|
float i = data[odd].Im;
|
|
|
|
float odduR = r * uR - i * uI;
|
|
float odduI = r * uI + i * uR;
|
|
|
|
r = data[even].Re;
|
|
i = data[even].Im;
|
|
|
|
data[even].Re = r + odduR;
|
|
data[even].Im = i + odduI;
|
|
|
|
data[odd].Re = r - odduR;
|
|
data[odd].Im = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT(ComplexF[] data, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
Fourier.FFT(data, data.Length, direction);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT(Complex[] data, int length, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < length)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(length) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
|
|
}
|
|
|
|
Fourier.SyncLookupTableLength(length);
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
double[] uRLookup = _uRLookup[level, signIndex];
|
|
double[] uILookup = _uILookup[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
double uR = uRLookup[j];
|
|
double uI = uILookup[j];
|
|
|
|
for (int even = j; even < length; even += N)
|
|
{
|
|
int odd = even + M;
|
|
|
|
double r = data[odd].Re;
|
|
double i = data[odd].Im;
|
|
|
|
double odduR = r * uR - i * uI;
|
|
double odduI = r * uI + i * uR;
|
|
|
|
r = data[even].Re;
|
|
i = data[even].Im;
|
|
|
|
data[even].Re = r + odduR;
|
|
data[even].Im = i + odduI;
|
|
|
|
data[odd].Re = r - odduR;
|
|
data[odd].Im = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT_Quick(Complex[] data, int length, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < length)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(length) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
|
|
}
|
|
|
|
Fourier.SyncLookupTableLength(length);
|
|
|
|
int ln = Fourier.Log2(length);
|
|
|
|
// reorder array
|
|
Fourier.ReorderArray(data);
|
|
|
|
// successive doubling
|
|
int N = 1;
|
|
int signIndex = (direction == FourierDirection.Forward) ? 0 : 1;
|
|
|
|
for (int level = 1; level <= ln; level++)
|
|
{
|
|
int M = N;
|
|
N <<= 1;
|
|
|
|
double[] uRLookup = _uRLookup[level, signIndex];
|
|
double[] uILookup = _uILookup[level, signIndex];
|
|
|
|
for (int j = 0; j < M; j++)
|
|
{
|
|
double uR = uRLookup[j];
|
|
double uI = uILookup[j];
|
|
|
|
for (int even = j; even < length; even += N)
|
|
{
|
|
int odd = even + M;
|
|
|
|
double r = data[odd].Re;
|
|
double i = data[odd].Im;
|
|
|
|
double odduR = r * uR - i * uI;
|
|
double odduI = r * uI + i * uR;
|
|
|
|
r = data[even].Re;
|
|
i = data[even].Im;
|
|
|
|
data[even].Re = r + odduR;
|
|
data[even].Im = i + odduI;
|
|
|
|
data[odd].Re = r - odduR;
|
|
data[odd].Im = i - odduI;
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D real-symmetric fast fourier transform.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="direction"></param>
|
|
public static void RFFT(float[] data, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
Fourier.RFFT(data, data.Length, direction);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 1D real-symmetric fast fourier transform.
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="length"></param>
|
|
/// <param name="direction"></param>
|
|
public static void RFFT(float[] data, int length, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < length)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be at least as large as 'data.Length' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(length) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("length", length, "must be a power of 2");
|
|
}
|
|
|
|
float c1 = 0.5f, c2;
|
|
float theta = (float)Math.PI / (length / 2);
|
|
|
|
if (direction == FourierDirection.Forward)
|
|
{
|
|
c2 = -0.5f;
|
|
FFT(data, length / 2, direction);
|
|
}
|
|
else
|
|
{
|
|
c2 = 0.5f;
|
|
theta = -theta;
|
|
}
|
|
|
|
float wtemp = (float)Math.Sin(0.5 * theta);
|
|
float wpr = -2 * wtemp * wtemp;
|
|
float wpi = (float)Math.Sin(theta);
|
|
float wr = 1 + wpr;
|
|
float wi = wpi;
|
|
|
|
// do / undo packing
|
|
for (int i = 1; i < length / 4; i++)
|
|
{
|
|
int a = 2 * i;
|
|
int b = length - 2 * i;
|
|
float h1r = c1 * (data[a] + data[b]);
|
|
float h1i = c1 * (data[a + 1] - data[b + 1]);
|
|
float h2r = -c2 * (data[a + 1] + data[b + 1]);
|
|
float h2i = c2 * (data[a] - data[b]);
|
|
data[a] = h1r + wr * h2r - wi * h2i;
|
|
data[a + 1] = h1i + wr * h2i + wi * h2r;
|
|
data[b] = h1r - wr * h2r + wi * h2i;
|
|
data[b + 1] = -h1i + wr * h2i + wi * h2r;
|
|
wr = (wtemp = wr) * wpr - wi * wpi + wr;
|
|
wi = wi * wpr + wtemp * wpi + wi;
|
|
}
|
|
|
|
if (direction == FourierDirection.Forward)
|
|
{
|
|
float hir = data[0];
|
|
data[0] = hir + data[1];
|
|
data[1] = hir - data[1];
|
|
}
|
|
else
|
|
{
|
|
float hir = data[0];
|
|
data[0] = c1 * (hir + data[1]);
|
|
data[1] = c1 * (hir - data[1]);
|
|
Fourier.FFT(data, length / 2, direction);
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 2D fast fourier transform on a data set of complex numbers (represented as pairs of floats)
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="xLength"></param>
|
|
/// <param name="yLength"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT2(float[] data, int xLength, int yLength, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < xLength * yLength * 2)
|
|
{
|
|
throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * 2' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(xLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(yLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
|
|
}
|
|
|
|
int xInc = 1;
|
|
int yInc = xLength;
|
|
|
|
if (xLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(xLength);
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
int xStart = y * yInc;
|
|
Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
|
|
}
|
|
}
|
|
|
|
if (yLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(yLength);
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int yStart = x * xInc;
|
|
Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 2D fast fourier transform on a data set of complex numbers
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="xLength"></param>
|
|
/// <param name="yLength"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT2(ComplexF[] data, int xLength, int yLength, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < xLength * yLength)
|
|
{
|
|
throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(xLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(yLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
|
|
}
|
|
|
|
int xInc = 1;
|
|
int yInc = xLength;
|
|
|
|
if (xLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(xLength);
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
int xStart = y * yInc;
|
|
Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
|
|
}
|
|
}
|
|
|
|
if (yLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(yLength);
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int yStart = x * xInc;
|
|
Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 2D fast fourier transform on a data set of complex numbers
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="xLength"></param>
|
|
/// <param name="yLength"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT2(Complex[] data, int xLength, int yLength, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < xLength * yLength)
|
|
{
|
|
throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(xLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(yLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
|
|
}
|
|
|
|
int xInc = 1;
|
|
int yInc = xLength;
|
|
|
|
if (xLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(xLength);
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
int xStart = y * yInc;
|
|
Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
|
|
}
|
|
}
|
|
|
|
if (yLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(yLength);
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int yStart = x * xInc;
|
|
Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 3D fast fourier transform on a data set of complex numbers
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="xLength"></param>
|
|
/// <param name="yLength"></param>
|
|
/// <param name="zLength"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT3(ComplexF[] data, int xLength, int yLength, int zLength, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < xLength * yLength * zLength)
|
|
{
|
|
throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(xLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(yLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(zLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("zLength", zLength, "must be a power of 2");
|
|
}
|
|
|
|
int xInc = 1;
|
|
int yInc = xLength;
|
|
int zInc = xLength * yLength;
|
|
|
|
if (xLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(xLength);
|
|
for (int z = 0; z < zLength; z++)
|
|
{
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
int xStart = y * yInc + z * zInc;
|
|
Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (yLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(yLength);
|
|
for (int z = 0; z < zLength; z++)
|
|
{
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int yStart = z * zInc + x * xInc;
|
|
Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (zLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(zLength);
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int zStart = y * yInc + x * xInc;
|
|
Fourier.LinearFFT_Quick(data, zStart, zInc, zLength, direction);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute a 3D fast fourier transform on a data set of complex numbers
|
|
/// </summary>
|
|
/// <param name="data"></param>
|
|
/// <param name="xLength"></param>
|
|
/// <param name="yLength"></param>
|
|
/// <param name="zLength"></param>
|
|
/// <param name="direction"></param>
|
|
public static void FFT3(Complex[] data, int xLength, int yLength, int zLength, FourierDirection direction)
|
|
{
|
|
if (data == null)
|
|
{
|
|
throw new ArgumentNullException("data");
|
|
}
|
|
if (data.Length < xLength * yLength * zLength)
|
|
{
|
|
throw new ArgumentOutOfRangeException("data.Length", data.Length, "must be at least as large as 'xLength * yLength * zLength' parameter");
|
|
}
|
|
if (Fourier.IsPowerOf2(xLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("xLength", xLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(yLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("yLength", yLength, "must be a power of 2");
|
|
}
|
|
if (Fourier.IsPowerOf2(zLength) == false)
|
|
{
|
|
throw new ArgumentOutOfRangeException("zLength", zLength, "must be a power of 2");
|
|
}
|
|
|
|
int xInc = 1;
|
|
int yInc = xLength;
|
|
int zInc = xLength * yLength;
|
|
|
|
if (xLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(xLength);
|
|
for (int z = 0; z < zLength; z++)
|
|
{
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
int xStart = y * yInc + z * zInc;
|
|
Fourier.LinearFFT_Quick(data, xStart, xInc, xLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (yLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(yLength);
|
|
for (int z = 0; z < zLength; z++)
|
|
{
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int yStart = z * zInc + x * xInc;
|
|
Fourier.LinearFFT_Quick(data, yStart, yInc, yLength, direction);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (zLength > 1)
|
|
{
|
|
Fourier.SyncLookupTableLength(zLength);
|
|
for (int y = 0; y < yLength; y++)
|
|
{
|
|
for (int x = 0; x < xLength; x++)
|
|
{
|
|
int zStart = y * yInc + x * xInc;
|
|
Fourier.LinearFFT_Quick(data, zStart, zInc, zLength, direction);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|