diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs index c9703e4b..a7832c85 100644 --- a/src/Numerics/IntegralTransforms/Fourier.cs +++ b/src/Numerics/IntegralTransforms/Fourier.cs @@ -28,6 +28,8 @@ // using System; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Storage; using MathNet.Numerics.Providers.FourierTransform; namespace MathNet.Numerics.IntegralTransforms @@ -50,6 +52,21 @@ namespace MathNet.Numerics.IntegralTransforms Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); } + public static void ForwardMultiDim(Complex[] samples, int[] dimensions) + { + Control.FourierTransformProvider.ForwardInplaceMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling); + } + + public static void Forward2D(Complex[] samplesRowWise, int rows, int columns) + { + ForwardMultiDim(samplesRowWise, new[] {rows, columns}); + } + + public static void Forward2D(Matrix samples) + { + Forward2D(samples, FourierOptions.Default); + } + /// /// Applies the forward Fast Fourier Transform (FFT) to arbitrary-length sample vectors. /// @@ -76,54 +93,159 @@ namespace MathNet.Numerics.IntegralTransforms } } + public static void ForwardMultiDim(Complex[] samples, int[] dimensions, FourierOptions options) + { + switch (options) + { + case FourierOptions.NoScaling: + case FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.ForwardInplaceMultidim(samples, dimensions, FourierTransformScaling.NoScaling); + break; + case FourierOptions.InverseExponent: + Control.FourierTransformProvider.BackwardInplaceMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.NoScaling: + case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.BackwardInplaceMultidim(samples, dimensions, FourierTransformScaling.NoScaling); + break; + default: + Control.FourierTransformProvider.ForwardInplaceMultidim(samples, dimensions, FourierTransformScaling.SymmetricScaling); + break; + } + } + + public static void Forward2D(Complex[] samplesRowWise, int rows, int columns, FourierOptions options) + { + ForwardMultiDim(samplesRowWise, new[] { rows, columns }, options); + } + + public static void Forward2D(Matrix samples, FourierOptions options) + { + // since dense matrix data is column major, we switch rows and columns + + var denseStorage = samples.Storage as DenseColumnMajorMatrixStorage; + if (denseStorage == null) + { + var samplesColumnMajor = samples.ToColumnWiseArray(); + ForwardMultiDim(samplesColumnMajor, new[] { samples.ColumnCount, samples.RowCount }, options); + denseStorage = new DenseColumnMajorMatrixStorage(samples.RowCount, samples.ColumnCount, samplesColumnMajor); + denseStorage.CopyToUnchecked(samples.Storage, ExistingData.Clear); + return; + } + + ForwardMultiDim(denseStorage.Data, new[] { samples.ColumnCount, samples.RowCount }, options); + } + /// /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. /// - /// Sample vector, where the FFT is evaluated in place. - public static void Inverse(Complex[] samples) + /// Sample vector, where the FFT is evaluated in place. + public static void Inverse(Complex[] spectrum) + { + Control.FourierTransformProvider.BackwardInplace(spectrum, FourierTransformScaling.SymmetricScaling); + } + + public static void InverseMultiDim(Complex[] spectrum, int[] dimensions) + { + Control.FourierTransformProvider.BackwardInplaceMultidim(spectrum, dimensions, FourierTransformScaling.SymmetricScaling); + } + + public static void Inverse2D(Complex[] spectrumRowWise, int rows, int columns) + { + InverseMultiDim(spectrumRowWise, new[] { rows, columns }); + } + + public static void Inverse2D(Matrix samples) { - Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.SymmetricScaling); + Inverse2D(samples, FourierOptions.Default); } /// /// Applies the inverse Fast Fourier Transform (iFFT) to arbitrary-length sample vectors. /// - /// Sample vector, where the FFT is evaluated in place. + /// Sample vector, where the FFT is evaluated in place. /// Fourier Transform Convention Options. - public static void Inverse(Complex[] samples, FourierOptions options) + public static void Inverse(Complex[] specrum, FourierOptions options) { switch (options) { case FourierOptions.NoScaling: - Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.NoScaling); + Control.FourierTransformProvider.BackwardInplace(specrum, FourierTransformScaling.NoScaling); break; case FourierOptions.AsymmetricScaling: - Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.BackwardScaling); + Control.FourierTransformProvider.BackwardInplace(specrum, FourierTransformScaling.BackwardScaling); break; case FourierOptions.InverseExponent: - Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); + Control.FourierTransformProvider.ForwardInplace(specrum, FourierTransformScaling.SymmetricScaling); break; case FourierOptions.InverseExponent | FourierOptions.NoScaling: - Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.NoScaling); + Control.FourierTransformProvider.ForwardInplace(specrum, FourierTransformScaling.NoScaling); break; case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling: - Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.ForwardScaling); + Control.FourierTransformProvider.ForwardInplace(specrum, FourierTransformScaling.ForwardScaling); break; default: - Control.FourierTransformProvider.BackwardInplace(samples, FourierTransformScaling.SymmetricScaling); + Control.FourierTransformProvider.BackwardInplace(specrum, FourierTransformScaling.SymmetricScaling); break; } } + public static void InverseMultiDim(Complex[] specrum, int[] dimensions, FourierOptions options) + { + switch (options) + { + case FourierOptions.NoScaling: + Control.FourierTransformProvider.BackwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.NoScaling); + break; + case FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.BackwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.BackwardScaling); + break; + case FourierOptions.InverseExponent: + Control.FourierTransformProvider.ForwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.SymmetricScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.NoScaling: + Control.FourierTransformProvider.ForwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.NoScaling); + break; + case FourierOptions.InverseExponent | FourierOptions.AsymmetricScaling: + Control.FourierTransformProvider.ForwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.ForwardScaling); + break; + default: + Control.FourierTransformProvider.BackwardInplaceMultidim(specrum, dimensions, FourierTransformScaling.SymmetricScaling); + break; + } + } + + public static void Inverse2D(Complex[] spectrumRowWise, int rows, int columns, FourierOptions options) + { + InverseMultiDim(spectrumRowWise, new[] { rows, columns }, options); + } + + public static void Inverse2D(Matrix spectrum, FourierOptions options) + { + // since dense matrix data is column major, we switch rows and columns + + var denseStorage = spectrum.Storage as DenseColumnMajorMatrixStorage; + if (denseStorage == null) + { + var samplesColumnMajor = spectrum.ToColumnWiseArray(); + InverseMultiDim(samplesColumnMajor, new[] { spectrum.ColumnCount, spectrum.RowCount }, options); + denseStorage = new DenseColumnMajorMatrixStorage(spectrum.RowCount, spectrum.ColumnCount, samplesColumnMajor); + denseStorage.CopyToUnchecked(spectrum.Storage, ExistingData.Clear); + return; + } + + InverseMultiDim(denseStorage.Data, new[] { spectrum.ColumnCount, spectrum.RowCount }, options); + } + /// /// Naive forward DFT, useful e.g. to verify faster algorithms. /// - /// Time-space sample vector. + /// Time-space sample vector. /// Fourier Transform Convention Options. /// Corresponding frequency-space vector. - public static Complex[] NaiveForward(Complex[] timeSpace, FourierOptions options) + public static Complex[] NaiveForward(Complex[] samples, FourierOptions options) { - var frequencySpace = Naive(timeSpace, SignByOptions(options)); + var frequencySpace = Naive(samples, SignByOptions(options)); ForwardScaleByOptions(options, frequencySpace); return frequencySpace; } @@ -131,12 +253,12 @@ namespace MathNet.Numerics.IntegralTransforms /// /// Naive inverse DFT, useful e.g. to verify faster algorithms. /// - /// Frequency-space sample vector. + /// Frequency-space sample vector. /// Fourier Transform Convention Options. /// Corresponding time-space vector. - public static Complex[] NaiveInverse(Complex[] frequencySpace, FourierOptions options) + public static Complex[] NaiveInverse(Complex[] spectrum, FourierOptions options) { - var timeSpace = Naive(frequencySpace, -SignByOptions(options)); + var timeSpace = Naive(spectrum, -SignByOptions(options)); InverseScaleByOptions(options, timeSpace); return timeSpace; } @@ -156,13 +278,13 @@ namespace MathNet.Numerics.IntegralTransforms /// /// Radix-2 inverse FFT for power-of-two sized sample vectors. /// - /// Sample vector, where the FFT is evaluated in place. + /// Sample vector, where the FFT is evaluated in place. /// Fourier Transform Convention Options. /// - public static void Radix2Inverse(Complex[] samples, FourierOptions options) + public static void Radix2Inverse(Complex[] spectrum, FourierOptions options) { - Radix2Parallel(samples, -SignByOptions(options)); - InverseScaleByOptions(options, samples); + Radix2Parallel(spectrum, -SignByOptions(options)); + InverseScaleByOptions(options, spectrum); } /// @@ -179,12 +301,12 @@ namespace MathNet.Numerics.IntegralTransforms /// /// Bluestein inverse FFT for arbitrary sized sample vectors. /// - /// Sample vector, where the FFT is evaluated in place. + /// Sample vector, where the FFT is evaluated in place. /// Fourier Transform Convention Options. - public static void BluesteinInverse(Complex[] samples, FourierOptions options) + public static void BluesteinInverse(Complex[] spectrum, FourierOptions options) { - Bluestein(samples, -SignByOptions(options)); - InverseScaleByOptions(options, samples); + Bluestein(spectrum, -SignByOptions(options)); + InverseScaleByOptions(options, spectrum); } /// diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index fe21d002..01fe4844 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -270,12 +270,16 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl public void ForwardInplaceMultidim(Complex[] complex, int[] dimensions, FourierTransformScaling scaling) { - throw new NotImplementedException(); + Kernel kernel = Configure(dimensions, scaling); + SafeNativeMethods.z_fft_forward(kernel.Handle, complex); + Release(kernel); } public void BackwardInplaceMultidim(Complex[] complex, int[] dimensions, FourierTransformScaling scaling) { - throw new NotImplementedException(); + Kernel kernel = Configure(dimensions, scaling); + SafeNativeMethods.z_fft_backward(kernel.Handle, complex); + Release(kernel); } static double ForwardScaling(FourierTransformScaling scaling, long length)