*	SPECTRUM.FOR  V 2.0

*	Routine to compute Amplitude, Auto and Cross-Spectral
*	Density Functions via Fast Fourier Transforms

*	David E. Hess
*	Fluid Flow Group - Process Measurements Division
*	Chemical Science and Technology Laboratory
*	National Institute of Standards and Technology
*	April 15, 1992

*	References: 1. Bendat, J.S. and Piersol, A.G. 1986 Random Data,
*	               2nd Ed., Wiley, New York.
*		    2. Harris, F.J. 1978 On the use of windows for
*		       harmonic analysis with the discrete fourier
*		       transform, Proceedings of the IEEE, 66, 51-83.
*		    3. Nuttall, A.H. 1981 Some windows with very good
*		       sidelobe behavior, IEEE Transactions on Acoustics,
*		       Speech and Signal Processing, ASSP-29, 84-91.
*	            4. Stearns, S.D. 1975 Digital Signal Analysis,
*	               Hayden Book Co., Rochelle Park, New Jersey.
*	            5. Sirivat, A. 1985 Statistical Package for the
*		       Analysis of Data (SPAD), private communication. 

*	   IFMAX  : maximum # of files that may be processed.
*	   JFMAX  : maximum # of files read from the command line.
*	   NMAX   : maximum # points per record.
*	   NSMAX  : maximum # of filter sections.
*          NUMCON : # of coefficients; order of polynomial + 1
*	   EPS    : smallest allowable number in output data.

*	       Phase Adjust Code parameters:

*	   NPTS  : # pts before current point to fit
*	   NSD   : # points in moving average of s.d. (odd # only)
*	   SDLIM : quit when s.d. theta .GT. SDLIM

*	This code is compiled using Microsoft Fortran V 5.1.  For full
*	optimization, at least 603K of RAM must be available when
*	compiled.  For additional information refer to documentation.

*	This program reads in a data file consisting of integers
*	created by an A/D sampling routine and transforms the data
*	to voltages or to some other variable.  This other variable
*	is computed using up to a fourth order polynomial; the
*	coefficients must be provided in either a sequential formatted
*	(Ascii) or a sequential binary file.  Alternatively, the routine
*	can read REAL*4 data which has already been transformed by some
*	other program.  The data can then be tapered in order to
*	reduce side-lobe leakage.  The tapering may employ a 50% overlap
*	and the user has a choice of six different window functions.
*	Then, the spectra are computed; the user may choose from :
*	amplitude, power or cross spectra.  The user may process
*	up to IFMAX files; the parameters chosen for the first file
*	are then applied to each of the IFMAX files.  Each record may
*	contain no more than NMAX points.  Note that in the case of two
*	channels of data per record, each channel may consist of up to
*	NMAX/2 points.  The speed of the routine is limited by the large
*	number of I/O operations that must take place.  Reduce # of points
*	per record and/or the # of records per file to increase the speed
*	of the routine.  Alternatively, specify in the configuration file
*	a RAM disk for storage of temporary files -- this also increases
*	speed.  The output is stored to a .PRN file with the same name
*	as the input data file in order that it may be imported into
*	LOTUS and plotted.  In order to avoid taking the logarithm of
*	zero, any number in the output data less than EPS is set equal
*	to EPS.  Then, the logarithm of EPS will be a finite number.
*	The order of the recursive digital filters is twice the number
*	of cascaded sections desired; the number of filter sections
*	must be less than or equal to NSMAX.  The routine will read
*	up to JFMAX filenames from the command line.

	INCLUDE		'SPECTRUM.FI'	! Interface statements
	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

*	Allocatable arrays may not be placed within COMMON blocks
*	and may not be dummy arguments within functions or subroutines.

	INTEGER*2 	IFIL[ALLOCATABLE](:),GAIN(0:7)
	INTEGER*2 	NDATA[ALLOCATABLE,HUGE](:)
	REAL*4		WIN[ALLOCATABLE,HUGE](:)
	REAL*4		FR[ALLOCATABLE,HUGE](:),FI[ALLOCATABLE,HUGE](:)
	REAL*4		FR2[ALLOCATABLE,HUGE](:),FI2[ALLOCATABLE,HUGE](:)
	REAL*4 		CONST [ALLOCATABLE](:,:),B (NUMCON)
	REAL*4		CONST2[ALLOCATABLE](:,:),B2(NUMCON)
	DATA 		TNAM,WNAM,PNAM,HNAM,SNAM,DNAM
     +                  /'Tukey','Welch','Parzen','Hanning',
     +			'Min Side Lobe','Max Decay'/

*	Initialize.

	PI=2.0*ASIN(1.0)
	PI2=2.0*PI
	COH=.FALSE.
	NOCOH=.TRUE.
	ROOT1=.TRUE.
	ROOT2=.TRUE.
	SUFX=.FALSE.

*	Code to read data from configuration file.

	INQUIRE (FILE='SPECTRUM.CFG',EXIST=EXISTS)
	IF (.NOT. EXISTS) GO TO 5

*	File exists. Get the info.

	CALL RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)

*	Code to perform command line parsing.

