diff --git a/src/Numerics/ComplexExtensions.cs b/src/Numerics/ComplexExtensions.cs index 78d46b92..532ea9dd 100644 --- a/src/Numerics/ComplexExtensions.cs +++ b/src/Numerics/ComplexExtensions.cs @@ -45,6 +45,16 @@ namespace MathNet.Numerics /// public static class ComplexExtensions { + /// + /// Gets the squared magnitude of the Complex number. + /// + /// The number to perfom this operation on. + /// The squared magnitude of the Complex number. + public static double MagnitudeSquared(this Complex32 complex) + { + return (complex.Real * complex.Real) + (complex.Imaginary * complex.Imaginary); + } + /// /// Gets the squared magnitude of the Complex number. /// diff --git a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs index 820c88e3..46c38d30 100644 --- a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs +++ b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs @@ -47,6 +47,38 @@ namespace MathNet.Numerics.IntegralTransforms /// const int BluesteinSequenceLengthThreshold = 46341; + /// + /// Generate the bluestein sequence for the provided problem size. + /// + /// Number of samples. + /// Bluestein sequence exp(I*Pi*k^2/N) + 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; + } + /// /// Generate the bluestein sequence for the provided problem size. /// @@ -79,6 +111,61 @@ namespace MathNet.Numerics.IntegralTransforms return sequence; } + /// + /// Convolution with the bluestein sequence (Parallel Version). + /// + /// Sample Vector. + 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]; + } + } + /// /// Convolution with the bluestein sequence (Parallel Version). /// @@ -134,6 +221,18 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Swap the real and imaginary parts of each sample. + /// + /// Sample Vector. + static void SwapRealImaginary(Complex32[] samples) + { + for (int i = 0; i < samples.Length; i++) + { + samples[i] = new Complex32(samples[i].Imaginary, samples[i].Real); + } + } + /// /// Swap the real and imaginary parts of each sample. /// @@ -146,6 +245,33 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Bluestein generic FFT for arbitrary sized sample vectors. + /// + /// Time-space sample vector. + /// Fourier series exponent sign. + 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); + } + } + /// /// Bluestein generic FFT for arbitrary sized sample vectors. /// diff --git a/src/Numerics/IntegralTransforms/Fourier.Naive.cs b/src/Numerics/IntegralTransforms/Fourier.Naive.cs index 609abbde..b03ebca7 100644 --- a/src/Numerics/IntegralTransforms/Fourier.Naive.cs +++ b/src/Numerics/IntegralTransforms/Fourier.Naive.cs @@ -41,6 +41,36 @@ namespace MathNet.Numerics.IntegralTransforms /// public static partial class Fourier { + /// + /// Naive generic DFT, useful e.g. to verify faster algorithms. + /// + /// Time-space sample vector. + /// Fourier series exponent sign. + /// Corresponding frequency-space vector. + 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; + } + /// /// Naive generic DFT, useful e.g. to verify faster algorithms. /// diff --git a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs index 98564ee5..0d029bff 100644 --- a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs +++ b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs @@ -70,6 +70,29 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Radix-2 Step Helper Method + /// + /// Sample vector. + /// Fourier series exponent sign. + /// Level Group Size. + /// Index inside of the level. + 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; + } + } + /// /// Radix-2 Step Helper Method /// @@ -93,6 +116,29 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Radix-2 generic FFT for power-of-two sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier series exponent sign. + /// + 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); + } + } + } + /// /// Radix-2 generic FFT for power-of-two sized sample vectors. /// @@ -116,6 +162,34 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier series exponent sign. + /// + 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); + } + }); + } + } + /// /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs index dcf8710f..52943aa4 100644 --- a/src/Numerics/IntegralTransforms/Fourier.cs +++ b/src/Numerics/IntegralTransforms/Fourier.cs @@ -44,6 +44,15 @@ namespace MathNet.Numerics.IntegralTransforms /// public static partial class Fourier { + /// + /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + public static void Forward(Complex32[] samples) + { + Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.SymmetricScaling); + } + /// /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. /// @@ -52,6 +61,32 @@ namespace MathNet.Numerics.IntegralTransforms { Control.FourierTransformProvider.Forward(samples, FourierTransformScaling.SymmetricScaling); } + + /// + /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + 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; + } + } /// /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. @@ -78,6 +113,37 @@ namespace MathNet.Numerics.IntegralTransforms break; } } + + /// + /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. + /// + /// Real part of the sample vector, where the FFT is evaluated in place. + /// Imaginary part of the sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + 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; + } + } /// /// 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; } } + + /// + /// 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. + /// + /// Data array of length N+2 (if N is even) or N+1 (if N is odd). + /// The number of samples. + /// Fourier Transform Convention Options. + 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; + } + } /// /// Packed Real-Complex forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. @@ -144,6 +244,36 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data. + /// + /// Sample data, where the FFT is evaluated in place. + /// + /// 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. + /// + /// Fourier Transform Convention Options. + 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; + } + } + /// /// Applies the forward Fast Fourier Transform (FFT) to multiple dimensional sample data. /// @@ -173,6 +303,19 @@ namespace MathNet.Numerics.IntegralTransforms break; } } + + /// + /// Applies the forward Fast Fourier Transform (FFT) to two dimensional sample data. + /// + /// Sample data, organized row by row, where the FFT is evaluated in place + /// The number of rows. + /// The number of columns. + /// Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments. + /// Fourier Transform Convention Options. + public static void Forward2D(Complex32[] samplesRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default) + { + ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options); + } /// /// 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); } + + /// + /// Applies the forward Fast Fourier Transform (FFT) to a two dimensional data in form of a matrix. + /// + /// Sample matrix, where the FFT is evaluated in place + /// Fourier Transform Convention Options. + public static void Forward2D(Matrix 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(samples.RowCount, samples.ColumnCount, columnMajorArray); + denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear); + } /// /// 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); } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. + /// + /// Spectrum data, where the iFFT is evaluated in place. + public static void Inverse(Complex32[] spectrum) + { + Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.SymmetricScaling); + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. /// @@ -224,6 +404,36 @@ namespace MathNet.Numerics.IntegralTransforms Control.FourierTransformProvider.Backward(spectrum, FourierTransformScaling.SymmetricScaling); } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. + /// + /// Spectrum data, where the iFFT is evaluated in place. + /// Fourier Transform Convention Options. + 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; + } + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. /// @@ -254,6 +464,37 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. + /// + /// Real part of the sample vector, where the iFFT is evaluated in place. + /// Imaginary part of the sample vector, where the iFFT is evaluated in place. + /// Fourier Transform Convention Options. + 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; + } + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. /// @@ -285,6 +526,42 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// 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. + /// + /// Data array of length N+2 (if N is even) or N+1 (if N is odd). + /// The number of samples. + /// Fourier Transform Convention Options. + 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; + } + } + /// /// 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 } } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data. + /// + /// Spectrum data, where the iFFT is evaluated in place. + /// + /// 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. + /// + /// Fourier Transform Convention Options. + 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; + } + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to multiple dimensional sample data. /// @@ -355,6 +666,19 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data. + /// + /// Sample data, organized row by row, where the iFFT is evaluated in place + /// The number of rows. + /// The number of columns. + /// Data available organized column by column instead of row by row can be processed directly by swapping the rows and columns arguments. + /// Fourier Transform Convention Options. + public static void Inverse2D(Complex32[] spectrumRowWise, int rows, int columns, FourierOptions options = FourierOptions.Default) + { + InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options); + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to two dimensional sample data. /// @@ -368,6 +692,34 @@ namespace MathNet.Numerics.IntegralTransforms InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options); } + /// + /// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix. + /// + /// Sample matrix, where the iFFT is evaluated in place + /// Fourier Transform Convention Options. + public static void Inverse2D(Matrix 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(spectrum.RowCount, spectrum.ColumnCount, columnMajorArray); + denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear); + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to a two dimensional data in form of a matrix. /// @@ -396,6 +748,19 @@ namespace MathNet.Numerics.IntegralTransforms denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear); } + /// + /// Naive forward DFT, useful e.g. to verify faster algorithms. + /// + /// Time-space sample vector. + /// Fourier Transform Convention Options. + /// Corresponding frequency-space vector. + public static Complex32[] NaiveForward(Complex32[] samples, FourierOptions options = FourierOptions.Default) + { + var frequencySpace = Naive(samples, SignByOptions(options)); + ForwardScaleByOptions(options, frequencySpace); + return frequencySpace; + } + /// /// Naive forward DFT, useful e.g. to verify faster algorithms. /// @@ -409,6 +774,19 @@ namespace MathNet.Numerics.IntegralTransforms return frequencySpace; } + /// + /// Naive inverse DFT, useful e.g. to verify faster algorithms. + /// + /// Frequency-space sample vector. + /// Fourier Transform Convention Options. + /// Corresponding time-space vector. + public static Complex32[] NaiveInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default) + { + var timeSpace = Naive(spectrum, -SignByOptions(options)); + InverseScaleByOptions(options, timeSpace); + return timeSpace; + } + /// /// Naive inverse DFT, useful e.g. to verify faster algorithms. /// @@ -422,6 +800,18 @@ namespace MathNet.Numerics.IntegralTransforms return timeSpace; } + /// + /// Radix-2 forward FFT for power-of-two sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + /// + public static void Radix2Forward(Complex32[] samples, FourierOptions options = FourierOptions.Default) + { + Radix2Parallel(samples, SignByOptions(options)); + ForwardScaleByOptions(options, samples); + } + /// /// Radix-2 forward FFT for power-of-two sized sample vectors. /// @@ -434,6 +824,18 @@ namespace MathNet.Numerics.IntegralTransforms ForwardScaleByOptions(options, samples); } + /// + /// Radix-2 inverse FFT for power-of-two sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + /// + public static void Radix2Inverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default) + { + Radix2Parallel(spectrum, -SignByOptions(options)); + InverseScaleByOptions(options, spectrum); + } + /// /// Radix-2 inverse FFT for power-of-two sized sample vectors. /// @@ -446,6 +848,17 @@ namespace MathNet.Numerics.IntegralTransforms InverseScaleByOptions(options, spectrum); } + /// + /// Bluestein forward FFT for arbitrary sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + public static void BluesteinForward(Complex32[] samples, FourierOptions options = FourierOptions.Default) + { + Bluestein(samples, SignByOptions(options)); + ForwardScaleByOptions(options, samples); + } + /// /// Bluestein forward FFT for arbitrary sized sample vectors. /// @@ -457,6 +870,17 @@ namespace MathNet.Numerics.IntegralTransforms ForwardScaleByOptions(options, samples); } + /// + /// Bluestein inverse FFT for arbitrary sized sample vectors. + /// + /// Sample vector, where the FFT is evaluated in place. + /// Fourier Transform Convention Options. + public static void BluesteinInverse(Complex32[] spectrum, FourierOptions options = FourierOptions.Default) + { + Bluestein(spectrum, -SignByOptions(options)); + InverseScaleByOptions(options, spectrum); + } + /// /// Bluestein inverse FFT for arbitrary sized sample vectors. /// @@ -479,6 +903,26 @@ namespace MathNet.Numerics.IntegralTransforms return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1; } + /// + /// Rescale FFT-the resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + 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; + } + } + /// /// Rescale FFT-the resulting vector according to the provided convention options. /// @@ -499,6 +943,30 @@ namespace MathNet.Numerics.IntegralTransforms } } + /// + /// Rescale the iFFT-resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + 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; + } + } + /// /// Rescale the iFFT-resulting vector according to the provided convention options. /// diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index 4809a664..86cc9a79 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -55,13 +55,19 @@ namespace MathNet.Numerics.Providers.FourierTransform /// 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); } } diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index f29e36a8..463e138c 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/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(); diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 4566d565..514d99f9 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/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); } diff --git a/src/UnitTests/IntegralTransformsTests/FourierTest.cs b/src/UnitTests/IntegralTransformsTests/FourierTest.cs index 17443794..0c5a8a3d 100644 --- a/src/UnitTests/IntegralTransformsTests/FourierTest.cs +++ b/src/UnitTests/IntegralTransformsTests/FourierTest.cs @@ -53,6 +53,41 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests return new ContinuousUniform(-1, 1, new System.Random(seed)); } + /// + /// Naive transforms real sine correctly. + /// + [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"); + } + } + } + /// /// Naive transforms real sine correctly. /// @@ -88,6 +123,20 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests } } + /// + /// Radix2XXX when not power of two throws ArgumentException. + /// + [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)); + } + /// /// Radix2XXX when not power of two throws ArgumentException. /// diff --git a/src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs b/src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs index 7db77ffc..f1ce9e62 100644 --- a/src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs +++ b/src/UnitTests/IntegralTransformsTests/InverseTransformTest.cs @@ -52,6 +52,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests return new ContinuousUniform(-1, 1, new System.Random(seed)); } + /// + /// Fourier naive is reversible. + /// + /// Fourier options. + [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); + } + /// /// Fourier naive is reversible. /// @@ -71,6 +90,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(samples, work, 12); } + /// + /// Fourier radix2xx is reversible. + /// + /// Fourier options. + [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); + } + /// /// Fourier radix2xx is reversible. /// @@ -90,6 +128,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(samples, work, 12); } + /// + /// Fourier bluestein is reversible. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein is reversible. /// @@ -109,6 +166,25 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(samples, work, 10); } + /// + /// Fourier bluestein is reversible. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein is reversible. /// @@ -147,6 +223,29 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(samples, work, 12); } + /// + /// Fourier default transform is reversible. + /// + [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); + } + /// /// Fourier default transform is reversible. /// diff --git a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs b/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs index 2fc25462..43235569 100644 --- a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs +++ b/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 naive, + Action 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 expected, + Action 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); } + /// + /// Fourier Radix2XX matches naive on real sine. + /// + /// Fourier options. + [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); + } + /// /// Fourier Radix2XX matches naive on real sine. /// @@ -103,6 +152,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 12, options, Fourier.NaiveInverse, Fourier.Radix2Inverse); } + /// + /// Fourier Radix2XX matches naive on random. + /// + /// Fourier options. + [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); + } + /// /// Fourier Radix2XX matches naive on random. /// @@ -118,6 +182,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.Radix2Inverse); } + /// + /// Fourier bluestein matches naive on real sine non-power of two. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein matches naive on real sine non-power of two. /// @@ -133,6 +212,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 12, options, Fourier.NaiveInverse, Fourier.BluesteinInverse); } + /// + /// Fourier bluestein matches naive on random power of two. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein matches naive on random power of two. /// @@ -148,6 +242,21 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse); } + /// + /// Fourier bluestein matches naive on random non-power of two. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein matches naive on random non-power of two. /// @@ -163,6 +272,24 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests Verify(samples, 10, options, Fourier.NaiveInverse, Fourier.BluesteinInverse); } + /// + /// Fourier bluestein matches providers on random power of two. + /// + /// Fourier options. + [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); + } + /// /// Fourier bluestein matches providers on random power of two. /// @@ -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() { diff --git a/src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs b/src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs index de5aaadb..b0f8189e 100644 --- a/src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs +++ b/src/UnitTests/IntegralTransformsTests/ParsevalTheoremTest.cs @@ -54,6 +54,29 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests return new ContinuousUniform(-1, 1, new System.Random(seed)); } + /// + /// Fourier default transform satisfies Parseval's theorem. + /// + /// Samples count. + [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); + } + /// /// Fourier default transform satisfies Parseval's theorem. ///