FFT: forward and backward scaling support, drop manual correction in parseval test

src/NativeProviders/MKL/fft.cpp

@@ -7,14 +7,20 @@
 #include "mkl_service.h"
 #include "mkl_dfti.h"
 
-template<typename Data, typename FFT>
-inline MKL_LONG fft_inplace(MKL_LONG n, Data x[], DFTI_CONFIG_VALUE precision, DFTI_CONFIG_VALUE domain, FFT fft)
+template<typename Data, typename Precision, typename FFT>
+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);
 	}
 }

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
     {
         /// <summary>
@@ -41,10 +48,10 @@ namespace MathNet.Numerics.Providers.FourierTransform
         /// </summary>
         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);
     }
 }

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;
+            }
+        }
     }
 }

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

@@ -26,6 +26,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>
 
+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;
+            }
+        }
     }
 }

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
 

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);
         }