#1




rlft3.f90 inverse transform
Hi all
I am trying to implement rlft3 subroutine to perform a inverse fourier transform. The size of my arrays are in powers of two, and my input complex array satisfies the required symmetry properties. However the output I obtain is not strictly a real array. I tried to fourier transform the output, and it gives the input array. So I thought may there is a mistake with how I set up my speq array. do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) speq(j,k)=data0(0,jc,kc) enddo enddo do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) enddo enddo function indc(i,N) ! implicit none integer:: indc,i,N ! if(i<=N/2) then indc=i1 else indc=i1N endif ! return end function indc I ensure symmetry in the code too. Is there is a mistake with what I am doing ? kindly help me in this regard Thanks John 
#2




Quote:
Anyhow, where did the "something" come from? That is, where and how did you get the data that comprises the 3D complex collection of data values that supposedly represents the FFT of a 3D real function? (See Footnote.) Maybe you can show us how you fill the data array that you are using with rlft3 from the complex data values. ... Quote:
Quote:
Quote:
Regards, Dave Footnote: In order to understand what the heck is going on with the data array and spec and speq, some people find it instructive to start with a real data array and perform the forward FFT. (Use a nontrivial something for which is easy to predict the results.) Look at the data array and look at speq after the call to rlft2 or rlft3 and make sure you can see exactly what happened. 
#3




