我的程序如下：

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

2、编制了主程序：MAIN.F90，内容如下：
主程序：
DIMENSION DATA(30000),F(15000),G(15000),R(15000)
DIMENSION FREQ(15000),TAU(15000)
!
open(5,file='ua0258.DAT')
open(6,file='OUTPUT.DAT')
!
!Do i=1,11
!read(5,*)
!enddo
!
READ(5,501) NN,DT,(DATA(M),M=1,NN)
CALL FPAC(NN,DATA,30000,DT,111,F,G,R,15000,NFOLD,DF)
!
DO 210 J=1,NFOLD
TAU(J)=REAL(J-1)*DT
FREQ(J)=J/(NFOLD*2)*DT
210 CONTINUE
!
WRITE (6,601) (MK-1,FREQ(MK),F(MK),G(MK),TAU(J),R(MK),MK=1,NFOLD)
STOP
!
501 FORMAT(T16,I5/T30,F5.3//8F15.11)
601 FORMAT(T35,'EXAMPLE WAVE'//T15,'--FOURIER, POWER SPECTRUM AND AUTOCORRELATION--'&
          ///T5,'M',TR3,'FREQUENCY',TR7,'FOURIER SPECTRUM',TR2,'POWER SPECTRUM',&
          TR2,'TIME LAG',TR2,'AUTOCORRELATION'//(T5,I5,TR5,F10.3,TR5,F10.3,TR5,F10.3,&
          TR5,F10.3,TR5,F10.3))
END

3、读入的数据文件的格式如下：
NO. OF POINTS:  24256 
EQUALLY SPACED INTERVALS OF:  0.005  SEC

-1.072970E-001 -1.072970E-001 -9.958957E-002 -9.958957E-002 -1.072970E-001 -1.034439E-001 -9.573565E-002 -1.034439E-001
-9.958957E-002 -1.034439E-001 -1.072970E-001 -9.958957E-002 -9.958957E-002 -1.111513E-001 -9.958957E-002 -9.958957E-002
-1.072970E-001 -1.034439E-001 -1.072970E-001 -1.072970E-001 -1.072970E-001 -9.958957E-002 -1.072970E-001 -1.111513E-001

4、现在出现的问题是，运行之后，报错显示主程序的第11行有错：input conversion error
出错报告截图：




转载本文请联系原作者获取授权，同时请注明本文来自黄艳科学网博客。
链接地址：https://blog.sciencenet.cn/blog-456941-455853.html

上一篇：阅读手记之五——地震动的谱分析入门
下一篇：阅读手记之七——地震动的谱分析入门