diff --git a/src/Managed.UnitTests/IntegralTransformsTests/DhtTest.cs b/src/Managed.UnitTests/IntegralTransformsTests/DhtTest.cs new file mode 100644 index 00000000..a0614695 --- /dev/null +++ b/src/Managed.UnitTests/IntegralTransformsTests/DhtTest.cs @@ -0,0 +1,93 @@ +// +// 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. +// + +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); + } + } +} diff --git a/src/Managed.UnitTests/Managed.UnitTests.csproj b/src/Managed.UnitTests/Managed.UnitTests.csproj index 589a731c..7d796f70 100644 --- a/src/Managed.UnitTests/Managed.UnitTests.csproj +++ b/src/Managed.UnitTests/Managed.UnitTests.csproj @@ -67,6 +67,7 @@ + diff --git a/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs b/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs new file mode 100644 index 00000000..c2e693cd --- /dev/null +++ b/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Naive.cs @@ -0,0 +1,94 @@ +// +// 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. +// + +namespace MathNet.Numerics.IntegralTransforms.Algorithms +{ + using System; + using Threading; + + /// + /// Fast (FHT) Implementation of the Discrete Hartley Transform (DHT). + /// + public partial class DiscreteHartleyTransform + { + /// + /// Naive generic DHT, useful e.g. to verify faster algorithms. + /// + /// Time-space sample vector. + /// Corresponding frequency-space vector. + 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; + } + + /// + /// Naive forward DHT, useful e.g. to verify faster algorithms. + /// + /// Time-space sample vector. + /// Hartley Transform Convention Options. + /// Corresponding frequency-space vector. + public double[] NaiveForward(double[] timeSpace, HartleyOptions options) + { + var frequencySpace = Naive(timeSpace); + ForwardScaleByOptions(options, frequencySpace); + return frequencySpace; + } + + /// + /// Naive inverse DHT, useful e.g. to verify faster algorithms. + /// + /// Frequency-space sample vector. + /// Hartley Transform Convention Options. + /// Corresponding time-space vector. + public double[] NaiveInverse(double[] frequencySpace, HartleyOptions options) + { + var timeSpace = Naive(frequencySpace); + InverseScaleByOptions(options, timeSpace); + return timeSpace; + } + } +} diff --git a/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Options.cs b/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Options.cs new file mode 100644 index 00000000..ddf94166 --- /dev/null +++ b/src/Managed/IntegralTransforms/Algorithms/DiscreteHartleyTransform.Options.cs @@ -0,0 +1,82 @@ +// +// 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. +// + +namespace MathNet.Numerics.IntegralTransforms.Algorithms +{ + using System; + + /// + /// Fast (FHT) Implementation of the Discrete Hartley Transform (DHT). + /// + public partial class DiscreteHartleyTransform + { + /// + /// Rescale FFT-the resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + 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; + } + } + + /// + /// Rescale the iFFT-resulting vector according to the provided convention options. + /// + /// Fourier Transform Convention Options. + /// Sample Vector. + 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; + } + } + } +} diff --git a/src/Managed/IntegralTransforms/HartleyOptions.cs b/src/Managed/IntegralTransforms/HartleyOptions.cs new file mode 100644 index 00000000..fd6f12e0 --- /dev/null +++ b/src/Managed/IntegralTransforms/HartleyOptions.cs @@ -0,0 +1,58 @@ +// +// 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. +// + +namespace MathNet.Numerics.IntegralTransforms +{ + using System; + + /// + /// Hartley Transform Convention + /// + [Flags] + public enum HartleyOptions + { + // FLAGS: + + /// + /// Only scale by 1/N in the inverse direction; No scaling in forward direction. + /// + AsymmetricScaling = 0x02, + + /// + /// Don't scale at all (neither on forward nor on inverse transformation). + /// + NoScaling = 0x04, + + // USABILITY POINTERS: + + /// + /// Universal; Symmetric scaling. + /// + Default = 0, + } +} diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index 069bd8f8..053c9833 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -56,6 +56,9 @@ + + + diff --git a/src/Native.UnitTests/Native.UnitTests.csproj b/src/Native.UnitTests/Native.UnitTests.csproj index 1d0c56ae..9c4f56bf 100644 --- a/src/Native.UnitTests/Native.UnitTests.csproj +++ b/src/Native.UnitTests/Native.UnitTests.csproj @@ -83,6 +83,9 @@ IntegralTransformsTests\DftTest.cs + + IntegralTransformsTests\DhtTest.cs + IntegrationTests\IntegrationTest.cs diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index d8d5ab04..cb10aad0 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -89,9 +89,18 @@ IntegralTransforms\Algorithms\DiscreteFourierTransform.RadixN.cs + + IntegralTransforms\Algorithms\DiscreteHartleyTransform.Naive.cs + + + IntegralTransforms\Algorithms\DiscreteHartleyTransform.Options.cs + IntegralTransforms\FourierOptions.cs + + IntegralTransforms\HartleyOptions.cs + IntegralTransforms\Transform.cs