us1zcljz
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

FFT: clearer internal structure

spatial
Christoph Ruegg 9 years ago
parent
7614be2ee6
commit
1f75cdf7e6
3 changed files with 163 additions and 89 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 77
      src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs
  2. 111
      src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs
  3. 64
      src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs

77
src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Bluestein.cs
View File

@ -28,7 +28,6 @@
using System; using System;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
using MathNet.Numerics.Properties;
using MathNet.Numerics.Threading; using MathNet.Numerics.Threading;
using Complex = System.Numerics.Complex; using Complex = System.Numerics.Complex;
@ -38,9 +37,9 @@ namespace MathNet.Numerics.Providers.FourierTransform
internal partial class ManagedFourierTransformProvider internal partial class ManagedFourierTransformProvider
{ {
/// <summary> /// <summary>
/// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1 /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1
/// will cause k*k in the Bluestein sequence to overflow (GH-286). /// will cause k*k in the Bluestein sequence to overflow (GH-286).
/// </summary> /// </summary>
const int BluesteinSequenceLengthThreshold = 46341; const int BluesteinSequenceLengthThreshold = 46341;
/// <summary> /// <summary>
@ -135,7 +134,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
b[i] = sequence[m - i]; b[i] = sequence[m - i];
} }
Radix2(b, -1); Radix2Forward(b);
}, },
() => () =>
{ {
@ -145,7 +144,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
a[i] = sequence[i].Conjugate() * samples[i]; a[i] = sequence[i].Conjugate() * samples[i];
} }
Radix2(a, -1); Radix2Forward(a);
}); });
for (int i = 0; i < a.Length; i++) for (int i = 0; i < a.Length; i++)
@ -153,7 +152,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
a[i] *= b[i]; a[i] *= b[i];
} }
Radix2Parallel(a, 1); Radix2InverseParallel(a);
var nbinv = 1.0f / m; var nbinv = 1.0f / m;
for (int i = 0; i < samples.Length; i++) for (int i = 0; i < samples.Length; i++)
@ -190,7 +189,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
b[i] = sequence[m - i]; b[i] = sequence[m - i];
} }
Radix2(b, -1); Radix2Forward(b);
}, },
() => () =>
{ {
@ -200,7 +199,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
a[i] = sequence[i].Conjugate() * samples[i]; a[i] = sequence[i].Conjugate() * samples[i];
} }
Radix2(a, -1); Radix2Forward(a);
}); });
for (int i = 0; i < a.Length; i++) for (int i = 0; i < a.Length; i++)
@ -208,7 +207,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
a[i] *= b[i]; a[i] *= b[i];
} }
Radix2Parallel(a, 1); Radix2InverseParallel(a);
var nbinv = 1.0 / m; var nbinv = 1.0 / m;
for (int i = 0; i < samples.Length; i++) for (int i = 0; i < samples.Length; i++)
@ -244,55 +243,37 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// <summary> /// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors. /// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary> /// </summary>
/// <param name="samples">Time-space sample vector.</param> private static void BluesteinForward(Complex[] samples)
/// <param name="exponentSign">Fourier series exponent sign.</param>
private static void Bluestein(Complex32[] samples, int exponentSign)
{ {
int n = samples.Length;
if (n.IsPowerOfTwo())
{
Radix2Parallel(samples, exponentSign);
return;
}
if (exponentSign == 1)
{
SwapRealImaginary(samples);
}
BluesteinConvolutionParallel(samples); BluesteinConvolutionParallel(samples);
if (exponentSign == 1)
{
SwapRealImaginary(samples);
}
} }
/// <summary> /// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors. /// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary> /// </summary>
/// <param name="samples">Time-space sample vector.</param> private static void BluesteinInverse(Complex[] spectrum)
/// <param name="exponentSign">Fourier series exponent sign.</param>
private static void Bluestein(Complex[] samples, int exponentSign)
{ {
int n = samples.Length; SwapRealImaginary(spectrum);
if (n.IsPowerOfTwo()) BluesteinConvolutionParallel(spectrum);
{ SwapRealImaginary(spectrum);
Radix2Parallel(samples, exponentSign); }
return;
}
if (exponentSign == 1)
{
SwapRealImaginary(samples);
}
/// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary>
private static void BluesteinForward(Complex32[] samples)
{
BluesteinConvolutionParallel(samples); BluesteinConvolutionParallel(samples);
}
if (exponentSign == 1) /// <summary>
{ /// Bluestein generic FFT for arbitrary sized sample vectors.
SwapRealImaginary(samples); /// </summary>
} private static void BluesteinInverse(Complex32[] spectrum)
{
SwapRealImaginary(spectrum);
BluesteinConvolutionParallel(spectrum);
SwapRealImaginary(spectrum);
} }
} }
} }

