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 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) 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 ) 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 PTCSPC( M,L,NF,MJ,P,FNC,AMP,ANG,COH ) C C ... Draws cross spectra and 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 AMP(II,I,J): Amplitude spectra C ANG(II,I,J): Phase spectra C COH(II,I,J): Coherency C IMPLICIT REAL*8(A-H,V-Z) IMPLICIT COMPLEX*16(O-U) CHARACTER TITLE*72, VNAME*8 DIMENSION P(0:NF,MJ,MJ), FNC(0:NF,MJ,MJ), FSP(0:200) DIMENSION AMP(0:NF,MJ,MJ), ANG(0:NF,MJ,MJ), COH(0:NF,MJ,MJ) DIMENSION VNAME(3), VALUE(3) COMMON /CMDATA/ TITLE C WX = 16.5/L -0.5 VNAME(1) = 'M = ' VNAME(2) = 'L = ' VALUE(1) = M VALUE(2) = L TITLE = 'CROSS SPECTRUM ANALYSIS' C CALL PLOTS CALL HEADER( 'PROGRAM 6.2: CROSS SPECRA',40,2,VNAME,VALUE ) CALL PLOT( 2.0,17.0-WX,-3 ) C DO 100 I=1,L IF( I.GE.2 ) CALL PLOT( -SNGL(WX+0.5)*(L-1),-SNGL(WX)-0.5,-3 ) C DO 100 J=1,L IF( J.GE.2 ) CALL PLOT( SNGL(WX)+0.5,0.0,-3 ) IPEN = 1 IF( I.EQ.J ) THEN DO 20 II=0,NF 20 FSP(II) = DLOG10( DREAL( P(II,I,I) ) ) CALL MAXMIN ( FSP(0),NF+1,YMIN,YMAX,DY ) CALL PLTBOX( WX,WX,0.0D0,DBLE(NF),10.0D0,YMIN,YMAX,DY,2 ) CALL PLOTY( FSP(0),NF+1,YMIN,YMAX,WX,WX,0,1 ) ELSE IF( I.LT.J ) THEN DO 30 II=0,NF 30 FSP(II) = DLOG10( AMP(II,I,J) ) CALL MAXMIN ( FSP(0),NF+1,YMIN,YMAX,DY ) CALL PLTBOX(WX,WX,0.D0,DBLE(NF),10.D0,YMIN,YMAX,DY,IPEN) CALL PLOTY( FSP(0),NF+1,YMIN,YMAX,WX,WX,0,1 ) ELSE YMAX = 3.1415926535D0 CALL PLTBOX(WX,WX,0.D0,DBLE(NF),10.D0,-YMAX,YMAX,YMAX,IPEN) CALL PLOTY( ANG(0,J,I),NF+1,-YMAX,YMAX,WX,WX,0,1 ) END IF END IF 100 CONTINUE C CALL PLOTI CALL HEADER( 'POWER SPECTRA AND COHERENCY',40,2,VNAME,VALUE ) CALL PLOT( 2.0,17.0-SNGL(WX),-3 ) C DO 200 I=1,L IF( I.GE.2 ) CALL PLOT( -SNGL(WX+0.5)*(L-I+1),-SNGL(WX)-0.5,-3 ) C DO 200 J=I,L IF( J.GE.2 ) CALL PLOT( SNGL(WX)+0.5,0.0,-3 ) IPEN = 1 IF( I.EQ.J ) THEN DO 110 II=0,NF 110 FSP(II) = DLOG10( DREAL( P(II,I,I) ) ) CALL MAXMIN ( FSP(0),NF+1,YMIN,YMAX,DY ) CALL PLTBOX( WX,WX,0.0D0,DBLE(NF),10.0D0,YMIN,YMAX,DY,2 ) CALL PLOTY( FSP(0),NF+1,YMIN,YMAX,WX,WX,0,1 ) ELSE CALL PLTBOX( WX,WX,0.D0,DBLE(NF),10.D0,0.D0,1.D0,0.5D0,IPEN) CALL PLOTY( COH(0,I,J),NF+1,0.0D0,1.0D0,WX,WX,0,1 ) END IF 200 CONTINUE C WX = 4.5D0 C DO 300 I=1,L IF( MOD(I,3).EQ.1 ) THEN CALL PLOTI CALL HEADER( 'PROGRAM 6.2: NOISE CONTRIBUTION ', * 40,2,VNAME,VALUE ) CALL SYMBOL( 2.5,17.5,0.25,'NOISE CONTRIBUTION',0.0,18 ) CALL SYMBOL(13.0,17.5,0.25,'RELATIVE NOISE CONTRIBUTION',0.,27) CALL PLOT( 2.5,17.0-SNGL(WX),-3 ) END IF IF( I.GE.2 ) CALL PLOT( -SNGL(WX)-3.0,-SNGL(WX)-1.0,-3 ) C CALL MAXMIN( FNC(0,I,L),NF+1,YMIN,YMAX,DY ) CALL AXISXY( 0.0D0,0.0D0,WX,WX,0.0D0,0.5D0,YMIN,YMAX, * 0.1D0,DY,0.2D0,1,1,1 ) DO 220 J=1,L 220 CALL PLOTY ( FNC(0,I,J),NF+1,YMIN,YMAX,WX,WX,IPOS,1 ) CALL PLOT( SNGL(WX)+3.0,0.0,-3 ) C CALL AXISXY( 0.0D0,0.0D0,2*WX,WX,0.0D0,0.5D0,0.0D0,1.0D0, * 0.1D0,0.2D0,0.2D0,1,1,1 ) DO 240 J=1,L DO 230 II=0,NF 230 FSP(II) = FNC(II,I,J)/FNC(II,I,L) 240 CALL PLOTY ( FSP(0),NF+1,0.0D0,1.0D0,2*WX,WX,IPOS,1 ) 300 CONTINUE C CALL PLOTE C 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 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 PLTBOX ( WX,WY,XMIN,XMAX,XSTEP,YMIN,YMAX,YSTEP,IPEN ) C C ... This subroutine draws box window ... C C Inputs: C WX,WY: width and height of the box window C XMIN,XMAX: C XSTEP: C YMIN,YMAX: C YSTEP: C IPEN: thickness or colour of the window C REAL*8 WX, WY, XMIN, XMAX, XSTEP, YMIN, YMAX, YSTEP C CALL NEWPEN( IPEN ) CALL PLOT( 0.0,0.0,3 ) CALL PLOT( SNGL(WX),0.0,2 ) CALL PLOT( SNGL(WX),SNGL(WY),2 ) CALL PLOT( 0.0,SNGL(WY),2 ) CALL PLOT( 0.0,0.0,2 ) C CALL NEWPEN(1) DX = XSTEP/(XMAX-XMIN) NX = 1.0D0/DX DX = DX*WX NX = ( XMAX-XMIN )/XSTEP DO 10 I=1,NX CALL PLOT( DX*I,0.0,3 ) 10 CALL PLOT( DX*I,0.1,2 ) DY = YSTEP/( YMAX-YMIN ) NY = 1.0D0/DY DY = DY*WY DO 20 I=1,NY CALL PLOT( 0.0,DY*I,3 ) 20 CALL PLOT( 0.1,DY*I,2 ) 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