Re: rlft3
Hi Dave
Thanks for your detailed reply. I ve attached my code for you to look at here is my code program fftcheck ! implicit none ! integer, parameter :: Nv=8 complex,dimension(Nv,Nv) :: speq complex,dimension(Nv/2:Nv/2,Nv/2:Nv/2,Nv/2:Nv/2) :: data0 complex,dimension(Nv/2,Nv,Nv) :: data double precision :: ran1,f1,f2 integer :: idum,i,j,k,indc,jc,kc ! call system_clock(idum) ! open(unit=220,file='input_speq.dat',status='unknow n') open(unit=221,file='outpu_real.dat',status='unknow n') open(unit=222,file='real_inver.dat',status='unknow n') open(unit=223,file='freq_nyqst.dat',status='unknow n') ! do i=Nv/2,Nv/2 do j=Nv/2,Nv/2 do k=Nv/2,Nv/2 f1=ran1(idum) f2=ran1(idum) data0(i,j,k)=cmplx(f1,f2) enddo enddo enddo ! ! Symmetry due to the realness of data in realspace ! do i=Nv/2,Nv/2 do j=Nv/2,Nv/2 do k=Nv/2,Nv/2 data0(i,j,k)=conjg(data0(i,j,k)) enddo enddo enddo ! ! Store the workspace arrays for fft ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) data(i,j,k)=data0(i1,jc,kc) enddo enddo enddo ! do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) speq(j,k)=data0(Nv/2,jc,kc) enddo enddo ! ! We try to retain symmetry here ! data(1,1,1)=cmplx(real(data(1,1,1)),0.0) ! global average field data(1,Nv/2+1,1)=cmplx(real(data(1,Nv/2+1,1)),0.0) data(1,1,Nv/2+1)=cmplx(real(data(1,1,Nv/2+1)),0.0) data(1,Nv/2+1,Nv/2+1)=cmplx(real(data(1,Nv/2+1,Nv/2+1)),0.0) speq(1,1)=cmplx(real(speq(1,1)),0.0) speq(Nv/2+1,1)=cmplx(real(speq(Nv/2+1,1)),0.0) speq(1,Nv/2+1)=cmplx(real(speq(1,Nv/2+1)),0.0) speq(Nv/2+1,Nv/2+1)=cmplx(real(speq(Nv/2+1,Nv/2+1)),0.0) ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(220,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc enddo enddo ! ! Inverse FFT ing the data ! call rlft3(data,speq,Nv,Nv,Nv,1) ! do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc enddo enddo ! ! FFT ! ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(221,'(2f16.8,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! call rlft3(data,speq,Nv,Nv,Nv,1) ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(222,'(2f16.8,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! stop ! end program fftcheck ! !************************************ !************************************ ! function indc(i,N) ! implicit none integer:: indc,i,N ! if(i<=N/2) then indc=i1 else indc=i1N endif ! return end function indc ! !************************************ !************************************ Thanks a lot for your help Cheers John 
#4




And to check if my rlft3.f90 works right, I took real data, fourier transformed it, Then inverse fourier transformed it again to check if things were fine, and it worked right
I ve attached the code here program fftcheck ! implicit none ! integer, parameter :: Nv=16 complex,dimension(Nv,Nv) :: speq complex,dimension(Nv/2:Nv/2,Nv/2:Nv/2,Nv/2:Nv/2) :: data0 complex,dimension(Nv/2,Nv,Nv) :: data double precision :: ran1,f1,f2 integer :: idum,i,j,k,indc,jc,kc,finder ! call system_clock(idum) ! write(*,*) 'Randomn number is ',idum ! open(unit=220,file='input_real.dat',status='unknow n') open(unit=221,file='outpu_speq.dat',status='unknow n') open(unit=222,file='speq_inver.dat',status='unknow n') open(unit=223,file='freq_nyqst.dat',status='unknow n') ! do i=Nv/2,Nv/2 do j=Nv/2,Nv/2 do k=Nv/2,Nv/2 f1=ran1(idum) f2=0.0d0 data0(i,j,k)=cmplx(f1,f2) enddo enddo enddo ! ! ! Store the workspace arrays for fft ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) data(i,j,k)=data0(i1,jc,kc) enddo enddo enddo ! ! ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(220,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! call rlft3(data,speq,Nv,Nv,Nv,1) ! do jc=1,Nv j=finder(jc,Nv) do kc=1,Nv k=finder(kc,Nv) if (data(1,jc,kc).ne.conjg(data(1,j,k))) write(*,*) 'Hello',jc,j,kc,k enddo enddo ! do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc enddo enddo ! ! Finally do the FFT ! ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(221,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! call rlft3(data,speq,Nv,Nv,Nv,1) ! do i=1,Nv/2 do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(222,'(2f16.5,3I8)') real(data(i,j,k)),aimag(data(i,j,k)),i1,jc,kc enddo enddo enddo ! do j=1,Nv jc=indc(j,Nv) do k=1,Nv kc=indc(k,Nv) write(223,'(2f16.5,2I8)') real(speq(j,k)),aimag(speq(j,k)),jc,kc enddo enddo ! stop ! end program fftcheck ! !************************************ !************************************ function indc(i,N) ! implicit none integer:: indc,i,N ! if(i<=N/2) then indc=i1 else indc=i1N endif ! return end function indc ! !************************************ !************************************ ! Thank you once again John 
#5




Here's what I was trying to suggest for starters: Do the forward transform by calling filling the 8x8x8 real data matrix and calling rlft3(,,,+1). Inspect the results (the transformed 4x8x8 complex data matrix and the 8x8 complex speq matrix).
By my reckoning, here is a list of members of the transformed data matrix and the speq matrix that satisfy the symmetry requirements. (That is, I show index values of pairs that are complex conjugates.) If the members of the matrices do not have the conjugate relationship indicated in the table below, the matrices might not represent the results of a 3d Fourier transformation of a matrix of real numbers (so trying rlft3(,,,1) won't give anything really meaningful in the context of what I perceive you are asking.) Code:
Number of complex data matrix values = 256 Here is a list of data matrix conjugates: * (1,1,1) (1,1,1) * (1,1,2) (1,1,8) (1,1,3) (1,1,7) (1,1,4) (1,1,6) * (1,1,5) (1,1,5) * (1,2,1) (1,8,1) (1,2,2) (1,8,8) (1,2,3) (1,8,7) (1,2,4) (1,8,6) (1,2,5) (1,8,5) (1,2,6) (1,8,4) (1,2,7) (1,8,3) (1,2,8) (1,8,2) (1,3,1) (1,7,1) (1,3,2) (1,7,8) (1,3,3) (1,7,7) (1,3,4) (1,7,6) (1,3,5) (1,7,5) (1,3,6) (1,7,4) (1,3,7) (1,7,3) (1,3,8) (1,7,2) (1,4,1) (1,6,1) (1,4,2) (1,6,8) (1,4,3) (1,6,7) (1,4,4) (1,6,6) (1,4,5) (1,6,5) (1,4,6) (1,6,4) (1,4,7) (1,6,3) (1,4,8) (1,6,2) * (1,5,1) (1,5,1) * (1,5,2) (1,5,8) (1,5,3) (1,5,7) (1,5,4) (1,5,6) * (1,5,5) (1,5,5) * Number of complex speq matrix values = 64 Here is a list of speq matrix conjugates: * (1,1) (1,1) * (1,2) (1,8) (1,3) (1,7) (1,4) (1,6) * (1,5) (1,5) * (2,1) (8,1) (2,2) (8,8) (2,3) (8,7) (2,4) (8,6) (2,5) (8,5) (2,6) (8,4) (2,7) (8,3) (2,8) (8,2) (3,1) (7,1) (3,2) (7,8) (3,3) (7,7) (3,4) (7,6) (3,5) (7,5) (3,6) (7,4) (3,7) (7,3) (3,8) (7,2) (4,1) (6,1) (4,2) (6,8) (4,3) (6,7) (4,4) (6,6) (4,5) (6,5) (4,6) (6,4) (4,7) (6,3) (4,8) (6,2) * (5,1) (5,1) * (5,2) (5,8) (5,3) (5,7) (5,4) (5,6) * (5,5) (5,5) * Of course a more formal statement and derivation of the conditions would be appropriate, but I think you might get started by looking at some actual results. Maybe you will test a few and then try some actual derivations. Try forcing those conjugate relationships onto your randomly generated matrices and see if rlft3(,,,1) followed by rlft3(,,,+1) show you what you want to see. Maybe there is more to the story (and maybe there is less). Bottom line: don't take my word for it. In other words: Your Mileage May Vary. Regards, Dave 
Thread Tools  
Display Modes  

