傅里叶变换--fortran

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(J.GT.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((M.GE.2).AND.(J.GT.M)) THEN

J = J - M

M = M / 2

GO TO 1

END IF

J = J + M

11 CONTINUE

MMAX = 2

2 IF(N.GT.MMAX) THEN

ISTEP = 2 * MMAX

THETA = 6.28318530717959D0 / (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) + SN GL(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 spli t c

c radix FFT routine presente

d by C.S. Burrus. Th

e origin al 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 --> fowar

d transform

c

c itype .eq. -1 --> backwar

d transfor m c

c

c

c The forwar

d transform computes

c

c X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N)

c

c

c

c The backwar

d transform computes

c

c x(j) = (1/N) * sum_{k=0}^{N-1} X(k)*exp(2ijk*pi/N)

c

c

c

c

c

c Requires standar

d FORTRAN functions - sin, cos

c

c

c

c Steve Kifowit, 9 July 1997

c

c

c

C-------------------------------------------------------------C

C A Duhamel-Hollman Split-Radix DIF FFT

C

C Reference: Electronics Letters, January 5, 1984

C

C Complex input and output in data arrays X and Y

C

C Length is N = 2**M

C

C

C

C C.S. Burrus Rice University

Dec 1984 C

C-------------------------------------------------------------C

c

SUBROUTINE SFFTEU( X, Y, N, M, ITYPE )

INTEGER N, M, ITYPE

REAL X(*), Y(*)

INTEGER I, J, K, N1, N2, N4, IS, ID, I0, I1, I2, I3

REAL TWOPI, E, A, A3, CC1, SS1, CC3, SS3

REAL R1, R2, S1, S2, S3, XT

INTRINSIC SIN, COS

PARAMETER ( TWOPI = 6.2831853071795864769 )

c

IF ( N .EQ. 1 ) RETURN

c

IF ( ITYPE .EQ. -1 ) THEN

DO 1, I = 1, N

Y(I) = - Y(I)

1 CONTINUE

ENDIF

c

N2 = 2 * N

DO 10, K = 1, M-1

N2 = N2 / 2

N4 = N2 / 4

E = TWOPI / N2

A = 0.0

DO 20, J = 1, N4

A3 = 3 * A

CC1 = COS( A )

SS1 = SIN( A )

CC3 = COS( A3 )