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
      DIMENSION  AR(MJ,NMAX), 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  ...  PRINT CHANGING SPECTRUM  ...
C
      CALL  prtvsp( AR,SIG2,MJ,M,N,NOBS,VAR )
      STOP
      E N D
      SUBROUTINE  prtvsp( A,SIG2,MJ,M,N,NOBS,VAR )
C
C  ...  COMPUTE AND print 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), VAR(3000)
      NF = 200
C
      write(6,620)  m, n, nobs
      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  ...  print out LOG-SPECTRUM  ...
C
      DO 10 I=0,NF
   10 SP(I) = SP(I) + DLOG10( VAR(J*NOBS-NOBS/2) )
      write(6,600)  j
      write(6,610)  (sp(i),i=0,nf)
  100 CONTINUE
C
      RETURN
  600 format( 1H ,'J =',I3 )
  610 format( 1H ,(10F8.3) )
  620 format( 1H ,'M =',I3,3X,'N =',I5,3X,'NOBS =',I3 )
      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