transforms: hartley - naive

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>

src/Managed.UnitTests/IntegralTransformsTests/DhtTest.cs

@@ -0,0 +1,93 @@
+// <copyright file="DhtTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 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.Collections.Generic;
+
+namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
+{
+    using System;
+    using System.Linq;
+    using MbUnit.Framework;
+    using IntegralTransforms;
+    using IntegralTransforms.Algorithms;
+    using Statistics;
+
+    [TestFixture]
+    public class DhtTest
+    {
+        private static Random _random = new Random();
+
+        private static double[] ProvideSamples(int count)
+        {
+            var samples = new double[count];
+            for (int i = 0; i < samples.Length; i++)
+            {
+                samples[i] = 1 - (2 * _random.NextDouble());
+            }
+
+            return samples;
+        }
+
+        [Test]
+        [Row(HartleyOptions.Default, FourierOptions.Default)]
+        [Row(HartleyOptions.AsymmetricScaling, FourierOptions.AsymmetricScaling)]
+        [Row(HartleyOptions.NoScaling, FourierOptions.NoScaling)]
+        public void NaiveMatchesDFT(HartleyOptions hartleyOptions, FourierOptions fourierOptions)
+        {
+            var samples = ProvideSamples(0x80);
+
+            var dht = new DiscreteHartleyTransform();
+            var hartley = dht.NaiveForward(samples, hartleyOptions);
+
+            var fourierComplex = Array.ConvertAll(samples, s => new Complex(s, s));
+            Transform.FourierForward(fourierComplex, fourierOptions);
+            var fourierReal = Array.ConvertAll(fourierComplex, s => s.Real);
+
+            AssertHelpers.AlmostEqualList(fourierReal, hartley, 1e-10);
+        }
+
+        [Test]
+        [Row(HartleyOptions.Default)]
+        [Row(HartleyOptions.AsymmetricScaling)]
+        public void NaiveIsReversible(HartleyOptions options)
+        {
+            var samples = ProvideSamples(0x80);
+            var work = new double[samples.Length];
+            samples.CopyTo(work, 0);
+
+            var dht = new DiscreteHartleyTransform();
+            work = dht.NaiveForward(work, options);
+
+            Assert.IsFalse(work.AlmostEqualListWithError(samples, 1e-10));
+
+            work = dht.NaiveInverse(work, options);
+
+            AssertHelpers.AlmostEqualList(samples, work, 1e-10);
+        }
+    }
+}

src/Managed.UnitTests/Managed.UnitTests.csproj

@@ -67,6 +67,7 @@
     <Compile Include="DistributionTests\Continuous\ContinuousUniformTests.cs" />
     <Compile Include="DistributionTests\Continuous\GammaTests.cs" />
     <Compile Include="DistributionTests\Continuous\NormalTests.cs" />
+    <Compile Include="IntegralTransformsTests\DhtTest.cs" />
     <Compile Include="IntegralTransformsTests\DftTest.cs" />
     <Compile Include="IntegrationTests\IntegrationTest.cs" />
     <Compile Include="InterpolationTests\InterpolationContract.cs" />

src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs

@@ -0,0 +1,94 @@
+// <copyright file="DiscreteHartleyTransform.Naive.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 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>
+
+namespace MathNet.Numerics.IntegralTransforms.Algorithms
+{
+    using System;
+    using Threading;
+
+    /// <summary>
+    /// Fast (FHT) Implementation of the Discrete Hartley Transform (DHT).
+    /// </summary>
+    public partial class DiscreteHartleyTransform
+    {
+        /// <summary>
+        /// Naive generic DHT, useful e.g. to verify faster algorithms.
+        /// </summary>
+        /// <param name="samples">Time-space sample vector.</param>
+        /// <returns>Corresponding frequency-space vector.</returns>
+        internal static double[] Naive(double[] samples)
+        {
+            double w0 = 2 * Constants.Pi / samples.Length;
+            var spectrum = new double[samples.Length];
+
+            Parallel.For(
+                0,
+                samples.Length,
+                k =>
+                {
+                    double wk = w0 * k;
+                    double sum = 0.0;
+                    for (int n = 0; n < samples.Length; n++)
+                    {
+                        double w = n * wk;
+                        sum += samples[n] * Constants.Sqrt2 * Math.Cos(w - Constants.PiOver4);
+                    }
+
+                    spectrum[k] = sum;
+                });
+
+            return spectrum;
+        }
+
+        /// <summary>
+        /// Naive forward DHT, useful e.g. to verify faster algorithms.
+        /// </summary>
+        /// <param name="timeSpace">Time-space sample vector.</param>
+        /// <param name="options">Hartley Transform Convention Options.</param>
+        /// <returns>Corresponding frequency-space vector.</returns>
+        public double[] NaiveForward(double[] timeSpace, HartleyOptions options)
+        {
+            var frequencySpace = Naive(timeSpace);
+            ForwardScaleByOptions(options, frequencySpace);
+            return frequencySpace;
+        }
+
+        /// <summary>
+        /// Naive inverse DHT, useful e.g. to verify faster algorithms.
+        /// </summary>
+        /// <param name="frequencySpace">Frequency-space sample vector.</param>
+        /// <param name="options">Hartley Transform Convention Options.</param>
+        /// <returns>Corresponding time-space vector.</returns>
+        public double[] NaiveInverse(double[] frequencySpace, HartleyOptions options)
+        {
+            var timeSpace = Naive(frequencySpace);
+            InverseScaleByOptions(options, timeSpace);
+            return timeSpace;
+        }
+    }
+}

