Hi, I'm using the FFTW package and am using the fortan wrappers to compute a 3D transform. This I seem to do successfully, however I can not seem to disentangle the raw output.>From studying the figure on page 8 of the fftw 3.1.2 manualthis is what I think is going on: For the output complex array I do the following shuffle: ------------------------------------------------- subroutine FFTWshuffle_R2C_3D(array) integer :: ngrid1,ngrid2,ngrid3 integer :: i,j,k,ii,jj,kk complex(DPC), dimension(:,:,:) :: array complex(DPC), allocatable, dimension(:,:,:) :: image ngrid1 = size(array,dim=1) ngrid2 = size(array,dim=2) ngrid3 = size(array,dim=3) allocate (image(ngrid1,ngrid2,ngrid3)) do j=1,ngrid2 do i=1,ngrid1 ii = i + ngrid1/2 jj = j + ngrid2/2 if (ii.gt.ngrid1) ii=ii-ngrid1 if (jj.gt.ngrid2) jj=jj-ngrid2 image(i,j,:)=array(ii,jj,:) enddo enddo array=image deallocate (image) end subroutine FFTWshuffle_R2C_3D ------------------------------------------ (I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1 and already contains just the +ve frequncies.) ---------This does not give me the correct answer, so any help as to what is going wrong would be welcome. Thanks Robert

# FFTW Fortran Real2Complex R2C output array structure? fftw_plan_dft_r2c_3d?

Started by ●June 12, 2007

Reply by ●June 13, 20072007-06-13

On Jun 12, 1:16 pm, "Robert S." <r...@astro.upenn.edu> wrote:> I'm using the FFTW package and am using the fortan wrappers > to compute a 3D transform. This I seem to do successfully, however > I can not seem to disentangle the raw output.>From studying the figure on page 8 of the fftw 3.1.2 manual > > [...] > > (I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1 > and already > contains just the +ve frequncies.)Remember that a Fortran array corresponds to a C array with *reversed* dimensions (column-major vs. row-major format). So, in the C FFTW interface, for an r2c (real-input) transform the last dimension is cut in ~half, but in the Fortran interface the *first* dimension is cut in ~half. This is discussed in the FFTW manual. See the third bullet point of: http://www.fftw.org/doc/Fortran_002dinterface-routines.html and the last example of: http://www.fftw.org/doc/Fortran-Examples.html Regards, Steven G. Johnson

Reply by ●June 13, 20072007-06-13

On Jun 13, 5:40 am, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:> On Jun 12, 1:16 pm, "Robert S." <r...@astro.upenn.edu> wrote: > > > I'm using the FFTW package and am using the fortan wrappers > > to compute a 3D transform. This I seem to do successfully, however > > I can not seem to disentangle the raw output.>From studying the figure on page 8 of the fftw 3.1.2 manual > > > [...] > > > (I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1 > > and already > > contains just the +ve frequncies.) > > Remember that a Fortran array corresponds to a C array with *reversed* > dimensions (column-major vs. row-major format). So, in the C FFTW > interface, for an r2c (real-input) transform the last dimension is cut > in ~half, but in the Fortran interface the *first* dimension is cut in > ~half. > > This is discussed in the FFTW manual. See the third bullet point of: > http://www.fftw.org/doc/Fortran_002dinterface-routines.html > and the last example of: > http://www.fftw.org/doc/Fortran-Examples.html > > Regards, > Steven G. JohnsonDear Steven, Many thanks! The routine works pefectly. Best, Robert AND just for completeness: The working frequency shuffle routine, based on Stephen's comments, is presented below. This places the zero frequency value at the position array(1, ng2/2 , ng3/2 ) ----------------------------------------------------------------------------------------------- subroutine FFTWshuffle_R2C_3D(array) integer :: ng1,ng2,ng3 integer :: j,k,jj,kk complex(DPC), dimension(:,:,:) :: array complex(DPC), allocatable, dimension(:,:,:) :: image ng1 = size(array,dim=1) ng2 = size(array,dim=2) ng3 = size(array,dim=3) allocate (image(ng1,ng2,ng3)) do k=1,ng3 do j=1,ng2 jj = j + ng2/2 + 1 kk = k + ng3/2 + 1 if (jj.gt.ng2) jj=jj-ng2 if (kk.gt.ng3) kk=kk-ng3 image(:,j,k)=array(:,jj,kk) enddo enddo array=image deallocate (image) end subroutine FFTWshuffle_R2C_3D -------------------------------------------------------------------