Skip to content

Implement managed ForwardMultidim using separable 1D FFT #1151

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bryancostanich wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Implement managed ForwardMultidim using separable 1D FFT #1151

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 113 additions & 0 deletions src/Numerics.Tests/IntegralTransformsTests/FourierTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,5 +167,118 @@ public void FourierDefaultTransformsRealSineCorrectly64()
}
}
}
[Test]
public void Forward2DAndInverse2DRoundtrip64()
{
const int rows = 4;
const int cols = 8;
var original = new Complex[rows * cols];
var rng = new System.Random(42);
for (int i = 0; i < original.Length; i++)
{
original[i] = new Complex(rng.NextDouble(), rng.NextDouble());
}

var samples = (Complex[])original.Clone();

Fourier.Forward2D(samples, rows, cols, FourierOptions.Default);

// Verify the transform actually changed the data
bool anyDifferent = false;
for (int i = 0; i < samples.Length; i++)
{
if ((samples[i] - original[i]).Magnitude > 1e-10)
{
anyDifferent = true;
break;
}
}
Assert.IsTrue(anyDifferent, "Forward2D should modify the data");

Fourier.Inverse2D(samples, rows, cols, FourierOptions.Default);

// Verify roundtrip recovers original data
for (int i = 0; i < samples.Length; i++)
{
Assert.AreEqual(original[i].Real, samples[i].Real, 1e-10, $"Real part mismatch at index {i}");
Assert.AreEqual(original[i].Imaginary, samples[i].Imaginary, 1e-10, $"Imaginary part mismatch at index {i}");
}
}

[Test]
public void Forward2DAndInverse2DRoundtrip32()
{
const int rows = 4;
const int cols = 8;
var original = new Complex32[rows * cols];
var rng = new System.Random(42);
for (int i = 0; i < original.Length; i++)
{
original[i] = new Complex32((float)rng.NextDouble(), (float)rng.NextDouble());
}

var samples = (Complex32[])original.Clone();

Fourier.Forward2D(samples, rows, cols, FourierOptions.Default);
Fourier.Inverse2D(samples, rows, cols, FourierOptions.Default);

// Verify roundtrip recovers original data
for (int i = 0; i < samples.Length; i++)
{
Assert.AreEqual(original[i].Real, samples[i].Real, 1e-4, $"Real part mismatch at index {i}");
Assert.AreEqual(original[i].Imaginary, samples[i].Imaginary, 1e-4, $"Imaginary part mismatch at index {i}");
}
}

[Test]
public void Forward2DKnownDCValue64()
{
// A constant 4x4 array should produce all energy in the DC bin
const int rows = 4;
const int cols = 4;
var samples = new Complex[rows * cols];
for (int i = 0; i < samples.Length; i++)
{
samples[i] = new Complex(1.0, 0.0);
}

Fourier.Forward2D(samples, rows, cols, FourierOptions.NoScaling);

// DC component should equal rows*cols = 16
Assert.AreEqual(16.0, samples[0].Real, 1e-10, "DC real");
Assert.AreEqual(0.0, samples[0].Imaginary, 1e-10, "DC imag");

// All other components should be zero
for (int i = 1; i < samples.Length; i++)
{
Assert.AreEqual(0.0, samples[i].Real, 1e-10, $"Non-DC real at {i}");
Assert.AreEqual(0.0, samples[i].Imaginary, 1e-10, $"Non-DC imag at {i}");
}
}

