Community
cancel
Showing results for
Did you mean: Beginner
179 Views

## 2d FFT CCE format in C

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
1 Solution Employee
179 Views

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

Thanks
Dima

5 Replies Employee
179 Views
Jozsef,

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

Todd Beginner
179 Views
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 Valued Contributor I
179 Views
Hi,

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. Employee
180 Views

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

Thanks
Dima Beginner
179 Views
Hi Dima,

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

Thanks a bunch,
Jozsef 