关于用小波分析的程序做周期分析。。。

殷2333

应气学渣一枚。。。老师让对气温等气象要素的年变化求气候倾向率和做周期分析。。。
     给了小波的程序，让求周期分析。。。下面是程序。。。
PROGRAM   WAVELET_MORLET
      PARAMETER(dt=1.)
        PARAMETER(MN=52,MM=25,INCRFC=1)
        PARAMETER(KVT=2)
      DOUBLE PRECISION  FF(MN),F(MN),PA,A,T,B,RMO,RIMO,A1(MN)
        DOUBLE PRECISION  RW(MN,MM),IW(MN,MM),W1(MN,MM),wx(mn,mm)
      REAL FF1(mN),W2(MN,MM),wx1(mn,mm),rw1(mn,mm),AA(MN,KVT)
     &,da(mm)
        PA=3.14159D0
       
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccc
ccc            REading the source       
        OPEN(20,FILE='f:\xb\0.txt')
        do j=1,mn
      READ(20,*)A1(J)
        enddo
        do j=1,mn
      A1(J)=A1(J)/10
        enddo
c      close(20)
c      DO I=1,MN
c        A1(I)=AA(1,I)
c        ENDDO
        WRITE(*,*) 'FINISHED READING DATA !'
      DO I=1,MN
       FF(I)=DBLE(A1(I))
      ENDDO
!-----------------------------------------
      CALL NORMAL(MN,FF,F)
!-----------------------------------------
!-----------------------------------------
**      CALL FILTER(F,MN,dt,PERC,PERL)
!-----------------------------------------
      DO IA=1,MM
                DO IB=1,MN
       A=DBLE(1.*IA*INCRFC)  ! Frequency/cycle parameter!不明白
       B=DBLE(IB-(MN+1)/2.)   ! Time parameter
       RW(IB,IA)=0.0d0
       IW(IB,IA)=0.0d0
       DO IT=1,MN

        T=DBLE(IT-(MN-1)/2-1)

        CALL RMORLET(T,A,B,PA,RMO)
        CALL RIMORLET(T,A,B,PA,RIMO)
       
      RW(IB,IA)=RW(IB,IA)+DBLE(F(IT))*RMO
      IW(IB,IA)=IW(IB,IA)-DBLE(F(IT))*RIMO
        ENDDO
      RW(IB,IA)=RW(IB,IA)/DSQRT(A)
      IW(IB,IA)=IW(IB,IA)/DSQRT(A)
      W1(IB,IA)=RW(IB,IA)**2+IW(IB,IA)**2
      W2(IB,IA)=sngl(W1(IB,IA))
        RW1(IB,IA)=sngl(RW(IB,IA))
      WX(IB,IA)=DATAN(IW(IB,IA)/RW(IB,IA))
        WX1(IB,IA)=SNGL(WX(IB,IA))
      ENDDO;ENDDO
*        do 111 ia=1,mm
*        var=0.0
*        do 112 ib=1,mn
*        var=var+w1(ib,ia)**2
             
*112        continue
*        write(*,*)ia
*        write(222)var
*111        continue

!---------------------------------------------------
ccccc小波的模方差贡献ccccc
        open(120,file='f:\xb\mo.grd',
     &        form='binary')
!        do i=1,mn
!     write(*,10) (w2(i,j),J=1,mm)
!        enddo

        do j=1,MM
                da(j)=0.0
                do i=1,MN
                        da(j)=da(j)+w1(i,j)
                enddo
        enddo
        do i=1,MM
        write(120) da(i)
        enddo
CCCCCCCCCCCCCCCCCCC写小波分析的不同时间周期的功率谱
!        open(122,file='g:\xb\xb\pu.txt',form='formatted')
       
        do j=1,MM
!        if((j.eq.5).or.(j.eq.9).or.(j.eq.15))then
       
!        write(122,'(47f6.2)') (iw(i,j),i=1,MN)
       
!     endif
        enddo