[Test]
public void Forward2DNonPowerOfTwoDimensions64()
{
// Test with non-power-of-two dimensions to exercise Bluestein path
const int rows = 3;
const int cols = 5;
var original = new Complex[rows * cols];
var rng = new System.Random(123);
for (int i = 0; i < original.Length; i++)
{
original[i] = new Complex(rng.NextDouble(), rng.NextDouble());
}

var samples = (Complex[])original.Clone();

Fourier.Forward2D(samples, rows, cols, FourierOptions.Default);
Fourier.Inverse2D(samples, rows, cols, FourierOptions.Default);

for (int i = 0; i < samples.Length; i++)
{
Assert.AreEqual(original[i].Real, samples[i].Real, 1e-10, $"Real part mismatch at index {i}");
Assert.AreEqual(original[i].Imaginary, samples[i].Imaginary, 1e-10, $"Imaginary part mismatch at index {i}");
}
}
}
}
184 changes: 180 additions & 4 deletions src/Numerics/Providers/FourierTransform/ManagedFourierTransformProvider.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -320,22 +320,198 @@ public void BackwardReal(double[] spectrum, int n, FourierTransformScaling scali

public void ForwardMultidim(Complex32[] samples, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
TransformMultidim(samples, dimensions, true);
ApplyMultidimScaling(samples, scaling, true);
}

public void ForwardMultidim(Complex[] samples, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
TransformMultidim(samples, dimensions, true);
ApplyMultidimScaling(samples, scaling, true);
}

public void BackwardMultidim(Complex32[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
TransformMultidim(spectrum, dimensions, false);
ApplyMultidimScaling(spectrum, scaling, false);
}

public void BackwardMultidim(Complex[] spectrum, int[] dimensions, FourierTransformScaling scaling)
{
throw new NotSupportedException();
TransformMultidim(spectrum, dimensions, false);
ApplyMultidimScaling(spectrum, scaling, false);
}

/// <summary>
/// Apply separable 1D FFT along each dimension of a row-major N-D array.
/// </summary>
static void TransformMultidim(Complex32[] samples, int[] dimensions, bool forward)
{
int totalLength = samples.Length;

for (int d = 0; d < dimensions.Length; d++)
{
int dimLength = dimensions[d];
int innerLength = 1;
for (int k = d + 1; k < dimensions.Length; k++)
{
innerLength *= dimensions[k];
}

int outerLength = totalLength / (dimLength * innerLength);

var slice = new Complex32[dimLength];

for (int outer = 0; outer < outerLength; outer++)
{
for (int inner = 0; inner < innerLength; inner++)
{
int baseIndex = outer * dimLength * innerLength + inner;

// Extract slice along dimension d
for (int i = 0; i < dimLength; i++)
{
slice[i] = samples[baseIndex + i * innerLength];
}

// Apply 1D FFT (no scaling - we scale once at the end)
if (forward)
{
if (dimLength.IsPowerOfTwo())
{
if (dimLength >= 1024)
Radix2ForwardParallel(slice);
else
Radix2Forward(slice);
}
else
{
BluesteinForward(slice);
}
}
else
{
if (dimLength.IsPowerOfTwo())
{
if (dimLength >= 1024)
Radix2InverseParallel(slice);
else
Radix2Inverse(slice);
}
else
{
BluesteinInverse(slice);
}
}

// Write back
for (int i = 0; i < dimLength; i++)
{
samples[baseIndex + i * innerLength] = slice[i];
}
}
}
}
}

/// <summary>
/// Apply separable 1D FFT along each dimension of a row-major N-D array.
/// </summary>
static void TransformMultidim(Complex[] samples, int[] dimensions, bool forward)
{
int totalLength = samples.Length;

for (int d = 0; d < dimensions.Length; d++)
{
int dimLength = dimensions[d];
int innerLength = 1;
for (int k = d + 1; k < dimensions.Length; k++)
{
innerLength *= dimensions[k];
}

int outerLength = totalLength / (dimLength * innerLength);

var slice = new Complex[dimLength];

for (int outer = 0; outer < outerLength; outer++)
{
for (int inner = 0; inner < innerLength; inner++)
{
int baseIndex = outer * dimLength * innerLength + inner;

// Extract slice along dimension d
for (int i = 0; i < dimLength; i++)
{
slice[i] = samples[baseIndex + i * innerLength];
}

// Apply 1D FFT (no scaling - we scale once at the end)
if (forward)
{
if (dimLength.IsPowerOfTwo())
{
if (dimLength >= 1024)
Radix2ForwardParallel(slice);
else
Radix2Forward(slice);
}
else
{
BluesteinForward(slice);
}
}
else
{
if (dimLength.IsPowerOfTwo())
{
if (dimLength >= 1024)
Radix2InverseParallel(slice);
else
Radix2Inverse(slice);
}
else
{
BluesteinInverse(slice);
}
}

// Write back
for (int i = 0; i < dimLength; i++)
{
samples[baseIndex + i * innerLength] = slice[i];
}
}
}
}
}

static void ApplyMultidimScaling(Complex32[] samples, FourierTransformScaling scaling, bool forward)
{
switch (scaling)
{
case FourierTransformScaling.SymmetricScaling:
HalfRescale(samples);
break;
case FourierTransformScaling.ForwardScaling when forward:
case FourierTransformScaling.BackwardScaling when !forward:
FullRescale(samples);
break;
}
}

static void ApplyMultidimScaling(Complex[] samples, FourierTransformScaling scaling, bool forward)
{
switch (scaling)
{
case FourierTransformScaling.SymmetricScaling:
HalfRescale(samples);
break;
case FourierTransformScaling.ForwardScaling when forward:
case FourierTransformScaling.BackwardScaling when !forward:
FullRescale(samples);
break;
}
}
}
}