Christoph Ruegg
17 years ago
10 changed files with 843 additions and 0 deletions
@ -0,0 +1,246 @@ |
|||
// <copyright file="DftTest.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.UnitTests.IntegralTransformsTests |
|||
{ |
|||
using System; |
|||
using MbUnit.Framework; |
|||
using IntegralTransforms; |
|||
using IntegralTransforms.Algorithms; |
|||
|
|||
[TestFixture] |
|||
public class DftTest |
|||
{ |
|||
private static Random _random = new Random(); |
|||
|
|||
private static Complex[] ProvideSamples(int count) |
|||
{ |
|||
var samples = new Complex[count]; |
|||
for (int i = 0; i < samples.Length; i++) |
|||
{ |
|||
samples[i] = Complex.WithRealImaginary( |
|||
1 - (2 * _random.NextDouble()), |
|||
1 - (2 * _random.NextDouble())); |
|||
} |
|||
|
|||
return samples; |
|||
} |
|||
|
|||
[Test] |
|||
public void NaiveTransformsRealSineCorrectly() |
|||
{ |
|||
var realSine = new Complex[16]; |
|||
for (int i = 0; i < realSine.Length; i++) |
|||
{ |
|||
realSine[i] = Math.Sin(i / 8.0 * Constants.Pi); |
|||
} |
|||
|
|||
// real-odd transforms to imaginary odd
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrum = dft.NaiveForward(realSine, FourierOptions.Matlab); |
|||
|
|||
// all real components must be zero
|
|||
foreach (var c in spectrum) |
|||
{ |
|||
Assert.AreApproximatelyEqual(0, c.Real, 1e-12, "real"); |
|||
} |
|||
|
|||
// all imaginary components except second and last musth be zero
|
|||
for(int i = 0; i<spectrum.Length; i++) |
|||
{ |
|||
if(i == 1) |
|||
{ |
|||
Assert.AreApproximatelyEqual(-8, spectrum[i].Imaginary, 1e-12, "imag second"); |
|||
} |
|||
else if (i == spectrum.Length - 1) |
|||
{ |
|||
Assert.AreApproximatelyEqual(8, spectrum[i].Imaginary, 1e-12, "imag last"); |
|||
} |
|||
else |
|||
{ |
|||
Assert.AreApproximatelyEqual(0, spectrum[i].Imaginary, 1e-12, "imag"); |
|||
} |
|||
} |
|||
} |
|||
|
|||
[Test] |
|||
[Row(FourierOptions.Default)] |
|||
[Row(FourierOptions.Matlab)] |
|||
public void NaiveIsReversible(FourierOptions options) |
|||
{ |
|||
var samples = ProvideSamples(0x80); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
work = dft.NaiveForward(work, options); |
|||
|
|||
Assert.IsFalse(work.AlmostEqualListWithError(samples, 1e-12)); |
|||
|
|||
work = dft.NaiveInverse(work, options); |
|||
|
|||
AssertHelpers.AlmostEqualList(samples, work, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
public void Radix2MatchesNaiveOnRealSine() |
|||
{ |
|||
var realSine = new Complex[16]; |
|||
for (int i = 0; i < realSine.Length; i++) |
|||
{ |
|||
realSine[i] = Math.Sin(i / 8.0 * Constants.Pi); |
|||
} |
|||
|
|||
// real-odd transforms to imaginary odd
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrumNaive = dft.NaiveForward(realSine, FourierOptions.Matlab); |
|||
|
|||
var spectrumRadix2 = new Complex[realSine.Length]; |
|||
realSine.CopyTo(spectrumRadix2, 0); |
|||
dft.Radix2Forward(spectrumRadix2, FourierOptions.Matlab); |
|||
|
|||
AssertHelpers.AlmostEqualList(spectrumNaive, spectrumRadix2, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
public void Radix2MatchesNaiveOnRandom() |
|||
{ |
|||
var samples = ProvideSamples(0x80); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrumNaive = dft.NaiveForward(samples, FourierOptions.Matlab); |
|||
dft.Radix2Forward(work, FourierOptions.Matlab); |
|||
|
|||
AssertHelpers.AlmostEqualList(spectrumNaive, work, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
[Row(FourierOptions.Default)] |
|||
[Row(FourierOptions.Matlab)] |
|||
public void Radix2IsReversible(FourierOptions options) |
|||
{ |
|||
var samples = ProvideSamples(0x8000); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
dft.Radix2Forward(work, options); |
|||
|
|||
Assert.IsFalse(work.AlmostEqualListWithError(samples, 1e-12)); |
|||
|
|||
dft.Radix2Inverse(work, options); |
|||
|
|||
AssertHelpers.AlmostEqualList(samples, work, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
public void Radix2ThrowsWhenNotPowerOfTwo() |
|||
{ |
|||
var samples = ProvideSamples(0x7F); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
|
|||
Assert.Throws( |
|||
typeof(ArgumentException), |
|||
() => dft.Radix2Forward(samples, FourierOptions.Default)); |
|||
|
|||
Assert.Throws( |
|||
typeof(ArgumentException), |
|||
() => dft.Radix2Inverse(samples, FourierOptions.Default)); |
|||
} |
|||
|
|||
[Test] |
|||
public void BluesteinMatchesNaiveOnRealSine() |
|||
{ |
|||
var realSine = new Complex[14]; |
|||
for (int i = 0; i < realSine.Length; i++) |
|||
{ |
|||
realSine[i] = Math.Sin(i / 7.0 * Constants.Pi); |
|||
} |
|||
|
|||
// real-odd transforms to imaginary odd
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrumNaive = dft.NaiveForward(realSine, FourierOptions.Matlab); |
|||
|
|||
var spectrumBluestein = new Complex[realSine.Length]; |
|||
realSine.CopyTo(spectrumBluestein, 0); |
|||
dft.BluesteinForward(spectrumBluestein, FourierOptions.Matlab); |
|||
|
|||
AssertHelpers.AlmostEqualList(spectrumNaive, spectrumBluestein, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
public void BluesteinMatchesNaiveOnRandomPowerOfTwo() |
|||
{ |
|||
var samples = ProvideSamples(0x80); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrumNaive = dft.NaiveForward(samples, FourierOptions.Matlab); |
|||
dft.BluesteinForward(work, FourierOptions.Matlab); |
|||
|
|||
AssertHelpers.AlmostEqualList(spectrumNaive, work, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
public void BluesteinMatchesNaiveOnRandomNonPowerOfTwo() |
|||
{ |
|||
var samples = ProvideSamples(0x7F); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
var spectrumNaive = dft.NaiveForward(samples, FourierOptions.Matlab); |
|||
dft.BluesteinForward(work, FourierOptions.Matlab); |
|||
|
|||
AssertHelpers.AlmostEqualList(spectrumNaive, work, 1e-12); |
|||
} |
|||
|
|||
[Test] |
|||
[Row(FourierOptions.Default)] |
|||
[Row(FourierOptions.Matlab)] |
|||
public void BluesteinIsReversible(FourierOptions options) |
|||
{ |
|||
var samples = ProvideSamples(0x7FFF); |
|||
var work = new Complex[samples.Length]; |
|||
samples.CopyTo(work, 0); |
|||
|
|||
var dft = new DiscreteFourierTransform(); |
|||
dft.BluesteinForward(work, options); |
|||
|
|||
Assert.IsFalse(work.AlmostEqualListWithError(samples, 1e-12)); |
|||
|
|||
dft.BluesteinInverse(work, options); |
|||
|
|||
AssertHelpers.AlmostEqualList(samples, work, 1e-12); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,175 @@ |
|||
// <copyright file="DiscreteFourierTransform.Bluestein.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 NumberTheory; |
|||
using Threading; |
|||
|
|||
/// <summary>
|
|||
/// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
|
|||
/// </summary>
|
|||
public partial class DiscreteFourierTransform |
|||
{ |
|||
/// <summary>
|
|||
/// Generate the bluestein sequence for the provided problem size.
|
|||
/// </summary>
|
|||
/// <param name="n">Number of samples.</param>
|
|||
/// <returns>Bluestein sequence exp(I*Pi*k^2/N)</returns>
|
|||
private static Complex[] BluesteinSequence(int n) |
|||
{ |
|||
double s = Constants.Pi / n; |
|||
var sequence = new Complex[n]; |
|||
|
|||
for (int k = 0; k < sequence.Length; k++) |
|||
{ |
|||
double t = s * (k * k); |
|||
sequence[k] = Complex.WithRealImaginary(Math.Cos(t), Math.Sin(t)); |
|||
} |
|||
|
|||
return sequence; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Convolution with the bluestein sequence.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample Vector.</param>
|
|||
private static void BluesteinConvolution(Complex[] samples) |
|||
{ |
|||
int n = samples.Length; |
|||
Complex[] sequence = BluesteinSequence(n); |
|||
|
|||
// Padding to power of two >= 2N–1 so we can apply Radix-2 FFT.
|
|||
int m = ((n << 1) - 1).CeilingToPowerOfTwo(); |
|||
Complex[] b = new Complex[m]; |
|||
Complex[] a = new Complex[m]; |
|||
|
|||
Parallel.Invoke( |
|||
() => |
|||
{ |
|||
// Build and transform padded sequence b_k = exp(I*Pi*k^2/N)
|
|||
for (int i = 0; i < n; i++) |
|||
{ |
|||
b[i] = sequence[i]; |
|||
} |
|||
|
|||
for (int i = m - n + 1; i < b.Length; i++) |
|||
{ |
|||
b[i] = sequence[m - i]; |
|||
} |
|||
|
|||
Radix2(b, -1); |
|||
}, |
|||
() => |
|||
{ |
|||
// Build and transform padded sequence a_k = x_k * exp(-I*Pi*k^2/N)
|
|||
for (int i = 0; i < samples.Length; i++) |
|||
{ |
|||
a[i] = sequence[i].Conjugate * samples[i]; |
|||
} |
|||
|
|||
Radix2(a, -1); |
|||
}); |
|||
|
|||
for (int i = 0; i < a.Length; i++) |
|||
{ |
|||
a[i] *= b[i]; |
|||
} |
|||
|
|||
Radix2(a, 1); |
|||
|
|||
var nbinv = 1.0 / m; |
|||
for (int i = 0; i < samples.Length; i++) |
|||
{ |
|||
samples[i] = nbinv * sequence[i].Conjugate * a[i]; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Swap the real and imaginary parts of each sample.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample Vector.</param>
|
|||
private static void SwapRealImaginary(Complex[] samples) |
|||
{ |
|||
for (int i = 0; i < samples.Length; i++) |
|||
{ |
|||
samples[i] = Complex.WithRealImaginary(samples[i].Imaginary, samples[i].Real); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bluestein 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>
|
|||
internal static void Bluestein(Complex[] samples, int exponentSign) |
|||
{ |
|||
int n = samples.Length; |
|||
if (n.IsPowerOfTwo()) |
|||
{ |
|||
Radix2(samples, exponentSign); |
|||
return; |
|||
} |
|||
|
|||
if (exponentSign == 1) |
|||
{ |
|||
SwapRealImaginary(samples); |
|||
} |
|||
|
|||
BluesteinConvolution(samples); |
|||
|
|||
if (exponentSign == 1) |
|||
{ |
|||
SwapRealImaginary(samples); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bluestein forward FFT for arbitrary sample vectors.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
public void BluesteinForward(Complex[] samples, FourierOptions options) |
|||
{ |
|||
Bluestein(samples, SignByOptions(options)); |
|||
ForwardScaleByOptions(options, samples); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Bluestein inverse FFT for arbitrary sample vectors.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
public void BluesteinInverse(Complex[] samples, FourierOptions options) |
|||
{ |
|||
Bluestein(samples, -SignByOptions(options)); |
|||
InverseScaleByOptions(options, samples); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,95 @@ |
|||
// <copyright file="DiscreteFourierTransform.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>
|
|||
/// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
|
|||
/// </summary>
|
|||
public partial class DiscreteFourierTransform |
|||
{ |
|||
/// <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) |
|||
{ |
|||
double w0 = exponentSign * 2 * Constants.Pi / samples.Length; |
|||
var spectrum = new Complex[samples.Length]; |
|||
|
|||
Parallel.For( |
|||
0, |
|||
samples.Length, |
|||
k => |
|||
{ |
|||
double wk = w0 * k; |
|||
Complex sum = Complex.Zero; |
|||
for (int n = 0; n < samples.Length; n++) |
|||
{ |
|||
double w = n * wk; |
|||
sum += samples[n] * Complex.WithRealImaginary(Math.Cos(w), Math.Sin(w)); |
|||
} |
|||
|
|||
spectrum[k] = sum; |
|||
}); |
|||
|
|||
return spectrum; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Naive forward DFT, useful e.g. to verify faster algorithms.
|
|||
/// </summary>
|
|||
/// <param name="timeSpace">Time-space sample vector.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
/// <returns>Corresponding frequency-space vector.</returns>
|
|||
public Complex[] NaiveForward(Complex[] timeSpace, FourierOptions options) |
|||
{ |
|||
var frequencySpace = Naive(timeSpace, SignByOptions(options)); |
|||
ForwardScaleByOptions(options, frequencySpace); |
|||
return frequencySpace; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Naive inverse DFT, useful e.g. to verify faster algorithms.
|
|||
/// </summary>
|
|||
/// <param name="frequencySpace">Frequency-space sample vector.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
/// <returns>Corresponding time-space vector.</returns>
|
|||
public Complex[] NaiveInverse(Complex[] frequencySpace, FourierOptions options) |
|||
{ |
|||
var timeSpace = Naive(frequencySpace, -SignByOptions(options)); |
|||
InverseScaleByOptions(options, timeSpace); |
|||
return timeSpace; |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,93 @@ |
|||
// <copyright file="DiscreteFourierTransform.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>
|
|||
/// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
|
|||
/// </summary>
|
|||
public partial class DiscreteFourierTransform |
|||
{ |
|||
/// <summary>
|
|||
/// Extract the exponent sign to be used in forward transforms according to the
|
|||
/// provided convention options.
|
|||
/// </summary>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
/// <returns>Fourier series exponent sign.</returns>
|
|||
private static int SignByOptions(FourierOptions options) |
|||
{ |
|||
return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1; |
|||
} |
|||
|
|||
/// <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(FourierOptions options, Complex[] samples) |
|||
{ |
|||
if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling || |
|||
(options & FourierOptions.AsymmetricScaling) == FourierOptions.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(FourierOptions options, Complex[] samples) |
|||
{ |
|||
if ((options & FourierOptions.NoScaling) == FourierOptions.NoScaling) |
|||
{ |
|||
return; |
|||
} |
|||
|
|||
var scalingFactor = 1.0 / samples.Length; |
|||
if ((options & FourierOptions.AsymmetricScaling) != FourierOptions.AsymmetricScaling) |
|||
{ |
|||
scalingFactor = Math.Sqrt(scalingFactor); |
|||
} |
|||
|
|||
for (int i = 0; i < samples.Length; i++) |
|||
{ |
|||
samples[i] *= scalingFactor; |
|||
} |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,137 @@ |
|||
// <copyright file="DiscreteFourierTransform.RadixN.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 NumberTheory; |
|||
using Properties; |
|||
|
|||
/// <summary>
|
|||
/// Complex Fast (FFT) Implementation of the Discrete Fourier Transform (DFT).
|
|||
/// </summary>
|
|||
public partial class DiscreteFourierTransform |
|||
{ |
|||
/// <summary>
|
|||
/// Radix-2 Reorder Helper Method
|
|||
/// </summary>
|
|||
/// <typeparam name="T">Sample type</typeparam>
|
|||
/// <param name="samples">Sample vector</param>
|
|||
private static void Radix2Reorder<T>(T[] samples) |
|||
{ |
|||
int j = 0; |
|||
for (int i = 0; i < samples.Length - 1; i++) |
|||
{ |
|||
if (i < j) |
|||
{ |
|||
T temp = samples[i]; |
|||
samples[i] = samples[j]; |
|||
samples[j] = temp; |
|||
} |
|||
|
|||
int m = samples.Length; |
|||
|
|||
do |
|||
{ |
|||
m >>= 1; |
|||
j ^= m; |
|||
} while ((j & m) == 0); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Radix-2 Step Helper Method
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector.</param>
|
|||
/// <param name="exponentSign">Fourier series exponent sign.</param>
|
|||
/// <param name="levelSize">Level Group Size.</param>
|
|||
/// <param name="k">Index inside of the level.</param>
|
|||
private static void Radix2Step(Complex[] samples, int exponentSign, int levelSize, int k) |
|||
{ |
|||
// Twiddle Factor
|
|||
double exponent = (exponentSign * k) * Constants.Pi / levelSize; |
|||
Complex w = Complex.WithRealImaginary(Math.Cos(exponent), Math.Sin(exponent)); |
|||
|
|||
int step = levelSize << 1; |
|||
for (int i = k; i < samples.Length; i += step) |
|||
{ |
|||
Complex ai = samples[i]; |
|||
Complex t = w * samples[i + levelSize]; |
|||
samples[i] = ai + t; |
|||
samples[i + levelSize] = ai - t; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Radix-2 generic FFT for power-of-two sample vectors.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
|
|||
/// <param name="exponentSign">Fourier series exponent sign.</param>
|
|||
/// <exception cref="ArgumentException"/>
|
|||
internal static void Radix2(Complex[] samples, int exponentSign) |
|||
{ |
|||
if (!samples.Length.IsPowerOfTwo()) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentPowerOfTwo); |
|||
} |
|||
|
|||
Radix2Reorder(samples); |
|||
for (int levelSize = 1; levelSize < samples.Length; levelSize *= 2) |
|||
{ |
|||
for (int k = 0; k <= levelSize - 1; k++) |
|||
{ |
|||
Radix2Step(samples, exponentSign, levelSize, k); |
|||
} |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Radix-2 forward FFT for power-of-two sample vectors.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
/// <exception cref="ArgumentException"/>
|
|||
public void Radix2Forward(Complex[] samples, FourierOptions options) |
|||
{ |
|||
Radix2(samples, SignByOptions(options)); |
|||
ForwardScaleByOptions(options, samples); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Radix-2 inverse FFT for power-of-two sample vectors.
|
|||
/// </summary>
|
|||
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param>
|
|||
/// <param name="options">Fourier Transform Convention Options.</param>
|
|||
/// <exception cref="ArgumentException"/>
|
|||
public void Radix2Inverse(Complex[] samples, FourierOptions options) |
|||
{ |
|||
Radix2(samples, -SignByOptions(options)); |
|||
InverseScaleByOptions(options, samples); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,73 @@ |
|||
// <copyright file="FourierOptions.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>
|
|||
/// Fourier Transform Convention
|
|||
/// </summary>
|
|||
[Flags] |
|||
public enum FourierOptions |
|||
{ |
|||
// FLAGS:
|
|||
|
|||
/// <summary>
|
|||
/// Inverse integrand exponent (forward: positive sign; inverse: negative sign).
|
|||
/// </summary>
|
|||
InverseExponent = 0x01, |
|||
|
|||
/// <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 and common exponent (used in Maple).
|
|||
/// </summary>
|
|||
Default = 0, |
|||
|
|||
/// <summary>
|
|||
/// Only scale by 1/N in the inverse direction; No scaling in forward direction (used in Matlab). [= AsymmetricScaling]
|
|||
/// </summary>
|
|||
Matlab = AsymmetricScaling, |
|||
|
|||
/// <summary>
|
|||
/// Inverse integrand exponent; No scaling at all (used in all Numerical Recipes based implementations). [= InverseExponent | NoScaling]
|
|||
/// </summary>
|
|||
NumericalRecipes = InverseExponent | NoScaling |
|||
} |
|||
} |
|||
Loading…
Reference in new issue