*	COLOR.FOR

*	Program to produce colored noise for spectrum test.

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

	IMPLICIT	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	PARAMETER 	(NUMO=2,NSMAX=20,NMAX=16384)
	INTEGER*2 	NDATA [ALLOCATABLE,HUGE] (:)
	INTEGER*2	GAIN(0:7)
	INTEGER*4 	IRSIZE,N,J,ITST,IDELTMS
	REAL*4 		X(NSMAX+1,5)
	REAL*4 		A(NSMAX),B(NSMAX),C(NSMAX)
	REAL*4		D(NSMAX),E(NSMAX),GRAF(2,20)
	CHARACTER*1	FIRST
	CHARACTER*4 	FLNM
	CHARACTER*8 	FILNAM

*	This routine uses a random number generator to produce uniform
*	deviates between 0 and 1.  These are then transformed to a 
*	sequence having a power density = 1 with power concentrated
*	between f1 and f2 Hz for a sampling interval of dt secs.

*	Initialization.

	ICHANS=1
	GAIN=0
	VOFST=20.0/4096.0

*	Get bandpass critical (-3 dB) points.

	WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
	READ (*,*) F1
	WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
	READ (*,*) F2

*	Get the sampling interval.

10	WRITE (*,'(/1X,A\)') 'Enter delta T secs. (1.0/Sample Rate) : '
	READ (*,*) DELT
	  XTST=SQRT(6.0/DELT)*0.5
	  ITST=INT4(XTST/VOFST)
	  IF (ITST .GT. 32767) THEN
	    WRITE (*,'(//1X,A,F7.6,A)')
     +        'This choice of delta t = ',DELT,' ,'
	    WRITE (*,'(1X,A,F6.2,A)')
     +        'leads to a maximum data value of ',XTST,' ,'
	    WRITE (*,'(1X,A,A,I6,A)')
     +        'which, when converted to an integer,',
     +        ' has a value of ',ITST,' .'
	    WRITE (*,'(1X,A/)')
     +        'This cannot be stored in the Integer*2 array NDATA.'
	    GO TO 10
	  ENDIF

*	Get order of filter.

20	WRITE (*,'(/1X,A\)') '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 20
	ENDIF

*	Get number of points per record.

30	WRITE (*,'(/1X,A\)') 'Enter # of points per record (N) : '
	READ (*,*) N
	IF (N .GT. NMAX) THEN
	  WRITE (*,'(1X,A,I5)') 'Reduce the # of points to ',NMAX
	  GO TO 30
	ENDIF

*	Allocate space for NDATA array.

	ALLOCATE (NDATA(N), STAT=IERR)
	IF (IERR .NE. 0)
     +          STOP 'Not enough storage for data.  Aborting ...'

*	Get number of records in file.

	WRITE (*,'(/1X,A\)') 'Enter # of records (NUMREC) : '
	READ (*,*) NUMREC

*	Get seed value for random number generator.

	WRITE (*,'(/1X,A\)') 'Enter a negative integer : '
	READ (*,*) IDUM

*	Get output file name.

40	WRITE (*,'(/1X,A\)') 'Enter output file name (4 chars) : '
	READ (*,'(A)') FLNM
	WRITE (*,'( )')

*	Convert to uppercase and check first character alphabetic.

	DO J=4,1,-1
	  FIRST=FLNM(J:J)
	  IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
	    IHOLD=ICHAR(FIRST)-32
	    FIRST=CHAR(IHOLD)
	    FLNM(J:J)=FIRST
	  ENDIF
	ENDDO
	IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
	  WRITE (*,'(/1X,A,A,A/1X,A,A,A/1X,A)') 
     +      'Filename ',FLNM,' began with',
     +      'the nonalphabetic character ',FIRST,'.',
     +      'Re-enter the filename correctly.'
	  GO TO 40
	ENDIF

	FILNAM=FLNM // '.DAT'
	IRSIZE=ICHANS*N*2
	IDELTMS=NINT(DELT*1.0E6)

*	Design the bandpass filter.

	CALL BPDES (F1,F2,DELT,NS,A,B,C,D,E,GRAF)

*	Initialize past values in filter to zero.

	X=0.0

