home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / cpm / fortran / spectral.lbr / SPECTRA3.FZR / SPECTRA3.FOR
Encoding:
Text File  |  1987-01-28  |  6.9 KB  |  232 lines
  1. C    SPECTRA3.FOR - CROSS PERIODOGRAM SMOOTHING FROM TWO SERIES
  2. C    TO ESTIMATE COHERENCY AND PHASE SPECTRUM.  COHERENCY
  3. C    GIVES SHARED HARMONICS AT SAME FREQUENCY AND PHASE SPECTRUM
  4. C    GIVES THE DIRECTION OF THIS ASSOCIATION.  SPECTRUM ESTIMATE
  5. C    IS COMPUTED BY REPEATED SMOOTHING WITH MODIFIED DANIELL
  6. C    WEIGHTS.  INPUT IS:
  7. C
  8. C    NOBS    -  INTEGER     - LENGTH OF ORIGINAL TIME SERIES
  9. C    NK    -  INTEGER     - NUMBER OF SMOOTHING PASSES
  10. C    K    -  INTEGER ARRAY - ARRAY OF SMOOTHING PARAMETERS
  11. C
  12. C    THE TRANSFORMS OF BOTH SERIES TO BE SMOOTHED AND THE 
  13. C    SPECTRUMS OF EACH ARE READ IN BY DATIN.  THE SMOOTHING OF
  14. C    BOTH SPECTRUMS (K,NK) MUST HAVE BEEN SAME FOR BOTH SERIES
  15. C    AND IDENTICAL TO THAT SPECIFIED HERE.
  16. C
  17.     DIMENSION TR1(513),TI1(513),TR2(513),TI2(513),K(10)
  18.     CHARACTER * 40 DF
  19.     DATA PI /3.141593/
  20.     WRITE(1,2050)
  21.  2050    FORMAT(/" SPECTRA3 - Cross-Spectrum by Smoothed Cross-Periodogram")
  22.  2001    WRITE(1) "Enter UPPERCASE Transform Input File Name for 1st File: "
  23.     READ(1) DF
  24.     IF(IOREAD(5,2,0,DF)) GOTO 2001
  25.     CALL DATIN(TR1,NPGM1,START,STEP,5)
  26.     CALL DATIN(TI1,NPGM1,START,STEP,5)
  27.     I=IOCLOS(5)
  28.  2002    WRITE(1) "Enter UPPERCASE Transform Input File Name for 2nd File: "
  29.     READ(1) DF
  30.     IF(IOREAD(5,2,0,DF)) GOTO 2002
  31.     CALL DATIN(TR2,NPGM2,START,STEP,5)
  32.     CALL DATIN(TI2,NPGM2,START,STEP,5)
  33.     I=IOCLOS(5)
  34.     IF(NPGM1.EQ.NPGM2) GOTO 20
  35.     WRITE(1,2) NPGM1,NPGM2
  36.     2    FORMAT(' *** ERROR - TRANSFORM LENGTHS ',2I5,' DIFFER.')
  37.     STOP
  38.    20    N=NPGM1
  39.     NP2=(N-1)*2
  40.     WRITE(1,7001)
  41.  7001    FORMAT(' Number of Observations in Original Time Series: '$
  42.         READ(1,7002) NOBS
  43.  7002    FORMAT(I0)
  44.     CON=FLOAT(NP2)**2/(PI*FLOAT(NOBS)*2.0)
  45.     WRITE(1,7003)
  46.  7003    FORMAT(' Number of Modified Daniell Passes to Make: '$
  47.         READ(1,7002) NK
  48.     DO 7005 I=1,NK
  49.     WRITE(1,7004) I
  50.  7004    FORMAT('    Weight (1/2 length of moving average) for pass ',I2,': '$
  51.  7005   READ(1,7002) K(I)
  52.     DO 30 I=1,N
  53.     CR=TR1(I)*TR2(I)+TI1(I)*TI2(I)
  54.     CI=TI1(I)*TR2(I)-TR1(I)*TI2(I)
  55.     TR1(I)=CR*CON
  56.    30    TI1(I)=CI*CON
  57.     WRITE(1,7006)
  58.  7006    FORMAT(' Smoothing Pass - '$
  59.     DO 40 I=1,NK
  60.     WRITE(1,7007) I
  61.  7007    FORMAT(' ',I2$
  62.     CALL MODDAN(TR1,TR2,N,K(I),1.0)
  63.    40    CALL MODDAN(TI1,TI2,N,K(I),-1.0)
  64.     WRITE(1,7008)
  65.  7008    FORMAT(1X)
  66.     CALL POLAR(TR1,TI1,N)
  67.  2003    WRITE(1) "Enter UPPERCASE Spectrum Input File Name for 1st File: "
  68.     READ(1) DF
  69.     IF(IOREAD(5,2,0,DF)) GOTO 2003
  70.     CALL DATIN(TR2,N1,START,STEP,5)
  71.     I=IOCLOS(5)
  72.  2004    WRITE(1) "Enter UPPERCASE Spectrum Input File Name for 2nd File: "
  73.     READ(1) DF
  74.     IF(IOREAD(5,2,0,DF)) GOTO 2003
  75.     CALL DATIN(TI2,N2,START,STEP,5)
  76.     I=IOCLOS(5)
  77.     IF((N1.EQ.N).AND.(N2.EQ.N)) GOTO 60
  78.     WRITE(1,3) N,N1,N2
  79.     3    FORMAT(' *** ERROR - TRANSFORM LENGTH ',I5,' NOT EQUAL TO '
  80.      *  'SPECTRUM LENGTHS ',2I5)
  81.     STOP
  82.    60    DO 50 I=1,N
  83.    50    TR1(I)=TR1(I)/SQRT(TR2(I)*TI2(I))
  84.     START=0.0
  85.     STEP=PI/FLOAT(N-1)
  86.  2052    WRITE(1) "Enter UPPERCASE File for Cross-Spectrum Output: "
  87.     READ(1) DF
  88.     IF(IOWRIT(8,2,0,DF)) GOTO 2052
  89.     DF="CROSS-SPECTRUM ESTIMATE"
  90.     CALL DATOUC(TR1,TI1,N,START,STEP,8,DF)
  91.     I=IOCLOS(8)
  92.     STOP
  93.     END
  94. C
  95.     FUNCTION EXTEND(X,I,N,SYM)
  96. C
  97. C    RETURNS THE ITH TERM IN SERIES X, EXTENDED IF NECESSARY WITH
  98. C    EVEN OR ODD SYMMETRY ACCORDING TO SIGN OF SYM (+1 OR -1).
  99. C    IF SYM=0 THE EXTENDED VALUE IS ZERO.
  100. C
  101.     DIMENSION X(N)
  102.     IF(N.GT.1) GOTO 10
  103.     WRITE(1,1) N
  104.     1    FORMAT(' *** ERROR VALUE OF N IN EXTEND IS ',I5)
  105.     STOP
  106.    10    J=I
  107.     CON=1
  108.    20    IF(J.GE.1) GOTO 30
  109.     J=2-J
  110.     CON=CON*SYM
  111.    30    IF(J.LE.N) GOTO 40
  112.     J=2*N-J
  113.     CON=CON*SYM
  114.     GOTO 20
  115.    40    EXTEND=X(J)*CON
  116.     RETURN
  117.     END
  118. C
  119.     SUBROUTINE MODDAN(X,Y,N,K,SYM)
  120. C
  121. C    MODIFIED DANIELL SMOOTHING TO SERIES X.  COMPUTES SPECTRUM
  122. C    ESTIMATE FROM PERIODOGRAMS.  ASSUMES SERIES IS SYMMETRIC ABOUT
  123. C    BOTH END VALUES TO PROVIDE END VALUES OF OUTPUT.  PARAMETERS
  124. C    ARE:
  125. C
  126. C    X    -  REAL ARRAY  -  THE SERIES    -  RETURNS SMOOTHED
  127. C    Y    -  REAL ARRAY  -  SCRATCH    -  RETURNS SERIES
  128. C    K    -  INTEGER     -  HALF-LENGTH    -  (UNCHANGED)
  129. C                  OF MOVING AVERAGE
  130. C    SYM    -  REAL        -  +1 OR -1    -  (UNCHANGED)
  131. C                  SYMMETRY INDICATOR
  132. C    N    -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
  133. C
  134.     DIMENSION X(N),Y(N)
  135.     DO 10 I=1,N
  136.    10    Y(I)=X(I)
  137.     IF(K.LE.0) RETURN
  138.     LIM=K-1
  139.     CON=1.0/FLOAT(2*K)
  140.     DO 20 I=1,N
  141.     X(I)=Y(I)
  142.     IF(LIM.EQ.0) GOTO 20
  143.     DO 30 J=1,LIM
  144.    30    X(I)=X(I)+EXTEND(Y,I-J,N,SYM)+EXTEND(Y,I+J,N,SYM)
  145.    20    X(I)=(X(I)+(EXTEND(Y,I-K,N,SYM)+EXTEND(Y,I+K,N,SYM))*0.5)*CON
  146.     RETURN
  147.     END
  148. C
  149.     SUBROUTINE POLAR(X,Y,N)
  150. C
  151. C    CONVERTS REAL AND IMAGINARY PARTS OF CROSS SPECTRUM TO
  152. C    MAGNITUDES AND PHASES.  CONVERSION IS DONE IN PLACE GIVEN:
  153. C
  154. C    X    -  REAL ARRAY  -  REAL PART    -  RETURNS MAGNITUDE
  155. C    Y    -  REAL ARRAY  -  IMAG PART    -  RETURNS PHASE
  156. C    N    -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
  157. C
  158.     DIMENSION X(N),Y(N)
  159.     DO 10 I=1,N
  160.     R=SQRT(X(I)**2+Y(I)**2)
  161.     PHI=ATAN2(Y(I),X(I))
  162.     X(I)=R
  163.    10    Y(I)=PHI
  164.     RETURN
  165.     END
  166. C
  167.     SUBROUTINE DATIN(X,N,START,STEP,M)
  168. C
  169. C    INPUT A SERIES OF VALUES IN RUNTIME FORMAT.  THE FIRST FOUR
  170. C    CARDS MUST BE (FREE FORMAT):
  171. C        1) A TITLE CARD (<72 COLS.)
  172. C        2) NUMBER OF CASES (I5)
  173. C        3) DATA FORMAT (<72 COLS.)
  174. C        4) START TIME VALUE AND STEP (2F10.5)
  175. C
  176. C    X    -  REAL ARRAY  -  RETURN THE SERIES
  177. C    N    -  INTEGER     -  RETURN SERIES LENGTH
  178. C    START    -  REAL           -  RETURN SERIES TIME VALUE AT 1ST POINT
  179. C    STEP    -  REAL        -  RETURN TIME INCREMENT BETWEEN POINTS
  180. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  181. C
  182.     DIMENSION X(513)
  183.     CHARACTER * 72 HEAD,FMT
  184.     WRITE(1,7001)
  185.  7001    FORMAT(' Reading data...')
  186.     READ(M,2) HEAD,N,FMT,START,STEP
  187.     2    FORMAT(A0/I5/A0/2F10.5)
  188.     WRITE(1,3) HEAD,N,FMT,START,STEP
  189.     3    FORMAT(5X,'The Input Data Header Reads:'/5X,A0/
  190.      *    5X,'Series Length Input is = ',I5/5X,'Format Of Input Data:'/
  191.      *    5X,A0/5X,'Time Origin = ',F11.5,'  Time Increment = ',F11.5)
  192.         READ(M,FMT) (X(I),I=1,N)
  193.     RETURN
  194.     END
  195. C
  196.     SUBROUTINE DATOUC(X,Y,N,START,STEP,M,HEAD)
  197. C
  198. C    OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN.
  199. C    UNLIKE DATOUT (1) FREQ VALUES, (2) PERIOD FREQ, (3) COHERENCY
  200. C    (4) SQUARED COHERENCY, (5) PHASE SPECTRUM AND (6) TIME LEAD
  201. C       AT FREQUENCEY ARE OUTPUT IN A FORMAT SO DATA MAY BE
  202. C    EDITED FOR PLOTTING SOFTWARE INPUT.
  203. C
  204. C    X    -  REAL ARRAY  -  COHERENCY
  205. C    X    -  REAL ARRAY  -  PHASE SPECTRUM
  206. C    N    -  INTEGER     -  SERIES LENGTH
  207. C    START    -  REAL           -  SERIES TIME VALUE AT 1ST POINT
  208. C    STEP    -  REAL        -  TIME INCREMENT BETWEEN POINTS
  209. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  210. C    HEAD    -  CHAR*40     -  HEADER TITLE
  211. C
  212.     DIMENSION X(N),Y(N)
  213.     CHARACTER * 40 HEAD
  214.     DATA PI /3.141593/
  215.     WRITE(1,7001)
  216.  7001    FORMAT(' Writing data...')
  217.     WRITE(M,1) HEAD,N,START,STEP
  218.     1    FORMAT(1X,A0/1X,I5/' *CANNOT BE REPROCESSED*'/1X,2F10.5)
  219.     Z=1024.0
  220.     DO 3 I=1,N
  221.         X2=X(I)**2
  222.         T=0.0
  223.         IF(START.GT.0.0) T=Y(I)/START
  224.     WRITE(M,2) START,Z,X(I),X2,Y(I),T
  225.     2    FORMAT(1X,F7.5,1X,F9.4,1X,F7.5,1X,F7.5,1X,F9.5,1X,F10.4)
  226.     START=START+STEP
  227.     3    Z=2.0*PI/START
  228.     RETURN
  229.     END
  230. D FOR DATIN.
  231. C    UNLIKE DATOUT (1) FREQ VALUES, (2) PERIOD FREQ, (3) COHERENCY
  232. C    (4) SQUARED COHERENCY, (5) PHASE SPEC