C PROGRAM 13.3 TVSPC C C ... CHANGING SPECTRUM ... C C INPUTS: C N: NUMBER OF POINTS C M: AR ORDER C K: TREND ORDER C NOBS: SAMPLING INTERVAL C IVAR=1: FOR VARIANCE CORRECTION C AR(I,J): TIME VARYING AR COEFFICIENT C VAR(I): TIME VARYING VARIANCE C PARAMETRS: C NMAX: ADJUSTABLE DIMENSION OF AR (NMAX >= N) C MJ: ADJUSTABLE DIMENSION OF AR (MJ >= M) C NF: NUMBER OF FREQUENCIES C @TEST.UNI: 5/16/89 1/07/91 9/09/92 C PARAMETER( NMAX=200,NF=200,MJ=10,IDEV1=2,IDEV2=3 ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72, VNAME*8 DIMENSION AR(MJ,NMAX), VNAME(5), VALUE(5), VAR(3000) COMMON /CMDATA/ TITLE DATA VAR /3000*1.0D0/ C OPEN( IDEV1,FILE='A:$LASERLIB$TEMP3.DAT' ) READ(IDEV1,1) TITLE READ(IDEV1,*) N, M, K, NOBS, IVAR READ(IDEV1,*) SIG2 DO 10 J=1,N 10 READ(IDEV1,*) (AR(I,J),I=1,M) CLOSE( IDEV1 ) NN = N*NOBS IF(IVAR.EQ.1 ) OPEN( IDEV2,FILE='A:$LASERLIB$TEMP4.DAT' ) IF(IVAR.EQ.1 ) READ( IDEV2,*) (VAR(I),I=1,NN) IF(IVAR.EQ.1 ) CLOSE( IDEV2 ) 1 FORMAT( A72 ) C C ... PLOT CHANGING SPECTRUM ... C C VNAME(1) = 'N =' VNAME(2) = 'M =' VNAME(3) = 'K =' VNAME(4) = 'NOBS =' VNAME(5) = 'SIG2 =' VALUE(1) = N VALUE(2) = M VALUE(3) = K VALUE(4) = NOBS VALUE(5) = SIG2 CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 12.2: CHANGING SPECTRUM',32,5,VNAME, * VALUE ) CALL PLOT( 3.0,2.0,-3 ) CALL PT3DSP( AR,SIG2,MJ,M,N,NOBS,VAR ) CALL PLOTE C call plot( 0.0,0.0,999 ) STOP E N D SUBROUTINE PT3DSP( A,SIG2,MJ,M,N,NOBS,VAR ) C C ... COMPUTE AND DRAW TIME VARYING SPECTRUM ... C C INPUTS: C A,SIG2,M: AR MODELS C N: DETA LENGTH C NOBS: LOCAL STATIONARY SPAN C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(MJ,N), SP(0:200), YH(2000), VAR(3000), B(20) ANGLE = 70.0 NF = 200 C WX = 5.0 WY = 3.0 WZ =10.0 X0 = 0.0 X1 = 0.5 DX = 0.5 Z0 = 0.0 Z1 = N*NOBS DZ = 500.0 K = 2 CALL ARMASP( A,M,B,0,SIG2,NF,SP ) CALL MAXMIN( SP(0),NF+1,Y0,Y1,DY ) C C ... DRAW 3D-AXIS ... C CALL PLOT3A( NF,M,WX,WY,WZ,ANGLE,X0,X1,DX,Y0,Y1,DY,Z0,Z1, * DZ,KN,K,YH ) DO 100 J=1,N C C ... SPECTRUM OF AR MODEL AT TIME J*NOBS ... C CALL ARMASP( A(1,J),M,B,0,SIG2,NF,SP ) C C ... DRAW SPECTRUM AND SHIFT THE ORIGIN ... C DO 10 I=0,NF 10 SP(I) = SP(I) + DLOG10( VAR(J*NOBS-NOBS/2) ) CALL PLOT3B( SP(0),NF+1,N,WX,WY,WZ,ANGLE,Y0,Y1,K,YH,ZS,ZC ) 100 CONTINUE C CALL PLOT( -SNGL(ZC),-SNGL(ZS),-3 ) C RETURN E N D SUBROUTINE PLOT3A( N,M,WX,WY,WZ,ANGLE,X0,X1,DX,Y0,Y1,DY,Z0,Z1, * DZ,KN,K,YH ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION YH(2000) C PI = 3.1415926535D0 THETA = ANGLE * ( PI/180.0 ) C CALL AXISXY( 0.0D0,0.0D0,WX,WY,X0,X1,Y0,Y1,DX,DY,0.2D0,0,5,5 ) CALL AXISZ( WX,WZ,THETA,Z0,Z1,DZ,0.2D0,1 ) ZS = WZ*DSIN(THETA) KN = K*(N-1) + 1 DO 30 I=1,2000 30 YH(I) = Y0 + 2*(ZS/N)/(WY/(Y1-Y0)) RETURN E N D SUBROUTINE PLOT3B( Y,N,M,WX,WY,WZ,ANGLE,Y0,Y1,K,YH,ZS,ZC ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N+1), YH(2000) C PI = 3.1415926535D0 THETA = ANGLE * ( PI/180.0 ) C ZS = WZ*DSIN(THETA) ZC = WZ*DCOS(THETA) KN = K*(N-1) + 1 SHIFT = (ZC/M) / (WX/(KN-1)) IS = SHIFT DS = SHIFT-IS KN1 = KN - IS DO 40 I = 1,KN1 L = I+IS 40 YH(I) = YH(L)+ ((YH(L+1)-YH(L))*DS) - (ZS/M)/(WY/(Y1-Y0)) DO 50 I=KN1+1,KN 50 YH(I) = Y0 DO 60 I=1,KN L = (I-1)/K + 1 IF( I.LT.KN ) W = Y(L) +(( Y(L+1) - Y(L) )*MOD(I-1,K))/K IF( I.EQ.KN ) W = Y(L) IF(L.EQ.N) W = Y(L) 60 YH(I) = DMAX1( YH(I),W) ZC1 = ZC/M ZS1 = ZS/M CALL PLOT( SNGL(ZC1),SNGL(ZS1),-3 ) CALL PLOTY( YH,KN,Y0,Y1,WX,WY,0,1 ) C RETURN E N D SUBROUTINE AXISZ( WX,WZ,THETA,Z0,Z1,DZ,DWC,IZ ) IMPLICIT REAL*8(A-H,O-Z) REAL*4 WC C CALL NEWPEN(2) CALL PLOT( SNGL(WX),0.0,3 ) CALL PLOT( SNGL(WX+WZ*DCOS(THETA)),SNGL(WZ*DSIN(THETA)),2 ) C WC = DWC NZ = (Z1-Z0)/DZ + 0.001 MZ = -DLOG10(DZ) IF ( MZ.LE.0) MZ = -1 D = WZ*DZ /( Z1-Z0) DO 10 I = 1,NZ+1 ZN = Z0 + DZ*(I-1) ZZ = D*(I-1) CALL PLOT(SNGL(WX+ZZ*DCOS(THETA)),SNGL(ZZ*DSIN(THETA)),3 ) CALL PLOT(SNGL(WX+ZZ*DCOS(THETA)+(WC/2)),SNGL(ZZ*DSIN(THETA)),2) IF( MOD(I,IZ).EQ.1 .OR. IZ.EQ.1 ) * CALL NUMBER( SNGL(WX+ZZ*DCOS(THETA)+2*WC),SNGL(ZZ*DSIN(THETA) * -(WC/2)),WC,ZN,0.0,MZ ) 10 CONTINUE CALL NEWPEN( 1 ) RETURN E N D SUBROUTINE AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE ARMASP( A,M,B,L,SIG2,NF,SP ) C C ... Logarithm of the power spectrum of the ARMA model ... C C Inputs: C M: AR order C L: MA order C A(I): AR coefficient C B(I): MA coefficient C SIG2: Innovation variance C NF: Number of frequencies C Output: C SP(I): Power spectrum (in log scale) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(M), B(L) DIMENSION SP(0:NF), H(0:500), FR(0:500), FI(0:500) C H(0) = 1.0D0 DO 10 I=1,M 10 H(I) = -A(I) C CALL FOURIE( H,M+1,NF+1,FR,FI ) C DO 20 I=0,NF 20 SP(I) = SIG2/( FR(I)**2 + FI(I)**2 ) C H(0) = 1.0D0 DO 30 I=1,L 30 H(I) = -B(I) CALL FOURIE( H,L+1,NF+1,FR,FI ) DO 40 I=0,NF 40 SP(I) = SP(I)*( FR(I)**2 + FI(I)**2 ) C DO 50 I=0,NF 50 SP(I) = DLOG10( SP(I) ) C RETURN E N D SUBROUTINE FOURIE( X,N,M,FC,FS ) C C ... Discrete Fourier transformation by Goertzel method ... C C Inputs: C X(I): data (I=1,N) C N: data length C M: number of Fourier components C FC(J): Fourier cosine transform (J=1,M) C FS(J): Fourier sine transform (J=1,M) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N), FC(M), FS(M) DATA PI/3.14159265358979D0/ C W = PI/(M-1) DO 20 I=1,M CI = DCOS(W*(I-1)) SI = DSIN(W*(I-1)) T1 = 0.0 T2 = 0.0 DO 10 J=N,2,-1 T0 = 2*CI*T1 - T2 + X(J) T2 = T1 10 T1 = T0 FC(I) = CI*T1 - T2 + X(1) 20 FS(I) = SI*T1 C RETURN E N D