src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Options.cs

@@ -0,0 +1,82 @@
+// <copyright file="DiscreteHartleyTransform.Options.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 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>
+
+namespace MathNet.Numerics.IntegralTransforms.Algorithms
+{
+    using System;
+
+    /// <summary>
+    /// Fast (FHT) Implementation of the Discrete Hartley Transform (DHT).
+    /// </summary>
+    public partial class DiscreteHartleyTransform
+    {
+        /// <summary>
+        /// Rescale FFT-the resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        private static void ForwardScaleByOptions(HartleyOptions options, double[] samples)
+        {
+            if ((options & HartleyOptions.NoScaling) == HartleyOptions.NoScaling ||
+                (options & HartleyOptions.AsymmetricScaling) == HartleyOptions.AsymmetricScaling)
+            {
+                return;
+            }
+
+            var scalingFactor = Math.Sqrt(1.0 / samples.Length);
+            for (int i = 0; i < samples.Length; i++)
+            {
+                samples[i] *= scalingFactor;
+            }
+        }
+
+        /// <summary>
+        /// Rescale the iFFT-resulting vector according to the provided convention options.
+        /// </summary>
+        /// <param name="options">Fourier Transform Convention Options.</param>
+        /// <param name="samples">Sample Vector.</param>
+        private static void InverseScaleByOptions(HartleyOptions options, double[] samples)
+        {
+            if ((options & HartleyOptions.NoScaling) == HartleyOptions.NoScaling)
+            {
+                return;
+            }
+
+            var scalingFactor = 1.0 / samples.Length;
+            if ((options & HartleyOptions.AsymmetricScaling) != HartleyOptions.AsymmetricScaling)
+            {
+                scalingFactor = Math.Sqrt(scalingFactor);
+            }
+
+            for (int i = 0; i < samples.Length; i++)
+            {
+                samples[i] *= scalingFactor;
+            }
+        }
+    }
+}

src/Managed/IntegralTransforms/HartleyOptions.cs

@@ -0,0 +1,58 @@
+// <copyright file="HartleyOptions.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://mathnet.opensourcedotnet.info
+//
+// Copyright (c) 2009 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>
+
+namespace MathNet.Numerics.IntegralTransforms
+{
+    using System;
+
+    /// <summary>
+    /// Hartley Transform Convention
+    /// </summary>
+    [Flags]
+    public enum HartleyOptions
+    {
+        // FLAGS:
+
+        /// <summary>
+        /// Only scale by 1/N in the inverse direction; No scaling in forward direction.
+        /// </summary>
+        AsymmetricScaling = 0x02,
+
+        /// <summary>
+        /// Don't scale at all (neither on forward nor on inverse transformation).
+        /// </summary>
+        NoScaling = 0x04,
+
+        // USABILITY POINTERS:
+
+        /// <summary>
+        /// Universal; Symmetric scaling.
+        /// </summary>
+        Default = 0,
+    }
+}

src/Managed/Managed.csproj

@@ -56,6 +56,9 @@
     <Compile Include="Distributions\IContinuousDistribution.cs" />
     <Compile Include="Distributions\IDiscreteDistribution.cs" />
     <Compile Include="Distributions\IDistribution.cs" />
+    <Compile Include="IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs" />
+    <Compile Include="IntegralTransforms\Algorithms\DiscreteHartleyTransform.Options.cs" />
+    <Compile Include="IntegralTransforms\HartleyOptions.cs" />
     <Compile Include="IPrecisionSupport.cs" />
     <Compile Include="IntegralTransforms\Algorithms\DiscreteFourierTransform.Options.cs" />
     <Compile Include="IntegralTransforms\Algorithms\DiscreteFourierTransform.Bluestein.cs" />

src/Native.UnitTests/Native.UnitTests.csproj

@@ -83,6 +83,9 @@
     <Compile Include="..\Managed.UnitTests\IntegralTransformsTests\DftTest.cs">
       <Link>IntegralTransformsTests\DftTest.cs</Link>
     </Compile>
+    <Compile Include="..\Managed.UnitTests\IntegralTransformsTests\DhtTest.cs">
+      <Link>IntegralTransformsTests\DhtTest.cs</Link>
+    </Compile>
     <Compile Include="..\Managed.UnitTests\IntegrationTests\IntegrationTest.cs">
       <Link>IntegrationTests\IntegrationTest.cs</Link>
     </Compile>

src/Native/Native.csproj

@@ -89,9 +89,18 @@
     <Compile Include="..\Managed\IntegralTransforms\Algorithms\DiscreteFourierTransform.RadixN.cs">
       <Link>IntegralTransforms\Algorithms\DiscreteFourierTransform.RadixN.cs</Link>
     </Compile>
+    <Compile Include="..\Managed\IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs">
+      <Link>IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs</Link>
+    </Compile>
+    <Compile Include="..\Managed\IntegralTransforms\Algorithms\DiscreteHartleyTransform.Options.cs">
+      <Link>IntegralTransforms\Algorithms\DiscreteHartleyTransform.Options.cs</Link>
+    </Compile>
     <Compile Include="..\Managed\IntegralTransforms\FourierOptions.cs">
       <Link>IntegralTransforms\FourierOptions.cs</Link>
     </Compile>
+    <Compile Include="..\Managed\IntegralTransforms\HartleyOptions.cs">
+      <Link>IntegralTransforms\HartleyOptions.cs</Link>
+    </Compile>
     <Compile Include="..\Managed\IntegralTransforms\Transform.cs">
       <Link>IntegralTransforms\Transform.cs</Link>
     </Compile>