CTITLEMFFT42 -- FAST FOURIER TRANSFORM 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR DAVE THOMPSON 00000200 CA DESIGNER DAVE THOMPSON 00000300 CA LANGUAGE FORTRAN 00000400 CA SYSTEM IBM AND CRAY 00000500 CA WRITTEN 00000600 CA REVISED 00000700 CA 00000800 CA 00000900 CA CALL MFFT42 (N2, A, B, ISN) 00001000 CA 00001100 CA IN N2 I4 2-EXPONENT FOR LENGTH OF DATA TO BE TRANSFORMED 00001200 CA IN A R4 ARRAY CONTAINING REAL COMPONENT OF INPUT DATA 00001300 CA CONTENTS OF A ARE REPLACED BY REAL COMPONENT OF 00001400 CA TRANSFORM ON OUTPUT 00001500 CA IN B R4 ARRAY CONTAINING IMAGINARY COMPONENT OF INPUT DATA00001600 CA CONTENTS OF B ARE REPLACED BY IMAGINARY COMPONENT 00001700 CA TRANSFORM ON OUTPUT. 00001800 CA IN ISN I4 = -1 FOR DIRECT TRANSFORM 00001900 CA = 1 FOR NORMALIZED INVERSE TRANSFORM 00002000 CA = 2 FOR UNNORMALIZED INVERSE TRANSFORM 00002100 CA 00002200 CA THIS SUBROUTINE DOES A FAST FOURIER TRANSFORM 00002300 CA 00002400 CA CAUTIONS: 00002500 CA A. RADIX-4 FFT WITH FINAL RADIX-2 PASS IF NECESSARY. 00002600 CA B. NO SIN-COS TABLES REQUIRED. 00002700 CA C. TRIG FUNCTION CALCULATION COMPENSATED FOR TRUNCATED 00002800 CA ARITHMETIC (SINGLETON). 00002900 CA 00003000 CA REFERENCES: 00003100 CA A. BASIC ALGORITHM DERIVED FROM SUBROUTINE FRXFM BY 00003200 CA LANGDON AND SANDE, PRINCETON, 1965. 00003300 CA B. TRIG FUNCTION GENERATION FROM SUBROUTINE FFT 00003400 CA BY SINGLETON, IEEE TRANS. AUDIO AND ELECTRO- 00003500 CA ACOUSTICS, JUNE 1969, P. 100. 00003600 CA C. REVERSE BINARY ORDERING FROM SUBROUTINE NLOGN 00003700 CA BY ROBINSON, 'MULTICHANNEL TIME SERIES 00003800 CA ANALYSIS WITH DIG. COMP. PROGS.', P. 63. 00003900 CA 00004000 CA 4. SPECIAL NOTES 00004100 CA A. EXTERNAL ROUTINES CALLED - NONE 00004200 CA B. CALLING PROGRAM DIMENSIONS - A(L),B(L) 00004300 C 00004400 C 00004500 C 00004600 SUBROUTINE MFFT42 (N2,A,B,ISN) 00004700 DIMENSION A(1),B(1),M(16) 00004800 L=2**N2 00004900 N4=N2/2 00005000 RADF=4.0*ATAN(1.0) 00005100 IF (ISN.LT.0) RADF=-RADF 00005200 IF (N4.EQ.0) GO TO 60 00005300 00005400 C DO RADIX-4 TRANSFORM 00005500 00005600 KS=L 00005700 DO 50 I=1,N4 00005800 LEN=KS 00005900 KS=KS/4 00006000 SD=RADF/FLOAT(LEN) 00006100 CD=2.0*SIN(SD)**2 00006200 SD=SIN(SD+SD) 00006300 C1=1.0 00006400 S1=0.0 00006500 DO 50 J=1,KS 00006600 IF (J.EQ.1) GO TO 10 00006700 C2=C1-(CD*C1+SD*S1) 00006800 S1=(SD*C1-CD*S1)+S1 00006900 C1=0.5/(C2*C2+S1*S1)+0.5 00007000 S1=C1*S1 00007100 C1=C1*C2 00007200 C2=C1*C1-S1*S1 00007300 S2=2.0*C1*S1 00007400 C3=C2*C1-S2*S1 00007500 S3=C2*S1+S2*C1 00007600 10 JL=J-LEN 00007700 DO 50 ISEQ=LEN,L,LEN 00007800 K1=ISEQ+JL 00007900 K2=K1+KS 00008000 K3=K2+KS 00008100 K4=K3+KS 00008200 AKP=A(K1)+A(K3) 00008300 AKM=A(K1)-A(K3) 00008400 AJP=A(K2)+A(K4) 00008500 AJM=A(K2)-A(K4) 00008600 A(K1)=AKP+AJP 00008700 AJP=AKP-AJP 00008800 BKP=B(K1)+B(K3) 00008900 BKM=B(K1)-B(K3) 00009000 BJP=B(K2)+B(K4) 00009100 BJM=B(K2)-B(K4) 00009200 B(K1)=BKP+BJP 00009300 BJP=BKP-BJP 00009400 IF (ISN.LT.0) GO TO 30 00009500 AKP=AKM-BJM 00009600 AKM=AKM+BJM 00009700 BKP=BKM+AJM 00009800 BKM=BKM-AJM 00009900 IF (J.EQ.1) GO TO 40 00010000 20 A(K3)=AKP*C1-BKP*S1 00010100 B(K3)=AKP*S1+BKP*C1 00010200 A(K2)=AJP*C2-BJP*S2 00010300 B(K2)=AJP*S2+BJP*C2 00010400 A(K4)=AKM*C3-BKM*S3 00010500 B(K4)=AKM*S3+BKM*C3 00010600 GO TO 50 00010700 30 AKP=AKM+BJM 00010800 AKM=AKM-BJM 00010900 BKP=BKM-AJM 00011000 BKM=BKM+AJM 00011100 IF (J.GT.1) GO TO 20 00011200 40 A(K3)=AKP 00011300 B(K3)=BKP 00011400 A(K2)=AJP 00011500 B(K2)=BJP 00011600 A(K4)=AKM 00011700 B(K4)=BKM 00011800 50 CONTINUE 00011900 60 IF (N2.EQ.N4+N4) GO TO 80 00012000 00012100 C DO FINAL RADIX-2 PASS IF REQUIRED 00012200 00012300 DO 70 J=1,L,2 00012400 T=A(J)+A(J+1) 00012500 A(J+1)=A(J)-A(J+1) 00012600 A(J)=T 00012700 T=B(J)+B(J+1) 00012800 B(J+1)=B(J)-B(J+1) 00012900 70 B(J)=T 00013000 00013100 C REARRANGE OUTPUT IN REVERSE BINARY ORDER 00013200 00013300 80 K=L 00013400 DO 90 I=1,N2 00013500 K=K/2 00013600 90 M(I)=K 00013700 K=0 00013800 DO 120 J=1,L 00013900 IF (K.LT.J) GO TO 100 00014000 T=A(J) 00014100 A(J)=A(K+1) 00014200 A(K+1)=T 00014300 T=B(J) 00014400 B(J)=B(K+1) 00014500 B(K+1)=T 00014600 100 DO 110 I=1,N2 00014700 II=I 00014800 IF (K.LT.M(I)) GO TO 120 00014900 110 K=K-M(I) 00015000 120 K=K+M(II) 00015100 IF (ISN.NE.1) RETURN 00015200 00015300 C NORMALIZE OUTPUT IF ISN = 1 00015400 00015500 FLI=1./FLOAT(L) 00015600 DO 130 I=1,L 00015700 A(I)=FLI*A(I) 00015800 130 B(I)=FLI*B(I) 00015900 RETURN 00016000 END 00016100