diff --git a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
index 820c88e3..192240f7 100644
--- a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
+++ b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
@@ -107,7 +107,7 @@ namespace MathNet.Numerics.IntegralTransforms
b[i] = sequence[m - i];
}
- Radix2(b, -1);
+ Radix2(b, false);
},
() =>
{
@@ -117,7 +117,7 @@ namespace MathNet.Numerics.IntegralTransforms
a[i] = sequence[i].Conjugate()*samples[i];
}
- Radix2(a, -1);
+ Radix2(a, false);
});
for (int i = 0; i < a.Length; i++)
@@ -125,7 +125,7 @@ namespace MathNet.Numerics.IntegralTransforms
a[i] *= b[i];
}
- Radix2Parallel(a, 1);
+ Radix2Parallel(a, true);
var nbinv = 1.0/m;
for (int i = 0; i < samples.Length; i++)
@@ -150,24 +150,24 @@ namespace MathNet.Numerics.IntegralTransforms
/// Bluestein generic FFT for arbitrary sized sample vectors.
///
/// Time-space sample vector.
- /// Fourier series exponent sign.
- internal static void Bluestein(Complex[] samples, int exponentSign)
+ /// Fourier series exponent sign: true for positive, false for negative.
+ internal static void Bluestein(Complex[] samples, bool positiveExponentSign)
{
int n = samples.Length;
if (n.IsPowerOfTwo())
{
- Radix2Parallel(samples, exponentSign);
+ Radix2Parallel(samples, positiveExponentSign);
return;
}
- if (exponentSign == 1)
+ if (positiveExponentSign)
{
SwapRealImaginary(samples);
}
BluesteinConvolutionParallel(samples);
- if (exponentSign == 1)
+ if (positiveExponentSign)
{
SwapRealImaginary(samples);
}
diff --git a/src/Numerics/IntegralTransforms/Fourier.Naive.cs b/src/Numerics/IntegralTransforms/Fourier.Naive.cs
index 609abbde..0e768201 100644
--- a/src/Numerics/IntegralTransforms/Fourier.Naive.cs
+++ b/src/Numerics/IntegralTransforms/Fourier.Naive.cs
@@ -45,11 +45,11 @@ namespace MathNet.Numerics.IntegralTransforms
/// Naive generic DFT, useful e.g. to verify faster algorithms.
///
/// Time-space sample vector.
- /// Fourier series exponent sign.
+ /// Fourier series exponent sign: true for positive, false for negative.
/// Corresponding frequency-space vector.
- internal static Complex[] Naive(Complex[] samples, int exponentSign)
+ internal static Complex[] Naive(Complex[] samples, bool positiveExponentSign)
{
- var w0 = exponentSign*Constants.Pi2/samples.Length;
+ var w0 = positiveExponentSign ? Constants.Pi2/samples.Length : -Constants.Pi2/samples.Length;
var spectrum = new Complex[samples.Length];
CommonParallel.For(0, samples.Length, (u, v) =>
diff --git a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs
index 98564ee5..e075a589 100644
--- a/src/Numerics/IntegralTransforms/Fourier.RadixN.cs
+++ b/src/Numerics/IntegralTransforms/Fourier.RadixN.cs
@@ -3,7 +3,7 @@
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
-// Copyright (c) 2009-2014 Math.NET
+// Copyright (c) 2009-2016 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@@ -97,9 +97,9 @@ namespace MathNet.Numerics.IntegralTransforms
/// Radix-2 generic FFT for power-of-two sized sample vectors.
///
/// Sample vector, where the FFT is evaluated in place.
- /// Fourier series exponent sign.
+ /// Fourier series exponent sign: true for positive, false for negative.
///
- internal static void Radix2(Complex[] samples, int exponentSign)
+ internal static void Radix2(Complex[] samples, bool positiveExponentSign)
{
if (!samples.Length.IsPowerOfTwo())
{
@@ -107,11 +107,21 @@ namespace MathNet.Numerics.IntegralTransforms
}
Radix2Reorder(samples);
- for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
+
+ for (int level = 0, levelSize = 1; levelSize < samples.Length; level++, levelSize *= 2)
{
- for (var k = 0; k < levelSize; k++)
+ Complex[] twiddles = TwiddleFactors(level);
+ for (int k = 0; k < levelSize; k++)
{
- Radix2Step(samples, exponentSign, levelSize, k);
+ Complex twiddle = positiveExponentSign ? twiddles[k] : twiddles[k].Conjugate();
+ int step = levelSize << 1;
+ for (int i = k; i < samples.Length; i += step)
+ {
+ Complex ai = samples[i];
+ Complex t = twiddle * samples[i + levelSize];
+ samples[i] = ai + t;
+ samples[i + levelSize] = ai - t;
+ }
}
}
}
@@ -120,9 +130,9 @@ namespace MathNet.Numerics.IntegralTransforms
/// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
///
/// Sample vector, where the FFT is evaluated in place.
- /// Fourier series exponent sign.
+ /// Fourier series exponent sign: true for positive, false for negative.
///
- internal static void Radix2Parallel(Complex[] samples, int exponentSign)
+ internal static void Radix2Parallel(Complex[] samples, bool positiveExponentSign)
{
if (!samples.Length.IsPowerOfTwo())
{
@@ -130,18 +140,52 @@ namespace MathNet.Numerics.IntegralTransforms
}
Radix2Reorder(samples);
- for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2)
+
+ for (int level = 0, levelSize = 1; levelSize < samples.Length; level++, levelSize *= 2)
{
- var size = levelSize;
+ Complex[] twiddles = TwiddleFactors(level);
+ var step = levelSize << 1;
+
+ int size = levelSize;
CommonParallel.For(0, size, 64, (u, v) =>
{
- for (int i = u; i < v; i++)
+ for (int k = u; k < v; k++)
{
- Radix2Step(samples, exponentSign, size, i);
+ Complex twiddle = twiddles[k];
+ double tr = twiddle.Real;
+ double ti = positiveExponentSign ? twiddle.Imaginary : -twiddle.Imaginary;
+ for (var i = k; i < samples.Length; i += step)
+ {
+ var ai = samples[i];
+ var b = samples[i + size];
+ double ttr = tr * b.Real - ti * b.Imaginary;
+ double tti = ti * b.Real + tr * b.Imaginary;
+ samples[i] = new Complex(ai.Real + ttr, ai.Imaginary + tti);
+ samples[i + size] = new Complex(ai.Real - ttr, ai.Imaginary - tti);
+ }
}
});
}
}
+
+ static readonly Complex[][] Twiddle = new Complex[64][];
+
+ static Complex[] TwiddleFactors(int level)
+ {
+ Complex[] tf = Twiddle[level];
+ if (tf == null)
+ {
+ tf = new Complex[level.PowerOfTwo()];
+ for (int i = 0; i < tf.Length; i++)
+ {
+ double exponent = i * Constants.Pi / tf.Length;
+ tf[i] = new Complex(Math.Cos(exponent), Math.Sin(exponent));
+ }
+ Twiddle[level] = tf;
+ }
+
+ return tf;
+ }
}
}
diff --git a/src/Numerics/IntegralTransforms/Fourier.cs b/src/Numerics/IntegralTransforms/Fourier.cs
index 2364f9ba..ada232dc 100644
--- a/src/Numerics/IntegralTransforms/Fourier.cs
+++ b/src/Numerics/IntegralTransforms/Fourier.cs
@@ -392,7 +392,7 @@ namespace MathNet.Numerics.IntegralTransforms
/// Corresponding frequency-space vector.
public static Complex[] NaiveForward(Complex[] samples, FourierOptions options = FourierOptions.Default)
{
- var frequencySpace = Naive(samples, SignByOptions(options));
+ var frequencySpace = Naive(samples, PositiveSignByOptions(options));
ForwardScaleByOptions(options, frequencySpace);
return frequencySpace;
}
@@ -405,7 +405,7 @@ namespace MathNet.Numerics.IntegralTransforms
/// Corresponding time-space vector.
public static Complex[] NaiveInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
{
- var timeSpace = Naive(spectrum, -SignByOptions(options));
+ var timeSpace = Naive(spectrum, !PositiveSignByOptions(options));
InverseScaleByOptions(options, timeSpace);
return timeSpace;
}
@@ -418,7 +418,7 @@ namespace MathNet.Numerics.IntegralTransforms
///
public static void Radix2Forward(Complex[] samples, FourierOptions options = FourierOptions.Default)
{
- Radix2Parallel(samples, SignByOptions(options));
+ Radix2Parallel(samples, PositiveSignByOptions(options));
ForwardScaleByOptions(options, samples);
}
@@ -430,7 +430,7 @@ namespace MathNet.Numerics.IntegralTransforms
///
public static void Radix2Inverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
{
- Radix2Parallel(spectrum, -SignByOptions(options));
+ Radix2Parallel(spectrum, !PositiveSignByOptions(options));
InverseScaleByOptions(options, spectrum);
}
@@ -441,7 +441,7 @@ namespace MathNet.Numerics.IntegralTransforms
/// Fourier Transform Convention Options.
public static void BluesteinForward(Complex[] samples, FourierOptions options = FourierOptions.Default)
{
- Bluestein(samples, SignByOptions(options));
+ Bluestein(samples, PositiveSignByOptions(options));
ForwardScaleByOptions(options, samples);
}
@@ -452,7 +452,7 @@ namespace MathNet.Numerics.IntegralTransforms
/// Fourier Transform Convention Options.
public static void BluesteinInverse(Complex[] spectrum, FourierOptions options = FourierOptions.Default)
{
- Bluestein(spectrum, -SignByOptions(options));
+ Bluestein(spectrum, !PositiveSignByOptions(options));
InverseScaleByOptions(options, spectrum);
}
@@ -462,9 +462,9 @@ namespace MathNet.Numerics.IntegralTransforms
///
/// Fourier Transform Convention Options.
/// Fourier series exponent sign.
- static int SignByOptions(FourierOptions options)
+ static bool PositiveSignByOptions(FourierOptions options)
{
- return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent ? 1 : -1;
+ return (options & FourierOptions.InverseExponent) == FourierOptions.InverseExponent;
}
///
diff --git a/src/UnitTests/IntegralTransformsTests/FourierTest.cs b/src/UnitTests/IntegralTransformsTests/FourierTest.cs
index 17443794..3f210364 100644
--- a/src/UnitTests/IntegralTransformsTests/FourierTest.cs
+++ b/src/UnitTests/IntegralTransformsTests/FourierTest.cs
@@ -98,8 +98,8 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Forward(samples, FourierOptions.Default));
Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Inverse(samples, FourierOptions.Default));
- Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2(samples, -1));
- Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, -1));
+ Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2(samples, false));
+ Assert.Throws(typeof (ArgumentException), () => Fourier.Radix2Parallel(samples, false));
}
}
}