C---------------------------------------------------------C
C     An Extended Split-Radix DIF FFT                     C
C     Complex input and output in data arrays X and Y     C
C     Length is N = 2**M                                  C
C                                                         C
C     D. Takahashi, University of Tsukuba, Sep. 2002      C
C                                                         C
C     Reference:                                          C
C         D. Takahashi, "An extended split-radix FFT      C
C         algorithm," IEEE Signal Processing Lett.,       C
C         vol. 8, pp. 145-147, May 2001.                  C
C---------------------------------------------------------C
      SUBROUTINE ESRFFT(X,Y,IBETA,N,M,IV)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 X(N),Y(N)
      INTEGER*4 IBETA(N/2)
      PARAMETER (C21=0.70710678118654752D0,
     +           TWOPI=6.28318530717958647D0)
C-------Beta vector table-------
      IBETA(1)=1
      IBETA(2)=0
      IBETA(3)=0
      IBETA(4)=0
      ID=1
      DO 30 K=3,M-1
        DO 20 J=1,4
          IS=(J+3)*ID+1
          DO 10 I=IS,IS+ID-1
            IBETA(I)=IBETA(I-IS+1)
   10     CONTINUE
   20   CONTINUE
        ID=2*ID
   30 CONTINUE
C
      IF (IV .EQ. -1) THEN
        DO 40 I=1,N
          Y(I)=-Y(I)
   40   CONTINUE
      END IF
C-------L shaped butterflies-------
      L=1
      N2=N
      DO 70 K=1,M-2
        N8=N2/8
        E=-TWOPI/DBLE(N2)
        A=0.0D0
        DO 60 J=1,N8
          A3=3.0D0*A
          A5=5.0D0*A
          A7=7.0D0*A
          CC1=DCOS(A)
          SS1=DSIN(A)
          CC3=DCOS(A3)
          SS3=DSIN(A3)
          CC5=DCOS(A5)
          SS5=DSIN(A5)
          CC7=DCOS(A7)
          SS7=DSIN(A7)
          A=DBLE(J)*E
          DO 50 I=1,L
            IF (IBETA(I) .EQ. 1) THEN
              I0=(I-1)*N2+J
              I1=I0+N8
              I2=I1+N8
              I3=I2+N8
              I4=I3+N8
              I5=I4+N8
              I6=I5+N8
              I7=I6+N8
              X0=X(I0)-X(I4)
              Y0=Y(I0)-Y(I4)
              X1=X(I1)-X(I5)
              Y1=Y(I1)-Y(I5)
              X2=Y(I2)-Y(I6)
              Y2=X(I6)-X(I2)
              X3=X(I3)-X(I7)
              Y3=Y(I3)-Y(I7)
              X(I0)=X(I0)+X(I4)
              Y(I0)=Y(I0)+Y(I4)
              X(I1)=X(I1)+X(I5)
              Y(I1)=Y(I1)+Y(I5)
              X(I2)=X(I2)+X(I6)
              Y(I2)=Y(I2)+Y(I6)
              X(I3)=X(I3)+X(I7)
              Y(I3)=Y(I3)+Y(I7)
              U0=X0+C21*(X1-X3)
              V0=Y0+C21*(Y1-Y3)
              U1=X0-C21*(X1-X3)
              V1=Y0-C21*(Y1-Y3)
              U2=X2+C21*(Y1+Y3)
              V2=Y2-C21*(X1+X3)
              U3=X2-C21*(Y1+Y3)
              V3=Y2+C21*(X1+X3)
              X(I4)=CC1*(U0+U2)-SS1*(V0+V2)
              Y(I4)=CC1*(V0+V2)+SS1*(U0+U2)
              X(I5)=CC5*(U1+U3)-SS5*(V1+V3)
              Y(I5)=CC5*(V1+V3)+SS5*(U1+U3)
              X(I6)=CC3*(U1-U3)-SS3*(V1-V3)
              Y(I6)=CC3*(V1-V3)+SS3*(U1-U3)
              X(I7)=CC7*(U0-U2)-SS7*(V0-V2)
              Y(I7)=CC7*(V0-V2)+SS7*(U0-U2)
            END IF
   50     CONTINUE
   60   CONTINUE
        L=2*L
        N2=N2/2
   70 CONTINUE
C-------Length four butterflies-------
      DO 80 I=1,N/4
        IF (IBETA(I) .EQ. 1) THEN
          I0=4*I-3
          I1=I0+1
          I2=I1+1
          I3=I2+1
          X0=X(I0)-X(I2)
          Y0=Y(I0)-Y(I2)
          X1=Y(I1)-Y(I3)
          Y1=X(I3)-X(I1)
          X(I0)=X(I0)+X(I2)
          Y(I0)=Y(I0)+Y(I2)
          X(I1)=X(I1)+X(I3)
          Y(I1)=Y(I1)+Y(I3)
          X(I2)=X0+X1
          Y(I2)=Y0+Y1
          X(I3)=X0-X1
          Y(I3)=Y0-Y1
        END IF
   80 CONTINUE
C-------Length two butterflies-------
      DO 90 I=1,N/2
        IF (IBETA(I) .EQ. 1) THEN
          I0=2*I-1
          I1=I0+1
          X0=X(I0)-X(I1)
          Y0=Y(I0)-Y(I1)
          X(I0)=X(I0)+X(I1)
          Y(I0)=Y(I0)+Y(I1)
          X(I1)=X0
          Y(I1)=Y0
        END IF
   90 CONTINUE
C-------Digit reverse counter-------
      J=1
      DO 110 I=1,N-1
        IF (I .LT. J) THEN
          X0=X(I)
          Y0=Y(I)
          X(I)=X(J)
          Y(I)=Y(J)
          X(J)=X0
          Y(J)=Y0
        END IF
        K=N/2
  100   IF (K .LT. J) THEN
          J=J-K
          K=K/2
          GO TO 100
        END IF
        J=J+K
  110 CONTINUE
C
      IF (IV .EQ. -1) THEN
        DO 120 I=1,N
          X(I)=X(I)/DBLE(N)
          Y(I)=-Y(I)/DBLE(N)
  120   CONTINUE
      END IF
      RETURN
      END