! ***************** ! * NFFT1.FOR * ! * PURSER 1994 * ! ***************** ! BLOCKDATA DAT235 COMMON/FFT235/ LN2,LN3,LN5 DATA LN2/0/,LN3/0/,LN5/0/ END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINES CFFT, DFFT C ! C ! Fourier analyze (CFFT) or synthesize (DFFT) a line of complex data C ! RB <-> Real part of data and transform C ! QB <-> Imaginary part of data and transform C ! N --> Number of complex data of series (product of 2's,3's,5's) C ! PERIOD --> Period of one cycle C ! W <-> Coefficient array of N reals (re)computed only when C ! JUMBLE(0)=0 C ! JUMBLE <-> Coefficient array of N integers (re)computed only when C ! JUMBLE(0) set to zero on input to this routine C ! C ! For related routines, see subroutine RUMBLE (in which the C ! coefficient arrays are initialized), XCFFT, YCFFT (matrix versions C ! of CFFT), XDFFT, YDFFT (matrix versions of DFFT), RFFT (Fourier C ! analysis of real data), HFFT (Fourier synthesis of real data), and C ! matrix counterparts, XRFFT, YRFFT, XHFFT, YHFFT. C !------------------------------------------------------------------------------C SUBROUTINE CFFT(RB,QB,N,PERIOD,W,JUMBLE) COMMON/FFT235/ LN2,LN3,LN5 DIMENSION RB(0:*),QB(0:*),W(0:*),JUMBLE(0:*) RFAC=PERIOD/N GOTO 300 ENTRY DFFT(RB,QB,N,PERIOD,W,JUMBLE) RFAC=1./PERIOD ! FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS: DO J=N/2+1,N-1 I=N-J T =RB(J) RB(J)=RB(I) RB(I)=T T =QB(J) QB(J)=QB(I) QB(I)=T ENDDO 300 IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N) IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W) NM=N-1 NH=N/2 ! SCALE AND PERMUTE THE DATA: RB(0)=RB(0)*RFAC QB(0)=QB(0)*RFAC DO I=1,NM J=JUMBLE(I) IF(J.GT.I)THEN T1=RB(I) RB(I)=RB(J)*RFAC RB(J)=T1 T1=QB(I) QB(I)=QB(J)*RFAC QB(J)=T1 ELSE RB(I)=RB(I)*RFAC QB(I)=QB(I)*RFAC ENDIF ENDDO ! TRANSFORM THE DATA: MA=1 MB=N ! RADIX 4 LS=0 DO L=2,LN2,2 MB=MB/4 LS=L MA4=MA*4 MB2=MB*2 DO J=0,MA-1 JMB=J*MB JMB2=J*MB2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 DO I=0,NM,MA4 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA R0=RB(K0) R1=RB(K1) R2=RB(K2) R3=RB(K3) Q0=QB(K0) Q1=QB(K1) Q2=QB(K2) Q3=QB(K3) T1=R3*RF3-Q3*QF3 ! Q13 Q3=R3*QF3+Q3*RF3 ! R13 R3=R2*RF1-Q2*QF1 ! R12 R2=R2*QF1+Q2*RF1 ! Q12 Q2=Q3-R2 ! R23 R2=Q3+R2 ! Q22 Q3=R3+T1 ! R22 T1=R3-T1 ! Q23 R3=R1*RF2-Q1*QF2 ! R11 Q1=R1*QF2+Q1*RF2 ! Q11 R1=R0-R3 ! R21 R0=R0+R3 ! R20 RB(K3)=R1-Q2 ! R3 RB(K1)=R1+Q2 ! R1 Q2=Q0+Q1 ! Q20 Q1=Q0-Q1 ! Q21 QB(K0)=Q2+R2 ! Q0 QB(K2)=Q2-R2 ! Q2 RB(K2)=R0-Q3 ! R2 RB(K0)=R0+Q3 ! R0 QB(K3)=Q1-T1 ! Q3 QB(K1)=Q1+T1 ! Q1 ENDDO ENDDO MA=MA4 ENDDO IF(LS.NE.LN2)THEN ! RADIX 2 MB=MB/2 MA2=MA*2 DO J=0,MA-1 JMB=J*MB DO I=0,NM,MA2 K0=J+I K1=K0+MA RF1=W(JMB) QF1=W(NH+JMB) R0=RB(K0) Q0=QB(K0) R1=RB(K1) Q1=QB(K1) T1=R1*QF1+Q1*RF1 ! Q11 Q1=R1*RF1-Q1*QF1 ! R11 RB(K1)=R0-Q1 ! R1 RB(K0)=R0+Q1 ! R0 QB(K1)=Q0-T1 ! Q1 QB(K0)=Q0+T1 ! Q0 ENDDO ENDDO MA=MA2 ENDIF ! RADIX 3 REP=-.5 REC=1.5 QEP=.5*SQRT(3.) DO L=1,LN3 MB=MB/3 MA3=MA*3 DO J=0,MA-1 JMB=J*MB RF1=W(JMB) QF1=W(NH+JMB) RF2=RF1*RF1-QF1*QF1 QF2=2*RF1*QF1 DO I=0,NM,MA3 K0=I+J K1=K0+MA K2=K1+MA R1=RB(K1) Q1=QB(K1) R2=RB(K2) Q2=QB(K2) T1=R2*QF2+Q2*RF2 ! R12 Q2=R2*RF2-Q2*QF2 ! Q12 R2=R1*QF1+Q1*RF1 ! Q11 R1=R1*RF1-Q1*QF1 ! R11 Q1=R2+T1 ! Q21 R2=(R2-T1)*QEP ! R22 T1=R1+Q2 ! R21 R1=(R1-Q2)*QEP ! Q22 RB(K0)=RB(K0)+T1 ! R0 QB(K0)=QB(K0)+Q1 ! Q0 T1=RB(K0)-T1*REC ! R21 Q1=QB(K0)-Q1*REC ! Q21 QB(K2)=Q1-R1 ! Q2 QB(K1)=Q1+R1 ! Q1 RB(K1)=T1-R2 ! R1 RB(K2)=T1+R2 ! R2 ENDDO ENDDO MA=MA3 ENDDO IF(LN5.GT.0)THEN ! RADIX 5 NZE=N/5 RZE=W(NZE) QZE=W(NH+NZE) RZC=1.-RZE RET=RZE*RZE-QZE*QZE QET=2*RZE*QZE REC=1.-RET DO L=1,LN5 MB=MB/5 MA5=MA*5 DO J=0,MA-1 JMB=J*MB JMB2=JMB*2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 RF4=RF2*RF2-QF2*QF2 QF4=2*RF2*QF2 DO I=0,NM,MA5 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA K4=K3+MA R1=RB(K1) R2=RB(K2) R3=RB(K3) R4=RB(K4) Q1=QB(K1) Q2=QB(K2) Q3=QB(K3) Q4=QB(K4) T1=R1*QF1+Q1*RF1 ! Q11 R1=R1*RF1-Q1*QF1 ! R11 Q1=R4*RF4-Q4*QF4 ! Q14 R4=R4*QF4+Q4*RF4 ! R14 Q4=R1-Q1 ! Q24 R1=R1+Q1 ! R21 Q1=T1+R4 ! Q21 R4=T1-R4 ! R24 T1=R3*RF3-Q3*QF3 ! Q13 R3=R3*QF3+Q3*RF3 ! R13 Q3=R2*QF2+Q2*RF2 ! Q12 R2=R2*RF2-Q2*QF2 ! R12 Q2=Q3+R3 ! Q22 R3=Q3-R3 ! R23 Q3=R2-T1 ! Q23 R2=R2+T1 ! R22 RB(K0)=RB(K0)+R1+R2 ! R0 QB(K0)=QB(K0)+Q1+Q2 ! Q0 T1= R4*QZE+R3*QET ! R34 R3= R3*QZE-R4*QET ! R33 R4=RB(K0)-R2*RZC-R1*REC ! R32 R1=RB(K0)-R1*RZC-R2*REC ! R31 RB(K2)=R4+R3 ! R2 RB(K3)=R4-R3 ! R3 RB(K4)=R1+T1 ! R4 RB(K1)=R1-T1 ! R1 T1=QB(K0)-Q1*RZC-Q2*REC ! Q31 Q2=QB(K0)-Q2*RZC-Q1*REC ! Q32 Q1= Q3*QZE-Q4*QET ! Q33 Q4= Q4*QZE+Q3*QET ! Q34 QB(K3)=Q2+Q1 ! Q3 QB(K2)=Q2-Q1 ! Q2 QB(K1)=T1+Q4 ! Q1 QB(K4)=T1-Q4 ! Q4 ENDDO ENDDO MA=MA5 ENDDO ENDIF RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINE GET235 C ! C ! Analyze the facorization of N in terms of 2's, 3's, and 5's and verify C ! that these are the only prime factors. Set LN2, LN3, LN5 respectively to C ! the number of powers of 2, 3, 5. C ! N --> Number of data along the line of Fourier transformation C !------------------------------------------------------------------------------C SUBROUTINE GET235(N) COMMON/FFT235/ LN2,LN3,LN5 NKT=N NM=N-1 NH=N/2 LN2=-1 400 NK=NKT NKT=NK/2 LN2=LN2+1 IF(2*NKT.EQ.NK)GOTO 400 NKT=NK LN3=-1 401 NK=NKT NKT=NK/3 LN3=LN3+1 IF(3*NKT.EQ.NK)GOTO 401 NKT=NK LN5=-1 402 NK=NKT NKT=NK/5 LN5=LN5+1 IF(5*NKT.EQ.NK)GOTO 402 IF(NK.NE.1)STOP'N CANNOT BE FACTORED INTO ONLY 2s, 3s, 5s' RETURN END !------------------------------------------------------------------------------C ! R.J.Purser, National Meteorological Center, Washington D.C. 1994 C ! SUBROUTINE RUMBLE C ! C ! Initialize coefficient arrays JUMBLE and TUMBLE for use with the C ! fast-Fourier-transform routines when the number N of data has the prime- C ! factorization 2**LN2*3**LN3*5**LN5 C ! C ! JUMBLE <-- Permutation of N data encoded as a sequence of transpositions C ! TUMBLE <-- Trigonometric coefficients for use in FFT. The first half are C ! the cosines, the second half the sines, of uniformly increasing C ! relevant angles C !------------------------------------------------------------------------------C SUBROUTINE RUMBLE(JUMBLE,TUMBLE) PARAMETER(ML=20) COMMON/FFT235/ LN2,LN3,LN5 DIMENSION ND(ML),MD(ML),JUMBLE(0:*),TUMBLE(0:*) LN=LN2+LN3+LN5 N=2**LN2*3**LN3*5**LN5 NM=N-1 NH=N/2 PI2ON=8.*ATAN(1.)/N DO I=0,NH-1 ANG=PI2ON*I TUMBLE(I)=COS(ANG) TUMBLE(I+NH)=SIN(ANG) ENDDO ID=1 IS=0 DO I=1,LN5 IS=IS+1 MD(IS)=ID ID=ID*5 ENDDO DO I=1,LN3 IS=IS+1 MD(IS)=ID ID=ID*3 ENDDO DO I=1,LN2 IS=IS+1 MD(IS)=ID ID=ID*2 ENDDO ID=1 DO I=1,LN2 ND(IS)=ID ID=ID*2 IS=IS-1 ENDDO DO I=1,LN3 ND(IS)=ID ID=ID*3 IS=IS-1 ENDDO DO I=1,LN5 ND(IS)=ID ID=ID*5 IS=IS-1 ENDDO JUMBLE(0)=N DO I=1,NM IR=I J=0 DO L=1,LN KD=IR/ND(L) IR=IR-KD*ND(L) J=J+KD*MD(L) ENDDO JUMBLE(I)=J ENDDO DO I=1,NM J=JUMBLE(I) IF(J.LT.I)THEN 400 J=JUMBLE(J) IF(J.LT.I)GOTO 400 JUMBLE(I)=J ENDIF ENDDO RETURN END SUBROUTINE XCFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE) COMMON/FFT235/ LN2,LN3,LN5 DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*) RFAC=PERIOD/N ! GOTO 300 ! ! ENTRY XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE) ! RFAC=1./PERIOD ! !C FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS: ! DO J=N/2+1,N-1 ! I=N-J ! DO IY=1,NY ! T =RB(J,IY) ! RB(J,IY)=RB(I,IY) ! RB(I,IY)=T ! T =QB(J,IY) ! QB(J,IY)=QB(I,IY) ! QB(I,IY)=T ! ENDDO ! ENDDO ! 300 IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N) IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W) NM=N-1 NH=N/2 ! SCALE AND PERMUTE THE DATA: DO IY=1,NY RB(0,IY)=RB(0,IY)*RFAC QB(0,IY)=QB(0,IY)*RFAC ENDDO DO I=1,NM J=JUMBLE(I) IF(J.GT.I)THEN DO IY=1,NY T=RB(I,IY) RB(I,IY)=RB(J,IY)*RFAC RB(J,IY)=T T=QB(I,IY) QB(I,IY)=QB(J,IY)*RFAC QB(J,IY)=T ENDDO ELSE DO IY=1,NY RB(I,IY)=RB(I,IY)*RFAC QB(I,IY)=QB(I,IY)*RFAC ENDDO ENDIF ENDDO ! TRANSFORM THE DATA: MA=1 MB=N ! RADIX 4 LS=0 DO L=2,LN2,2 MB=MB/4 LS=L MA4=MA*4 MB2=MB*2 DO J=0,MA-1 JMB=J*MB JMB2=J*MB2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 DO I=0,NM,MA4 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA DO IY=1,NY T =RB(K3,IY)*RF3-QB(K3,IY)*QF3 ! Q13 QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3 ! R13 RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1 ! R12 RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1 ! Q12 QB(K2,IY)=QB(K3,IY)-RB(K2,IY) ! R23 RB(K2,IY)=QB(K3,IY)+RB(K2,IY) ! Q22 QB(K3,IY)=RB(K3,IY)+T ! R22 T =RB(K3,IY)-T ! Q23 RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2 ! R11 QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2 ! Q11 RB(K1,IY)=RB(K0,IY)-RB(K3,IY) ! R21 RB(K0,IY)=RB(K0,IY)+RB(K3,IY) ! R20 RB(K3,IY)=RB(K1,IY)-QB(K2,IY) ! R3 RB(K1,IY)=RB(K1,IY)+QB(K2,IY) ! R1 QB(K2,IY)=QB(K0,IY)+QB(K1,IY) ! Q20 QB(K1,IY)=QB(K0,IY)-QB(K1,IY) ! Q21 QB(K0,IY)=QB(K2,IY)+RB(K2,IY) ! Q0 QB(K2,IY)=QB(K2,IY)-RB(K2,IY) ! Q2 RB(K2,IY)=RB(K0,IY)-QB(K3,IY) ! R2 RB(K0,IY)=RB(K0,IY)+QB(K3,IY) ! R0 QB(K3,IY)=QB(K1,IY)-T ! Q3 QB(K1,IY)=QB(K1,IY)+T ! Q1 ENDDO ENDDO ENDDO MA=MA4 ENDDO IF(LS.NE.LN2)THEN ! RADIX 2 MB=MB/2 MA2=MA*2 DO J=0,MA-1 JMB=J*MB DO I=0,NM,MA2 K0=J+I K1=K0+MA RF1=W(JMB) QF1=W(NH+JMB) DO IY=1,NY T =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 RB(K1,IY)=RB(K0,IY)-QB(K1,IY) ! R1 RB(K0,IY)=RB(K0,IY)+QB(K1,IY) ! R0 QB(K1,IY)=QB(K0,IY)-T ! Q1 QB(K0,IY)=QB(K0,IY)+T ! Q0 ENDDO ENDDO ENDDO MA=MA2 ENDIF ! RADIX 3 REP=-.5 REC=1.5 QEP=.5*SQRT(3.) DO L=1,LN3 MB=MB/3 MA3=MA*3 DO J=0,MA-1 JMB=J*MB RF1=W(JMB) QF1=W(NH+JMB) RF2=RF1*RF1-QF1*QF1 QF2=2*RF1*QF1 DO I=0,NM,MA3 K0=I+J K1=K0+MA K2=K1+MA DO IY=1,NY T =RB(K2,IY)*QF2+QB(K2,IY)*RF2 ! R12 QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2 ! Q12 RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 QB(K1,IY)=RB(K2,IY)+T ! Q21 RB(K2,IY)=(RB(K2,IY)-T)*QEP ! R22 T =RB(K1,IY)+QB(K2,IY) ! R21 RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP ! Q22 RB(K0,IY)=RB(K0,IY)+T ! R0 QB(K0,IY)=QB(K0,IY)+QB(K1,IY) ! Q0 T =RB(K0,IY)-T*REC ! R21 QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC ! Q21 QB(K2,IY)=QB(K1,IY)-RB(K1,IY) ! Q2 QB(K1,IY)=QB(K1,IY)+RB(K1,IY) ! Q1 RB(K1,IY)=T -RB(K2,IY) ! R1 RB(K2,IY)=T +RB(K2,IY) ! R2 ENDDO ENDDO ENDDO MA=MA3 ENDDO IF(LN5.GT.0)THEN ! RADIX 5 NZE=N/5 RZE=W(NZE) QZE=W(NH+NZE) RZC=1.-RZE RET=RZE*RZE-QZE*QZE QET=2*RZE*QZE REC=1.-RET DO L=1,LN5 MB=MB/5 MA5=MA*5 DO J=0,MA-1 JMB=J*MB JMB2=JMB*2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 RF4=RF2*RF2-QF2*QF2 QF4=2*RF2*QF2 DO I=0,NM,MA5 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA K4=K3+MA DO IY=1,NY T =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4 ! Q14 RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4 ! R14 QB(K4,IY)=RB(K1,IY)-QB(K1,IY) ! Q24 RB(K1,IY)=RB(K1,IY)+QB(K1,IY) ! R21 QB(K1,IY)=T +RB(K4,IY) ! Q21 RB(K4,IY)=T -RB(K4,IY) ! R24 T =RB(K3,IY)*RF3-QB(K3,IY)*QF3 ! Q13 RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3 ! R13 QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2 ! Q12 RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2 ! R12 QB(K2,IY)=QB(K3,IY)+RB(K3,IY) ! Q22 RB(K3,IY)=QB(K3,IY)-RB(K3,IY) ! R23 QB(K3,IY)=RB(K2,IY)-T ! Q23 RB(K2,IY)=RB(K2,IY)+T ! R22 RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY) ! R0 QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY) ! Q0 T =RB(K4,IY)*QZE+RB(K3,IY)*QET ! R34 RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET ! R33 RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32 RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31 RB(K2,IY)=RB(K4,IY)+RB(K3,IY) ! R2 RB(K3,IY)=RB(K4,IY)-RB(K3,IY) ! R3 RB(K4,IY)=RB(K1,IY)+T ! R4 RB(K1,IY)=RB(K1,IY)-T ! R1 T =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31 QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32 QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET ! Q33 QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET ! Q34 QB(K3,IY)=QB(K2,IY)+QB(K1,IY) ! Q3 QB(K2,IY)=QB(K2,IY)-QB(K1,IY) ! Q2 QB(K1,IY)=T +QB(K4,IY) ! Q1 QB(K4,IY)=T -QB(K4,IY) ! Q4 ENDDO ENDDO ENDDO MA=MA5 ENDDO ENDIF RETURN END SUBROUTINE XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE) COMMON/FFT235/ LN2,LN3,LN5 DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*) RFAC=1./PERIOD ! FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS: DO J=N/2+1,N-1 I=N-J DO IY=1,NY T =RB(J,IY) RB(J,IY)=RB(I,IY) RB(I,IY)=T T =QB(J,IY) QB(J,IY)=QB(I,IY) QB(I,IY)=T ENDDO ENDDO 300 IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N) IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W) NM=N-1 NH=N/2 ! SCALE AND PERMUTE THE DATA: DO IY=1,NY RB(0,IY)=RB(0,IY)*RFAC QB(0,IY)=QB(0,IY)*RFAC ENDDO DO I=1,NM J=JUMBLE(I) IF(J.GT.I)THEN DO IY=1,NY T=RB(I,IY) RB(I,IY)=RB(J,IY)*RFAC RB(J,IY)=T T=QB(I,IY) QB(I,IY)=QB(J,IY)*RFAC QB(J,IY)=T ENDDO ELSE DO IY=1,NY RB(I,IY)=RB(I,IY)*RFAC QB(I,IY)=QB(I,IY)*RFAC ENDDO ENDIF ENDDO ! TRANSFORM THE DATA: MA=1 MB=N ! RADIX 4 LS=0 DO L=2,LN2,2 MB=MB/4 LS=L MA4=MA*4 MB2=MB*2 DO J=0,MA-1 JMB=J*MB JMB2=J*MB2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 DO I=0,NM,MA4 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA DO IY=1,NY T =RB(K3,IY)*RF3-QB(K3,IY)*QF3 ! Q13 QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3 ! R13 RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1 ! R12 RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1 ! Q12 QB(K2,IY)=QB(K3,IY)-RB(K2,IY) ! R23 RB(K2,IY)=QB(K3,IY)+RB(K2,IY) ! Q22 QB(K3,IY)=RB(K3,IY)+T ! R22 T =RB(K3,IY)-T ! Q23 RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2 ! R11 QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2 ! Q11 RB(K1,IY)=RB(K0,IY)-RB(K3,IY) ! R21 RB(K0,IY)=RB(K0,IY)+RB(K3,IY) ! R20 RB(K3,IY)=RB(K1,IY)-QB(K2,IY) ! R3 RB(K1,IY)=RB(K1,IY)+QB(K2,IY) ! R1 QB(K2,IY)=QB(K0,IY)+QB(K1,IY) ! Q20 QB(K1,IY)=QB(K0,IY)-QB(K1,IY) ! Q21 QB(K0,IY)=QB(K2,IY)+RB(K2,IY) ! Q0 QB(K2,IY)=QB(K2,IY)-RB(K2,IY) ! Q2 RB(K2,IY)=RB(K0,IY)-QB(K3,IY) ! R2 RB(K0,IY)=RB(K0,IY)+QB(K3,IY) ! R0 QB(K3,IY)=QB(K1,IY)-T ! Q3 QB(K1,IY)=QB(K1,IY)+T ! Q1 ENDDO ENDDO ENDDO MA=MA4 ENDDO IF(LS.NE.LN2)THEN ! RADIX 2 MB=MB/2 MA2=MA*2 DO J=0,MA-1 JMB=J*MB DO I=0,NM,MA2 K0=J+I K1=K0+MA RF1=W(JMB) QF1=W(NH+JMB) DO IY=1,NY T =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 RB(K1,IY)=RB(K0,IY)-QB(K1,IY) ! R1 RB(K0,IY)=RB(K0,IY)+QB(K1,IY) ! R0 QB(K1,IY)=QB(K0,IY)-T ! Q1 QB(K0,IY)=QB(K0,IY)+T ! Q0 ENDDO ENDDO ENDDO MA=MA2 ENDIF ! RADIX 3 REP=-.5 REC=1.5 QEP=.5*SQRT(3.) DO L=1,LN3 MB=MB/3 MA3=MA*3 DO J=0,MA-1 JMB=J*MB RF1=W(JMB) QF1=W(NH+JMB) RF2=RF1*RF1-QF1*QF1 QF2=2*RF1*QF1 DO I=0,NM,MA3 K0=I+J K1=K0+MA K2=K1+MA DO IY=1,NY T =RB(K2,IY)*QF2+QB(K2,IY)*RF2 ! R12 QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2 ! Q12 RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 QB(K1,IY)=RB(K2,IY)+T ! Q21 RB(K2,IY)=(RB(K2,IY)-T)*QEP ! R22 T =RB(K1,IY)+QB(K2,IY) ! R21 RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP ! Q22 RB(K0,IY)=RB(K0,IY)+T ! R0 QB(K0,IY)=QB(K0,IY)+QB(K1,IY) ! Q0 T =RB(K0,IY)-T*REC ! R21 QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC ! Q21 QB(K2,IY)=QB(K1,IY)-RB(K1,IY) ! Q2 QB(K1,IY)=QB(K1,IY)+RB(K1,IY) ! Q1 RB(K1,IY)=T -RB(K2,IY) ! R1 RB(K2,IY)=T +RB(K2,IY) ! R2 ENDDO ENDDO ENDDO MA=MA3 ENDDO IF(LN5.GT.0)THEN ! RADIX 5 NZE=N/5 RZE=W(NZE) QZE=W(NH+NZE) RZC=1.-RZE RET=RZE*RZE-QZE*QZE QET=2*RZE*QZE REC=1.-RET DO L=1,LN5 MB=MB/5 MA5=MA*5 DO J=0,MA-1 JMB=J*MB JMB2=JMB*2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 RF4=RF2*RF2-QF2*QF2 QF4=2*RF2*QF2 DO I=0,NM,MA5 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA K4=K3+MA DO IY=1,NY T =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11 RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11 QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4 ! Q14 RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4 ! R14 QB(K4,IY)=RB(K1,IY)-QB(K1,IY) ! Q24 RB(K1,IY)=RB(K1,IY)+QB(K1,IY) ! R21 QB(K1,IY)=T +RB(K4,IY) ! Q21 RB(K4,IY)=T -RB(K4,IY) ! R24 T =RB(K3,IY)*RF3-QB(K3,IY)*QF3 ! Q13 RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3 ! R13 QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2 ! Q12 RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2 ! R12 QB(K2,IY)=QB(K3,IY)+RB(K3,IY) ! Q22 RB(K3,IY)=QB(K3,IY)-RB(K3,IY) ! R23 QB(K3,IY)=RB(K2,IY)-T ! Q23 RB(K2,IY)=RB(K2,IY)+T ! R22 RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY) ! R0 QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY) ! Q0 T =RB(K4,IY)*QZE+RB(K3,IY)*QET ! R34 RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET ! R33 RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32 RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31 RB(K2,IY)=RB(K4,IY)+RB(K3,IY) ! R2 RB(K3,IY)=RB(K4,IY)-RB(K3,IY) ! R3 RB(K4,IY)=RB(K1,IY)+T ! R4 RB(K1,IY)=RB(K1,IY)-T ! R1 T =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31 QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32 QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET ! Q33 QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET ! Q34 QB(K3,IY)=QB(K2,IY)+QB(K1,IY) ! Q3 QB(K2,IY)=QB(K2,IY)-QB(K1,IY) ! Q2 QB(K1,IY)=T +QB(K4,IY) ! Q1 QB(K4,IY)=T -QB(K4,IY) ! Q4 ENDDO ENDDO ENDDO MA=MA5 ENDDO ENDIF RETURN END SUBROUTINE YCFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE) COMMON/FFT235/ LN2,LN3,LN5 DIMENSION RB(NCB,0:*),QB(NCB,0:*),W(0:*),JUMBLE(0:*) RFAC=PERIOD/N GOTO 300 ENTRY YDFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE) RFAC=1./PERIOD ! FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS: DO J=N/2+1,N-1 I=N-J DO IX=1,NX,NDX T =RB(IX,J) RB(IX,J)=RB(IX,I) RB(IX,I)=T T =QB(IX,J) QB(IX,J)=QB(IX,I) QB(IX,I)=T ENDDO ENDDO 300 IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N) IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W) NM=N-1 NH=N/2 ! SCALE AND PERMUTE THE DATA: DO IX=1,NX,NDX RB(IX,0)=RB(IX,0)*RFAC QB(IX,0)=QB(IX,0)*RFAC ENDDO DO I=1,NM J=JUMBLE(I) IF(J.GT.I)THEN DO IX=1,NX,NDX T=RB(IX,I) RB(IX,I)=RB(IX,J)*RFAC RB(IX,J)=T T=QB(IX,I) QB(IX,I)=QB(IX,J)*RFAC QB(IX,J)=T ENDDO ELSE DO IX=1,NX,NDX RB(IX,I)=RB(IX,I)*RFAC QB(IX,I)=QB(IX,I)*RFAC ENDDO ENDIF ENDDO ! TRANSFORM THE DATA: MA=1 MB=N ! RADIX 4 LS=0 DO L=2,LN2,2 MB=MB/4 LS=L MA4=MA*4 MB2=MB*2 DO J=0,MA-1 JMB=J*MB JMB2=J*MB2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 DO I=0,NM,MA4 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA DO IX=1,NX,NDX T =RB(IX,K3)*RF3-QB(IX,K3)*QF3 ! Q13 QB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3 ! R13 RB(IX,K3)=RB(IX,K2)*RF1-QB(IX,K2)*QF1 ! R12 RB(IX,K2)=RB(IX,K2)*QF1+QB(IX,K2)*RF1 ! Q12 QB(IX,K2)=QB(IX,K3)-RB(IX,K2) ! R23 RB(IX,K2)=QB(IX,K3)+RB(IX,K2) ! Q22 QB(IX,K3)=RB(IX,K3)+T ! R22 T =RB(IX,K3)-T ! Q23 RB(IX,K3)=RB(IX,K1)*RF2-QB(IX,K1)*QF2 ! R11 QB(IX,K1)=RB(IX,K1)*QF2+QB(IX,K1)*RF2 ! Q11 RB(IX,K1)=RB(IX,K0)-RB(IX,K3) ! R21 RB(IX,K0)=RB(IX,K0)+RB(IX,K3) ! R20 RB(IX,K3)=RB(IX,K1)-QB(IX,K2) ! R3 RB(IX,K1)=RB(IX,K1)+QB(IX,K2) ! R1 QB(IX,K2)=QB(IX,K0)+QB(IX,K1) ! Q20 QB(IX,K1)=QB(IX,K0)-QB(IX,K1) ! Q21 QB(IX,K0)=QB(IX,K2)+RB(IX,K2) ! Q0 QB(IX,K2)=QB(IX,K2)-RB(IX,K2) ! Q2 RB(IX,K2)=RB(IX,K0)-QB(IX,K3) ! R2 RB(IX,K0)=RB(IX,K0)+QB(IX,K3) ! R0 QB(IX,K3)=QB(IX,K1)-T ! Q3 QB(IX,K1)=QB(IX,K1)+T ! Q1 ENDDO ENDDO ENDDO MA=MA4 ENDDO IF(LS.NE.LN2)THEN ! RADIX 2 MB=MB/2 MA2=MA*2 DO J=0,MA-1 JMB=J*MB DO I=0,NM,MA2 K0=J+I K1=K0+MA RF1=W(JMB) QF1=W(NH+JMB) DO IX=1,NX,NDX T =RB(IX,K1)*QF1+QB(IX,K1)*RF1 ! Q11 QB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1 ! R11 RB(IX,K1)=RB(IX,K0)-QB(IX,K1) ! R1 RB(IX,K0)=RB(IX,K0)+QB(IX,K1) ! R0 QB(IX,K1)=QB(IX,K0)-T ! Q1 QB(IX,K0)=QB(IX,K0)+T ! Q0 ENDDO ENDDO ENDDO MA=MA2 ENDIF ! RADIX 3 REP=-.5 REC=1.5 QEP=.5*SQRT(3.) DO L=1,LN3 MB=MB/3 MA3=MA*3 DO J=0,MA-1 JMB=J*MB RF1=W(JMB) QF1=W(NH+JMB) RF2=RF1*RF1-QF1*QF1 QF2=2*RF1*QF1 DO I=0,NM,MA3 K0=I+J K1=K0+MA K2=K1+MA T =RB(IX,K2)*QF2+QB(IX,K2)*RF2 ! R12 QB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2 ! Q12 RB(IX,K2)=RB(IX,K1)*QF1+QB(IX,K1)*RF1 ! Q11 RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1 ! R11 QB(IX,K1)=RB(IX,K2)+T ! Q21 RB(IX,K2)=(RB(IX,K2)-T )*QEP ! R22 T =RB(IX,K1)+QB(IX,K2) ! R21 RB(IX,K1)=(RB(IX,K1)-QB(IX,K2))*QEP ! Q22 RB(IX,K0)=RB(IX,K0)+T ! R0 QB(IX,K0)=QB(IX,K0)+QB(IX,K1) ! Q0 T =RB(IX,K0)-T *REC ! R21 QB(IX,K1)=QB(IX,K0)-QB(IX,K1)*REC ! Q21 QB(IX,K2)=QB(IX,K1)-RB(IX,K1) ! Q2 QB(IX,K1)=QB(IX,K1)+RB(IX,K1) ! Q1 RB(IX,K1)=T -RB(IX,K2) ! R1 RB(IX,K2)=T +RB(IX,K2) ! R2 ENDDO ENDDO MA=MA3 ENDDO IF(LN5.GT.0)THEN ! RADIX 5 NZE=N/5 RZE=W(NZE) QZE=W(NH+NZE) RZC=1.-RZE RET=RZE*RZE-QZE*QZE QET=2*RZE*QZE REC=1.-RET DO L=1,LN5 MB=MB/5 MA5=MA*5 DO J=0,MA-1 JMB=J*MB JMB2=JMB*2 RF1=W(JMB) QF1=W(NH+JMB) RF2=W(JMB2) QF2=W(NH+JMB2) RF3=RF1*RF2-QF1*QF2 QF3=RF1*QF2+QF1*RF2 RF4=RF2*RF2-QF2*QF2 QF4=2*RF2*QF2 DO I=0,NM,MA5 K0=I+J K1=K0+MA K2=K1+MA K3=K2+MA K4=K3+MA T =RB(IX,K1)*QF1+QB(IX,K1)*RF1 ! Q11 RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1 ! R11 QB(IX,K1)=RB(IX,K4)*RF4-QB(IX,K4)*QF4 ! Q14 RB(IX,K4)=RB(IX,K4)*QF4+QB(IX,K4)*RF4 ! R14 QB(IX,K4)=RB(IX,K1)-QB(IX,K1) ! Q24 RB(IX,K1)=RB(IX,K1)+QB(IX,K1) ! R21 QB(IX,K1)=T +RB(IX,K4) ! Q21 RB(IX,K4)=T -RB(IX,K4) ! R24 T =RB(IX,K3)*RF3-QB(IX,K3)*QF3 ! Q13 RB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3 ! R13 QB(IX,K3)=RB(IX,K2)*QF2+QB(IX,K2)*RF2 ! Q12 RB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2 ! R12 QB(IX,K2)=QB(IX,K3)+RB(IX,K3) ! Q22 RB(IX,K3)=QB(IX,K3)-RB(IX,K3) ! R23 QB(IX,K3)=RB(IX,K2)-T ! Q23 RB(IX,K2)=RB(IX,K2)+T ! R22 RB(IX,K0)=RB(IX,K0)+RB(IX,K1)+RB(IX,K2) ! R0 QB(IX,K0)=QB(IX,K0)+QB(IX,K1)+QB(IX,K2) ! Q0 T =RB(IX,K4)*QZE+RB(IX,K3)*QET ! R34 RB(IX,K3)=RB(IX,K3)*QZE-RB(IX,K4)*QET ! R33 RB(IX,K4)=RB(IX,K0)-RB(IX,K2)*RZC-RB(IX,K1)*REC ! R32 RB(IX,K1)=RB(IX,K0)-RB(IX,K1)*RZC-RB(IX,K2)*REC ! R31 RB(IX,K2)=RB(IX,K4)+RB(IX,K3) ! R2 RB(IX,K3)=RB(IX,K4)-RB(IX,K3) ! R3 RB(IX,K4)=RB(IX,K1)+T ! R4 RB(IX,K1)=RB(IX,K1)-T ! R1 T =QB(IX,K0)-QB(IX,K1)*RZC-QB(IX,K2)*REC ! Q31 QB(IX,K2)=QB(IX,K0)-QB(IX,K2)*RZC-QB(IX,K1)*REC ! Q32 QB(IX,K1)=QB(IX,K3)*QZE-QB(IX,K4)*QET ! Q33 QB(IX,K4)=QB(IX,K4)*QZE+QB(IX,K3)*QET ! Q34 QB(IX,K3)=QB(IX,K2)+QB(IX,K1) ! Q3 QB(IX,K2)=QB(IX,K2)-QB(IX,K1) ! Q2 QB(IX,K1)=T +QB(IX,K4) ! Q1 QB(IX,K4)=T -QB(IX,K4) ! Q4 ENDDO ENDDO MA=MA5 ENDDO ENDIF RETURN END