mfftp2.f

来自「FFT源代码汇集（各种FFT源代码）」· F 代码 · 共 62 行

      SUBROUTINE MFFTP2(INDX,I2,IW,NM)
*
*     THIS ROUTINE COMPUTES TWO INDEX TABLES FOR THE
*     PERMUTATION DUE TO REPRESENTATION INVERSION ("BIT REVERSAL")
*     THE INDEX TABLES ARE STORED ONE AFTER THE OTHER, AND ARE
*     RECIPROCAL.
*
*     WARNING: BY "BIT REVERSAL" WE MEAN THE SHUFFLING OF INDEXES
*              AS REQUIRED BY IN-PLACE FFT ALGORITHMS, REGARDLESS
*              OF WHAT IS THEIR RADIX (I.E. 2,3,5).
*              THE SHUFFLING IS EFFECTIVELY AN 'INTEGER BIT REVERSAL'
*              ONLY IN CASE OF RADIX-2 ALGORITHMS.
*
      INTEGER INDX(0:*),I2(0:*),IW(14)
      PARAMETER (MAXFAC=3)
      INTEGER FACTOR(MAXFAC)
      DATA (FACTOR(I),I=1,MAXFAC)/2,3,5/
      DO 1 I=0,NM-1
        INDX(I)=I
1     CONTINUE
      IFAC=IW(14)
      LX=1
      MX=NM
      DO 30 IFAC=1,IFAC
        DO 10 I=1,IW(IFAC)-1,2
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(INDX,I2,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(I2,INDX,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
10      CONTINUE
 
*...  IF W(IFAC)ODD,THEN
        IF (I.EQ.IW(IFAC))THEN
          MX=MX/FACTOR(IFAC)
          CALL MFFTP3(INDX,I2,MX,FACTOR(IFAC),LX)
          LX=LX*FACTOR(IFAC)
          DO 20 I=0,NM-1
            INDX(I)=I2(I)
20        CONTINUE
        ENDIF
30    CONTINUE
*...     INVERSE PERMUTATION
         DO 40 I=0,NM-1
           I2(I)=INDX(I)
           INDX(I+NM)=I
40       CONTINUE
         DO 59 I=1,NM-3
51         IF(I2(I).NE.I)THEN
             IDEST=I2(I)
             IT=INDX(I+NM)
             INDX(I+NM)=INDX(IDEST+NM)
             INDX(IDEST+NM)=IT
             I2(I)=I2(IDEST)
             I2(IDEST)=IDEST
           GOTO 51
           ENDIF
59       CONTINUE
 
      END