FFTW vs fourn subroutine result does not match
I'm converting/rewriting a old Fortran codebase to modern one. One of the segment of the code base uses fourn subroutine (from the Numerical receipies book) for FFT purpose. But when I'm trying to do that exact thing with the FFTW library it does not produce same result. I'm confused here. You can find the code the input data here : https://github.com/Koushikphy/fft_test/tree/master/notworking
The code that uses fourn:
program test
implicit none
integer, parameter :: n=65536
complex(kind=8) ::inp(n) = 0.0d0
real(kind=8) :: sn, urt(2*n)
integer :: i, ii
sn = 1.0d0/sqrt(real(n,kind=8))
do i=1,9070
read(75,'(i4, 2f20.16)') ii, inp(i)
enddo
do i=1,n
urt(2*i-1)= real(inp(i))
urt(2*i) = aimag(inp(i))
enddo
! forward
call fourn(urt,[n],1,1)
do i=1,n
write(201,'(i4, 2f20.16)')i, urt(2*i-1), urt(2*i)
enddo
end program test
SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN)
INTEGER ISIGN,NDIM,NN(NDIM)
! C REAL DATA(*)
DOUBLE PRECISION DATA(*)
INTEGER I1,I2,I2REV,I3,I3REV,IBIT,IDIM,IFP1,IFP2,IP1,IP2,IP3,K1,K2,N,NPREV,NREM,NTOT
! C REAL TEMPI,TEMPR
DOUBLE PRECISION TEMPI,TEMPR
DOUBLE PRECISION THETA,WI,WPI,WPR,WR,WTEMP
NTOT=1
DO IDIM=1,NDIM
NTOT=NTOT*NN(IDIM)
ENDDO
NPREV=1
DO IDIM=1,NDIM
N=NN(IDIM)
NREM=NTOT/(N*NPREV)
IP1=2*NPREV
IP2=IP1*N
IP3=IP2*NREM
I2REV=1
DO I2=1,IP2,IP1
IF (I2.LT.I2REV) THEN
DO I1=I2,I2+IP1-2,2
DO I3=I1,IP3,IP2
I3REV=I2REV+I3-I2
TEMPR=DATA(I3)
TEMPI=DATA(I3+1)
DATA(I3)=DATA(I3REV)
DATA(I3+1)=DATA(I3REV+1)
DATA(I3REV)=TEMPR
DATA(I3REV+1)=TEMPI
ENDDO
ENDDO
ENDIF
IBIT=IP2/2
1 IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN
I2REV=I2REV-IBIT
IBIT=IBIT/2
GOTO 1
ENDIF
I2REV=I2REV+IBIT
ENDDO
IFP1=IP1
2 IF (IFP1.LT.IP2) THEN
IFP2=2*IFP1
THETA=ISIGN*6.28318530717959D0/(IFP2/IP1)
WPR=-2.0D0*SIN(0.5D0*THETA)**2
WPI=SIN(THETA)
WR=1.0D0
WI=0.0D0
DO I3=1,IFP1,IP1
DO I1=I3,I3+IP1-2,2
DO I2=I1,IP3,IFP2
K1=I2
K2=K1+IFP1
TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1)
TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2)
DATA(K2)=DATA(K1)-TEMPR
DATA(K2+1)=DATA(K1+1)-TEMPI
DATA(K1)=DATA(K1)+TEMPR
DATA(K1+1)=DATA(K1+1)+TEMPI
ENDDO
ENDDO
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
ENDDO
IFP1=IFP2
GOTO 2
ENDIF
NPREV=N*NPREV
ENDDO
RETURN
END
The code that uses fftw:
program test
implicit none
integer, parameter :: n=65536
complex(kind=8) :: inp(n)=0.0d0
integer(kind=8) :: plan
real(kind=8) :: sn
integer :: i, ii
sn = 1.0d0/sqrt(real(n,kind=8))
call dfftw_plan_dft_1d(plan,n,inp,inp,-1,0) !forward plan
do i=1,9070
read(75,'(i4, 2f20.16)') ii, inp(i)
enddo
! forward transform
call dfftw_execute_dft(plan, inp, inp)
do i =1,n
write(101,'(i4, 2f20.16)') i, inp(i)
enddo
end program test
And the input file fort.75 can be found here https://github.com/Koushikphy/fft_test/blob/master/notworking/fort.75
For, test I have also done a test with a different input where I've done FFT for a sin data, where the result match perfectly (https://github.com/Koushikphy/fft_test/tree/master/working).
The fftw approach
program test
implicit none
integer, parameter :: n=65536
real, parameter :: pi = 4.0*atan(1.0)
complex(kind=8), dimension(n) :: x,y,grid,sin2y,out
integer(kind=8) :: pForward, pBackward
real(kind=8) :: sn
integer :: i
sn = 1.0d0/sqrt(real(n,kind=8))
call dfftw_plan_dft_1d(pForward,n,x,y,-1,0) !forward plan
call dfftw_plan_dft_1d(pBackward,n,x,y,+1,0)! backward plan
grid = [(i*2*pi/n, i=1,n)]
sin2y = sin(2*grid)
!actual data
write(100,'(2f20.16)')sin2y
! forward transform
call dfftw_execute_dft(pForward, sin2y, out)
out = out*sn
write(101,'(2f20.16)') out
! backward transform
call dfftw_execute_dft(pBackward, out, sin2y)
sin2y = sin2y*sn
write(102,'(2f20.16)') sin2y
end program test
and the fourn approach
program test
implicit none
integer, parameter :: n=8192
real, parameter :: pi = 4.0*atan(1.0)
complex(kind=8), dimension(n) ::grid,sin2y
real(kind=8) :: sn, urt(2*n)
integer :: i, nn(1)
sn = 1.0d0/sqrt(real(n,kind=8))
grid = [(i*2*pi/n, i=1,n)]
sin2y = sin(2*grid)
!actual data
write(200,'(2f20.16)')sin2y
do i=1,n
urt(2*i-1)= real(sin2y(i))
urt(2*i) = aimag(sin2y(i))
enddo
nn = n
! forward
call fourn(urt,nn,1,1)
urt = urt*sn
do i=1,n
write(201,'(2f20.16)')urt(2*i-1:2*i)
enddo
!backward
call fourn(urt,nn,1,-1)
urt = urt*sn
do i=1,n
write(202,'(2f20.16)')urt(2*i-1:2*i)
enddo
end program test
Can anyone tell me what am I doing wrong here?
It looks like the problem is due to the different definition of discrete FFT in FFTW and Numerical Recipes. Specifically, according to the manual page, the "forward" FFT in FFTW are defined as
(which corresponds to FFTW_FORWARD = -1 as defined in fftw3.f). On the other hand, the "forward" FFT in Numerical Recipes (NR) seems to be defined with exp(+i ...), according to Sec.12.4: "FFT in Two or More Dimensions" and eq.(12.4.1) (in the "NR in Fortran" book). The header part of fourn() says:
Replaces data by its ndim-dimensional discrete Fourier transform, if
isignis input as 1. (...snip...) Ifisignis input as −1, data is replaced by its inverse transform times the product of the lengths of all dimensions.
so it seems that the "forward" transform in FFTW corresponds to the "backward" transform in NR. Since both FFTW and fourn() do not normalize the results in any step, I think we can just change the isign from 1 to -1 to compare the results:
...
urt(:) = 0 !<--- clear the entire urt(:) by zero...
do i=1,n
urt(2*i-1) = real(inp(i))
urt(2*i) = aimag(inp(i))
enddo
! forward
!call fourn(urt, [n], 1, 1) !<--- uses exp(+i ...) for "forward" transform in NR
call fourn(urt, [n], 1, -1) !<--- uses exp(-i ...) for "backward" transform in NR
...
Then, both the codes give the same result for the input file fort.75 (with the Re-Im curves shown below, which match between NR and FFTW).
For the second code using the sin data, the results of the "forward" transform are different between FFTW and NR (i.e. complex-conjugate to each other), while they become identical if we flip isign in fourn() (as expected).