diff --git a/src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs b/src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs index 78bd4b59..86535ab0 100644 --- a/src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs +++ b/src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs @@ -40,11 +40,8 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestFixture, Category("FFT")] public class FourierTest { - /// - /// Naive transforms real sine correctly. - /// [Test] - public void NaiveTransformsRealSineCorrectly32() + public void ReferenceDftTransformsRealSineCorrectly32() { var samples = Generate.PeriodicMap(16, w => new Complex32((float)Math.Sin(w), 0), 16, 1.0, Constants.Pi2); @@ -75,11 +72,8 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests } } - /// - /// Naive transforms real sine correctly. - /// [Test] - public void NaiveTransformsRealSineCorrectly() + public void ReferenceDftTransformsRealSineCorrectly64() { var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2); @@ -109,5 +103,69 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests } } } + + [Test] + public void FourierDefaultTransformsRealSineCorrectly32() + { + var samples = Generate.PeriodicMap(16, w => new Complex32((float)Math.Sin(w), 0), 16, 1.0, Constants.Pi2); + + // real-odd transforms to imaginary odd + Fourier.Forward(samples, FourierOptions.Matlab); + + // all real components must be zero + foreach (var c in samples) + { + Assert.AreEqual(0, c.Real, 1e-6, "real"); + } + + // all imaginary components except second and last musth be zero + for (var i = 0; i < samples.Length; i++) + { + if (i == 1) + { + Assert.AreEqual(-8, samples[i].Imaginary, 1e-12, "imag second"); + } + else if (i == samples.Length - 1) + { + Assert.AreEqual(8, samples[i].Imaginary, 1e-12, "imag last"); + } + else + { + Assert.AreEqual(0, samples[i].Imaginary, 1e-6, "imag"); + } + } + } + + [Test] + public void FourierDefaultTransformsRealSineCorrectly64() + { + var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2); + + // real-odd transforms to imaginary odd + Fourier.Forward(samples, FourierOptions.Matlab); + + // all real components must be zero + foreach (var c in samples) + { + Assert.AreEqual(0, c.Real, 1e-12, "real"); + } + + // all imaginary components except second and last musth be zero + for (var i = 0; i < samples.Length; i++) + { + if (i == 1) + { + Assert.AreEqual(-8, samples[i].Imaginary, 1e-12, "imag second"); + } + else if (i == samples.Length - 1) + { + Assert.AreEqual(8, samples[i].Imaginary, 1e-12, "imag last"); + } + else + { + Assert.AreEqual(0, samples[i].Imaginary, 1e-12, "imag"); + } + } + } } } diff --git a/src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs b/src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs index 2b274dea..812d34a4 100644 --- a/src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs +++ b/src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs @@ -54,7 +54,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Fourier options. [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] - public void FourierNaiveIsReversible32(FourierOptions options) + public void ReferenceDftIsReversible32(FourierOptions options) { var samples = Generate.RandomComplex32(0x80, GetUniform(1)); var work = new Complex32[samples.Length]; @@ -73,7 +73,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Fourier options. [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] - public void FourierNaiveIsReversible(FourierOptions options) + public void ReferenceDftIsReversible64(FourierOptions options) { var samples = Generate.RandomComplex(0x80, GetUniform(1)); var work = new Complex[samples.Length]; @@ -111,7 +111,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Fourier options. [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] - public void FourierRadix2IsReversible(FourierOptions options) + public void FourierRadix2IsReversible64(FourierOptions options) { var samples = Generate.RandomComplex(0x8000, GetUniform(1)); var work = new Complex[samples.Length]; @@ -149,7 +149,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Fourier options. [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] - public void FourierBluesteinIsReversible(FourierOptions options) + public void FourierBluesteinIsReversible64(FourierOptions options) { var samples = Generate.RandomComplex(0x7FFF, GetUniform(1)); var work = new Complex[samples.Length]; @@ -187,7 +187,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Fourier options. [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] - public void FourierRealIsReversible(FourierOptions options) + public void FourierRealIsReversible64(FourierOptions options) { var samples = Generate.Random(0x7FFF, GetUniform(1)); var work = new double[samples.Length+2]; @@ -206,7 +206,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Hartley options. [TestCase(HartleyOptions.Default)] [TestCase(HartleyOptions.AsymmetricScaling)] - public void HartleyNaiveIsReversible(HartleyOptions options) + public void HartleyNaiveIsReversible64(HartleyOptions options) { var samples = Generate.Random(0x80, GetUniform(1)); var work = new double[samples.Length]; @@ -218,51 +218,5 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests work = Hartley.NaiveInverse(work, options); 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. - /// - [Test] - public void FourierDefaultTransformIsReversible() - { - var samples = Generate.RandomComplex(0x7FFF, GetUniform(1)); - var work = new Complex[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); - } } } diff --git a/src/Numerics.Tests/IntegralTransformsTests/MatchingNaiveTransformTest.cs b/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs similarity index 68% rename from src/Numerics.Tests/IntegralTransformsTests/MatchingNaiveTransformTest.cs rename to src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs index 665763ea..a114c368 100644 --- a/src/Numerics.Tests/IntegralTransformsTests/MatchingNaiveTransformTest.cs +++ b/src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs @@ -40,7 +40,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Matching Naive transform tests. /// [TestFixture, Category("FFT")] - public class MatchingNaiveTransformTest + public class MatchingReferenceTransformTest { /// /// Continuous uniform distribution. @@ -86,6 +86,42 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests AssertHelpers.AlmostEqual(spectrumExpected, spectrumActual, maximumErrorDecimalPlaces); } + static void Verify( + Complex32[] samples, + int maximumErrorDecimalPlaces, + FourierTransformScaling 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 Verify( + Complex[] samples, + int maximumErrorDecimalPlaces, + FourierTransformScaling options, + Action expected, + Action actual) + { + var spectrumExpected = new Complex[samples.Length]; + samples.CopyTo(spectrumExpected, 0); + expected(spectrumExpected, options); + + var spectrumActual = new Complex[samples.Length]; + samples.CopyTo(spectrumActual, 0); + actual(spectrumActual, options); + + AssertHelpers.AlmostEqual(spectrumExpected, spectrumActual, maximumErrorDecimalPlaces); + } + /// /// Fourier Radix2XX matches naive on real sine. /// @@ -93,10 +129,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierRadix2MatchesNaive_RealSine_32(FourierOptions options) + public void FourierRadix2MatchesReferenceRealSine32(FourierOptions options) { var samples = Generate.PeriodicMap(16, w => new Complex32((float)Math.Sin(w), 0), 16, 1.0, Constants.Pi2); - Verify(samples, 6, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 6, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -108,10 +143,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierRadix2MatchesNaive_RealSine(FourierOptions options) + public void FourierRadix2MatchesReferenceRealSine64(FourierOptions options) { var samples = Generate.PeriodicMap(16, w => new Complex(Math.Sin(w), 0), 16, 1.0, Constants.Pi2); - Verify(samples, 12, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 12, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -123,10 +157,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierRadix2MatchesNaive_Random_32(FourierOptions options) + public void FourierRadix2MatchesReferenceRandom32(FourierOptions options) { var samples = Generate.RandomComplex32(0x80, GetUniform(1)); - Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -138,10 +171,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierRadix2MatchesNaive_Random(FourierOptions options) + public void FourierRadix2MatchesReferenceRandom64(FourierOptions options) { var samples = Generate.RandomComplex(0x80, GetUniform(1)); - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -153,10 +185,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_RealSine_Arbitrary_32(FourierOptions options) + public void FourierBluesteinMatchesReferenceRealSineArbitrary32(FourierOptions options) { var samples = Generate.PeriodicMap(14, w => new Complex32((float)Math.Sin(w), 0), 14, 1.0, Constants.Pi2); - Verify(samples, 6, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 6, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -168,10 +199,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_RealSine_Arbitrary(FourierOptions options) + public void FourierBluesteinMatchesReferenceRealSineArbitrary64(FourierOptions options) { var samples = Generate.PeriodicMap(14, w => new Complex(Math.Sin(w), 0), 14, 1.0, Constants.Pi2); - Verify(samples, 12, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 12, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -183,10 +213,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_Random_PowerOfTwo_32(FourierOptions options) + public void FourierBluesteinMatchesReferenceRandomPowerOfTwo32(FourierOptions options) { var samples = Generate.RandomComplex32(0x80, GetUniform(1)); - Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -198,10 +227,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_Random_PowerOfTwo(FourierOptions options) + public void FourierBluesteinMatchesReferenceRandomPowerOfTwo64(FourierOptions options) { var samples = Generate.RandomComplex(0x80, GetUniform(1)); - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -213,10 +241,9 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_Random_Arbitrary_32(FourierOptions options) + public void FourierBluesteinMatchesReferenceRandomArbitrary32(FourierOptions options) { var samples = Generate.RandomComplex32(0x7F, GetUniform(1)); - Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 5, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } @@ -228,50 +255,13 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.Default)] [TestCase(FourierOptions.Matlab)] [TestCase(FourierOptions.NumericalRecipes)] - public void FourierBluesteinMatchesNaive_Random_Arbitrary(FourierOptions options) + public void FourierBluesteinMatchesReferenceRandomArbitrary64(FourierOptions options) { var samples = Generate.RandomComplex(0x7F, GetUniform(1)); - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Inverse, Fourier.Inverse); } - /// - /// 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)); - - Verify(samples, 5, options, Fourier.Forward, Fourier.Forward); - Verify(samples, 5, options, Fourier.Inverse, Fourier.Inverse); - } - - /// - /// 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(FourierOptions options) - { - var samples = Generate.RandomComplex(0x7F, GetUniform(1)); - - Verify(samples, 10, options, Fourier.Forward, Fourier.Forward); - Verify(samples, 10, options, Fourier.Inverse, Fourier.Inverse); - } - [TestCase(FourierOptions.Default, 128)] [TestCase(FourierOptions.Default, 129)] [TestCase(FourierOptions.NoScaling, 128)] @@ -305,7 +295,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests [TestCase(FourierOptions.NoScaling, 129)] [TestCase(FourierOptions.AsymmetricScaling, 128)] [TestCase(FourierOptions.AsymmetricScaling, 129)] - public void RealMatchesComplex(FourierOptions options, int n) + public void RealMatchesComplex64(FourierOptions options, int n) { var real = Generate.Random(n.IsEven() ? n + 2 : n + 1, GetUniform(1)); real[n] = 0d; @@ -327,119 +317,95 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests } [Test] - public void ProviderMatchesManagedProvider_PowerOfTwo_Large_32() + public void ProviderMatchesManagedProviderPowerOfTwoLarge32() { // 65536 = 2^16 var samples = Generate.RandomComplex32(65536, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 5, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test] - public void ProviderMatchesManagedProvider_PowerOfTwo_Large() + public void ProviderMatchesManagedProviderPowerOfTwoLarge64() { // 65536 = 2^16 var samples = Generate.RandomComplex(65536, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 10, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test] - public void ProviderMatchesManagedProvider_Arbitrary_Large_32() + public void ProviderMatchesManagedProviderArbitraryLarge32() { // 30870 = 2*3*3*5*7*7*7 var samples = Generate.RandomComplex32(30870, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 5, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test] - public void ProviderMatchesManagedProvider_Arbitrary_Large() + public void ProviderMatchesManagedProviderArbitraryLarge64() { // 30870 = 2*3*3*5*7*7*7 var samples = Generate.RandomComplex(30870, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 10, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test] - public void ProviderMatchesManagedProvider_Arbitrary_Large_GH286_32() + public void ProviderMatchesManagedProviderArbitraryLarge32_GH286() { var samples = Generate.RandomComplex32(46500, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 5, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 5, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test] - public void ProviderMatchesManagedProvider_Arbitrary_Large_GH286() + public void ProviderMatchesManagedProviderArbitraryLarge64_GH286() { var samples = Generate.RandomComplex(46500, GetUniform(1)); - var managed = FourierTransformControl.CreateManaged(); - - Verify(samples, 10, FourierOptions.NoScaling, (s, o) => managed.Forward(s, FourierTransformScaling.NoScaling), (s, o) => FourierTransformControl.Provider.Forward(s, FourierTransformScaling.NoScaling)); + Verify(samples, 10, FourierTransformScaling.NoScaling, FourierTransformControl.CreateManaged().Forward, FourierTransformControl.Provider.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_PowerOfTwo_Large_32() + public void AlgorithmsMatchReferencePowerOfTwoLarge32() { // 65536 = 2^16 - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex32(65536, GetUniform(1)); - - Verify(samples, 3, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 3, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_PowerOfTwo_Large() + public void AlgorithmsMatchReferencePowerOfTwoLarge64() { // 65536 = 2^16 - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex(65536, GetUniform(1)); - - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 10, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_Arbitrary_Large_32() + public void AlgorithmsMatchReferenceArbitraryLarge32() { // 30870 = 2*3*3*5*7*7*7 - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex32(30870, GetUniform(1)); - - Verify(samples, 4, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 4, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_Arbitrary_Large() + public void AlgorithmsMatchReferenceArbitraryLarge64() { // 30870 = 2*3*3*5*7*7*7 - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex(30870, GetUniform(1)); - - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 10, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_Arbitrary_Large_GH286_32() + public void AlgorithmsMatchReferenceArbitraryLarge32_GH286() { - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex32(46500, GetUniform(1)); - - Verify(samples, 4, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 4, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } [Test, Explicit("Long-Running")] - public void AlgorithmsMatchNaive_Arbitrary_Large_GH286() + public void AlgorithmsMatchReferenceArbitraryLarge64_GH286() { - const FourierOptions options = FourierOptions.NoScaling; var samples = Generate.RandomComplex(46500, GetUniform(1)); - - Verify(samples, 10, options, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); + Verify(samples, 10, FourierOptions.NoScaling, ReferenceDiscreteFourierTransform.Forward, Fourier.Forward); } } } diff --git a/src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs b/src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs index c61f8185..ebf41cd3 100644 --- a/src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs +++ b/src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs @@ -56,19 +56,15 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Samples count. [TestCase(0x1000)] [TestCase(0x7FF)] - public void FourierDefaultTransformSatisfiesParsevalsTheorem32(int count) + public void ReferenceDftSatisfiesParsevalsTheorem32(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(); + var spectrum = new Complex32[samples.Length]; + samples.CopyTo(spectrum, 0); + ReferenceDiscreteFourierTransform.Forward(spectrum); + var frequencySpaceEnergy = (from s in spectrum select s.MagnitudeSquared()).Mean(); Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-7); } @@ -79,19 +75,53 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Samples count. [TestCase(0x1000)] [TestCase(0x7FF)] - public void FourierDefaultTransformSatisfiesParsevalsTheorem(int count) + public void ReferenceDftSatisfiesParsevalsTheorem64(int count) { var samples = Generate.RandomComplex(count, GetUniform(1)); + var timeSpaceEnergy = (from s in samples select s.MagnitudeSquared()).Mean(); + + var spectrum = new Complex[samples.Length]; + samples.CopyTo(spectrum, 0); + ReferenceDiscreteFourierTransform.Forward(spectrum); + var frequencySpaceEnergy = (from s in spectrum select s.MagnitudeSquared()).Mean(); + + Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); + } + /// + /// 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 Complex[samples.Length]; - samples.CopyTo(work, 0); + var spectrum = new Complex32[samples.Length]; + samples.CopyTo(spectrum, 0); + Fourier.Forward(spectrum); + var frequencySpaceEnergy = (from s in spectrum select s.MagnitudeSquared()).Mean(); + + Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-7); + } - // Default -> Symmetric Scaling - Fourier.Forward(work); + /// + /// Fourier default transform satisfies Parseval's theorem. + /// + /// Samples count. + [TestCase(0x1000)] + [TestCase(0x7FF)] + public void FourierDefaultTransformSatisfiesParsevalsTheorem64(int count) + { + var samples = Generate.RandomComplex(count, GetUniform(1)); + var timeSpaceEnergy = (from s in samples select s.MagnitudeSquared()).Mean(); - var frequencySpaceEnergy = (from s in work select s.MagnitudeSquared()).Mean(); + var spectrum = new Complex[samples.Length]; + samples.CopyTo(spectrum, 0); + Fourier.Forward(spectrum); + var frequencySpaceEnergy = (from s in spectrum select s.MagnitudeSquared()).Mean(); Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); } @@ -102,19 +132,16 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests /// Samples count. [TestCase(0x40)] [TestCase(0x1F)] - public void HartleyDefaultNaiveSatisfiesParsevalsTheorem(int count) + public void HartleyDefaultNaiveSatisfiesParsevalsTheorem64(int count) { var samples = Generate.Random(count, GetUniform(1)); - var timeSpaceEnergy = (from s in samples select s*s).Mean(); - var work = new double[samples.Length]; - samples.CopyTo(work, 0); - - // Default -> Symmetric Scaling - work = Hartley.NaiveForward(work, HartleyOptions.Default); + var spectrum = new double[samples.Length]; + samples.CopyTo(spectrum, 0); + spectrum = Hartley.NaiveForward(spectrum, HartleyOptions.Default); - var frequencySpaceEnergy = (from s in work select s*s).Mean(); + var frequencySpaceEnergy = (from s in spectrum select s*s).Mean(); Assert.AreEqual(timeSpaceEnergy, frequencySpaceEnergy, 1e-12); } } diff --git a/src/Numerics/IntegralTransforms/Fourier.Naive.cs b/src/Numerics/IntegralTransforms/Fourier.Naive.cs deleted file mode 100644 index 28490d7a..00000000 --- a/src/Numerics/IntegralTransforms/Fourier.Naive.cs +++ /dev/null @@ -1,101 +0,0 @@ -// -// Math.NET Numerics, part of the Math.NET Project -// http://numerics.mathdotnet.com -// http://github.com/mathnet/mathnet-numerics -// -// Copyright (c) 2009-2014 Math.NET -// -// Permission is hereby granted, free of charge, to any person -// obtaining a copy of this software and associated documentation -// files (the "Software"), to deal in the Software without -// restriction, including without limitation the rights to use, -// copy, modify, merge, publish, distribute, sublicense, and/or sell -// copies of the Software, and to permit persons to whom the -// Software is furnished to do so, subject to the following -// conditions: -// -// The above copyright notice and this permission notice shall be -// included in all copies or substantial portions of the Software. -// -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -// OTHER DEALINGS IN THE SOFTWARE. -// - -using System; -using System.Numerics; -using MathNet.Numerics.Threading; - -namespace MathNet.Numerics.IntegralTransforms -{ - /// - /// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT). - /// - 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. - /// - /// Time-space sample vector. - /// Fourier series exponent sign. - /// Corresponding frequency-space vector. - internal static Complex[] Naive(Complex[] samples, int exponentSign) - { - var w0 = exponentSign*Constants.Pi2/samples.Length; - var spectrum = new Complex[samples.Length]; - - CommonParallel.For(0, samples.Length, (u, v) => - { - for (int i = u; i < v; i++) - { - var wk = w0*i; - var sum = Complex.Zero; - for (var n = 0; n < samples.Length; n++) - { - var w = n*wk; - sum += samples[n]*new Complex(Math.Cos(w), Math.Sin(w)); - } - - spectrum[i] = sum; - } - }); - - return spectrum; - } - } -} diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs index 34396e73..ea45243f 100644 --- a/src/Numerics/IntegralTransforms/Fourier.cs +++ b/src/Numerics/IntegralTransforms/Fourier.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // -// Copyright (c) 2009-2016 Math.NET +// Copyright (c) 2009-2018 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation