小波分析实部图无法显示长周期

自由行

画出来的小波分析实部图没有显示出长周期，不知道是怎么回事，是在fortran里面需要该参数吗，比较着急想请教各位大神有何方法喃。太感谢了。贴出fortran源程序,不确定能不能在红色标注那里更改mm的值，自己试过更改反而出图更有问题，很是奇怪，怎么才能显示完她的长周期
CCCC**************************************************************CCC      !极端降水量
CCCC       MORLET WAVELET ANALYSIS FOR SPECTRUM DISTRIBUATION     CCC
CCCC                       OF TIME SERIES              
CCC
CCCC                ******PARAMETER  TABLE******                  CCC
CCCC PA===> CONSTANT PI=3.1425926;                                CCC
CCCC MN===> THE LENGTH OF TIME SERIES ; MM===> WINDOWS LENGTH     CCC
CCCC  F===> PRIMATIVE DATA ; RW===> WINDOWS BLOCK                 CCC  
CCCC INCRFC===> FREQUENCY/CYCLE OF INTERVAL                       CCC
CCCC  normal===>将资料标准化；RMORLET===>小波变换（实部）         CCC
CCCC                          RIORLET===>小波变换（虚部）         CCC
CCCC 本程序使用的MOVLET 小波变换                                  CCC
CCCC______________________________________________________________CCC   
CCCC______________________________________________________________CCC
CCC                                                               CCC
CCC THIS PROGRAM IS FROM Xv Jianjun,NANJING UNIVERSITY,01/05/1996 CCC
CCC             NOTED BY yrl,NIM,16/3/2002           CCC
CCCC                  COPYRIGHT RESERVED  !!!!!!                  CCC
CCCC**************************************************************CCC
      PROGRAM   WAVELET_MORLET
      PARAMETER(dt=1.)
PARAMETER(MN=55,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='R95p.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='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='pattern1.grd',form='binary')
write(119) ((w2(i,j),I=1,mn),J=1,mm)
cccccc小波的实部ccccc  
open(188,file='real-j1.grd',form='binary')
write(188) ((rw1(i,j),J=1,mm),I=1,mn)
cccccc小波的位相ccccc
      open(189,file='phase1.grd',form='binary')
      write(189) ((wx1(i,j),I=1,mn),j=1,mm)
cccccc小波的虚部ccccc     
open(121,file='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
            

极地冷高压

居然没人回复，楼主图画的漂亮，水一下

自由行

极地冷高压 发表于 2016-10-23 15:13
居然没人回复，楼主图画的漂亮，水一下

ai，谢谢，您知道这是为什么吗