搜我FFT_code_in_Fortran

2020年11月25日发(作者：何光)
快速傅立叶变换(FFT)的FORTRAN程序代码公开！

虽然现在网上的FF T程序很多，但是有的结果是错误的，比如说，经过一次正反变换后，
不能得到原来的数列，所以，本人 测试了多个程序，提供几个能用的，给大家，希望大家支
持！

下面这个短小精悍，速度快，大牛牛写的！

SUBROUTINE FOUR1(DATA,NN,ISIGN)
! ISIGN: -1:反变换 1： 正变换

REAL*8 WR, WI, WPR, WPI, WTEMP, THETA
DIMENSION DATA(2*NN)

N = 2*NN
J = 1
DO 11 I = 1, N, 2
IF(.I) THEN
TEMPR = DATA(J)
TEMPI = DATA(J+1)
DATA(J) = DATA(I)
DATA(J+1) = DATA(I+1)
DATA(I) = TEMPR
DATA(I+1) = TEMPI
END IF
M = N / 2
1 IF((.2).AND.(.M)) THEN
J = J - M
M = M / 2
GO TO 1
END IF
J = J + M
11 CONTINUE
MMAX = 2
2 IF() THEN
ISTEP = 2 * MMAX
THETA = 6.283D0 / (ISIGN*MMAX)
WPR = -2.D0 * DSIN(0.5D0*THETA)**2
WPI = DSIN(THETA)
WR = 1.D0
WI = 0.D0
DO 13 M = 1, MMAX, 2
DO 12 I = M, N, ISTEP
J = I + MMAX
TEMPR = SNGL(WR) * DATA(J) - SNGL(WI) * DATA(J+1)
TEMPI = SNGL(WR) * DATA(J+1) + SNGL(WI) * DATA(J)
DATA(J) = DATA(I) - TEMPR
DATA(J+1) = DATA(I+1) - TEMPI
DATA(I) = DATA(I) + TEMPR
DATA(I+1) = DATA(I+1) + TEMPI
12 CONTINUE
WTEMP = WR
WR = WR * WPR - WI * WPI + WR
WI = WI * WPR + WTEMP * WPI + WI
13 CONTINUE
MMAX = ISTEP
GO TO 2
END IF

RETURN
END

这个程序也很不错！

c----------------- --------------------------------------------c
c c
c Subroutine sffteu( x, y, n, m, itype ) c
c c
c This routine is a slight modification of a complex split c
c radix FFT routine presented by C.S. Burrus. The original c
c program header is shown below. c
c c
c Arguments: c
c x - real array containing real parts of transform c
c sequence (in/out) c
c y - real array containing imag parts of transform c
c sequence (in/out) c
c n - integer length of transform (in) c
c m - integer such that n = 2**m (in) c
c itype - integer job specifier (in) c
c itype .ne. -1 --> foward transform c
c itype .eq. -1 --> backward transform c
c c
c The forward transform computes c
c X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N) c
c c
c The backward transform computes c
c x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N) c