diff --git a/src/NativeProviders/MKL/fft.cpp b/src/NativeProviders/MKL/fft.cpp index c35985b6..bc92dfb4 100644 --- a/src/NativeProviders/MKL/fft.cpp +++ b/src/NativeProviders/MKL/fft.cpp @@ -7,14 +7,20 @@ #include "mkl_service.h" #include "mkl_dfti.h" -template -inline MKL_LONG fft_inplace(MKL_LONG n, Data x[], DFTI_CONFIG_VALUE precision, DFTI_CONFIG_VALUE domain, FFT fft) +template +inline MKL_LONG fft_1d_inplace(const MKL_LONG n, Data x[], const Precision forward_scale, const Precision backward_scale, const DFTI_CONFIG_VALUE precision, const DFTI_CONFIG_VALUE domain, FFT fft) { MKL_LONG status = 0; DFTI_DESCRIPTOR_HANDLE descriptor = 0; status = DftiCreateDescriptor(&descriptor, precision, domain, 1, n); if (0 != status) goto failed; + status = DftiSetValue(descriptor, DFTI_FORWARD_SCALE, forward_scale); + if (0 != status) goto failed; + + status = DftiSetValue(descriptor, DFTI_BACKWARD_SCALE, backward_scale); + if (0 != status) goto failed; + status = DftiCommitDescriptor(descriptor); if (0 != status) goto failed; @@ -32,23 +38,23 @@ failed: extern "C" { - DLLEXPORT MKL_LONG z_fft_forward_inplace(MKL_LONG n, MKL_Complex16 x[]) + DLLEXPORT MKL_LONG z_fft_forward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) { - return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); + return fft_1d_inplace(n, x, scaling, 1.0, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG c_fft_forward_inplace(MKL_LONG n, MKL_Complex8 x[]) + DLLEXPORT MKL_LONG c_fft_forward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) { - return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); + return fft_1d_inplace(n, x, scaling, 1.0f, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeForward); } - DLLEXPORT MKL_LONG z_fft_backward_inplace(MKL_LONG n, MKL_Complex16 x[]) + DLLEXPORT MKL_LONG z_fft_backward_inplace(const MKL_LONG n, const double scaling, MKL_Complex16 x[]) { - return fft_inplace(n, x, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); + return fft_1d_inplace(n, x, 1.0, scaling, DFTI_DOUBLE, DFTI_COMPLEX, DftiComputeBackward); } - DLLEXPORT MKL_LONG c_fft_backward_inplace(MKL_LONG n, MKL_Complex8 x[]) + DLLEXPORT MKL_LONG c_fft_backward_inplace(const MKL_LONG n, const float scaling, MKL_Complex8 x[]) { - return fft_inplace(n, x, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); + return fft_1d_inplace(n, x, 1.0f, scaling, DFTI_SINGLE, DFTI_COMPLEX, DftiComputeBackward); } } diff --git a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs index ba52fc2f..2e4dea88 100644 --- a/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/IFourierTransformProvider.cs @@ -34,6 +34,13 @@ namespace MathNet.Numerics.Providers.FourierTransform using Complex = System.Numerics.Complex; #endif + public enum FourierTransformScaling : int + { + NoScaling = 0, + SymmetricScaling = 1, + AsymmetricScaling = 2 + } + public interface IFourierTransformProvider { /// @@ -41,10 +48,10 @@ namespace MathNet.Numerics.Providers.FourierTransform /// void InitializeVerify(); - void ForwardInplace(Complex[] complex); - void BackwardInplace(Complex[] complex); + void ForwardInplace(Complex[] complex, FourierTransformScaling scaling); + void BackwardInplace(Complex[] complex, FourierTransformScaling scaling); - Complex[] Forward(Complex[] complexTimeSpace); - Complex[] Backward(Complex[] complexFrequenceSpace); + Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling); + Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling); } } diff --git a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs index 14a37198..696f8f95 100644 --- a/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs @@ -37,30 +37,44 @@ namespace MathNet.Numerics.Providers.FourierTransform { } - public virtual void ForwardInplace(Complex[] complex) + public virtual void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinForward(complex, FourierOptions.Default); + Fourier.BluesteinForward(complex, Options(scaling)); } - public virtual void BackwardInplace(Complex[] complex) + public virtual void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) { - Fourier.BluesteinInverse(complex, FourierOptions.Default); + Fourier.BluesteinInverse(complex, Options(scaling)); } - public virtual Complex[] Forward(Complex[] complexTimeSpace) + public virtual Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexTimeSpace.Length]; complexTimeSpace.Copy(work); - ForwardInplace(work); + ForwardInplace(work, scaling); return work; } - public virtual Complex[] Backward(Complex[] complexFrequenceSpace) + public virtual Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexFrequenceSpace.Length]; complexFrequenceSpace.Copy(work); - BackwardInplace(work); + BackwardInplace(work, scaling); return work; } + + private FourierOptions Options(FourierTransformScaling scaling) + { + switch (scaling) + { + case FourierTransformScaling.NoScaling: + return FourierOptions.NoScaling; + case FourierTransformScaling.AsymmetricScaling: + return FourierOptions.AsymmetricScaling; + case FourierTransformScaling.SymmetricScaling: + default: + return FourierOptions.Default; + } + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs index 49f57311..c3c66484 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/MklFourierTransformProvider.cs @@ -26,6 +26,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using System.Numerics; namespace MathNet.Numerics.Providers.FourierTransform.Mkl @@ -36,30 +37,54 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl { } - public void ForwardInplace(Complex[] complex) + public void ForwardInplace(Complex[] complex, FourierTransformScaling scaling) { - SafeNativeMethods.z_fft_forward_inplace(complex.Length, complex); + SafeNativeMethods.z_fft_forward_inplace(complex.Length, ForwardScaling(scaling, complex.Length), complex); } - public void BackwardInplace(Complex[] complex) + public void BackwardInplace(Complex[] complex, FourierTransformScaling scaling) { - SafeNativeMethods.z_fft_backward_inplace(complex.Length, complex); + SafeNativeMethods.z_fft_backward_inplace(complex.Length, BackwardScaling(scaling, complex.Length), complex); } - public Complex[] Forward(Complex[] complexTimeSpace) + public Complex[] Forward(Complex[] complexTimeSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexTimeSpace.Length]; complexTimeSpace.Copy(work); - ForwardInplace(work); + ForwardInplace(work, scaling); return work; } - public Complex[] Backward(Complex[] complexFrequenceSpace) + public Complex[] Backward(Complex[] complexFrequenceSpace, FourierTransformScaling scaling) { Complex[] work = new Complex[complexFrequenceSpace.Length]; complexFrequenceSpace.Copy(work); - BackwardInplace(work); + BackwardInplace(work, scaling); return work; } + + private double ForwardScaling(FourierTransformScaling scaling, int length) + { + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + return Math.Sqrt(1.0/length); + default: + return 1.0; + } + } + + private double BackwardScaling(FourierTransformScaling scaling, int length) + { + switch (scaling) + { + case FourierTransformScaling.SymmetricScaling: + return Math.Sqrt(1.0/length); + case FourierTransformScaling.AsymmetricScaling: + return 1.0/length; + default: + return 1.0; + } + } } } diff --git a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs index 27cfcb5a..02be1430 100644 --- a/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs +++ b/src/Numerics/Providers/FourierTransform/Mkl/SafeNativeMethods.cs @@ -82,16 +82,16 @@ namespace MathNet.Numerics.Providers.FourierTransform.Mkl #region FFT [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long z_fft_forward_inplace(long n, [In, Out] Complex[] x); + internal static extern long z_fft_forward_inplace(long n, double scaling, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_forward_inplace(long n, [In, Out] Complex32[] x); + internal static extern long c_fft_forward_inplace(long n, float scaling, [In, Out] Complex32[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long z_fft_backward_inplace(long n, [In, Out] Complex[] x); + internal static extern long z_fft_backward_inplace(long n, double scaling, [In, Out] Complex[] x); [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern long c_fft_backward_inplace(long n, [In, Out] Complex32[] x); + internal static extern long c_fft_backward_inplace(long n, float scaling, [In, Out] Complex32[] x); #endregion FFT diff --git a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs index 588768b5..91e67507 100644 --- a/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs +++ b/src/UnitTests/FourierTransformProviderTests/FourierTransformProviderTests.cs @@ -29,6 +29,7 @@ using System; using MathNet.Numerics.Distributions; +using MathNet.Numerics.Providers.FourierTransform; using MathNet.Numerics.Statistics; using NUnit.Framework; @@ -54,7 +55,7 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests // real-odd transforms to imaginary odd samples.Copy(spectrum); - Control.FourierTransformProvider.ForwardInplace(spectrum); + Control.FourierTransformProvider.ForwardInplace(spectrum, FourierTransformScaling.AsymmetricScaling); // all real components must be zero foreach (var c in spectrum) @@ -87,15 +88,8 @@ namespace MathNet.Numerics.UnitTests.FourierTransformProviderTests var samples = Generate.RandomComplex(count, GetUniform(1)); var timeSpaceEnergy = Generate.Map(samples, s => s.MagnitudeSquared()).Mean(); - var work = new Complex[samples.Length]; - samples.Copy(work); - - Control.FourierTransformProvider.ForwardInplace(work); - - var frequencySpaceEnergy = Generate.Map(work, s => s.MagnitudeSquared()).Mean(); - - // TODO: normalize scaling - this should instead be controllable, not needed by default - frequencySpaceEnergy /= count; + Control.FourierTransformProvider.ForwardInplace(samples, FourierTransformScaling.SymmetricScaling); + var frequencySpaceEnergy = Generate.Map(samples, s => s.MagnitudeSquared()).Mean(); Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); }