Decomposition of 2D FFT (and inverse fft) with fftw in fortran77

Question

[edit 1] Added figures to show the original data and the obtained data

[edit 2] I found my mistakes, I used fftw_measure instead of fftw_estimate in the call of dfftw_plan_many_dft

[edit 3] corrected a typo in the code (replace u with u2d in dfftw_execute_dft_r2c )

I'm trying perform the 2D fft of an array using multiple 1D fft instead of using the 2D fft function already existing in the fftw library. And subsequently, I need to perform the inverse 2D fft. The reason for doing that is that (in the future) my array will be too big to be loaded in one go to perform the 2D fft.

The 1st draft of my code looks more or less like this at the moment:

  double precision u2d(nx,ny),u2d2(nx,ny)
  double complex qhat2d(nx/2+1,ny),qhat2d2(nx/2+1,ny)
  integer N(1)
  integer howmany, idist, odist, istride, ostride
  integer inembed(2), onembed(2)
  integer rank

  ! some function to read the data into u2d

  ! perform x-fft
  N(1) = NX
  howmany = NY
  inembed(1) = NX
  inembed(2) = NY
  istride = 1
  idist = NX
  ostride = 1
  odist = (NX/2+1)
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  write(*,*) 'u', u2d(1,1)
  CALL dfftw_plan_many_dft_r2c(PLAN,rank,N(1),howmany,
&                              u2d,inembed,
&                              istride,idist,
&                              qhat2d,onembed,
&                              ostride,odist,FFTW_ESTIMATE) ! 
  CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d) ! x-fft
  CALL dfftw_destroy_plan(PLAN)


  ! perform y-fft
  N(1) = NY
  howmany = (NX/2+1)
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = (NX/2+1)
  idist = 1
  ostride = (NX/2+1)
  odist = 1
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
&                              qhat2d,inembed,
&                              istride,idist,
&                              qhat2d2,onembed,
&                              ostride,odist,FFTW_FORWARD,
&                              FFTW_MEASURE) ! 
  CALL dfftw_execute_dft(PLAN,qhat2d,qhat2d2) ! y-fft
  CALL dfftw_destroy_plan(PLAN)

  ! normally here, perform some filtering operation
  ! but at the moment, I do nothing

  ! perform inv-y-fft
  N(1) = NY
  howmany = (NX/2+1)
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = (NX/2+1)
  idist = 1
  ostride = (NX/2+1)
  odist = 1
  onembed(1) = (NX/2+1)
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany,
 &                             qhat2d2,inembed,
 &                              istride,idist,
 &                              qhat2d,onembed,
 &                              ostride,odist,FFTW_BACKWARD,
 &                              FFTW_MEASURE) ! 
  CALL dfftw_execute_dft(PLAN,qhat2d2,qhat2d) ! inv-y-fft
  CALL dfftw_destroy_plan(PLAN)

  ! perform inv-x-fft
  N(1) = NX ! I'm not too sure about this value here
  howmany = NY
  inembed(1) = (NX/2+1)
  inembed(2) = NY
  istride = 1
  idist = (NX/2+1)
  ostride = 1
  odist = NX
  onembed(1) = NX
  onembed(2) = NY
  rank = 1
  CALL dfftw_plan_many_dft_c2r(PLAN,rank,N(1),howmany,
&                              qhat2d,inembed,
&                              istride,idist,
&                              u2d2,onembed,
&                              ostride,odist,FFTW_ESTIMATE) ! 
  CALL dfftw_execute_dft_c2r(PLAN,qhat2d,u2d2) ! x-fft
  CALL dfftw_destroy_plan(PLAN)
  write(*,*) 'u2d2', u2d2(1,1)

  do i=1,nx
   do j=1,ny
    u2d2(i,j) = u2d2(i,j) / (nx*ny)
   enddo
  enddo
  write(*,*) 'u2d2', u2d2(1,1) ! here the values u2d2(1,1) is different from u2d(1,1)

  ! some action to write u2d2 to file
  end

I was expecting u2d and u2d2 to be the same but I obtain relatively different values. Did I make a mistake somewhere?

The original and results are shown hereunder. The shape looks similar but the values are relatively different (the minimum and maximum for example).

Original field

Field obtained after fft and i-fft