*	Write the data in the form of binary numbers to a data
*	file that may be read by the spectrum routine.

	OPEN (NUMO,FILE=FILNAM,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
     +        FORM='BINARY')
	WRITE (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
	WRITE (NUMO) (GAIN(I),I=0,7)

*	Put message on screen.

	WRITE (*,'(/////////////////////16X,
     +           ''C O L O R E D   N O I S E   U T I L I T Y'')')
	WRITE (*,'(/25X,''Creating '',A,'' now.''/)') FILNAM 

*	Create NUMREC records of the sequence.

	DO IB=1,NUMREC

	  IF (IB .EQ. 1) ITOSS=1000
	  IF (IB .NE. 1) ITOSS=0

*	  Generate the sequence NDATA(J).

	  DO J=1,N+ITOSS

	    X(1,5)=SQRT(6.0/DELT)*(RAN1(IDUM)-0.5)

*	    Go thru NS sections of filter.

	    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

*	    Push down past values in filter.

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

*	    Throw away ITOSS values to eliminate startup transient.

	    IF (IB .EQ. 1) THEN
	      IF (J .LE. ITOSS) CYCLE
	      NDATA(J-ITOSS)=INT2(X(NS+1,5)/VOFST)
	      CYCLE
	    ENDIF

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

	    NDATA(J)=INT2(X(NS+1,5)/VOFST)

	  ENDDO

*	  Display record number message.

	  IF (IB .EQ. 1) THEN
	    WRITE (*,50) IB
50	    FORMAT (25X,'Writing Record ',I4.4)
	  ELSE
	    WRITE (*,60) IB
60	    FORMAT ('+',24X,'Writing Record ',I4.4)
	  ENDIF

*	  Create file for time series.

	  IF (IB .EQ. 1) THEN
	    OPEN (3,FILE='COLRTIM.PRN',STATUS='UNKNOWN')
	      WRITE (3,'(G15.5,2X,G15.5)')
     +          (FLOAT(J-1)*0.01,FLOAT(NDATA(J))*VOFST,J=1,N)
	    CLOSE (3)
	  ENDIF

*	  Save output to a file.

	  WRITE (NUMO) (NDATA(J),J=1,N)

	ENDDO

	CLOSE (NUMO,STATUS='KEEP')

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

*	Bandpass filter routine.   (BPDES.FOR)

*	Bandpass Butterworth Digital Filter Design Subroutine

*	Inputs are passband (-3 dB) frequencies F1 and F2 in Hz,
*			Sampling interval T in seconds, and
*			number NS of filter sections.

*	Outputs are NS sets of filter coefficients, i.e.,
*			A(K) thru E(K) for K = 1 thru NS, and
*			20 pairs of frequency and power gain, i.e.,
*			GRAF(1,K) thru GRAF(2,K) FOR K = 1 thru 20.

*	Note that A(K) thru E(K) as well as GRAF(2,20) must be
*			dimensioned in the calling program.

*	The digital filter has NS sections in cascade.  The Kth
*			section has the transfer function:

*				A(K)*(Z**4-2*Z**2+1)
*		H(Z) =  ------------------------------------
*			Z**4+B(K)*Z**3+C(K)*Z**2+D(K)*Z+E(K)

*	Thus, if F(M) and G(M) are the input and output of the 
*			Kth section at time M*T, then

*		G(M)=A(K)*(F(M)-2*F(M-2)+F(M-4))-B(K)*G(M-1)
*			-C(K)*G(M-2)-D(K)*G(M-3)-E(K)*G(M-4)

*	Do not use this routine for lowpass design.  Use LPDES.
*	Do not use this routine for highpass design.  Use HPDES.

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

	IMPLICIT 	REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
	REAL*4 		A(*),B(*),C(*),D(*),E(*),GRAF(2,20)

*	PI=2.0*ASIN(1.0)
	PI=3.1415926536

	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 150 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
150	CONTINUE

	DO 170 J=1,2
	  DO 160 I=1,10
	    K=I*(2-J)+(21-I)*(J-1)
	    GRAF(2,K)=0.01+0.98*FLOAT(I-1)/9.0
	    X=(1.0/GRAF(2,K)-1.0)**(1.0/FLOAT(4*NS))
	    SQ=SQRT(WC*WC*X*X+4.0*W1*W2)
	    GRAF(1,K)=ABS(ATAN(0.5*(WC*X+FLOAT(2*J-3)*SQ)))/(PI*T)
160	  CONTINUE
170	CONTINUE

	RETURN
	END

	FUNCTION RAN1(IDUM)

*	Returns a uniform random deviate between 0.0 and 1.0.  Set
*	IDUM to any negative value to initialize or reinitialize the
*	sequence. Taken from Numerical Recipes, p.196-197.

	INTEGER*2	IDUM
	REAL*4		R(97)
	PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
	PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
	PARAMETER (M3=243000,IA3=4561,IC3=51349)
	DATA		IFF /0/

	IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
	  IFF=1
	  IX1=MOD(IC1-IDUM,M1)
	  IX1=MOD(IA1*IX1+IC1,M1)
	  IX2=MOD(IX1,M2)
	  IX1=MOD(IA1*IX1+IC1,M1)
	  IX3=MOD(IX1,M3)
	  DO J=1,97
	    IX1=MOD(IA1*IX1+IC1,M1)
	    IX2=MOD(IA2*IX2+IC2,M2)
	    R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
	  ENDDO
	  IDUM=1
	ENDIF

	IX1=MOD(IA1*IX1+IC1,M1)
	IX2=MOD(IA2*IX2+IC2,M2)
	IX3=MOD(IA3*IX3+IC3,M3)
	J=1+(97*IX3)/M3
	IF(J.GT.97.OR.J.LT.1)PAUSE
	RAN1=R(J)
	R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

	RETURN
	END
