home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / MISC / MATH / SPECTR20.ZIP / COLOR.FOR < prev    next >
Encoding:
Text File  |  1992-04-22  |  8.1 KB  |  326 lines
  1. *    COLOR.FOR
  2.  
  3. *    Program to produce colored noise for spectrum test.
  4.  
  5. *    David E. Hess
  6. *    Fluid Flow Group - Process Measurements Division
  7. *    Chemical Science and Technology Laboratory
  8. *    National Institute of Standards and Technology
  9. *    April 15, 1992
  10.  
  11.     IMPLICIT    REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  12.     PARAMETER     (NUMO=2,NSMAX=20,NMAX=16384)
  13.     INTEGER*2     NDATA [ALLOCATABLE,HUGE] (:)
  14.     INTEGER*2    GAIN(0:7)
  15.     INTEGER*4     IRSIZE,N,J,ITST,IDELTMS
  16.     REAL*4         X(NSMAX+1,5)
  17.     REAL*4         A(NSMAX),B(NSMAX),C(NSMAX)
  18.     REAL*4        D(NSMAX),E(NSMAX),GRAF(2,20)
  19.     CHARACTER*1    FIRST
  20.     CHARACTER*4     FLNM
  21.     CHARACTER*8     FILNAM
  22.  
  23. *    This routine uses a random number generator to produce uniform
  24. *    deviates between 0 and 1.  These are then transformed to a 
  25. *    sequence having a power density = 1 with power concentrated
  26. *    between f1 and f2 Hz for a sampling interval of dt secs.
  27.  
  28. *    Initialization.
  29.  
  30.     ICHANS=1
  31.     GAIN=0
  32.     VOFST=20.0/4096.0
  33.  
  34. *    Get bandpass critical (-3 dB) points.
  35.  
  36.     WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
  37.     READ (*,*) F1
  38.     WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
  39.     READ (*,*) F2
  40.  
  41. *    Get the sampling interval.
  42.  
  43. 10    WRITE (*,'(/1X,A\)') 'Enter delta T secs. (1.0/Sample Rate) : '
  44.     READ (*,*) DELT
  45.       XTST=SQRT(6.0/DELT)*0.5
  46.       ITST=INT4(XTST/VOFST)
  47.       IF (ITST .GT. 32767) THEN
  48.         WRITE (*,'(//1X,A,F7.6,A)')
  49.      +        'This choice of delta t = ',DELT,' ,'
  50.         WRITE (*,'(1X,A,F6.2,A)')
  51.      +        'leads to a maximum data value of ',XTST,' ,'
  52.         WRITE (*,'(1X,A,A,I6,A)')
  53.      +        'which, when converted to an integer,',
  54.      +        ' has a value of ',ITST,' .'
  55.         WRITE (*,'(1X,A/)')
  56.      +        'This cannot be stored in the Integer*2 array NDATA.'
  57.         GO TO 10
  58.       ENDIF
  59.  
  60. *    Get order of filter.
  61.  
  62. 20    WRITE (*,'(/1X,A\)') 'Enter # of filter sections to cascade : '
  63.     READ (*,*) NS
  64.     IF (NS .GT. NSMAX) THEN
  65.       WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
  66.       GO TO 20
  67.     ENDIF
  68.  
  69. *    Get number of points per record.
  70.  
  71. 30    WRITE (*,'(/1X,A\)') 'Enter # of points per record (N) : '
  72.     READ (*,*) N
  73.     IF (N .GT. NMAX) THEN
  74.       WRITE (*,'(1X,A,I5)') 'Reduce the # of points to ',NMAX
  75.       GO TO 30
  76.     ENDIF
  77.  
  78. *    Allocate space for NDATA array.
  79.  
  80.     ALLOCATE (NDATA(N), STAT=IERR)
  81.     IF (IERR .NE. 0)
  82.      +          STOP 'Not enough storage for data.  Aborting ...'
  83.  
  84. *    Get number of records in file.
  85.  
  86.     WRITE (*,'(/1X,A\)') 'Enter # of records (NUMREC) : '
  87.     READ (*,*) NUMREC
  88.  
  89. *    Get seed value for random number generator.
  90.  
  91.     WRITE (*,'(/1X,A\)') 'Enter a negative integer : '
  92.     READ (*,*) IDUM
  93.  
  94. *    Get output file name.
  95.  
  96. 40    WRITE (*,'(/1X,A\)') 'Enter output file name (4 chars) : '
  97.     READ (*,'(A)') FLNM
  98.     WRITE (*,'( )')
  99.  
  100. *    Convert to uppercase and check first character alphabetic.
  101.  
  102.     DO J=4,1,-1
  103.       FIRST=FLNM(J:J)
  104.       IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
  105.         IHOLD=ICHAR(FIRST)-32
  106.         FIRST=CHAR(IHOLD)
  107.         FLNM(J:J)=FIRST
  108.       ENDIF
  109.     ENDDO
  110.     IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
  111.       WRITE (*,'(/1X,A,A,A/1X,A,A,A/1X,A)') 
  112.      +      'Filename ',FLNM,' began with',
  113.      +      'the nonalphabetic character ',FIRST,'.',
  114.      +      'Re-enter the filename correctly.'
  115.       GO TO 40
  116.     ENDIF
  117.  
  118.     FILNAM=FLNM // '.DAT'
  119.     IRSIZE=ICHANS*N*2
  120.     IDELTMS=NINT(DELT*1.0E6)
  121.  
  122. *    Design the bandpass filter.
  123.  
  124.     CALL BPDES (F1,F2,DELT,NS,A,B,C,D,E,GRAF)
  125.  
  126. *    Initialize past values in filter to zero.
  127.  
  128.     X=0.0
  129.  
  130. *    Write the data in the form of binary numbers to a data
  131. *    file that may be read by the spectrum routine.
  132.  
  133.     OPEN (NUMO,FILE=FILNAM,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
  134.      +        FORM='BINARY')
  135.     WRITE (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
  136.     WRITE (NUMO) (GAIN(I),I=0,7)
  137.  
  138. *    Put message on screen.
  139.  
  140.     WRITE (*,'(/////////////////////16X,
  141.      +           ''C O L O R E D   N O I S E   U T I L I T Y'')')
  142.     WRITE (*,'(/25X,''Creating '',A,'' now.''/)') FILNAM 
  143.  
  144. *    Create NUMREC records of the sequence.
  145.  
  146.     DO IB=1,NUMREC
  147.  
  148.       IF (IB .EQ. 1) ITOSS=1000
  149.       IF (IB .NE. 1) ITOSS=0
  150.  
  151. *      Generate the sequence NDATA(J).
  152.  
  153.       DO J=1,N+ITOSS
  154.  
  155.         X(1,5)=SQRT(6.0/DELT)*(RAN1(IDUM)-0.5)
  156.  
  157. *        Go thru NS sections of filter.
  158.  
  159.         DO I=1,NS
  160.           TEMP=A(I)*(X(I,5)-2.0*X(I,3)+X(I,1))-B(I)*X(I+1,4)
  161.           X(I+1,5)=TEMP-C(I)*X(I+1,3)-D(I)*X(I+1,2)-E(I)*X(I+1,1)
  162.         ENDDO
  163.  
  164. *        Push down past values in filter.
  165.  
  166.         DO I=1,NS+1
  167.           DO K=1,4
  168.             X(I,K)=X(I,K+1)
  169.           ENDDO
  170.         ENDDO
  171.  
  172. *        Throw away ITOSS values to eliminate startup transient.
  173.  
  174.         IF (IB .EQ. 1) THEN
  175.           IF (J .LE. ITOSS) CYCLE
  176.           NDATA(J-ITOSS)=INT2(X(NS+1,5)/VOFST)
  177.           CYCLE
  178.         ENDIF
  179.  
  180. *        Set NDATA(J) equal to the output of the filter and continue.
  181.  
  182.         NDATA(J)=INT2(X(NS+1,5)/VOFST)
  183.  
  184.       ENDDO
  185.  
  186. *      Display record number message.
  187.  
  188.       IF (IB .EQ. 1) THEN
  189.         WRITE (*,50) IB
  190. 50        FORMAT (25X,'Writing Record ',I4.4)
  191.       ELSE
  192.         WRITE (*,60) IB
  193. 60        FORMAT ('+',24X,'Writing Record ',I4.4)
  194.       ENDIF
  195.  
  196. *      Create file for time series.
  197.  
  198.       IF (IB .EQ. 1) THEN
  199.         OPEN (3,FILE='COLRTIM.PRN',STATUS='UNKNOWN')
  200.           WRITE (3,'(G15.5,2X,G15.5)')
  201.      +          (FLOAT(J-1)*0.01,FLOAT(NDATA(J))*VOFST,J=1,N)
  202.         CLOSE (3)
  203.       ENDIF
  204.  
  205. *      Save output to a file.
  206.  
  207.       WRITE (NUMO) (NDATA(J),J=1,N)
  208.  
  209.     ENDDO
  210.  
  211.     CLOSE (NUMO,STATUS='KEEP')
  212.  
  213.     WRITE (*,'( )')
  214.     STOP '                    Program terminated successfully.'
  215.     END
  216.  
  217. *    Bandpass filter routine.   (BPDES.FOR)
  218.  
  219. *    Bandpass Butterworth Digital Filter Design Subroutine
  220.  
  221. *    Inputs are passband (-3 dB) frequencies F1 and F2 in Hz,
  222. *            Sampling interval T in seconds, and
  223. *            number NS of filter sections.
  224.  
  225. *    Outputs are NS sets of filter coefficients, i.e.,
  226. *            A(K) thru E(K) for K = 1 thru NS, and
  227. *            20 pairs of frequency and power gain, i.e.,
  228. *            GRAF(1,K) thru GRAF(2,K) FOR K = 1 thru 20.
  229.  
  230. *    Note that A(K) thru E(K) as well as GRAF(2,20) must be
  231. *            dimensioned in the calling program.
  232.  
  233. *    The digital filter has NS sections in cascade.  The Kth
  234. *            section has the transfer function:
  235.  
  236. *                A(K)*(Z**4-2*Z**2+1)
  237. *        H(Z) =  ------------------------------------
  238. *            Z**4+B(K)*Z**3+C(K)*Z**2+D(K)*Z+E(K)
  239.  
  240. *    Thus, if F(M) and G(M) are the input and output of the 
  241. *            Kth section at time M*T, then
  242.  
  243. *        G(M)=A(K)*(F(M)-2*F(M-2)+F(M-4))-B(K)*G(M-1)
  244. *            -C(K)*G(M-2)-D(K)*G(M-3)-E(K)*G(M-4)
  245.  
  246. *    Do not use this routine for lowpass design.  Use LPDES.
  247. *    Do not use this routine for highpass design.  Use HPDES.
  248.  
  249.     SUBROUTINE BPDES (F1,F2,T,NS,A,B,C,D,E,GRAF)
  250.  
  251.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  252.     REAL*4         A(*),B(*),C(*),D(*),E(*),GRAF(2,20)
  253.  
  254. *    PI=2.0*ASIN(1.0)
  255.     PI=3.1415926536
  256.  
  257.     W1=SIN(F1*PI*T)/COS(F1*PI*T)
  258.     W2=SIN(F2*PI*T)/COS(F2*PI*T)
  259.     WC=W2-W1
  260.     Q=WC*WC+2.0*W1*W2
  261.     S=W1*W1*W2*W2
  262.  
  263.     DO 150 K=1,NS
  264.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  265.       P=-2.0*WC*CS
  266.       R=P*W1*W2
  267.       X=1.0+P+Q+R+S
  268.       A(K)=WC*WC/X
  269.       B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
  270.       C(K)=(6.0-2.0*Q+6.0*S)/X
  271.       D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
  272.       E(K)=(1.0-P+Q-R+S)/X
  273. 150    CONTINUE
  274.  
  275.     DO 170 J=1,2
  276.       DO 160 I=1,10
  277.         K=I*(2-J)+(21-I)*(J-1)
  278.         GRAF(2,K)=0.01+0.98*FLOAT(I-1)/9.0
  279.         X=(1.0/GRAF(2,K)-1.0)**(1.0/FLOAT(4*NS))
  280.         SQ=SQRT(WC*WC*X*X+4.0*W1*W2)
  281.         GRAF(1,K)=ABS(ATAN(0.5*(WC*X+FLOAT(2*J-3)*SQ)))/(PI*T)
  282. 160      CONTINUE
  283. 170    CONTINUE
  284.  
  285.     RETURN
  286.     END
  287.  
  288.     FUNCTION RAN1(IDUM)
  289.  
  290. *    Returns a uniform random deviate between 0.0 and 1.0.  Set
  291. *    IDUM to any negative value to initialize or reinitialize the
  292. *    sequence. Taken from Numerical Recipes, p.196-197.
  293.  
  294.     INTEGER*2    IDUM
  295.     REAL*4        R(97)
  296.     PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
  297.     PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
  298.     PARAMETER (M3=243000,IA3=4561,IC3=51349)
  299.     DATA        IFF /0/
  300.  
  301.     IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
  302.       IFF=1
  303.       IX1=MOD(IC1-IDUM,M1)
  304.       IX1=MOD(IA1*IX1+IC1,M1)
  305.       IX2=MOD(IX1,M2)
  306.       IX1=MOD(IA1*IX1+IC1,M1)
  307.       IX3=MOD(IX1,M3)
  308.       DO J=1,97
  309.         IX1=MOD(IA1*IX1+IC1,M1)
  310.         IX2=MOD(IA2*IX2+IC2,M2)
  311.         R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
  312.       ENDDO
  313.       IDUM=1
  314.     ENDIF
  315.  
  316.     IX1=MOD(IA1*IX1+IC1,M1)
  317.     IX2=MOD(IA2*IX2+IC2,M2)
  318.     IX3=MOD(IA3*IX3+IC3,M3)
  319.     J=1+(97*IX3)/M3
  320.     IF(J.GT.97.OR.J.LT.1)PAUSE
  321.     RAN1=R(J)
  322.     R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
  323.  
  324.     RETURN
  325.     END
  326.