- 积分
- 852
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-12-4
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
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=66,MM=33,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='qv.txt')
READ(20,*)(A1(J),J=1,MN)
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) !
B=DBLE(IB-(MN+1)/2.)
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:\paper\fangcha.txt',
& form='formatted')
! 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='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='mod.grd',form='binary')
write(119) ((w2(i,j),I=1,mn),J=1,mm)
cccccc小波的实部ccccc
open(188,file='real.grd',form='binary')
write(188) ((rw1(i,j),I=1,mn),J=1,mm)
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[img=0,1]file:///C:\Users\Administrator\AppData\Roaming\Tencent\Users\479266557\QQ\WinTemp\RichOle\~%~8[D6UPJV~UJOQF{LGA$S.png[/img]
END
这是小波分析 的程序只要运行就会出现如图情况
想问问呢哪里出错了呢?
|
-
|