- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi all,
I am doing a very similar test to the dfti-examples.
Take example real_3d_cce_double_ex2.c:
(1) It does a forward and a backward fft and it checks the error on the real domain to see if we get back the original data.
Now I'm trying to do the opposite with the same fft settings:
(2) I first do a backward and then a forward fft and I check the error on the complex domain to see if I can get back the original
(complex) data.
However, (1) works (i.e. I get error ~1.0e-15), but in case (2) the error is large.
Details:
The code below uses the following constants:
My settings for (both forward and backward) fft are as follows:
Settings for only forward fft:
Settings for only backward fft:
The code I use to initialize, transform and then test the complex data is as follows:
QUESTION: Are there any particular settings that need to be done as a setup to the FFT routines if one wants to apply backward fft first? What am I missing?
Thanks,
Jozsef
I am doing a very similar test to the dfti-examples.
Take example real_3d_cce_double_ex2.c:
(1) It does a forward and a backward fft and it checks the error on the real domain to see if we get back the original data.
Now I'm trying to do the opposite with the same fft settings:
(2) I first do a backward and then a forward fft and I check the error on the complex domain to see if I can get back the original
(complex) data.
However, (1) works (i.e. I get error ~1.0e-15), but in case (2) the error is large.
Details:
The code below uses the following constants:
[cpp]const MKL_LONG N = 32; const MKL_LONG LENGTHS[3] = {N,N,N}; const MKL_LONG PS[4] = {0,N*N,N,1}; const MKL_LONG SS[4] = {0,N*(N/2+1),N/2+1,1}; const MKL_LONG PO = N*N*N; const MKL_LONG SO = N*N*(N/2+1);
const MKL_LONG n = 1;[/cpp]
My settings for (both forward and backward) fft are as follows:
[cpp]if ( status = DftiCreateDescriptor(fft, DFTI_DOUBLE, DFTI_REAL, 3, LENGTHS) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_PLACEMENT, DFTI_NOT_INPLACE) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_PACKED_FORMAT, DFTI_CCE_FORMAT) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_NUMBER_OF_TRANSFORMS, n) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_FORWARD_SCALE, 1.0) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(*fft, DFTI_BACKWARD_SCALE, 1.0/(double)PO) ) ERR(DftiErrorMessage(status));[/cpp]
Settings for only forward fft:
[cpp] // specifically initialize for forward (physical to spectral) transform if ( status = DftiSetValue(fft, DFTI_INPUT_DISTANCE, PO) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_OUTPUT_DISTANCE, SO) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_INPUT_STRIDES, PS) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_OUTPUT_STRIDES, SS) ) ERR(DftiErrorMessage(status)); if ( status = DftiCommitDescriptor(fft) ) ERR(DftiErrorMessage(status)); // perform forward FFT if ( status = DftiComputeForward(fft,u,U) ) ERR(DftiErrorMessage(status));[/cpp]
Settings for only backward fft:
[cpp] // specifically initialize for backward (spectral to physical) transform if ( status = DftiSetValue(fft, DFTI_INPUT_DISTANCE, SO) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_OUTPUT_DISTANCE, PO) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_INPUT_STRIDES, SS) ) ERR(DftiErrorMessage(status)); if ( status = DftiSetValue(fft, DFTI_OUTPUT_STRIDES, PS) ) ERR(DftiErrorMessage(status)); if ( status = DftiCommitDescriptor(fft) ) ERR(DftiErrorMessage(status)); // perform backward FFT if ( status = DftiComputeBackward(fft,U,u) ) ERR(DftiErrorMessage(status));[/cpp]
The code I use to initialize, transform and then test the complex data is as follows:
[cpp] MKL_LONG ki; double maxerr; double complex *P, *expP; double *p; // allocate memory if ( !(P = (double complex*)malloc(SO*sizeof(double complex))) ) ERR("Can't allocate memory!"); if ( !(exP = (double complex*)malloc(SO*sizeof(double complex))) ) ERR("Can't allocate memory!"); if ( !(p = (double*)malloc(PO*sizeof(double))) ) ERR("Can't allocate memory!"); // initialize complex data for ( ki = 0; ki < SO; ki++ ) P[ki] = some complex random number // inverse & forward fft do backward FFT P -> p do forward FFT p -> exP // test error maxerr = 0.0; for ( ki = 0; ki < SO; ki++ ) { if (fabs(creal(P[ki]-exP[ki])) > maxerr) maxerr = fabs(creal(P[ki]-exP[ki])); if (fabs(cimag(P[ki]-exP[ki])) > maxerr) maxerr = fabs(cimag(P[ki]-exP[ki])); } printf("maxerr = %g\n",maxerr);[/cpp]
QUESTION: Are there any particular settings that need to be done as a setup to the FFT routines if one wants to apply backward fft first? What am I missing?
Thanks,
Jozsef
Link Copied
3 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Jozsef,
Probably, it is input data issue. Complex-to-real transform expects the input is conjugate-even [that isP(n1,n2,n3)==conj(P(N-n1,N-n2,N-n3)) ]. Ifthe input is not conjugate-even,the transform produces junk, because it doesn't attempt to 'symmetrize' the input.
Thanks
Dima
Probably, it is input data issue. Complex-to-real transform expects the input is conjugate-even [that isP(n1,n2,n3)==conj(P(N-n1,N-n2,N-n3)) ]. Ifthe input is not conjugate-even,the transform produces junk, because it doesn't attempt to 'symmetrize' the input.
Thanks
Dima
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
But I thought the conjugate-symmetric part is not stored. Is it?
J
J
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Most of it is not. However, some symmetry should be present (e.g. Im(Y(0,0,0)==0).
Dima
Dima
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page