Browse Source

Single precision Fourier Transform support (#481)

v3
AlexHild 10 years ago
committed by Christoph Ruegg
parent
commit
aa83f4862e
  1. 10
      src/Numerics/ComplexExtensions.cs
  2. 126
      src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
  3. 30
      src/Numerics/IntegralTransforms/Fourier.Naive.cs
  4. 74
      src/Numerics/IntegralTransforms/Fourier.RadixN.cs
  5. 468
      src/Numerics/IntegralTransforms/Fourier.cs
  6. 10
      src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs
  7. 104
      src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs
  8. 119
      src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs
  9. 49
      src/UnitTests/IntegralTransformsTests/FourierTest.cs
  10. 99
      src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs
  11. 222
      src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs
  12. 23
      src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs

10
src/Numerics/ComplexExtensions.cs

@ -45,6 +45,16 @@ namespace MathNet.Numerics
/// </summary>
public static class ComplexExtensions
{
/// <summary>
/// Gets the squared magnitude of the <c>Complex</c> number.
/// </summary>
/// <param name="complex">The <see cref="Complex32"/> number to perfom this operation on.</param>
/// <returns>The squared magnitude of the <c>Complex</c> number.</returns>
public static double MagnitudeSquared(this Complex32 complex)
{
return (complex.Real * complex.Real) + (complex.Imaginary * complex.Imaginary);
}
/// <summary>
/// Gets the squared magnitude of the <c>Complex</c> number.
/// </summary>

126
src/Numerics/IntegralTransforms/Fourier.Bluestein.cs

@ -47,6 +47,38 @@ namespace MathNet.Numerics.IntegralTransforms
/// </summary>
const int BluesteinSequenceLengthThreshold = 46341;
/// <summary>
/// Generate the bluestein sequence for the provided problem size.
/// </summary>
/// <param name="n">Number of samples.</param>
/// <returns>Bluestein sequence exp(I*Pi*k^2/N)</returns>
static Complex32[] BluesteinSequence32(int n)
{
double s = Constants.Pi / n;
var sequence = new Complex32[n];
// TODO: benchmark whether the second variation is significantly
// faster than the former one. If not just use the former one always.
if (n > BluesteinSequenceLengthThreshold)
{
for (int k = 0; k < sequence.Length; k++)
{
double t = (s * k) * k;
sequence[k] = new Complex32((float)Math.Cos(t), (float)Math.Sin(t));
}
}
else
{
for (int k = 0; k < sequence.Length; k++)
{
double t = s * (k * k);
sequence[k] = new Complex32((float)Math.Cos(t), (float)Math.Sin(t));
}
}
return sequence;
}
/// <summary>
/// Generate the bluestein sequence for the provided problem size.
/// </summary>
@ -79,6 +111,61 @@ namespace MathNet.Numerics.IntegralTransforms
return sequence;
}
/// <summary>
/// Convolution with the bluestein sequence (Parallel Version).
/// </summary>
/// <param name="samples">Sample Vector.</param>
static void BluesteinConvolutionParallel(Complex32[] samples)
{
int n = samples.Length;
Complex32[] sequence = BluesteinSequence32(n);
// Padding to power of two >= 2N–1 so we can apply Radix-2 FFT.
int m = ((n << 1) - 1).CeilingToPowerOfTwo();
var b = new Complex32[m];
var a = new Complex32[m];
CommonParallel.Invoke(
() =>
{
// Build and transform padded sequence b_k = exp(I*Pi*k^2/N)
for (int i = 0; i < n; i++)
{
b[i] = sequence[i];
}
for (int i = m - n + 1; i < b.Length; i++)
{
b[i] = sequence[m - i];
}
Radix2(b, -1);
},
() =>
{
// Build and transform padded sequence a_k = x_k * exp(-I*Pi*k^2/N)
for (int i = 0; i < samples.Length; i++)
{
a[i] = sequence[i].Conjugate() * samples[i];
}
Radix2(a, -1);
});
for (int i = 0; i < a.Length; i++)
{
a[i] *= b[i];
}
Radix2Parallel(a, 1);
var nbinv = 1.0f / m;
for (int i = 0; i < samples.Length; i++)
{
samples[i] = nbinv * sequence[i].Conjugate() * a[i];
}
}
/// <summary>
/// Convolution with the bluestein sequence (Parallel Version).
/// </summary>
@ -134,6 +221,18 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Swap the real and imaginary parts of each sample.
/// </summary>
/// <param name="samples">Sample Vector.</param>
static void SwapRealImaginary(Complex32[] samples)
{
for (int i = 0; i < samples.Length; i++)
{
samples[i] = new Complex32(samples[i].Imaginary, samples[i].Real);
}
}
/// <summary>
/// Swap the real and imaginary parts of each sample.
/// </summary>
@ -146,6 +245,33 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary>
/// <param name="samples">Time-space sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
internal static void Bluestein(Complex32[] samples, int exponentSign)
{
int n = samples.Length;
if (n.IsPowerOfTwo())
{
Radix2Parallel(samples, exponentSign);
return;
}
if (exponentSign == 1)
{
SwapRealImaginary(samples);
}
BluesteinConvolutionParallel(samples);
if (exponentSign == 1)
{
SwapRealImaginary(samples);
}
}
/// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary>

30
src/Numerics/IntegralTransforms/Fourier.Naive.cs

@ -41,6 +41,36 @@ namespace MathNet.Numerics.IntegralTransforms
/// </summary>
public static partial class Fourier
{
/// <summary>
/// Naive generic DFT, useful e.g. to verify faster algorithms.
/// </summary>
/// <param name="samples">Time-space sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <returns>Corresponding frequency-space vector.</returns>
internal static Complex32[] Naive(Complex32[] samples, int exponentSign)
{
var w0 = exponentSign * Constants.Pi2 / samples.Length;
var spectrum = new Complex32[samples.Length];
CommonParallel.For(0, samples.Length, (u, v) =>
{
for (int i = u; i < v; i++)
{
var wk = w0 * i;
var sum = Complex32.Zero;
for (var n = 0; n < samples.Length; n++)
{
var w = n * wk;
sum += samples[n] * new Complex32((float)Math.Cos(w), (float)Math.Sin(w));
}
spectrum[i] = sum;
}
});
return spectrum;
}
/// <summary>
/// Naive generic DFT, useful e.g. to verify faster algorithms.
/// </summary>

74
src/Numerics/IntegralTransforms/Fourier.RadixN.cs

@ -70,6 +70,29 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Radix-2 Step Helper Method
/// </summary>
/// <param name="samples">Sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <param name="levelSize">Level Group Size.</param>
/// <param name="k">Index inside of the level.</param>
static void Radix2Step(Complex32[] samples, int exponentSign, int levelSize, int k)
{
// Twiddle Factor
var exponent = (exponentSign * k) * Constants.Pi / levelSize;
var w = new Complex32((float)Math.Cos(exponent), (float)Math.Sin(exponent));
var step = levelSize << 1;
for (var i = k; i < samples.Length; i += step)
{
var ai = samples[i];
var t = w * samples[i + levelSize];
samples[i] = ai + t;
samples[i + levelSize] = ai - t;
}
}
/// <summary>
/// Radix-2 Step Helper Method
/// </summary>
@ -93,6 +116,29 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
internal static void Radix2(Complex32[] samples, int exponentSign)
{
if (!samples.Length.IsPowerOfTwo())
{
throw new ArgumentException(Resources.ArgumentPowerOfTwo);
}
Radix2Reorder(samples);
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
{
for (var k = 0; k < levelSize; k++)
{
Radix2Step(samples, exponentSign, levelSize, k);
}
}
}
/// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
@ -116,6 +162,34 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
internal static void Radix2Parallel(Complex32[] samples, int exponentSign)
{
if (!samples.Length.IsPowerOfTwo())
{
throw new ArgumentException(Resources.ArgumentPowerOfTwo);
}
Radix2Reorder(samples);
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
{
var size = levelSize;
CommonParallel.For(0, size, 64, (u, v) =>
{
for (int i = u; i < v; i++)
{
Radix2Step(samples, exponentSign, size, i);
}
});
}
}
/// <summary>
/// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary>

468
src/Numerics/IntegralTransforms/Fourier.cs

@ -44,6 +44,15 @@ namespace MathNet.Numerics.IntegralTransforms
/// </summary>
public static partial class Fourier
{
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
public static void Forward(Complex32[] samples)
{
Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.SymmetricScaling);
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
/// </summary>
@ -52,6 +61,32 @@ namespace MathNet.Numerics.IntegralTransforms
{
Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.SymmetricScaling);
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Forward(Complex32[] samples, FourierOptions options)
{
switch (options)
{
case FourierOptions.NoScaling:
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.NoScaling);
break;
case FourierOptions.InverseExponent:
Control.FourierTransformProvider.Backward(samples, FourierTransformScaling.SymmetricScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.NoScaling:
case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.Backward(samples, FourierTransformScaling.NoScaling);
break;
default:
Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
@ -78,6 +113,37 @@ namespace MathNet.Numerics.IntegralTransforms
break;
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="real">Real part of the sample vector, where the FFT is evaluated in place.</param>
/// <param name="imaginary">Imaginary part of the sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Forward(float[] real, float[] imaginary, FourierOptions options = FourierOptions.Default)
{
if (real.Length != imaginary.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
// TODO: consider to support this natively by the provider, without the need for copying
// TODO: otherwise, consider ArrayPool
Complex32[] data = new Complex32[real.Length];
for (int i = 0; i < data.Length; i++)
{
data[i] = new Complex32(real[i], imaginary[i]);
}
Forward(data, options);
for (int i = 0; i < data.Length; i++)
{
real[i] = data[i].Real;
imaginary[i] = data[i].Imaginary;
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
@ -109,6 +175,40 @@ namespace MathNet.Numerics.IntegralTransforms
imaginary[i] = data[i].Imaginary;
}
}
/// <summary>
/// Packed Real-Complex forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
/// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
/// the spectrum can be fully reconstructed form the positive frequencies only (first half).
/// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
/// </summary>
/// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
/// <param name="n">The number of samples.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void ForwardReal(float[] data, int n, FourierOptions options = FourierOptions.Default)
{
int length = n.IsEven() ? n + 2 : n + 1;
if (data.Length < length)
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, length));
}
if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
{
throw new NotSupportedException();
}
switch (options)
{
case FourierOptions.NoScaling:
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.ForwardReal(data, n, FourierTransformScaling.NoScaling);
break;
default:
Control.FourierTransformProvider.ForwardReal(data, n, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Packed Real-Complex forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors.
@ -144,6 +244,36 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data.
/// </summary>
/// <param name="samples">Sample data, where the FFT is evaluated in place.</param>
/// <param name="dimensions">
/// The data size per dimension. The first dimension is the major one.
/// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
/// </param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void ForwardMultiDim(Complex32[] samples, int[] dimensions, FourierOptions options = FourierOptions.Default)
{
switch (options)
{
case FourierOptions.NoScaling:
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.ForwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
break;
case FourierOptions.InverseExponent:
Control.FourierTransformProvider.BackwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.NoScaling:
case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.BackwardMultidim(samples, dimensions, FourierTransformScaling.NoScaling);
break;
default:
Control.FourierTransformProvider.ForwardMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data.
/// </summary>
@ -173,6 +303,19 @@ namespace MathNet.Numerics.IntegralTransforms
break;
}
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to two dimensional sample data.
/// </summary>
/// <param name="samplesRowWise">Sample data, organized row by row, where the FFT is evaluated in place</param>
/// <param name="rows">The number of rows.</param>
/// <param name="columns">The number of columns.</param>
/// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Forward2D(Complex32[] samplesRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
{
ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options);
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to two dimensional sample data.
@ -186,6 +329,34 @@ namespace MathNet.Numerics.IntegralTransforms
{
ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options);
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to a two dimensional data in form of a matrix.
/// </summary>
/// <param name="samples">Sample matrix, where the FFT is evaluated in place</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Forward2D(Matrix<Complex32> samples, FourierOptions options = FourierOptions.Default)
{
var rowMajorArray = samples.AsRowMajorArray();
if (rowMajorArray != null)
{
ForwardMultiDim(rowMajorArray, new[] { samples.RowCount, samples.ColumnCount }, options);
return;
}
var columnMajorArray = samples.AsColumnMajorArray();
if (columnMajorArray != null)
{
ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
return;
}
// Fall Back
columnMajorArray = samples.ToColumnMajorArray();
ForwardMultiDim(columnMajorArray, new[] { samples.ColumnCount, samples.RowCount }, options);
var denseStorage = new DenseColumnMajorMatrixStorage<Complex32>(samples.RowCount, samples.ColumnCount, columnMajorArray);
denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear);
}
/// <summary>
/// Applies the forward Fast Fourier Transform (FFT) to a two dimensional data in form of a matrix.
@ -215,6 +386,15 @@ namespace MathNet.Numerics.IntegralTransforms
denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
public static void Inverse(Complex32[] spectrum)
{
Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
@ -224,6 +404,36 @@ namespace MathNet.Numerics.IntegralTransforms
Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Inverse(Complex32[] spectrum, FourierOptions options)
{
switch (options)
{
case FourierOptions.NoScaling:
Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.NoScaling);
break;
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.BackwardScaling);
break;
case FourierOptions.InverseExponent:
Control.FourierTransformProvider.Forward(spectrum, FourierTransformScaling.SymmetricScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.NoScaling:
Control.FourierTransformProvider.Forward(spectrum, FourierTransformScaling.NoScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.Forward(spectrum, FourierTransformScaling.ForwardScaling);
break;
default:
Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
@ -254,6 +464,37 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
/// <param name="real">Real part of the sample vector, where the iFFT is evaluated in place.</param>
/// <param name="imaginary">Imaginary part of the sample vector, where the iFFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Inverse(float[] real, float[] imaginary, FourierOptions options = FourierOptions.Default)
{
if (real.Length != imaginary.Length)
{
throw new ArgumentException(Resources.ArgumentArraysSameLength);
}
// TODO: consider to support this natively by the provider, without the need for copying
// TODO: otherwise, consider ArrayPool
Complex32[] data = new Complex32[real.Length];
for (int i = 0; i < data.Length; i++)
{
data[i] = new Complex32(real[i], imaginary[i]);
}
Inverse(data, options);
for (int i = 0; i < data.Length; i++)
{
real[i] = data[i].Real;
imaginary[i] = data[i].Imaginary;
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// </summary>
@ -285,6 +526,42 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Packed Real-Complex inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
/// the spectrum can be fully reconstructed form the positive frequencies only (first half).
/// The data array needs to be N+2 (if N is even) or N+1 (if N is odd) long in order to support such a packed spectrum.
/// </summary>
/// <param name="data">Data array of length N+2 (if N is even) or N+1 (if N is odd).</param>
/// <param name="n">The number of samples.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void InverseReal(float[] data, int n, FourierOptions options = FourierOptions.Default)
{
int length = n.IsEven() ? n + 2 : n + 1;
if (data.Length < length)
{
throw new ArgumentException(string.Format(Resources.ArrayTooSmall, length));
}
if ((options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent)
{
throw new NotSupportedException();
}
switch (options)
{
case FourierOptions.NoScaling:
Control.FourierTransformProvider.BackwardReal(data, n, FourierTransformScaling.NoScaling);
break;
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.BackwardReal(data, n, FourierTransformScaling.BackwardScaling);
break;
default:
Control.FourierTransformProvider.BackwardReal(data, n, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Packed Real-Complex inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors.
/// Since for real-valued time samples the complex spectrum is conjugate-even (symmetry),
@ -321,6 +598,40 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data.
/// </summary>
/// <param name="spectrum">Spectrum data, where the iFFT is evaluated in place.</param>
/// <param name="dimensions">
/// The data size per dimension. The first dimension is the major one.
/// For example, with two dimensions "rows" and "columns" the samples are assumed to be organized row by row.
/// </param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void InverseMultiDim(Complex32[] spectrum, int[] dimensions, FourierOptions options = FourierOptions.Default)
{
switch (options)
{
case FourierOptions.NoScaling:
Control.FourierTransformProvider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
break;
case FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.BackwardScaling);
break;
case FourierOptions.InverseExponent:
Control.FourierTransformProvider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.NoScaling:
Control.FourierTransformProvider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.NoScaling);
break;
case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling:
Control.FourierTransformProvider.ForwardMultidim(spectrum, dimensions, FourierTransformScaling.ForwardScaling);
break;
default:
Control.FourierTransformProvider.BackwardMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling);
break;
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data.
/// </summary>
@ -355,6 +666,19 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data.
/// </summary>
/// <param name="spectrumRowWise">Sample data, organized row by row, where the iFFT is evaluated in place</param>
/// <param name="rows">The number of rows.</param>
/// <param name="columns">The number of columns.</param>
/// <remarks>Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments.</remarks>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Inverse2D(Complex32[] spectrumRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default)
{
InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data.
/// </summary>
@ -368,6 +692,34 @@ namespace MathNet.Numerics.IntegralTransforms
InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix.
/// </summary>
/// <param name="spectrum">Sample matrix, where the iFFT is evaluated in place</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void Inverse2D(Matrix<Complex32> spectrum, FourierOptions options = FourierOptions.Default)
{
var rowMajorArray = spectrum.AsRowMajorArray();
if (rowMajorArray != null)
{
InverseMultiDim(rowMajorArray, new[] { spectrum.RowCount, spectrum.ColumnCount }, options);
return;
}
var columnMajorArray = spectrum.AsColumnMajorArray();
if (columnMajorArray != null)
{
InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
return;
}
// Fall Back
columnMajorArray = spectrum.ToColumnMajorArray();
InverseMultiDim(columnMajorArray, new[] { spectrum.ColumnCount, spectrum.RowCount }, options);
var denseStorage = new DenseColumnMajorMatrixStorage<Complex32>(spectrum.RowCount, spectrum.ColumnCount, columnMajorArray);
denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear);
}
/// <summary>
/// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix.
/// </summary>
@ -396,6 +748,19 @@ namespace MathNet.Numerics.IntegralTransforms
denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear);
}
/// <summary>
/// Naive forward DFT, useful e.g. to verify faster algorithms.
/// </summary>
/// <param name="samples">Time-space sample vector.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <returns>Corresponding frequency-space vector.</returns>
public static Complex32[] NaiveForward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
{
var frequencySpace = Naive(samples, SignByOptions(options));
ForwardScaleByOptions(options, frequencySpace);
return frequencySpace;
}
/// <summary>
/// Naive forward DFT, useful e.g. to verify faster algorithms.
/// </summary>
@ -409,6 +774,19 @@ namespace MathNet.Numerics.IntegralTransforms
return frequencySpace;
}
/// <summary>
/// Naive inverse DFT, useful e.g. to verify faster algorithms.
/// </summary>
/// <param name="spectrum">Frequency-space sample vector.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <returns>Corresponding time-space vector.</returns>
public static Complex32[] NaiveInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
{
var timeSpace = Naive(spectrum, -SignByOptions(options));
InverseScaleByOptions(options, timeSpace);
return timeSpace;
}
/// <summary>
/// Naive inverse DFT, useful e.g. to verify faster algorithms.
/// </summary>
@ -422,6 +800,18 @@ namespace MathNet.Numerics.IntegralTransforms
return timeSpace;
}
/// <summary>
/// Radix-2 forward FFT for power-of-two sized sample vectors.
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <exception cref="ArgumentException"/>
public static void Radix2Forward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
{
Radix2Parallel(samples, SignByOptions(options));
ForwardScaleByOptions(options, samples);
}
/// <summary>
/// Radix-2 forward FFT for power-of-two sized sample vectors.
/// </summary>
@ -434,6 +824,18 @@ namespace MathNet.Numerics.IntegralTransforms
ForwardScaleByOptions(options, samples);
}
/// <summary>
/// Radix-2 inverse FFT for power-of-two sized sample vectors.
/// </summary>
/// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <exception cref="ArgumentException"/>
public static void Radix2Inverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
{
Radix2Parallel(spectrum, -SignByOptions(options));
InverseScaleByOptions(options, spectrum);
}
/// <summary>
/// Radix-2 inverse FFT for power-of-two sized sample vectors.
/// </summary>
@ -446,6 +848,17 @@ namespace MathNet.Numerics.IntegralTransforms
InverseScaleByOptions(options, spectrum);
}
/// <summary>
/// Bluestein forward FFT for arbitrary sized sample vectors.
/// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void BluesteinForward(Complex32[] samples, FourierOptions options = FourierOptions.Default)
{
Bluestein(samples, SignByOptions(options));
ForwardScaleByOptions(options, samples);
}
/// <summary>
/// Bluestein forward FFT for arbitrary sized sample vectors.
/// </summary>
@ -457,6 +870,17 @@ namespace MathNet.Numerics.IntegralTransforms
ForwardScaleByOptions(options, samples);
}
/// <summary>
/// Bluestein inverse FFT for arbitrary sized sample vectors.
/// </summary>
/// <param name="spectrum">Sample vector, where the FFT is evaluated in place.</param>
/// <param name="options">Fourier Transform Convention Options.</param>
public static void BluesteinInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default)
{
Bluestein(spectrum, -SignByOptions(options));
InverseScaleByOptions(options, spectrum);
}
/// <summary>
/// Bluestein inverse FFT for arbitrary sized sample vectors.
/// </summary>
@ -479,6 +903,26 @@ namespace MathNet.Numerics.IntegralTransforms
return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1;
}
/// <summary>
/// Rescale FFT-the resulting vector according to the provided convention options.
/// </summary>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <param name="samples">Sample Vector.</param>
static void ForwardScaleByOptions(FourierOptions options, Complex32[] samples)
{
if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling ||
(options & FourierOptions.AsymmetricScaling) == FourierOptions.AsymmetricScaling)
{
return;
}
var scalingFactor = (float)Math.Sqrt(1.0 / samples.Length);
for (int i = 0; i < samples.Length; i++)
{
samples[i] *= scalingFactor;
}
}
/// <summary>
/// Rescale FFT-the resulting vector according to the provided convention options.
/// </summary>
@ -499,6 +943,30 @@ namespace MathNet.Numerics.IntegralTransforms
}
}
/// <summary>
/// Rescale the iFFT-resulting vector according to the provided convention options.
/// </summary>
/// <param name="options">Fourier Transform Convention Options.</param>
/// <param name="samples">Sample Vector.</param>
static void InverseScaleByOptions(FourierOptions options, Complex32[] samples)
{
if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling)
{
return;
}
var scalingFactor = (float)1.0 / samples.Length;
if ((options & FourierOptions.AsymmetricScaling) != FourierOptions.AsymmetricScaling)
{
scalingFactor = (float)Math.Sqrt(scalingFactor);
}
for (int i = 0; i < samples.Length; i++)
{
samples[i] *= scalingFactor;
}
}
/// <summary>
/// Rescale the iFFT-resulting vector according to the provided convention options.
/// </summary>

10
src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs

@ -55,13 +55,19 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// </summary>
void InitializeVerify();
void Forward(Complex32[] samples, FourierTransformScaling scaling);
void Forward(Complex[] samples, FourierTransformScaling scaling);
void Backward(Complex32[] spectrum, FourierTransformScaling scaling);
void Backward(Complex[] spectrum, FourierTransformScaling scaling);
void ForwardReal(float[] samples, int n, FourierTransformScaling scaling);
void ForwardReal(double[] samples, int n, FourierTransformScaling scaling);
void BackwardReal(float[] spectrum, int n, FourierTransformScaling scaling);
void BackwardReal(double[] spectrum, int n, FourierTransformScaling scaling);
void ForwardMultidim(Complex32[] samples, int[] dimensions, FourierTransformScaling scaling);
void ForwardMultidim(Complex[] samples, int[] dimensions, FourierTransformScaling scaling);
void BackwardMultidim(Complex32[] spectrum, int[] dimensions, FourierTransformScaling scaling);
void BackwardMultidim(Complex[] spectrum, int[] dimensions, FourierTransformScaling scaling);
}
}

