Fairfield33
Programmer
This document on Integer FFT's is not compiling for me. I am a beginner in Fortran and am having some problems. I would greatly appreciate any Help with this document.
SUBROUTINE FASTI (IRE,IMG,ISIZE)
C PROVIDE INTEGER WORD LENGTH AND PIE TO AVAILABLE PRECISION
DATA NBITTS/16/,PIE/3.14159/
DIMENSION IRE (1) ,IMG(1)
IBIG=2**(NBITS-2)
IBG2=IBIG/2
XBIG=FLOAT(IBIG)
N=IABS(ISIZE)
C SET UP INITIAL TRIG FUNCTIONS
Z=PIE/FLOAT(N)
IAS=SIN(Z)*XBIG+0.5
IAS=COS(Z)*XBIG+0.5
C SET UP INTIAL VALUES OF TRANSORM SPLIT
IFA = N/2
C IF AN INVERSE TRANSFORM, CONJUGATE THE DATA
IF(ISIZE)1,3,3
1 DO 2 K=1,N
2 IMG(K)=-IMG(K)
3 I TIME=0
4 IFAB=IFA*2
ITIME=ITIME+1
C DOUBLE THE BASIC ANGLE
IBSIN = 2*MPY(IAS,IAC)
IBCOS=MPY(IAS,IAS)
IAS=IBSIN
IAC=(IBG2-IBCOS)*2
CALL COMP(IAS,IAC)
C MAIN CALCLATION OF RADIX 2 TRANSFORMS
ISW=0
ICW=IBIG
DO 11 LA = 1,IFA
DO 9 IP=LA,N,IFAB
IQ=IP+IFA
C EXTRACT DATA AND SCALE
K1=IRE(IP)
K2=TRE(IQ)
L1=IMG(IP)
L2=IMG(IQ)
IF(ISIZE)6,18,5
5 K1=K1/2
K2=K2/2
L1=L1/2
L2=L2/2
C FLUTTER BY
6 IRE(IP)=K1+K2
K1=K1-K2
IMG(IP)=L1+L2
L1=L1-L2
IF(LA-1)18,7,8
7 IRE(IQ)=K1
IMG(IQ)=L1
GO TO 9
C MULTIPLY BY TWIDDLE FACTORS IF REQUIRED
8 IRE(IQ)=MPY(K1,ICW)+MPY(L1,ISW)
IMG(IQ)=MPY(L1,ICW)+MPY(K1,ISW)
9 CONTINUE
IF(LA-IFA)10,11,18
C CALCULATE A NEW SET OF TWIDDLE FACTORS
10 IT=ICW-MPY(IBSIN,ISW)-2*MPY(IBCOS,ICW)
ISW=ISW+MPY(IBSIN,ICW)-2*MPY(IBCOS,ISW)
ICW=IT
CALL COMP(ISW,ICW)
11 CONTINUE
IFA=IFA/2
IF(IFA)12,12,4
C TRANSFORM COMPLETE - CONJUGATE RESULT IF INVERSE
12 IF(ISIZE)13,18,15
13 DO 14 K=1,N
14 IMG(K)=-IMG(K)
15 DO 17 J20=1,N
IZ=J20-1
CALL BITRV(IZ,II,ITIME)
IF(II-J20)17,17,16
16 II=II+1
IT=IRE(II)
IRE(II)=IRE(J20)
IRE(J20)=IT
IT=IMG(II)
IMG(II)=IMG(J20)
IMG(J20)=IT
17 CONTINUE
18 RETURN
END
FUNCTION MPY(I,J)
C FUNCTION TO MULTIPLY FIXED POINT FRACTIONS
C THIS VERSION INCORPORATES ROUNDING
C
C PROVIDE CONSTANTS FOR DECOMPOSISTION FO PRODUCT
DATA IMAD,IMID/64,128/
ILEF=I/IMID
IRIT=I-ILEF*IMID
JLEF=J/IMID
JRIT=J-JLEF*IMID
C DO NOT REARRANGE THE FOLLOWING MULTIPLICATION
II=(IRIT*JRIT)/IMID+IRIT*JLEF+JRIT*ILEF
MPY= (II+ISIGN(IMAD,II))/IMID+ILEF*JLEF
RETURN
END
FUNCTION MPY(I,J)
C
C FUNCTION TO MULITPLY FIXED POINT FRACTIONS
C THIS VERSION DOES NOT ROUND
C
C PROVED CONSTANTS FOR DECOMPOSITION OF PRODUCT
DATA IMID/128/
ILEF=I/IMID
IRIT=I-ILEF*IMID
JLEF=J/IMID
JRIT=J-JLEF*IMID
C DO NOT REARRANGE THE FOLLOWING MULTIPLICATION
MPY=((IRIT*JRIT)/IMID+IRIT*JLEF+JRIT*ILEF)/IMID+ILEF*JLEF
RETURN
END
SUBROUTINE BITRV(I,J,NBITS)
C
C SUBROUTINE WHICH RETURNS NBITS BITS OF I REVERSED IN J.
C
IZ=I
IJ=0
DO 1 K=1,NBITS
C DO NOT REARRANGE THE FOLLOWING STATEMENT
IJ=IJ+IJ+IZ/2*2
1 IZ=IZ/2
J=IJ
RETURN
END
SUBROUTINE FASTI (IRE,IMG,ISIZE)
C PROVIDE INTEGER WORD LENGTH AND PIE TO AVAILABLE PRECISION
DATA NBITTS/16/,PIE/3.14159/
DIMENSION IRE (1) ,IMG(1)
IBIG=2**(NBITS-2)
IBG2=IBIG/2
XBIG=FLOAT(IBIG)
N=IABS(ISIZE)
C SET UP INITIAL TRIG FUNCTIONS
Z=PIE/FLOAT(N)
IAS=SIN(Z)*XBIG+0.5
IAS=COS(Z)*XBIG+0.5
C SET UP INTIAL VALUES OF TRANSORM SPLIT
IFA = N/2
C IF AN INVERSE TRANSFORM, CONJUGATE THE DATA
IF(ISIZE)1,3,3
1 DO 2 K=1,N
2 IMG(K)=-IMG(K)
3 I TIME=0
4 IFAB=IFA*2
ITIME=ITIME+1
C DOUBLE THE BASIC ANGLE
IBSIN = 2*MPY(IAS,IAC)
IBCOS=MPY(IAS,IAS)
IAS=IBSIN
IAC=(IBG2-IBCOS)*2
CALL COMP(IAS,IAC)
C MAIN CALCLATION OF RADIX 2 TRANSFORMS
ISW=0
ICW=IBIG
DO 11 LA = 1,IFA
DO 9 IP=LA,N,IFAB
IQ=IP+IFA
C EXTRACT DATA AND SCALE
K1=IRE(IP)
K2=TRE(IQ)
L1=IMG(IP)
L2=IMG(IQ)
IF(ISIZE)6,18,5
5 K1=K1/2
K2=K2/2
L1=L1/2
L2=L2/2
C FLUTTER BY
6 IRE(IP)=K1+K2
K1=K1-K2
IMG(IP)=L1+L2
L1=L1-L2
IF(LA-1)18,7,8
7 IRE(IQ)=K1
IMG(IQ)=L1
GO TO 9
C MULTIPLY BY TWIDDLE FACTORS IF REQUIRED
8 IRE(IQ)=MPY(K1,ICW)+MPY(L1,ISW)
IMG(IQ)=MPY(L1,ICW)+MPY(K1,ISW)
9 CONTINUE
IF(LA-IFA)10,11,18
C CALCULATE A NEW SET OF TWIDDLE FACTORS
10 IT=ICW-MPY(IBSIN,ISW)-2*MPY(IBCOS,ICW)
ISW=ISW+MPY(IBSIN,ICW)-2*MPY(IBCOS,ISW)
ICW=IT
CALL COMP(ISW,ICW)
11 CONTINUE
IFA=IFA/2
IF(IFA)12,12,4
C TRANSFORM COMPLETE - CONJUGATE RESULT IF INVERSE
12 IF(ISIZE)13,18,15
13 DO 14 K=1,N
14 IMG(K)=-IMG(K)
15 DO 17 J20=1,N
IZ=J20-1
CALL BITRV(IZ,II,ITIME)
IF(II-J20)17,17,16
16 II=II+1
IT=IRE(II)
IRE(II)=IRE(J20)
IRE(J20)=IT
IT=IMG(II)
IMG(II)=IMG(J20)
IMG(J20)=IT
17 CONTINUE
18 RETURN
END
FUNCTION MPY(I,J)
C FUNCTION TO MULTIPLY FIXED POINT FRACTIONS
C THIS VERSION INCORPORATES ROUNDING
C
C PROVIDE CONSTANTS FOR DECOMPOSISTION FO PRODUCT
DATA IMAD,IMID/64,128/
ILEF=I/IMID
IRIT=I-ILEF*IMID
JLEF=J/IMID
JRIT=J-JLEF*IMID
C DO NOT REARRANGE THE FOLLOWING MULTIPLICATION
II=(IRIT*JRIT)/IMID+IRIT*JLEF+JRIT*ILEF
MPY= (II+ISIGN(IMAD,II))/IMID+ILEF*JLEF
RETURN
END
FUNCTION MPY(I,J)
C
C FUNCTION TO MULITPLY FIXED POINT FRACTIONS
C THIS VERSION DOES NOT ROUND
C
C PROVED CONSTANTS FOR DECOMPOSITION OF PRODUCT
DATA IMID/128/
ILEF=I/IMID
IRIT=I-ILEF*IMID
JLEF=J/IMID
JRIT=J-JLEF*IMID
C DO NOT REARRANGE THE FOLLOWING MULTIPLICATION
MPY=((IRIT*JRIT)/IMID+IRIT*JLEF+JRIT*ILEF)/IMID+ILEF*JLEF
RETURN
END
SUBROUTINE BITRV(I,J,NBITS)
C
C SUBROUTINE WHICH RETURNS NBITS BITS OF I REVERSED IN J.
C
IZ=I
IJ=0
DO 1 K=1,NBITS
C DO NOT REARRANGE THE FOLLOWING STATEMENT
IJ=IJ+IJ+IZ/2*2
1 IZ=IZ/2
J=IJ
RETURN
END