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