5	NUMARGS=NARGS()
	IF (NUMARGS .EQ. 1) THEN
	  NUMARGS=NUMARGS-1
	  GO TO 10
	ELSE
	  NUMARGS=NUMARGS-1
	  IF (NUMARGS .GT. JFMAX) THEN
	    WRITE (*,'(/1X,A/1X,I1,A/1X,A)')
     +      'The routine is currently set to recognize up to ',
     +      JFMAX,' filenames on the command line and will',
     +      'proceed with these filenames only.'
	    NUMARGS=JFMAX
	  ENDIF
	  DO I=1,NUMARGS
	    CALL GETARG (I,BUFFER,STATUS)
	    DO J=4,1,-1
	      FIRST=BUFFER(J:J)
	      IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
	        IHOLD=ICHAR(FIRST)-32
	        FIRST=CHAR(IHOLD)
	        BUFFER(J:J)=FIRST
	      ENDIF
	    ENDDO
	    IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
	      WRITE (*,'(/1X,A,I1,A/1X,A,A,A/)') 
     +       'Filename ',I,' on the command line began with',
     +       'the nonalphabetic character ',FIRST,'.'
	      STOP 'Re-enter the command line correctly'
	    ENDIF
	    IF (.NOT. ROOT1) THEN
	      ISL=INDEX (PTH,'\',.TRUE.)
	      FLNM(I)=PTH(:ISL) // BUFFER // '.DAT'
	    ELSE
	      FLNM(I)=BUFFER // '.DAT'
	    ENDIF
	    INQUIRE (FILE=FLNM(I),EXIST=EXISTS)
	    IF (.NOT. EXISTS) THEN
	      IPD=INDEX (FLNM(I),'.',.TRUE.)
	      WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
     +                                FLNM(I)(:IPD+3)
	      STOP 'Re-enter the command line correctly'
	    ENDIF
	  ENDDO
	  CONSEC=.TRUE.
	  ARB=.FALSE.
	  IFSTRT=1
	  IFSTP=NUMARGS
	  IFTOT=NUMARGS
	  GO TO 40
	ENDIF

*	This controls the values that are fed to IF.

10	WRITE (*,'(/1X,A,A\)') 'Process files (C)onsecutively',
     +                         ' or in an (A)rbitrary order : '
	READ (*,'(A)') INP
	WRITE (*,'( )')
	IF (INP .EQ. 'c') INP = 'C'
	IF (INP .EQ. 'a') INP = 'A'
	CONSEC=(INP .EQ. 'C')
	ARB=(INP .EQ. 'A')
	IF (.NOT. CONSEC .AND. .NOT. ARB) GO TO 10

*	These statements will allow IF to be set equal to values
*	stored in the IFIL array in the loop below.

	IF (ARB) THEN
20	  WRITE (*,'(1X,A\)') 'Enter number of files to process : '
	  READ (*,*) IFTOT
	  WRITE (*,'( )')
	  IF (IFTOT .GT. IFMAX) THEN
	    WRITE (*,'(1X,''# files must be .LE. '',I3,''.''//)') IFMAX
	    GO TO 20
	  ENDIF
	  ALLOCATE (IFIL(IFTOT), STAT=IERR)  ! Allocate space for IFIL
	  IF (IERR .NE. 0)
     +       STOP 'Problem allocating storage for IFIL.  Aborting ...'
	  IFILMIN=IFMAX
	  IFILMAX=1
	  DO I=1,IFTOT
	    WRITE (*,'(1X,''Enter last digits of file '',I2,'' : ''\)') I
	    READ (*,*) IFIL(I)
	    IFILMIN=MIN0(IFIL(I),IFILMIN)  ! Minimum of IFIL array
	    IFILMAX=MAX0(IFIL(I),IFILMAX)  ! Maximum of IFIL array
	  ENDDO
	  WRITE (*,'( )')
	  IFSTRT=1
	  IFSTP=IFTOT
	  WRITE (*,'(1X,A)') 'Files are : '
	  WRITE (*,'(13X,10(I2,2X))') (IFIL(I),I=1,IFTOT)
	  WRITE (*,'( )')
	ELSE
25	  WRITE (*,'(1X,A\)') 'Enter starting & ending data file # : '
	  READ (*,*) IFSTRT,IFSTP
	  IF (IFSTP .LT. IFSTRT) THEN
	    WRITE (*,'(1X,A)')
     +             'Ending data file # >= Starting data file #.'
	    GO TO 25
	  ENDIF
	  WRITE (*,'( )')
	  IFTOT=IFSTP-IFSTRT+1
	ENDIF

*	Get the first letter.

30	WRITE (*,'(1X,A\)') 'Enter first letter of data file : '
	READ (*,'(A)') FIRST
	IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
	  IHOLD=ICHAR(FIRST)-32
	  FIRST=CHAR(IHOLD)
	ENDIF
	IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
	  WRITE (*,'(1X,A/)') 'Enter an alphabetic character (A-Z).'
	  GO TO 30
	ENDIF

*	Assign coefficient data file names.

40	IF (.NOT. ROOT1) THEN
	  ISL=INDEX (PTH,'\',.TRUE.)
	  XNAM =PTH(:ISL) // FIRST // 'CON.DAT'
	  YNAM =PTH(:ISL) // FIRST // 'CON.PRN'
	  XNAM2=PTH(:ISL) // FIRST // 'CON2.DAT'
	  YNAM2=PTH(:ISL) // FIRST // 'CON2.PRN'
	ELSE
	  XNAM=FIRST//'CON.DAT'
	  YNAM=FIRST//'CON.PRN'
	  XNAM2=FIRST//'CON2.DAT'
	  YNAM2=FIRST//'CON2.PRN'
	ENDIF

*	Transform to voltages, other quantity or read Real data ?

50	WRITE (*,'(/1X,A/1X,A\)')
     +    'Transform into (V)oltages or into (O)ther quantity',
     +    'or read (R)eal data which is already transformed : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'v') INP = 'V'
	IF (INP .EQ. 'o') INP = 'O'
	IF (INP .EQ. 'r') INP = 'R'
	VOLT=(INP .EQ. 'V')
	OTHER=(INP .EQ. 'O')
	REALDAT=(INP .EQ. 'R')
	IF (.NOT. VOLT .AND. .NOT. OTHER .AND.
     +      .NOT. REALDAT) GO TO 50
	IF (REALDAT) GO TO 70

*	Get the voltage transformation factor.

60	WRITE (*,'(/1X,A\)')
     +   'Alternate voltage transformation factor (Y/N) : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'y') INP = 'Y'
	IF (INP .EQ. 'n') INP = 'N'
	ALTVOLT=(INP .EQ. 'Y')
	USEGAIN=(INP .EQ. 'N')
	IF (.NOT. ALTVOLT .AND. .NOT. USEGAIN) GO TO 60

	IF (ALTVOLT) THEN
	  WRITE (*,'(/1X,A\)') 'Enter voltage transformation factor : '
	  READ (*,*) VOFST1
	ENDIF

*	Remove the mean value from the file ?

70	WRITE (*,'(/1X,A,A\)') 'Remove (M)ean value,',
     +         ' (D)etrend or (N)either : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'm') INP = 'M'
	IF (INP .EQ. 'd') INP = 'D'
	IF (INP .EQ. 'n') INP = 'N'
	RMEAN=(INP .EQ. 'M')
	DTREND=(INP .EQ. 'D')
	NOOP=(INP .EQ. 'N')
	IF (.NOT. RMEAN .AND. .NOT. DTREND .AND. .NOT. NOOP) GO TO 70
	WRITE (*,'( )')
	IF (VOLT .OR. REALDAT) GO TO 90

*	Determine the type of data file for coefficients.

80	WRITE (*,'(1X,A,A\)') 'Read from (A)scii or',
     +                        ' (B)inary data file : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'a') INP = 'A'
	IF (INP .EQ. 'b') INP = 'B'
	ASCII=(INP .EQ. 'A')
	BINRY=(INP .EQ. 'B')
	IF (.NOT. ASCII .AND. .NOT. BINRY) GO TO 80
	WRITE (*,'( )')

*	Allocate space for CONST array.

	IF (CONSEC) ALLOCATE (CONST(IFSTRT:IFSTP,NUMCON), STAT=IERR)
	IF (ARB) ALLOCATE (CONST(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
	IF (IERR .NE. 0)
     +     STOP 'Problem allocating storage for CONST.  Aborting ...'

*	Open the formatted file containing the transformation
*	constants and fill the array.

	IF (ASCII) THEN

	  INQUIRE (FILE=YNAM,EXIST=EXISTS)
	  IF (.NOT. EXISTS) THEN
	    IPD=INDEX (YNAM,'.',.TRUE.)
	    WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
     +                              YNAM(:IPD+3)
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

	  OPEN (NUMI,FILE=YNAM,STATUS='OLD')
	    READ (NUMI,*) (B(I),I=1,5)
	  CLOSE (NUMI,ERR=400)

	  DO I=IFSTRT,IFSTP
	    IF (CONSEC) IF=I
	    IF (ARB) IF=IFIL(I)
	    DO J=1,5
	      CONST(IF,J)=B(J)
	    ENDDO
	  ENDDO

	ELSE IF (BINRY) THEN

*	  Open the binary file containing the transformation
*	  constants and fill the array.

	  INQUIRE (FILE=XNAM,EXIST=EXISTS)
	  IF (.NOT. EXISTS) THEN
	    IPD=INDEX (XNAM,'.',.TRUE.)
	    WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
     +                              XNAM(:IPD+3)
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

	  OPEN (NUMI,FILE=XNAM,STATUS='OLD',ACCESS='SEQUENTIAL',
     +          FORM='BINARY')
	    READ (NUMI) NUMSETS,NSTART
	    IF (
     +      (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
     +      .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
     +                         NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
	      WRITE (*,'(/1X,A/1X,A/1X,A/)')
     +	      'Starting or ending data file #''s stored in coefficient',
     +        'data file for channel 1 do not match values entered',
     +        'into program.  This must be corrected before proceeding.'
	      CLOSE (NUMI)
	      STOP '                   Program terminated unsuccessfully.'
	    ENDIF	      
	    READ (NUMI)
     +         ((CONST(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
	  CLOSE (NUMI,ERR=400)

	ENDIF

*	Filter the data ?

90	WRITE (*,'(1X,A/20X,A\)')
     +    'Filter the data using (L)owpass, (H)ighpass,',
     +                      '(B)andpass, Band(S)top or (N)one : '
	READ (*,'(A)') FTYPE
	IF (FTYPE .EQ. 'l') FTYPE = 'L'
	IF (FTYPE .EQ. 'h') FTYPE = 'H'
	IF (FTYPE .EQ. 'b') FTYPE = 'B'
	IF (FTYPE .EQ. 's') FTYPE = 'S'
	IF (FTYPE .EQ. 'n') FTYPE = 'N'
	LPASS=(FTYPE .EQ. 'L')
	HPASS=(FTYPE .EQ. 'H')
	BPASS=(FTYPE .EQ. 'B')
	BSTOP=(FTYPE .EQ. 'S')
	NOFIL=(FTYPE .EQ. 'N')
	IF (.NOT. LPASS .AND. .NOT. HPASS .AND. .NOT. BPASS .AND.
     +      .NOT. BSTOP .AND. .NOT. NOFIL) GO TO 90
	IF (NOFIL) GO TO 110

*	Get the cutoff frequencies.

	FC=0.0
	F1=0.0
	F2=0.0

	IF (LPASS .OR. HPASS) THEN
	  WRITE (*,'(/1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
	  READ (*,*) FC
	ELSE IF (BPASS .OR. BSTOP) THEN
	  WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
	  READ (*,*) F1
	  WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
	  READ (*,*) F2
	ENDIF

100	WRITE (*,'(/1X,A/5X,A/1X,A\)')
     +             'The order of the filter will be twice',
     +             'the number of cascaded sections.',
     +             'Enter # of filter sections to cascade : '
	READ (*,*) NS
	IF (NS .GT. NSMAX) THEN
	  WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
	  GO TO 100
	ENDIF

*	Two channels of data ?

110	WRITE (*,'(/1X,A\)') 'How many channels of data (1 or 2) : '
	READ (*,*) ICHANS
	WRITE (*,'( )')
	ONECHAN=(ICHANS .EQ. 1)
	TWOCHAN=(ICHANS .EQ. 2)
	IF (.NOT. ONECHAN .AND. .NOT. TWOCHAN) GO TO 110
	IF (ONECHAN) GO TO 170
	IF (TWOCHAN .AND. REALDAT) GO TO 130

*	Transform second channel to voltages or to some other quantity ?

120	WRITE (*,'(1X,A,A\)') 'Transform 2nd channel into',
     +                ' (V)oltages or into (O)ther quantity : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'v') INP = 'V'
	IF (INP .EQ. 'o') INP = 'O'
	VOLT2=(INP .EQ. 'V')
	OTHER2=(INP .EQ. 'O')
	IF (.NOT. VOLT2 .AND. .NOT. OTHER2) GO TO 120
	WRITE (*,'( )')

*	Get the voltage transformation factor for channel 2.

125	WRITE (*,'(1X,A\)')
     +   'Alternate voltage transformation factor (Y/N) : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'y') INP = 'Y'
	IF (INP .EQ. 'n') INP = 'N'
	ALTVOLT2=(INP .EQ. 'Y')
	USEGAIN2=(INP .EQ. 'N')
	IF (.NOT. ALTVOLT2 .AND. .NOT. USEGAIN2) GO TO 125
	WRITE (*,'( )')

	IF (ALTVOLT2) THEN
	  WRITE (*,'(1X,A\)') 'Enter voltage transformation factor : '
	  READ (*,*) VOFST2
	  WRITE (*,'( )')
	ENDIF

*	Remove the mean value from the second channel of data ?

130	WRITE (*,'(1X,A,A\)') '2nd channel: remove (M)ean value,',
     +         ' (D)etrend or (N)either : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'm') INP = 'M'
	IF (INP .EQ. 'd') INP = 'D'
	IF (INP .EQ. 'n') INP = 'N'
	RMEAN2=(INP .EQ. 'M')
	DTREND2=(INP .EQ. 'D')
	NOOP2=(INP .EQ. 'N')
	IF (.NOT. RMEAN2 .AND. .NOT. DTREND2 .AND.
     +      .NOT. NOOP2) GO TO 130
	WRITE (*,'( )')
	IF (VOLT2 .OR. REALDAT) GO TO 150

*	Determine type of data file for coefficients for second channel.

140	WRITE (*,'(1X,A,A\)') '2nd channel: read from (A)scii or',
     +                        ' (B)inary data file : '
	READ (*,'(A)') INP
	IF (INP .EQ. 'a') INP = 'A'
	IF (INP .EQ. 'b') INP = 'B'
	ASCII2=(INP .EQ. 'A')
	BINRY2=(INP .EQ. 'B')
	IF (.NOT. ASCII2 .AND. .NOT. BINRY2) GO TO 140
	WRITE (*,'( )')

*	Allocate space for CONST2 array.

	IF (CONSEC) ALLOCATE (CONST2(IFSTRT:IFSTP,NUMCON), STAT=IERR)
	IF (ARB) ALLOCATE (CONST2(IFILMIN:IFILMAX,NUMCON), STAT=IERR)
	IF (IERR .NE. 0)
     +     STOP 'Problem allocating storage for CONST2.  Aborting ...'

*	Open the formatted file containing the transformation
*	constants for the second channel and fill the array.

	IF (ASCII2) THEN

	  INQUIRE (FILE=YNAM2,EXIST=EXISTS)
	  IF (.NOT. EXISTS) THEN
	    IPD=INDEX (YNAM2,'.',.TRUE.)
	    WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
     +                              YNAM2(:IPD+3)
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

	  OPEN (NUMI,FILE=YNAM2,STATUS='OLD')
	    READ (NUMI,*) (B2(I),I=1,5)
	  CLOSE (NUMI,ERR=400)

	  DO I=IFSTRT,IFSTP
	    IF (CONSEC) IF=I
	    IF (ARB) IF=IFIL(I)
	    DO J=1,5
	      CONST2(IF,J)=B2(J)
	    ENDDO
	  ENDDO

	ELSE IF (BINRY2) THEN

*	  Open the binary file containing the transformation
*	  constants for the second channel and fill the array.

	  INQUIRE (FILE=XNAM2,EXIST=EXISTS)
	  IF (.NOT. EXISTS) THEN
	    IPD=INDEX (XNAM2,'.',.TRUE.)
	    WRITE (*,'(/1X,A,A/)') 'Cannot find file : ',
     +                              XNAM2(:IPD+3)
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

	  OPEN (NUMI,FILE=XNAM2,STATUS='OLD',ACCESS='SEQUENTIAL',
     +          FORM='BINARY')
	    READ (NUMI) NUMSETS,NSTART
	    IF (
     +      (CONSEC .AND. (NSTART .NE. IFSTRT .OR. NUMSETS .NE. IFTOT))
     +      .OR. (ARB .AND. (NSTART .NE. IFILMIN .OR.
     +                         NSTART+NUMSETS-1 .NE. IFILMAX)) ) THEN
	      WRITE (*,'(/1X,A/1X,A/1X,A/)')
     +	      'Starting or ending data file #''s stored in coefficient',
     +        'data file for channel 2 do not match values entered',
     +        'into program.  This must be corrected before proceeding.'
	      CLOSE (NUMI)
	      STOP '                   Program terminated unsuccessfully.'
	    ENDIF	      
	    READ (NUMI)
     +        ((CONST2(NSTART+II-1,JJ),JJ=1,NUMCON),II=1,NUMSETS)
	  CLOSE (NUMI,ERR=400)

	ENDIF

*	Filter the data for the second channel ?

150	WRITE (*,'(1X,A/20X,A\)')
     +    'Filter the data using (L)owpass, (H)ighpass,',
     +                      '(B)andpass, Band(S)top or (N)one : '
	READ (*,'(A)') FTYP2
	IF (FTYP2 .EQ. 'l') FTYP2 = 'L'
	IF (FTYP2 .EQ. 'h') FTYP2 = 'H'
	IF (FTYP2 .EQ. 'b') FTYP2 = 'B'
	IF (FTYP2 .EQ. 's') FTYP2 = 'S'
	IF (FTYP2 .EQ. 'n') FTYP2 = 'N'
	LPASS2=(FTYP2 .EQ. 'L')
	HPASS2=(FTYP2 .EQ. 'H')
	BPASS2=(FTYP2 .EQ. 'B')
	BSTOP2=(FTYP2 .EQ. 'S')
	NOFIL2=(FTYP2 .EQ. 'N')
	IF (.NOT. LPASS2 .AND. .NOT. HPASS2 .AND. .NOT. BPASS2 .AND.
     +      .NOT. BSTOP2 .AND. .NOT. NOFIL2) GO TO 150
	WRITE (*,'( )')
	IF (NOFIL2) GO TO 170

*	Get the cutoff frequencies.

	FC2=0.0
	F12=0.0
	F22=0.0

	IF (LPASS2 .OR. HPASS2) THEN
	  WRITE (*,'(1X,A\)') 'Enter the (-3 dB) frequency cutoff : '
	  READ (*,*) FC2
	ELSE IF (BPASS2 .OR. BSTOP2) THEN
	  WRITE (*,'(1X,A\)') 'Enter the low frequency cutoff : '
	  READ (*,*) F12
	  WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
	  READ (*,*) F22
	ENDIF

160	WRITE (*,'(/1X,A\)') 'Enter # of filter sections to cascade : '
	READ (*,*) NS2
	WRITE (*,'( )')
	IF (NS2 .GT. NSMAX) THEN
	  WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
	  GO TO 160
	ENDIF

*	Determine the kind of spectrum.

170	WRITE (*,'(1X,A\)')
     +            '(A)mplitude, (P)ower or (C)ross spectrum : '
	READ (*,'(A)') INP
	WRITE (*,'( )')
	IF (INP .EQ. 'a') INP='A'
	IF (INP .EQ. 'p') INP='P'
	IF (INP .EQ. 'c') INP='C'
	AMP=(INP .EQ. 'A')
	POWER=(INP .EQ. 'P')
	CROSS=(INP .EQ. 'C')
	IF (.NOT. AMP .AND. .NOT. POWER .AND. .NOT. CROSS) GO TO 170
	IF (CROSS .AND. ONECHAN) THEN
	  WRITE (*,'(1X,A,A/)') 'Cross Spectrum requires',
     +                          ' two channels of data.'
	  GO TO 170
	ENDIF

*	Calculate the coherence function ?

180	IF (CROSS) THEN
	  WRITE (*,'(1X,A/1X,A//1X,A\)')
     +    'Coherence function requires calculation',
     +    'of both cross and autospectral density functions.',
     +    'Calculate the coherence function gamma (f) (Y/N) : '
	  READ (*,'(A)') INP
	  IF (INP .EQ. 'y') INP = 'Y'
	  IF (INP .EQ. 'n') INP = 'N'
	  COH=(INP .EQ. 'Y')
	  NOCOH=(INP .EQ. 'N')
	  IF (.NOT. COH .AND. .NOT. NOCOH) GO TO 180
	  WRITE (*,'( )')
	ENDIF

*	Taper the data or not ?

190	WRITE (*,'(1X,A\)') 'Taper the data (Y/N) : '
	READ (*,'(A)') INP
	WRITE (*,'( )')
	IF (INP .EQ. 'y') INP='Y'
	IF (INP .EQ. 'n') INP='N'
	TAPER=(INP .EQ. 'Y')
	NOTAPER=(INP .EQ. 'N')
	IF (.NOT. TAPER .AND. .NOT. NOTAPER) GO TO 190

*	Choose the window function.

	IF (TAPER) THEN
200	  WRITE (*,'(1X,A/1X,A\)')
     +	       '(H)anning, (P)arzen, (T)ukey, (W)elch,',
     +         'Min (S)ide Lobe or Max (D)ecay window : '
	  READ (*,'(A)') INP
	  WRITE (*,'( )')
	  IF (INP .EQ. 'h') INP='H'
	  IF (INP .EQ. 'p') INP='P'
	  IF (INP .EQ. 't') INP='T'
	  IF (INP .EQ. 'w') INP='W'
	  IF (INP .EQ. 's') INP='S'
	  IF (INP .EQ. 'd') INP='D'
	  HANN =(INP .EQ. 'H')
	  PARZ =(INP .EQ. 'P')
	  TUKEY=(INP .EQ. 'T')
	  WELCH=(INP .EQ. 'W')
	  MINSL=(INP .EQ. 'S')
	  MXDEC=(INP .EQ. 'D')
	  IF (.NOT. HANN  .AND. .NOT. PARZ  .AND.
     +        .NOT. WELCH .AND. .NOT. TUKEY .AND.
     +        .NOT. MINSL .AND. .NOT. MXDEC) GO TO 200
	ENDIF

*	Overlap the data or not ?

	IF (TAPER) THEN
205	  WRITE (*,'(1X,A\)') 'Use overlap (Y/N) : '
	  READ (*,'(A)') INP
	  WRITE (*,'( )')
	  IF (INP .EQ. 'y') INP='Y'
	  IF (INP .EQ. 'n') INP='N'
	  OVRLP  =(INP .EQ. 'Y')
	  NOOVRLP=(INP .EQ. 'N')
	  IF (.NOT. OVRLP .AND. .NOT. NOOVRLP) GO TO 205
	ENDIF

*	Get the time just before starting.

	CALL GETTIM (IHR,IMIN,ISEC,I100TH)

*	Do the calculations for each of the desired files.

	DO KOUNT=IFSTRT,IFSTP
	  IF (CONSEC) IF=KOUNT
	  IF (ARB) IF=IFIL(KOUNT)

*	  Assign filenames.

	  IF (NUMARGS .GT. 0) THEN
	    FIN=FLNM(IF)
	  ELSE
	    IF (.NOT. ROOT1) THEN
	      ISL=INDEX (PTH,'\',.TRUE.)
	      WRITE (FIN,'(A,A1,I3.3,A4)') PTH(:ISL),FIRST,IF,'.DAT'
	    ELSE
	      WRITE (FIN,'(A1,I3.3,A4)') FIRST,IF,'.DAT'
	    ENDIF
	    INQUIRE (FILE=FIN,EXIST=EXISTS)
	    IF (.NOT. EXISTS) THEN
	      IPD=INDEX (FIN,'.',.TRUE.)
	      WRITE (*,'(25X,A,A,A/)') 'Cannot find file : ',
     +                                  FIN(:IPD+3),'.'
	      CYCLE
	    ENDIF
	  ENDIF

	  IPD=INDEX (FIN,'.',.TRUE.)
	  IF (ROOT2) THEN
	    TDOUT =FIN(IPD-4:IPD-1) // '.TRF'
	    TDOUT2=FIN(IPD-4:IPD-1) // '.TR2'
	    TOUT  =FIN(IPD-4:IPD-1) // '.TAP'
	    TOUT2 =FIN(IPD-4:IPD-1) // '.TP2'
	  ELSE IF (.NOT. ROOT2) THEN
	    ISL=INDEX (TPTH,'\',.TRUE.)
	    TDOUT =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TRF'
	    TDOUT2=TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TR2'
	    TOUT  =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TAP'
	    TOUT2 =TPTH(:ISL) // FIN(IPD-4:IPD-1) // '.TP2'
	  ENDIF

	  IF (.NOT. SUFX) THEN
	    IF (POWER) POUT=FIN(:IPD-1) // 'P.PRN'
	    IF (AMP)   POUT=FIN(:IPD-1) // 'A.PRN'
	    IF (CROSS) POUT=FIN(:IPD-1) // 'C.PRN'
	    EOUT  =FIN(:IPD-1) // 'E.PRN'
	  ELSE
	    IF (POWER) POUT=FIN(:IPD-1) // 'P' // SUFFIX
	    IF (AMP)   POUT=FIN(:IPD-1) // 'A' // SUFFIX
	    IF (CROSS) POUT=FIN(:IPD-1) // 'C' // SUFFIX
	    EOUT  =FIN(:IPD-1) // 'E' // SUFFIX
	  ENDIF

*	  Put message on screen.

	  IPD=INDEX (FIN,'.',.TRUE.)
	  IF (KOUNT .EQ. IFSTRT) THEN
	    WRITE (*,'(/////////////////////////13X,
     +    ''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'')')
	    WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
	  ELSE
	    WRITE (*,'(/25X,''Processing '',A,'' now.''/)') FIN(:IPD+3)
	  ENDIF

*	  Open the data file and input the data.

	  OPEN (NUMI,FILE=FIN,STATUS='OLD',ACCESS='SEQUENTIAL',
     +          FORM='BINARY')

	  READ (NUMI,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
	  READ (NUMI,ERR=400) (GAIN(I),I=0,7)

*	  Calculate # points per record and adjust record size for
*	  the file containing transformed data.

	  IF (REALDAT) THEN
	    N=IRSIZE/4
	    IRSIZ2=IRSIZE
	    IRSIZE=IRSIZE/2
	  ELSE
	    N=IRSIZE/2
	    IRSIZ2=2*IRSIZE
	  ENDIF
	  DELT=FLOAT(IDELTMS)/1000000.0
	  FN=FLOAT(N)

*	  N less than or equal to NMAX error checking.

	  IF (ONECHAN) NTST=NMAX
	  IF (TWOCHAN) NTST=NMAX/2
	  IF (N .GT. NMAX) THEN
	    WRITE (*,'(15X,A,I5,A/15X,A,I5/)')
     +         'This data file contains ',N,' data points.',
     +         '# data points per record per channel must be < ',NTST
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

*	  Test to make sure that N is a power of 2.

	  ITST=NINT(ALOG10(FN)/ALOG10(2.0))
	  ITST2=INT(2**ITST)-N

	  IF (ITST2 .NE. 0) THEN
	    WRITE (*,'(20X,A,I5,A/20X,A/)') 'This data file contains ',
     +          N,' data points.','# data points must be a power of 2.'
	    STOP '                   Program terminated unsuccessfully.'
	  ENDIF

*	  Allocate space for various arrays.

	  IF (REALDAT) THEN
	    ALLOCATE (WIN(N), FR (N), FI (N),
     +                                  FR2(N), FI2(N), STAT=IERR)
	    IF (IERR .NE. 0)
     +              STOP 'Not enough storage for data.  Aborting ...'
	  ELSE
	    ALLOCATE (NDATA(N), WIN(N), FR (N), FI (N),
     +                                  FR2(N), FI2(N), STAT=IERR)
	    IF (IERR .NE. 0)
     +              STOP 'Not enough storage for data.  Aborting ...'
	  ENDIF

*	  Test to be sure data file really has one channel.

	  IF (ONECHAN .AND. ICHANS .NE. 1) THEN
	    WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
     +      ' file indicates two channels.',
     +      'You have specified one channel.  Please correct problem.'
	    STOP '                       Program terminated unsuccessfully
     +.'
	  ENDIF

*	  Test to be sure data file really has two channels.

	  IF (TWOCHAN .AND. ICHANS .NE. 2) THEN
	    WRITE (*,'(13X,A,A/13X,A)') 'The header of this data',
     +      ' file indicates only one channel.',
     +      'You have specified two channels.  Please correct problem.'
	    STOP '                       Program terminated unsuccessfully
     +.'
	  ENDIF

*	  Estimate amount of temporary space needed.

	  IF (REALDAT) THEN
	    NBYTE=IRSIZ2*NUMREC
	  ELSE
	    NBYTE=IRSIZE*NUMREC
	  ENDIF

	  IF (.NOT. REALDAT) TRFFIL=2*NBYTE
	  IF (REALDAT) TRFFIL=NBYTE
	  IF (NOTAPER) TAPFIL=0
	  IF (TAPER .AND. OVRLP) TAPFIL=2*TRFFIL
	  IF (TAPER .AND. NOOVRLP) TAPFIL=TRFFIL

	  TEMP_SPACE=TRFFIL+TAPFIL
	  TSPCMB=FLOAT(TEMP_SPACE)/1048576.0

	  WRITE (*,'(13X,A,I8,A,F6.2,A/)')
     +      'Temporary space required = ',TEMP_SPACE,
     +      ' bytes (',TSPCMB,' Mb).'

*	  Perform a write test to determine if space is available
*	  to hold the temporary data files.

	  OPEN (NUMT, FILE=TDOUT, FORM='BINARY',RECL=1,
     +	    ACCESS='DIRECT', STATUS='UNKNOWN')
	    IREC=TEMP_SPACE+(HEADER_SIZE*2)-1
	    WRITE (NUMT, REC=IREC, ERR=390) #FF
	  CLOSE (NUMT, STATUS='DELETE')

*	  Determine the voltage transformation factor.

	  IF (USEGAIN)  VOFST1=(20.0/4096.0)/(2**GAIN(0))
	  IF (USEGAIN2) VOFST2=(20.0/4096.0)/(2**GAIN(1))

*	  Open the output file to hold the transformed data.

	  OPEN (NUMO,FILE=TDOUT,STATUS='UNKNOWN',
     +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)
	  IF (TWOCHAN) OPEN (NUMO2,FILE=TDOUT2,STATUS='UNKNOWN',
     +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)

*	  Write file header.

	  IF (ONECHAN) THEN
	    WRITE (NUMO,ERR=400) ICHANS,IRSIZ2,NUMREC,IDELTMS
	  ELSE IF (TWOCHAN) THEN
	    WRITE (NUMO,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
	    WRITE (NUMO2,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
	  ENDIF

*	  Transformation procedure for one channel of data.

	  IF (ONECHAN) THEN

	    IF (OTHER) THEN
	      DO J=1,NUMCON
	        B(J)=CONST(IF,J)
	      ENDDO
	    ELSE
	      B=0.0
	    ENDIF

	    IF (TRANS_ONE_CHAN (NDATA,FR,FI,B) .NE. 0) STOP
     +        '                      Program terminated unsuccessfully.'

*	    Close the input data file and proceed
*	    to the tapering operation.

	    CLOSE (NUMI,STATUS='KEEP',ERR=400)

*	  Transformation procedure for two channels of data.

	  ELSE IF (TWOCHAN) THEN

	    IF (OTHER) THEN
	      DO J=1,NUMCON
	        B(J)=CONST(IF,J)
	      ENDDO
	    ELSE
	      B=0.0
	    ENDIF

	    IF (OTHER2) THEN
	      DO J=1,NUMCON
	        B2(J)=CONST2(IF,J)
	      ENDDO
	    ELSE
	      B2=0.0
	    ENDIF

	    IF (TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)
     +        .NE. 0) STOP
     +        '                      Program terminated unsuccessfully.'

*	    Close the input data file and proceed
*	    to the tapering operation.

	    CLOSE (NUMI,STATUS='KEEP',ERR=400)

	  ENDIF

*	  Messages concerning mean removal, detrending or filtering.

	  MSG=.FALSE.
	  IF (TWOCHAN) THEN
	    IF (RMEAN) WRITE (*,'(25X,A,A)') 'Mean value removed',
     +                       ' from channel 1.'
	    IF (DTREND) WRITE (*,'(25X,A)') 'Channel 1 detrended.'
	    IF (LPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
     +                                    ' was lowpass filtered.'
	    IF (HPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
     +                                    ' was highpass filtered.'
	    IF (BPASS) WRITE (*,'(25X,A,A)') 'Channel 1',
     +                                    ' was bandpass filtered.'
	    IF (BSTOP) WRITE (*,'(25X,A,A)') 'Channel 1',
     +                                    ' was bandstop filtered.'
	    IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
     +                                'Order of filter = ',2*NS,'.'
	    IF (RMEAN2) WRITE (*,'(25X,A,A)') 'Mean value removed',
     +                       ' from channel 2.'
	    IF (DTREND2) WRITE (*,'(25X,A)') 'Channel 2 detrended.'
	    IF (LPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
     +                                    ' was lowpass filtered.'
	    IF (HPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
     +                                    ' was highpass filtered.'
	    IF (BPASS2) WRITE (*,'(25X,A,A)') 'Channel 2',
     +                                    ' was bandpass filtered.'
	    IF (BSTOP2) WRITE (*,'(25X,A,A)') 'Channel 2',
     +                                    ' was bandstop filtered.'
	    IF (.NOT. NOFIL2) WRITE (*,'(25X,A,I2,A)')
     +                               'Order of filter = ',2*NS2,'.'
	    IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
	    IF (RMEAN2 .OR. DTREND2 .OR. .NOT. NOFIL) MSG=.TRUE.
	  ELSE
	    IF (RMEAN) WRITE (*,'(25X,A)') 'Mean value removed.'
	    IF (DTREND) WRITE (*,'(25X,A)') 'File was detrended.'
	    IF (LPASS) WRITE (*,'(25X,A)') 'File was lowpass filtered.'
	    IF (HPASS) WRITE (*,'(25X,A)') 'File was highpass filtered.'
	    IF (BPASS) WRITE (*,'(25X,A)') 'File was bandpass filtered.'
	    IF (BSTOP) WRITE (*,'(25X,A)') 'File was bandstop filtered.'
	    IF (.NOT. NOFIL) WRITE (*,'(25X,A,I2,A)')
     +                                     'Order of filter = ',2*NS,'.'
	    IF (RMEAN .OR. DTREND .OR. .NOT. NOFIL) MSG=.TRUE.
	  ENDIF

*	  Read the file header.

	  REWIND (NUMO)
	  READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
	  IF (TWOCHAN) THEN
	    REWIND (NUMO2)
	    READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
	  ENDIF

	  N=IRSIZE/4
	  N2=N/2
	  NREC=NUMREC             ! original numrec before overlap
	  NREC2=2*NUMREC
	  K=INT4(FLOAT(N)/10.0)

*	  No Tapering applied.

	  IF (NOTAPER) THEN
	    IF (MSG) THEN
	      WRITE (*,'(/25X,A/25X,I4.4,A)') 'No tapering applied.  ',
     +               NUMREC, ' records will be used.'
	    ELSE
	      WRITE (*,'(25X,A/25X,I4.4,A)') 'No tapering applied.  ',
     +               NUMREC, ' records will be used.'
	    ENDIF
	    GO TO 300
	  ENDIF

*	  Compute the scale factor.

	  IF (WIN_SCALE (WIN) .NE. 0) STOP
     +      '                      Program terminated unsuccessfully.'

*	  Open the data files to hold tapered data.

	  OPEN (NUMT,FILE=TOUT,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
     +          FORM='BINARY',ERR=400)

	  IF (TWOCHAN) OPEN (NUMT2,FILE=TOUT2,STATUS='UNKNOWN',
     +          ACCESS='SEQUENTIAL',FORM='BINARY',ERR=400)

	  IF (NOOVRLP) THEN

*	    No overlap - use NUMREC records.

	    WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NUMREC,IDELTMS
	    IF (TWOCHAN) WRITE (NUMT2,ERR=400)
     +                   ICHANS,IRSIZE,NUMREC,IDELTMS
	  ELSE

*	    50% overlap - use NREC2 records.

	    WRITE (NUMT,ERR=400) ICHANS,IRSIZE,NREC2,IDELTMS
	    IF (TWOCHAN) WRITE (NUMT2,ERR=400)
     +                   ICHANS,IRSIZE,NREC2,IDELTMS
	  ENDIF

	  IF (MSG) WRITE (*,'( )')
	  IF (HANN)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // HNAM // '.'
	  IF (WELCH)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // WNAM // '.'
	  IF (PARZ)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // PNAM // '.'
	  IF (TUKEY)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // TNAM // '.'
	  IF (MINSL)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // SNAM // '.'
	  IF (MXDEC)
     +  	WRITE (*,'(25X,A)') 'Taper applied = ' // DNAM // '.'

	  IF (OVRLP) THEN
	    WRITE (*,'(25X,A/25X,I4.4,A/)') '50% overlap used.',
     +	              NREC2,' records will be written.'
	  ELSE
	    WRITE (*,'(25X,A/25X,I4.4,A/)') 'No overlap used.',
     +	              NUMREC,' records will be written.'
	  ENDIF

*	  Do the tapering now.

	  IF (TAPER_DATA (FR,FR2,FI,FI2,WIN) .NE. 0) STOP
     +      '                      Program terminated unsuccessfully.'

*	  Throw away transformed data file and prepare tapered data file.

	  CLOSE (NUMO,STATUS='DELETE',ERR=400)
	  IF (TWOCHAN) CLOSE (NUMO2,STATUS='DELETE',ERR=400)
	  REWIND (NUMT)
	  READ (NUMT) ICHANS,IRSIZE,NUMREC,IDELTMS
	  IF (TWOCHAN) THEN
	    REWIND (NUMT2)
	    READ (NUMT2) ICHANS,IRSIZE,NUMREC,IDELTMS
	  ENDIF

*	  Begin spectrum computations.

300	  IF (NOTAPER) THEN
	    INUM=NUMO
	    INUM2=NUMO2
	  ELSE IF (TAPER) THEN
	    INUM=NUMT
	    INUM2=NUMT2
	  ENDIF

*	  Test for odd number of records.  This should only occur
*	  at this point for the one channel, no taper or one channel,
*	  taper and no overlap cases.

	  SIZTST=FLOAT(NUMREC)	  
	  REMAIN=AMOD(SIZTST,2.0)
	  EVENUM=(REMAIN .EQ. 0.0)
	  ODDNUM=(REMAIN .NE. 0.0)

*	  Initialization.

	  N=IRSIZE/4
	  N2=N/2
	  NP1=N+1
	  IOD=INT4(ALOG10(FLOAT(N))/ALOG10(2.0)+.5)
	  NRECD2=NUMREC/2
	  IF (ODDNUM) NRECD2=(NUMREC+1)/2
	  IF (OVRLP) NUMREC=NUMREC-1
	  XMEAN=0.0
	  YMEAN=0.0
	  DO I=1,N
	    FR2(I)=0.0
	    FI2(I)=0.0
	    WIN(I)=0.0
	  ENDDO

	  IF (ONECHAN) THEN

*	    Spectrum computations for one channel.

	    IF (SPEC_ONE_CHAN (FR,FI,FR2,FI2) .NE. 0) STOP
     +        '                      Program terminated unsuccessfully.'

*	    Deallocate memory in various arrays.

	    IF (REALDAT) THEN
	      DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
	      IF (IERR .NE. 0)
     +          STOP 'Problem attempting to deallocate.  Aborting ...'
	    ELSE
	      DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
	      IF (IERR .NE. 0)
     +          STOP 'Problem attempting to deallocate.  Aborting ...'
	    ENDIF

	    CYCLE

	  ELSE IF (TWOCHAN) THEN

*	    Spectrum computations for two channels.

	    IF (SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN) .NE. 0) STOP
     +        '                      Program terminated unsuccessfully.'

*	    Deallocate memory in various arrays.

	    IF (REALDAT) THEN
	      DEALLOCATE (WIN, FR, FI, FR2, FI2, STAT=IERR)
	      IF (IERR .NE. 0)
     +          STOP 'Problem attempting to deallocate.  Aborting ...'
	    ELSE
	      DEALLOCATE (NDATA, WIN, FR, FI, FR2, FI2, STAT=IERR)
	      IF (IERR .NE. 0)
     +          STOP 'Problem attempting to deallocate.  Aborting ...'
	    ENDIF

	  ENDIF

	ENDDO

*	Get the time just before ending.

	CALL GETTIM (IHR2,IMIN2,ISEC2,I100TH2)

	WRITE (*,'(/25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
     +    'Starting time : ',IHR,':',IMIN,':',ISEC,'.',I100TH

	WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
     +    'Ending time   : ',IHR2,':',IMIN2,':',ISEC2,'.',I100TH2

*	Calculate elapsed time - hundredths of a second.

	IF (I100TH2-I100TH .LT. 0) THEN
	  I100TH2=I100TH2+100
	  I100H=I100TH2-I100TH
	  ISEC2=ISEC2-1
	ELSE
	  I100H=I100TH2-I100TH
	ENDIF

*	Seconds.

	IF (ISEC2-ISEC .LT. 0) THEN
	  ISEC2=ISEC2+60
	  ISECH=ISEC2-ISEC
	  IMIN2=IMIN2-1
	ELSE
	  ISECH=ISEC2-ISEC
	ENDIF

*	Minutes.

	IF (IMIN2-IMIN .LT. 0) THEN
	  IMIN2=IMIN2+60
	  IMINH=IMIN2-IMIN
	  IHR2=IHR2-1
	ELSE
	  IMINH=IMIN2-IMIN
	ENDIF

*	Hours.

	IF (IHR2-IHR .LT. 0) THEN
	  IHR2=IHR2+24
	  IHRH=IHR2-IHR
	ELSE
	  IHRH=IHR2-IHR
	ENDIF

*	Display elapsed time.

	WRITE (*,'(25X,A,I2.2,A,I2.2,A,I2.2,A,I2.2)')
     +    'Elapsed time  : ',IHRH,':',IMINH,':',ISECH,'.',I100H

	WRITE (*,'( )')
	STOP '                        Program terminated successfully.'

*	I/O Error checking.

390	WRITE (*,'(/17X,A/)')
     +    'Not enough space available for temporary files.'

	STOP '                      Program terminated unsuccessfully.'

400	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +    'No such device or out of space !'
	STOP '                      Program terminated unsuccessfully.'
	END

	FUNCTION TRANS_ONE_CHAN (NDATA,FR,FI,B)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	INTEGER*2 	NDATA[HUGE](*)
	REAL*4		FR[HUGE](*)
	REAL*4		FI[HUGE](*)
	REAL*4 		B(*)

*	Initialization.

	TRANS_ONE_CHAN=0
	SUM1UT=0.0
	SUM2UT=0.0

*	Perform mean and RMS value calculations for each record
*	in this data file.

	DO IB=1,NUMREC

*	  Fill the array with the current record.

	  IF (REALDAT) THEN
	    READ (NUMI) (FR(I),I=1,N)
	    DO I=1,N
	      FI(I)=FR(I)*FR(I)
	    ENDDO
	  ELSE
	    READ (NUMI) (NDATA(I),I=1,N)
	  ENDIF

*	  Transform the data into voltages or velocities.
*	  Calculate squared values for RMS calculations.

	  IF (REALDAT) GO TO 10
	  DO I=1,N
	    V=FLOAT(NDATA(I))*VOFST1
	    IF (VOLT) THEN
	      FR(I)=V
	      FI(I)=FR(I)*FR(I)
	      CYCLE
	    ENDIF
	    FR(I)=B(5)
	    DO J=4,1,-1
	      FR(I)=FR(I)*V+B(J)
	    ENDDO
	    FI(I)=FR(I)*FR(I)
	  ENDDO

*	  Determine mean and RMS values for each record.

10	  SUM1U=ADDSUM(FR,N)
	  SUM1UT=SUM1UT+SUM1U
	  SUM2U=ADDSUM(FI,N)
	  SUM2UT=SUM2UT+SUM2U
	  UBAR=SUM1U/FN
	  URMS=SQRT(SUM2U/FN-UBAR**2)

*	  Display output for each record.

	  IF (IB .EQ. 1) THEN
	    WRITE (*,20) IB,UBAR,URMS
20	    FORMAT (16X,'Record ',I4.4,4X,'Ave = ',G11.5,
     +              2X,'Rms = ',G11.5)
	  ELSE
	    WRITE (*,30) IB,UBAR,URMS
30	    FORMAT ('+',15X,'Record ',I4.4,4X,'Ave = ',G11.5,
     +              2X,'Rms = ',G11.5)
	  ENDIF

*	  Save transformed values for this record.

	  IF (.NOT. NOFIL)
     +        CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
	  IF (RMEAN) THEN
	    DO I=1,N
	      FR(I)=FR(I)-UBAR
	    ENDDO
	  ENDIF
	  IF (DTREND) CALL DTRD (FR,N,DELT)
	  WRITE (NUMO,ERR=40) (FR(I), I=1,N)

	ENDDO

*	All records of this data file have been processed
*	and stored. Compute mean and rms values for this
*	data file and continue.

	TOTAL=FLOAT(NUMREC*N)
	AVEU=SUM1UT/TOTAL
	RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)

*	Display the results for this file.

	WRITE (*,'(/16X,A,1X,A,G11.5,2X,A,G11.5/)') 'File Totals : ',
     +         'Ave = ',AVEU,'Rms = ',RMSU

	RETURN

40	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +                      'No such device or out of space !'
	TRANS_ONE_CHAN=1
	RETURN
	END

	FUNCTION TRANS_TWO_CHAN (NDATA,FR,FI,FR2,FI2,WIN,B,B2)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	INTEGER*2 	NDATA[HUGE](*)
	REAL*4		FR[HUGE](*)
	REAL*4		FR2[HUGE](*)
	REAL*4		FI[HUGE](*)
	REAL*4		FI2[HUGE](*)
	REAL*4		WIN[HUGE](*)
	REAL*4		B(*),B2(*)

*	Initialization.

	TRANS_TWO_CHAN=0
	SUM1UT=0.0
	SUM2UT=0.0
	SUM1VT=0.0
	SUM2VT=0.0

	N=N/2
	FN=FLOAT(N)

*	Perform mean and RMS value calculations for each record
*	in this data file.

	DO IB=1,NUMREC

*	  Fill the array with the current record.

	  IF (REALDAT) THEN
	    READ (NUMI) (WIN(I),I=1,2*N)
	    DO II=1,2*N,2
	      I=(II+1)/2
	      FR(I)=WIN(II)
	      FR2(I)=WIN(II+1)
	    ENDDO
	    DO I=1,N
	      FI(I)=FR(I)*FR(I)
	      FI2(I)=FR2(I)*FR2(I)
	    ENDDO
	  ELSE
	    READ (NUMI) (NDATA(I),I=1,2*N)
	  ENDIF

*	  Transform the data into voltages or velocities.
*	  Calculate squared values for RMS calculations.

	  IF (REALDAT) GO TO 20
	  DO II=1,2*N,2
	    I=(II+1)/2
	    V=FLOAT(NDATA(II))*VOFST1
	    V2=FLOAT(NDATA(II+1))*VOFST2
	    IF (VOLT) THEN
	      FR(I)=V
	      FI(I)=FR(I)*FR(I)
	      GO TO 10
	    ENDIF
	    FR(I)=B(5)
	    DO J=4,1,-1
	      FR(I)=FR(I)*V+B(J)
	    ENDDO
	    FI(I)=FR(I)*FR(I)
10	    IF (VOLT2) THEN
	      FR2(I)=V2
	      FI2(I)=FR2(I)*FR2(I)
	      CYCLE
	    ENDIF
	    FR2(I)=B2(5)
	    DO J=4,1,-1
	      FR2(I)=FR2(I)*V2+B2(J)
	    ENDDO
	    FI2(I)=FR2(I)*FR2(I)
	  ENDDO

*	  Determine mean and RMS values for each record.

20	  SUM1U=ADDSUM(FR,N)
	  SUM1UT=SUM1UT+SUM1U
	  SUM2U=ADDSUM(FI,N)
	  SUM2UT=SUM2UT+SUM2U
	  UBAR=SUM1U/FN
	  URMS=SQRT(SUM2U/FN-UBAR**2)

	  SUM1V=ADDSUM(FR2,N)
	  SUM1VT=SUM1VT+SUM1V
	  SUM2V=ADDSUM(FI2,N)
	  SUM2VT=SUM2VT+SUM2V
	  VBAR=SUM1V/FN
	  VRMS=SQRT(SUM2V/FN-VBAR**2)

*	  Display output for each record.

	  IF (IB .EQ. 1) THEN
	    WRITE (*,30) IB,UBAR,URMS,VBAR,VRMS
30	    FORMAT (1X,'Rec ',I4.4,2X,'Ave= ',G11.5,1X,
     +       'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
	  ELSE
	    WRITE (*,40) IB,UBAR,URMS,VBAR,VRMS
40	    FORMAT ('+','Rec ',I4.4,2X,'Ave= ',G11.5,1X,
     +       'Rms= ',G11.5,1X,'Ave2= ',G11.5,1X,'Rms2= ',G11.5)
	  ENDIF

*	  Save transformed values for this record.

	  IF (.NOT. NOFIL)
     +        CALL FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)
	  IF (RMEAN) THEN
	    DO I=1,N
	      FR(I)=FR(I)-UBAR
	    ENDDO
	  ENDIF
	  IF (DTREND) CALL DTRD (FR,N,DELT)
	  WRITE (NUMO,ERR=50) (FR(I), I=1,N)

	  IF (.NOT. NOFIL2)
     +        CALL FILTER (FTYP2,FC2,F12,F22,NS2,DELT,N,FR2)
	  IF (RMEAN2) THEN
	    DO I=1,N
	      FR2(I)=FR2(I)-VBAR
	    ENDDO
	  ENDIF
	  IF (DTREND2) CALL DTRD (FR2,N,DELT)
	  WRITE (NUMO2,ERR=50) (FR2(I), I=1,N)

	ENDDO

*	All records of this data file have been processed
*	and stored. Compute mean and rms values for this
*	data file and continue.

	TOTAL=FLOAT(NUMREC*N)
	AVEU=SUM1UT/TOTAL
	RMSU=SQRT(SUM2UT/TOTAL-AVEU**2)

	AVEV=SUM1VT/TOTAL
	RMSV=SQRT(SUM2VT/TOTAL-AVEV**2)

*	Display the results for this file.

	WRITE (*,'(/1X,A,2(A,G11.5,1X,A,G11.5,1X)/)')
     +         'File :    ','Ave= ',AVEU,'Rms= ',RMSU,
     +                    'Ave2= ',AVEV,'Rms2= ',RMSV

	RETURN

50	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +                      'No such device or out of space !'
	TRANS_TWO_CHAN=1
	RETURN
	END

	FUNCTION WIN_SCALE (WIN)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	REAL*4	WIN[HUGE](*)
	DATA	A0,A1,A2,A3 /0.355768,0.487396,0.144232,0.012604/
	DATA	C0,C1,C2,C3 /0.3125,0.46875,0.1875,0.03125/

*	Initialization.

	WIN_SCALE=0
	S=0.0
	TNM1=FLOAT(N-1)

*	Compute the scale factor.

	DO I=1,N
	  TIM1=FLOAT(I-1)
	  IF (HANN)
     +      WIN(I)=(1.0-COS(PI2*TIM1/TNM1))/2.0
	  IF (WELCH)
     +	    WIN(I)=1.0-((TIM1-TNM1/2.0)/(TNM1/2.0))**2
	  IF (PARZ)
     +	    WIN(I)=1.0-ABS((TIM1-TNM1/2.0)/(TNM1/2.0))
	  IF (TUKEY) THEN
	    IF (I .LE. K .OR. I .GE. N-K)
     +        WIN(I)=(1.0-COS(10.0*PI*TIM1/TNM1))/2.0
	    IF (I .GT. K .AND. I .LT. N-K) WIN(I)=1.0
	  ENDIF
	  IF (MINSL)
     +       WIN(I)=A0-A1*COS(PI2*TIM1/TNM1)
     +       +A2*COS(2.0*PI2*TIM1/TNM1)-A3*COS(3.0*PI2*TIM1/TNM1)
	  IF (MXDEC)
     +       WIN(I)=C0-C1*COS(PI2*TIM1/TNM1)
     +       +C2*COS(2.0*PI2*TIM1/TNM1)-C3*COS(3.0*PI2*TIM1/TNM1)
	  S=S+WIN(I)*WIN(I)
	ENDDO

	S=S/FLOAT(N)
	S=SQRT(1.0/S)

	RETURN

40	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +                      'No such device or out of space !'
	WIN_SCALE=1
	RETURN
	END

	FUNCTION TAPER_DATA (FR,FR2,FI,FI2,WIN)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	REAL*4		FR[HUGE](*)
	REAL*4		FR2[HUGE](*)
	REAL*4		FI[HUGE](*)
	REAL*4		FI2[HUGE](*)
	REAL*4		WIN[HUGE](*)

*	Initialization.

	TAPER_DATA=0
	IBO=0

*	Do the tapering now.

	DO IB=1,NUMREC+1

	  IF (IB .EQ. NUMREC+1 .AND. NOOVRLP) RETURN

*	  Reposition input file if last time thru loop.

	  IF (IB .EQ. NUMREC+1 .AND. OVRLP) THEN
	    REWIND (NUMO)
	    READ (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
	    IF (TWOCHAN) THEN
	      REWIND (NUMO2)
	      READ (NUMO2) ICHANS,IRSIZE,NUMREC,IDELTMS
	    ENDIF
	  ENDIF

*	  Get the next N points.

	  READ (NUMO) (FR(I),I=1,N)
	  IF (TWOCHAN) READ (NUMO2) (FR2(I),I=1,N)
	  IF (IB .EQ. 1 .OR. NOOVRLP) GO TO 10

*	  Operate on the first half of FR and store in FI.

	  DO I=N2+1,N
	    FI(I)=FR(I-N2)*S*WIN(I)
	    IF (TWOCHAN) FI2(I)=FR2(I-N2)*S*WIN(I)
	  ENDDO

*	  Write contents of FI to output file.

	  IBO=IBO+1
	  WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO

	  WRITE (NUMT,ERR=20) (FI(I),I=1,N)
	  IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FI2(I),I=1,N)
	  IF (IB .EQ. NUMREC+1) RETURN

*	  Operate on all of FR.

10	  DO I=1,N

*	    Operate on the second half of FR and store in FI.

	    IF (I .LE. N2 .AND. OVRLP) THEN
	      FI(I)=FR(N2+I)*S*WIN(I)
	      IF (TWOCHAN) FI2(I)=FR2(N2+I)*S*WIN(I)
	    ENDIF

*	    Now operate on all of FR.

	    FR(I)=FR(I)*S*WIN(I)
	    IF (TWOCHAN) FR2(I)=FR2(I)*S*WIN(I)

	  ENDDO

*	  Write contents of FR to output file.

	  IBO=IBO+1
	  IF (IBO .EQ. 1) THEN
	    WRITE (*,'(25X,''Tapering record '',I4.4)') IBO
	  ELSE
	    WRITE (*,'(''+'',24X,''Tapering record '',I4.4)') IBO
	  ENDIF

	  WRITE (NUMT,ERR=20) (FR(I),I=1,N)
	  IF (TWOCHAN) WRITE (NUMT2,ERR=20) (FR2(I),I=1,N)

	ENDDO

	RETURN

20	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +                      'No such device or out of space !'
	TAPER_DATA=1
	RETURN
	END

	FUNCTION SPEC_ONE_CHAN (FR,FI,FR2,FI2)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	REAL*4		FR[HUGE](*)
	REAL*4		FI[HUGE](*)
	REAL*4		FR2[HUGE](*)
	REAL*4		FI2[HUGE](*)

*	Initialization.

	SPEC_ONE_CHAN=0

*	Spectrum computations for one channel.

	DO IB=1,NRECD2

	  IF (IB .EQ. 1) THEN
	    WRITE (*,'(/25X,''Transforming records '',I4.4,
     +                 '' and '',I4.4)') 2*IB-1,2*IB
	  ELSE
	    WRITE (*,'(''+'',24X,''Transforming records '',I4.4,
     +                 '' and '',I4.4)') 2*IB-1,2*IB
	  ENDIF

	  IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2)
     +      WRITE (*,'(/25X,A,I4.4,A)')
     +      'Record ',2*IB,' was discarded.'

*	  Process two records simultaneously. This increases speed but
*	  requires an even number of records in the one channel case.

	  READ (INUM) (FR(I),I=1,N)

	  IF (ODDNUM .AND. IB .EQ. NRECD2) THEN
	    DO I=1,N
	      FI(I)=FR(I)
	    ENDDO
	  ELSE
	    READ (INUM) (FI(I),I=1,N)
	  ENDIF

	  CALL FFT (FR,FI,IOD)

	  IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) FI(1)=0.0

	  XMEAN=XMEAN+FR(1)
	  YMEAN=YMEAN+FI(1)

          DO K=1,N2

	    AK=FR(K+1)
	    BK=FI(K+1)
	    ANK=FR(NP1-K)
	    BNK=FI(NP1-K)

*	    Unscramble X and Y.

	    XR=AK+ANK
	    XI=BK-BNK
	    YR=BNK+BK
	    YI=ANK-AK

	    IF (POWER) THEN
	      B1K=XR*XR+XI*XI
	      B2K=YR*YR+YI*YI
	    ELSE
	      B1K=SQRT(XR*XR+XI*XI)
	      B2K=SQRT(YR*YR+YI*YI)
	    ENDIF

	    IF ( (OVRLP .OR. ODDNUM) .AND. IB .EQ. NRECD2) B2K=0.0

	    FR2(K)=FR2(K)+B1K
	    FI2(K)=FI2(K)+B2K

	  ENDDO
	ENDDO

	CLOSE (INUM,STATUS='DELETE',ERR=10)

	IB=IB-1		! Decrement added count in loop variable

*	Put results into one buffer.

	XMEAN=XMEAN+YMEAN
	DO I=1,N
	  FR2(I)=FR2(I)+FI2(I)
	ENDDO
	IB=IB*2
	IF (OVRLP .OR. ODDNUM) IB=IB-1	! Remove bad record.

*	Finished.

	FN=FLOAT(N)
	TOTAL=FN*FLOAT(IB)
	XMEAN=XMEAN/TOTAL
	DT=DELT
	DELF=1.0/(FN*DT)
	IF (POWER) SCALE=DT/(2.0*TOTAL)
	IF (AMP) THEN
	  SCALE=DT/(2.0*FLOAT(IB))
	  XMEAN=XMEAN*TOTAL*SCALE*2.0
	ENDIF

	DO I=1,N
	  FR2(I)=FR2(I)*SCALE
	ENDDO

	IF (POWER) XRMS=SQRT(DELF*ADDSUM(FR2,N2))

*	Put results into output [.PRN] file.

	IPD=INDEX (POUT,'.',.TRUE.)
	WRITE (*,'(/25X,''Writing '',A,'' now.'')') POUT(:IPD+3)

	OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=10)

	IF (POWER) THEN
	  WRITE (NUMP,'(F11.5,2X,F11.5/)',ERR=10) XMEAN,XRMS
	  DO I=1,N2
	    IF (FR2(I) .LT. EPS) FR2(I)=EPS
	    WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=10)
     +          DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I))
	  ENDDO
	ELSE
	  WRITE (NUMP,'(G15.5,2X,G15.5,2X,G15.5)',ERR=10)
     +           EPS,ALOG10(EPS),XMEAN
	  WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5)',ERR=10)
     +          (DELF*I,ALOG10(DELF*I),FR2(I), I=1,N2)
	ENDIF

	CLOSE (NUMP,STATUS='KEEP',ERR=10)

*	Display results on the screen.

	IF (POWER) THEN
	   WRITE (*,'(/29X,A,10X,A/21X,1P,2G15.5)')
     +                                    'mean','rms',XMEAN,XRMS
	   PERR = 1.0/SQRT(FLOAT(NREC))  ! use numrec before overlap
	   IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
	   WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
	ENDIF

	RETURN

10	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +                      'No such device or out of space !'
	SPEC_ONE_CHAN=1
	RETURN
	END

	FUNCTION SPEC_TWO_CHAN (FR,FI,FR2,FI2,WIN)

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	INCLUDE		'SPECTRUM.FN'	! Common statements

	REAL*4		FR[HUGE](*)
	REAL*4		FI[HUGE](*)
	REAL*4		FR2[HUGE](*)
	REAL*4		FI2[HUGE](*)
	REAL*4		WIN[HUGE](*)

*	Initialization.

	SPEC_TWO_CHAN=0

*	Spectrum computations for two channels.

	DO IB=1,NUMREC

	  IF (IB .EQ. 1) THEN
	    WRITE (*,'(/25X,''Transforming record '',I4.4)') IB
	  ELSE
	    WRITE (*,'(''+'',24X,''Transforming record '',I4.4)') IB
	  ENDIF

	  IF (OVRLP .AND. IB .EQ. NUMREC)
     +      WRITE (*,'(/25X,A,I4.4,A)')
     +      'Record ',IB+1, ' was discarded for each channel.'

	  READ (INUM) (FR(I),I=1,N)
	  READ (INUM2) (FI(I),I=1,N)

	  CALL FFT (FR,FI,IOD)

	  XMEAN=XMEAN+FR(1)
	  YMEAN=YMEAN+FI(1)

          DO K=1,N2

	    AK=FR(K+1)
	    BK=FI(K+1)
	    ANK=FR(NP1-K)
	    BNK=FI(NP1-K)

*	    Unscramble X and Y.

	    XR=AK+ANK
	    XI=BK-BNK
	    YR=BNK+BK
	    YI=ANK-AK

	    IF (POWER) THEN
	      B1K=XR*XR+XI*XI
	      B2K=YR*YR+YI*YI
	    ELSE IF (CROSS) THEN
	      B1K=XR*YR+XI*YI
	      B2K=XI*YR-XR*YI
	      IF (COH) THEN
	        B3K=XR*XR+XI*XI
	        B4K=YR*YR+YI*YI
	      ENDIF
	    ELSE
	      B1K=SQRT(XR*XR+XI*XI)
	      B2K=SQRT(YR*YR+YI*YI)
	    ENDIF

	    FR2(K)=FR2(K)+B1K
	    FI2(K)=FI2(K)+B2K

	    IF (COH) THEN
	      WIN(K)=WIN(K)+B3K
	      WIN(N2+K)=WIN(N2+K)+B4K
	    ENDIF

	  ENDDO
	ENDDO

	CLOSE (INUM,STATUS='DELETE',ERR=20)
	CLOSE (INUM2,STATUS='DELETE',ERR=20)

	IB=IB-1		! Decrement added count in loop variable

*	Finished.

	FN=FLOAT(N)
	TOTAL=FN*FLOAT(IB)
	XMEAN=XMEAN/TOTAL
	YMEAN=YMEAN/TOTAL
	DT=DELT
	DELF=1.0/(FN*DT)
	IF (POWER .OR. CROSS) SCALE=DT/(2.0*TOTAL)
	IF (AMP) THEN
	  SCALE=DT/(2.0*FLOAT(IB))
	  XMEAN=XMEAN*TOTAL*SCALE*2.0
	  YMEAN=YMEAN*TOTAL*SCALE*2.0
	ENDIF

	DO I=1,N
	  FR2(I)=FR2(I)*SCALE
	  FI2(I)=FI2(I)*SCALE
	ENDDO

	IF (COH) THEN
	  DO I=1,N
	    WIN(I)=WIN(I)*SCALE
	  ENDDO
	ENDIF

	IF (POWER) THEN
	  XRMS=SQRT(DELF*ADDSUM(FR2,N2))
	  YRMS=SQRT(DELF*ADDSUM(FI2,N2))
	ENDIF

	IF (CROSS) THEN
	  DEGRAD=360.0/PI2
	  S2=SQRT(2.0)
	  NR=NREC                ! use numrec before overlap
	  SNR=SQRT(FLOAT(NR))
	  S2NR=SQRT(FLOAT(NR)*2.0)
	  RXY0=DELF*ADDSUM(FR2,N2)
	  DO I=1,N2
	    FII=FLOAT(I)
	    GXYMAG=SQRT(FR2(I)**2+FI2(I)**2)
	    PHASE=ATAN2(FI2(I),FR2(I))
	    FR2(I)=GXYMAG
	    FI2(I)=PHASE
	    FI2(N2+I)=PHASE

*	    If coherence is calculated then compute errors associated
*	    with cross spectrum estimates.

	    IF (COH) THEN
	      GAMMA=FR2(I)**2/(WIN(I)*WIN(N2+I))
	      WIN(I)=GAMMA
	      FR(I)=1.0/(SQRT(WIN(I))*SNR)                  ! e [ |Gxy(f)| ]
	      FR(N2+I)=SQRT(1.0-WIN(I))/(SQRT(WIN(I))*S2NR) ! s.d. [ Oxy(f) ]
	      FR(N2+I)=FR(N2+I)*DEGRAD                      ! s.d. [O] in deg
	      FI(I)=S2*(1.0-WIN(I))/(SQRT(WIN(I))*SNR)      ! e [ gxy2(f) ]
	      IF (TAPER .AND. NOOVRLP) THEN
	        FR(I)=FR(I)*S2
	        FR(N2+I)=FR(N2+I)*S2
	        FI(I)=FI(I)*S2
	      ENDIF
	    ENDIF
	  ENDDO

*	  Apply a phase shift to straighten out Theta only if
*	  coherence is calculated -- need s.d. estimate.

	  IF (COH) THEN

	    NSD2=(NSD-1)/2

*	    Check on proper relationship of NPTS and NSD.

	    IF (NSD2+1 .GT. NPTS) THEN
	      WRITE (*,'(/21X,A,A/21X,A)') 'NSD too big or NPTS',
     +        ' too small.','No adjustment of phase angle performed.'
	      GO TO 10
	    ENDIF

*	    Get first NPTS points. Pass this to fit subroutine.

	    DO I=1,NPTS
	      X(I)=DELF*I
	      Y(I)=FI2(N2+I)
	    ENDDO

*	    Sum first NSD points to check on average s.d. of theta.

	    SDSUM=0.0
	    DO I=1,NSD
	      SDSUM=SDSUM+FR(N2+NPTS-NSD2-1+I)
	    ENDDO

	    DO I=NPTS+1,N2

*	      Check on moving average of s.d. of theta.
*	      If too high -- quit.

	      SDSUM=SDSUM+FR(N2+I+NSD2)
	      SDSUM=SDSUM-FR(N2+I-NSD2-1)
	      SDAVG=SDSUM/NSD
	      IF (SDAVG .GT. SDLIM) GO TO 10

*	      Perform linear least squares fit on NPTS points before
*	      the current point. Guess value of current point from fit.

	      CALL LSF (X,Y,NPTS,XA1,XB2)
	      TEST=XA1+XB2*DELF*I

*	      XLO & XHI are current point guess +/- pi.

	      XLO=TEST-PI
	      XHI=TEST+PI

*	      If actual value is .GT. pi away from guessed value then
*	      add 2*pi to or subtract 2*pi from actual value as required.

	      DO WHILE (FI2(N2+I) .LT. XLO)
	        FI2(N2+I)=FI2(N2+I)+PI2
	      ENDDO

	      DO WHILE (FI2(N2+I) .GT. XHI)
	        FI2(N2+I)=FI2(N2+I)-PI2
	      ENDDO

*	      Prepare set of NPTS points before next current point.

	      DO J=1,NPTS
	        X(J)=DELF*(I-NPTS+J)
	        Y(J)=FI2(N2+I-NPTS+J)
	      ENDDO

	    ENDDO
	  ENDIF
	ENDIF

*	Put results into output [.PRN] file.

10	IPD=INDEX (POUT,'.',.TRUE.)
	WRITE (*,'(/25X,''Writing '',A,'' now.''/)') POUT(:IPD+3)

	OPEN (NUMP,FILE=POUT,STATUS='UNKNOWN',ERR=20)

	IF (POWER) THEN
	  WRITE (NUMP,'(4(F11.5,2X))',ERR=20) XMEAN,XRMS,YMEAN,YRMS
	  WRITE (NUMP,'( )',ERR=20)
	  DO I=1,N2
	    IF (FR2(I) .LT. EPS) FR2(I)=EPS
	    IF (FI2(I) .LT. EPS) FI2(I)=EPS
	    WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,
     +        F11.5,1X,G15.5,1X,F11.5)',ERR=20)
     +        DELF*I,ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),
     +        FI2(I),ALOG10(FI2(I))
	  ENDDO
	ELSE IF (CROSS) THEN
	  IF (NOCOH) THEN
	    WRITE (NUMP,'(G15.5)',ERR=20) RXY0
	    WRITE (NUMP,'( )',ERR=20)
	    WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,
     +        1X,F11.5)',ERR=20) (DELF*I,ALOG10(DELF*I),
     +        FR2(I),ALOG10(FR2(I)),FI2(I),I=1,N2)
	  ELSE IF (COH) THEN
	    WRITE (NUMP,'(G15.5)',ERR=20) RXY0
	    WRITE (NUMP,'( )',ERR=20)
	    WRITE (NUMP,'(F11.5,1X,F11.5,1X,G15.5,1X,F11.5,1X,
     +        F11.5,1X,F11.5,1X,G15.5)',ERR=20) (DELF*I,
     +        ALOG10(DELF*I),FR2(I),ALOG10(FR2(I)),FI2(I),FI2(N2+I),
     +        WIN(I),I=1,N2)
	  ENDIF
	ELSE
	  WRITE (NUMP,'(4(G15.5,2X))',ERR=20)
     +      EPS,ALOG10(EPS),XMEAN,YMEAN
	  WRITE (NUMP,'(F11.5,2X,G15.5,2X,G15.5,2X,G15.5)',ERR=20)
     +      (DELF*I,ALOG10(DELF*I),FR2(I),FI2(I),I=1,N2)
	ENDIF

	CLOSE (NUMP,STATUS='KEEP',ERR=20)

*	Display results on the screen.

	IF (POWER) THEN
	   WRITE (*,'(29X,A,10X,A,2(/21X,1P,2G15.5))')
     +       'mean','rms',XMEAN,XRMS,YMEAN,YRMS
	   PERR = 1.0/SQRT(FLOAT(NREC))  ! use numrec before overlap
	   IF (TAPER .AND. NOOVRLP) PERR=PERR*SQRT(2.0)
	   WRITE (*,'(/25X,A,F6.5)') 'norm. random error = ',PERR
	ENDIF

	IF (CROSS) WRITE (*,'(21X,A,1P,G15.5)')
     +    'Rxy(0) = E [x(t) y(t)] = ',RXY0

*	If coherence was calculated then store
*	error estimates into separate error file.

	IF (COH) THEN
	  IPD=INDEX (EOUT,'.',.TRUE.)
	  WRITE (*,'(/25X,''Writing '',A,'' now.'')') EOUT(:IPD+3)

	  OPEN (NUME,FILE=EOUT,STATUS='UNKNOWN',ERR=20)

	  WRITE (NUME,'(5(F11.5,2X))',ERR=20) (DELF*I,
     +      ALOG10(DELF*I),FR(I),FR(N2+I),FI(I),I=1,N2)

	  CLOSE (NUME,STATUS='KEEP',ERR=20)
	ENDIF

	RETURN

20	WRITE (*,'(/13X,A,A/)') 'Run-time I/O error : ',
     +    'No such device or out of space !'
	SPEC_TWO_CHAN=1
	RETURN
	END

********************************************************
	SUBROUTINE FFT (FR,FI,K)
*                                                      *	
*	Fast Fourier Transform of a                    *
*	Complex Time Series                            *
*	Features : time decomposition (odd-even)       *
*	         : input bit reversal                  *
*	         : computation 'in place'              *
*	         : N= number of points                 *
*	         : FR+j*FI - complex series            *
*	Source   : Stearns 1975                        *
********************************************************

	IMPLICIT	INTEGER*4 (I-N)
	REAL*4		FR[HUGE](*)
	REAL*4		FI[HUGE](*)

	N=2**K
	MR=0
	NN=N-1

*	Input bit reversal.

	DO M=1,NN
	  L=N
10	  L=L/2
	  IF(MR+L .GT. NN) GO TO 10
	  MR=MOD(MR,L)+L
	  IF (MR.LE. M) CYCLE
	  TR=FR(M+1)
	  FR(M+1)=FR(MR+1)
	  FR(MR+1)=TR
	  TI=FI(M+1)
	  FI(M+1)=FI(MR+1)
	  FI(MR+1)=TI
	ENDDO

*	Compute the coefficients.

	L=1
	PI=2.0*ASIN(1.0)
20	IF (L.GE.N) RETURN
	ISTEP=2*L
	EL=L

	DO M=1,L
	  A=PI*FLOAT(1-M)/EL
	  WR=COS(A)
	  WI=SIN(A)
	  DO I=M,N,ISTEP
	    J=I+L
	    TR=WR*FR(J)-WI*FI(J)
	    TI=WR*FI(J)+WI*FR(J)
	    FR(J)=FR(I)-TR
	    FI(J)=FI(I)-TI
	    FR(I)=FR(I)+TR
	    FI(I)=FI(I)+TI
	  ENDDO
	ENDDO

	L=ISTEP
	GO TO 20
	END

	SUBROUTINE RDCFG (ROOT1,ROOT2,SUFX,PTH,TPTH,SUFFIX)

*	Routine to read config file - spectrum.cfg

	PARAMETER	(NUM=2)
	LOGICAL*1	ROOT1,ROOT2,SUFX
	CHARACTER*1	FIRST
	CHARACTER*4	SUFFIX
	CHARACTER*72	FLINE,PTH,TPTH

*	File exists. Get various paths.

	KOUNT=0
	OPEN (NUM,FILE='SPECTRUM.CFG',STATUS='OLD')
	  IF (EOF(NUM)) GO TO 30
10	  READ (NUM,'(A)',ERR=30) FLINE

*	  Check for a comment line.

	  IF (FLINE (1:1) .EQ. '*') GO TO 10
	  KOUNT=KOUNT+1

*	  Convert to upper case characters if necessary.

	  DO J=72,1,-1
	    FIRST=FLINE(J:J)
	    IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
	      IHOLD=ICHAR(FIRST)-32
	      FIRST=CHAR(IHOLD)
	      FLINE(J:J)=FIRST
	    ENDIF
	  ENDDO

*	  Assign appropriate paths.

	  IF (KOUNT .EQ. 1) THEN
	    IF (FLINE(:4) .EQ. 'ROOT') THEN
	      ROOT1=.TRUE.
	    ELSE
	      ROOT1=.FALSE.
	      ISL=INDEX (FLINE,'\',.TRUE.)
	      PTH=FLINE(:ISL)
	    ENDIF
	  ELSE IF (KOUNT .EQ. 2) THEN
	    IF (FLINE(:4) .EQ. 'ROOT') THEN
	      ROOT2=.TRUE.
	    ELSE
	      ROOT2=.FALSE.
	      ISL=INDEX (FLINE,'\',.TRUE.)
	      TPTH=FLINE(:ISL)
	    ENDIF
	  ELSE IF (KOUNT .EQ. 3) THEN
	    SUFX=.TRUE.
	    SUFFIX=FLINE(:4)
	  ENDIF
	  GO TO 10

30	CLOSE (NUM)
	RETURN
	END

	SUBROUTINE LSF(X,Y,N,A,B)

	IMPLICIT REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	REAL*4 X(N),Y(N)

*	Initialization.

	XS=0.0
	X2S=0.0
	YS=0.0
	XYS=0.0

*	Compute sums.

	DO I=1,N
	  XS=XS+X(I)
	  X2S=X2S+X(I)*X(I)
	  YS=YS+Y(I)
	  XYS=XYS+X(I)*Y(I)
	ENDDO

*	Compute coefficients and guess next value.

	DEN=N*X2S-XS*XS
	A=(X2S*YS-XS*XYS)/DEN
	B=(N*XYS-XS*YS)/DEN

	RETURN
	END

	FUNCTION ADDSUM (Y,N)

*	SUMMING of elements of real*4 array Y. This function uses
*	a tree summming process to minimize truncation error for
*	large N(>500). A Real*4 and Integer*2 buffer of size log2(n)
*	is required. This version received from Dr. Anuvat Sirivat.

	REAL*4      Y[HUGE](*)
	REAL*4      B(16),ADDSUM
	INTEGER*2   W(16)
	INTEGER*4   N

	K=1
	W(1)=0
	B(1)=0.0

	DO I=1,N
	  B(K)=B(K)+Y(I)
	  W(K)=W(K)+1
10	  K1=K-1

	  IF(K1.LE.0) GO TO 20
	  IF(W(K1).GT.W(K)) GO TO 20

	  B(K1)=B(K1)+B(K)
	  W(K1)=W(K1)+W(K)
	  K=K-1
	  GO TO 10

20	  IF(W(K).LT.2) CYCLE

	  K=K+1
	  W(K)=0
	  B(K)=0.0
	ENDDO

30	K2=K/2
	IF(K2.EQ.0) GO TO 40

	J=1
	DO I=1,K2
	  B(I)=B(J)+B(J+1)
	  J=J+2
	ENDDO

        IF(J.EQ.K) B(1)=B(1)+B(K)
        K=K2
	GO TO 30

40	ADDSUM=B(1)
	RETURN
	END

	SUBROUTINE DTRD (FR,N,DELT)

*	This routine will remove a linear mean trend from the
*	original data by means of least squares fit. This is a
*	modified version of a routine received from Dr. Sirivat.

	IMPLICIT REAL*4 (A-H,O-Z), INTEGER*4 (I-N)
	INTEGER*4	N
	REAL*4		DELT
	REAL*4		FR[HUGE](*)

*	Initialization.

	FN=FLOAT(N)
	YSR=0.0
	XYSR=0.0

*	Compute required sums for least squares fit.

	DO I=1,N
	  FII=FLOAT(I)
	  YSR=YSR+FR(I)
	  XYSR=XYSR+FII*FR(I)
	ENDDO

	XS=FN*(FN+1.0)/2.0*DELT
	X2S=FN*(FN+1.0)*(2.0*FN+1.0)/6.0*DELT*DELT
	XYSR=XYSR*DELT

	DER=FN*X2S-XS*XS

*	Compute coefficients and remove trend.

	B0R=(YSR*X2S-XYSR*XS)/DER
	B1R=(FN*XYSR-XS*YSR)/DER

 	DO I=1,N
	  FII=FLOAT(I)
	  FR(I)=FR(I)-B0R-B1R*FII*DELT
	ENDDO

	RETURN
	END

	SUBROUTINE FILTER (FTYPE,FC,F1,F2,NS,DELT,N,FR)

*	Routine to filter a data file.  (FILTER.FOR)

*	This routine will filter a data file using any one of
*	four different types of IIR recursive digital filters:
*	Lowpass, Highpass, Bandpass, and Bandstop. The routines
*	are modified versions of those given in Stearns, 1975.

	PARAMETER 	(NSMAX=20)
	LOGICAL*1 	LPASS,HPASS,BPASS,BSTOP
	INTEGER*2 	NS
	INTEGER*4 	N,J
	REAL*4 		FR[HUGE](*)
	REAL*4 		X(NSMAX+1,5)
	REAL*4 		A(NSMAX),B(NSMAX),C(NSMAX),D(NSMAX),E(NSMAX),G
	CHARACTER*1 	FTYPE

*	Initialization.

	X=0.0
	G=0.0

	LPASS=(FTYPE .EQ. 'L')
	HPASS=(FTYPE .EQ. 'H')
	BPASS=(FTYPE .EQ. 'B')
	BSTOP=(FTYPE .EQ. 'S')

*	Design the desired filter.

	IF (LPASS .OR. HPASS) KHI=3
	IF (BPASS .OR. BSTOP) KHI=5

	IF (LPASS) CALL LP (FC,DELT,NS,A,B,C)
	IF (HPASS) CALL HP (FC,DELT,NS,A,B,C)
	IF (BPASS) CALL BP (F1,F2,DELT,NS,A,B,C,D,E)
	IF (BSTOP) CALL BS (F1,F2,DELT,NS,A,B,C,D,E,G)
	G2=2.0*G
	GG=2.0+G*G

*	Filter the array FR(J).

	DO J=1,N

*	  Go thru NS sections of filter.

	  IF (LPASS) THEN
	    X(1,3)=FR(J)
	    DO I=1,NS
	      TEMP=A(I)*(X(I,3)+2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
	      X(I+1,3)=TEMP-C(I)*X(I+1,1)
	    ENDDO
	  ELSE IF (HPASS) THEN
	    X(1,3)=FR(J)
	    DO I=1,NS
	      TEMP=A(I)*(X(I,3)-2.0*X(I,2)+X(I,1))-B(I)*X(I+1,2)
	      X(I+1,3)=TEMP-C(I)*X(I+1,1)
	    ENDDO
	  ELSE IF (BPASS) THEN
	    X(1,5)=FR(J)
	    DO I=1,NS
	      TEMP=A(I)*(X(I,5)-2.0*X(I,3)+X(I,1))-B(I)*X(I+1,4)
	      X(I+1,5)=TEMP-C(I)*X(I+1,3)-D(I)*X(I+1,2)-E(I)*X(I+1,1)
	    ENDDO
	  ELSE IF (BSTOP) THEN
	    X(1,5)=FR(J)
	    DO I=1,NS
	      TEMP=A(I)*(X(I,5)+G2*X(I,4)+GG*X(I,3)+G2*X(I,2)+X(I,1))
	      X(I+1,5)=TEMP-B(I)*X(I+1,4)-C(I)*X(I+1,3)-D(I)*X(I+1,2)
     +                 -E(I)*X(I+1,1)
	    ENDDO
	  ENDIF

*	  Push down past values in filter.

	  DO I=1,NS+1
	    DO K=1,KHI-1
	      X(I,K)=X(I,K+1)
	    ENDDO
	  ENDDO

*	  Set FR(J) equal to the output of the filter and continue.

	  FR(J)=X(NS+1,KHI)

	ENDDO

	RETURN
	END

	SUBROUTINE LP (FC,T,NS,A,B,C)

*	Lowpass Butterworth Digital Filter 
*	Design Subroutine (LP.FOR)

	INTEGER*2 	NS
	REAL*4 		A(*),B(*),C(*)

	PI=2.0*ASIN(1.0)

	WCP=SIN(FC*PI*T)/COS(FC*PI*T)

	DO K=1,NS
	  CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
	  X=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
	  A(K)=WCP*WCP*X
	  B(K)=2.0*(WCP*WCP-1.0)*X
	  C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*X
	ENDDO

	RETURN
	END

	SUBROUTINE HP (FC,T,NS,A,B,C)

*	Highpass Butterworth Digital Filter 
*	Design Subroutine (HP.FOR)

	INTEGER*2 	NS
	REAL*4 		A(*),B(*),C(*)

	PI=2.0*ASIN(1.0)

	WCP=SIN(FC*PI*T)/COS(FC*PI*T)

	DO K=1,NS
	  CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
	  A(K)=1.0/(1.0+WCP*WCP-2.0*WCP*CS)
	  B(K)=2.0*(WCP*WCP-1.0)*A(K)
	  C(K)=(1.0+WCP*WCP+2.0*WCP*CS)*A(K)
	ENDDO

	RETURN
	END

	SUBROUTINE BP (F1,F2,T,NS,A,B,C,D,E)

*	Bandpass Butterworth Digital Filter 
*	Design Subroutine (BP.FOR)

	INTEGER*2 	NS
	REAL*4 		A(*),B(*),C(*),D(*),E(*)

	PI=2.0*ASIN(1.0)

	W1=SIN(F1*PI*T)/COS(F1*PI*T)
	W2=SIN(F2*PI*T)/COS(F2*PI*T)
	WC=W2-W1
	Q=WC*WC+2.0*W1*W2
	S=W1*W1*W2*W2

	DO K=1,NS
	  CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
	  P=-2.0*WC*CS
	  R=P*W1*W2
	  X=1.0+P+Q+R+S
	  A(K)=WC*WC/X
	  B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
	  C(K)=(6.0-2.0*Q+6.0*S)/X
	  D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
	  E(K)=(1.0-P+Q-R+S)/X
	ENDDO

	RETURN
	END

	SUBROUTINE BS (F1,F2,T,NS,A,B,C,D,E,G)

*	Bandstop Butterworth Digital Filter 
*	Design Subroutine (BS.FOR)

	INTEGER*2 	NS
	REAL*4 		A(*),B(*),C(*),D(*),E(*),G

	PI=2.0*ASIN(1.0)

	W1=SIN(F1*PI*T)/COS(F1*PI*T)
	W2=SIN(F2*PI*T)/COS(F2*PI*T)
	WC=W2-W1
	XK=1.0+W1*W2
	G=(2.0*XK-4.0)/XK
	Q=WC*WC+2.0*W1*W2
	S=W1*W1*W2*W2

	DO K=1,NS
	  CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
	  P=-2.0*WC*CS
	  R=P*W1*W2
	  X=1.0+P+Q+R+S
	  A(K)=XK*XK/X
	  B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
	  C(K)=(6.0-2.0*Q+6.0*S)/X
	  D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
	  E(K)=(1.0-P+Q-R+S)/X
	ENDDO

	RETURN
	END