FFTW does [not normalize the values](http://www.fftw.org/faq/section3.html#whyscaled). Therefore, there will be a factor of the array length if you transform forward and backwards on the same data. — Alexander Vogt, May 23 '16 at 18:27
See, e.g., https://stackoverflow.com/questions/3721125/fftw-inverse-of-forward-fft-not-equal-to-original-function/7871634#7871634 — Alexander Vogt, May 23 '16 at 18:28
or https://stackoverflow.com/questions/4855958/normalising-fft-data-fftw — Alexander Vogt, May 23 '16 at 18:29
or https://stackoverflow.com/questions/30402282/fftw3-inverse-transform-not-work/30402957#30402957 — Alexander Vogt, May 23 '16 at 18:30
In addition to other comments I stress that you must show what is wrong. You must show which values you get and why you believe they are wrong. — Vladimir F Героям слава, May 23 '16 at 21:21
@AlexanderVogt My original data u2d is real so I thought that using r2c was the correct way of doing the first fft. Furthermore at the end I do renormalize my data by dividing the result of the 2 inverse fft by (nx*ny). — Anh Khoa, May 24 '16 at 09:04
@roygvib Actually, a duplicate? http://stackoverflow.com/questions/36885526/when-using-r2c-and-c2r-fftw-in-fortran-are-the-forward-and-backward-dimensions?rq=1 The title is absolutely wrong, however. — Vladimir F Героям слава, May 24 '16 at 17:58
@VladimirF I initially thought so, but if the intermediate Fourier transform is not re-used afterward (even "filtered" before entering the inverse transform), the code seems to work (by assuming that FFTW_ESTIMATE keeps the input array intact). So hmm, I'm not very sure... (one thing sure is that, I would avoid mixing plan creation and actual transform at the same time, though.) — roygvib, May 24 '16 at 18:07

score 1 · Answer 1 · answered May 24 '16 at 11:09

1

I found my mistakes, I used fftw_measure instead of fftw_estimate in the call of dfftw_plan_many_dft.

Correcting that gives me the appropriate original field.

answered May 24 '16 at 11:09

Anh Khoa

23
4

1

This shouldn't matter at all. What could matter is FFTW_PRESERVE_INPUT, but only if you don't save the input somewhere else. – Vladimir F Героям слава May 24 '16 at 11:23
@VladimirF Checking the manual page http://www.fftw.org/doc/Planner-Flags.html, it actually seems that if one uses FFTW_MEASURE, the input array will be overwritten when creating a plan. Very confusing interface... – roygvib May 24 '16 at 16:55
@roygvib Yes, you should not just take the input array as `intent(in)`, it is not such. Not even during the transforms. You should use ` FFTW_PRESERVE_INPUT` then. I strongly suggest to save the input values somewhere else if you need them later for any purpose. It is not save enough to just use `FFTW_ESTIMATE` as the answer tries to suggest. That may hide the problem, but it can appear at any time when the estimate selects an algorithm which does not preserve the input. The input will then get overwritten during a transform instead of during the planning. Therefore this answer is incorrect. – Vladimir F Героям слава May 24 '16 at 17:10
Yeah, and more generally, I suggest separating the creation of plans and actual transforms somehow. In my case, I create a module that handles this internally for the user. But I think the answer itself is correct, because the coding "style" is a matter of taste (to some extent). Nevertheless, I think it should be very useful if a different answer is also added as an alternative (safer) style... – roygvib May 24 '16 at 17:21
@roygvib It is not correct, you can use `FFTW_ESTIMATE` and still get your input overwritten. The only luck is that the planner chose an algorithm that doesn't do the overwriting, but that is pure luck. Next time, with the same code on a different machine, it can fail again. – Vladimir F Героям слава May 24 '16 at 17:43
Although he has `CALL dfftw_execute_dft_r2c(PLAN,u,qhat2d)`, why not `CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d)`? Where is `u` defined at all? I suspect it is a type and my conclusion is still correct. – Vladimir F Героям слава May 24 '16 at 17:46
Ah, "u" should be a typo of u2d. The code worked on my computer with that modification (while changing FFTW_MEASURE to FFTW_ESTIMATE). BTW, is there really some case where FFTW_ESTIMATE can change the input array? Because I want to know more on this, I would really appreciate any additional info... – roygvib May 24 '16 at 17:51
@roygvib not during planning, I thing. But during the transform it could happen. Hard to reproduce, because you don't know what the estimation chosses, but this guy http://stackoverflow.com/questions/29322026/fortran-complex-to-real-fftw-issue?rq=1 did use `FFTW_ESTIMATE`. And actually, your answer here http://stackoverflow.com/questions/36885526/when-using-r2c-and-c2r-fftw-in-fortran-are-the-forward-and-backward-dimensions?rq=1 has encountered the same issue, or not? – Vladimir F Героям слава May 25 '16 at 10:16
[@VladimirF Hi, I will check this more later and get back. (RE the linked Q/A, I remember I got a wrong result if I use c2r thing.)] – roygvib May 25 '16 at 20:23

Vladimir F Героям слава · Accepted Answer · 2016-05-24T20:19:47.033

To clear the confusion. What happens is that the FFTW c2r and r2c routines are not guaranteed to preserve the input arrays. They can overwrite the result with garbage.

Now, you can be lucky and just use FFTW_ESTIMATE instead of FFTW_MEASURE, but it is not a good idea.FFTW_MEASURE tries many algorithms and is therefore likely to try also one which does not preserve the input. FFTW_ESTIMATE will not try to compute anything and will not overwrite the input.

The problem is that your input can be overwritten when performing the transform (executing the plan). You will only be lucky when FFT_ESTIMATE selects an algorithm which preserves the input for you. It is a lottery.

To be sure the input gets preserved you should use FFTW_INPUT_PRESERVE in addition to either FFT_ESTIMATE or FFTW_MEASURE.

You can also not use it and instead save the input somewhere. I think that is often better because FFTW_INPUT_PRESERVE can (quite likely) select a slower algorithm.

Decomposition of 2D FFT (and inverse fft) with fftw in fortran77

2 Answers2