website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/addnoise.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   addnoise.f
  23 c
  24 c   add Gaussian / uniform white noise
  25 c   author T. Schreiber (1998)
  26 c===========================================================================
  27
  28       parameter(nx=1000000)
  29       character*72 file, fout
  30       dimension x(nx)
  31       external rand
  32       data eps/0./, frac/0./, iuni/0/
  33       data iverb/1/
  34
  35       call whatido("add Gaussian/uniform noise",iverb)
  36       eps=fcan("r",eps)
  37       frac=fcan("v",frac)
  38       nmaxx=ican("l",nx)
  39       nexcl=ican("x",0)
  40       jcol=ican("c",0)
  41       if(lopt("u",1).eq.1) iuni=1
  42       isout=igetout(fout,iverb)
  43       if(eps.eq.0.and.frac.eq.0.) call usage()
  44
  45       if(lopt("0",1).eq.1.and.eps.gt.0) then
  46          if(isout.eq.1) fout="0_noisy"
  47          do 10 n=1,nmaxx
  48             if(iuni.eq.1) then
  49                x(n)=rand(0.0)*eps
  50             else
  51                x(n)=rgauss(0.0,eps)
  52             endif
  53  10         continue
  54          call writefile(nmaxx,x,fout,iverb)
  55          stop
  56       endif
  57
  58       do 20 ifi=1,nstrings()
  59          call nthstring(ifi,file)
  60          nmax=nmaxx
  61          call readfile(nmax,x,nexcl,jcol,file,iverb)
  62          if(file.eq."-") file="stdin"
  63          call rms(nmax,x,sc,sd)
  64          if(frac.gt.0) eps=sd*frac
  65          if(iuni.eq.1) then
  66             if(iv_io(iverb).eq.1) write(istderr(),*)
  67      .      "adding uniform noise in [0,", eps,"]"
  68          else
  69             if(iv_io(iverb).eq.1) write(istderr(),*)
  70      .      "adding Gaussian noise of amplitude", eps
  71          endif
  72          if(sd.gt.0.and.iv_io(iverb).eq.1) write(istderr(),*)
  73      .      "that is",eps/sd,"* rms of data"
  74          do 30 n=1,nmax
  75             if(iuni.eq.1) then
  76                x(n)=x(n)+rand(0.0)*eps
  77             else
  78                x(n)=x(n)+rgauss(0.0,eps)
  79             endif
  80  30         continue
  81          if(isout.eq.1) call addsuff(fout,file,"_noisy")
  82  20      call writefile(nmax,x,fout,iverb)
  83       end
  84
  85       subroutine usage()
  86 c usage message
  87
  88       call whatineed(
  89      .   "[-r# | -v#] [-u -0 -o outfile -l# -x# -c# -V# -h] file(s)")
  90       call ptext("either -r or -v must be present")
  91       call popt("r","absolute noise level")
  92       call popt("v","same as fraction of standard deviation")
  93       call popt("u","add uniform noise (default Gaussian)")
  94       call popt("0","do not read input, just issue random numbers")
  95       call popt("l","number of values to be read (all)")
  96       call popt("x","number of values to be skipped (0)")
  97       call popt("c","column to be read (1 or file,#)")
  98       call pout("file_noisy")
  99       call pall()
 100       stop
 101       end
 102