A simple 2D fft computation question.

A simple 2D fft computation question.

Hi All,

I just begin to try to use MKL for my 2D FFT computation. But I am not able to figure it out with the provided sample code. My code is very simple (as following). Loading data from a file, then I have a 8x8 2D matrix.But the result is not the same as the outcome of Matlab.

Can someone give me some help why my code is incorrect?

Thanks

int main(int args, char *argv[])
{
int m = 8+2, n = 8+2; //8 by 8 matrix
float *input = (float *)malloc(sizeof(float) * m * n); // 10 = (8/2+1)*2
int i = 0, j = 0;
DFTI_DESCRIPTOR *my_desc1_handle;
double data = 0.0;
long status;
long l[2];
l[0] = 8; l[1] = 8;

//Load data
FILE *stream_input;

if( (stream_input = fopen( "C:\temp\im_cp_MKL_part_fft.txt", "r" )) == NULL ) //data file
return -1;

memset(input, 0, sizeof(float) * m * n); //reset the memory to be all zeros.

//obtain the input data
for(i = 0; i < (m-2)*(n-2); i++)
{
fscanf(stream_input,"%lf", &data);
input[i] = data;
}

fclose(stream_input);

status = DftiCreateDescriptor( &my_desc1_handle, DFTI_SINGLE,DFTI_REAL, 2, l);
status = DftiCommitDescriptor( my_desc1_handle);
status = DftiComputeForward( my_desc1_handle, input);
status = DftiFreeDescriptor(&my_desc1_handle);

free(input);
return 0;
}

2 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

Hello,

DftiComputeForward/DftiComputeBackward functions fetch input and store output data according to stride parameters, DFTI_INPUT_STRIDES and DFTI_OUTPUT_STRIDES. If not set explicitly before the call to DftiCommitDescriptor, the strides will have default values which are {0, N, 1} for two-dimensional MxN transform in C/C++ (row-major layout) and {0,1,M} in Fortran (column-major layout).

In the example provided, real-valued input element X(r,c) will be fetched from location X[r*N + c*1] which is likely what was expected. For output, however, element Y(r,c) should be stored into extended matrix of size (M+2)x(N+2), so the strides should be set explicitly to {0,N+2,1} through the call of DftiSetValue(hand, DFTI_OUTPUT_STRIDES,...).

Once the output strides are set as described, the out-of-place transform will work correctly. However, the default is in-place computation, so the out-of-place should be set explicitly through the call of DftiSetValue(hand,DFTI_PLACEMENT, DFTI_NOT_INPLACE).

If, however, in-place transform is needed, be aware that the in-place transform disregards the output strides, so one has to set the input strides {0,N+2,1}, through the call of DftiSetValue(hand,DFTI_INPUT_STRIDES,...). With this, however, the data should also be placed with the stride, so the read loop should be modified appropriately.

This is one part of the story, the other part is how to interpret the result which is stored by default in 2D CCS packed format. The format is rather tricky to comprehend. Instead of explaining it, may I suggest using CCE format instead. This can be achieved through the call of DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX) ?

Hope this helps.

Thanks

Dima

Leave a Comment

Please sign in to add a comment. Not a member? Join today