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) + 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
c                                                             c
c                                                             c
c  Requires standard 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 )
     SS3 = SIN( A3 )
     A = J * E
     IS = J
     ID = 2 * N2
 40         DO 30, I0 = IS, N-1, ID
        I1 = I0 + N4
        I2 = I1 + N4
        I3 = I2 + N4
        R1 = X(I0) - X(I2)
        X(I0) = X(I0) + X(I2)
        R2 = X(I1) - X(I3)
        X(I1) = X(I1) + X(I3)
        S1 = Y(I0) - Y(I2)
        Y(I0) = Y(I0) + Y(I2)
        S2 = Y(I1) - Y(I3)
        Y(I1) = Y(I1) + Y(I3)
        S3 = R1 - S2
        R1 = R1 + S2
        S2 = R2 - S1
        R2 = R2 + S1
        X(I2) = R1 * CC1 - S2 * SS1
        Y(I2) = - S2 * CC1 - R1 * SS1
        X(I3) = S3 * CC3 + R2 * SS3
        Y(I3) = R2 * CC3 - S3 * SS3
 30         CONTINUE
     IS = 2 * ID - N2 + J
     ID = 4 * ID
     IF ( IS .LT. N ) GOTO 40
 20      CONTINUE
 10   CONTINUE
c
C--------LAST STAGE, LENGTH-2 BUTTERFLY ----------------------C
c
      IS = 1
      ID = 4
 50   DO 60, I0 = IS, N, ID
  I1 = I0 + 1
  R1 = X(I0)
  X(I0) = R1 + X(I1)
  X(I1) = R1 - X(I1)
  R1 = Y(I0)
  Y(I0) = R1 + Y(I1)
  Y(I1) = R1 - Y(I1)
 60   CONTINUE
      IS = 2 * ID - 1
      ID = 4 * ID
      IF ( IS .LT. N ) GOTO 50
c
C-------BIT REVERSE COUNTER-----------------------------------C
c
 100  J = 1
      N1 = N - 1
      DO 104, I = 1, N1
  IF ( I .GE. J ) GOTO 101
  XT = X(J)
  X(J) = X(I)
  X(I) = XT
  XT = Y(J)
  Y(J) = Y(I)
  Y(I) = XT
 101     K = N / 2
 102     IF ( K .GE. J ) GOTO 103
  J = J - K
  K = K / 2
  GOTO 102
 103     J = J + K
 104  CONTINUE
c
      IF ( ITYPE .EQ. -1 ) THEN
  DO 2, I = 1, N
     X(I) = X(I) / N
     Y(I) = - Y(I) / N
 2       CONTINUE
      ENDIF
c
      RETURN
c
c ... End of subroutine SFFTEU ...
c
      END
!------------------------------------------------------------------

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。