การสลายตัวของ 2D FFT (และ inverse fft) ด้วย fftw ใน fortran77

[แก้ไข 1] เพิ่มตัวเลขเพื่อแสดงข้อมูลต้นฉบับและข้อมูลที่ได้รับ

[แก้ไข 2] ฉันพบข้อผิดพลาด ฉันใช้ fftw_measure แทน fftw_estimate ในการเรียก dfftw_plan_many_dft

[แก้ไข 3] แก้ไขการพิมพ์ผิดในรหัส (แทนที่ u ด้วย u2d ใน dfftw_execute_dft_r2c )

ฉันกำลังลองทำ 2D fft ของอาร์เรย์โดยใช้ 1D fft หลายอัน แทนที่จะใช้ฟังก์ชัน 2D fft ที่มีอยู่แล้วในไลบรารี fftw และต่อมา ฉันต้องทำ fft 2D แบบผกผัน เหตุผลในการทำเช่นนั้นคือ (ในอนาคต) อาเรย์ของฉันจะใหญ่เกินกว่าจะโหลดได้ในคราวเดียวเพื่อทำ 2D fft

โค้ดฉบับร่างที่ 1 ของฉันมีลักษณะเช่นนี้ไม่มากก็น้อยในขณะนี้:

  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

ฉันคาดหวังว่า u2d และ u2d2 จะเหมือนกัน แต่ฉันได้รับค่าที่แตกต่างกันค่อนข้างมาก ฉันทำผิดพลาดที่ไหนสักแห่งหรือไม่?

ต้นฉบับและผลลัพธ์แสดงไว้ด้านล่างนี้ รูปร่างดูคล้ายกันแต่ค่าค่อนข้างแตกต่างกัน (เช่น ค่าต่ำสุดและค่าสูงสุด)

ช่องต้นฉบับ

ฟิลด์ที่ได้รับหลังจาก fft และ i-fft