CTITLEMFORAX -- DOUBLE PRECISION FFT 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. D. THOMPSON 00000020 CA DESIGNER D. D. THOMPSON 00000030 CA LANGUAGE FORTRAN 77 00000040 CA SYSTEM IBM & CRAY 00000041 CA WRITTEN 1972 00000050 C REVISED MO-DA-YR 00000060 C REVISED 05-20-85 TWH. ADAPT TO RUN ON IBM & 00000061 C CRAY SYSTEMS. 00000062 C 00000070 CA 00000080 CA 00000090 CA CALL MFORAX (N,X,Y,ISIGN) 00000100 CA INPUT N = LOG BASE 2 OF THE FFT LENGTH I4 00000110 CA IN/OUT X = REAL DATA ARRAY R8 00000120 CA IN/OUT Y = IMAG. DATA ARRAY R8 00000130 CA INPUT ISIGN = INPUT FLAG I4 00000140 CA -1 FOR DIRECT TRANSFORM 00000150 CA +1 FOR NORMALIZED INVERSE TRANSFORM 00000160 CA +2 FOR UNNORMALIZED INVERSE TRANSFORM 00000170 CA 00000180 CA 00000190 CA THIS ROUTINE PERFORMS EITHER FORWARD OR REVERSE 00000200 CA FFT IN PLACE WITH DOUBLE PRECISION. 00000210 CA 00000220 CA 00000230 C 00000240 C SUBROUTINES CALLED: NONE 00000250 C 00000260 C EJECT 00000270 C 00000280 SUBROUTINE MFORAX(N,X,Y,ISIGN) 00000290 CAEND 00000300 C 00000310 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00000320 DIMENSIONX(1), Y(1), M( 00000330 *16) 00000340 LX=2**N 00000350 FLXI=1.D0/LX 00000360 FP=DATAN(1.D0)*FLXI*8.D0 00000370 IF(ISIGN.EQ.-1) FP=-FP 00000380 C 00000390 DO 10 00000400 * I=1,N 00000410 C 00000420 10 M(I)=2**(N-I) 00000430 C 00000440 DO 40 00000450 * L=1,N 00000460 NB=2**(L-1) 00000470 LB=LX/NB 00000480 LBH=LB/2 00000490 K=0 00000500 C 00000510 DO 40 00000520 * IB=1,NB 00000530 FK=K 00000540 V=FK*FP 00000550 W1=DCOS(V) 00000560 W2=DSIN(V) 00000570 IS=LB*(IB-1) 00000580 C 00000590 DO 20 00000600 * I=1,LBH 00000610 J=IS+I 00000620 JH=J+LBH 00000630 Q1=X(JH)*W1-Y(JH)*W2 00000640 Q2=X(JH)*W2+Y(JH)*W1 00000650 X(JH)=X(J)-Q1 00000660 Y(JH)=Y(J)-Q2 00000670 X(J)=X(J)+Q1 00000680 C 00000690 20 Y(J)=Y(J)+Q2 00000700 C 00000710 DO 30 00000720 * I=2,N 00000730 II=I 00000740 IF (K.LT.M(I)) GO TO 40 00000750 C 00000760 30 K=K-M(I) 00000770 C 00000780 40 K=K+M(II) 00000790 K=0 00000800 C 00000810 DO 70 00000820 * J=1,LX 00000830 IF (K.LT.J) GO TO 50 00000840 H1=X(J) 00000850 H2=Y(J) 00000860 X(J)=X(K+1) 00000870 Y(J)=Y(K+1) 00000880 X(K+1)=H1 00000890 Y(K+1)=H2 00000900 C 00000910 50 DO 60 00000920 * I=1,N 00000930 II=I 00000940 IF (K.LT.M(I)) GO TO 70 00000950 C 00000960 60 K=K-M(I) 00000970 C 00000980 70 K=K+M(II) 00000990 C 00001000 IF(ISIGN.EQ.1) GO TO 80 00001010 RETURN 00001020 C 00001030 80 DO 90 00001040 * I=1,LX 00001050 X(I)=FLXI*X(I) 00001060 C 00001070 90 Y(I)=FLXI*Y(I) 00001080 C 00001090 RETURN 00001100 END 00001110