傲娇师兄开贴帮师妹调程序（暑假归来，继续）

jrzhu1108 · 发表于 2014-8-26 20:54:51

{:eb513:}{:eb513:}{:eb513:}

2006202788 · 发表于 2014-8-27 08:40:55

谢谢分享不胜感激

roland · 发表于 2014-8-28 18:41:07

太好了，谢谢分享

邗珏 · 发表于 2014-8-30 21:53:39

请问楼主对配套魏凤英老师书的最大熵谱的fortran程序有研究吗？
用配套的数据，程序也会出错，数组溢出了貌似
能请楼主帮忙看看吗？

MESA.FOR (2.93 KB, 下载次数: 0)

邗珏 · 发表于 2014-8-30 21:59:22

菜鸟，不懂怎么把附件设成免费唉
我还是直接帖程序吧

C    THIS IS A PROGRAM OF MAXIMUM ENTROPY SPECTRUM
C    ANALYSIS (MESA) FOR A TIME SERIES
PROGRAM MESA
PARAMETER (N=44, M=22)
DIMENSION X(N),T(N),SE(N,M)
DIMENSION B1(N),B2(N),A(N),AA(N),FI(M),SIGMA(N)
CHARACTER*9 FNAME1,FNAME2
C    ***************************************************
C    * N: SAMPLE SIZE OF THE TIME SERIES             *
C    * M: MAXIMUM LAG TIME LENGTH                   *
C    * X(N): TIME SERIES                            *
C    * T(N):LENGTH OF PERIOD                         *
C    * SE(N,M): MAXIMUM ENTROPY SPECTRUM VALUES       *
C    ***************************************************
WRITE(*,*) 'PLEASE ENTER THE DATA FILE NAME:'
READ(*,'(A)') FNAME1
WRITE(*,*) 'PLEASE ENTER THE OUTPUT FILE NAME:'
READ(*,'(A)') FNAME2
NW2=1
OPEN (7,FILE=FNAME1)
READ (7,*,END=9)(X(I),I=1,N)
9    CONTINUE
CLOSE (7)
CALL SMESA(X,N,M,NW2,K0,T,SE,FI,SIGMA,A,AA,B1,B2)
OPEN (8,FILE=FNAME2)
WRITE(8,80)
  80 FORMAT(20X,'*******MAXIMUM ENTROPY SPECTRUM ANALYSIS*******'/)
WRITE (8,90)
WRITE (8,100) (L,T(L),SE(L,K0), L=1,M)
CLOSE(8)
  90 FORMAT (1X,7X,1HL,4X,1HT,10X,2HSE/,
   &       8X,23H-----------------------)
100 FORMAT (1X,I8,F7.2,E15.7)
STOP
END

CCC------------------------------END------------------------------------CCC

SUBROUTINE SMESA(X,N,K,NW,K0,T,SE,FI,SIGMA,A,AA,B1,B2)
DIMENSION X(N),T(N),SE(N,K)
DIMENSION B1(N),B2(N),A(N),AA(N),FI(K),SIGMA(N)
C
C----------------------------------------------------------------------C
C
IF (NW.EQ.0) GOTO 10
WRITE(*,'(/9X,1Hk,6X,5HFI(K))')
10 S=0.
DO 50 I=1,N
50 S=S+X(I)
S=S/N
DO 60 I=1,N
X(I)=X(I)-S
60 CONTINUE
S=0.
DO 11 I=1,N
11 S=S+X(I)*X(I)
SIGMA(1)=S/N
M=1
B1(1)=X(1)
B2(N-1)=X(N)
DO 12 J=2,N-1
B1(J)=X(J)
12 B2(J-1)=X(J)
IF (NW.EQ.0) GOTO 301
WRITE(*,'(I10,F12.3)') M, SIGMA(1)
  301 GOTO 101
  102 M=M+1
DO 13 J=1,M-1
13 AA(J)=A(J)
DO 14 J=1,N-M
B1(J)=B1(J)-AA(M-1)*B2(J)
14 B2(J)=B2(J+1)-AA(M-1)*B1(J+1)
  101 COM=0.
DEN=0.
DO 15 J=1,N-M
COM=COM+B1(J)*B2(J)
15 DEN=DEN+B1(J)**2+B2(J)**2
A(M)=2.0*COM/DEN
SIGMA(M+1)=SIGMA(M)*(1-A(M)**2)
IF(M.EQ.1)GOTO 102
  103 DO 16 J=1,M-1
A(J)=AA(J)-A(M)*AA(M-J)
16 CONTINUE
FI(M) = (FLOAT(N+M)/FLOAT(N-M))*SIGMA(M+1)
IF (FI(M-1).lt.FI(M-2).and.FI(M-1).lt.FI(M)) then
K0=M-1
GOTO 555
   ELSE
      CONTINUE
   ENDIF
IF (NW.EQ.0) GOTO 300
WRITE(*,'(I10,F12.3)') M, FI(m)
  300 IF (M.EQ.K) GOTO 555
DO 19 I=1,K
F1=0.
F2=0.
DO 18 J=1,M
E=2.0*3.1416*FLOAT(J)*FLOAT(I)/FLOAT(N)
F1=F1+A(J)*COS(E)
18 F2=F2+A(J)*SIN(E)
SE(I,M) = SIGMA(M+1)/((1.0-F1)**2+F2**2)
19 CONTINUE
IF (M.LT.K) GOTO 102
  555 DO 150 I=1,K
T(I) = 2.*K/float(I)
  150 CONTINUE
IF (NW.EQ.0) GOTO 999
WRITE (*,'(/,12x,4Hk0 =,i4/)') K0
  999 RETURN
END
CCC-----------------------------END----------------------------------

邗珏 · 发表于 2014-8-30 22:00:34

运行程序出现的错误截图也一并贴上

运行程序出先的错误

邗珏 · 发表于 2014-8-31 23:26:08

今天又仔细看了，
似乎问题出在103上面一行: IF(M.EQ.1)GOTO 102
这样M=2的时候运行到下面这行
IF (FI(M-1).lt.FI(M-2).and.FI(M-1).lt.FI(M)) then
FI(M-2)就是FI(0)了，显然是超出数组界限了
把 IF(M.EQ.1)GOTO 102改成 IF(M.LT.3)GOTO 102可以解决问题
不知道对不对

虫儿飞 · 发表于 2014-9-2 15:51:56

邗珏发表于 2014-8-31 23:26
今天又仔细看了，
似乎问题出在103上面一行: IF(M.EQ.1)GOTO 102
这样M=2的时候运行到下面这行

好几天没登陆了
你为了让程序能够运行这样修改是没有问题的，但是是不是符合算法我还要再翻一下书。

邗珏 · 发表于 2014-9-2 23:47:35

虫儿飞发表于 2014-9-2 15:51
好几天没登陆了
你为了让程序能够运行这样修改是没有问题的，但是是不是符合算法我还要再翻一下书。

先多谢楼主

wggyhq · 发表于 2014-9-3 17:35:06

这个师兄太好了

		自动登录	找回密码
密码			立即注册

[经验总结] 傲娇师兄开贴帮师妹调程序（暑假归来，继续）