Two-Dimensional Mixed Radix Mass Storage Fourier Transform
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: FFT2D
C AUTHOR: RICHARD C. SINGLETON
C SRI INTERNATIONAL, MENLO PARK, CALIFORNIA 94025
C INPUT: NONE
C-----------------------------------------------------------------------
C
COMMON /FFTCOM/ LUN
C
C SET UP DYNAMIC STORAGE ALLOCATION, AS IN PORT
C
COMMON /CSTAK/ DSTAK(2500)
C
COMMON /MASS/ RMAS(54450)
C
C NOTE: MINICOMPUTERS MAY NOT BE ABLE TO ACCOMMODATE
C THE LARGE ARRAY USED IN THIS TEST PROGRAM TO SIMULATE
C MASS STORAGE. IN SUCH CASES, TRUE MASS STORAGE SHOULD
C BE USED.
C
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
REAL RSTAK(5000)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (DSTAK(1),RSTAK(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(3),LUSED)
C
C SET UP MACHINE CONSTANTS, USING PORT SUBPROGRAMS
C
IOUTD = I1MACH(2)
C
C THIS PROGRAM USES RANDOM DISK STORAGE, WITH A LOGICAL RECORD SIZE
C OF M1*M2. THE TOTAL FILE IS TREATED AS IF IT WERE AN ARRAY OF
C SIZE N1=K1*M1 BY N2=K2*M2.
C N1 MUST BE AN INTEGER NUMBER OF RECORDS, THUS K1 MUST BE A MULTIPLE
C OF M2. M1 MUST ALSO BE A MULTIPLE OF 2.
C
C USES SUBROUTINES FFT, FFTMX, REALT, ISTKGT AND ISTKRL FROM SEC.1.4
C
C USES DONALD FRASER'S SUBROUTINES MFREAD AND MFWRIT, WITH THE MASS
C STORAGE UNIT NUMBER LUN COMMUNICATED THROUGH LABELED COMMON FFTCOM.
C
LUN = 2
C
C THE FOLLOWING INPUT PARAMETERS SPECIFY A 240 BY 225 REAL FFT
C
J1 = 1
K2 = 15
I1 = 8
M2 = 15
C
K1 = J1*M2
KT = K1*K2
M1 = 2*I1
MT = M1*M2
N1 = K1*M1
N2 = K2*M2
IB = MAX0((N1+2)*M2,N2*M1)
C
WRITE (IOUTD,9999)
9999 FORMAT (///12X, 39HTEST OF TWO-DIMENSIONAL MIXED RADIX FFT)
WRITE (IOUTD,9998) MT, IB, N1, N2
9998 FORMAT (/14H RECORD SIZE =, I4, 6X, 13HBUFFER SIZE =, I5, I12,
* 3H BY, I5, 10H TRANSFORM)
C
C SET UP WORKING STORAGE FOR INPUT/OUTPUT
C
IA = ISTKGT(IB,3)
C
C WRITE A TEST FILE
C
IB = IA + MT - 1
L = 0
DO 20 J=1,KT
DO 10 K=IA,IB
RSTAK(K) = L
L = L + 1
10 CONTINUE
CALL MFWRIT(RSTAK(IA), MT, J)
20 CONTINUE
C
C TRANSFORM FROM TIME TO FREQUENCY
C
CALL FFT2T(RSTAK(IA), J1, K2, I1, M2, 1)
C
C HERE WE ARE IN THE FREQUENCY DOMAIN --- ANY FILTERING TO DO?
C
C TRANSFORM BACK FROM FREQUENCY TO TIME DOMAIN
C
CALL FFT2I(RSTAK(IA), J1, K2, I1, M2, 1)
C
C NOW BACK IN THE TIME DOMAIN...WILL CHECK TO SEE IF WE RETURNED SAFELY
C
L = 0
SSA = 0.0
DO 40 J=1,KT
CALL MFREAD(RSTAK(IA), MT, J)
DO 30 K=IA,IB
SSA = (RSTAK(K)-FLOAT(L))**2 + SSA
L = L + 1
30 CONTINUE
40 CONTINUE
C
C COMPARE ERROR SUM OF SQUARES WITH SUM OF SQUARES ABOUT MEAN FOR DATA
C
SSB = FLOAT(L)*FLOAT(L+1)*FLOAT(L-1)/12.0
SSA = SQRT(SSA/SSB)
WRITE (IOUTD,9997) SSA
9997 FORMAT (/44H RMS ERROR NORM FOR TRANSFORM-INVERSE PAIR =, E10.3)
C
WRITE (IOUTD,9996) LUSED
9996 FORMAT (/21H MAXIMUM STACK SIZE =, I6)
C
C DE-ALLOCATE WORKING STORAGE...AS A FINAL CHECK ON ARRAY BOUNDS
C
IF (LOUT.NE.0) CALL ISTKRL(LOUT)
C
STOP
END
C
C-----------------------------------------------------------------------
C BLOCK DATA --- INITIALIZES LABELED COMMON
C-----------------------------------------------------------------------
C
BLOCK DATA
C
COMMON /CSTAK/ DSTAK(2500)
C
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
INTEGER ISIZE(5)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(2),LNOW)
EQUIVALENCE (ISTAK(3),LUSED)
EQUIVALENCE (ISTAK(4),LMAX)
EQUIVALENCE (ISTAK(5),LBOOK)
EQUIVALENCE (ISTAK(6),ISIZE(1))
C
DATA ISIZE(1), ISIZE(2), ISIZE(3), ISIZE(4), ISIZE(5) /1,1,1,2,2/
DATA LOUT, LNOW, LUSED, LMAX, LBOOK /0,10,10,5000,10/
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFT2T
C COMPUTES TWO-DIMENSIONAL FOURIER TRANSFORM FOR REAL OR COMPLEX DATA,
C IN TWO PASSES OF A MASS STORAGE UNIT LUN
C-----------------------------------------------------------------------
C
SUBROUTINE FFT2T(A, JJ, KK, LL, MM, IR)
C
C THE PARAMETER IR CONTROLS THE CALLING OF REALT. IF IR .NE. 0,
C COMPUTES AN N1 BY N2 FOURIER TRANSFORM OF REAL DATA STORED IN
C RECORDS OF SIZE M1*M2 ON A MASS STORAGE FILE LUN. THE FILE
C NUMBER LUN MUST BE STORED BY THE USER IN COMMON FFTCOM.
C
C ARRAY A IS THE INPUT/OUTPUT ARRAY FOR THE MASS STORAGE FILE LUN,
C DIMENSIONED MAX0(N1*M2,N2*M1) FOR A COMPLEX TRANSFORM
C OR MAX0((N1+2)*M2,N2*M1) FOR A REAL TRANSFORM.
C
C THE PARAMETERS M1, M2, N1 AND N2 ARE COMMUNICATED AS FOLLOWS:
C
C M1=IABS(LL)*2 I.E., M1 MUST BE EVEN
C M2=IABS(MM)
C N1=IABS(JJ)*M1*M2 I.E., N1 MUST BE A MULTIPLE OF M1*M2
C N2=IABS(KK)*M2
C
C ON ENTRY, THE REAL INPUT VALUES ARE ASSUMED TO BE STORED IN RECORDS
C OF LENGTH M1*M2, ARRANGED IN SEQUENCE WITH N1 VALUES FOR THE
C FIRST COLUMN, N1 VALUES FOR THE SECOND COLUMN, ETC. THE COLUMNS
C ARE TRANSFORMED, THEN THE DATA ARE PERMUTED TO M1 BY M2 SEGMENTS,
C SO THAT THE FILE CAN BE READ BY ROWS TO TRANSFORM THE SECOND
C DIMENSION.
C
C ON EXIT, THE FOURIER COEFFICIENTS ARE ON MASS STORAGE FILE LUN,
C WHICH IS LEFT PARTITIONED IN M1 BY M2 SUB-MATRICES...SO THAT
C THE RESULTS CAN BE READ EITHER BY ROW OR COLUMN. THE STRUCTURE OF
C THE FIRST SEGMENT, FOR EXAMPLE, IS AS FOLLOWS:
C
C 0,0 0,1 0,2 ... 0,M2-1
C 1,0 1,1 1,2 ... 1,M2-1
C 2,0 2,1 2,2 ... 2,M2-1
C . . . .
C . . . .
C . . . .
C M1-1,0 M1-1,1 M1-1,2 ... M1-1,M2-1
C
C THE COSINE AND SINE COEFFICIENTS ALTERNATE DOWN COLUMNS, EXCEPT
C THAT THE FOLDING FREQUENCY COSINE COEFFICIENT VALUES (ALTERNATING
C WITH ZERO SINE COEFFICIENT VALUES) ARE STORED AT THE END OF THE
C MAIN FILE IN RECORDS OF LENGTH M1*M2, WITH A POSSIBLE SHORT
C RECORD TO BRING THE TOTAL LENGTH OF THE ADDED RECORDS TO 2*N2.
C
C TO DO AN N1/2 BY N2 COMPLEX FOURIER TRANSFORM, WHERE THE INPUT
C VALUES ARE ARRANGED WITH REAL AND IMAGINARY COMPONENTS ALTERN-
C ATING, CALL THIS SUBROUTINE WITH IR=0.
C
COMMON /FFTCOM/ LUN
DIMENSION A(1)
C
COMMON /CSTAK/ DSTAK(2500)
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
REAL RSTAK(5000)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (DSTAK(1),RSTAK(1))
C
IF (JJ*KK*LL*MM.NE.0) GO TO 10
IERR = I1MACH(4)
WRITE (IERR,9999) JJ, KK, LL, MM
9999 FORMAT (33H ERROR - ZERO IN FFT2T PARAMETERS, 4I9)
RETURN
C
10 J1 = IABS(JJ)
K2 = IABS(KK)
M1H = IABS(LL)
M1 = 2*M1H
M2 = IABS(MM)
M2C = 2*M2
K1 = J1*M2
KT = K1*K2
MT = M1*M2
N1H = K1*M1H
N1 = K1*M1
N2 = K2*M2
L1 = K1*MT
L2 = K2*MT
JM = J1*MT
IF (IR.EQ.0) GO TO 20
C
C SET UP WORKING STORAGE FOR FOLDING FREQUENCY COEFFICIENTS
C
IC = ISTKGT(2*N2,3)
JC = IC
C
20 DO 60 J=1,KT,K1
I = J
DO 30 L=1,L1,MT
CALL MFREAD(A(L), MT, I)
I = I + 1
30 CONTINUE
CALL FFT(A, A(2), M2, N1H, 1, -2)
IF (IR.EQ.0) GO TO 40
CALL REALT(A, A(2), M2, N1H, 1, -2)
CALL XFR(A(L1+1), RSTAK(JC), M2C)
JC = JC + M2C
C
C THE FOLLOWING SECTION PARTITIONS THE N1 BY N2 MASS STORAGE INTO
C M1 BY M2 SUB-MATRICES.
C
40 CALL TRNSP(A, N1, M1, M2)
I = J
L = 1
50 CALL MFWRIT(A(L), MT, I)
I = I + 1
L = L + JM
IF (L.LT.L1) GO TO 50
L = L - L1 + MT
IF (L.LT.JM) GO TO 50
C
60 CONTINUE
C
DO 90 J=1,K1
I = J
DO 70 L=1,L2,MT
CALL MFREAD(A(L), MT, I)
I = I + K1
70 CONTINUE
CALL FFT(A, A(2), 1, N2, M1H, -2)
I = J
DO 80 L=1,L2,MT
CALL MFWRIT(A(L), MT, I)
I = I + K1
80 CONTINUE
90 CONTINUE
C
IF (IR.EQ.0) RETURN
CALL FFT(RSTAK(IC), RSTAK(IC+1), 1, N2, 1, -2)
J = KT + 1
JC = IC
K = 2*N2
100 CALL MFWRIT(RSTAK(JC), MIN0(K,MT), J)
J = J + 1
JC = JC + MT
K = K - MT
IF (K.GT.0) GO TO 100
CALL ISTKRL(1)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFT2I
C COMPUTES TWO-DIMENSIONAL FOURIER TRANSFORM INVERSE FOR REAL OR
C COMPLEX DATA, IN TWO PASSES OF A MASS STORAGE UNIT LUN
C-----------------------------------------------------------------------
C
SUBROUTINE FFT2I(A, JJ, KK, LL, MM, IR)
C
C THE PARAMETER IR CONTROLS THE CALLING OF REALT. IF IR .NE. 0,
C COMPUTES AN INVERSE FOURIER TRANSFORM, USING AS INPUT AN N1 BY N2
C FILE OF FOURIER COEFFICIENTS ON MASS STORAGE FILE LUN. THE FILE
C NUMBER LUN MUST BE STORED BY THE USER IN COMMON FFTCOM.
C
C ARRAY A IS THE INPUT/OUTPUT ARRAY FOR THE MASS STORAGE FILE LUN,
C DIMENSIONED MAX0(N1*M2,N2*M1) FOR A COMPLEX TRANSFORM
C OR MAX0((N1+2)*M2,N2*M1) FOR A REAL TRANSFORM.
C
C THE PARAMETERS M1, M2, N1 AND N2 ARE COMMUNICATED AS FOLLOWS:
C
C M1=IABS(LL)*2 I.E., M1 MUST BE EVEN
C M2=IABS(MM)
C N1=IABS(JJ)*M1*M2 I.E., N1 MUST BE A MULTIPLE OF M1*M2
C N2=IABS(KK)*M2
C
C ON ENTRY, THE FOURIER COEFFICIENTS ARE ON MASS STORAGE FILE LUN,
C WHICH IS PARTITIONED IN M1 BY M2 SUB-MATRICES...SO THAT
C THE RESULTS CAN BE READ EITHER BY ROW OR COLUMN. THE STRUCTURE OF
C THE FIRST SEGMENT, FOR EXAMPLE, IS AS FOLLOWS:
C
C 0,0 0,1 0,2 ... 0,M2-1
C 1,0 1,1 1,2 ... 1,M2-1
C 2,0 2,1 2,2 ... 2,M2-1
C . . . .
C . . . .
C . . . .
C M1-1,0 M1-1,1 M1-1,2 ... M1-1,M2-1
C
C THE COSINE AND SINE COEFFICIENTS ALTERNATE DOWN COLUMNS, EXCEPT
C THAT THE FOLDING FREQUENCY COSINE COEFFICIENT VALUES (ALTERNATING
C WITH ZERO SINE COEFFICIENT VALUES) ARE STORED AT THE END OF THE
C MAIN FILE IN RECORDS OF LENGTH M1*M2, WITH A POSSIBLE SHORT
C RECORD TO BRING THE TOTAL LENGTH OF THE ADDED RECORDS TO 2*N2.
C
C ON EXIT, THE REAL RESULT VALUES ARE STORED ON THE FILE LUN IN RECORDS
C OF LENGTH M1*M2, ARRANGED IN SEQUENCE WITH N1 VALUES FOR THE
C FIRST COLUMN, N1 VALUES FOR THE SECOND COLUMN, ETC.
C
C THE FOLLOWING PAIR OF CALLS PRODUCES AN IDENTITY TRANSFORMATION:
C
C CALL FFT2D(A,JJ,KK,LL,MM,IR)
C CALL FFT2I(A,JJ,KK,LL,MM,IR)
C
C TO DO AN N1/2 BY N2 COMPLEX INVERSE FOURIER TRANSFORM, WHERE THE
C INPUT VALUES ARE ARRANGED WITH COSINE AND SINE COEFFICIENTS
C ALTERNATING, CALL THIS SUBROUTINE WITH IR=0.
C
COMMON /FFTCOM/ LUN
DIMENSION A(1)
C
COMMON /CSTAK/ DSTAK(2500)
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
REAL RSTAK(5000)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (DSTAK(1),RSTAK(1))
C
IF (JJ*KK*LL*MM.NE.0) GO TO 10
IERR = I1MACH(4)
WRITE (IERR,9999) JJ, KK, LL, MM
9999 FORMAT (33H ERROR - ZERO IN FFT2I PARAMETERS, 4I9)
RETURN
C
10 J1 = IABS(JJ)
K2 = IABS(KK)
M1H = IABS(LL)
M1 = 2*M1H
M2 = IABS(MM)
M2C = 2*M2
K1 = J1*M2
KT = K1*K2
MT = M1*M2
N1H = K1*M1H
N1 = K1*M1
N2 = K2*M2
L1 = K1*MT
L2 = K2*MT
KM = M2*MT
IF (IR.EQ.0) GO TO 30
C
C SET UP WORKING STORAGE FOR FOLDING FREQUENCY COEFFICIENTS, AND
C RETRIEVE THEM FROM THE END OF THE MAIN FILE.
C
IC = ISTKGT(2*N2,3)
J = KT + 1
JC = IC
K = 2*N2
20 CALL MFREAD(RSTAK(JC), MIN0(K,MT), J)
J = J + 1
JC = JC + MT
K = K - MT
IF (K.GT.0) GO TO 20
CALL FFT(RSTAK(IC), RSTAK(IC+1), 1, N2, 1, 2)
JC = IC
C
30 DO 60 J=1,K1
I = J
DO 40 L=1,L2,MT
CALL MFREAD(A(L), MT, I)
I = I + K1
40 CONTINUE
CALL FFT(A, A(2), 1, N2, M1H, 2)
I = J
DO 50 L=1,L2,MT
CALL MFWRIT(A(L), MT, I)
I = I + K1
50 CONTINUE
60 CONTINUE
C
DO 100 J=1,KT,K1
C
C THE FOLLOWING SECTION RESTORES THE ORIGINAL ORDER OF THE N1 BY N2
C MASS STORAGE, ELIMINATING THE M1 BY M2 SUB-MATRICES.
C
I = J
L = 1
70 CALL MFREAD(A(L), MT, I)
I = I + 1
L = L + KM
IF (L.LT.L1) GO TO 70
L = L - L1 + MT
IF (L.LT.KM) GO TO 70
CALL TRNSP(A, N1, M1, M2)
C
IF (IR.EQ.0) GO TO 80
CALL XFR(RSTAK(JC), A(L1+1), M2C)
JC = JC + M2C
CALL REALT(A, A(2), M2, N1H, 1, 2)
80 CALL FFT(A, A(2), M2, N1H, 1, 2)
I = J
DO 90 L=1,L1,MT
CALL MFWRIT(A(L), MT, I)
I = I + 1
90 CONTINUE
100 CONTINUE
C
IF (IR.NE.0) CALL ISTKRL(1)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: XFR
C MOVES A(J) TO B(J) FOR J=1,2,...,N
C-----------------------------------------------------------------------
C
SUBROUTINE XFR(A, B, N)
C
DIMENSION A(1), B(1)
DO 10 J=1,N
B(J) = A(J)
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRNSP
C SELF-INVERSE PERMUTATION OF N1*M2 DATA VALUES IN ARRAY A
C-----------------------------------------------------------------------
C
SUBROUTINE TRNSP(A, N1, M1, M2)
C
C SUPPOSE ARRAY A IS THE INPUT/OUTPUT BUFFER FOR A MASS STORE,
C WHERE THE RECORD SIZE IS M1*M2. SUPPOSE THIS MASS STORE HOLDS
C A RECTANGULAR ARRAY, STORED BY COLUMNS OF LENGTH N1, WHERE N1
C IS A MULTIPLE OF M1*M2 AND THE NUMBER OF ROWS IS A MULTIPLE
C OF M2. THIS SUBROUTINE REARRANGES THE BUFFER SO THAT ALL
C ELEMENTS BELONGING TO A ROW ARE IN THE SAME RECORD SEGMENT.
C WHILE THE SEGMENTS ARE NOT IN NORMAL ROW ORDER, THIS ORDER
C CAN BE OBTAINED BY REORDERING THE OUTPUT, AS IN 'FFT2T'...THEN
C BEFORE CALLING THIS SUBROUTINE A SECOND TIME TO RESTORE THE
C ORIGINAL ORDER, THE SEGMENTS MUST BE REORDERED ON INPUT, AS
C IN 'FFT2I'.
C
DIMENSION A(1)
NB = N1*M2
MT = M1*M2
MM = MT - M1
J = 0
K = 0
10 J = J + M1
K = K + N1
IF (J.GE.K) GO TO 10
20 CALL EXCH(A(J), A(K), M1)
J = J + M1
K = K + N1
IF (K.LT.NB) GO TO 20
K = K - NB + MT
IF (K.LT.N1) GO TO 10
K = K - N1 + M1
IF (K.NE.MM) GO TO 10
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: EXCH
C EXCHANGES A(J) AND B(J) FOR J=2,3,...,N+1
C-----------------------------------------------------------------------
C
SUBROUTINE EXCH(A, B, N)
C
DIMENSION A(1), B(1)
J = 1
10 J = J + 1
T = A(J)
A(J) = B(J)
B(J) = T
IF (J.LE.N) GO TO 10
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MFREAD
C READS RECORD JB FROM MASS STORE TO BUFA, NB REAL VALUES
C SOURCE: DONALD FRASER, OPTIMIZED MASS STORAGE FFT PROGRAM, APPENDIX C
C-----------------------------------------------------------------------
C
SUBROUTINE MFREAD(BUFA, NB, JB)
REAL BUFA(NB)
COMMON /FFTCOM/ LUN
C
C READ(LUN'JB) BUFA
C
COMMON /MASS/ RMAS(54450)
J = (JB-1)*NB + 1
CALL XFR(RMAS(J), BUFA, NB)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: MFWRIT
C WRITES RECORD JB FROM BUFA TO MASS STORE, NB REAL VALUES
C SOURCE: DONALD FRASER, OPTIMIZED MASS STORAGE FFT PROGRAM, APPENDIX C
C-----------------------------------------------------------------------
C
SUBROUTINE MFWRIT(BUFA, NB, JB)
REAL BUFA(NB)
COMMON /FFTCOM/ LUN
C
C WRITE(LUN'JB) BUFA
C
COMMON /MASS/ RMAS(54450)
J = (JB-1)*NB + 1
CALL XFR(BUFA, RMAS(J), NB)
RETURN
END
Return to Main Software Page