! It runs and computes correctly although it is not in any sense ! sophisticated. However, it is free, and useful if you are in a hurry ! and don't want to think things out for yourself. ! Arthur Wouk (wouk@brl-vgr) ! ******************* ! Fast Fourier Transform ! ******************* ! Code converted using TO_F90 by Alan Miller ! Date: 2003-07-09 Time: 17:02:48 SUBROUTINE fft(x, n, k) ! FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X ! (A COMPLEX ARRAY OF DIMENSION N). ! SOURCE: Joel H. Ferziger; Numerical methods for engineering applications. ! 2nd edition, Wiley (1998). ! X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM. ! N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<32769). ! K = 1 FOR FORWARD TRANSFORM. ! K = -1 FOR INVERSE TRANSFORM. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) INTEGER, INTENT(IN) :: n COMPLEX (dp), INTENT(IN OUT) :: x(n) INTEGER, INTENT(IN) :: k INTEGER :: sby2, s REAL (dp) :: ang, re, im COMPLEX (dp) :: xtemp, t, u(16), v, w LOGICAL :: NEW INTEGER :: i, INDEX, itemp, j, jndex, l, l2n, p, q, stage REAL (dp) :: pi2 = 6.2831853071796_dp, gain = 1.0_dp INTEGER :: no = 0, ko = 0 ! TEST FIRST CALL? NEW = (no /= n) IF (NEW) THEN ! IF FIRST CALL COMPUTE LOG2 (N). l2n = 0 no = 1 10 l2n = l2n + 1 no = no + no IF (no < n) GO TO 10 gain = 1. / n ang = pi2 * gain re = COS(ang) im = SIN(ang) END IF ! COMPUTE COMPLEX EXPONENTIALS IF NOT FIRST CALL IF (.NOT.(.NOT.NEW .AND. k*ko >= 1)) THEN u(1) = CMPLX(re, -SIGN(im, DBLE(k)), KIND=dp) DO i = 2, l2n u(i) = u(i-1) * u(i-1) END DO ! The code from netlib cntained the line: ! K0 = K ko = k END IF ! MAIN LOOP sby2 = n DO stage = 1, l2n v = u(stage) w = (1.0_dp, 0.0_dp) s = sby2 sby2 = s / 2 DO l = 1, sby2 DO i = 1, n, s p = i + l - 1 q = p + sby2 t = x(p) + x(q) x(q) = (x(p)-x(q)) * w x(p) = t END DO w = w * v END DO END DO ! REORDER THE ELEMENTS BY BIT REVERSAL DO i = 1, n INDEX = i - 1 jndex = 0 DO j = 1, l2n jndex = jndex + jndex itemp = INDEX / 2 IF (itemp+itemp /= INDEX) jndex = jndex + 1 INDEX = itemp END DO j = jndex + 1 IF (j >= i) THEN xtemp = x(j) x(j) = x(i) x(i) = xtemp END IF END DO ! FORWARD TRANSFORM DONE IF (k > 0) RETURN ! INVERSE TRANSFORM x(1:n) = x(1:n) * gain RETURN END SUBROUTINE fft