1、编制子程序:FFT.F90——用于快速傅里叶变换 FPAC.F90——用于求傅里叶谱,功率谱和自相关系数 FPAC子程序: !****************************************************************** ! SUBROUTINE FOR FOURIER, POWER SPECTRA AND ZUTOCORRELATION !****************************************************************** ! SUBROUTINE FPAC(N,X,ND1,DT,IND,F,G,R,ND2,NFOLD,DF) COMPLEX C(32678) DIMENSION X(ND1),F(ND2),G(ND2),R(ND2) ! ! INITIALIZATION ! DO 110 M=1,N C(M)=CMPLX(X(M),0.) 110 CONTINUE NT=2 120 IF(NT.GE.N) GO TO 130 NT=NT*2 GO TO 120 130 IF(NT.EQ.N) GO TO 150 DO 140 M=N+1,NT C(M)=(0.,0.) 140 CONTINUE 150 NFOLD=NT/2+1 T=REAL(NT)*DT DF=1./T ! ! FOURIER TRANSFORM ! CALL FAST(NT,C,32678,-1) ! ! FOURIER SPECTRUM ! IF(IND.EQ.001) GO TO 180 DO 160 K=1,NFOLD F(K)=CABS(C(K))*DT 160 CONTINUE IF(IND.EQ.100) RETURN ! ! POWER SPECTRUM ! IF(IND.EQ.101) GO TO 180 G(1)=F(1)**2/T DO 170 K=2,NFOLD-1 G(K)=2.*F(K)**2/T 170 CONTINUE G(NFOLD)=F(NFOLD)**2/T IF(MOD(IND,10).EQ.0) RETURN ! ! AUTOCORRELATION ! 180 DO 190 K=1,NT C(K)=C(K)*CONJG(C(K)) 190 CONTINUE CALL FAST(NT,C,32678,+1) R0=REAL(C(1)) DO 200 J=1,NFOLD R(J)=REAL(C(J))/R0 200 CONTINUE RETURN END
FFT子程序:
!************************************************* ! SUBROUTINE FOR FAST FOURIER TRANSFORM !************************************************* ! ! ! SUBROUTINE FAST(N,X,ND,IND) ! COMPLEX X(ND),TEMP,THETA ! J=1 DO 140 I=1,N IF(I.GE.J) GO TO 110 TEMP=X(J) X(J)=X(I) X(I)=TEMP 110 M=N/2 120 IF(J.LE.M) GO TO 130 J=J-M M=M/2 IF(M.GE.2) GO TO 120 130 J=J+M 140 CONTINUE KMAX=1 150 IF(KMAX.GE.N) RETURN ISTEP=KMAX*2 DO 170 K=1,KMAX THETA=CMPLX(0.0,3.141593*REAL(IND*(K-1)))/REAL(KMAX) DO 160 I=K,N,ISTEP J=I+KMAX TEMP=X(J)*CEXP(THETA) X(J)=X(I)-TEMP X(I)=X(I)+TEMP 160 CONTINUE 170 CONTINUE KMAX=ISTEP GO TO 150 END