Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6956 Discussions

FFT real to complex transform gives different result than Matlab

Zhu_W_Intel
Employee
669 Views
I wonder whether I did something wrong here.
Here is my example: a simple 2*3 matrix
float* pWeightMatrixfft = SYS_MALLOC(float, 6);

pWeightMatrixfft[0] = 1.0f;
pWeightMatrixfft[1] = 2.0f;
pWeightMatrixfft[2] = 3.0f;
pWeightMatrixfft[3] = 4.0f;
pWeightMatrixfft[4] = 5.0f;
pWeightMatrixfft[5] = 6.0f;

float* pWeightMatrixfftout = SYS_MALLOC(float, 12);
memset(pWeightMatrixfftout, 0, 12*sizeof(float));

DFTI_DESCRIPTOR_HANDLE my_desc1_handled =0;
l[0] = 2;
l[1] = 3;
status = DftiCreateDescriptor( &my_desc1_handled, DFTI_SINGLE, DFTI_REAL, 2, l);
status = DftiSetValue( my_desc1_handled, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
status = DftiCommitDescriptor( my_desc1_handled);
status = DftiComputeForward( my_desc1_handled, pWeightMatrixfft, pWeightMatrixfftout );
status = DftiFreeDescriptor(&my_desc1_handled)

The output is : 21 0 -3 -9 0 0 9 0 0 0 0 0

But if I ran in Matlab:
x = [1 2 3
4 5 6]

x =

1 2 3
4 5 6

>> fft2(x)

ans =

21.0000 -3.0000 + 1.7321i -3.0000 - 1.7321i
-9.0000 0 0


Did I do anything wrong? If I run with DFTI_COMPLEX when the input is a complex array with all imag part to be 0, I can get the same result as in matlab.

Please help.
0 Kudos
7 Replies
John_S_18
Beginner
669 Views
Zhu:
How did you set your output strides?

DftiSetValue(my_desc1_handled,DFTI_OUTPUT_STRIDES,????);

-Ozzer





0 Kudos
Zhu_W_Intel
Employee
669 Views
Quoting - ozzer



I didn't set. How should I set?
0 Kudos
John_S_18
Beginner
669 Views
I am just a user like you, and so I will refer you to the specifics of FFT's in Chapter 11 of the MKL manual. There are many details about the layout of input and output data. Look for the "Strides" section. I recommend that you invest some time and read Chapter 11 very carefully, especially when transforming real-valued data out-of-place and when dealing with the packed output format.

You have a lead now, and you're welcome to follow it. Good luck, Zhu!

-Ozzer


Quoting - Zhu Wang (Intel)

I didn't set. How should I set?

0 Kudos
Zhu_W_Intel
Employee
669 Views
Quoting - ozzer


Thanks for the useful suggestion. I read the manual, but still don't get it.
I can get the imag part now, but the order the pack format doesn't make sense to me.

I would appreciate some examples which can output similiar results as if I use DFTI_COMPLEX.
0 Kudos
Gennady_F_Intel
Moderator
669 Views
Quoting - ozzer



I didn't set. How should I set?

zhu, specially for such cases, MKL package conatins a lot of DFTI examples (C-interfaces). Please look at the dirdftcsource dir.
0 Kudos
Zhu_W_Intel
Employee
669 Views

zhu, specially for such cases, MKL package conatins a lot of DFTI examples (C-interfaces). Please look at the dirdftcsource dir.

I saw those examples. But I don't see anyone compare the Complex FFT vs. Real FFT.
For example, the matlab will return:
x =

1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16

>> fft2(x)

ans =

1.0e+02 *

1.3600 -0.0800 + 0.0800i -0.0800 -0.0800 - 0.0800i
-0.3200 + 0.3200i 0 0 0
-0.3200 0 0 0
-0.3200 - 0.3200i 0 0 0


But if I run Real FFT in MKL,
I got numbers like:
136 -8 8 -8
-0.539062 0.841543 0 0
-32 32 0 0
-0.841543 -0.539062 0 0
-32 0 0 0
0.539062 -0.841543 0 0
-32 -32 0 0

I don't know how the -0.539062 0.841543 come from.

Would appreciate any help.
0 Kudos
Evgueni_P_Intel
Employee
669 Views
Hi Zhu,

By default real-to-complex FFTs store the transformed data in the CCS format (see the default value for DFTI_PACKED_FORMATp.2999 MKLReference Manual).
This format is described in detail for 2D FFTs onp. 3004.
Notice that, for a 2D FFT of MxN real numbers, the transformed data consist of (M+2)x(N+2) real numbers.
Furthermore, these transformed data are stored in a not so straightforward way :-)
We need to allocatemore memory (not 12 but 20 floats in the very first example in this thread.)

By default all FFTs are in-place and only input strides have a default value.
If we ask for a not-in-place FFT, we must set the output strides ourselves.
See example dftc/real_2d_ccs_single_ex2.c for an example that is closest to your use case.
For the first example in this thread, the output strides must be {0, 5, 1}.

Evgueni.
0 Kudos
Reply