2
我想計算3D FFT使用Intel MKL的數組，其具有約300×200×200元素。這種3D數組存儲爲double類型的縱列方式一維數組：使用帶有零填充的英特爾MKL的3D FFT
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;
}
我要輸入陣列上執行not-in-place FFT和防止屏幕修改（我需要在我的代碼後使用它）然後進行反向計算in-place。我也想要零填充。
我的問題是：
- 我怎麼能執行零填充？
- 在計算中包含零填充時，我應該如何處理
FFT函數使用的數組大小？ - 如何取出零填充結果並獲得實際結果？
這裏是我對這個問題的嘗試，我會絕對感謝的任何評論，建議或提示。
#include <stdio.h>
#include "mkl.h"
int max(int a, int b, int c)
{
int m = a;
(m < b) && (m = b);
(m < c) && (m = c);
return m;
}
void FFT3D_R2C(// Real to Complex 3D FFT.
double *in, int nRowsIn , int nColsIn , int nHeightsIn ,
double *out)
{
int n = max(nRowsIn , nColsIn , nHeightsIn );
// Round up to the next highest power of 2.
unsigned int N = (unsigned int) n; // compute the next highest power of 2 of 32-bit n.
N--;
N |= N >> 1;
N |= N >> 2;
N |= N >> 4;
N |= N >> 8;
N |= N >> 16;
N++;
/* Strides describe data layout in real and conjugate-even domain. */
MKL_LONG rs[4], cs[4];
// DFTI descriptor.
DFTI_DESCRIPTOR_HANDLE fft_desc = 0;
// Variables needed for out-of-place computations.
MKL_Complex16 *in_fft = new MKL_Complex16 [ N*N*N ];
MKL_Complex16 *out_fft = new MKL_Complex16 [ N*N*N ];
double *out_ZeroPadded = new double [ N*N*N ];
/* Compute strides */
rs[3] = 1; cs[3] = 1;
rs[2] = (N/2+1)*2; cs[2] = (N/2+1);
rs[1] = N*(N/2+1)*2; cs[1] = N*(N/2+1);
rs[0] = 0; cs[0] = 0;
// Create DFTI descriptor.
MKL_LONG sizes[] = { N, N, N };
DftiCreateDescriptor(&fft_desc, DFTI_DOUBLE, DFTI_REAL, 3, sizes);
// Configure DFTI descriptor.
DftiSetValue(fft_desc, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
DftiSetValue(fft_desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE ); // Out-of-place transformation.
DftiSetValue(fft_desc, DFTI_INPUT_STRIDES , rs );
DftiSetValue(fft_desc, DFTI_OUTPUT_STRIDES , cs );
DftiCommitDescriptor(fft_desc);
DftiComputeForward (fft_desc, in , in_fft );
// Change strides to compute backward transform.
DftiSetValue (fft_desc, DFTI_INPUT_STRIDES , cs);
DftiSetValue (fft_desc, DFTI_OUTPUT_STRIDES, rs);
DftiCommitDescriptor(fft_desc);
DftiComputeBackward (fft_desc, out_fft, out_ZeroPadded);
// Printing the zero padded 3D FFT result.
for(long long i = 0; i < (long long)N*N*N; i++)
printf("%f\n", out_ZeroPadded[i]);
/* I don't know how to take out the zero padded results and
save the actual result in the variable named "out" */
DftiFreeDescriptor (&fft_desc);
delete[] in_fft;
delete[] out_ZeroPadded ;
}
int main()
{
int n = 10;
double *a = new double [n*n*n]; // This array is real.
double *afft = new double [n*n*n];
// Fill the array with some 'real' numbers.
for(int i = 0; i < n*n*n; i++)
a[ i ] = 1.0;
// Calculate FFT.
FFT3D_R2C(a, n, n, n, afft);
printf("FFT results:\n");
for(int i = 0; i < n*n*n; i++)
printf("%15.8f\n", afft[i]);
delete[] a;
delete[] afft;
return 0;
}