From 520a32743d44be7f6705b7cb654af4ff8af9d4ec Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Fri, 13 Dec 2013 21:30:59 +0100 Subject: [PATCH] Generate: port signal generation routines from neodym --- src/Numerics/Generate.cs | 361 +++++++++++++++++++++++++++++++++ src/Numerics/Numerics.csproj | 1 + src/UnitTests/GenerateTests.cs | 161 +++++++++++++++ src/UnitTests/UnitTests.csproj | 1 + 4 files changed, 524 insertions(+) create mode 100644 src/Numerics/Generate.cs create mode 100644 src/UnitTests/GenerateTests.cs diff --git a/src/Numerics/Generate.cs b/src/Numerics/Generate.cs new file mode 100644 index 00000000..3665ba51 --- /dev/null +++ b/src/Numerics/Generate.cs @@ -0,0 +1,361 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 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.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.Random; + +namespace MathNet.Numerics +{ + public static class Generate + { + /// + /// Generate a linearly spaced sample vector of the given length between the specified values (inclusive). + /// Equivalent to MATLAB linspace but with the length as first instead of last argument. + /// + public static double[] LinearSpaced(int length, double start, double stop) + { + if (length <= 0) return new double[0]; + if (length == 1) return new[] { stop }; + + double step = (stop - start)/(length - 1); + + var data = new double[length]; + for (int i = 0; i < data.Length; i++) + { + data[i] = start + i*step; + } + data[data.Length - 1] = stop; + return data; + } + + /// + /// Generate a base 10 logarithmically spaced sample vector of the given length between the specified decade exponents (inclusive). + /// Equivalent to MATLAB logspace but with the length as first instead of last argument. + /// + public static double[] LogSpaced(int length, double startExponent, double stopExponent) + { + if (length <= 0) return new double[0]; + if (length == 1) return new[] { Math.Pow(10, stopExponent) }; + + double step = (stopExponent - startExponent)/(length - 1); + + var data = new double[length]; + for (int i = 0; i < data.Length; i++) + { + data[i] = Math.Pow(10, startExponent + i*step); + } + data[data.Length - 1] = Math.Pow(10, stopExponent); + return data; + } + + /// + /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and step 1. + /// Equivalent to MATLAB colon operator (:). + /// + public static double[] LinearRange(int start, int stop) + { + if (start == stop) return new double[] { start }; + if (start < stop) + { + var data = new double[stop - start + 1]; + for (int i = 0; i < data.Length; i++) + { + data[i] = start + i; + } + return data; + } + else + { + var data = new double[start - stop + 1]; + for (int i = 0; i < data.Length; i++) + { + data[i] = start - i; + } + return data; + } + } + + /// + /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and the provide step. + /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step. + /// Equivalent to MATLAB double colon operator (::). + /// + public static double[] LinearRange(int start, int step, int stop) + { + if (start == stop) return new double[] { start }; + if (start < stop && step < 0 || start > stop && step > 0 || step == 0d) + { + return new double[0]; + } + + var data = new double[(stop - start)/step + 1]; + for (int i = 0; i < data.Length; i++) + { + data[i] = start + i*step; + } + return data; + } + + /// + /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and the provide step. + /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step. + /// Equivalent to MATLAB double colon operator (::). + /// + public static double[] LinearRange(double start, double step, double stop) + { + if (start == stop) return new double[] { start }; + if (start < stop && step < 0 || start > stop && step > 0 || step == 0d) + { + return new double[0]; + } + + var data = new double[(int)Math.Floor((stop - start)/step + 1d)]; + for (int i = 0; i < data.Length; i++) + { + data[i] = start + i*step; + } + return data; + } + + /// + /// Create a Sine sample vector. + /// + /// The number of samples to generate. + /// Samples per unit. + /// Frequency in samples per unit. + /// The maximal reached peak. + /// The mean, or dc part, of the signal. + /// Optional phase offset. + /// Optional delay, relative to the phase. + public static double[] Sinusoidal(int length, double samplingRate, double frequency, double amplitude, double mean = 0.0, double phase = 0.0, int delay = 0) + { + double step = frequency/samplingRate*Constants.Pi2; + phase = (phase - delay*step)%Constants.Pi2; + + var data = new double[length]; + for (int i = 0; i < length; i++) + { + data[i] = mean + amplitude*Math.Sin(phase + i*step); + } + return data; + } + + /// + /// Create an infinite Sine sample sequence. + /// + /// Samples per unit. + /// Frequency in samples per unit. + /// The maximal reached peak. + /// The mean, or dc part, of the signal. + /// Optional phase offset. + /// Optional delay, relative to the phase. + public static IEnumerable SinusoidalSequence(double samplingRate, double frequency, double amplitude, double mean = 0.0, double phase = 0.0, int delay = 0) + { + double step = frequency/samplingRate*Constants.Pi2; + phase = (phase - delay*step)%Constants.Pi2; + + while (true) + { + for (int i = 0; i < 1000; i++) + { + yield return mean + amplitude*Math.Sin(phase + i*step); + } + phase = (phase + 1000*step)%Constants.Pi2; + } + } + + /// + /// Create a Heaviside Step sample vector. + /// + /// The number of samples to generate. + /// The maximal reached peak. + /// Offset to the time axis. + public static double[] Step(int length, double amplitude, int delay) + { + var data = new double[length]; + for (int i = Math.Max(0, delay); i < data.Length; i++) + { + data[i] = amplitude; + } + return data; + } + + /// + /// Create an infinite Heaviside Step sample sequence. + /// + /// The maximal reached peak. + /// Offset to the time axis. + public static IEnumerable StepSequence(double amplitude, int delay) + { + for (int i = 0; i < delay; i++) + { + yield return 0d; + } + + while (true) + { + yield return amplitude; + } + } + + /// + /// Create a Dirac Delta Impulse sample vector. + /// + /// The number of samples to generate. + /// impulse sequence period. -1 for single impulse only. + /// The maximal reached peak. + /// Offset to the time axis. Zero or positive. + public static double[] Impulse(int length, int period, double amplitude, int delay) + { + var data = new double[length]; + if (period <= 0) + { + if (delay >= 0 && delay < length) + { + data[delay] = amplitude; + } + } + else + { + delay = ((delay%period) + period)%period; + while (delay < length) + { + data[delay] = amplitude; + delay += period; + } + } + return data; + } + + + /// + /// Create a Dirac Delta Impulse sample vector. + /// + /// impulse sequence period. -1 for single impulse only. + /// The maximal reached peak. + /// Offset to the time axis. Zero or positive. + public static IEnumerable ImpulseSequence(int period, double amplitude, int delay) + { + if (period <= 0) + { + for (int i = 0; i < delay; i++) + { + yield return 0d; + } + + yield return amplitude; + + while (true) + { + yield return 0d; + } + } + else + { + delay = ((delay%period) + period)%period; + + for (int i = 0; i < delay; i++) + { + yield return 0d; + } + + while (true) + { + yield return amplitude; + + for (int i = 1; i < period; i++) + { + yield return 0d; + } + } + } + } + + /// + /// Create random samples. + /// + public static double[] Random(int length, IContinuousDistribution distribution) + { + return distribution.Samples().Take(length).ToArray(); + } + + /// + /// Create an infinite random sample sequence. + /// + public static IEnumerable Random(IContinuousDistribution distribution) + { + return distribution.Samples(); + } + + /// + /// Create samples with independent amplitudes of normal distribution and a flat spectral density. + /// + public static double[] WhiteGaussianNoise(int length, double mean, double standardDeviation) + { + return Normal.Samples(MersenneTwister.Default, mean, standardDeviation).Take(length).ToArray(); + } + + /// + /// Create an infinite sample sequence with independent amplitudes of normal distribution and a flat spectral density. + /// + public static IEnumerable WhiteGaussianNoiseSequence(double mean, double standardDeviation) + { + return Normal.Samples(MersenneTwister.Default, mean, standardDeviation); + } + + /// + /// Create skew alpha stable samples. + /// + /// The number of samples to generate. + /// Stability alpha-parameter of the stable distribution + /// Skewness beta-parameter of the stable distribution + /// Scale c-parameter of the stable distribution + /// Location mu-parameter of the stable distribution + public static double[] StableNoise(int length, double alpha, double beta, double scale, double location) + { + return Stable.Samples(MersenneTwister.Default, alpha, beta, scale, location).Take(length).ToArray(); + } + + /// + /// Create skew alpha stable samples. + /// + /// Stability alpha-parameter of the stable distribution + /// Skewness beta-parameter of the stable distribution + /// Scale c-parameter of the stable distribution + /// Location mu-parameter of the stable distribution + public static IEnumerable StableNoiseSequence(double alpha, double beta, double scale, double location) + { + return Stable.Samples(MersenneTwister.Default, alpha, beta, scale, location); + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 684c5d2b..97f7ef12 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -87,6 +87,7 @@ + diff --git a/src/UnitTests/GenerateTests.cs b/src/UnitTests/GenerateTests.cs new file mode 100644 index 00000000..80e97483 --- /dev/null +++ b/src/UnitTests/GenerateTests.cs @@ -0,0 +1,161 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 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.Linq; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests +{ + [TestFixture] + public class GenerateTests + { + [Test] + public void LinearSpaced() + { + Assert.That(Generate.LinearSpaced(0, 0d, 2d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearSpaced(1, 0d, 2d), Is.EqualTo(new[] { 2d }).AsCollection); + Assert.That(Generate.LinearSpaced(2, 0d, 2d), Is.EqualTo(new[] { 0d, 2d }).AsCollection); + Assert.That(Generate.LinearSpaced(3, 0d, 2d), Is.EqualTo(new[] { 0d, 1d, 2d }).AsCollection); + Assert.That(Generate.LinearSpaced(4, 0d, 2d), Is.EqualTo(new[] { 0d, 2d/3d, 4d/3d, 2d }).Within(1e-12).AsCollection); + + Assert.That(Generate.LinearSpaced(0, 2d, 0d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearSpaced(1, 2d, 0d), Is.EqualTo(new[] { 0d }).AsCollection); + Assert.That(Generate.LinearSpaced(2, 2d, 0d), Is.EqualTo(new[] { 2d, 0d }).AsCollection); + Assert.That(Generate.LinearSpaced(3, 2d, 0d), Is.EqualTo(new[] { 2d, 1d, 0d }).AsCollection); + Assert.That(Generate.LinearSpaced(4, 2d, 0d), Is.EqualTo(new[] { 2d, 4d/3d, 2d/3d, 0d }).Within(1e-12).AsCollection); + } + + [Test] + public void LogSpaced() + { + Assert.That(Generate.LogSpaced(0, 0d, 2d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LogSpaced(1, 0d, 2d), Is.EqualTo(new[] { 100.0 }).AsCollection); + Assert.That(Generate.LogSpaced(2, 0d, 2d), Is.EqualTo(new[] { 1.0, 100.0 }).AsCollection); + Assert.That(Generate.LogSpaced(3, 0d, 2d), Is.EqualTo(new[] { 1.0, 10.0, 100.0 }).AsCollection); + Assert.That(Generate.LogSpaced(4, 0d, 2d), Is.EqualTo(new[] { 1.0, Math.Pow(10.0, 2.0/3.0), Math.Pow(10.0, 4.0/3.0), 100.0 }).Within(1e-12).AsCollection); + + Assert.That(Generate.LogSpaced(0, 2d, 0d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LogSpaced(1, 2d, 0d), Is.EqualTo(new[] { 1.0 }).AsCollection); + Assert.That(Generate.LogSpaced(2, 2d, 0d), Is.EqualTo(new[] { 100.0, 1.0 }).AsCollection); + Assert.That(Generate.LogSpaced(3, 2d, 0d), Is.EqualTo(new[] { 100.0, 10.0, 1.0 }).AsCollection); + Assert.That(Generate.LogSpaced(4, 2d, 0d), Is.EqualTo(new[] { 100.0, Math.Pow(10.0, 4.0/3.0), Math.Pow(10, 2.0/3.0), 1.0 }).Within(1e-12).AsCollection); + + Assert.That(Generate.LogSpaced(5, -2d, 2d), Is.EqualTo(new[] { 0.01, 0.1, 1.0, 10.0, 100.0 }).AsCollection); + Assert.That(Generate.LogSpaced(5, 2d, -2d), Is.EqualTo(new[] { 100.0, 10.0, 1.0, 0.1, 0.01 }).AsCollection); + } + + [Test] + public void LinearRange() + { + Assert.That(Generate.LinearRange(1, 1), Is.EqualTo(new[] { 1d }).AsCollection); + + Assert.That(Generate.LinearRange(1, 3), Is.EqualTo(new[] { 1d, 2d, 3d }).AsCollection); + Assert.That(Generate.LinearRange(-1, -3), Is.EqualTo(new[] { -1d, -2d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(-3, -1), Is.EqualTo(new[] { -3d, -2d, -1d }).AsCollection); + Assert.That(Generate.LinearRange(1, -2), Is.EqualTo(new[] { 1d, 0d, -1d, -2d }).AsCollection); + } + + [Test] + public void LinearRangeStep() + { + Assert.That(Generate.LinearRange(1, 1, 1), Is.EqualTo(new[] { 1d }).AsCollection); + + Assert.That(Generate.LinearRange(1, -1, 2), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(2, 1, 1), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(1, 0, 2), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(2, 0, 1), Is.EqualTo(new double[0]).AsCollection); + + Assert.That(Generate.LinearRange(1, 1, 5), Is.EqualTo(new[] { 1d, 2d, 3d, 4d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1, 2, 5), Is.EqualTo(new[] { 1d, 3d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1, 2, 6), Is.EqualTo(new[] { 1d, 3d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1, 2, 4), Is.EqualTo(new[] { 1d, 3d }).AsCollection); + + Assert.That(Generate.LinearRange(1, -1, -3), Is.EqualTo(new[] { 1d, 0d, -1d, -2d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1, -2, -3), Is.EqualTo(new[] { 1d, -1d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1, -2, -4), Is.EqualTo(new[] { 1d, -1d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1, -2, -2), Is.EqualTo(new[] { 1d, -1d }).AsCollection); + } + + [Test] + public void LinearRangeFloatingPoint() + { + Assert.That(Generate.LinearRange(1d, 1d, 1d), Is.EqualTo(new[] { 1d }).AsCollection); + + Assert.That(Generate.LinearRange(1d, -1d, 2d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(2d, 1d, 1d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(1d, 0d, 2d), Is.EqualTo(new double[0]).AsCollection); + Assert.That(Generate.LinearRange(2d, 0d, 1d), Is.EqualTo(new double[0]).AsCollection); + + Assert.That(Generate.LinearRange(1d, 1d, 5d), Is.EqualTo(new[] { 1d, 2d, 3d, 4d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1d, 2d, 5d), Is.EqualTo(new[] { 1d, 3d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1d, 2d, 6d), Is.EqualTo(new[] { 1d, 3d, 5d }).AsCollection); + Assert.That(Generate.LinearRange(1d, 2d, 4d), Is.EqualTo(new[] { 1d, 3d }).AsCollection); + Assert.That(Generate.LinearRange(1d, 1.5d, 5d), Is.EqualTo(new[] { 1d, 2.5d, 4d }).AsCollection); + Assert.That(Generate.LinearRange(1d, 1.5d, 6.5d), Is.EqualTo(new[] { 1d, 2.5d, 4d, 5.5d }).AsCollection); + + Assert.That(Generate.LinearRange(1d, -1d, -3d), Is.EqualTo(new[] { 1d, 0d, -1d, -2d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1d, -2d, -3d), Is.EqualTo(new[] { 1d, -1d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1d, -2d, -4d), Is.EqualTo(new[] { 1d, -1d, -3d }).AsCollection); + Assert.That(Generate.LinearRange(1d, -2d, -2d), Is.EqualTo(new[] { 1d, -1d }).AsCollection); + Assert.That(Generate.LinearRange(1d, -1.5d, -3d), Is.EqualTo(new[] { 1d, -0.5d, -2d }).AsCollection); + Assert.That(Generate.LinearRange(1d, -1.5d, -3.5d), Is.EqualTo(new[] { 1d, -0.5d, -2d, -3.5 }).AsCollection); + Assert.That(Generate.LinearRange(1d, -1.5d, -4d), Is.EqualTo(new[] { 1d, -0.5d, -2d, -3.5 }).AsCollection); + } + + [Test] + public void SinusoidalConsistentWithSequence() + { + Assert.That( + Generate.SinusoidalSequence(32, 2, 5, 1, 0.5, -6).Take(1000).ToArray(), + Is.EqualTo(Generate.Sinusoidal(1000, 32, 2, 5, 1, 0.5, -6)).AsCollection); + } + + [Test] + public void StepConsistentWithSequence() + { + Assert.That( + Generate.StepSequence(5, 40).Take(1000).ToArray(), + Is.EqualTo(Generate.Step(1000, 5, 40)).AsCollection); + } + + [Test] + public void ImpulseConsistentWithSequence() + { + Assert.That( + Generate.ImpulseSequence(0, 5, 40).Take(1000).ToArray(), + Is.EqualTo(Generate.Impulse(1000, 0, 5, 40)).AsCollection); + + Assert.That( + Generate.ImpulseSequence(100, 5, 40).Take(1000).ToArray(), + Is.EqualTo(Generate.Impulse(1000, 100, 5, 40)).AsCollection); + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 56026b9c..491bf401 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -133,6 +133,7 @@ +