c2fft.f

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

      SUBROUTINE C2FFT(C,ID,NM,NN,WM,WN,ISIG,IORD,IWORK,IERR)
****PURPOSE:
*       THIS ROUTINE PERFORMS A 2-DIMENSIONAL COMPLEX FOURIER TRANSFORM,
*       OF ORDER NM*NN.
****USAGE:
*       THE USER IS EXPECTED TO PROVIDE THE DATA IN A 2-DIMENSIONAL
*       COMPLEX ARRAY C, DIMENSIONED IN THE CALLING PROGRAM C(ID,NN);
*       ID CAN BE DIFFERENT FROM NM, AND IT IS RECOMMENDED THAT IT IS
*       CHOSEN EQUAL TO NM+1 IF NM IS EVEN OR TO NM IF NM IS ODD.
*       THE ELEMENTS C(K,*),NM<K<=ID MUST BE ZEROED. THE ROUTINE IS
*       INTENDED FOR REPEATED USAGE, THUS SEPARATE SET-UP AND
*       OPERATING CALLS ARE AVAILABLE : THE USER SHOULD ALWAYS PERFORM
*       A SET-UP CALL ( ISIG=0 ) PASSING THE CHOSEN PARAMETERS, BEFORE
*       PERFORMING THE ACTUAL TRANSFORM ( ISIG= +1 OR -1 ); THE USER CAN
*       CHOOSE WHETHER TO OBTAIN THE RESULTS OF THE DIRECT TRANSFORM
*       IN NATURAL ORDER (ISIG=-1,IORD=1) OR LEAVE THEM IN THE
*      'BIT-REVERSED' ORDER( ISIG=-1,IORD=0); THIS CHOICE SAVES
*       SOME COMPUTER TIME, AND IT IS RECOMMENDED IN CASES DISCUSSED
*       IN THE LONG WRITE-UP. ANALOGOUSLY,THE INVERSE TRANSFORM ACCEPTS
*       INPUT ( PLEASE NOTE| ) DATA IN NATURAL ORDER ( ISIG=1,IORD=1),
*       OR DATA ALREADY SUBJECTED TO A BIT-REVERSAL PERMUTATION( ISIG=1
*       IORD=0 ).
****PARAMETERS :
*       INPUT :
*       C : ARRAY TO BE TRANSFORMED; DECLARED COMPLEX C(ID,NN) IN THE
*           CALLING PROGRAM;
*       ID : FIRST DIMENSION OF C IN THE CALLING PROGRAM
*       ISIG : OPTION FLAG : ISIG=0 : SET-UP RUN, C NOT USED
*                            ISIG=-1: DIRECT TRANSFORM
*                            ISIG=+1: INVERSE TRANSFORM
*       WM,WN : INTEGER ARRAYS , USED TO HOST TABLES FOR THE TRANSFORM;
*               DIMENSIONED IN THE CALLING PROGRAM AT LEAST (4*NM+14)
*               AND (4*NN+14) RESPECTIVELY; IF ISIG.NE.0 ,THEY ARE
*               ASSUMED TO HAVE BEEN SET BY A PREVIOUS CALL WITH ISIG=0
*               (ALL OTHER ARGUMENTS UNCHANGED), AND NEVER HAVE BEEN
*               MODIFIED SINCE THEN;
*               IF NM=NN, WM AND WN DO NOT NEED TO BE DISTINCT;
*       NM : ORDER OF THE TRANSFORM ALONG THE COLUMNS OF C
*       NN : ORDER OF THE TRANSFORM ALONG THE ROWS OF C
*       IORD : OPTION FLAG : =1 : OUTPUT IN NATURAL ORDER (ISIG=-1)
*                                 INPUT IN NATURAL ORDER  (ISIG=+1)
*                            =0 : OUTPUT IN BIT-REVERSED ORDER(ISIG=-1)
*                                 INPUT IN BIT-REVERSED ORDER(ISIG=+1)
*       IWORK : INTEGER ARRAY, USED AS WORK AREA FOR REORDERING IF IORD=
*                1, IT MUST BE AT LEAST MAX(NM,NN) WORDS LONG.
*                IT IS UNUSED IF IORD=0.
*
****OUTPUT :
*       C : TRANSFORMED ARRAY
*       WM, WN : ONLY IF ISIG=0, WM AND WN ARE FILLED WITH THE
*               APPROPRIATE  TABLES (TWIDDLE FACTORS).
*       IWORK : UNDEFINED
*       IERR  : ERROR CODE : =0 : SUCCESSFUL
*                            =1 : DATA DIMENSIONS ARE NOT CORRECT
*                            =2 : PRIME FACTORS DIFFERENT FROM 2,3,5
*                                 ARE PRESENT IN DATA DIMENSIONS
*                            =3 : TABLES NOT CORRECTLY INITIALIZED
*
      COMPLEX C(*)
      INTEGER WM(-14:*),WN(-14:*)
      INTEGER IWORK(*)
      INTEGER IDERR,FACERR,TBERR
      PARAMETER (IDERR=1,FACERR=2,TBERR=3)
      IF(ID.LT.NM)THEN
        IERR=IDERR
        RETURN
      ENDIF
      IERR=0
      IF(ISIG.EQ.0)THEN
        CALL MFFTP(NM,WM,0,IERR)
        IF(IERR.NE.0)RETURN
        IF(NN.NE.NM)THEN
          CALL MFFTP(NN,WN,0,IERR)
          IF(IERR.NE.0)RETURN
        ELSE
          CALL MFFTZ0(WM,1,4*NM+14,WN,1)
        ENDIF
        RETURN
 
      ELSE IF(ISIG.GT.0)THEN
 
        IF(IORD.NE.0)THEN
         CALL MFFTOV(C,1,ID,NM,NN,WM(3*NM),IWORK)
         CALL MFFTOV(C,ID,1,NN,NM,WN(3*NN),IWORK)
        ENDIF
 
        CALL MFFTIV(C,1,ID,NM,NN,WM,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTIV(C,ID,1,NN,NM,WN,IERR)
        IF(IERR.NE.0)RETURN
 
      ELSE
        CALL MFFTDV(C,1,ID,NM,NN,WM,IERR)
        IF(IERR.NE.0)RETURN
        CALL MFFTDV(C,ID,1,NN,NM,WN,IERR)
        IF(IERR.NE.0)RETURN
 
        IF(IORD.NE.0)THEN
         CALL MFFTOV(C,1,ID,NM,NN,WM(2*NM),IWORK)
         CALL MFFTOV(C,ID,1,NN,NM,WN(2*NN),IWORK)
        ENDIF
      ENDIF
 
      END