Search code examples
c#mathfftcomplex-numbersmanaged-cuda

Is it normal for complex array fft-ifft pair radically change values on each iteration?


So simple 3D fft-ifft pair on array of complex values code here:

using System;
using System.Diagnostics;
using System.Numerics;
using System.Runtime.InteropServices;
using ManagedCuda.BasicTypes;
using ManagedCuda;
using ManagedCuda.CudaFFT;
using ManagedCuda.VectorTypes;

class ManagedCuda3DComplexFFT
{
    public static class Extensions
    {
        public delegate void ActOn3DArrayElement(ref cuFloatComplex element, int i, int j, int k);

        public static void ForEach(ref cuFloatComplex[,,] arr, ActOn3DArrayElement act)
        {
            for (int i = 0; i < arr.GetLength(0); i++)
            {
                for (int j = 0; j < arr.GetLength(1); j++)
                {
                    for (int k = 0; k < arr.GetLength(2); k++)
                    {
                        act(ref arr[i, j, k], i, j, k);
                    }
                }
            }   
        }
        public static string ToString(ref cuFloatComplex[,,] arr)
        {
            string result = "";
            for (int i = 0; i < arr.GetLength(0); i++)
            {
                result += "\n" + i.ToString() + " {[\n";
                for (int j = 0; j < arr.GetLength(1); j++)
                {
                    result +=  "[";
                    for (int k = 0; k < arr.GetLength(2); k++)
                    {
                        result += arr[i, j, k].ToString() + "; ";
                    }
                    result += "]\n";
                }
                result += "]};";
            }

            return result;
        }
    }


    static void Main(string[] args)
    {
        const int DIM = 2;

        CudaContext ctx = new CudaContext(0);

        var dx = new CudaDeviceVariable<cuFloatComplex>( DIM*DIM*DIM);
        var data = new cuFloatComplex[DIM,DIM,DIM];
        var fftdData = new cuFloatComplex[DIM,DIM,DIM];
        var ifftdData = new cuFloatComplex[DIM,DIM,DIM];
        Extensions.ForEach(ref data, (ref cuFloatComplex element, int i, int j, int k) =>
        {
            element.real = 1;
            element.imag = 1;
        } );
        Console.WriteLine("source: " + Extensions.ToString(ref data));
        dx.CopyToDevice(data);
        var plan = new CudaFFTPlan3D(DIM, DIM, DIM, cufftType.C2C);
        for (int z = 0; z < 3; z++)
        {
            Console.WriteLine("iteration == " + z);
            plan.Exec(dx.DevicePointer, TransformDirection.Forward); // now inside dx.DevicePointer fftd data 
            dx.CopyToHost(fftdData);
            plan.Exec(dx.DevicePointer, TransformDirection.Inverse); // now inside dx.DevicePointer ifftd data 
            dx.CopyToHost(ifftdData);

            Console.WriteLine("fft: " + Extensions.ToString(ref fftdData));
            Console.WriteLine("ifft: " + Extensions.ToString(ref ifftdData));
            Console.WriteLine("==========================");

        }

        return;
}

to return changing values example output:

source: 
0 {[
[1 + 1i; 1 + 1i; ]
[1 + 1i; 1 + 1i; ]
]};
1 {[
[1 + 1i; 1 + 1i; ]
[1 + 1i; 1 + 1i; ]
]};
iteration == 0
fft: 
0 {[
[8 + 8i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
1 {[
[0 + 0i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
ifft: 
0 {[
[8 + 8i; 8 + 8i; ]
[8 + 8i; 8 + 8i; ]
]};
1 {[
[8 + 8i; 8 + 8i; ]
[8 + 8i; 8 + 8i; ]
]};
==========================
iteration == 1
fft: 
0 {[
[64 + 64i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
1 {[
[0 + 0i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
ifft: 
0 {[
[64 + 64i; 64 + 64i; ]
[64 + 64i; 64 + 64i; ]
]};
1 {[
[64 + 64i; 64 + 64i; ]
[64 + 64i; 64 + 64i; ]
]};
==========================
iteration == 2
fft: 
0 {[
[512 + 512i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
1 {[
[0 + 0i; 0 + 0i; ]
[0 + 0i; 0 + 0i; ]
]};
ifft: 
0 {[
[512 + 512i; 512 + 512i; ]
[512 + 512i; 512 + 512i; ]
]};
1 {[
[512 + 512i; 512 + 512i; ]
[512 + 512i; 512 + 512i; ]
]};
==========================

Solution

  • You have a non-normalized FFT/IFFT pair. Each iteration the values of the array grow by a factor of N (the number of samples in the array).

    Usually it is the inverse transform that has a 1/N factor in it, but some implementations put that factor in the forward transform, or split it by adding a factor 1/sqrt(N) to both transforms.

    Some implementations just want to be fast and leave out the normalization to save N multiplications. You seem to have one of these implementations. You will have to add the normalization yourself.