home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / cpm / languags / fortran / spectral.lbr / SPECTRA2.FZR / SPECTRA2.FOR
Text File  |  1987-01-27  |  5KB  |  168 lines
  1. C    SPECTRA2.FOR - PERIODOGRAM SMOOTHING TO OBTAIN SPECTRUM
  2. C    ESTIMATE FROM OUTPUT OF SPECTRAL.FOR.  SPECTRUM ESTIMATE
  3. C    IS COMPUTED BY REPEATED SMOOTHING WITH MODIFIED DANIELL
  4. C    WEIGHTS.  INPUT IS:
  5. C
  6. C    NOBS    -  INTEGER     - LENGTH OF ORIGINAL TIME SERIES
  7. C    NK    -  INTEGER     - NUMBER OF SMOOTHING PASSES
  8. C    K    -  INTEGER ARRAY - ARRAY OF SMOOTHING PARAMETERS
  9. C
  10. C    THE PERIODOGRAM TO BE SMOOTHED IS READ IN BY DATIN.
  11. C
  12.     DIMENSION X(513),Y(513),K(10)
  13.     CHARACTER * 40 DF
  14.     DATA PI /3.141593/
  15.     WRITE(1,2050)
  16.  2050    FORMAT(/" SPECTRA2 - Spectrum Estimates by Smoothed Periodogram")
  17.  2001    WRITE(1) "Enter UPPERCASE Periodogram Input File Name: "
  18.     READ(1) DF
  19.     IF(IOREAD(5,2,0,DF)) GOTO 2001
  20.     CALL DATIN(X,NPGM,START,STEP,5)
  21.     NP2=(NPGM-1)*2
  22.     WRITE(1,7001)
  23.  7001    FORMAT(' Number of Observations in Original Time Series: '$
  24.         READ(1,7002) NOBS
  25.  7002    FORMAT(I0)
  26.     WRITE(1,7003)
  27.  7003    FORMAT(' Number of Modified Daniell Passes to Make: '$
  28.         READ(1,7002) NK
  29.     DO 7005 I=1,NK
  30.     WRITE(1,7004) I
  31.  7004    FORMAT('    Weight (1/2 length of moving average) for pass ',I2,': '$
  32.  7005   READ(1,7002) K(I)
  33.     CON=FLOAT(NOBS)/(8.0*PI)
  34.     DO 10 I=1,NPGM
  35.    10    X(I)=X(I)*CON
  36.     WRITE(1,7006)
  37.  7006    FORMAT(' Smoothing Pass - '$
  38.     DO 50 I=1,NK
  39.     WRITE(1,7007) I
  40.  7007    FORMAT(' ',I2$
  41.    50    CALL MODDAN(X,Y,NPGM,K(I),1.0)
  42.     WRITE(1,7008)
  43.  7008    FORMAT(1X)
  44.     START=0.0
  45.     STEP=2.0*PI/FLOAT(NP2)
  46.  2052    WRITE(1) "Enter UPPERCASE File for Spectrum Output: "
  47.     READ(1) DF
  48.     IF(IOWRIT(8,2,0,DF)) GOTO 2052
  49.     DF="SPECTRUM ESTIMATE"
  50.     CALL DATOUP(X,NPGM,START,STEP,8,DF)
  51.     I=IOCLOS(8)
  52.     STOP
  53.     END
  54. C
  55.     FUNCTION EXTEND(X,I,N,SYM)
  56. C
  57. C    RETURNS THE ITH TERM IN SERIES X, EXTENDED IF NECESSARY WITH
  58. C    EVEN OR ODD SYMMETRY ACCORDING TO SIGN OF SYM (+1 OR -1).
  59. C    IF SYM=0 THE EXTENDED VALUE IS ZERO.
  60. C
  61.     DIMENSION X(N)
  62.     IF(N.GT.1) GOTO 10
  63.     WRITE(1,1) N
  64.     1    FORMAT(' *** ERROR VALUE OF N IN EXTEND IS ',I5)
  65.     STOP
  66.    10    J=I
  67.     CON=1
  68.    20    IF(J.GE.1) GOTO 30
  69.     J=2-J
  70.     CON=CON*SYM
  71.    30    IF(J.LE.N) GOTO 40
  72.     J=2*N-J
  73.     CON=CON*SYM
  74.     GOTO 20
  75.    40    EXTEND=X(J)*CON
  76.     RETURN
  77.     END
  78. C
  79.     SUBROUTINE MODDAN(X,Y,N,K,SYM)
  80. C
  81. C    MODIFIED DANIELL SMOOTHING TO SERIES X.  COMPUTES SPECTRUM
  82. C    ESTIMATE FROM PERIODOGRAMS.  ASSUMES SERIES IS SYMMETRIC ABOUT
  83. C    BOTH END VALUES TO PROVIDE END VALUES OF OUTPUT.  PARAMETERS
  84. C    ARE:
  85. C
  86. C    X    -  REAL ARRAY  -  THE SERIES    -  RETURNS SMOOTHED
  87. C    Y    -  REAL ARRAY  -  SCRATCH    -  RETURNS SERIES
  88. C    K    -  INTEGER     -  HALF-LENGTH    -  (UNCHANGED)
  89. C                  OF MOVING AVERAGE
  90. C    SYM    -  REAL        -  +1 OR -1    -  (UNCHANGED)
  91. C                  SYMMETRY INDICATOR
  92. C    N    -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
  93. C
  94.     DIMENSION X(N),Y(N)
  95.     DO 10 I=1,N
  96.    10    Y(I)=X(I)
  97.     IF(K.LE.0) RETURN
  98.     LIM=K-1
  99.     CON=1.0/FLOAT(2*K)
  100.     DO 20 I=1,N
  101.     X(I)=Y(I)
  102.     IF(LIM.EQ.0) GOTO 20
  103.     DO 30 J=1,LIM
  104.    30    X(I)=X(I)+EXTEND(Y,I-J,N,SYM)+EXTEND(Y,I+J,N,SYM)
  105.    20    X(I)=(X(I)+(EXTEND(Y,I-K,N,SYM)+EXTEND(Y,I+K,N,SYM))*0.5)*CON
  106.     RETURN
  107.     END
  108. C
  109.     SUBROUTINE DATIN(X,N,START,STEP,M)
  110. C
  111. C    INPUT A SERIES OF VALUES IN RUNTIME FORMAT.  THE FIRST FOUR
  112. C    CARDS MUST BE (FREE FORMAT):
  113. C        1) A TITLE CARD (<72 COLS.)
  114. C        2) NUMBER OF CASES (I5)
  115. C        3) DATA FORMAT (<72 COLS.)
  116. C        4) START TIME VALUE AND STEP (2F10.5)
  117. C
  118. C    X    -  REAL ARRAY  -  RETURN THE SERIES
  119. C    N    -  INTEGER     -  RETURN SERIES LENGTH
  120. C    START    -  REAL           -  RETURN SERIES TIME VALUE AT 1ST POINT
  121. C    STEP    -  REAL        -  RETURN TIME INCREMENT BETWEEN POINTS
  122. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  123. C
  124.     DIMENSION X(513)
  125.     CHARACTER * 72 HEAD,FMT
  126.     WRITE(1,7001)
  127.  7001    FORMAT(' Reading data...')
  128.     READ(M,2) HEAD,N,FMT,START,STEP
  129.     2    FORMAT(A0/I5/A0/2F10.5)
  130.     WRITE(1,3) HEAD,N,FMT,START,STEP
  131.     3    FORMAT(5X,'The Input Data Header Reads:'/5X,A0/
  132.      *    5X,'Series Length Input is = ',I5/5X,'Format Of Input Data:'/
  133.      *    5X,A0/5X,'Time Origin = ',F11.5,'  Time Increment = ',F11.5)
  134.         READ(M,FMT) (X(I),I=1,N)
  135.     RETURN
  136.     END
  137. C
  138.     SUBROUTINE DATOUP(X,N,START,STEP,M,HEAD)
  139. C
  140. C    OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN.
  141. C    UNLIKE DATOUT (1) FREQ VALUES, (2) PERIOD FREQ,  AND
  142. C    (3) LOG OF SPECTRUM ARE OUTPUT WITH (4) SPECTRUM
  143. C    SO DATA MAY BE EDITED FOR PLOTTING SOFTWARE INPUT.
  144. C
  145. C    X    -  REAL ARRAY  -  THE SERIES
  146. C    N    -  INTEGER     -  SERIES LENGTH
  147. C    START    -  REAL           -  SERIES TIME VALUE AT 1ST POINT
  148. C    STEP    -  REAL        -  TIME INCREMENT BETWEEN POINTS
  149. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  150. C    HEAD    -  CHAR*40     -  HEADER TITLE
  151. C
  152.     DIMENSION X(N)
  153.     CHARACTER * 40 HEAD
  154.     DATA PI /3.141593/
  155.     WRITE(1,7001)
  156.  7001    FORMAT(' Writing data...')
  157.     WRITE(M,1) HEAD,N,START,STEP
  158.     1    FORMAT(1X,A0/1X,I5/' (28X,E15.8)'/1X,2F10.5)
  159.     Y=1024.0
  160.     DO 3 I=1,N
  161.     Z=ALOG(X(I))
  162.     WRITE(M,2) START,Y,Z,X(I)
  163.     2    FORMAT(1X,F7.5,1X,F9.4,1X,F9.5,1X,E15.8)
  164.     START=START+STEP
  165.     3    Y=2.0*PI/START
  166.     RETURN
  167.     END
  168.