111
src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.Radix2.cs
View File

@ -28,7 +28,6 @@
using System; using System;
using System.Runtime.CompilerServices; using System.Runtime.CompilerServices;
using MathNet.Numerics.Properties;
using MathNet.Numerics.Threading; using MathNet.Numerics.Threading;
using Complex = System.Numerics.Complex; using Complex = System.Numerics.Complex;
@ -120,22 +119,29 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// <summary> /// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors. /// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary> /// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param> private static void Radix2Forward(Complex32[] data)
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
private static void Radix2(Complex32[] samples, int exponentSign)
{ {
if (!samples.Length.IsPowerOfTwo()) Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
throw new ArgumentException(Resources.ArgumentPowerOfTwo); for (var k = 0; k < levelSize; k++)
{
Radix2Step(data, -1, levelSize, k);
}
} }
}
Radix2Reorder(samples); /// <summary>
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) /// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
private static void Radix2Forward(Complex[] data)
{
Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
for (var k = 0; k < levelSize; k++) for (var k = 0; k < levelSize; k++)
{ {
Radix2Step(samples, exponentSign, levelSize, k); Radix2Step(data, -1, levelSize, k);
} }
} }
} }
@ -143,22 +149,29 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// <summary> /// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors. /// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary> /// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param> private static void Radix2Inverse(Complex32[] data)
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
private static void Radix2(Complex[] samples, int exponentSign)
{ {
if (!samples.Length.IsPowerOfTwo()) Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
throw new ArgumentException(Resources.ArgumentPowerOfTwo); for (var k = 0; k < levelSize; k++)
{
Radix2Step(data, 1, levelSize, k);
}
} }
}
Radix2Reorder(samples); /// <summary>
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) /// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
private static void Radix2Inverse(Complex[] data)
{
Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
for (var k = 0; k < levelSize; k++) for (var k = 0; k < levelSize; k++)
{ {
Radix2Step(samples, exponentSign, levelSize, k); Radix2Step(data, 1, levelSize, k);
} }
} }
} }
@ -166,18 +179,30 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// <summary> /// <summary>
/// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary> /// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param> private static void Radix2ForwardParallel(Complex32[] data)
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
private static void Radix2Parallel(Complex32[] samples, int exponentSign)
{ {
if (!samples.Length.IsPowerOfTwo()) Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
throw new ArgumentException(Resources.ArgumentPowerOfTwo); var size = levelSize;
CommonParallel.For(0, size, 64, (u, v) =>
{
for (int i = u; i < v; i++)
{
Radix2Step(data, -1, size, i);
}
});
} }
}
Radix2Reorder(samples); /// <summary>
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary>
private static void Radix2ForwardParallel(Complex[] data)
{
Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
var size = levelSize; var size = levelSize;
@ -185,7 +210,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
{ {
for (int i = u; i < v; i++) for (int i = u; i < v; i++)
{ {
Radix2Step(samples, exponentSign, size, i); Radix2Step(data, -1, size, i);
} }
}); });
} }
@ -194,18 +219,30 @@ namespace MathNet.Numerics.Providers.FourierTransform
/// <summary> /// <summary>
/// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version). /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary> /// </summary>
/// <param name="samples">Sample vector, where the FFT is evaluated in place.</param> private static void Radix2InverseParallel(Complex32[] data)
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <exception cref="ArgumentException"/>
private static void Radix2Parallel(Complex[] samples, int exponentSign)
{ {
if (!samples.Length.IsPowerOfTwo()) Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
throw new ArgumentException(Resources.ArgumentPowerOfTwo); var size = levelSize;
CommonParallel.For(0, size, 64, (u, v) =>
{
for (int i = u; i < v; i++)
{
Radix2Step(data, 1, size, i);
}
});
} }
}
Radix2Reorder(samples); /// <summary>
for (var levelSize = 1; levelSize < samples.Length; levelSize *= 2) /// Radix-2 generic FFT for power-of-two sample vectors (Parallel Version).
/// </summary>
private static void Radix2InverseParallel(Complex[] data)
{
Radix2Reorder(data);
for (var levelSize = 1; levelSize < data.Length; levelSize *= 2)
{ {
var size = levelSize; var size = levelSize;
@ -213,7 +250,7 @@ namespace MathNet.Numerics.Providers.FourierTransform
{ {
for (int i = u; i < v; i++) for (int i = u; i < v; i++)
{ {
Radix2Step(samples, exponentSign, size, i); Radix2Step(data, 1, size, i);
} }
}); });
} }