104
src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs

@ -58,6 +58,23 @@ namespace MathNet.Numerics.Providers.FourierTransform
return "Managed";
}
public void Forward(Complex32[] samples, FourierTransformScaling scaling)
{
switch (scaling)
{
case FourierTransformScaling.SymmetricScaling:
Fourier.BluesteinForward(samples, FourierOptions.Default);
break;
case FourierTransformScaling.ForwardScaling:
// Only backward scaling can be expressed with options, hence the double-inverse
Fourier.BluesteinInverse(samples, FourierOptions.AsymmetricScaling | FourierOptions.InverseExponent);
break;
default:
Fourier.BluesteinForward(samples, FourierOptions.NoScaling);
break;
}
}
public void Forward(Complex[] samples, FourierTransformScaling scaling)
{
switch (scaling)
@ -75,6 +92,22 @@ namespace MathNet.Numerics.Providers.FourierTransform
}
}
public void Backward(Complex32[] spectrum, FourierTransformScaling scaling)
{
switch (scaling)
{
case FourierTransformScaling.SymmetricScaling:
Fourier.BluesteinInverse(spectrum, FourierOptions.Default);
break;
case FourierTransformScaling.BackwardScaling:
Fourier.BluesteinInverse(spectrum, FourierOptions.AsymmetricScaling);
break;
default:
Fourier.BluesteinInverse(spectrum, FourierOptions.NoScaling);
break;
}
}
public void Backward(Complex[] spectrum, FourierTransformScaling scaling)
{
switch (scaling)
@ -91,6 +124,37 @@ namespace MathNet.Numerics.Providers.FourierTransform
}
}
public void ForwardReal(float[] samples, int n, FourierTransformScaling scaling)
{
// TODO: backport proper, optimized implementation from Iridium
Complex32[] data = new Complex32[n];
for (int i = 0; i < data.Length; i++)
{
data[i] = new Complex32(samples[i], 0.0f);
}
Forward(data, scaling);
samples[0] = data[0].Real;
samples[1] = 0f;
for (int i = 1, j = 2; i < data.Length / 2; i++)
{
samples[j++] = data[i].Real;
samples[j++] = data[i].Imaginary;
}
if (n.IsEven())
{
samples[n] = data[data.Length / 2].Real;
samples[n + 1] = 0f;
}
else
{
samples[n - 1] = data[data.Length / 2].Real;
samples[n] = data[data.Length / 2].Imaginary;
}
}
public void ForwardReal(double[] samples, int n, FourierTransformScaling scaling)
{
// TODO: backport proper, optimized implementation from Iridium
@ -122,6 +186,36 @@ namespace MathNet.Numerics.Providers.FourierTransform
}
}
public void BackwardReal(float[] spectrum, int n, FourierTransformScaling scaling)
{
// TODO: backport proper, optimized implementation from Iridium
Complex32[] data = new Complex32[n];
data[0] = new Complex32(spectrum[0], 0f);
for (int i = 1, j = 2; i < data.Length / 2; i++)
{
data[i] = new Complex32(spectrum[j++], spectrum[j++]);
data[data.Length - i] = data[i].Conjugate();
}
if (n.IsEven())
{
data[data.Length / 2] = new Complex32(spectrum[n], 0f);
}
else
{
data[data.Length / 2] = new Complex32(spectrum[n - 1], spectrum[n]);
data[data.Length / 2 + 1] = data[data.Length / 2].Conjugate();
}
Backward(data, scaling);
for (int i = 0; i < data.Length; i++)
{
spectrum[i] = data[i].Real;
}
spectrum[n] = 0f;
}
public void BackwardReal(double[] spectrum, int n, FourierTransformScaling scaling)
{
// TODO: backport proper, optimized implementation from Iridium
@ -152,11 +246,21 @@ namespace MathNet.Numerics.Providers.FourierTransform
spectrum[n] = 0d;
}
public void ForwardMultidim(Complex32[] samples, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
}
public void ForwardMultidim(Complex[] samples, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
}
public void BackwardMultidim(Complex32[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
}
public void BackwardMultidim(Complex[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();

119
src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs

@ -43,6 +43,7 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
public int[] Dimensions;
public FourierTransformScaling Scaling;
public bool Real;
public bool Single;
}
Kernel _kernel;
@ -143,7 +144,7 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
return MklProvider.Describe();
}
Kernel Configure(int length, FourierTransformScaling scaling, bool real)
Kernel Configure(int length, FourierTransformScaling scaling, bool real, bool single)
{
Kernel kernel = Interlocked.Exchange(ref _kernel, null);
@ -153,36 +154,54 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
{
Dimensions = new[] {length},
Scaling = scaling,
Real = real
Real = real,
Single = single
};
if (real) SafeNativeMethods.d_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
else SafeNativeMethods.z_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
if (single)
{
if (real) SafeNativeMethods.s_fft_create(out kernel.Handle, length, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
else SafeNativeMethods.c_fft_create(out kernel.Handle, length, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
}
else
{
if (real) SafeNativeMethods.d_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
else SafeNativeMethods.z_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
}
return kernel;
}
if (kernel.Dimensions.Length != 1 || kernel.Dimensions[0] != length || kernel.Scaling != scaling || kernel.Real != real)
if (kernel.Dimensions.Length != 1 || kernel.Dimensions[0] != length || kernel.Scaling != scaling || kernel.Real != real || kernel.Single != single)
{
SafeNativeMethods.x_fft_free(ref kernel.Handle);
if (real) SafeNativeMethods.d_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
else SafeNativeMethods.z_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
if (single)
{
if (real) SafeNativeMethods.s_fft_create(out kernel.Handle, length, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
else SafeNativeMethods.c_fft_create(out kernel.Handle, length, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
}
else
{
if (real) SafeNativeMethods.d_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
else SafeNativeMethods.z_fft_create(out kernel.Handle, length, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
}
kernel.Dimensions = new[] {length};
kernel.Scaling = scaling;
kernel.Real = real;
kernel.Single = single;
return kernel;
}
return kernel;
}
Kernel Configure(int[] dimensions, FourierTransformScaling scaling)
Kernel Configure(int[] dimensions, FourierTransformScaling scaling, bool single)
{
if (dimensions.Length == 1)
{
return Configure(dimensions[0], scaling, false);
return Configure(dimensions[0], scaling, false, single);
}
Kernel kernel = Interlocked.Exchange(ref _kernel, null);
@ -193,7 +212,8 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
{
Dimensions = dimensions,
Scaling = scaling,
Real = false
Real = false,
Single = single
};
long length = 1;
@ -202,11 +222,19 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
length *= dimensions[i];
}
SafeNativeMethods.z_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
if (single)
{
SafeNativeMethods.c_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
}
else
{
SafeNativeMethods.z_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
}
return kernel;
}
bool mismatch = kernel.Dimensions.Length != dimensions.Length || kernel.Scaling != scaling || kernel.Real != false;
bool mismatch = kernel.Dimensions.Length != dimensions.Length || kernel.Scaling != scaling || kernel.Real != false || kernel.Single != single;
if (!mismatch)
{
for (int i = 0; i < dimensions.Length; i++)
@ -228,11 +256,20 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
}
SafeNativeMethods.x_fft_free(ref kernel.Handle);
SafeNativeMethods.z_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
if (single)
{
SafeNativeMethods.c_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, (float)ForwardScaling(scaling, length), (float)BackwardScaling(scaling, length));
}
else
{
SafeNativeMethods.z_fft_create_multidim(out kernel.Handle, dimensions.Length, dimensions, ForwardScaling(scaling, length), BackwardScaling(scaling, length));
}
kernel.Dimensions = dimensions;
kernel.Scaling = scaling;
kernel.Real = false;
kernel.Single = single;
return kernel;
}
@ -248,46 +285,90 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl
}
}
public void Forward(Complex32[] samples, FourierTransformScaling scaling)
{
Kernel kernel = Configure(samples.Length, scaling, false, true);
SafeNativeMethods.c_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void Forward(Complex[] samples, FourierTransformScaling scaling)
{
Kernel kernel = Configure(samples.Length, scaling, false);
Kernel kernel = Configure(samples.Length, scaling, false, false);
SafeNativeMethods.z_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void Backward(Complex32[] spectrum, FourierTransformScaling scaling)
{
Kernel kernel = Configure(spectrum.Length, scaling, false, true);
SafeNativeMethods.c_fft_backward(kernel.Handle, spectrum);
Release(kernel);
}
public void Backward(Complex[] spectrum, FourierTransformScaling scaling)
{
Kernel kernel = Configure(spectrum.Length, scaling, false);
Kernel kernel = Configure(spectrum.Length, scaling, false, false);
SafeNativeMethods.z_fft_backward(kernel.Handle, spectrum);
Release(kernel);
}
public void ForwardReal(float[] samples, int n, FourierTransformScaling scaling)
{
Kernel kernel = Configure(n, scaling, true, true);
SafeNativeMethods.s_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void ForwardReal(double[] samples, int n, FourierTransformScaling scaling)
{
Kernel kernel = Configure(n, scaling, true);
Kernel kernel = Configure(n, scaling, true, false);
SafeNativeMethods.d_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void BackwardReal(float[] spectrum, int n, FourierTransformScaling scaling)
{
Kernel kernel = Configure(n, scaling, true, true);
SafeNativeMethods.s_fft_backward(kernel.Handle, spectrum);
Release(kernel);
spectrum[n] = 0f;
}
public void BackwardReal(double[] spectrum, int n, FourierTransformScaling scaling)
{
Kernel kernel = Configure(n, scaling, true);
Kernel kernel = Configure(n, scaling, true, false);
SafeNativeMethods.d_fft_backward(kernel.Handle, spectrum);
Release(kernel);
spectrum[n] = 0d;
}
public void ForwardMultidim(Complex32[] samples, int[] dimensions, FourierTransformScaling scaling)
{
Kernel kernel = Configure(dimensions, scaling, true);
SafeNativeMethods.c_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void ForwardMultidim(Complex[] samples, int[] dimensions, FourierTransformScaling scaling)
{
Kernel kernel = Configure(dimensions, scaling);
Kernel kernel = Configure(dimensions, scaling, false);
SafeNativeMethods.z_fft_forward(kernel.Handle, samples);
Release(kernel);
}
public void BackwardMultidim(Complex32[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
Kernel kernel = Configure(dimensions, scaling, true);
SafeNativeMethods.c_fft_backward(kernel.Handle, spectrum);
Release(kernel);
}
public void BackwardMultidim(Complex[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
Kernel kernel = Configure(dimensions, scaling);
Kernel kernel = Configure(dimensions, scaling, false);
SafeNativeMethods.z_fft_backward(kernel.Handle, spectrum);
Release(kernel);
}

49
src/UnitTests/IntegralTransformsTests/FourierTest.cs

@ -53,6 +53,41 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
return new ContinuousUniform(-1, 1, new System.Random(seed));
}
/// <summary>
/// Naive transforms real sine correctly.
/// </summary>
[Test]
public void NaiveTransformsRealSineCorrectly32()
{
var samples = Generate.PeriodicMap(16, w => new Complex32((float)Math.Sin(w), 0), 16, 1.0, Constants.Pi2);
// real-odd transforms to imaginary odd
var spectrum = Fourier.NaiveForward(samples, FourierOptions.Matlab);
// all real components must be zero
foreach (var c in spectrum)
{
Assert.AreEqual(0, c.Real, 1e-6, "real");
}
// all imaginary components except second and last musth be zero
for (var i = 0; i < spectrum.Length; i++)
{
if (i == 1)
{
Assert.AreEqual(-8, spectrum[i].Imaginary, 1e-12, "imag second");
}
else if (i == spectrum.Length - 1)
{
Assert.AreEqual(8, spectrum[i].Imaginary, 1e-12, "imag last");
}
else
{
Assert.AreEqual(0, spectrum[i].Imaginary, 1e-6, "imag");
}
}
}
/// <summary>
/// Naive transforms real sine correctly.
/// </summary>
@ -88,6 +123,20 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
}
}
/// <summary>
/// Radix2XXX when not power of two throws <c>ArgumentException</c>.
/// </summary>
[Test]
public void Radix2ThrowsWhenNotPowerOfTwo32()
{
var samples = Generate.RandomComplex32(0x7F, GetUniform(1));
Assert.Throws(typeof(ArgumentException), () => Fourier.Radix2Forward(samples, FourierOptions.Default));
Assert.Throws(typeof(ArgumentException), () => Fourier.Radix2Inverse(samples, FourierOptions.Default));
Assert.Throws(typeof(ArgumentException), () => Fourier.Radix2(samples, -1));
Assert.Throws(typeof(ArgumentException), () => Fourier.Radix2Parallel(samples, -1));
}
/// <summary>
/// Radix2XXX when not power of two throws <c>ArgumentException</c>.
/// </summary>

99
src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs

@ -52,6 +52,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
return new ContinuousUniform(-1, 1, new System.Random(seed));
}
/// <summary>
/// Fourier naive is reversible.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
public void FourierNaiveIsReversible32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x80, GetUniform(1));
var work = new Complex32[samples.Length];
samples.CopyTo(work, 0);
work = Fourier.NaiveForward(work, options);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
work = Fourier.NaiveInverse(work, options);
AssertHelpers.AlmostEqual(samples, work, 11);
}
/// <summary>
/// Fourier naive is reversible.
/// </summary>
@ -71,6 +90,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(samples, work, 12);
}
/// <summary>
/// Fourier radix2xx is reversible.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
public void FourierRadix2IsReversible32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x8000, GetUniform(1));
var work = new Complex32[samples.Length];
samples.CopyTo(work, 0);
Fourier.Radix2Forward(work, options);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
Fourier.Radix2Inverse(work, options);
AssertHelpers.AlmostEqual(samples, work, 12);
}
/// <summary>
/// Fourier radix2xx is reversible.
/// </summary>
@ -90,6 +128,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(samples, work, 12);
}
/// <summary>
/// Fourier bluestein is reversible.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
public void FourierBluesteinIsReversible32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x7FFF, GetUniform(1));
var work = new Complex32[samples.Length];
samples.CopyTo(work, 0);
Fourier.BluesteinForward(work, options);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
Fourier.BluesteinInverse(work, options);
AssertHelpers.AlmostEqual(samples, work, 10);
}
/// <summary>
/// Fourier bluestein is reversible.
/// </summary>
@ -109,6 +166,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(samples, work, 10);
}
/// <summary>
/// Fourier bluestein is reversible.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
public void FourierRealIsReversible32(FourierOptions options)
{
var samples = Generate.RandomSingle(0x7FFF, GetUniform(1));
var work = new float[samples.Length + 2];
samples.CopyTo(work, 0);
Fourier.ForwardReal(work, samples.Length, options);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
Fourier.InverseReal(work, samples.Length, options);
AssertHelpers.AlmostEqual(samples, work, 5);
}
/// <summary>
/// Fourier bluestein is reversible.
/// </summary>
@ -147,6 +223,29 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(samples, work, 12);
}
/// <summary>
/// Fourier default transform is reversible.
/// </summary>
[Test]
public void FourierDefaultTransformIsReversible32()
{
var samples = Generate.RandomComplex32(0x7FFF, GetUniform(1));
var work = new Complex32[samples.Length];
samples.CopyTo(work, 0);
Fourier.Forward(work);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
Fourier.Inverse(work);
AssertHelpers.AlmostEqual(samples, work, 10);
Fourier.Inverse(work, FourierOptions.Default);
Assert.IsFalse(work.ListAlmostEqual(samples, 6));
Fourier.Forward(work, FourierOptions.Default);
AssertHelpers.AlmostEqual(samples, work, 10);
}
/// <summary>
/// Fourier default transform is reversible.
/// </summary>

