2d FFT CCE format in C

[Jozsef]

I wonder how the frequencies are associated to the output matrix elements after a 2d FFT in CCE format in C.

As I understand, for a real-to-complex FFT, I get the output as an array of y[N/2+1] complex numbers. Now, there is a conjugate-even symmetry in both dimensions. Accordingly, in C, only N/2+1 elements are stored in the N direction. My question is: how should the symmetry be understood in the M direction or in other words, how are the positive and negative frequencies are associated to the complex coefficients in the output array y?

An example might make it more clear. In 1D, the loop

int k;
double Y = 0.0;

for ( k = -N/2; k < N/2+N%2; k++ )
Y += y[abs(k)*2+0]*cos(2.0*PI*j/n*k) - SGN(k)*y[abs(k)*2+1]*sin(2.0*PI*j/n*k);

seems to correctly sum up in Y the approximate Fourier series, evaluated at 2.0*PI*j/n ,using the complex coefficients (stored in y in CCS format) after a 1D FFT.

The question is how can one construct a similar loop in 2D?

Thanks,
Jozsef

[Todd_R_Intel]

Jozsef,

Have you seen the 2D example code in the examples/fftc directory?

Todd

[Jozsef]

Hi Todd,

yes, I have looked at them. (The examples/dftc folder.) They do forward and backward transform. I don't see how this is helpful in how the data should be interpreted.

Thanks,
Jozsef

[barragan_villanueva_]

Hi,

In order to help you it needs to specify your function.
Currently, I have just a guess that you want after forward DFT-transform to restore all data manually doing backward transform. Am I right?
Before you implemented such function for 1D CCE. Maybe to simplify your task let us try to implement such function for 2D complex-to-complex transform first. And after that just correct it for 2D CCE.

[Dmitry_B_Intel · Accepted solution]

Hi Jozsef,

Multidimensional real transform N1-by-N2-by-N3-... produces conjugate-even sequence y with the following symmetry property (MKL Reference Manual, subsection Forward domain of transform):

y(n1,n2,n3,...) = conj( y(N1-n1,N2-n2,N3-n3,...) ).

For 2D case M-by-N, we have:

y(m,n) = conj( y(M-m,N-n) )

Now, how the y(m,n) for the M-by-N transform are stored in C/C++?
Let us use notation y(m,n), for whichwe define:

#define y(m,n) y

This definition may be useful when the dimensions of the data are not known at compile time, for example:

complex *y = (complex*)malloc( sizeof(complex) * M * N ); //M,N are not known until run-time
#define y(m,n) y[(m)*N + (n)] // row-major layout

For the full 2d-sequence, (m,n) would run from (0,0) to (M-1,N-1). Due to symmetry property, only one half of the sequence is stored, where (m,n) run from (0,0) to (M-1,N/2). Thus if one needs y(m,k) with k > N/2, he should use y(m,k) = conj( y(m,N-k) ).

A mapping of the conjugate-even sequence to a form where frequences span from negative to positive with zero frequency in the middle employs the fact that, according to definition of DFT,
y(m,k) = y(M*l+m,N*k+n) for any integers k,l. Thus we may chose this mapping to full sequence y:

m = -(M-1)/2...M/2
n = -(N-1)/2...N/2
y_shifted(m,n) =
y(m,n) where 0<=m<=M/2, 0<=n<=N/2
y(m+M,n) where -(M-1)/2<=m<0, 0<=n<=N/2
y(m,n+N) where 0<=m<=M/2, -(N-1)/2<=n<0
y(m+M,n+N) where -(M-1)/2<=m<0, -(N-1)/2<=n<0

This definition doesn't suit our purpose because it addresses items with n>N/2, that is the half that is not stored. One can map the definition of y_shifted to the computed CCE data that spans index range (0,0)...(M-1,N/2) using the conjugate-even symmetry:

m = -(M-1)/2...M/2
n = -(N-1)/2...N/2
y_shifted(m,n) =
y(m,n) where 0<=m<=M/2, 0<=n<=N/2
y(m+M,n) where -(M-1)/2<=m<0, 0<=n<=N/2
conj( y(0,-n) ) where 0=m, -(N-1)/2<=n<0
conj( y(M-m,-n) ) where 0<=M/2, -(N-1)/2<=n<0
conj( y(-m,-n) ) where -(M-1)/2<=m<0, -(N-1)/2<=n<0

I hope this will help you.

Thanks
Dima

[Jozsef]

Hi Dima,

that's exactly what I was interested in. Works perfect.

Thanks a bunch,
Jozsef