如何在C#中实现类似MATLAB fftn()的3D-FFT与3D-IFFT？

我需要在C#中实现和MATLAB fftn() 功能一致的3D-FFT与3D-IFFT，但找不到可用的免费类库，恳请帮助。

我试过Math.net、WFFT.Net，都无法正常使用；ILNumerics可行但并非免费软件。自己通过1D-FFT拼接实现了3D-FFT和3D-IFFT算法，但速度太慢：处理64×64×1000点的3D-FFT耗时2秒，而MATLAB仅需0.2秒。以下是我的实现代码：

static Complex[,,] Fft3D(float[,,] data)
{
    int M = data.GetLength(0);
    int N = data.GetLength(1);
    int K = data.GetLength(2);
    Complex[,,] data3D = new Complex[M, N, K];
    for (int i = 0; i < M; i++)
    {
        for (int j = 0; j < N; j++)
        {
            for (int k = 0; k < K; k++)
            {
                data3D[i, j, k] = new Complex(data[i, j, k], 0.0);
            }
        }
    }          
    for (int i = 0; i < M; i++)
    {
        for (int j = 0; j < N; j++)
        {                   
            Complex[] data1D = new Complex[K];
            for(int k = 0; k < K; k++)
            {
                data1D[k] = data3D[i, j, k];
            }
            Fourier.Forward(data1D, FourierOptions.Matlab);
            for (int k = 0; k < K; k++)
            {
                data3D[i, j, k] = data1D[k];
            }
        }
    }
    for (int i = 0; i < M; i++)
    {
        for (int k = 0; k < K; k++)
        {
            Complex[] data1D = new Complex[N];
            for (int j = 0; j < N; j++)
            {
                data1D[j] = data3D[i, j, k];
            }
            Fourier.Forward(data1D, FourierOptions.Matlab);
            for (int j = 0; j < N; j++)
            {
                data3D[i, j, k] = data1D[j];
            }
        }
    }
    for (int j = 0; j < N; j++)
    {
        for (int k = 0; k < K; k++)
        {
            Complex[] data1D = new Complex[M];
            for (int i = 0; i < M; i++)
            {
                data1D[i] = data3D[i, j, k];
            }
            Fourier.Forward(data1D, FourierOptions.Matlab);
            for (int i = 0; i < M; i++)
            {
                data3D[i, j, k] = data1D[i];
            }
        }
    }
    return data3D;
}
static float[,,] Ifft3D(Complex[,,] data)
{
    int M = data.GetLength(0);
    int N = data.GetLength(1);
    int K = data.GetLength(2);
    float[,,] data3D = new float[M, N, K];
    Complex[,,] dataCopy = (Complex[,,])data.Clone();
    for (int i = 0; i < M; i++)
    {
        for (int j = 0; j < N; j++)
        {
            Complex[] data1D = new Complex[K];
            for (int k = 0; k < K; k++)
            {
                data1D[k] = dataCopy[i, j, k];
            }
            Fourier.Inverse(data1D, FourierOptions.Matlab);
            for (int k = 0; k < K; k++)
            {
                dataCopy[i, j, k] = data1D[k];
            }
        }
    }
    for (int i = 0; i < M; i++)
    {
        for (int k = 0; k < K; k++)
        {
            Complex[] data1D = new Complex[N];
            for (int j = 0; j < N; j++)
            {
                data1D[j] = dataCopy[i, j, k];
            }
            Fourier.Inverse(data1D, FourierOptions.Matlab);
            for (int j = 0; j < N; j++)
            {
                dataCopy[i, j, k] = data1D[j];
            }
        }
    }
    for (int j = 0; j < N; j++)
    {
        for (int k = 0; k < K; k++)
        {
            Complex[] data1D = new Complex[M];
            for (int i = 0; i < M; i++)
            {
                data1D[i] = dataCopy[i, j, k];
            }
            Fourier.Inverse(data1D, FourierOptions.Matlab);
            for (int i = 0; i < M; i++)
            {
                dataCopy[i, j, k] = data1D[i];
            }
        }
    }
    for (int i = 0; i < M; i++)
    {
        for (int j = 0; j < N; j++)
        {
            for (int k = 0; k < K; k++)
            {
                data3D[i, j, k] = (float)dataCopy[i, j, k].Real;
            }
        }
    }
    return data3D;
}

内容的提问来源于stack exchange，提问作者ZhouY