222
src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs

@ -54,6 +54,22 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
return new ContinuousUniform(-1, 1, new System.Random(seed));
}
static void Verify(
Complex32[] samples,
int maximumErrorDecimalPlaces,
FourierOptions options,
Func<Complex32[], FourierOptions, Complex32[]> naive,
Action<Complex32[], FourierOptions> fast)
{
var spectrumNaive = naive(samples, options);
var spectrumFast = new Complex32[samples.Length];
samples.CopyTo(spectrumFast, 0);
fast(spectrumFast, options);
AssertHelpers.AlmostEqual(spectrumNaive, spectrumFast, maximumErrorDecimalPlaces);
}
static void Verify(
Complex[] samples,
int maximumErrorDecimalPlaces,
@ -70,6 +86,24 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(spectrumNaive, spectrumFast, maximumErrorDecimalPlaces);
}
static void VerifyInplace(
Complex32[] samples,
int maximumErrorDecimalPlaces,
FourierOptions options,
Action<Complex32[], FourierOptions> expected,
Action<Complex32[], FourierOptions> actual)
{
var spectrumExpected = new Complex32[samples.Length];
samples.CopyTo(spectrumExpected, 0);
expected(spectrumExpected, options);
var spectrumActual = new Complex32[samples.Length];
samples.CopyTo(spectrumActual, 0);
actual(spectrumActual, options);
AssertHelpers.AlmostEqual(spectrumExpected, spectrumActual, maximumErrorDecimalPlaces);
}
static void VerifyInplace(
Complex[] samples,
int maximumErrorDecimalPlaces,
@ -88,6 +122,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
AssertHelpers.AlmostEqual(spectrumExpected, spectrumActual, maximumErrorDecimalPlaces);
}
/// <summary>
/// Fourier Radix2XX matches naive on real sine.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
[TestCase(FourierOptions.NumericalRecipes)]
public void FourierRadix2MatchesNaive_RealSine_32(FourierOptions options)
{
var samples = Generate.PeriodicMap(16, w => new Complex32((float)Math.Sin(w), 0), 16, 1.0, Constants.Pi2);
Verify(samples, 12, options, Fourier.NaiveForward, Fourier.Radix2Forward);
Verify(samples, 12, options, Fourier.NaiveInverse, Fourier.Radix2Inverse);
}
/// <summary>
/// Fourier Radix2XX matches naive on real sine.
/// </summary>
@ -103,6 +152,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 12, options, Fourier.NaiveInverse, Fourier.Radix2Inverse);
}
/// <summary>
/// Fourier Radix2XX matches naive on random.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
[TestCase(FourierOptions.NumericalRecipes)]
public void FourierRadix2MatchesNaive_Random_32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x80, GetUniform(1));
Verify(samples, 10, options, Fourier.NaiveForward, Fourier.Radix2Forward);
Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.Radix2Inverse);
}
/// <summary>
/// Fourier Radix2XX matches naive on random.
/// </summary>
@ -118,6 +182,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.Radix2Inverse);
}
/// <summary>
/// Fourier bluestein matches naive on real sine non-power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
[TestCase(FourierOptions.NumericalRecipes)]
public void FourierBluesteinMatchesNaive_RealSine_Arbitrary_32(FourierOptions options)
{
var samples = Generate.PeriodicMap(14, w => new Complex32((float)Math.Sin(w), 0), 14, 1.0, Constants.Pi2);
Verify(samples, 11, options, Fourier.NaiveForward, Fourier.BluesteinForward);
Verify(samples, 11, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches naive on real sine non-power of two.
/// </summary>
@ -133,6 +212,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 12, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches naive on random power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
[TestCase(FourierOptions.NumericalRecipes)]
public void FourierBluesteinMatchesNaive_Random_PowerOfTwo_32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x80, GetUniform(1));
Verify(samples, 10, options, Fourier.NaiveForward, Fourier.BluesteinForward);
Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches naive on random power of two.
/// </summary>
@ -148,6 +242,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches naive on random non-power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.Matlab)]
[TestCase(FourierOptions.NumericalRecipes)]
public void FourierBluesteinMatchesNaive_Random_Arbitrary_32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x7F, GetUniform(1));
Verify(samples, 9, options, Fourier.NaiveForward, Fourier.BluesteinForward);
Verify(samples, 9, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches naive on random non-power of two.
/// </summary>
@ -163,6 +272,24 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches providers on random power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[TestCase(FourierOptions.Default)]
[TestCase(FourierOptions.NoScaling)]
[TestCase(FourierOptions.AsymmetricScaling)]
[TestCase(FourierOptions.InverseExponent)]
[TestCase(FourierOptions.InverseExponent | FourierOptions.NoScaling)]
[TestCase(FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling)]
public void FourierBluesteinMatchesProvider_Random_Arbitrary_32(FourierOptions options)
{
var samples = Generate.RandomComplex32(0x7F, GetUniform(1));
VerifyInplace(samples, 10, options, Fourier.Forward, Fourier.BluesteinForward);
VerifyInplace(samples, 10, options, Fourier.Inverse, Fourier.BluesteinInverse);
}
/// <summary>
/// Fourier bluestein matches providers on random power of two.
/// </summary>
@ -181,6 +308,33 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
VerifyInplace(samples, 10, options, Fourier.Inverse, Fourier.BluesteinInverse);
}
[TestCase(FourierOptions.Default, 128)]
[TestCase(FourierOptions.Default, 129)]
[TestCase(FourierOptions.NoScaling, 128)]
[TestCase(FourierOptions.NoScaling, 129)]
[TestCase(FourierOptions.AsymmetricScaling, 128)]
[TestCase(FourierOptions.AsymmetricScaling, 129)]
public void RealMatchesComplex32(FourierOptions options, int n)
{
var real = Generate.RandomSingle(n.IsEven() ? n + 2 : n + 1, GetUniform(1));
real[n] = 0f;
if (n.IsEven()) real[n + 1] = 0f;
var complex = new Complex32[n];
for (int i = 0; i < complex.Length; i++)
{
complex[i] = new Complex32(real[i], 0f);
}
Fourier.Forward(complex, options);
Fourier.ForwardReal(real, n, options);
int m = (n + 1) / 2;
for (int i = 0, j = 0; i < m; i++)
{
AssertHelpers.AlmostEqual(complex[i], new Complex32(real[j++], real[j++]), 6);
}
}
[TestCase(FourierOptions.Default, 128)]
[TestCase(FourierOptions.Default, 129)]
[TestCase(FourierOptions.NoScaling, 128)]
@ -208,6 +362,16 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
}
}
[Test]
public void AlgorithmsMatchProvider_PowerOfTwo_Large_32()
{
// 65536 = 2^16
var samples = Generate.RandomComplex32(65536, GetUniform(1));
VerifyInplace(samples, 10, FourierOptions.NoScaling, (s, o) => Control.FourierTransformProvider.Forward(s, FourierTransformScaling.NoScaling), Fourier.Radix2Forward);
VerifyInplace(samples, 10, FourierOptions.NoScaling, (s, o) => Control.FourierTransformProvider.Forward(s, FourierTransformScaling.NoScaling), Fourier.BluesteinForward);
}
[Test]
public void AlgorithmsMatchProvider_PowerOfTwo_Large()
{
@ -218,6 +382,20 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
VerifyInplace(samples, 10, FourierOptions.NoScaling, (s, o) => Control.FourierTransformProvider.Forward(s, FourierTransformScaling.NoScaling), Fourier.BluesteinForward);
}
[Test]
public void AlgorithmsMatchProvider_Arbitrary_Large_32()
{
// 30870 = 2*3*3*5*7*7*7
const FourierOptions options = FourierOptions.NoScaling;
var samples = Generate.RandomComplex32(30870, GetUniform(1));
var provider = new Complex32[samples.Length];
samples.Copy(provider);
Control.FourierTransformProvider.Forward(provider, FourierTransformScaling.NoScaling);
Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward);
}
[Test]
public void AlgorithmsMatchProvider_Arbitrary_Large()
{
@ -232,6 +410,19 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward);
}
[Test]
public void AlgorithmsMatchProvider_Arbitrary_Large_GH286_32()
{
const FourierOptions options = FourierOptions.NoScaling;
var samples = Generate.RandomComplex32(46500, GetUniform(1));
var provider = new Complex32[samples.Length];
samples.Copy(provider);
Control.FourierTransformProvider.Forward(provider, FourierTransformScaling.NoScaling);
Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward);
}
[Test]
public void AlgorithmsMatchProvider_Arbitrary_Large_GH286()
{
@ -245,6 +436,18 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, (a, b) => provider, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_PowerOfTwo_Large_32()
{
// 65536 = 2^16
const FourierOptions options = FourierOptions.NoScaling;
var samples = Generate.RandomComplex32(65536, GetUniform(1));
var naive = Fourier.NaiveForward(samples, options);
Verify(samples, 3, options, (a, b) => naive, Fourier.Radix2Forward);
Verify(samples, 3, options, (a, b) => naive, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_PowerOfTwo_Large()
{
@ -257,6 +460,16 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, (a, b) => naive, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_Arbitrary_Large_32()
{
// 30870 = 2*3*3*5*7*7*7
const FourierOptions options = FourierOptions.NoScaling;
var samples = Generate.RandomComplex32(30870, GetUniform(1));
Verify(samples, 4, options, Fourier.NaiveForward, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_Arbitrary_Large()
{
@ -268,6 +481,15 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Verify(samples, 10, options, (a, b) => naive, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_Arbitrary_Large_GH286_32()
{
const FourierOptions options = FourierOptions.NoScaling;
var samples = Generate.RandomComplex32(46500, GetUniform(1));
Verify(samples, 4, options, Fourier.NaiveForward, Fourier.BluesteinForward);
}
[Test, Explicit("Long-Running")]
public void AlgorithmsMatchNaive_Arbitrary_Large_GH286()
{

23
src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs

@ -54,6 +54,29 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
return new ContinuousUniform(-1, 1, new System.Random(seed));
}
/// <summary>
/// Fourier default transform satisfies Parseval's theorem.
/// </summary>
/// <param name="count">Samples count.</param>
[TestCase(0x1000)]
[TestCase(0x7FF)]
public void FourierDefaultTransformSatisfiesParsevalsTheorem32(int count)
{
var samples = Generate.RandomComplex32(count, GetUniform(1));
var timeSpaceEnergy = (from s in samples select s.MagnitudeSquared()).Mean();
var work = new Complex32[samples.Length];
samples.CopyTo(work, 0);
// Default -> Symmetric Scaling
Fourier.Forward(work);
var frequencySpaceEnergy = (from s in work select s.MagnitudeSquared()).Mean();
Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-7);
}
/// <summary>
/// Fourier default transform satisfies Parseval's theorem.
/// </summary>

Loading…
Cancel
Save