!2009-1-3  快速付利葉變換子程序的編寫
!電子工程與光電技術學院
program fft      
real,parameter::pi=3.1415926535898 	  
!real::DATA1(32)=(/1,0,3,0,0,0,0,0,0,2,1,0,3,2,1,0,1,1,3,0,3,0,1,1,3,0,3,3,1,1,2,2/)*pi/2;

integer,parameter::N1=32;
real data1(N1*4)

integer::ISIGN=1

do i=1,n1*2,2
data1(i)=1.0/((i+1)/2)
data1(i+1)=0
enddo
do i=2*n1+1,4*n1
data1(i)=0
enddo

print*,DATA1
call FOUR1(DATA1,64,ISIGN)

print*,DATA1

contains
	  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 subroutine
end