C PROGRAM 6.2 MARSPC C C ... This program computes cross spectra and noise contribution ... C C Inputs: C IDEV: Input device C M: AR order C L: Dimension of the time series C E(I,J): Innovation covariance matrix C A(II,I,J): AR coefficient matrices C Parametrs: C NF: Number of frequencies C MJ: Adjustable dimension C MJ1: Adjustable dimension C Written by G.K. PARAMETER( NF=200,MJ=3,MJ1=7,IDEV=1 ) IMPLICIT REAL*8(A-H) IMPLICIT COMPLEX*16(O-Z) DIMENSION A(MJ,MJ,MJ1), E(MJ,MJ) DIMENSION P(0:NF,MJ,MJ), FNC(0:NF,MJ,MJ), COH(0:NF,MJ,MJ) DIMENSION AMP(0:NF,MJ,MJ), ANG(0:NF,MJ,MJ) C open( IDEV,file='ar2.dat' ) READ(IDEV,*) M, L DO 10 I=1,L 10 READ(IDEV,*) (E(I,J),J=1,L) DO 20 II=1,M DO 20 I=1,L 20 READ(IDEV,*) (A(I,J,II),J=1,L) close( IDEV ) C C ... Cross spectrum ... C CALL MARSPC( M,L,A,E,NF,MJ,MJ1,P,FNC,AMP,ANG,COH ) C C ... Print out plot cross spectra, coherency and noise C contribution ... C CALL PRMSPC( M,L,NF,MJ,P,FNC,COH ) C CALL PTCSPC( M,L,NF,MJ,P,FNC,AMP,ANG,COH ) C STOP E N D SUBROUTINE PRMSPC( M,L,NF,MJ,P,FNC,COH ) C C ... Print out cross spectra and relative noise contribution ... C C Inputs: C M: AR order C L: Dimension of the time series C NF: Number of frequencies C MJ: Adjustable dimension of P and FNC C P(II,I,J): Cross spectrum of Y(I) and Y(J) at II-th frequency C FNC(II,I,J): Noise contribution of Y(J) to Y(I) C COH(II,I,J): Simple coherency C IMPLICIT COMPLEX*16(O-Z) IMPLICIT REAL*8(A-H) DIMENSION P(0:NF,MJ,MJ), FNC(0:NF,MJ,MJ), FRNC(0:200) DIMENSION COH(0:NF,MJ,MJ) C WRITE(6,600) WRITE(6,610) M, L, NF WRITE(6,620) DO 10 I=1,L DO 10 J=1,L WRITE(6,630) I, J 10 WRITE(6,640) (P(II,I,J),II=0,NF) WRITE(6,650) DO 20 I=1,L-1 DO 20 J=I+1,L WRITE(6,630) I, J 20 WRITE(6,660) (COH(II,I,J),II=0,NF) WRITE(6,670) DO 40 I=1,L DO 40 J=1,L WRITE(6,630) I, J DO 30 II=0,NF IF(J.EQ.1) FRNC(II) = FNC(II,I,J)/FNC(II,I,L) 30 IF(J.NE.1) FRNC(II) = (FNC(II,I,J)-FNC(II,I,J-1))/FNC(II,I,L) 40 WRITE(6,660) (FRNC(II),II=0,NF) C 600 FORMAT( 1H ,'PROGRAM 6.2: SPECTRUM ANALYSIS BY MAR MODEL' ) 610 FORMAT( 1H ,'M(ORDER) =',I3,5X,'L(DIMENSION) =',I2,5X, * 'NF(NUMBER OF FREQUENCIES =',I3 ) 620 FORMAT( 1H ,'CROSS SPECTRUM' ) 630 FORMAT( 1H ,'*** I =',I2,2X,'J =',I2,' ***' ) 640 FORMAT( 2(' (',D13.5,',',D13.5,')',4X) ) 650 FORMAT( 1H ,'COHERENCY' ) 660 FORMAT( 1H ,10F7.4 ) 670 FORMAT( 1H ,'RELATIVE NOISE CONTRIBUTION' ) RETURN E N D SUBROUTINE MARSPC( M,L,A,E,NF,MJ,MJ1,P,FNC,AMP,ANG,COH ) C C ... Cross spectra, amplitude, phase and noise contribution ... C C Inputs: C M: AR order C L: Dimension of the time series C E(I,J): Innovation covariance matrix C A(I,J,II): AR coefficient matrices C NF: Number of frequencies C MJ: Adjustable dimension C MJ1: Adjustable dimension C Outputs: C P(II,I,J): Cross spectra C FNC(II,I,J): Noise contribution of Y(J) to Y(I) C COH(II,I,J): Simple coherency C AMP(II,I,J): Amplitude spectra C ANG(II,I,J): Phase spectra C IMPLICIT REAL*8(A-H) IMPLICIT COMPLEX*16(O-Z) DIMENSION A(MJ,MJ,MJ1), E(MJ,MJ), C(0:20) DIMENSION P(0:NF,MJ,MJ), FNC(0:NF,MJ,MJ), COH(0:NF,MJ,MJ) DIMENSION AMP(0:NF,MJ,MJ), ANG(0:NF,MJ,MJ) DIMENSION ZA(0:200,10,10), ZB(10,10) DIMENSION FC(0:200), FS(0:200), WRK(10,10) C DO 30 I=1,L DO 30 J=1,L C(0) = 0.0D0 IF( I.EQ.J ) C(0) = 1.0D0 DO 10 II=1,M 10 C(II) = -A(I,J,II) C CALL FOURIE( C(0),M+1,NF+1,FC(0),FS(0) ) C DO 20 II=0,NF 20 ZA(II,I,J) = DCMPLX( FC(II),FS(II) ) 30 CONTINUE C DO 200 II=0,NF DO 40 I=1,L DO 40 J=1,L 40 ZB(I,J) = ZA(II,I,J) C CALL CINV( ZB,ZDET,L,10 ) C DO 60 I=1,L DO 60 J=1,L SUM = (0.0D0,0.0D0) DO 50 IJ=1,L 50 SUM = SUM + ZB(I,IJ)*E(IJ,J) 60 WRK(I,J) = SUM DO 80 I=1,L DO 80 J=1,L SUM = 0.0D0 DO 70 IJ=1,L 70 SUM = SUM + WRK(I,IJ)*DCONJG( ZB(J,IJ) ) 80 P(II,I,J) = SUM C C ... Amplitude and phase ... C DO 90 I=1,L-1 DO 90 J=I+1,L AMP(II,I,J) = DSQRT( DREAL(P(II,I,J))**2 + DIMAG(P(II,I,J))**2) ANG(II,I,J) = DATAN( DIMAG(P(II,I,J))/DREAL(P(II,I,J)) ) IF( DIMAG(P(II,I,J)).GT.0.0D0.AND.DREAL(P(II,I,J)).LT.0.0D0 ) * ANG(II,I,J) = ANG(II,I,J) + 3.1415926535D0 IF( DIMAG(P(II,I,J)).LT.0.0D0.AND.DREAL(P(II,I,J)).LT.0.0D0 ) * ANG(II,I,J) = ANG(II,I,J) - 3.1415926535D0 90 CONTINUE C C ... Simple coherency ... C DO 100 I=1,L-1 DO 100 J=I+1,L COH(II,I,J) = (DREAL(P(II,I,J))**2 + DIMAG(P(II,I,J))**2)/ * (P(II,I,I)*P(II,J,J)) 100 CONTINUE C C ... Power contribution ... C DO 120 I=1,L FSUM = 0.0D0 DO 110 J=1,L FSUM = FSUM + ZB(I,J)*DCONJG( ZB(I,J) )*E(J,J) 110 FNC(II,I,J) = FSUM 120 CONTINUE C 200 CONTINUE C RETURN E N D SUBROUTINE CINV( X,DET,M,MJ ) C C ... Inverse and determinant of a complex matrix X ... C C INPUTS: C X: M*M SQUARE MATRIX C M: DIMENSION OF X C MJ: ABSOLUTE DIMENSION OF X IN THE MAIN PROGRAM C C OUTPUTS: C X: INVERSE OF X C DET: DETERMINANT OF X C IMPLICIT COMPLEX*16 (A-H,O-Z) DIMENSION X(MJ,MJ), IND(100) C DET = 1.0D0 DO 60 L=1,M XMAX = 0.1D-10 IMAX = 0 DO 10 I=L,M IF( CDABS( X(I,L) ).GT.CDABS( XMAX ) ) THEN XMAX = X(I,L) IMAX = I END IF 10 CONTINUE IND(L) = IMAX IF( IMAX.NE.L ) THEN IF( IMAX.LE.0 ) THEN DET = 0.0D0 RETURN END IF C C ROW INTERCHANGE DO 20 J=1,M XTEMP = X(IMAX,J) X(IMAX,J) = X(L,J) 20 X(L,J) = XTEMP DET = -DET END IF DET = DET*XMAX X(L,L) = 1.0D0 DO 30 J=1,M 30 X(L,J) = X(L,J)/XMAX DO 50 I=1,M IF(I.NE.L) THEN XTEMP =X(I,L) X(I,L) = 0.0D0 DO 40 J=1,M 40 X(I,J) = X(I,J) - XTEMP*X(L,J) END IF 50 CONTINUE 60 CONTINUE C C COLUMN INTERCHANGE IF( M.GT.1 ) THEN DO 110 J=1,M-1 JJ = IND(M-J) IF( JJ.NE.M-J ) THEN DO 100 I=1,M XTEMP = X(I,JJ) X(I,JJ) = X(I,M-J) 100 X(I,M-J) = XTEMP END IF 110 CONTINUE END IF 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