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. Johnson
Dear 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
-------------------------------------------------------------------
Reply by Steven G. Johnson●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 Robert S.●June 12, 20072007-06-12
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 manual
this 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