- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

As I understand, for a real-to-complex FFT, I get the output as an array of 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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

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

I hope this will help you.

Thanks

Dima

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Todd

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

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

I hope this will help you.

Thanks

Dima

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Thanks a bunch,

Jozsef

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page