program noise ! this program produces white noise with a selected rms value. ! the noise is very close to gaussian. implicit real*8 (a - h, o - z) character*40 output write (*,'(5x, a)') 'enter filespec of output file.' read (*, '(a)') output open (unit = 1, file = output, status = 'unknown') ! input starting words for pseudo-random sequence (2 arbitrarily chosen ! integers); number of noise points to calculate; rms value of generated ! noise (it will approach this for a large number of generated points). write (*,'(5x, a)') 'enter idum1, idum2, npoints, rms' read (*, *) idum1, idum2, npoints, rms idum1 = - abs (idum1) idum2 = - abs (idum2) sumnoise = 0. do i = 1, npoints call gauss0 (idum1, idum2, rms, xnoise) write (1, '(i6, 1pe13.5)') i, xnoise sumnoise = sumnoise + xnoise**2 end do write (*, *) 'noise rms = ', sqrt (sumnoise / npoints) close (unit = 1) stop end ! subroutine gauss0 (idum1, idum2, sdev, x) ! generates normally distributed random number, zero-mean ! idum1 (input/output) must contain a negative number on the first entry ! idum2 (input/output) must contain a negative number on the first entry ! sdev (input) standard deviation of desired distribution ! x (output) the value of the computed normal variable ! box-muller method (polar algorithm): ! the algorithm is taken from donald e. knuth (standford university) ! "the art of computer programming", volume 2 / seminumerical algorithms, ! second edition, addison-wesley, isbn 0-201-03822-6, page 117 ! 09-jun-88: first implementation implicit real*8 (a - h, o - z) save 10 v1 = 2.0 * (ran3 (idum1)) - 1.0 v2 = 2.0 * (ran3 (idum2)) - 1.0 s = v1 * v1 + v2 * v2 if (.not. (s .lt. 1.0)) go to 10 ! select one of the two possible solutions if (ran3 (idum1) .gt. 0.5) then x = sdev * v1 * dsqrt (- 2.0 * dlog (s) / s) else x = sdev * v2 * dsqrt (- 2.0 * dlog (s) / s) end if return end ! function ran3 (idum) implicit real*8 (a - h, o - z) parameter (mbig = 1000000000, mseed = 161803398) parameter (mz = 0, fac = 1. / mbig) dimension ma (55) data iff /0/ save if (idum .lt. 0 .or. iff .eq. 0) then iff = 1 mj = mseed - iabs (idum) mj = mod (mj, mbig) ma (55) = mj mk = 1 do i = 1, 54 ii = mod (21 * i, 55) ma (ii) = mk mk = mj - mk if (mk .lt. mz) mk = mk + mbig mj = ma (ii) end do do k = 1, 4 do i = 1, 55 ma (i) = ma (i) - ma (1 + mod (i + 30, 55)) if (ma (i) .lt. mz) ma (i) = ma (i) + mbig end do end do inext = 0 inextp = 31 idum = 1 end if inext = inext + 1 if (inext .eq. 56) inext = 1 inextp = inextp + 1 if (inextp .eq. 56) inextp = 1 mj = ma (inext) - ma (inextp) if (mj .lt. mz) mj = mj + mbig ma (inext) = mj ran3 = mj * fac return end