Browse Source

FFT: cleanup tests

spatial
Christoph Ruegg 9 years ago
parent
commit
521fe92aba
  1. 74
      src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs
  2. 58
      src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs
  3. 178
      src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs
  4. 73
      src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs
  5. 101
      src/Numerics/IntegralTransforms/Fourier.Naive.cs
  6. 2
      src/Numerics/IntegralTransforms/Fourier.cs

74
src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs

@ -40,11 +40,8 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
[TestFixture, Category("FFT")]
public class FourierTest
{
/// <summary>
/// Naive transforms real sine correctly.
/// </summary>
[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
}
}
/// <summary>
/// Naive transforms real sine correctly.
/// </summary>
[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");
}
}
}
}
}

58
src/Numerics.Tests/IntegralTransformsTests/InverseTransformTest.cs

@ -54,7 +54,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
/// <param name="options">Fourier options.</param>
[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
/// <param name="options">Fourier options.</param>
[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
/// <param name="options">Fourier options.</param>
[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
/// <param name="options">Fourier options.</param>
[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
/// <param name="options">Fourier options.</param>
[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
/// <param name="options">Hartley options.</param>
[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);
}
/// <summary>
/// Fourier default transform is reversible.
/// </summary>
[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);
}
/// <summary>
/// Fourier default transform is reversible.
/// </summary>
[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);
}
}
}

178
src/Numerics.Tests/IntegralTransformsTests/MatchingNaiveTransformTest.cs → src/Numerics.Tests/IntegralTransformsTests/MatchingReferenceTransformTest.cs

@ -40,7 +40,7 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
/// Matching Naive transform tests.
/// </summary>
[TestFixture, Category("FFT")]
public class MatchingNaiveTransformTest
public class MatchingReferenceTransformTest
{
/// <summary>
/// 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<Complex32[], FourierTransformScaling> expected,
Action<Complex32[], FourierTransformScaling> 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<Complex[], FourierTransformScaling> expected,
Action<Complex[], FourierTransformScaling> 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);
}
/// <summary>
/// Fourier Radix2XX matches naive on real sine.
/// </summary>
@ -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);
}
/// <summary>
/// Fourier bluestein matches providers on random power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[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);
}
/// <summary>
/// Fourier bluestein matches providers on random power of two.
/// </summary>
/// <param name="options">Fourier options.</param>
[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);
}
}
}

73
src/Numerics.Tests/IntegralTransformsTests/ParsevalTheoremTest.cs

@ -56,19 +56,15 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
/// <param name="count">Samples count.</param>
[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
/// <param name="count">Samples count.</param>
[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);
}
/// <summary>
/// Fourier default transform satisfies Parseval's theorem.
/// </summary>
/// <param name="count">Samples count.</param>
[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);
/// <summary>
/// Fourier default transform satisfies Parseval's theorem.
/// </summary>
/// <param name="count">Samples count.</param>
[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
/// <param name="count">Samples count.</param>
[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);
}
}

101
src/Numerics/IntegralTransforms/Fourier.Naive.cs

@ -1,101 +0,0 @@
// <copyright file="Fourier.Naive.cs" company="Math.NET">
// 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.
// </copyright>
using System;
using System.Numerics;
using MathNet.Numerics.Threading;
namespace MathNet.Numerics.IntegralTransforms
{
/// <summary>
/// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
/// </summary>
public static partial class Fourier
{
/// <summary>
/// Naive generic DFT, useful e.g. to verify faster algorithms.
/// </summary>
/// <param name="samples">Time-space sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <returns>Corresponding frequency-space vector.</returns>
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;
}
/// <summary>
/// Naive generic DFT, useful e.g. to verify faster algorithms.
/// </summary>
/// <param name="samples">Time-space sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <returns>Corresponding frequency-space vector.</returns>
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;
}
}
}

2
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

Loading…
Cancel
Save