FFTW 与 fourn 子程序结果不匹配
FFTW vs fourn subroutine result does not match
我 converting/rewriting 从旧的 Fortran 代码库到现代代码库。代码库的一部分使用 fourn 子程序（来自 Numerical receipies 书）用于 FFT 目的。但是当我试图用 FFTW 库做同样的事情时，它不会产生相同的结果。我在这里很困惑。您可以在此处找到输入数据的代码：https://github.com/Koushikphy/fft_test/tree/master/notworking
使用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
使用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
输入文件fort.75可以在这里找到https://github.com/Koushikphy/fft_test/blob/master/notworking/fort.75
对于，测试我还用不同的输入做了一个测试，我对 sin 数据做了 FFT，结果完全匹配 (https://github.com/Koushikphy/fft_test/tree/master/working)。
fftw 方法
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
和fourn方法
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
谁能告诉我我做错了什么？
看起来问题是由于FFTW和Numerical Recipes中离散FFT的定义不同所致。具体来说，根据manual page，FFTW中的“前向”FFT定义为
（对应于 fftw3.f 中定义的 FFTW_FORWARD = -1）。另一方面，根据第 12.4 节：“二维或多维 FFT”和等式（12.4. 1）（在“Fortran 中的 NR”一书中）。 fourn() 的 header part 表示：
Replaces data by its ndim-dimensional discrete Fourier transform,
if isign is input as 1. (...snip...)
If isign is input as −1, data is replaced by its
inverse transform times the product of the lengths of all dimensions.
所以看起来 FFTW 中的“前向”变换对应于 NR 中的“后向”变换。由于 FFTW 和 fourn() 都没有在任何步骤中对结果进行归一化，我认为我们可以将 isign 从 1 更改为 -1 来比较结果：
...
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
...
然后，两个代码对输入文件给出相同的结果fort.75（Re-Im 曲线如下所示，在 NR 和 FFTW 之间匹配）。
对于使用 sin 数据的第二个代码，“正向”变换的结果在 FFTW 和 NR 之间是不同的（即彼此复共轭），而如果我们翻转它们会变得相同 isign在 fourn()（如预期）。
我 converting/rewriting 从旧的 Fortran 代码库到现代代码库。代码库的一部分使用 fourn 子程序（来自 Numerical receipies 书）用于 FFT 目的。但是当我试图用 FFTW 库做同样的事情时，它不会产生相同的结果。我在这里很困惑。您可以在此处找到输入数据的代码：https://github.com/Koushikphy/fft_test/tree/master/notworking
使用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
使用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
输入文件fort.75可以在这里找到https://github.com/Koushikphy/fft_test/blob/master/notworking/fort.75
对于，测试我还用不同的输入做了一个测试，我对 sin 数据做了 FFT，结果完全匹配 (https://github.com/Koushikphy/fft_test/tree/master/working)。
fftw 方法
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
和fourn方法
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
谁能告诉我我做错了什么？
看起来问题是由于FFTW和Numerical Recipes中离散FFT的定义不同所致。具体来说，根据manual page，FFTW中的“前向”FFT定义为
（对应于 fftw3.f 中定义的 FFTW_FORWARD = -1）。另一方面，根据第 12.4 节：“二维或多维 FFT”和等式（12.4. 1）（在“Fortran 中的 NR”一书中）。 fourn() 的 header part 表示：
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.
所以看起来 FFTW 中的“前向”变换对应于 NR 中的“后向”变换。由于 FFTW 和 fourn() 都没有在任何步骤中对结果进行归一化，我认为我们可以将 isign 从 1 更改为 -1 来比较结果：
...
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
...
然后，两个代码对输入文件给出相同的结果fort.75（Re-Im 曲线如下所示，在 NR 和 FFTW 之间匹配）。
对于使用 sin 数据的第二个代码，“正向”变换的结果在 FFTW 和 NR 之间是不同的（即彼此复共轭），而如果我们翻转它们会变得相同 isign在 fourn()（如预期）。