Program 12: Prime Factor FFT Algorithm

Prime Factor FFT Algorithm Below is the Fortran code for a Prime-Factor Algorithm (PFA) FFT allowing factors of the length of 2, 3, 4, 5, and 7. It is followed by an unscrambler.
C---------------------------------------------------
C
C   A PRIME FACTOR FFT PROGRAM WITH GENERAL MODULES
C   COMPLEX INPUT DATA IN ARRAYS  X AND Y
C   COMPLEX OUTPUT IN  A AND B
C   LENGTH  N  WITH  M  FACTORS IN ARRAY  NI
C     N = NI(1)*NI(2)* ... *NI(M)
C   UNSCRAMBLING CONSTANT  UNSC
C     UNSC = N/NI(1) + N/NI(2) +...+ N/NI(M), MOD N
C      C. S. BURRUS, RICE UNIVERSITY, JAN 1987
C
C--------------------------------------------------
C
    SUBROUTINE PFA(X,Y,N,M,NI,A,B,UNSC)
C
    INTEGER  NI(4), I(16), UNSC
        REAL X(1), Y(1), A(1), B(1)
C
    DATA  C31, C32  / -0.86602540,-1.50000000 /
    DATA  C51, C52  /  0.95105652,-1.53884180 /
    DATA  C53, C54  / -0.36327126, 0.55901699 /
    DATA  C55       / -1.25  /
    DATA  C71, C72  / -1.16666667,-0.79015647 /
    DATA  C73, C74  /  0.055854267, 0.7343022 /
    DATA  C75, C76  /  0.44095855,-0.34087293 /
    DATA  C77, C78  /  0.53396936, 0.87484229 /
C
C-----------------NESTED LOOPS----------------------
C
    DO 10 K=1, M
          N1 = NI(K)
      N2 = N/N1
      DO 15 J=1, N, N1
             IT   = J
             DO 30 L=1, N1
                I(L) = IT
            A(L) = X(IT)
            B(L) = Y(IT)
                IT = IT + N2
                IF (IT.GT.N)  IT = IT - N
   30        CONTINUE
             GOTO (20,102,103,104,105,20,107), N1
C
C----------------WFTA N=2--------------------------------
C
  102   R1    = A(1)
    A(1)  = R1 + A(2)
    A(2)  = R1 - A(2)
C
    R1   = B(1)
    B(1) = R1 + B(2)
    B(2) = R1 - B(2)
C
    GOTO 20
C----------------WFTA N=3--------------------------------
C
  103   R2 = (A(2) - A(3)) * C31
    R1 =  A(2) + A(3)
    A(1)= A(1) + R1
    R1  = A(1) + R1 * C32
C
    S2 = (B(2) - B(3)) * C31
    S1 =  B(2) + B(3)
    B(1)= B(1) + S1
    S1  = B(1) + S1 * C32
C
    A(2) = R1 - S2
    A(3) = R1 + S2
    B(2) = S1 + R2
    B(3) = S1 - R2
C
    GOTO 20
C
C----------------WFTA N=4---------------------------------
C
  104   R1 = A(1) + A(3)
    T1 = A(1) - A(3)
    R2 = A(2) + A(4)
    A(1) = R1 + R2
    A(3) = R1 - R2
C
    R1 = B(1) + B(3)
    T2 = B(1) - B(3)
    R2 = B(2) + B(4)
    B(1) = R1 + R2
    B(3) = R1 - R2
C
    R1 = A(2) - A(4)
    R2 = B(2) - B(4)
C
    A(2) = T1 + R2
    A(4) = T1 - R2
    B(2) = T2 - R1
    B(4) = T2 + R1
C
    GOTO 20
C
C----------------WFTA N=5--------------------------------
C
  105   R1 = A(2) + A(5)
    R4 = A(2) - A(5)
    R3 = A(3) + A(4)
    R2 = A(3) - A(4)
