home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / cpm / fortran / spectral.lbr / SPECTRAL.FZR / SPECTRAL.FOR
Encoding:
Text File  |  1987-01-28  |  6.9 KB  |  270 lines
  1. C    SPECTRAL.FOR - COMPUTATION OF THE DISCRETE FOURIER TRANSFORM
  2. C    AND THE PERIODOGRAM OF A TIME SERIES.  THE PERIODOGRAM IS
  3. C    COMPUTED AS THE SQUARED AMPLITUDE R(W)**2 RATHER THAN THE
  4. C    ACTUAL PERIODOGRAM (N/(8*PI))R(W)**2 FOR COMPUTATIONAL USE
  5. C    IN SUBSEQUENT PROCESSING.
  6. C
  7.     DIMENSION X(1024),Y(1024)
  8.     CHARACTER * 40 DF
  9.     WRITE(1,2050)
  10.  2050    FORMAT(/" SPECTRAL - Fourier Transform And Periodogram")
  11.  2001    WRITE(1) "Enter UPPERCASE Input File Name: "
  12.     READ(1) DF
  13.     IF(IOREAD(5,2,0,DF)) GOTO 2001
  14.     M=5
  15.     DO 10 I=1,1024
  16.     X(I)=0.0
  17.    10    Y(I)=0.0
  18.     CALL DATIN(X,N,START,STEP,M)
  19.     WRITE(1,7001) 
  20.  7001    FORMAT(' Order of Trend to Remove (0=Mean, 1=Linear): '$
  21.         READ(1,7002) K
  22.  7002    FORMAT(I0)
  23.     CALL DETRND(X,N,K)
  24.     WRITE(1,7003)
  25.  7003    FORMAT(' Proportion of Data to Taper: '$
  26.         READ(1,7004) P
  27.  7004    FORMAT(F0.0)
  28.     CALL TAPER(X,N,P)
  29.     NP2=1
  30.    20    NP2=NP2*2
  31.     IF(NP2.LT.N) GOTO 20
  32.     INV=0
  33.     CALL FT01A(X,Y,NP2,INV)
  34.     IF(INV.EQ.0) GOTO 30
  35.     WRITE(1) "*** ERROR IN SUBROUTINE FT01A ***"
  36.     GOTO 99
  37.    30    NPGM=NP2/2+1
  38.  2052    WRITE(1) "Enter UPPERCASE File for Transform Output: "
  39.     READ(1) DF
  40.     IF(IOWRIT(8,2,0,DF)) GOTO 2052
  41.     DF="REAL PART OF TRANSFORM"
  42.     CALL DATOUT(X,NPGM,0.0,1.0,8,DF)
  43.     DF="IMAGINARY PART OF TRANSFORM"
  44.     CALL DATOUT(Y,NPGM,0.0,1.0,8,DF)
  45.         INV=IOCLOS(8)
  46.     CON=(2.0*FLOAT(NP2)/FLOAT(N))**2
  47.     DO 40 I=1,NPGM
  48.    40    X(I)=(X(I)**2+Y(I)**2)*CON
  49.  2053    WRITE(1) "Enter UPPERCASE File for Periodogram Output: "
  50.     READ(1) DF
  51.     IF(IOWRIT(8,2,0,DF)) GOTO 2053
  52.     DF="SQUARED AMPLITUDE"
  53.     CALL DATOUT(X,NPGM,0.0,1.0,8,DF)
  54.         INV=IOCLOS(8)
  55.    99    STOP
  56.     END
  57. C
  58.     SUBROUTINE FT01A(XR,XI,N,INV)
  59. C
  60. C    THIS SUBROUTINE IMPLEMENTS THE SANDE-TUKEY RADIX-2 FAST FOURIER
  61. C    TRANSFORM.  EITHER THE DIRECT OR INVERSE TRANSFORM MAY BE COMP-
  62. C    UTED.  PARAMETERS ARE:
  63. C
  64. C    XR  -  REAL ARRAY  -  REAL PART OF THE  -  RETURNS REAL PART
  65. C                  SERIES           OF TRANSFORM
  66. C    XI  -  REAL ARRAY  -  IMAGINARY PART OF -  RETURNS IMAGINARY
  67. C                  THE SERIES       PART OF TRANSFORM
  68. C    N   -  INTEGER     -  LENGTH OF SERIES  -  (SAME)
  69. C    INV -  INTEGER     -  0 FOR DIRECT      -  RETURNS ERROR CODE
  70. C                  1 FOR INVERSE        (-1 IF ERROR)
  71. C
  72. C    DIRECT IS:
  73. C           N
  74. C    (1/N) SUM X(T)*EXP(-2*PI*I*(T-1)*(J-1)/N), J=1,...,N
  75. C             T=1
  76. C
  77. C    INVERSE IS:
  78. C           N
  79. C          SUM X(T)*EXP( 2*PI*I*(T-1)*(J-1)/N), J=1,...,N
  80. C             T=1
  81. C
  82.     DIMENSION XR(N),XI(N),UR(15),UI(15)
  83.     LOGICAL FIRST
  84.     DATA FIRST /.TRUE./
  85.     IF(.NOT.FIRST) GOTO 120
  86.     UR(1)=0.0
  87.     UI(1)=1.0
  88.     DO 110 I=2,15
  89.     UR(I)=SQRT((1.0+UR(I-1))/2.0)
  90.   110    UI(I)=UI(I-1)/(2.0*UR(I))
  91.     FIRST=.FALSE.
  92.   120    IF((N.GT.0).AND.(FLOAT(N).LE.2.0**16)) GOTO 130
  93.     INV=-1
  94.     RETURN
  95.   130    N0=1
  96.     II=0
  97.   140    N0=N0+N0
  98.     II=II+1
  99.     IF(N0.LT.N) GOTO 140
  100.     I1=N0/2
  101.     I3=1
  102.     I0=II
  103.     WRITE(1,7001)
  104.  7001    FORMAT(' FFTing data pass - '$
  105.     DO 260 I4=1,II
  106.     WRITE(1,7002) I4
  107.  7002    FORMAT(' ',I2$
  108.     DO 250 K=1,I1
  109.     WR=1.0
  110.     WI=0.0
  111.     KK=K-1
  112.     DO 230 I=1,I0
  113.     IF(KK.EQ.0) GOTO 240
  114.     IF(MOD(KK,2).EQ.0) GOTO 230
  115.     J0=I0-I
  116.     WS=WR*UR(J0)-WI*UI(J0)
  117.     WI=WR*UI(J0)+WI*UR(J0)
  118.     WR=WS
  119.   230    KK=KK/2
  120.   240    IF(INV.EQ.0) WI=-WI
  121.     L=K
  122.     DO 250 J=1,I3
  123.     L1=L+I1
  124.     ZR=XR(L)+XR(L1)
  125.     ZI=XI(L)+XI(L1)
  126.     Z=WR*(XR(L)-XR(L1))-WI*(XI(L)-XI(L1))
  127.     XI(L1)=WR*(XI(L)-XI(L1))+WI*(XR(L)-XR(L1))
  128.     XR(L1)=Z
  129.     XR(L)=ZR
  130.     XI(L)=ZI
  131.   250    L=L1+I1
  132.     I0=I0-1
  133.     I3=I3+I3
  134.   260    I1=I1/2
  135.     WRITE(1,7003)
  136.  7003    FORMAT('  *')
  137.     UM=1.0
  138.     IF(INV.EQ.0) UM=1.0/FLOAT(N0)
  139.     DO 310 J=1,N0
  140.     K=0
  141.     J1=J-1
  142.     DO 320 I=1,II
  143.     K=2*K+MOD(J1,2)
  144.   320    J1=J1/2
  145.     K=K+1
  146.     IF(K.LT.J) GOTO 310
  147.     ZR=XR(J)
  148.     ZI=XI(J)
  149.     XR(J)=XR(K)*UM
  150.     XI(J)=XI(K)*UM
  151.     XR(K)=ZR*UM
  152.     XI(K)=ZI*UM
  153.   310    CONTINUE
  154.     RETURN
  155.     END
  156. C
  157.     SUBROUTINE DETRND(X,N,K)
  158. C
  159. C    COMPUTES AND SUBTRACTS EITHER THE SERIES MEAN OR A LINEAR
  160. C    LEAST SQUARES TREND FROM THE SERIES DATA, E.G. POLYMONIAL
  161. C    OF DEGREE K
  162. C
  163. C    X  -  REAL ARRAY  -  THE INPUT SERIES -  RETURNS DETRENDED
  164. C    N  -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
  165. C    K  -  INTEGER     -  DEGREE OF POLYN  -  (UNCHANGED)
  166. C
  167.     DIMENSION X(N)
  168.     WRITE(1,7001)
  169.  7001    FORMAT(' Detrending data...')
  170.     SUMX=0.0
  171.     DO 20 I=1,N
  172.    20    SUMX=SUMX+X(I)
  173.     XBAR=SUMX/FLOAT(N)
  174.     WRITE(1,7002) XBAR
  175.  7002    FORMAT(5X,'Mean Removed = ',F13.5)
  176.     DO 30 I=1,N
  177.    30    X(I)=X(I)-XBAR
  178.     IF(K.LE.0) RETURN
  179.     TBAR=FLOAT(N+1)/2.0
  180.     SUMTT=(FLOAT(N)*(FLOAT(N)**2-1))/12.0
  181.     SUMTX=0.0
  182.     DO 40 I=1,N
  183.    40    SUMTX=SUMTX+X(I)*(FLOAT(I)-TBAR)
  184.     BETA=SUMTX/SUMTT
  185.     WRITE(1,7003) BETA
  186.  7003    FORMAT(5X,'Slope Removed = ',F13.5)
  187.     DO 50 I=1,N
  188.    50    X(I)=X(I)-BETA*(FLOAT(I)-TBAR)
  189.     RETURN
  190.     END
  191. C
  192.     SUBROUTINE TAPER(X,N,P)
  193. C
  194. C    APPLIES SPLIT-COSINE-BELL TAPERING TO THE TIME SERIES.  THE
  195. C    ARGUMENT P IS THE TOTAL PROPORTION OF THE SERIES TO TAPER.
  196. C    TUKEY SUGGESTS .1 TO .2, .25 IS THE DEFAULT IN MANY PACKAGES.
  197. C
  198. C    X  -  REAL ARRAY  -  THE INPUT SERIES -  RETURNS TAPERED
  199. C    N  -  INTEGER     -  SERIES LENGTH    -  (UNCHANGED)
  200. C    P  -  REAL        -  PROP. TO  TAPER  -  (UNCHANGED)
  201. C
  202.     DATA PI /3.141593/
  203.     DIMENSION X(N)
  204.     IF((P.LE.0.0).OR.(P.GE.1.0)) RETURN
  205.     WRITE(1,7001)
  206.  7001    FORMAT(' Tapering data...')
  207.     M=INT(P*FLOAT(N)+0.5)/2
  208.     DO 10 I=1,M
  209.     WEIGHT=0.5-0.5*COS(PI*(FLOAT(I)-0.5)/FLOAT(M))
  210.     X(I)=X(I)*WEIGHT
  211.    10    X(N+1-I)=X(N+1-I)*WEIGHT
  212.     RETURN
  213.     END
  214. C
  215.     SUBROUTINE DATIN(X,N,START,STEP,M)
  216. C
  217. C    INPUT A SERIES OF VALUES IN RUNTIME FORMAT.  THE FIRST FOUR
  218. C    CARDS MUST BE (FREE FORMAT):
  219. C        1) A TITLE CARD (<72 COLS.)
  220. C        2) NUMBER OF CASES (I5)
  221. C        3) DATA FORMAT (<72 COLS.)
  222. C        4) START TIME VALUE AND STEP (2F10.5)
  223. C
  224. C    X    -  REAL ARRAY  -  RETURN THE SERIES
  225. C    N    -  INTEGER     -  RETURN SERIES LENGTH
  226. C    START    -  REAL           -  RETURN SERIES TIME VALUE AT 1ST POINT
  227. C    STEP    -  REAL        -  RETURN TIME INCREMENT BETWEEN POINTS
  228. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  229. C
  230.     DIMENSION X(1024)
  231.     CHARACTER * 72 HEAD,FMT
  232.     WRITE(1,7001)
  233.  7001    FORMAT(' Reading data...')
  234.     READ(M,2) HEAD,N,FMT,START,STEP
  235.     2    FORMAT(A0/I5/A0/2F10.5)
  236.     WRITE(1,3) HEAD,N,FMT,START,STEP
  237.     3    FORMAT(5X,'The Input Data Header Reads:'/5X,A0/
  238.      *    5X,'Series Length Input is = ',I5/5X,'Format Of Input Data:'/
  239.      *    5X,A0/5X,'Time Origin = ',F11.5,'  Time Increment = ',F11.5)
  240.         READ(M,FMT) (X(I),I=1,N)
  241.     RETURN
  242.     END
  243. C
  244.     SUBROUTINE DATOUT(X,N,START,STEP,M,HEAD)
  245. C
  246. C    OUTPUT A SERIES OF VALUES IN FORMAT REQUIRED FOR DATIN.
  247. C
  248. C    X    -  REAL ARRAY  -  THE SERIES
  249. C    N    -  INTEGER     -  SERIES LENGTH
  250. C    START    -  REAL           -  SERIES TIME VALUE AT 1ST POINT
  251. C    STEP    -  REAL        -  TIME INCREMENT BETWEEN POINTS
  252. C    M    -  INTEGER     -  UNIT NUMBER (UNCHANGED)
  253. C    HEAD    -  CHAR*40     -  HEADER TITLE
  254. C
  255.     DIMENSION X(N)
  256.     CHARACTER * 40 HEAD
  257.     WRITE(1,7001)
  258.  7001    FORMAT(' Writing data...')
  259.     WRITE(M,1) HEAD,N,START,STEP
  260.     1    FORMAT(1X,A0/1X,I5/' (5E15.8)'/1X,2F10.5)
  261.     WRITE(M,2) (X(I),I=1,N)
  262.     2    FORMAT(1X,5E15.8)
  263.     RETURN
  264.     END
  265. D,FMT
  266.     WRITE(1,7001)
  267.  7001    FORMAT(' Reading data...')
  268.     READ(M,2) HEAD,N,FMT,START,STEP
  269.     
  270.