64
src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs
View File

@ -65,7 +65,21 @@ namespace MathNet.Numerics.Providers.FourierTransform
public void Forward(Complex32[] samples, FourierTransformScaling scaling) public void Forward(Complex32[] samples, FourierTransformScaling scaling)
{ {
Bluestein(samples, -1); if (samples.Length.IsPowerOfTwo())
{
if (samples.Length >= 1024)
{
Radix2ForwardParallel(samples);
}
else
{
Radix2Forward(samples);
}
}
else
{
BluesteinForward(samples);
}
switch (scaling) switch (scaling)
{ {
@ -84,7 +98,21 @@ namespace MathNet.Numerics.Providers.FourierTransform
public void Forward(Complex[] samples, FourierTransformScaling scaling) public void Forward(Complex[] samples, FourierTransformScaling scaling)
{ {
Bluestein(samples, -1); if (samples.Length.IsPowerOfTwo())
{
if (samples.Length >= 1024)
{
Radix2ForwardParallel(samples);
}
else
{
Radix2Forward(samples);
}
}
else
{
BluesteinForward(samples);
}
switch (scaling) switch (scaling)
{ {
@ -103,7 +131,21 @@ namespace MathNet.Numerics.Providers.FourierTransform
public void Backward(Complex32[] spectrum, FourierTransformScaling scaling) public void Backward(Complex32[] spectrum, FourierTransformScaling scaling)
{ {
Bluestein(spectrum, 1); if (spectrum.Length.IsPowerOfTwo())
{
if (spectrum.Length >= 1024)
{
Radix2InverseParallel(spectrum);
}
else
{
Radix2Inverse(spectrum);
}
}
else
{
BluesteinInverse(spectrum);
}
switch (scaling) switch (scaling)
{ {
@ -122,7 +164,21 @@ namespace MathNet.Numerics.Providers.FourierTransform
public void Backward(Complex[] spectrum, FourierTransformScaling scaling) public void Backward(Complex[] spectrum, FourierTransformScaling scaling)
{ {
Bluestein(spectrum, 1); if (spectrum.Length.IsPowerOfTwo())
{
if (spectrum.Length >= 1024)
{
Radix2InverseParallel(spectrum);
}
else
{
Radix2Inverse(spectrum);
}
}
else
{
BluesteinInverse(spectrum);
}
switch (scaling) switch (scaling)
{ {

Write Preview
Loading…
Cancel
Save