mfftds.f

FFT源代码汇集（各种FFT源代码）

      SUBROUTINE MFFTDS(C,IMSX,IVSX,IESX,NMX,NVX,NEX,TABLES,IERR)
 
*   PURPOSE:
*       THE SAME AS MFFTDM.    IT IS A VARIANT OF MFFTDM
*       AIMED AT OPTIMAL PERFORMANCE ON MATRICES HAVING THE
*       FIRST DIMENSION SMALLER THAN 64. IT REQUIRES THAT MFFTP HAS
*       BEEN CALLED WITH ID .NE. 0
*  IMPLEMENTATION:
*       THE TRANSFORMATION IS IMPLEMENTED THROUGH REPEATED CALLS TO THE
*       "BUTTERFLY" TRANSFORMATION MFFT?8; PARAMETERS OF THE "BUTTERFLY"
*       ARE COMMUNICATED THROUGH THE COMMON BLOCK MFFTPA.
*  ARGUMENTS:
*    INPUT  :
*       C : ARRAY TO BE TRANSFORMED.
*       IMSX,IVSX,IESX,NMX,NVX,NEX: THESE ARGUMENTS DEFINE THE STRUCTURE
*           OF C ACCORDING TO THE DEFINITIONS ABOVE. THEY ARE UNCHANGED
*           ON OUTPUT
*       TABLES : ARRAY PREPARED BY MFFTP. IT  IS NOT CHANGED ON OUTPUT.
*                IT SHOULD BE DECLARED INTEGER TABLES(4*NM+14);
*                IT MUST BE INITIALIZED BY MFFTP BEFORE USAGE.
*    OUTPUT:
*       C : TRANSFORM OF THE ORIGINAL ARRAY; "BIT REVERSED" ORDER
*       IERR : ERROR CODE : =0 : SUCCESSFUL
*                         : =3 :  'TABLES' NOT CORRECTLY INITIALIZED
      COMPLEX C(*)
      INTEGER TABLES(-14:*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      COMMON /MFFTPA/  IMS,IVS,IES,NM,NV,NE,MX,LX,MLIM,MSTEP,LLIM,LSTEP,
     $ NUSTEP,IVLIM,ILIM,MD2LIM,LD2LIM
      SAVE IBASE
*  LOADING THE COMMON BLOCK : CONSTANTS
      IMS=IMSX
      IVS=IVSX
      NM=NMX
      NV=NVX
      NE=NEX
      IVLIM=(NV-1)*IVS
      MSTEP=IMS
*  LOADING THE COMMON BLOCK : ITERATION-DEPENDENT QUANTITIES: INITIALIZA
      LX=1
      MX=NM
 
      IBASE=2*NM
*  SELECT THE HIGHEST FACTOR OF NM
      IFAC=TABLES(-1)
      GOTO(200,300,500)IFAC
      IERR=TBERR
      RETURN
*...  RADIX 5 LOOP
500   CONTINUE
      DO 510 IM=1,TABLES(-12)
        MX=MX/5
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        LSTEP=NUSTEP*5
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTC8(C,TABLES(2*IBASE))
        LX=LX*5
        IBASE=IBASE+NUSTEP*4
510   CONTINUE
*..  RADIX 3 LOOP
300   CONTINUE
      DO 310 IM=1,TABLES(-13)
        MX=MX/3
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        LSTEP=NUSTEP*3
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTB8(C,TABLES(2*IBASE))
        LX=LX*3
        IBASE=IBASE+NUSTEP*2
310   CONTINUE
*..  RADIX 2 LOOP
200   CONTINUE
      DO 210 IM=1,TABLES(-14)
        MX=MX/2
        NUSTEP=MX*MSTEP
        ILIM=NUSTEP-1
        MD2LIM=NUSTEP/2-MSTEP
        LSTEP=NUSTEP*2
        LLIM=NM*MSTEP-LSTEP
        CALL MFFTA8(C,TABLES(2*IBASE))
        LX=LX+LX
        IBASE=IBASE+NUSTEP
210   CONTINUE
 
      END