CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C C
C Subroutine FASTF C
C C
C Radix 4 complex discrete fast fourier transform C
C C
C XREAL = Array which on input contains REAL part of data for transformation C
C and on output gives REAL part of the result, type REAL, dimension C
C IABS(ISIZE) or greater. C
C C
C XIMAG = Array which on input contains IMAGINARY part of data for C
C transformation and on output gives imaginary part of the result, type real, C
C dimension iabs(isize) or greater. C
C C
C ISIZE = INTEGER variable or constant specifying length and type of C
C transform. The length is IABS(ISIZE) which must be a power of 2. Minimum C
C length 4, maximum 2**20. If ISIZE is positive the forward transform is C
C calculated, if negative the inverse transform is found. C
C C
C This subroutine was written by D.Monro. Imperial College, London. 1971 C
C C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE FASTF(XREAL,XIMAG,ISIZE)
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION XREAL(32384),XIMAG(32384)
INTEGER L(20)
EQUIVALENCE(L1,L(1)),(L2,L(2)),(L3,L(3)),(L4,L(4)),
1 (L5,L(5)),(L6,L(6)),(L7,L(7)),(L8,L(8)),(L9,L(9)),
2 (L10,L(10)),(L11,L(11)),(L12,L(12)),(L13,L(13)),
3 (L14,L(14)),(L15,L(15)),(L16,L(16)),(L17,L(17)),
4 (L18,L(18)),(L19,L(19)),(L20,L(20))
PIE=4.*ATAN(1.)
N=IABS(ISIZE)
IF(N-4)24,1,1
C SET UP INITIAL VALUES OF TRANSFORM SPLIT
1 IFACA=N/4
IF(ISIZE)2,4,4
C IF THIS IS TO BE AN INVERSE TRANSFORM, CONJUGATE THE DATA
2 DO 3 K=1,N
3 XIMAG(K)=-XIMAG(K)
4 ITIME=0
5 IFCAB=IFACA*4
ITIME=ITIME+2
C DO THE TRANSFORMS REQUIRED BY THIS STAGE
Z=PIE/FLOAT(IFCAB)
BCOS=-2.*(SIN(Z)**2)
BSIN=SIN(2.*Z)
CW1=1.
SW1=0.
DO 10 LITLA=1,IFACA
DO 8 I0=LITLA,N,IFCAB
C THIS IS THE MAIN CALCULATION OF RADIX 4 TRANSFORMS
I1=I0+IFACA
I2=I1+IFACA
I3=I2+IFACA
XS0=XREAL(I0)+XREAL(I2)
XS1=XREAL(I0)-XREAL(I2)
YS0=XIMAG(I0)+XIMAG(I2)
YS1=XIMAG(I0)-XIMAG(I2)
XS2=XREAL(I1)+XREAL(I3)
XS3=XREAL(I1)-XREAL(I3)
YS2=XIMAG(I1)+XIMAG(I3)
YS3=XIMAG(I1)-XIMAG(I3)
XREAL(I0)=XS0+XS2
XIMAG(I0)=YS0+YS2
X1=XS1+YS3
Y1=YS1-XS3
X2=XS0-XS2
Y2=YS0-YS2
X3=XS1-YS3
Y3=YS1+XS3
IF(LITLA-1)24,6,7
6 XREAL(I2)=X1
XIMAG(I2)=Y1
XREAL(I1)=X2
XIMAG(I1)=Y2
XREAL(I3)=X3
XIMAG(I3)=Y3
GO TO 8
C MULTIPLY BY TWIDDLE FACTORS IF REQUIRED
7 XREAL(I2)=X1*CW1+Y1*SW1
XIMAG(I2)=Y1*CW1-X1*SW1
XREAL(I1)=X2*CW2+Y2*SW2
XIMAG(I1)=Y2*CW2-X2*SW2
XREAL(I3)=X3*CW3+Y3*SW3
XIMAG(I3)=Y3*CW3-X3*SW3
8 CONTINUE
IF(LITLA-IFACA)9,10,24
C CALCULATE A NEW SET OF TWIDDLE FACTORS
9 Z=CW1*BCOS-SW1*BSIN+CW1
SW1=BCOS*SW1+BSIN*CW1+SW1
TEMPR=1.5-0.5*(Z*Z+SW1*SW1)
CW1=Z*TEMPR
SW1=SW1*TEMPR
CW2=CW1*CW1-SW1*SW1
SW2=2.*CW1*SW1
CW3=CW1*CW2-SW1*SW2
SW3=CW1*SW2+CW2*SW1
10 CONTINUE
IF(IFACA-1)14,14,11
C SET UP THE TRANSFORM SPLIT FOR THE NEXT STAGE
11 IFACA=IFACA/4
IF(IFACA)24,12,5
C THIS IS THE CALCULATION OF A RADIX TWO STAGE
12 DO 13 K=1,N,2
TEMPR=XREAL(K)+XREAL(K+1)
XREAL(K+1)=XREAL(K)-XREAL(K+1)
XREAL(K)=TEMPR
TEMPR=XIMAG(K)+XIMAG(K+1)
XIMAG(K+1)=XIMAG(K)-XIMAG(K+1)
13 XIMAG(K)=TEMPR
ITIME=ITIME+1
14 IF(ISIZE)15,19,17
C IF THIS WAS AN INVERSE TRANSFORM, CONJUGATE THE RESULT
15 DO 16 K=1,N
16 XIMAG(K)=-XIMAG(K)
GO TO 19
C IF THIS WAS A FORWARD TRANSFORM, SCALE THE RESULT
17 Z=1./FLOAT(N)
DO 18 K=1,N
XREAL(K)=XREAL(K)*Z
18 XIMAG(K)=XIMAG(K)*Z
C UNSCRAMBLE THE RESULT
19 I1=20-ITIME
DO 20 K=1,I1
20 L(K)=1
II=1
I1=I1+1
DO 21 K=I1,20
II=II*2
21 L(K)=II
II=1
DO 23 J1=1,L1
DO 23 J2=J1,L2,L1
DO 23 J3=J2,L3,L2
DO 23 J4=J3,L4,L3
DO 23 J5=J4,L5,L4
DO 23 J6=J5,L6,L5
DO 23 J7=J6,L7,L6
DO 23 J8=J7,L8,L7
DO 23 J9=J8,L9,L8
DO 23 J10=J9,L10,L9
DO 23 J11=J10,L11,L10
DO 23 J12=J11,L12,L11
DO 23 J13=J12,L13,L12
DO 23 J14=J13,L14,L13
DO 23 J15=J14,L15,L14
DO 23 J16=J15,L16,L15
DO 23 J17=J16,L17,L16
DO 23 J18=J17,L18,L17
DO 23 J19=J18,L19,L18
DO 23 J20=J19,L20,L19
IF(II-J20)22,23,23
22 TEMPR=XREAL(II)
XREAL(II)=XREAL(J20)
XREAL(J20)=TEMPR
TEMPR=XIMAG(II)
XIMAG(II)=XIMAG(J20)
XIMAG(J20)=TEMPR
23 II=II+1
24 RETURN
END