ccccc小波的模ccccc
      open(119,file='f:\xb\pattern1.grd',form='binary')
        write(119) ((w2(i,j),I=1,mn),J=1,mm)
cccccc小波的实部ccccc  
        open(188,file='f:\xb\real-j1.grd',form='binary')
        write(188) ((rw1(i,j),J=1,mm),I=1,mn)
cccccc小波的位相ccccc
      open(189,file='f:\xb\phase1.grd',form='binary')
      write(189) ((wx1(i,j),I=1,mn),j=1,mm)
cccccc小波的虚部ccccc     
        open(121,file='f:\xb\imaginary1.grd',form='binary')
        write(121) ((iw(i,j),I=1,mn),J=1,mm)

        close(119)
        close(188)
        close(189)
        close(121)
!---------------------------------------------------
      END
!---------------------------------------------
        SUBROUTINE RMORLET(T,A,B,PA,RMO)
       DOUBLE PRECISION Q,TT,PA,T,A,B,RMO
       TT=(T-B)/A
   !    Q=DBLE(1.0)
       Q=-((ABS(TT))**2)
       Q=Q/2           
       RMO=DCOS(2*PA*TT)*DEXP(Q)
        RETURN
       END
!---------------------------------------------
         SUBROUTINE  RIMORLET(T,A,B,PA,RIMO)
       DOUBLE PRECISION Q,TT,T,A,B,pa,RIMO
       TT=(T-B)/A
   !    Q=DBLE(1.0)
       Q=-((ABS(TT))**2)
       Q=Q/2.0
       RIMO=DSIN(2*PA*TT)*DEXP(Q)
        RETURN
       END
!---------------------------------------------
       SUBROUTINE NORMAL(N,X,XX)
       DOUBLE PRECISION X(N),XX(N),PX,S
       PX=0.d0
       DO I=1,N
        PX=PX+X(I)
       ENDDO
         PX=PX/FLOAT(N)
       S=0.d0
       DO I=1,N
        S=S+(X(I)-PX)**2
         ENDDO
       S=DSQRT(S/FLOAT(N))
       DO I=1,N
        XX(I)=(X(I)-PX)/S
       ENDDO
       END
!--------------------------------------------
      SUBROUTINE FILTER(AA,ITX,DT,PERC,PERL)
      DOUBLE PRECISION AA(ITX),WK1(ITX),WK2(ITX)
        DOUBLE PRECISION PI2,WO,WOS,W1,W2,DW,DWS,B1,B2,B3
        PI2=2.0*3.1415926d0
        WO=PI2/PERC
        WOS=WO*WO
        W1=PI2/PERL
        W2=WOS/W1
        DW=2.0*ABS((SIN(W1*DT))/(1.0+COS(W1*DT))-(SIN(W2*DT))/(1.0+
     &   COS(W2*DT)))
      DWS=(4.0*SIN(W1*DT)*SIN(W2*DT))/((1.0+COS(W1*DT))*(1.0+
     &    COS(W2*DT)))
      B3=(2.0*DW)/(4.0+2.0*DW+DWS)
        B2=(4.0-2.0*DW+DWS)/(4.0+2.0*DW+DWS)
        B1=(2.0*(DWS-4.0))/(4.0+2.0*DW+DWS)
      ITXM=ITX-1
        DO IT=1,ITX
       WK1(IT)=0.0
        ENDDO
        DO IT=3,ITX
       WK1(IT)=B3*(AA(IT)-AA(IT-2))-B1*WK1(IT-1)-B2*WK1(IT-2)
        ENDDO
        WK2(ITX-1)=WK1(ITX-1)
        WK2(ITX)=WK1(ITX)
        WK1(1)=AA(1)
        WK1(2)=AA(2)
        DO IT=2,ITXM
       II=ITX-IT
       WK2(II)=B3*(WK1(II)-WK1(II+2))-B1*WK2(II+1)-B2*WK2(II+2)
        ENDDO
        DO IT=1,ITX
       AA(IT)=WK2(IT)
        ENDDO
        END
            

发现其中让打开多个文件。。。不知道是什么意思？