使用Intel MKL进行3D卷积

Question

我已经编写了一个使用Intel MKL计算3D卷积的C/C++代码，数组大约有300×200×200个元素。我想应用一个大小为3×3×3或5×5×5的核。输入的3D数组和核都具有实值。

这个3D数组以列优先方式存储为double类型的1D数组。类似地，核也是double类型并以列优先方式保存。例如，

for( int k = 0; k < nk; k++ ) // Loop through the height.
    for( int j = 0; j < nj; j++ ) // Loop through the rows.
        for( int i = 0; i < ni; i++ ) // Loop through the columns.
        {
            ijk = i + ni * j + ni * nj * k;
            my3Darray[ ijk ] = 1.0;
        }

为了进行卷积计算，我想在输入数组和内核上执行非就地FFT，并防止它们被修改（我需要稍后在我的代码中使用它们），然后进行就地反向计算。

当我将我的代码得到的结果与MATLAB得到的结果进行比较时，它们非常不同。是否有人可以帮助我解决这个问题？我的代码缺少什么？

这是我使用的MATLAB代码：

a = ones( 10, 10, 10 );
kernel = ones( 3, 3, 3 );
aconvolved = convn( a, kernel, 'same' );

这是我的C/C++代码：

#include <stdio.h>
#include "mkl.h"

void Conv3D(
    double *in, double *ker, double *out,
    int nRows, int nCols, int nHeights)
{

    int NI = nRows;
    int NJ = nCols;
    int NK = nHeights;

    double *in_fft  = new double [NI*NJ*NK];
    double *ker_fft = new double [NI*NJ*NK];

    DFTI_DESCRIPTOR_HANDLE fft_desc = 0;
    MKL_LONG sizes[]   = { NK, NJ, NI };
    MKL_LONG strides[] = { 0, NJ*NI, NI, 1 };

    DftiCreateDescriptor( &fft_desc, DFTI_DOUBLE, DFTI_REAL, 3, sizes     );
    DftiSetValue        (  fft_desc, DFTI_PLACEMENT     , DFTI_NOT_INPLACE);   // Out-of-place computation.
    DftiSetValue        (  fft_desc, DFTI_INPUT_STRIDES , strides         );
    DftiSetValue        (  fft_desc, DFTI_OUTPUT_STRIDES, strides         );
    DftiSetValue        (  fft_desc, DFTI_BACKWARD_SCALE, 1/NI/NJ/NK      );
    DftiCommitDescriptor(  fft_desc );

    DftiComputeForward  (  fft_desc, in , in_fft  );
    DftiComputeForward  (  fft_desc, ker, ker_fft );

    for (long long i = 0; i < (long long)NI*NJ*NK; ++i )
        out[i] = in_fft[i]*ker_fft[i];

    // In-place computation.
    DftiSetValue        (  fft_desc, DFTI_PLACEMENT, DFTI_INPLACE );
    DftiCommitDescriptor(  fft_desc      );
    DftiComputeBackward (  fft_desc, out );

    DftiFreeDescriptor  ( &fft_desc );

    delete[] in_fft;
    delete[] ker_fft;

}

int main(int argc, char* argv[])
{
    int n = 10;
    int nkernel = 3;

    double *a          = new double [n*n*n]; // This array is real.
    double *aconvolved = new double [n*n*n]; // The convolved array is also real.
    double *kernel     = new double [nkernel*nkernel*nkernel]; // kernel is real.

    // Fill the array with some 'real' numbers.
    for( int i = 0; i < n*n*n; i++ )
        a[ i ] = 1.0;

    // Fill the kernel with some 'real' numbers.
    for( int i = 0; i < nkernel*nkernel*nkernel; i++ )
        kernel[ i ] = 1.0;

    // Calculate the convolution.
    Conv3D( a, kernel, aconvolved, n, n, n );

    printf("Convolved:\n");
    for( int i = 0; i < n*n*n; i++ )
        printf( "%15.8f\n", aconvolved[i] );

    delete[] a;
    delete[] kernel;
    delete[] aconvolved;

    return 0;
}

- afp_2008

---

我喜欢你的问题，但是我无法帮助你。然而，如果你可以在你的结果旁边提供期望的结果，那么其他人可能会更容易帮助你。 - nonsensickle

Answer 1

如果使用实值频率数据进行FFT，就无法进行逆向FFT（只能获取幅度）。正向FFT需要输出复杂数据。这可以通过将DFTI_FORWARD_DOMAIN设置为DFTI_COMPLEX来完成。

DftiCreateDescriptor( &fft_desc, DFTI_DOUBLE, DFTI_COMPLEX, 3, sizes     );

这样做会隐式地将后向域设置为复杂。

您还需要一个复杂的数据类型。可能是这样的：

MKL_Complex16* in_fft  = new MKL_Complex16[NI*NJ*NK];

这意味着你需要将实部和虚部都乘以相应的值：

for (size_t i = 0; i < (size_t)NI*NJ*NK; ++i) {
    out_fft[i].real = in_fft[i].real * ker_fft[i].real;
    out_fft[i].imag = in_fft[i].imag * ker_fft[i].imag;
}

反向FFT的输出也是复数，假设您的输入数据是实数，则可以直接获取`.real`组件作为结果。这意味着您需要一个临时的复杂数组（如上面的`out_fft`）。

另外请注意，为了避免伪像，您希望每个维度上的fft大小至少为M+N-1。通常，您会选择下一个最高的2的幂以提高速度。

我强烈建议您首先在MATLAB中使用FFT实现它。有许多这样的实现可用（例如example），但我建议您从基础开始自己制作一个简单的函数。

- chappjc