home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / s / spectr20.zip / SPECTRUM.FOR < prev   
Text File  |  1992-04-22  |  66KB  |  2,514 lines
  1. *    SPECTRUM.FOR  V 2.0
  2.  
  3. *    Routine to compute Amplitude, Auto and Cross-Spectral
  4. *    Density Functions via Fast Fourier Transforms
  5.  
  6. *    David E. Hess
  7. *    Fluid Flow Group - Process Measurements Division
  8. *    Chemical Science and Technology Laboratory
  9. *    National Institute of Standards and Technology
  10. *    April 15, 1992
  11.  
  12. *    References: 1. Bendat, J.S. and Piersol, A.G. 1986 Random Data,
  13. *                   2nd Ed., Wiley, New York.
  14. *            2. Harris, F.J. 1978 On the use of windows for
  15. *               harmonic analysis with the discrete fourier
  16. *               transform, Proceedings of the IEEE, 66, 51-83.
  17. *            3. Nuttall, A.H. 1981 Some windows with very good
  18. *               sidelobe behavior, IEEE Transactions on Acoustics,
  19. *               Speech and Signal Processing, ASSP-29, 84-91.
  20. *                4. Stearns, S.D. 1975 Digital Signal Analysis,
  21. *                   Hayden Book Co., Rochelle Park, New Jersey.
  22. *                5. Sirivat, A. 1985 Statistical Package for the
  23. *               Analysis of Data (SPAD), private communication. 
  24.  
  25. *       IFMAX  : maximum # of files that may be processed.
  26. *       JFMAX  : maximum # of files read from the command line.
  27. *       NMAX   : maximum # points per record.
  28. *       NSMAX  : maximum # of filter sections.
  29. *          NUMCON : # of coefficients; order of polynomial + 1
  30. *       EPS    : smallest allowable number in output data.
  31.  
  32. *           Phase Adjust Code parameters:
  33.  
  34. *       NPTS  : # pts before current point to fit
  35. *       NSD   : # points in moving average of s.d. (odd # only)
  36. *       SDLIM : quit when s.d. theta .GT. SDLIM
  37.  
  38. *    This code is compiled using Microsoft Fortran V 5.1.  For full
  39. *    optimization, at least 603K of RAM must be available when
  40. *    compiled.  For additional information refer to documentation.
  41.  
  42. *    This program reads in a data file consisting of integers
  43. *    created by an A/D sampling routine and transforms the data
  44. *    to voltages or to some other variable.  This other variable
  45. *    is computed using up to a fourth order polynomial; the
  46. *    coefficients must be provided in either a sequential formatted
  47. *    (Ascii) or a sequential binary file.  Alternatively, the routine
  48. *    can read REAL*4 data which has already been transformed by some
  49. *    other program.  The data can then be tapered in order to
  50. *    reduce side-lobe leakage.  The tapering may employ a 50% overlap
  51. *    and the user has a choice of six different window functions.
  52. *    Then, the spectra are computed; the user may choose from :
  53. *    amplitude, power or cross spectra.  The user may process
  54. *    up to IFMAX files; the parameters chosen for the first file
  55. *    are then applied to each of the IFMAX files.  Each record may
  56. *    contain no more than NMAX points.  Note that in the case of two
  57. *    channels of data per record, each channel may consist of up to
  58. *    NMAX/2 points.  The speed of the routine is limited by the large
  59. *    number of I/O operations that must take place.  Reduce # of points
  60. *    per record and/or the # of records per file to increase the speed
  61. *    of the routine.  Alternatively, specify in the configuration file
  62. *    a RAM disk for storage of temporary files -- this also increases
  63. *    speed.  The output is stored to a .PRN file with the same name
  64. *    as the input data file in order that it may be imported into
  65. *    LOTUS and plotted.  In order to avoid taking the logarithm of
  66. *    zero, any number in the output data less than EPS is set equal
  67. *    to EPS.  Then, the logarithm of EPS will be a finite number.
  68. *    The order of the recursive digital filters is twice the number
  69. *    of cascaded sections desired; the number of filter sections
  70. *    must be less than or equal to NSMAX.  The routine will read
  71. *    up to JFMAX filenames from the command line.
  72.  
  73.     INCLUDE        'SPECTRUM.FI'    ! Interface statements
  74.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  75.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  76.  
  77. *    Allocatable arrays may not be placed within COMMON blocks
  78. *    and may not be dummy arguments within functions or subroutines.
  79.  
  80.     INTEGER*2     IFIL[ALLOCATABLE](:),GAIN(0:7)
  81.     INTEGER*2     NDATA[ALLOCATABLE,HUGE](:)
  82.     REAL*4        WIN[ALLOCATABLE,HUGE](:)
  83.     REAL*4        FR[ALLOCATABLE,HUGE](:),FI[ALLOCATABLE,HUGE](:)
  84.     REAL*4        FR2[ALLOCATABLE,HUGE](:),FI2[ALLOCATABLE,HUGE](:)
  85.     REAL*4         CONST [ALLOCATABLE](:,:),B (NUMCON)
  86.     REAL*4        CONST2[ALLOCATABLE](:,:),B2(NUMCON)
  87.     DATA         TNAM,WNAM,PNAM,HNAM,SNAM,DNAM
  88.      +                  /'Tukey','Welch','Parzen','Hanning',
  89.      +            'Min Side Lobe','Max Decay'/
  90.  
  91. *    Initialize.
  92.  
  93.     PI=2.0*ASIN(1.0)
  94.     PI2=2.0*PI
  95.     COH=.FALSE.
  96.     NOCOH=.TRUE.
  97.     ROOT1=.TRUE.
  98.     ROOT2=.TRUE.
  99.     SUFX=.FALSE.
  100.  
  101. *    Code to read data from configuration file.
  102.  
  103.     INQUIRE (FILE='SPECTRUM.CFG',EXIST=EXISTS)
  104.     IF (.NOT. EXISTS) GO TO 5
  105.  
  106. *    File exists. Get the info.
  107.  
  108.     CALL RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)
  109.  
  110. *    Code to perform command line parsing.
  111.  
  112. 5    NUMARGS=NARGS()
  113.     IF (NUMARGS .EQ. 1) THEN
  114.       NUMARGS=NUMARGS-1
  115.       GO TO 10
  116.     ELSE
  117.       NUMARGS=NUMARGS-1
  118.       IF (NUMARGS .GT. JFMAX) THEN
  119.         WRITE (*,'(/1X,A/1X,I1,A/1X,A)')
  120.      +      'The routine is currently set to recognize up to ',
  121.      +      JFMAX,' filenames on the command line and will',
  122.      +      'proceed with these filenames only.'
  123.         NUMARGS=JFMAX
  124.       ENDIF
  125.       DO I=1,NUMARGS
  126.         CALL GETARG (I,BUFFER,STATUS)
  127.         DO J=4,1,-1
  128.           FIRST=BUFFER(J:J)
  129.           IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
  130.             IHOLD=ICHAR(FIRST)-32
  131.             FIRST=CHAR(IHOLD)
  132.             BUFFER(J:J)=FIRST
  133.           ENDIF
  134.         ENDDO
  135.         IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
  136.           WRITE (*,'(/1X,A,I1,A/1X,A,A,A/)') 
  137.      +       'Filename ',I,' on the command line began with',
  138.      +       'the nonalphabetic character ',FIRST,'.'
  139.           STOP 'Re-enter the command line correctly'
  140.         ENDIF
  141.         IF (.NOT. ROOT1) THEN
  142.           ISL=INDEX (PTH,'\',.TRUE.)
  143.           FLNM(I)=PTH(:ISL) // BUFFER // '.DAT'
  144.         ELSE
  145.           FLNM(I)=BUFFER // '.DAT'
  146.         ENDIF
  147.         INQUIRE (FILE=FLNM(I),EXIST=EXISTS)
  148.         IF (.NOT. EXISTS) THEN
  149.           IPD=INDEX (FLNM(I),'.',.TRUE.)
  150.           WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
  151.      +                                FLNM(I)(:IPD+3)
  152.           STOP 'Re-enter the command line correctly'
  153.         ENDIF
  154.       ENDDO
  155.       CONSEC=.TRUE.
  156.       ARB=.FALSE.
  157.       IFSTRT=1
  158.       IFSTP=NUMARGS
  159.       IFTOT=NUMARGS
  160.       GO TO 40
  161.     ENDIF
  162.  
  163. *    This controls the values that are fed to IF.
  164.  
  165. 10    WRITE (*,'(/1X,A,A\)') 'Process files (C)onsecutively',
  166.      +                         ' or in an (A)rbitrary order : '
  167.     READ (*,'(A)') INP
  168.     WRITE (*,'( )')
  169.     IF (INP .EQ. 'c') INP = 'C'
  170.     IF (INP .EQ. 'a') INP = 'A'
  171.     CONSEC=(INP .EQ. 'C')
  172.     ARB=(INP .EQ. 'A')
  173.     IF (.NOT. CONSEC .AND. .NOT. ARB) GO TO 10
  174.  
  175. *    These statements will allow IF to be set equal to values
  176. *    stored in the IFIL array in the loop below.
  177.  
  178.     IF (ARB) THEN
  179. 20      WRITE (*,'(1X,A\)') 'Enter number of files to process : '
  180.       READ (*,*) IFTOT
  181.       WRITE (*,'( )')
  182.       IF (IFTOT .GT. IFMAX) THEN
  183.         WRITE (*,'(1X,''# files must be .LE. '',I3,''.''//)') IFMAX
  184.         GO TO 20
  185.       ENDIF
  186.       ALLOCATE (IFIL(IFTOT), STAT=IERR)  ! Allocate space for IFIL
  187.       IF (IERR .NE. 0)
  188.      +       STOP 'Problem allocating storage for IFIL.  Aborting ...'
  189.       IFILMIN=IFMAX
  190.       IFILMAX=1
  191.       DO I=1,IFTOT
  192.         WRITE (*,'(1X,''Enter last digits of file '',I2,'' : ''\)') I
  193.         READ (*,*) IFIL(I)
  194.         IFILMIN=MIN0(IFIL(I),IFILMIN)  ! Minimum of IFIL array
  195.         IFILMAX=MAX0(IFIL(I),IFILMAX)  ! Maximum of IFIL array
  196.       ENDDO
  197.       WRITE (*,'( )')
  198.       IFSTRT=1
  199.       IFSTP=IFTOT
  200.       WRITE (*,'(1X,A)') 'Files are : '
  201.       WRITE (*,'(13X,10(I2,2X))') (IFIL(I),I=1,IFTOT)
  202.       WRITE (*,'( )')
  203.     ELSE
  204. 25      WRITE (*,'(1X,A\)') 'Enter starting & ending data file # : '
  205.       READ (*,*) IFSTRT,IFSTP
  206.       IF (IFSTP .LT. IFSTRT) THEN
  207.         WRITE (*,'(1X,A)')
  208.      +             'Ending data file # >= Starting data file #.'
  209.         GO TO 25
  210.       ENDIF
  211.       WRITE (*,'( )')
  212.       IFTOT=IFSTP-IFSTRT+1
  213.     ENDIF
  214.  
  215. *    Get the first letter.
  216.  
  217. 30    WRITE (*,'(1X,A\)') 'Enter first letter of data file : '
  218.     READ (*,'(A)') FIRST
  219.     IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
  220.       IHOLD=ICHAR(FIRST)-32
  221.       FIRST=CHAR(IHOLD)
  222.     ENDIF
  223.     IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
  224.       WRITE (*,'(1X,A/)') 'Enter an alphabetic character (A-Z).'
  225.       GO TO 30
  226.     ENDIF
  227.  
  228. *    Assign coefficient data file names.
  229.  
  230. 40    IF (.NOT. ROOT1) THEN
  231.       ISL=INDEX (PTH,'\',.TRUE.)
  232.       XNAM =PTH(:ISL) // FIRST // 'CON.DAT'
  233.       YNAM =PTH(:ISL) // FIRST // 'CON.PRN'
  234.       XNAM2=PTH(:ISL) // FIRST // 'CON2.DAT'
  235.       YNAM2=PTH(:ISL) // FIRST // 'CON2.PRN'
  236.     ELSE
  237.       XNAM=FIRST//'CON.DAT'
  238.       YNAM=FIRST//'CON.PRN'
  239.       XNAM2=FIRST//'CON2.DAT'
  240.       YNAM2=FIRST//'CON2.PRN'
  241.     ENDIF
  242.  
  243. *    Transform to voltages, other quantity or read Real data ?
  244.  
  245. 50    WRITE (*,'(/1X,A/1X,A\)')
  246.      +    'Transform into (V)oltages or into (O)ther quantity',
  247.      +    'or read (R)eal data which is already transformed : '
  248.     READ (*,'(A)') INP
  249.     IF (INP .EQ. 'v') INP = 'V'
  250.     IF (INP .EQ. 'o') INP = 'O'
  251.     IF (INP .EQ. 'r') INP = 'R'
  252.     VOLT=(INP .EQ. 'V')
  253.     OTHER=(INP .EQ. 'O')
  254.     REALDAT=(INP .EQ. 'R')
  255.     IF (.NOT. VOLT .AND. .NOT. OTHER .AND.
  256.      +      .NOT. REALDAT) GO TO 50
  257.     IF (REALDAT) GO TO 70
  258.  
  259. *    Get the voltage transformation factor.
  260.  
  261. 60    WRITE (*,'(/1X,A\)')
  262.      +   'Alternate voltage transformation factor (Y/N) : '
  263.     READ (*,'(A)') INP
  264.     IF (INP .EQ. 'y') INP = 'Y'
  265.     IF (INP .EQ. 'n') INP = 'N'
  266.     ALTVOLT=(INP .EQ. 'Y')
  267.     USEGAIN=(INP .EQ. 'N')
  268.     IF (.NOT. ALTVOLT .AND. .NOT. USEGAIN) GO TO 60
  269.  
  270.     IF (ALTVOLT) THEN
  271.       WRITE (*,'(/1X,A\)') 'Enter voltage transformation factor : '
  272.       READ (*,*) VOFST1
  273.     ENDIF
  274.  
  275. *    Remove the mean value from the file ?
  276.  
  277. 70    WRITE (*,'(/1X,A,A\)') 'Remove (M)ean value,',
  278.      +         ' (D)etrend or (N)either : '
  279.     READ (*,'(A)') INP
  280.     IF (INP .EQ. 'm') INP = 'M'
  281.     IF (INP .EQ. 'd') INP = 'D'
  282.     IF (INP .EQ. 'n') INP = 'N'
  283.     RMEAN=(INP .EQ. 'M')
  284.     DTREND=(INP .EQ. 'D')
  285.     NOOP=(INP .EQ. 'N')
  286.     IF (.NOT. RMEAN .AND. .NOT. DTREND .AND. .NOT. NOOP) GO TO 70
  287.     WRITE (*,'( )')
  288.     IF (VOLT .OR. REALDAT) GO TO 90
  289.  
  290. *    Determine the type of data file for coefficients.
  291.  
  292. 80    WRITE (*,'(1X,A,A\)') 'Read from (A)scii or',
  293.      +                        ' (B)inary data file : '
  294.     READ (*,'(A)') INP
  295.     IF (INP .EQ. 'a') INP = 'A'
  296.     IF (INP .EQ. 'b') INP = 'B'
  297.     ASCII=(INP .EQ. 'A')
  298.     BINRY=(INP .EQ. 'B')
  299.     IF (.NOT. ASCII .AND. .NOT. BINRY) GO TO 80
  300.     WRITE (*,'( )')
  301.  
  302. *    Allocate space for CONST array.
  303.  
  304.     IF (CONSEC) ALLOCATE (CONST(IFSTRT:IFSTP,NUMCON), STAT=IERR)
  305.     IF (ARB) ALLOCATE (CONST(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
  306.     IF (IERR .NE. 0)
  307.      +     STOP 'Problem allocating storage for CONST.  Aborting ...'
  308.  
  309. *    Open the formatted file containing the transformation
  310. *    constants and fill the array.
  311.  
  312.     IF (ASCII) THEN
  313.  
  314.       INQUIRE (FILE=YNAM,EXIST=EXISTS)
  315.       IF (.NOT. EXISTS) THEN
  316.         IPD=INDEX (YNAM,'.',.TRUE.)
  317.         WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
  318.      +                              YNAM(:IPD+3)
  319.         STOP '                   Program terminated unsuccessfully.'
  320.       ENDIF
  321.  
  322.       OPEN (NUMI,FILE=YNAM,STATUS='OLD')
  323.         READ (NUMI,*) (B(I),I=1,5)
  324.       CLOSE (NUMI,ERR=400)
  325.  
  326.       DO I=IFSTRT,IFSTP
  327.         IF (CONSEC) IF=I
  328.         IF (ARB) IF=IFIL(I)
  329.         DO J=1,5
  330.           CONST(IF,J)=B(J)
  331.         ENDDO
  332.       ENDDO
  333.  
  334.     ELSE IF (BINRY) THEN
  335.  
  336. *      Open the binary file containing the transformation
  337. *      constants and fill the array.
  338.  
  339.       INQUIRE (FILE=XNAM,EXIST=EXISTS)
  340.       IF (.NOT. EXISTS) THEN
  341.         IPD=INDEX (XNAM,'.',.TRUE.)
  342.         WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
  343.      +                              XNAM(:IPD+3)
  344.         STOP '                   Program terminated unsuccessfully.'
  345.       ENDIF
  346.  
  347.       OPEN (NUMI,FILE=XNAM,STATUS='OLD',ACCESS='SEQUENTIAL',
  348.      +          FORM='BINARY')
  349.         READ (NUMI) NUMSETS,NSTART
  350.         IF (
  351.      +      (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
  352.      +      .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
  353.      +                         NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
  354.           WRITE (*,'(/1X,A/1X,A/1X,A/)')
  355.      +          'Starting or ending data file #''s stored in coefficient',
  356.      +        'data file for channel 1 do not match values entered',
  357.      +        'into program.  This must be corrected before proceeding.'
  358.           CLOSE (NUMI)
  359.           STOP '                   Program terminated unsuccessfully.'
  360.         ENDIF          
  361.         READ (NUMI)
  362.      +         ((CONST(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
  363.       CLOSE (NUMI,ERR=400)
  364.  
  365.     ENDIF
  366.  
  367. *    Filter the data ?
  368.  
  369. 90    WRITE (*,'(1X,A/20X,A\)')
  370.      +    'Filter the data using (L)owpass, (H)ighpass,',
  371.      +                      '(B)andpass, Band(S)top or (N)one : '
  372.     READ (*,'(A)') FTYPE
  373.     IF (FTYPE .EQ. 'l') FTYPE = 'L'
  374.     IF (FTYPE .EQ. 'h') FTYPE = 'H'
  375.     IF (FTYPE .EQ. 'b') FTYPE = 'B'
  376.     IF (FTYPE .EQ. 's') FTYPE = 'S'
  377.     IF (FTYPE .EQ. 'n') FTYPE = 'N'
  378.     LPASS=(FTYPE .EQ. 'L')
  379.     HPASS=(FTYPE .EQ. 'H')
  380.     BPASS=(FTYPE .EQ. 'B')
  381.     BSTOP=(FTYPE .EQ. 'S')
  382.     NOFIL=(FTYPE .EQ. 'N')
  383.     IF (.NOT. LPASS .AND. .NOT. HPASS .AND. .NOT. BPASS .AND.
  384.      +      .NOT. BSTOP .AND. .NOT. NOFIL) GO TO 90
  385.     IF (NOFIL) GO TO 110
  386.  
  387. *    Get the cutoff frequencies.
  388.  
  389.     FC=0.0
  390.     F1=0.0
  391.     F2=0.0
  392.  
  393.     IF (LPASS .OR. HPASS) THEN
  394.       WRITE (*,'(/1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
  395.       READ (*,*) FC
  396.     ELSE IF (BPASS .OR. BSTOP) THEN
  397.       WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
  398.       READ (*,*) F1
  399.       WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
  400.       READ (*,*) F2
  401.     ENDIF
  402.  
  403. 100    WRITE (*,'(/1X,A/5X,A/1X,A\)')
  404.      +             'The order of the filter will be twice',
  405.      +             'the number of cascaded sections.',
  406.      +             'Enter # of filter sections to cascade : '
  407.     READ (*,*) NS
  408.     IF (NS .GT. NSMAX) THEN
  409.       WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
  410.       GO TO 100
  411.     ENDIF
  412.  
  413. *    Two channels of data ?
  414.  
  415. 110    WRITE (*,'(/1X,A\)') 'How many channels of data (1 or 2) : '
  416.     READ (*,*) ICHANS
  417.     WRITE (*,'( )')
  418.     ONECHAN=(ICHANS .EQ. 1)
  419.     TWOCHAN=(ICHANS .EQ. 2)
  420.     IF (.NOT. ONECHAN .AND. .NOT. TWOCHAN) GO TO 110
  421.     IF (ONECHAN) GO TO 170
  422.     IF (TWOCHAN .AND. REALDAT) GO TO 130
  423.  
  424. *    Transform second channel to voltages or to some other quantity ?
  425.  
  426. 120    WRITE (*,'(1X,A,A\)') 'Transform 2nd channel into',
  427.      +                ' (V)oltages or into (O)ther quantity : '
  428.     READ (*,'(A)') INP
  429.     IF (INP .EQ. 'v') INP = 'V'
  430.     IF (INP .EQ. 'o') INP = 'O'
  431.     VOLT2=(INP .EQ. 'V')
  432.     OTHER2=(INP .EQ. 'O')
  433.     IF (.NOT. VOLT2 .AND. .NOT. OTHER2) GO TO 120
  434.     WRITE (*,'( )')
  435.  
  436. *    Get the voltage transformation factor for channel 2.
  437.  
  438. 125    WRITE (*,'(1X,A\)')
  439.      +   'Alternate voltage transformation factor (Y/N) : '
  440.     READ (*,'(A)') INP
  441.     IF (INP .EQ. 'y') INP = 'Y'
  442.     IF (INP .EQ. 'n') INP = 'N'
  443.     ALTVOLT2=(INP .EQ. 'Y')
  444.     USEGAIN2=(INP .EQ. 'N')
  445.     IF (.NOT. ALTVOLT2 .AND. .NOT. USEGAIN2) GO TO 125
  446.     WRITE (*,'( )')
  447.  
  448.     IF (ALTVOLT2) THEN
  449.       WRITE (*,'(1X,A\)') 'Enter voltage transformation factor : '
  450.       READ (*,*) VOFST2
  451.       WRITE (*,'( )')
  452.     ENDIF
  453.  
  454. *    Remove the mean value from the second channel of data ?
  455.  
  456. 130    WRITE (*,'(1X,A,A\)') '2nd channel: remove (M)ean value,',
  457.      +         ' (D)etrend or (N)either : '
  458.     READ (*,'(A)') INP
  459.     IF (INP .EQ. 'm') INP = 'M'
  460.     IF (INP .EQ. 'd') INP = 'D'
  461.     IF (INP .EQ. 'n') INP = 'N'
  462.     RMEAN2=(INP .EQ. 'M')
  463.     DTREND2=(INP .EQ. 'D')
  464.     NOOP2=(INP .EQ. 'N')
  465.     IF (.NOT. RMEAN2 .AND. .NOT. DTREND2 .AND.
  466.      +      .NOT. NOOP2) GO TO 130
  467.     WRITE (*,'( )')
  468.     IF (VOLT2 .OR. REALDAT) GO TO 150
  469.  
  470. *    Determine type of data file for coefficients for second channel.
  471.  
  472. 140    WRITE (*,'(1X,A,A\)') '2nd channel: read from (A)scii or',
  473.      +                        ' (B)inary data file : '
  474.     READ (*,'(A)') INP
  475.     IF (INP .EQ. 'a') INP = 'A'
  476.     IF (INP .EQ. 'b') INP = 'B'
  477.     ASCII2=(INP .EQ. 'A')
  478.     BINRY2=(INP .EQ. 'B')
  479.     IF (.NOT. ASCII2 .AND. .NOT. BINRY2) GO TO 140
  480.     WRITE (*,'( )')
  481.  
  482. *    Allocate space for CONST2 array.
  483.  
  484.     IF (CONSEC) ALLOCATE (CONST2(IFSTRT:IFSTP,NUMCON), STAT=IERR)
  485.     IF (ARB) ALLOCATE (CONST2(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
  486.     IF (IERR .NE. 0)
  487.      +     STOP 'Problem allocating storage for CONST2.  Aborting ...'
  488.  
  489. *    Open the formatted file containing the transformation
  490. *    constants for the second channel and fill the array.
  491.  
  492.     IF (ASCII2) THEN
  493.  
  494.       INQUIRE (FILE=YNAM2,EXIST=EXISTS)
  495.       IF (.NOT. EXISTS) THEN
  496.         IPD=INDEX (YNAM2,'.',.TRUE.)
  497.         WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
  498.      +                              YNAM2(:IPD+3)
  499.         STOP '                   Program terminated unsuccessfully.'
  500.       ENDIF
  501.  
  502.       OPEN (NUMI,FILE=YNAM2,STATUS='OLD')
  503.         READ (NUMI,*) (B2(I),I=1,5)
  504.       CLOSE (NUMI,ERR=400)
  505.  
  506.       DO I=IFSTRT,IFSTP
  507.         IF (CONSEC) IF=I
  508.         IF (ARB) IF=IFIL(I)
  509.         DO J=1,5
  510.           CONST2(IF,J)=B2(J)
  511.         ENDDO
  512.       ENDDO
  513.  
  514.     ELSE IF (BINRY2) THEN
  515.  
  516. *      Open the binary file containing the transformation
  517. *      constants for the second channel and fill the array.
  518.  
  519.       INQUIRE (FILE=XNAM2,EXIST=EXISTS)
  520.       IF (.NOT. EXISTS) THEN
  521.         IPD=INDEX (XNAM2,'.',.TRUE.)
  522.         WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
  523.      +                              XNAM2(:IPD+3)
  524.         STOP '                   Program terminated unsuccessfully.'
  525.       ENDIF
  526.  
  527.       OPEN (NUMI,FILE=XNAM2,STATUS='OLD',ACCESS='SEQUENTIAL',
  528.      +          FORM='BINARY')
  529.         READ (NUMI) NUMSETS,NSTART
  530.         IF (
  531.      +      (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
  532.      +      .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
  533.      +                         NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
  534.           WRITE (*,'(/1X,A/1X,A/1X,A/)')
  535.      +          'Starting or ending data file #''s stored in coefficient',
  536.      +        'data file for channel 2 do not match values entered',
  537.      +        'into program.  This must be corrected before proceeding.'
  538.           CLOSE (NUMI)
  539.           STOP '                   Program terminated unsuccessfully.'
  540.         ENDIF          
  541.         READ (NUMI)
  542.      +        ((CONST2(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
  543.       CLOSE (NUMI,ERR=400)
  544.  
  545.     ENDIF
  546.  
  547. *    Filter the data for the second channel ?
  548.  
  549. 150    WRITE (*,'(1X,A/20X,A\)')
  550.      +    'Filter the data using (L)owpass, (H)ighpass,',
  551.      +                      '(B)andpass, Band(S)top or (N)one : '
  552.     READ (*,'(A)') FTYP2
  553.     IF (FTYP2 .EQ. 'l') FTYP2 = 'L'
  554.     IF (FTYP2 .EQ. 'h') FTYP2 = 'H'
  555.     IF (FTYP2 .EQ. 'b') FTYP2 = 'B'
  556.     IF (FTYP2 .EQ. 's') FTYP2 = 'S'
  557.     IF (FTYP2 .EQ. 'n') FTYP2 = 'N'
  558.     LPASS2=(FTYP2 .EQ. 'L')
  559.     HPASS2=(FTYP2 .EQ. 'H')
  560.     BPASS2=(FTYP2 .EQ. 'B')
  561.     BSTOP2=(FTYP2 .EQ. 'S')
  562.     NOFIL2=(FTYP2 .EQ. 'N')
  563.     IF (.NOT. LPASS2 .AND. .NOT. HPASS2 .AND. .NOT. BPASS2 .AND.
  564.      +      .NOT. BSTOP2 .AND. .NOT. NOFIL2) GO TO 150
  565.     WRITE (*,'( )')
  566.     IF (NOFIL2) GO TO 170
  567.  
  568. *    Get the cutoff frequencies.
  569.  
  570.     FC2=0.0
  571.     F12=0.0
  572.     F22=0.0
  573.  
  574.     IF (LPASS2 .OR. HPASS2) THEN
  575.       WRITE (*,'(1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
  576.       READ (*,*) FC2
  577.     ELSE IF (BPASS2 .OR. BSTOP2) THEN
  578.       WRITE (*,'(1X,A\)') 'Enter the low frequency cutoff : '
  579.       READ (*,*) F12
  580.       WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
  581.       READ (*,*) F22
  582.     ENDIF
  583.  
  584. 160    WRITE (*,'(/1X,A\)') 'Enter # of filter sections to cascade : '
  585.     READ (*,*) NS2
  586.     WRITE (*,'( )')
  587.     IF (NS2 .GT. NSMAX) THEN
  588.       WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
  589.       GO TO 160
  590.     ENDIF
  591.  
  592. *    Determine the kind of spectrum.
  593.  
  594. 170    WRITE (*,'(1X,A\)')
  595.      +            '(A)mplitude, (P)ower or (C)ross spectrum : '
  596.     READ (*,'(A)') INP
  597.     WRITE (*,'( )')
  598.     IF (INP .EQ. 'a') INP='A'
  599.     IF (INP .EQ. 'p') INP='P'
  600.     IF (INP .EQ. 'c') INP='C'
  601.     AMP=(INP .EQ. 'A')
  602.     POWER=(INP .EQ. 'P')
  603.     CROSS=(INP .EQ. 'C')
  604.     IF (.NOT. AMP .AND. .NOT. POWER .AND. .NOT. CROSS) GO TO 170
  605.     IF (CROSS .AND. ONECHAN) THEN
  606.       WRITE (*,'(1X,A,A/)') 'Cross Spectrum requires',
  607.      +                          ' two channels of data.'
  608.       GO TO 170
  609.     ENDIF
  610.  
  611. *    Calculate the coherence function ?
  612.  
  613. 180    IF (CROSS) THEN
  614.       WRITE (*,'(1X,A/1X,A//1X,A\)')
  615.      +    'Coherence function requires calculation',
  616.      +    'of both cross and autospectral density functions.',
  617.      +    'Calculate the coherence function gamma (f) (Y/N) : '
  618.       READ (*,'(A)') INP
  619.       IF (INP .EQ. 'y') INP = 'Y'
  620.       IF (INP .EQ. 'n') INP = 'N'
  621.       COH=(INP .EQ. 'Y')
  622.       NOCOH=(INP .EQ. 'N')
  623.       IF (.NOT. COH .AND. .NOT. NOCOH) GO TO 180
  624.       WRITE (*,'( )')
  625.     ENDIF
  626.  
  627. *    Taper the data or not ?
  628.  
  629. 190    WRITE (*,'(1X,A\)') 'Taper the data (Y/N) : '
  630.     READ (*,'(A)') INP
  631.     WRITE (*,'( )')
  632.     IF (INP .EQ. 'y') INP='Y'
  633.     IF (INP .EQ. 'n') INP='N'
  634.     TAPER=(INP .EQ. 'Y')
  635.     NOTAPER=(INP .EQ. 'N')
  636.     IF (.NOT. TAPER .AND. .NOT. NOTAPER) GO TO 190
  637.  
  638. *    Choose the window function.
  639.  
  640.     IF (TAPER) THEN
  641. 200      WRITE (*,'(1X,A/1X,A\)')
  642.      +           '(H)anning, (P)arzen, (T)ukey, (W)elch,',
  643.      +         'Min (S)ide Lobe or Max (D)ecay window : '
  644.       READ (*,'(A)') INP
  645.       WRITE (*,'( )')
  646.       IF (INP .EQ. 'h') INP='H'
  647.       IF (INP .EQ. 'p') INP='P'
  648.       IF (INP .EQ. 't') INP='T'
  649.       IF (INP .EQ. 'w') INP='W'
  650.       IF (INP .EQ. 's') INP='S'
  651.       IF (INP .EQ. 'd') INP='D'
  652.       HANN =(INP .EQ. 'H')
  653.       PARZ =(INP .EQ. 'P')
  654.       TUKEY=(INP .EQ. 'T')
  655.       WELCH=(INP .EQ. 'W')
  656.       MINSL=(INP .EQ. 'S')
  657.       MXDEC=(INP .EQ. 'D')
  658.       IF (.NOT. HANN  .AND. .NOT. PARZ  .AND.
  659.      +        .NOT. WELCH .AND. .NOT. TUKEY .AND.
  660.      +        .NOT. MINSL .AND. .NOT. MXDEC) GO TO 200
  661.     ENDIF
  662.  
  663. *    Overlap the data or not ?
  664.  
  665.     IF (TAPER) THEN
  666. 205      WRITE (*,'(1X,A\)') 'Use overlap (Y/N) : '
  667.       READ (*,'(A)') INP
  668.       WRITE (*,'( )')
  669.       IF (INP .EQ. 'y') INP='Y'
  670.       IF (INP .EQ. 'n') INP='N'
  671.       OVRLP  =(INP .EQ. 'Y')
  672.       NOOVRLP=(INP .EQ. 'N')
  673.       IF (.NOT. OVRLP .AND. .NOT. NOOVRLP) GO TO 205
  674.     ENDIF
  675.  
  676. *    Get the time just before starting.
  677.  
  678.     CALL GETTIM (IHR,IMIN,ISEC,I100TH)
  679.  
  680. *    Do the calculations for each of the desired files.
  681.  
  682.     DO KOUNT=IFSTRT,IFSTP
  683.       IF (CONSEC) IF=KOUNT
  684.       IF (ARB) IF=IFIL(KOUNT)
  685.  
  686. *      Assign filenames.
  687.  
  688.       IF (NUMARGS .GT. 0) THEN
  689.         FIN=FLNM(IF)
  690.       ELSE
  691.         IF (.NOT. ROOT1) THEN
  692.           ISL=INDEX (PTH,'\',.TRUE.)
  693.           WRITE (FIN,'(A,A1,I3.3,A4)') PTH(:ISL),FIRST,IF,'.DAT'
  694.         ELSE
  695.           WRITE (FIN,'(A1,I3.3,A4)') FIRST,IF,'.DAT'
  696.         ENDIF
  697.         INQUIRE (FILE=FIN,EXIST=EXISTS)
  698.         IF (.NOT. EXISTS) THEN
  699.           IPD=INDEX (FIN,'.',.TRUE.)
  700.           WRITE (*,'(25X,A,A,A/)') 'Cannot find file : ',
  701.      +                                  FIN(:IPD+3),'.'
  702.           CYCLE
  703.         ENDIF
  704.       ENDIF
  705.  
  706.       IPD=INDEX (FIN,'.',.TRUE.)
  707.       IF (ROOT2) THEN
  708.         TDOUT =FIN(IPD-4:IPD-1) // '.TRF'
  709.         TDOUT2=FIN(IPD-4:IPD-1) // '.TR2'
  710.         TOUT  =FIN(IPD-4:IPD-1) // '.TAP'
  711.         TOUT2 =FIN(IPD-4:IPD-1) // '.TP2'
  712.       ELSE IF (.NOT. ROOT2) THEN
  713.         ISL=INDEX (TPTH,'\',.TRUE.)
  714.         TDOUT =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TRF'
  715.         TDOUT2=TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TR2'
  716.         TOUT  =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TAP'
  717.         TOUT2 =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TP2'
  718.       ENDIF
  719.  
  720.       IF (.NOT. SUFX) THEN
  721.         IF (POWER) POUT=FIN(:IPD-1) // 'P.PRN'
  722.         IF (AMP)   POUT=FIN(:IPD-1) // 'A.PRN'
  723.         IF (CROSS) POUT=FIN(:IPD-1) // 'C.PRN'
  724.         EOUT  =FIN(:IPD-1) // 'E.PRN'
  725.       ELSE
  726.         IF (POWER) POUT=FIN(:IPD-1) // 'P' // SUFFIX
  727.         IF (AMP)   POUT=FIN(:IPD-1) // 'A' // SUFFIX
  728.         IF (CROSS) POUT=FIN(:IPD-1) // 'C' // SUFFIX
  729.         EOUT  =FIN(:IPD-1) // 'E' // SUFFIX
  730.       ENDIF
  731.  
  732. *      Put message on screen.
  733.  
  734.       IPD=INDEX (FIN,'.',.TRUE.)
  735.       IF (KOUNT .EQ. IFSTRT) THEN
  736.         WRITE (*,'(/////////////////////////13X,
  737.      +    ''S P E C T R U M   C A L C U L A T I O N   R O U T I N E'')')
  738.         WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
  739.       ELSE
  740.         WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
  741.       ENDIF
  742.  
  743. *      Open the data file and input the data.
  744.  
  745.       OPEN (NUMI,FILE=FIN,STATUS='OLD',ACCESS='SEQUENTIAL',
  746.      +          FORM='BINARY')
  747.  
  748.       READ (NUMI,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
  749.       READ (NUMI,ERR=400) (GAIN(I),I=0,7)
  750.  
  751. *      Calculate # points per record and adjust record size for
  752. *      the file containing transformed data.
  753.  
  754.       IF (REALDAT) THEN
  755.         N=IRSIZE/4
  756.         IRSIZ2=IRSIZE
  757.         IRSIZE=IRSIZE/2
  758.       ELSE
  759.         N=IRSIZE/2
  760.         IRSIZ2=2*IRSIZE
  761.       ENDIF
  762.       DELT=FLOAT(IDELTMS)/1000000.0
  763.       FN=FLOAT(N)
  764.  
  765. *      N less than or equal to NMAX error checking.
  766.  
  767.       IF (ONECHAN) NTST=NMAX
  768.       IF (TWOCHAN) NTST=NMAX/2
  769.       IF (N .GT. NMAX) THEN
  770.         WRITE (*,'(15X,A,I5,A/15X,A,I5/)')
  771.      +         'This data file contains ',N,' data points.',
  772.      +         '# data points per record per channel must be < ',NTST
  773.         STOP '                   Program terminated unsuccessfully.'
  774.       ENDIF
  775.  
  776. *      Test to make sure that N is a power of 2.
  777.  
  778.       ITST=NINT(ALOG10(FN)/ALOG10(2.0))
  779.       ITST2=INT(2**ITST)-N
  780.  
  781.       IF (ITST2 .NE. 0) THEN
  782.         WRITE (*,'(20X,A,I5,A/20X,A/)') 'This data file contains ',
  783.      +          N,' data points.','# data points must be a power of 2.'
  784.         STOP '                   Program terminated unsuccessfully.'
  785.       ENDIF
  786.  
  787. *      Allocate space for various arrays.
  788.  
  789.       IF (REALDAT) THEN
  790.         ALLOCATE (WIN(N), FR (N), FI (N),
  791.      +                                  FR2(N), FI2(N), STAT=IERR)
  792.         IF (IERR .NE. 0)
  793.      +              STOP 'Not enough storage for data.  Aborting ...'
  794.       ELSE
  795.         ALLOCATE (NDATA(N), WIN(N), FR (N), FI (N),
  796.      +                                  FR2(N), FI2(N), STAT=IERR)
  797.         IF (IERR .NE. 0)
  798.      +              STOP 'Not enough storage for data.  Aborting ...'
  799.       ENDIF
  800.  
  801. *      Test to be sure data file really has one channel.
  802.  
  803.       IF (ONECHAN .AND. ICHANS .NE. 1) THEN
  804.         WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
  805.      +      ' file indicates two channels.',
  806.      +      'You have specified one channel.  Please correct problem.'
  807.         STOP '                       Program terminated unsuccessfully
  808.      +.'
  809.       ENDIF
  810.  
  811. *      Test to be sure data file really has two channels.
  812.  
  813.       IF (TWOCHAN .AND. ICHANS .NE. 2) THEN
  814.         WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
  815.      +      ' file indicates only one channel.',
  816.      +      'You have specified two channels.  Please correct problem.'
  817.         STOP '                       Program terminated unsuccessfully
  818.      +.'
  819.       ENDIF
  820.  
  821. *      Estimate amount of temporary space needed.
  822.  
  823.       IF (REALDAT) THEN
  824.         NBYTE=IRSIZ2*NUMREC
  825.       ELSE
  826.         NBYTE=IRSIZE*NUMREC
  827.       ENDIF
  828.  
  829.       IF (.NOT. REALDAT) TRFFIL=2*NBYTE
  830.       IF (REALDAT) TRFFIL=NBYTE
  831.       IF (NOTAPER) TAPFIL=0
  832.       IF (TAPER .AND. OVRLP) TAPFIL=2*TRFFIL
  833.       IF (TAPER .AND. NOOVRLP) TAPFIL=TRFFIL
  834.  
  835.       TEMP_SPACE=TRFFIL+TAPFIL
  836.       TSPCMB=FLOAT(TEMP_SPACE)/1048576.0
  837.  
  838.       WRITE (*,'(13X,A,I8,A,F6.2,A/)')
  839.      +      'Temporary space required = ',TEMP_SPACE,
  840.      +      ' bytes (',TSPCMB,' Mb).'
  841.  
  842. *      Perform a write test to determine if space is available
  843. *      to hold the temporary data files.
  844.  
  845.       OPEN (NUMT, FILE=TDOUT, FORM='BINARY',RECL=1,
  846.      +        ACCESS='DIRECT', STATUS='UNKNOWN')
  847.         IREC=TEMP_SPACE+(HEADER_SIZE*2)-1
  848.         WRITE (NUMT, REC=IREC, ERR=390) #FF
  849.       CLOSE (NUMT, STATUS='DELETE')
  850.  
  851. *      Determine the voltage transformation factor.
  852.  
  853.       IF (USEGAIN)  VOFST1=(20.0/4096.0)/(2**GAIN(0))
  854.       IF (USEGAIN2) VOFST2=(20.0/4096.0)/(2**GAIN(1))
  855.  
  856. *      Open the output file to hold the transformed data.
  857.  
  858.       OPEN (NUMO,FILE=TDOUT,STATUS='UNKNOWN',
  859.      +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
  860.       IF (TWOCHAN) OPEN (NUMO2,FILE=TDOUT2,STATUS='UNKNOWN',
  861.      +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
  862.  
  863. *      Write file header.
  864.  
  865.       IF (ONECHAN) THEN
  866.         WRITE (NUMO,ERR=400) ICHANS,IRSIZ2,NUMREC,IDELTMS
  867.       ELSE IF (TWOCHAN) THEN
  868.         WRITE (NUMO,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
  869.         WRITE (NUMO2,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
  870.       ENDIF
  871.  
  872. *      Transformation procedure for one channel of data.
  873.  
  874.       IF (ONECHAN) THEN
  875.  
  876.         IF (OTHER) THEN
  877.           DO J=1,NUMCON
  878.             B(J)=CONST(IF,J)
  879.           ENDDO
  880.         ELSE
  881.           B=0.0
  882.         ENDIF
  883.  
  884.         IF (TRANS_ONE_CHAN (NDATA,FR,FI,B) .NE. 0) STOP
  885.      +        '                      Program terminated unsuccessfully.'
  886.  
  887. *        Close the input data file and proceed
  888. *        to the tapering operation.
  889.  
  890.         CLOSE (NUMI,STATUS='KEEP',ERR=400)
  891.  
  892. *      Transformation procedure for two channels of data.
  893.  
  894.       ELSE IF (TWOCHAN) THEN
  895.  
  896.         IF (OTHER) THEN
  897.           DO J=1,NUMCON
  898.             B(J)=CONST(IF,J)
  899.           ENDDO
  900.         ELSE
  901.           B=0.0
  902.         ENDIF
  903.  
  904.         IF (OTHER2) THEN
  905.           DO J=1,NUMCON
  906.             B2(J)=CONST2(IF,J)
  907.           ENDDO
  908.         ELSE
  909.           B2=0.0
  910.         ENDIF
  911.  
  912.         IF (TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)
  913.      +        .NE. 0) STOP
  914.      +        '                      Program terminated unsuccessfully.'
  915.  
  916. *        Close the input data file and proceed
  917. *        to the tapering operation.
  918.  
  919.         CLOSE (NUMI,STATUS='KEEP',ERR=400)
  920.  
  921.       ENDIF
  922.  
  923. *      Messages concerning mean removal, detrending or filtering.
  924.  
  925.       MSG=.FALSE.
  926.       IF (TWOCHAN) THEN
  927.         IF (RMEAN) WRITE (*,'(25X,A,A)') 'Mean value removed',
  928.      +                       ' from channel 1.'
  929.         IF (DTREND) WRITE (*,'(25X,A)') 'Channel 1 detrended.'
  930.         IF (LPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
  931.      +                                    ' was lowpass filtered.'
  932.         IF (HPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
  933.      +                                    ' was highpass filtered.'
  934.         IF (BPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
  935.      +                                    ' was bandpass filtered.'
  936.         IF (BSTOP) WRITE (*,'(25X,A,A)') 'Channel 1',
  937.      +                                    ' was bandstop filtered.'
  938.         IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
  939.      +                                'Order of filter = ',2*NS,'.'
  940.         IF (RMEAN2) WRITE (*,'(25X,A,A)') 'Mean value removed',
  941.      +                       ' from channel 2.'
  942.         IF (DTREND2) WRITE (*,'(25X,A)') 'Channel 2 detrended.'
  943.         IF (LPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
  944.      +                                    ' was lowpass filtered.'
  945.         IF (HPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
  946.      +                                    ' was highpass filtered.'
  947.         IF (BPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
  948.      +                                    ' was bandpass filtered.'
  949.         IF (BSTOP2) WRITE (*,'(25X,A,A)') 'Channel 2',
  950.      +                                    ' was bandstop filtered.'
  951.         IF (.NOT. NOFIL2) WRITE (*,'(25X,A,I2,A)')
  952.      +                               'Order of filter = ',2*NS2,'.'
  953.         IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
  954.         IF (RMEAN2 .OR. DTREND2 .OR. .NOT. NOFIL) MSG=.TRUE.
  955.       ELSE
  956.         IF (RMEAN) WRITE (*,'(25X,A)') 'Mean value removed.'
  957.         IF (DTREND) WRITE (*,'(25X,A)') 'File was detrended.'
  958.         IF (LPASS) WRITE (*,'(25X,A)') 'File was lowpass filtered.'
  959.         IF (HPASS) WRITE (*,'(25X,A)') 'File was highpass filtered.'
  960.         IF (BPASS) WRITE (*,'(25X,A)') 'File was bandpass filtered.'
  961.         IF (BSTOP) WRITE (*,'(25X,A)') 'File was bandstop filtered.'
  962.         IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
  963.      +                                     'Order of filter = ',2*NS,'.'
  964.         IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
  965.       ENDIF
  966.  
  967. *      Read the file header.
  968.  
  969.       REWIND (NUMO)
  970.       READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
  971.       IF (TWOCHAN) THEN
  972.         REWIND (NUMO2)
  973.         READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
  974.       ENDIF
  975.  
  976.       N=IRSIZE/4
  977.       N2=N/2
  978.       NREC=NUMREC             ! original numrec before overlap
  979.       NREC2=2*NUMREC
  980.       K=INT4(FLOAT(N)/10.0)
  981.  
  982. *      No Tapering applied.
  983.  
  984.       IF (NOTAPER) THEN
  985.         IF (MSG) THEN
  986.           WRITE (*,'(/25X,A/25X,I4.4,A)') 'No tapering applied.  ',
  987.      +               NUMREC, ' records will be used.'
  988.         ELSE
  989.           WRITE (*,'(25X,A/25X,I4.4,A)') 'No tapering applied.  ',
  990.      +               NUMREC, ' records will be used.'
  991.         ENDIF
  992.         GO TO 300
  993.       ENDIF
  994.  
  995. *      Compute the scale factor.
  996.  
  997.       IF (WIN_SCALE (WIN) .NE. 0) STOP
  998.      +      '                      Program terminated unsuccessfully.'
  999.  
  1000. *      Open the data files to hold tapered data.
  1001.  
  1002.       OPEN (NUMT,FILE=TOUT,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
  1003.      +          FORM='BINARY',ERR=400)
  1004.  
  1005.       IF (TWOCHAN) OPEN (NUMT2,FILE=TOUT2,STATUS='UNKNOWN',
  1006.      +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
  1007.  
  1008.       IF (NOOVRLP) THEN
  1009.  
  1010. *        No overlap - use NUMREC records.
  1011.  
  1012.         WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
  1013.         IF (TWOCHAN) WRITE (NUMT2,ERR=400)
  1014.      +                   ICHANS,IRSIZE,NUMREC,IDELTMS
  1015.       ELSE
  1016.  
  1017. *        50% overlap - use NREC2 records.
  1018.  
  1019.         WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NREC2,IDELTMS
  1020.         IF (TWOCHAN) WRITE (NUMT2,ERR=400)
  1021.      +                   ICHANS,IRSIZE,NREC2,IDELTMS
  1022.       ENDIF
  1023.  
  1024.       IF (MSG) WRITE (*,'( )')
  1025.       IF (HANN)
  1026.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // HNAM // '.'
  1027.       IF (WELCH)
  1028.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // WNAM // '.'
  1029.       IF (PARZ)
  1030.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // PNAM // '.'
  1031.       IF (TUKEY)
  1032.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // TNAM // '.'
  1033.       IF (MINSL)
  1034.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // SNAM // '.'
  1035.       IF (MXDEC)
  1036.      +      WRITE (*,'(25X,A)') 'Taper applied = ' // DNAM // '.'
  1037.  
  1038.       IF (OVRLP) THEN
  1039.         WRITE (*,'(25X,A/25X,I4.4,A/)') '50% overlap used.',
  1040.      +                  NREC2,' records will be written.'
  1041.       ELSE
  1042.         WRITE (*,'(25X,A/25X,I4.4,A/)') 'No overlap used.',
  1043.      +                  NUMREC,' records will be written.'
  1044.       ENDIF
  1045.  
  1046. *      Do the tapering now.
  1047.  
  1048.       IF (TAPER_DATA (FR,FR2,FI,FI2,WIN) .NE. 0) STOP
  1049.      +      '                      Program terminated unsuccessfully.'
  1050.  
  1051. *      Throw away transformed data file and prepare tapered data file.
  1052.  
  1053.       CLOSE (NUMO,STATUS='DELETE',ERR=400)
  1054.       IF (TWOCHAN) CLOSE (NUMO2,STATUS='DELETE',ERR=400)
  1055.       REWIND (NUMT)
  1056.       READ (NUMT) ICHANS,IRSIZE,NUMREC,IDELTMS
  1057.       IF (TWOCHAN) THEN
  1058.         REWIND (NUMT2)
  1059.         READ (NUMT2) ICHANS,IRSIZE,NUMREC,IDELTMS
  1060.       ENDIF
  1061.  
  1062. *      Begin spectrum computations.
  1063.  
  1064. 300      IF (NOTAPER) THEN
  1065.         INUM=NUMO
  1066.         INUM2=NUMO2
  1067.       ELSE IF (TAPER) THEN
  1068.         INUM=NUMT
  1069.         INUM2=NUMT2
  1070.       ENDIF
  1071.  
  1072. *      Test for odd number of records.  This should only occur
  1073. *      at this point for the one channel, no taper or one channel,
  1074. *      taper and no overlap cases.
  1075.  
  1076.       SIZTST=FLOAT(NUMREC)      
  1077.       REMAIN=AMOD(SIZTST,2.0)
  1078.       EVENUM=(REMAIN .EQ. 0.0)
  1079.       ODDNUM=(REMAIN .NE. 0.0)
  1080.  
  1081. *      Initialization.
  1082.  
  1083.       N=IRSIZE/4
  1084.       N2=N/2
  1085.       NP1=N+1
  1086.       IOD=INT4(ALOG10(FLOAT(N))/ALOG10(2.0)+.5)
  1087.       NRECD2=NUMREC/2
  1088.       IF (ODDNUM) NRECD2=(NUMREC+1)/2
  1089.       IF (OVRLP) NUMREC=NUMREC-1
  1090.       XMEAN=0.0
  1091.       YMEAN=0.0
  1092.       DO I=1,N
  1093.         FR2(I)=0.0
  1094.         FI2(I)=0.0
  1095.         WIN(I)=0.0
  1096.       ENDDO
  1097.  
  1098.       IF (ONECHAN) THEN
  1099.  
  1100. *        Spectrum computations for one channel.
  1101.  
  1102.         IF (SPEC_ONE_CHAN (FR,FI,FR2,FI2) .NE. 0) STOP
  1103.      +        '                      Program terminated unsuccessfully.'
  1104.  
  1105. *        Deallocate memory in various arrays.
  1106.  
  1107.         IF (REALDAT) THEN
  1108.           DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
  1109.           IF (IERR .NE. 0)
  1110.      +          STOP 'Problem attempting to deallocate.  Aborting ...'
  1111.         ELSE
  1112.           DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
  1113.           IF (IERR .NE. 0)
  1114.      +          STOP 'Problem attempting to deallocate.  Aborting ...'
  1115.         ENDIF
  1116.  
  1117.         CYCLE
  1118.  
  1119.       ELSE IF (TWOCHAN) THEN
  1120.  
  1121. *        Spectrum computations for two channels.
  1122.  
  1123.         IF (SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN) .NE. 0) STOP
  1124.      +        '                      Program terminated unsuccessfully.'
  1125.  
  1126. *        Deallocate memory in various arrays.
  1127.  
  1128.         IF (REALDAT) THEN
  1129.           DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
  1130.           IF (IERR .NE. 0)
  1131.      +          STOP 'Problem attempting to deallocate.  Aborting ...'
  1132.         ELSE
  1133.           DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
  1134.           IF (IERR .NE. 0)
  1135.      +          STOP 'Problem attempting to deallocate.  Aborting ...'
  1136.         ENDIF
  1137.  
  1138.       ENDIF
  1139.  
  1140.     ENDDO
  1141.  
  1142. *    Get the time just before ending.
  1143.  
  1144.     CALL GETTIM (IHR2,IMIN2,ISEC2,I100TH2)
  1145.  
  1146.     WRITE (*,'(/25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
  1147.      +    'Starting time : ',IHR,':',IMIN,':',ISEC,'.',I100TH
  1148.  
  1149.     WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
  1150.      +    'Ending time   : ',IHR2,':',IMIN2,':',ISEC2,'.',I100TH2
  1151.  
  1152. *    Calculate elapsed time - hundredths of a second.
  1153.  
  1154.     IF (I100TH2-I100TH .LT. 0) THEN
  1155.       I100TH2=I100TH2+100
  1156.       I100H=I100TH2-I100TH
  1157.       ISEC2=ISEC2-1
  1158.     ELSE
  1159.       I100H=I100TH2-I100TH
  1160.     ENDIF
  1161.  
  1162. *    Seconds.
  1163.  
  1164.     IF (ISEC2-ISEC .LT. 0) THEN
  1165.       ISEC2=ISEC2+60
  1166.       ISECH=ISEC2-ISEC
  1167.       IMIN2=IMIN2-1
  1168.     ELSE
  1169.       ISECH=ISEC2-ISEC
  1170.     ENDIF
  1171.  
  1172. *    Minutes.
  1173.  
  1174.     IF (IMIN2-IMIN .LT. 0) THEN
  1175.       IMIN2=IMIN2+60
  1176.       IMINH=IMIN2-IMIN
  1177.       IHR2=IHR2-1
  1178.     ELSE
  1179.       IMINH=IMIN2-IMIN
  1180.     ENDIF
  1181.  
  1182. *    Hours.
  1183.  
  1184.     IF (IHR2-IHR .LT. 0) THEN
  1185.       IHR2=IHR2+24
  1186.       IHRH=IHR2-IHR
  1187.     ELSE
  1188.       IHRH=IHR2-IHR
  1189.     ENDIF
  1190.  
  1191. *    Display elapsed time.
  1192.  
  1193.     WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
  1194.      +    'Elapsed time  : ',IHRH,':',IMINH,':',ISECH,'.',I100H
  1195.  
  1196.     WRITE (*,'( )')
  1197.     STOP '                        Program terminated successfully.'
  1198.  
  1199. *    I/O Error checking.
  1200.  
  1201. 390    WRITE (*,'(/17X,A/)')
  1202.      +    'Not enough space available for temporary files.'
  1203.  
  1204.     STOP '                      Program terminated unsuccessfully.'
  1205.  
  1206. 400    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1207.      +    'No such device or out of space !'
  1208.     STOP '                      Program terminated unsuccessfully.'
  1209.     END
  1210.  
  1211.     FUNCTION TRANS_ONE_CHAN (NDATA,FR,FI,B)
  1212.  
  1213.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1214.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1215.  
  1216.     INTEGER*2     NDATA[HUGE](*)
  1217.     REAL*4        FR[HUGE](*)
  1218.     REAL*4        FI[HUGE](*)
  1219.     REAL*4         B(*)
  1220.  
  1221. *    Initialization.
  1222.  
  1223.     TRANS_ONE_CHAN=0
  1224.     SUM1UT=0.0
  1225.     SUM2UT=0.0
  1226.  
  1227. *    Perform mean and RMS value calculations for each record
  1228. *    in this data file.
  1229.  
  1230.     DO IB=1,NUMREC
  1231.  
  1232. *      Fill the array with the current record.
  1233.  
  1234.       IF (REALDAT) THEN
  1235.         READ (NUMI) (FR(I),I=1,N)
  1236.         DO I=1,N
  1237.           FI(I)=FR(I)*FR(I)
  1238.         ENDDO
  1239.       ELSE
  1240.         READ (NUMI) (NDATA(I),I=1,N)
  1241.       ENDIF
  1242.  
  1243. *      Transform the data into voltages or velocities.
  1244. *      Calculate squared values for RMS calculations.
  1245.  
  1246.       IF (REALDAT) GO TO 10
  1247.       DO I=1,N
  1248.         V=FLOAT(NDATA(I))*VOFST1
  1249.         IF (VOLT) THEN
  1250.           FR(I)=V
  1251.           FI(I)=FR(I)*FR(I)
  1252.           CYCLE
  1253.         ENDIF
  1254.         FR(I)=B(5)
  1255.         DO J=4,1,-1
  1256.           FR(I)=FR(I)*V+B(J)
  1257.         ENDDO
  1258.         FI(I)=FR(I)*FR(I)
  1259.       ENDDO
  1260.  
  1261. *      Determine mean and RMS values for each record.
  1262.  
  1263. 10      SUM1U=ADDSUM(FR,N)
  1264.       SUM1UT=SUM1UT+SUM1U
  1265.       SUM2U=ADDSUM(FI,N)
  1266.       SUM2UT=SUM2UT+SUM2U
  1267.       UBAR=SUM1U/FN
  1268.       URMS=SQRT(SUM2U/FN-UBAR**2)
  1269.  
  1270. *      Display output for each record.
  1271.  
  1272.       IF (IB .EQ. 1) THEN
  1273.         WRITE (*,20) IB,UBAR,URMS
  1274. 20        FORMAT (16X,'Record ',I4.4,4X,'Ave = ',G11.5,
  1275.      +              2X,'Rms = ',G11.5)
  1276.       ELSE
  1277.         WRITE (*,30) IB,UBAR,URMS
  1278. 30        FORMAT ('+',15X,'Record ',I4.4,4X,'Ave = ',G11.5,
  1279.      +              2X,'Rms = ',G11.5)
  1280.       ENDIF
  1281.  
  1282. *      Save transformed values for this record.
  1283.  
  1284.       IF (.NOT. NOFIL)
  1285.      +        CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
  1286.       IF (RMEAN) THEN
  1287.         DO I=1,N
  1288.           FR(I)=FR(I)-UBAR
  1289.         ENDDO
  1290.       ENDIF
  1291.       IF (DTREND) CALL DTRD (FR,N,DELT)
  1292.       WRITE (NUMO,ERR=40) (FR(I), I=1,N)
  1293.  
  1294.     ENDDO
  1295.  
  1296. *    All records of this data file have been processed
  1297. *    and stored. Compute mean and rms values for this
  1298. *    data file and continue.
  1299.  
  1300.     TOTAL=FLOAT(NUMREC*N)
  1301.     AVEU=SUM1UT/TOTAL
  1302.     RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)
  1303.  
  1304. *    Display the results for this file.
  1305.  
  1306.     WRITE (*,'(/16X,A,1X,A,G11.5,2X,A,G11.5/)') 'File Totals : ',
  1307.      +         'Ave = ',AVEU,'Rms = ',RMSU
  1308.  
  1309.     RETURN
  1310.  
  1311. 40    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1312.      +                      'No such device or out of space !'
  1313.     TRANS_ONE_CHAN=1
  1314.     RETURN
  1315.     END
  1316.  
  1317.     FUNCTION TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)
  1318.  
  1319.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1320.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1321.  
  1322.     INTEGER*2     NDATA[HUGE](*)
  1323.     REAL*4        FR[HUGE](*)
  1324.     REAL*4        FR2[HUGE](*)
  1325.     REAL*4        FI[HUGE](*)
  1326.     REAL*4        FI2[HUGE](*)
  1327.     REAL*4        WIN[HUGE](*)
  1328.     REAL*4        B(*),B2(*)
  1329.  
  1330. *    Initialization.
  1331.  
  1332.     TRANS_TWO_CHAN=0
  1333.     SUM1UT=0.0
  1334.     SUM2UT=0.0
  1335.     SUM1VT=0.0
  1336.     SUM2VT=0.0
  1337.  
  1338.     N=N/2
  1339.     FN=FLOAT(N)
  1340.  
  1341. *    Perform mean and RMS value calculations for each record
  1342. *    in this data file.
  1343.  
  1344.     DO IB=1,NUMREC
  1345.  
  1346. *      Fill the array with the current record.
  1347.  
  1348.       IF (REALDAT) THEN
  1349.         READ (NUMI) (WIN(I),I=1,2*N)
  1350.         DO II=1,2*N,2
  1351.           I=(II+1)/2
  1352.           FR(I)=WIN(II)
  1353.           FR2(I)=WIN(II+1)
  1354.         ENDDO
  1355.         DO I=1,N
  1356.           FI(I)=FR(I)*FR(I)
  1357.           FI2(I)=FR2(I)*FR2(I)
  1358.         ENDDO
  1359.       ELSE
  1360.         READ (NUMI) (NDATA(I),I=1,2*N)
  1361.       ENDIF
  1362.  
  1363. *      Transform the data into voltages or velocities.
  1364. *      Calculate squared values for RMS calculations.
  1365.  
  1366.       IF (REALDAT) GO TO 20
  1367.       DO II=1,2*N,2
  1368.         I=(II+1)/2
  1369.         V=FLOAT(NDATA(II))*VOFST1
  1370.         V2=FLOAT(NDATA(II+1))*VOFST2
  1371.         IF (VOLT) THEN
  1372.           FR(I)=V
  1373.           FI(I)=FR(I)*FR(I)
  1374.           GO TO 10
  1375.         ENDIF
  1376.         FR(I)=B(5)
  1377.         DO J=4,1,-1
  1378.           FR(I)=FR(I)*V+B(J)
  1379.         ENDDO
  1380.         FI(I)=FR(I)*FR(I)
  1381. 10        IF (VOLT2) THEN
  1382.           FR2(I)=V2
  1383.           FI2(I)=FR2(I)*FR2(I)
  1384.           CYCLE
  1385.         ENDIF
  1386.         FR2(I)=B2(5)
  1387.         DO J=4,1,-1
  1388.           FR2(I)=FR2(I)*V2+B2(J)
  1389.         ENDDO
  1390.         FI2(I)=FR2(I)*FR2(I)
  1391.       ENDDO
  1392.  
  1393. *      Determine mean and RMS values for each record.
  1394.  
  1395. 20      SUM1U=ADDSUM(FR,N)
  1396.       SUM1UT=SUM1UT+SUM1U
  1397.       SUM2U=ADDSUM(FI,N)
  1398.       SUM2UT=SUM2UT+SUM2U
  1399.       UBAR=SUM1U/FN
  1400.       URMS=SQRT(SUM2U/FN-UBAR**2)
  1401.  
  1402.       SUM1V=ADDSUM(FR2,N)
  1403.       SUM1VT=SUM1VT+SUM1V
  1404.       SUM2V=ADDSUM(FI2,N)
  1405.       SUM2VT=SUM2VT+SUM2V
  1406.       VBAR=SUM1V/FN
  1407.       VRMS=SQRT(SUM2V/FN-VBAR**2)
  1408.  
  1409. *      Display output for each record.
  1410.  
  1411.       IF (IB .EQ. 1) THEN
  1412.         WRITE (*,30) IB,UBAR,URMS,VBAR,VRMS
  1413. 30        FORMAT (1X,'Rec ',I4.4,2X,'Ave= ',G11.5,1X,
  1414.      +       'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
  1415.       ELSE
  1416.         WRITE (*,40) IB,UBAR,URMS,VBAR,VRMS
  1417. 40        FORMAT ('+','Rec ',I4.4,2X,'Ave= ',G11.5,1X,
  1418.      +       'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
  1419.       ENDIF
  1420.  
  1421. *      Save transformed values for this record.
  1422.  
  1423.       IF (.NOT. NOFIL)
  1424.      +        CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
  1425.       IF (RMEAN) THEN
  1426.         DO I=1,N
  1427.           FR(I)=FR(I)-UBAR
  1428.         ENDDO
  1429.       ENDIF
  1430.       IF (DTREND) CALL DTRD (FR,N,DELT)
  1431.       WRITE (NUMO,ERR=50) (FR(I), I=1,N)
  1432.  
  1433.       IF (.NOT. NOFIL2)
  1434.      +        CALL FILTER (FTYP2,FC2,F12,F22,NS2,DELT,N,FR2)
  1435.       IF (RMEAN2) THEN
  1436.         DO I=1,N
  1437.           FR2(I)=FR2(I)-VBAR
  1438.         ENDDO
  1439.       ENDIF
  1440.       IF (DTREND2) CALL DTRD (FR2,N,DELT)
  1441.       WRITE (NUMO2,ERR=50) (FR2(I), I=1,N)
  1442.  
  1443.     ENDDO
  1444.  
  1445. *    All records of this data file have been processed
  1446. *    and stored. Compute mean and rms values for this
  1447. *    data file and continue.
  1448.  
  1449.     TOTAL=FLOAT(NUMREC*N)
  1450.     AVEU=SUM1UT/TOTAL
  1451.     RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)
  1452.  
  1453.     AVEV=SUM1VT/TOTAL
  1454.     RMSV=SQRT(SUM2VT/TOTAL-AVEV**2)
  1455.  
  1456. *    Display the results for this file.
  1457.  
  1458.     WRITE (*,'(/1X,A,2(A,G11.5,1X,A,G11.5,1X)/)')
  1459.      +         'File :    ','Ave= ',AVEU,'Rms= ',RMSU,
  1460.      +                    'Ave2= ',AVEV,'Rms2= ',RMSV
  1461.  
  1462.     RETURN
  1463.  
  1464. 50    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1465.      +                      'No such device or out of space !'
  1466.     TRANS_TWO_CHAN=1
  1467.     RETURN
  1468.     END
  1469.  
  1470.     FUNCTION WIN_SCALE (WIN)
  1471.  
  1472.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1473.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1474.  
  1475.     REAL*4    WIN[HUGE](*)
  1476.     DATA    A0,A1,A2,A3 /0.355768,0.487396,0.144232,0.012604/
  1477.     DATA    C0,C1,C2,C3 /0.3125,0.46875,0.1875,0.03125/
  1478.  
  1479. *    Initialization.
  1480.  
  1481.     WIN_SCALE=0
  1482.     S=0.0
  1483.     TNM1=FLOAT(N-1)
  1484.  
  1485. *    Compute the scale factor.
  1486.  
  1487.     DO I=1,N
  1488.       TIM1=FLOAT(I-1)
  1489.       IF (HANN)
  1490.      +      WIN(I)=(1.0-COS(PI2*TIM1/TNM1))/2.0
  1491.       IF (WELCH)
  1492.      +        WIN(I)=1.0-((TIM1-TNM1/2.0)/(TNM1/2.0))**2
  1493.       IF (PARZ)
  1494.      +        WIN(I)=1.0-ABS((TIM1-TNM1/2.0)/(TNM1/2.0))
  1495.       IF (TUKEY) THEN
  1496.         IF (I .LE. K .OR. I .GE. N-K)
  1497.      +        WIN(I)=(1.0-COS(10.0*PI*TIM1/TNM1))/2.0
  1498.         IF (I .GT. K .AND. I .LT. N-K) WIN(I)=1.0
  1499.       ENDIF
  1500.       IF (MINSL)
  1501.      +       WIN(I)=A0-A1*COS(PI2*TIM1/TNM1)
  1502.      +       +A2*COS(2.0*PI2*TIM1/TNM1)-A3*COS(3.0*PI2*TIM1/TNM1)
  1503.       IF (MXDEC)
  1504.      +       WIN(I)=C0-C1*COS(PI2*TIM1/TNM1)
  1505.      +       +C2*COS(2.0*PI2*TIM1/TNM1)-C3*COS(3.0*PI2*TIM1/TNM1)
  1506.       S=S+WIN(I)*WIN(I)
  1507.     ENDDO
  1508.  
  1509.     S=S/FLOAT(N)
  1510.     S=SQRT(1.0/S)
  1511.  
  1512.     RETURN
  1513.  
  1514. 40    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1515.      +                      'No such device or out of space !'
  1516.     WIN_SCALE=1
  1517.     RETURN
  1518.     END
  1519.  
  1520.     FUNCTION TAPER_DATA (FR,FR2,FI,FI2,WIN)
  1521.  
  1522.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1523.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1524.  
  1525.     REAL*4        FR[HUGE](*)
  1526.     REAL*4        FR2[HUGE](*)
  1527.     REAL*4        FI[HUGE](*)
  1528.     REAL*4        FI2[HUGE](*)
  1529.     REAL*4        WIN[HUGE](*)
  1530.  
  1531. *    Initialization.
  1532.  
  1533.     TAPER_DATA=0
  1534.     IBO=0
  1535.  
  1536. *    Do the tapering now.
  1537.  
  1538.     DO IB=1,NUMREC+1
  1539.  
  1540.       IF (IB .EQ. NUMREC+1 .AND. NOOVRLP) RETURN
  1541.  
  1542. *      Reposition input file if last time thru loop.
  1543.  
  1544.       IF (IB .EQ. NUMREC+1 .AND. OVRLP) THEN
  1545.         REWIND (NUMO)
  1546.         READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
  1547.         IF (TWOCHAN) THEN
  1548.           REWIND (NUMO2)
  1549.           READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
  1550.         ENDIF
  1551.       ENDIF
  1552.  
  1553. *      Get the next N points.
  1554.  
  1555.       READ (NUMO) (FR(I),I=1,N)
  1556.       IF (TWOCHAN) READ (NUMO2) (FR2(I),I=1,N)
  1557.       IF (IB .EQ. 1 .OR. NOOVRLP) GO TO 10
  1558.  
  1559. *      Operate on the first half of FR and store in FI.
  1560.  
  1561.       DO I=N2+1,N
  1562.         FI(I)=FR(I-N2)*S*WIN(I)
  1563.         IF (TWOCHAN) FI2(I)=FR2(I-N2)*S*WIN(I)
  1564.       ENDDO
  1565.  
  1566. *      Write contents of FI to output file.
  1567.  
  1568.       IBO=IBO+1
  1569.       WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO
  1570.  
  1571.       WRITE (NUMT,ERR=20) (FI(I),I=1,N)
  1572.       IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FI2(I),I=1,N)
  1573.       IF (IB .EQ. NUMREC+1) RETURN
  1574.  
  1575. *      Operate on all of FR.
  1576.  
  1577. 10      DO I=1,N
  1578.  
  1579. *        Operate on the second half of FR and store in FI.
  1580.  
  1581.         IF (I .LE. N2 .AND. OVRLP) THEN
  1582.           FI(I)=FR(N2+I)*S*WIN(I)
  1583.           IF (TWOCHAN) FI2(I)=FR2(N2+I)*S*WIN(I)
  1584.         ENDIF
  1585.  
  1586. *        Now operate on all of FR.
  1587.  
  1588.         FR(I)=FR(I)*S*WIN(I)
  1589.         IF (TWOCHAN) FR2(I)=FR2(I)*S*WIN(I)
  1590.  
  1591.       ENDDO
  1592.  
  1593. *      Write contents of FR to output file.
  1594.  
  1595.       IBO=IBO+1
  1596.       IF (IBO .EQ. 1) THEN
  1597.         WRITE (*,'(25X,''Tapering record '',I4.4)') IBO
  1598.       ELSE
  1599.         WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO
  1600.       ENDIF
  1601.  
  1602.       WRITE (NUMT,ERR=20) (FR(I),I=1,N)
  1603.       IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FR2(I),I=1,N)
  1604.  
  1605.     ENDDO
  1606.  
  1607.     RETURN
  1608.  
  1609. 20    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1610.      +                      'No such device or out of space !'
  1611.     TAPER_DATA=1
  1612.     RETURN
  1613.     END
  1614.  
  1615.     FUNCTION SPEC_ONE_CHAN (FR,FI,FR2,FI2)
  1616.  
  1617.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1618.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1619.  
  1620.     REAL*4        FR[HUGE](*)
  1621.     REAL*4        FI[HUGE](*)
  1622.     REAL*4        FR2[HUGE](*)
  1623.     REAL*4        FI2[HUGE](*)
  1624.  
  1625. *    Initialization.
  1626.  
  1627.     SPEC_ONE_CHAN=0
  1628.  
  1629. *    Spectrum computations for one channel.
  1630.  
  1631.     DO IB=1,NRECD2
  1632.  
  1633.       IF (IB .EQ. 1) THEN
  1634.         WRITE (*,'(/25X,''Transforming records '',I4.4,
  1635.      +                 '' and '',I4.4)') 2*IB-1,2*IB
  1636.       ELSE
  1637.         WRITE (*,'(''+'',24X,''Transforming records '',I4.4,
  1638.      +                 '' and '',I4.4)') 2*IB-1,2*IB
  1639.       ENDIF
  1640.  
  1641.       IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2)
  1642.      +      WRITE (*,'(/25X,A,I4.4,A)')
  1643.      +      'Record ',2*IB,' was discarded.'
  1644.  
  1645. *      Process two records simultaneously. This increases speed but
  1646. *      requires an even number of records in the one channel case.
  1647.  
  1648.       READ (INUM) (FR(I),I=1,N)
  1649.  
  1650.       IF (ODDNUM .AND. IB .EQ. NRECD2) THEN
  1651.         DO I=1,N
  1652.           FI(I)=FR(I)
  1653.         ENDDO
  1654.       ELSE
  1655.         READ (INUM) (FI(I),I=1,N)
  1656.       ENDIF
  1657.  
  1658.       CALL FFT (FR,FI,IOD)
  1659.  
  1660.       IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) FI(1)=0.0
  1661.  
  1662.       XMEAN=XMEAN+FR(1)
  1663.       YMEAN=YMEAN+FI(1)
  1664.  
  1665.           DO K=1,N2
  1666.  
  1667.         AK=FR(K+1)
  1668.         BK=FI(K+1)
  1669.         ANK=FR(NP1-K)
  1670.         BNK=FI(NP1-K)
  1671.  
  1672. *        Unscramble X and Y.
  1673.  
  1674.         XR=AK+ANK
  1675.         XI=BK-BNK
  1676.         YR=BNK+BK
  1677.         YI=ANK-AK
  1678.  
  1679.         IF (POWER) THEN
  1680.           B1K=XR*XR+XI*XI
  1681.           B2K=YR*YR+YI*YI
  1682.         ELSE
  1683.           B1K=SQRT(XR*XR+XI*XI)
  1684.           B2K=SQRT(YR*YR+YI*YI)
  1685.         ENDIF
  1686.  
  1687.         IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) B2K=0.0
  1688.  
  1689.         FR2(K)=FR2(K)+B1K
  1690.         FI2(K)=FI2(K)+B2K
  1691.  
  1692.       ENDDO
  1693.     ENDDO
  1694.  
  1695.     CLOSE (INUM,STATUS='DELETE',ERR=10)
  1696.  
  1697.     IB=IB-1        ! Decrement added count in loop variable
  1698.  
  1699. *    Put results into one buffer.
  1700.  
  1701.     XMEAN=XMEAN+YMEAN
  1702.     DO I=1,N
  1703.       FR2(I)=FR2(I)+FI2(I)
  1704.     ENDDO
  1705.     IB=IB*2
  1706.     IF (OVRLP .OR. ODDNUM) IB=IB-1    ! Remove bad record.
  1707.  
  1708. *    Finished.
  1709.  
  1710.     FN=FLOAT(N)
  1711.     TOTAL=FN*FLOAT(IB)
  1712.     XMEAN=XMEAN/TOTAL
  1713.     DT=DELT
  1714.     DELF=1.0/(FN*DT)
  1715.     IF (POWER) SCALE=DT/(2.0*TOTAL)
  1716.     IF (AMP) THEN
  1717.       SCALE=DT/(2.0*FLOAT(IB))
  1718.       XMEAN=XMEAN*TOTAL*SCALE*2.0
  1719.     ENDIF
  1720.  
  1721.     DO I=1,N
  1722.       FR2(I)=FR2(I)*SCALE
  1723.     ENDDO
  1724.  
  1725.     IF (POWER) XRMS=SQRT(DELF*ADDSUM(FR2,N2))
  1726.  
  1727. *    Put results into output [.PRN] file.
  1728.  
  1729.     IPD=INDEX (POUT,'.',.TRUE.)
  1730.     WRITE (*,'(/25X,''Writing '',A,'' now.'')') POUT(:IPD+3)
  1731.  
  1732.     OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=10)
  1733.  
  1734.     IF (POWER) THEN
  1735.       WRITE (NUMP,'(F11.5,2X,F11.5/)',ERR=10) XMEAN,XRMS
  1736.       DO I=1,N2
  1737.         IF (FR2(I) .LT. EPS) FR2(I)=EPS
  1738.         WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=10)
  1739.      +          DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I))
  1740.       ENDDO
  1741.     ELSE
  1742.       WRITE (NUMP,'(G15.5,2X,G15.5,2X,G15.5)',ERR=10)
  1743.      +           EPS,ALOG10(EPS),XMEAN
  1744.       WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5)',ERR=10)
  1745.      +          (DELF*I,ALOG10(DELF*I),FR2(I), I=1,N2)
  1746.     ENDIF
  1747.  
  1748.     CLOSE (NUMP,STATUS='KEEP',ERR=10)
  1749.  
  1750. *    Display results on the screen.
  1751.  
  1752.     IF (POWER) THEN
  1753.        WRITE (*,'(/29X,A,10X,A/21X,1P,2G15.5)')
  1754.      +                                    'mean','rms',XMEAN,XRMS
  1755.        PERR = 1.0/SQRT(FLOAT(NREC))  ! use numrec before overlap
  1756.        IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
  1757.        WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
  1758.     ENDIF
  1759.  
  1760.     RETURN
  1761.  
  1762. 10    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  1763.      +                      'No such device or out of space !'
  1764.     SPEC_ONE_CHAN=1
  1765.     RETURN
  1766.     END
  1767.  
  1768.     FUNCTION SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN)
  1769.  
  1770.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  1771.     INCLUDE        'SPECTRUM.FN'    ! Common statements
  1772.  
  1773.     REAL*4        FR[HUGE](*)
  1774.     REAL*4        FI[HUGE](*)
  1775.     REAL*4        FR2[HUGE](*)
  1776.     REAL*4        FI2[HUGE](*)
  1777.     REAL*4        WIN[HUGE](*)
  1778.  
  1779. *    Initialization.
  1780.  
  1781.     SPEC_TWO_CHAN=0
  1782.  
  1783. *    Spectrum computations for two channels.
  1784.  
  1785.     DO IB=1,NUMREC
  1786.  
  1787.       IF (IB .EQ. 1) THEN
  1788.         WRITE (*,'(/25X,''Transforming record '',I4.4)') IB
  1789.       ELSE
  1790.         WRITE (*,'(''+'',24X,''Transforming record '',I4.4)') IB
  1791.       ENDIF
  1792.  
  1793.       IF (OVRLP .AND. IB .EQ. NUMREC)
  1794.      +      WRITE (*,'(/25X,A,I4.4,A)')
  1795.      +      'Record ',IB+1, ' was discarded for each channel.'
  1796.  
  1797.       READ (INUM) (FR(I),I=1,N)
  1798.       READ (INUM2) (FI(I),I=1,N)
  1799.  
  1800.       CALL FFT (FR,FI,IOD)
  1801.  
  1802.       XMEAN=XMEAN+FR(1)
  1803.       YMEAN=YMEAN+FI(1)
  1804.  
  1805.           DO K=1,N2
  1806.  
  1807.         AK=FR(K+1)
  1808.         BK=FI(K+1)
  1809.         ANK=FR(NP1-K)
  1810.         BNK=FI(NP1-K)
  1811.  
  1812. *        Unscramble X and Y.
  1813.  
  1814.         XR=AK+ANK
  1815.         XI=BK-BNK
  1816.         YR=BNK+BK
  1817.         YI=ANK-AK
  1818.  
  1819.         IF (POWER) THEN
  1820.           B1K=XR*XR+XI*XI
  1821.           B2K=YR*YR+YI*YI
  1822.         ELSE IF (CROSS) THEN
  1823.           B1K=XR*YR+XI*YI
  1824.           B2K=XI*YR-XR*YI
  1825.           IF (COH) THEN
  1826.             B3K=XR*XR+XI*XI
  1827.             B4K=YR*YR+YI*YI
  1828.           ENDIF
  1829.         ELSE
  1830.           B1K=SQRT(XR*XR+XI*XI)
  1831.           B2K=SQRT(YR*YR+YI*YI)
  1832.         ENDIF
  1833.  
  1834.         FR2(K)=FR2(K)+B1K
  1835.         FI2(K)=FI2(K)+B2K
  1836.  
  1837.         IF (COH) THEN
  1838.           WIN(K)=WIN(K)+B3K
  1839.           WIN(N2+K)=WIN(N2+K)+B4K
  1840.         ENDIF
  1841.  
  1842.       ENDDO
  1843.     ENDDO
  1844.  
  1845.     CLOSE (INUM,STATUS='DELETE',ERR=20)
  1846.     CLOSE (INUM2,STATUS='DELETE',ERR=20)
  1847.  
  1848.     IB=IB-1        ! Decrement added count in loop variable
  1849.  
  1850. *    Finished.
  1851.  
  1852.     FN=FLOAT(N)
  1853.     TOTAL=FN*FLOAT(IB)
  1854.     XMEAN=XMEAN/TOTAL
  1855.     YMEAN=YMEAN/TOTAL
  1856.     DT=DELT
  1857.     DELF=1.0/(FN*DT)
  1858.     IF (POWER .OR. CROSS) SCALE=DT/(2.0*TOTAL)
  1859.     IF (AMP) THEN
  1860.       SCALE=DT/(2.0*FLOAT(IB))
  1861.       XMEAN=XMEAN*TOTAL*SCALE*2.0
  1862.       YMEAN=YMEAN*TOTAL*SCALE*2.0
  1863.     ENDIF
  1864.  
  1865.     DO I=1,N
  1866.       FR2(I)=FR2(I)*SCALE
  1867.       FI2(I)=FI2(I)*SCALE
  1868.     ENDDO
  1869.  
  1870.     IF (COH) THEN
  1871.       DO I=1,N
  1872.         WIN(I)=WIN(I)*SCALE
  1873.       ENDDO
  1874.     ENDIF
  1875.  
  1876.     IF (POWER) THEN
  1877.       XRMS=SQRT(DELF*ADDSUM(FR2,N2))
  1878.       YRMS=SQRT(DELF*ADDSUM(FI2,N2))
  1879.     ENDIF
  1880.  
  1881.     IF (CROSS) THEN
  1882.       DEGRAD=360.0/PI2
  1883.       S2=SQRT(2.0)
  1884.       NR=NREC                ! use numrec before overlap
  1885.       SNR=SQRT(FLOAT(NR))
  1886.       S2NR=SQRT(FLOAT(NR)*2.0)
  1887.       RXY0=DELF*ADDSUM(FR2,N2)
  1888.       DO I=1,N2
  1889.         FII=FLOAT(I)
  1890.         GXYMAG=SQRT(FR2(I)**2+FI2(I)**2)
  1891.         PHASE=ATAN2(FI2(I),FR2(I))
  1892.         FR2(I)=GXYMAG
  1893.         FI2(I)=PHASE
  1894.         FI2(N2+I)=PHASE
  1895.  
  1896. *        If coherence is calculated then compute errors associated
  1897. *        with cross spectrum estimates.
  1898.  
  1899.         IF (COH) THEN
  1900.           GAMMA=FR2(I)**2/(WIN(I)*WIN(N2+I))
  1901.           WIN(I)=GAMMA
  1902.           FR(I)=1.0/(SQRT(WIN(I))*SNR)                  ! e [ |Gxy(f)| ]
  1903.           FR(N2+I)=SQRT(1.0-WIN(I))/(SQRT(WIN(I))*S2NR) ! s.d. [ Oxy(f) ]
  1904.           FR(N2+I)=FR(N2+I)*DEGRAD                      ! s.d. [O] in deg
  1905.           FI(I)=S2*(1.0-WIN(I))/(SQRT(WIN(I))*SNR)      ! e [ gxy2(f) ]
  1906.           IF (TAPER .AND. NOOVRLP) THEN
  1907.             FR(I)=FR(I)*S2
  1908.             FR(N2+I)=FR(N2+I)*S2
  1909.             FI(I)=FI(I)*S2
  1910.           ENDIF
  1911.         ENDIF
  1912.       ENDDO
  1913.  
  1914. *      Apply a phase shift to straighten out Theta only if
  1915. *      coherence is calculated -- need s.d. estimate.
  1916.  
  1917.       IF (COH) THEN
  1918.  
  1919.         NSD2=(NSD-1)/2
  1920.  
  1921. *        Check on proper relationship of NPTS and NSD.
  1922.  
  1923.         IF (NSD2+1 .GT. NPTS) THEN
  1924.           WRITE (*,'(/21X,A,A/21X,A)') 'NSD too big or NPTS',
  1925.      +        ' too small.','No adjustment of phase angle performed.'
  1926.           GO TO 10
  1927.         ENDIF
  1928.  
  1929. *        Get first NPTS points. Pass this to fit subroutine.
  1930.  
  1931.         DO I=1,NPTS
  1932.           X(I)=DELF*I
  1933.           Y(I)=FI2(N2+I)
  1934.         ENDDO
  1935.  
  1936. *        Sum first NSD points to check on average s.d. of theta.
  1937.  
  1938.         SDSUM=0.0
  1939.         DO I=1,NSD
  1940.           SDSUM=SDSUM+FR(N2+NPTS-NSD2-1+I)
  1941.         ENDDO
  1942.  
  1943.         DO I=NPTS+1,N2
  1944.  
  1945. *          Check on moving average of s.d. of theta.
  1946. *          If too high -- quit.
  1947.  
  1948.           SDSUM=SDSUM+FR(N2+I+NSD2)
  1949.           SDSUM=SDSUM-FR(N2+I-NSD2-1)
  1950.           SDAVG=SDSUM/NSD
  1951.           IF (SDAVG .GT. SDLIM) GO TO 10
  1952.  
  1953. *          Perform linear least squares fit on NPTS points before
  1954. *          the current point. Guess value of current point from fit.
  1955.  
  1956.           CALL LSF (X,Y,NPTS,XA1,XB2)
  1957.           TEST=XA1+XB2*DELF*I
  1958.  
  1959. *          XLO & XHI are current point guess +/- pi.
  1960.  
  1961.           XLO=TEST-PI
  1962.           XHI=TEST+PI
  1963.  
  1964. *          If actual value is .GT. pi away from guessed value then
  1965. *          add 2*pi to or subtract 2*pi from actual value as required.
  1966.  
  1967.           DO WHILE (FI2(N2+I) .LT. XLO)
  1968.             FI2(N2+I)=FI2(N2+I)+PI2
  1969.           ENDDO
  1970.  
  1971.           DO WHILE (FI2(N2+I) .GT. XHI)
  1972.             FI2(N2+I)=FI2(N2+I)-PI2
  1973.           ENDDO
  1974.  
  1975. *          Prepare set of NPTS points before next current point.
  1976.  
  1977.           DO J=1,NPTS
  1978.             X(J)=DELF*(I-NPTS+J)
  1979.             Y(J)=FI2(N2+I-NPTS+J)
  1980.           ENDDO
  1981.  
  1982.         ENDDO
  1983.       ENDIF
  1984.     ENDIF
  1985.  
  1986. *    Put results into output [.PRN] file.
  1987.  
  1988. 10    IPD=INDEX (POUT,'.',.TRUE.)
  1989.     WRITE (*,'(/25X,''Writing '',A,'' now.''/)') POUT(:IPD+3)
  1990.  
  1991.     OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=20)
  1992.  
  1993.     IF (POWER) THEN
  1994.       WRITE (NUMP,'(4(F11.5,2X))',ERR=20) XMEAN,XRMS,YMEAN,YRMS
  1995.       WRITE (NUMP,'( )',ERR=20)
  1996.       DO I=1,N2
  1997.         IF (FR2(I) .LT. EPS) FR2(I)=EPS
  1998.         IF (FI2(I) .LT. EPS) FI2(I)=EPS
  1999.         WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,
  2000.      +        F11.5,1X,G15.5,1X,F11.5)',ERR=20)
  2001.      +        DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),
  2002.      +        FI2(I),ALOG10(FI2(I))
  2003.       ENDDO
  2004.     ELSE IF (CROSS) THEN
  2005.       IF (NOCOH) THEN
  2006.         WRITE (NUMP,'(G15.5)',ERR=20) RXY0
  2007.         WRITE (NUMP,'( )',ERR=20)
  2008.         WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,
  2009.      +        1X,F11.5)',ERR=20) (DELF*I,ALOG10(DELF*I),
  2010.      +        FR2(I),ALOG10(FR2(I)),FI2(I),I=1,N2)
  2011.       ELSE IF (COH) THEN
  2012.         WRITE (NUMP,'(G15.5)',ERR=20) RXY0
  2013.         WRITE (NUMP,'( )',ERR=20)
  2014.         WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,1X,
  2015.      +        F11.5,1X,F11.5,1X,G15.5)',ERR=20) (DELF*I,
  2016.      +        ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),FI2(I),FI2(N2+I),
  2017.      +        WIN(I),I=1,N2)
  2018.       ENDIF
  2019.     ELSE
  2020.       WRITE (NUMP,'(4(G15.5,2X))',ERR=20)
  2021.      +      EPS,ALOG10(EPS),XMEAN,YMEAN
  2022.       WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=20)
  2023.      +      (DELF*I,ALOG10(DELF*I),FR2(I),FI2(I),I=1,N2)
  2024.     ENDIF
  2025.  
  2026.     CLOSE (NUMP,STATUS='KEEP',ERR=20)
  2027.  
  2028. *    Display results on the screen.
  2029.  
  2030.     IF (POWER) THEN
  2031.        WRITE (*,'(29X,A,10X,A,2(/21X,1P,2G15.5))')
  2032.      +       'mean','rms',XMEAN,XRMS,YMEAN,YRMS
  2033.        PERR = 1.0/SQRT(FLOAT(NREC))  ! use numrec before overlap
  2034.        IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
  2035.        WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
  2036.     ENDIF
  2037.  
  2038.     IF (CROSS) WRITE (*,'(21X,A,1P,G15.5)')
  2039.      +    'Rxy(0) = E [x(t) y(t)] = ',RXY0
  2040.  
  2041. *    If coherence was calculated then store
  2042. *    error estimates into separate error file.
  2043.  
  2044.     IF (COH) THEN
  2045.       IPD=INDEX (EOUT,'.',.TRUE.)
  2046.       WRITE (*,'(/25X,''Writing '',A,'' now.'')') EOUT(:IPD+3)
  2047.  
  2048.       OPEN (NUME,FILE=EOUT,STATUS='UNKNOWN',ERR=20)
  2049.  
  2050.       WRITE (NUME,'(5(F11.5,2X))',ERR=20) (DELF*I,
  2051.      +      ALOG10(DELF*I),FR(I),FR(N2+I),FI(I),I=1,N2)
  2052.  
  2053.       CLOSE (NUME,STATUS='KEEP',ERR=20)
  2054.     ENDIF
  2055.  
  2056.     RETURN
  2057.  
  2058. 20    WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
  2059.      +    'No such device or out of space !'
  2060.     SPEC_TWO_CHAN=1
  2061.     RETURN
  2062.     END
  2063.  
  2064. ********************************************************
  2065.     SUBROUTINE FFT (FR,FI,K)
  2066. *                                                      *    
  2067. *    Fast Fourier Transform of a                    *
  2068. *    Complex Time Series                            *
  2069. *    Features : time decomposition (odd-even)       *
  2070. *             : input bit reversal                  *
  2071. *             : computation 'in place'              *
  2072. *             : N= number of points                 *
  2073. *             : FR+j*FI - complex series            *
  2074. *    Source   : Stearns 1975                        *
  2075. ********************************************************
  2076.  
  2077.     IMPLICIT    INTEGER*4 (I-N)
  2078.     REAL*4        FR[HUGE](*)
  2079.     REAL*4        FI[HUGE](*)
  2080.  
  2081.     N=2**K
  2082.     MR=0
  2083.     NN=N-1
  2084.  
  2085. *    Input bit reversal.
  2086.  
  2087.     DO M=1,NN
  2088.       L=N
  2089. 10      L=L/2
  2090.       IF(MR+L .GT. NN) GO TO 10
  2091.       MR=MOD(MR,L)+L
  2092.       IF (MR.LE. M) CYCLE
  2093.       TR=FR(M+1)
  2094.       FR(M+1)=FR(MR+1)
  2095.       FR(MR+1)=TR
  2096.       TI=FI(M+1)
  2097.       FI(M+1)=FI(MR+1)
  2098.       FI(MR+1)=TI
  2099.     ENDDO
  2100.  
  2101. *    Compute the coefficients.
  2102.  
  2103.     L=1
  2104.     PI=2.0*ASIN(1.0)
  2105. 20    IF (L.GE.N) RETURN
  2106.     ISTEP=2*L
  2107.     EL=L
  2108.  
  2109.     DO M=1,L
  2110.       A=PI*FLOAT(1-M)/EL
  2111.       WR=COS(A)
  2112.       WI=SIN(A)
  2113.       DO I=M,N,ISTEP
  2114.         J=I+L
  2115.         TR=WR*FR(J)-WI*FI(J)
  2116.         TI=WR*FI(J)+WI*FR(J)
  2117.         FR(J)=FR(I)-TR
  2118.         FI(J)=FI(I)-TI
  2119.         FR(I)=FR(I)+TR
  2120.         FI(I)=FI(I)+TI
  2121.       ENDDO
  2122.     ENDDO
  2123.  
  2124.     L=ISTEP
  2125.     GO TO 20
  2126.     END
  2127.  
  2128.     SUBROUTINE RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)
  2129.  
  2130. *    Routine to read config file - spectrum.cfg
  2131.  
  2132.     PARAMETER    (NUM=2)
  2133.     LOGICAL*1    ROOT1,ROOT2,SUFX
  2134.     CHARACTER*1    FIRST
  2135.     CHARACTER*4    SUFFIX
  2136.     CHARACTER*72    FLINE,PTH,TPTH
  2137.  
  2138. *    File exists. Get various paths.
  2139.  
  2140.     KOUNT=0
  2141.     OPEN (NUM,FILE='SPECTRUM.CFG',STATUS='OLD')
  2142.       IF (EOF(NUM)) GO TO 30
  2143. 10      READ (NUM,'(A)',ERR=30) FLINE
  2144.  
  2145. *      Check for a comment line.
  2146.  
  2147.       IF (FLINE (1:1) .EQ. '*') GO TO 10
  2148.       KOUNT=KOUNT+1
  2149.  
  2150. *      Convert to upper case characters if necessary.
  2151.  
  2152.       DO J=72,1,-1
  2153.         FIRST=FLINE(J:J)
  2154.         IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
  2155.           IHOLD=ICHAR(FIRST)-32
  2156.           FIRST=CHAR(IHOLD)
  2157.           FLINE(J:J)=FIRST
  2158.         ENDIF
  2159.       ENDDO
  2160.  
  2161. *      Assign appropriate paths.
  2162.  
  2163.       IF (KOUNT .EQ. 1) THEN
  2164.         IF (FLINE(:4) .EQ. 'ROOT') THEN
  2165.           ROOT1=.TRUE.
  2166.         ELSE
  2167.           ROOT1=.FALSE.
  2168.           ISL=INDEX (FLINE,'\',.TRUE.)
  2169.           PTH=FLINE(:ISL)
  2170.         ENDIF
  2171.       ELSE IF (KOUNT .EQ. 2) THEN
  2172.         IF (FLINE(:4) .EQ. 'ROOT') THEN
  2173.           ROOT2=.TRUE.
  2174.         ELSE
  2175.           ROOT2=.FALSE.
  2176.           ISL=INDEX (FLINE,'\',.TRUE.)
  2177.           TPTH=FLINE(:ISL)
  2178.         ENDIF
  2179.       ELSE IF (KOUNT .EQ. 3) THEN
  2180.         SUFX=.TRUE.
  2181.         SUFFIX=FLINE(:4)
  2182.       ENDIF
  2183.       GO TO 10
  2184.  
  2185. 30    CLOSE (NUM)
  2186.     RETURN
  2187.     END
  2188.  
  2189.     SUBROUTINE LSF(X,Y,N,A,B)
  2190.  
  2191.     IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  2192.     REAL*4 X(N),Y(N)
  2193.  
  2194. *    Initialization.
  2195.  
  2196.     XS=0.0
  2197.     X2S=0.0
  2198.     YS=0.0
  2199.     XYS=0.0
  2200.  
  2201. *    Compute sums.
  2202.  
  2203.     DO I=1,N
  2204.       XS=XS+X(I)
  2205.       X2S=X2S+X(I)*X(I)
  2206.       YS=YS+Y(I)
  2207.       XYS=XYS+X(I)*Y(I)
  2208.     ENDDO
  2209.  
  2210. *    Compute coefficients and guess next value.
  2211.  
  2212.     DEN=N*X2S-XS*XS
  2213.     A=(X2S*YS-XS*XYS)/DEN
  2214.     B=(N*XYS-XS*YS)/DEN
  2215.  
  2216.     RETURN
  2217.     END
  2218.  
  2219.     FUNCTION ADDSUM (Y,N)
  2220.  
  2221. *    SUMMING of elements of real*4 array Y. This function uses
  2222. *    a tree summming process to minimize truncation error for
  2223. *    large N(>500). A Real*4 and Integer*2 buffer of size log2(n)
  2224. *    is required. This version received from Dr. Anuvat Sirivat.
  2225.  
  2226.     REAL*4      Y[HUGE](*)
  2227.     REAL*4      B(16),ADDSUM
  2228.     INTEGER*2   W(16)
  2229.     INTEGER*4   N
  2230.  
  2231.     K=1
  2232.     W(1)=0
  2233.     B(1)=0.0
  2234.  
  2235.     DO I=1,N
  2236.       B(K)=B(K)+Y(I)
  2237.       W(K)=W(K)+1
  2238. 10      K1=K-1
  2239.  
  2240.       IF(K1.LE.0) GO TO 20
  2241.       IF(W(K1).GT.W(K)) GO TO 20
  2242.  
  2243.       B(K1)=B(K1)+B(K)
  2244.       W(K1)=W(K1)+W(K)
  2245.       K=K-1
  2246.       GO TO 10
  2247.  
  2248. 20      IF(W(K).LT.2) CYCLE
  2249.  
  2250.       K=K+1
  2251.       W(K)=0
  2252.       B(K)=0.0
  2253.     ENDDO
  2254.  
  2255. 30    K2=K/2
  2256.     IF(K2.EQ.0) GO TO 40
  2257.  
  2258.     J=1
  2259.     DO I=1,K2
  2260.       B(I)=B(J)+B(J+1)
  2261.       J=J+2
  2262.     ENDDO
  2263.  
  2264.         IF(J.EQ.K) B(1)=B(1)+B(K)
  2265.         K=K2
  2266.     GO TO 30
  2267.  
  2268. 40    ADDSUM=B(1)
  2269.     RETURN
  2270.     END
  2271.  
  2272.     SUBROUTINE DTRD (FR,N,DELT)
  2273.  
  2274. *    This routine will remove a linear mean trend from the
  2275. *    original data by means of least squares fit. This is a
  2276. *    modified version of a routine received from Dr. Sirivat.
  2277.  
  2278.     IMPLICIT REAL*4 (A-H,O-Z), INTEGER*4 (I-N)
  2279.     INTEGER*4    N
  2280.     REAL*4        DELT
  2281.     REAL*4        FR[HUGE](*)
  2282.  
  2283. *    Initialization.
  2284.  
  2285.     FN=FLOAT(N)
  2286.     YSR=0.0
  2287.     XYSR=0.0
  2288.  
  2289. *    Compute required sums for least squares fit.
  2290.  
  2291.     DO I=1,N
  2292.       FII=FLOAT(I)
  2293.       YSR=YSR+FR(I)
  2294.       XYSR=XYSR+FII*FR(I)
  2295.     ENDDO
  2296.  
  2297.     XS=FN*(FN+1.0)/2.0*DELT
  2298.     X2S=FN*(FN+1.0)*(2.0*FN+1.0)/6.0*DELT*DELT
  2299.     XYSR=XYSR*DELT
  2300.  
  2301.     DER=FN*X2S-XS*XS
  2302.  
  2303. *    Compute coefficients and remove trend.
  2304.  
  2305.     B0R=(YSR*X2S-XYSR*XS)/DER
  2306.     B1R=(FN*XYSR-XS*YSR)/DER
  2307.  
  2308.      DO I=1,N
  2309.       FII=FLOAT(I)
  2310.       FR(I)=FR(I)-B0R-B1R*FII*DELT
  2311.     ENDDO
  2312.  
  2313.     RETURN
  2314.     END
  2315.  
  2316.     SUBROUTINE FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
  2317.  
  2318. *    Routine to filter a data file.  (FILTER.FOR)
  2319.  
  2320. *    This routine will filter a data file using any one of
  2321. *    four different types of IIR recursive digital filters:
  2322. *    Lowpass, Highpass, Bandpass, and Bandstop. The routines
  2323. *    are modified versions of those given in Stearns, 1975.
  2324.  
  2325.     PARAMETER     (NSMAX=20)
  2326.     LOGICAL*1     LPASS,HPASS,BPASS,BSTOP
  2327.     INTEGER*2     NS
  2328.     INTEGER*4     N,J
  2329.     REAL*4         FR[HUGE](*)
  2330.     REAL*4         X(NSMAX+1,5)
  2331.     REAL*4         A(NSMAX),B(NSMAX),C(NSMAX),D(NSMAX),E(NSMAX),G
  2332.     CHARACTER*1     FTYPE
  2333.  
  2334. *    Initialization.
  2335.  
  2336.     X=0.0
  2337.     G=0.0
  2338.  
  2339.     LPASS=(FTYPE .EQ. 'L')
  2340.     HPASS=(FTYPE .EQ. 'H')
  2341.     BPASS=(FTYPE .EQ. 'B')
  2342.     BSTOP=(FTYPE .EQ. 'S')
  2343.  
  2344. *    Design the desired filter.
  2345.  
  2346.     IF (LPASS .OR. HPASS) KHI=3
  2347.     IF (BPASS .OR. BSTOP) KHI=5
  2348.  
  2349.     IF (LPASS) CALL LP (FC,DELT,NS,A,B,C)
  2350.     IF (HPASS) CALL HP (FC,DELT,NS,A,B,C)
  2351.     IF (BPASS) CALL BP (F1,F2,DELT,NS,A,B,C,D,E)
  2352.     IF (BSTOP) CALL BS (F1,F2,DELT,NS,A,B,C,D,E,G)
  2353.     G2=2.0*G
  2354.     GG=2.0+G*G
  2355.  
  2356. *    Filter the array FR(J).
  2357.  
  2358.     DO J=1,N
  2359.  
  2360. *      Go thru NS sections of filter.
  2361.  
  2362.       IF (LPASS) THEN
  2363.         X(1,3)=FR(J)
  2364.         DO I=1,NS
  2365.           TEMP=A(I)*(X(I,3)+2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
  2366.           X(I+1,3)=TEMP-C(I)*X(I+1,1)
  2367.         ENDDO
  2368.       ELSE IF (HPASS) THEN
  2369.         X(1,3)=FR(J)
  2370.         DO I=1,NS
  2371.           TEMP=A(I)*(X(I,3)-2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
  2372.           X(I+1,3)=TEMP-C(I)*X(I+1,1)
  2373.         ENDDO
  2374.       ELSE IF (BPASS) THEN
  2375.         X(1,5)=FR(J)
  2376.         DO I=1,NS
  2377.           TEMP=A(I)*(X(I,5)-2.0*X(I,3)+X(I,1))-B(I)*X(I+1,4)
  2378.           X(I+1,5)=TEMP-C(I)*X(I+1,3)-D(I)*X(I+1,2)-E(I)*X(I+1,1)
  2379.         ENDDO
  2380.       ELSE IF (BSTOP) THEN
  2381.         X(1,5)=FR(J)
  2382.         DO I=1,NS
  2383.           TEMP=A(I)*(X(I,5)+G2*X(I,4)+GG*X(I,3)+G2*X(I,2)+X(I,1))
  2384.           X(I+1,5)=TEMP-B(I)*X(I+1,4)-C(I)*X(I+1,3)-D(I)*X(I+1,2)
  2385.      +                 -E(I)*X(I+1,1)
  2386.         ENDDO
  2387.       ENDIF
  2388.  
  2389. *      Push down past values in filter.
  2390.  
  2391.       DO I=1,NS+1
  2392.         DO K=1,KHI-1
  2393.           X(I,K)=X(I,K+1)
  2394.         ENDDO
  2395.       ENDDO
  2396.  
  2397. *      Set FR(J) equal to the output of the filter and continue.
  2398.  
  2399.       FR(J)=X(NS+1,KHI)
  2400.  
  2401.     ENDDO
  2402.  
  2403.     RETURN
  2404.     END
  2405.  
  2406.     SUBROUTINE LP (FC,T,NS,A,B,C)
  2407.  
  2408. *    Lowpass Butterworth Digital Filter 
  2409. *    Design Subroutine (LP.FOR)
  2410.  
  2411.     INTEGER*2     NS
  2412.     REAL*4         A(*),B(*),C(*)
  2413.  
  2414.     PI=2.0*ASIN(1.0)
  2415.  
  2416.     WCP=SIN(FC*PI*T)/COS(FC*PI*T)
  2417.  
  2418.     DO K=1,NS
  2419.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  2420.       X=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
  2421.       A(K)=WCP*WCP*X
  2422.       B(K)=2.0*(WCP*WCP-1.0)*X
  2423.       C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*X
  2424.     ENDDO
  2425.  
  2426.     RETURN
  2427.     END
  2428.  
  2429.     SUBROUTINE HP (FC,T,NS,A,B,C)
  2430.  
  2431. *    Highpass Butterworth Digital Filter 
  2432. *    Design Subroutine (HP.FOR)
  2433.  
  2434.     INTEGER*2     NS
  2435.     REAL*4         A(*),B(*),C(*)
  2436.  
  2437.     PI=2.0*ASIN(1.0)
  2438.  
  2439.     WCP=SIN(FC*PI*T)/COS(FC*PI*T)
  2440.  
  2441.     DO K=1,NS
  2442.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  2443.       A(K)=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
  2444.       B(K)=2.0*(WCP*WCP-1.0)*A(K)
  2445.       C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*A(K)
  2446.     ENDDO
  2447.  
  2448.     RETURN
  2449.     END
  2450.  
  2451.     SUBROUTINE BP (F1,F2,T,NS,A,B,C,D,E)
  2452.  
  2453. *    Bandpass Butterworth Digital Filter 
  2454. *    Design Subroutine (BP.FOR)
  2455.  
  2456.     INTEGER*2     NS
  2457.     REAL*4         A(*),B(*),C(*),D(*),E(*)
  2458.  
  2459.     PI=2.0*ASIN(1.0)
  2460.  
  2461.     W1=SIN(F1*PI*T)/COS(F1*PI*T)
  2462.     W2=SIN(F2*PI*T)/COS(F2*PI*T)
  2463.     WC=W2-W1
  2464.     Q=WC*WC+2.0*W1*W2
  2465.     S=W1*W1*W2*W2
  2466.  
  2467.     DO K=1,NS
  2468.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  2469.       P=-2.0*WC*CS
  2470.       R=P*W1*W2
  2471.       X=1.0+P+Q+R+S
  2472.       A(K)=WC*WC/X
  2473.       B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
  2474.       C(K)=(6.0-2.0*Q+6.0*S)/X
  2475.       D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
  2476.       E(K)=(1.0-P+Q-R+S)/X
  2477.     ENDDO
  2478.  
  2479.     RETURN
  2480.     END
  2481.  
  2482.     SUBROUTINE BS (F1,F2,T,NS,A,B,C,D,E,G)
  2483.  
  2484. *    Bandstop Butterworth Digital Filter 
  2485. *    Design Subroutine (BS.FOR)
  2486.  
  2487.     INTEGER*2     NS
  2488.     REAL*4         A(*),B(*),C(*),D(*),E(*),G
  2489.  
  2490.     PI=2.0*ASIN(1.0)
  2491.  
  2492.     W1=SIN(F1*PI*T)/COS(F1*PI*T)
  2493.     W2=SIN(F2*PI*T)/COS(F2*PI*T)
  2494.     WC=W2-W1
  2495.     XK=1.0+W1*W2
  2496.     G=(2.0*XK-4.0)/XK
  2497.     Q=WC*WC+2.0*W1*W2
  2498.     S=W1*W1*W2*W2
  2499.  
  2500.     DO K=1,NS
  2501.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  2502.       P=-2.0*WC*CS
  2503.       R=P*W1*W2
  2504.       X=1.0+P+Q+R+S
  2505.       A(K)=XK*XK/X
  2506.       B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
  2507.       C(K)=(6.0-2.0*Q+6.0*S)/X
  2508.       D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
  2509.       E(K)=(1.0-P+Q-R+S)/X
  2510.     ENDDO
  2511.  
  2512.     RETURN
  2513.     END
  2514.