PROGRAM FFT PARAMETER (NDIMX=512) COMPLEX H1(NDIMX), P1(NDIMX), HOUT(NDIMX) !COMPLEX COEF(NDIMX+15) REAL TMAX, PI, DT, DW INTEGER I, SIGN, ERROR OPEN (UNIT = 20, File = "FFT.txt", STATUS = "UNKNOWN", IOSTAT = ERROR) IF (ERROR /= 0) THEN PRINT *, "***ERROR IN OPENING FILE***", Error STOP END IF OPEN (UNIT = 30, File = "FFT-1.txt", STATUS = "UNKNOWN", IOSTAT = ERROR) IF (ERROR /= 0) THEN PRINT *, "***ERROR IN OPENING FILE***", Error STOP END IF PI = 4. * ATAN(1.) ! Pi TMAX = 10. ! Time interval being sampled DT = TMAX / REAL(NDIMX) ! delta t for time sampling DW = 2. * PI / TMAX ! delta omega for frequency grid ! Sample input function to be Fourier transformed DO I = 1, NDIMX TIME = (I - 1) * DT H1(I) = EXP(- Time * Time / (5.)) * SIN(40. * TIME) P1(I) = H1(I) END DO ! Go to Fourier space SIGN = -1 CALL FFTBYHAND(SIGN, H1, HOUT) ! PRINT THE SPECTRUM DO I = NDIMX / 2 + 1, NDIMX ! here the negative momenta WRITE(*,*) (I - 1. - NDIMX) * DW, H1(I), ABS(HOUT(I)) WRITE(20,*) (I - 1. - NDIMX) * DW, H1(I), ABS(HOUT(I)) END DO DO I = 1, NDIMX / 2 ! here the positive momenta WRITE(*,*) (I - 1.) * DW, H1(I), ABS(HOUT(I)) WRITE(20,*) (I - 1.) * DW, H1(I), ABS(HOUT(I)) WRITE(30,*) (I - 1.) * DW, REAL(H1(I)), REAL(P1(I)) END DO ! Go back to real space SIGN = 1 CALL FFTBYHAND(Sign, HOUT, H1) ! Test the FFT ! compare f(t) with the forward and backward transform DO I = 1, NDIMX IF(ABS((P1(I) - H1(I)) / P1(I)) > 1E-6) THEN WRITE (*,*) "WARNING SOMETHING SCREWED UP" END IF END DO STOP ! ------------------------------------------------------------------------------------------------------- CONTAINS ! ------------------------------------------------------------------------------------------------------- ! Subroutine for slow fourier transform SUBROUTINE FFTBYHAND(SIGN, HIN, HOUT) PARAMETER (NDIMX = 512) COMPLEX E1, SUM, TPII, HIN(NDIMX), HOUT(NDIMX) INTEGER I, I1, SIGN TPII = 2. * (4. * ATAN(1.0)) * (0.,1.) ! 2PI DO I = 1, NDIMX SUM = (0., 0.) DO I1 = 1, NDIMX E1 = EXP(SIGN * TPII * (I1 - I) * (I - 1.) / (NDIMX + 0.)) SUM = SUM + E1 * HIN(I1) END DO HOUT(I) = SUM IF (SIGN == 1) HOUT(I) = SUM / (NDIMX + 0.) ! Normalization END DO RETURN END SUBROUTINE END PROGRAM