1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
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.
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.
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
21 c===========================================================================
23 c simple nonlinear noise reduction
24 c see H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
25 c University Press (1997,2004)
26 c author T. Schreiber (1998)
27 c===========================================================================
29 dimension x(nx), x0(nx), xc(nx)
30 character*72 file, fout
31 data eps/0./, frac/0./, imax/1/
34 call whatido("simple nonlinear noise reduction",iverb)
42 isout=igetout(fout,iverb)
43 if(eps.eq.0.and.frac.eq.0.) call usage()
45 call nthstring(1,file)
46 call readfile(nmax,x,nexcl,jcol,file,iverb)
47 if(file.eq."-") file="stdin"
48 if(isout.eq.1) call addsuff(fout,file,"_l")
49 call rms(nmax,x,sc,sd)
50 if(frac.gt.0) eps=sd*frac
54 call nrlazy(nmax,x,xc,m,eps)
55 if(fout.ne." ".or.isout.eq.1.or.it.eq.imax) then
56 if(isout.eq.1) call suffix(fout,"c")
57 call outfile(fout,iunit,iverb)
59 30 write(iunit,*) xc(n), x0(n)-xc(n)
60 if(iunit.ne.istdout()) close(iunit)
61 if(iv_io(iverb).eq.1) call writereport(nmax,fout)
65 eps=eps+(xc(n)-x(n))**2
69 if(iv_io(iverb).eq.1) write(istderr(),*)
70 . 'Zero correction, finished'
73 20 if(iv_io(iverb).eq.1) write(istderr(),*)
74 . 'New diameter of neighbourhoods is ', eps
81 . "-m# [-r# | -v#] [-i# -o outfile -l# -x# -c# -V# -h] file")
82 call ptext("either -r or -v must be present")
83 call popt("m","embedding dimension")
84 call popt("r","absolut radius of neighbourhoods")
85 call popt("v","same as fraction of standard deviation")
86 call popt("i","number of iterations (1)")
87 call popt("l","number of values to be read (all)")
88 call popt("x","number of values to be skipped (0)")
89 call popt("c","column to be read (1 or file,#)")
90 call pout("file_lc, file_lcc (etc.)")
95 subroutine nrlazy(nmax,y,yc,m,eps)
96 parameter(im=100,ii=100000000,nx=1000000)
97 dimension y(nmax),yc(nmax),jh(0:im*im),jpntr(nx),nlist(nx)
99 if(nmax.gt.nx) stop "nrlazy: make nx larger."
100 call base(nmax,y,1,m,jh,jpntr,eps)
104 call neigh(nmax,y,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
107 30 av=av+y(nlist(nn)-(m-1)/2) ! average middle coordinate
108 20 yc(n-(m-1)/2)=av/nfound