C
    T = (R1 - R3) * C54
    R1 = R1 + R3
    A(1) = A(1) + R1
    R1   = A(1) + R1 * C55
C
    R3 = R1 - T
    R1 = R1 + T
C
    T = (R4 + R2) * C51
    R4 =  T + R4 * C52
    R2 =  T + R2 * C53
C
    S1 = B(2) + B(5)
    S4 = B(2) - B(5)
    S3 = B(3) + B(4)
    S2 = B(3) - B(4)
C
    T = (S1 - S3) * C54
    S1 = S1 + S3
    B(1) = B(1) + S1
    S1   = B(1) + S1 * C55
C
    S3 =  S1 - T
    S1 =  S1 + T
C
    T = (S4 + S2) * C51
    S4 =  T + S4 * C52
    S2 =  T + S2 * C53
C
    A(2) = R1 + S2
    A(5) = R1 - S2
    A(3) = R3 - S4
    A(4) = R3 + S4
C
    B(2) = S1 - R2
    B(5) = S1 + R2
    B(3) = S3 + R4
    B(4) = S3 - R4
C
    GOTO 20
C-----------------WFTA N=7--------------------------
C
  107   R1 = A(2) + A(7)
    R6 = A(2) - A(7)
    S1 = B(2) + B(7)
    S6 = B(2) - B(7)
    R2 = A(3) + A(6)
    R5 = A(3) - A(6)
    S2 = B(3) + B(6)
    S5 = B(3) - B(6)
    R3 = A(4) + A(5)
    R4 = A(4) - A(5)
    S3 = B(4) + B(5)
    S4 = B(4) - B(5)
C
    T3 = (R1 - R2) * C74
    T  = (R1 - R3) * C72
    R1 = R1 + R2 + R3
    A(1) = A(1) + R1
    R1   = A(1) + R1 * C71
    R2 =(R3 - R2) * C73
    R3 = R1 - T + R2
    R2 = R1 - R2 - T3
    R1 = R1 + T + T3
    T = (R6 - R5) * C78
    T3 =(R6 + R4) * C76
    R6 =(R6 + R5 - R4) * C75
    R5 =(R5 + R4) * C77
    R4 = R6 - T3 + R5
    R5 = R6 - R5 - T
    R6 = R6 + T3 + T
C
    T3 = (S1 - S2) * C74
    T  = (S1 - S3) * C72
    S1 =  S1 + S2 + S3
    B(1) = B(1) + S1
    S1   = B(1) + S1 * C71
    S2 =(S3 - S2) * C73
    S3 = S1 - T  + S2
    S2 = S1 - S2 - T3
    S1 = S1 + T  + T3
    T  = (S6 - S5) * C78
    T3 = (S6 + S4) * C76
    S6 = (S6 + S5 - S4) * C75
    S5 = (S5 + S4) * C77
    S4 = S6 - T3 + S5
    S5 = S6 - S5 - T
    S6 = S6 + T3 + T
C
    A(2) = R3 + S4
    A(7) = R3 - S4
    A(3) = R1 + S6
    A(6) = R1 - S6
    A(4) = R2 - S5
    A(5) = R2 + S5
    B(4) = S2 + R5
    B(5) = S2 - R5
    B(2) = S3 - R4
    B(7) = S3 + R4
    B(3) = S1 - R6
    B(6) = S1 + R6
C
   20        IT   = J
             DO 31 L=1, N1
                I(L) = IT
            X(IT) = A(L)
            Y(IT) = B(L)
                IT = IT + N2
                IF (IT.GT.N)  IT = IT - N
   31        CONTINUE
   15     CONTINUE
   10   CONTINUE
C
C-----------------UNSCRAMBLING----------------------
C
    L = 1
    DO 2 K=1, N
       A(K) = X(L)
       B(K) = Y(L)
           L = L + UNSC
           IF (L.GT.N)  L = L - N
    2   CONTINUE
    RETURN
    END
C