Я конвертирую / переписываю старую кодовую базу Fortran на современную. Один из сегментов кодовой базы использует подпрограмму fourn (из книги числовых квитанций) для целей БПФ. Но когда я пытаюсь сделать то же самое с библиотекой 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, где результат полностью совпадает (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
Может ли кто-нибудь сказать мне, что я здесь делаю не так?