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===========================================================================
22 c nonlinear noise reduction
23 c see H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
24 c University Press (1997,2004)
25 c authors T. Schreiber, H. Kantz, R. Hegger (1998) based on earlier versions
26 c===========================================================================
28 dimension x(nx), x0(nx), xc(nx)
29 character*72 file, fout
33 call whatido("nonlinear noise reduction (see also: noise)",iverb)
42 isout=igetout(fout,iverb)
44 call nthstring(1,file)
45 call readfile(nmax,x,nexcl,jcol,file,iverb)
46 if(file.eq."-") file="stdin"
47 if(isout.eq.1) call addsuff(fout,file,"_")
52 call clean(nmax,x,xc,m,kmin,nq,eps,iverb)
53 if(fout.ne." ".or.isout.eq.1.or.it.eq.imax) then
54 if(isout.eq.1) call suffix(fout,"c")
55 call outfile(fout,iunit,iverb)
57 30 write(iunit,*) xc(n), x0(n)-xc(n)
58 if(iunit.ne.istdout()) close(iunit)
59 if(iv_io(iverb).eq.1) call writereport(nmax,fout)
63 eps=eps+(xc(n)-x(n))**2
66 20 if(iv_io(iverb).eq.1)
67 . write(istderr(),*) 'New diameter of neighbourhoods is ', eps
74 . "-m# -q# -r# -k# [-i# -o outfile -l# -x# -c# -V# -h] file")
75 call popt("m","embedding dimension")
76 call popt("q","dimension of manifold")
77 call popt("r","radius of neighbourhoods")
78 call popt("k","minimal number of neighbours")
79 call popt("i","number of iterations (1)")
80 call popt("l","number of values to be read (all)")
81 call popt("x","number of values to be skipped (0)")
82 call popt("c","column to be read (1 or file,#)")
83 call pout("file_c, file_cc (etc.)")
85 call ptext("Verbosity levels (add what you want):")
86 call ptext(" 1 = input/output" )
87 call ptext(" 2 = state of neighbour search")
92 subroutine clean(nmax,y,yc,m,kmin,nq,d,iverb)
93 parameter(im=100,ii=100000000,nx=100000,mm=15,small=0.0001)
94 dimension y(nmax),yc(nmax),r(mm),ju(nx),c(mm,mm),cm(mm),
95 . jh(0:im*im),jpntr(nx),nlist(nx), zcm(mm,nx)
97 if(nmax.gt.nx.or.m.gt.mm) stop "clean: make mm/nx larger."
98 sr=2*small+m-2 ! ${\rm tr}(1/r)=1$
101 10 if(i.eq.m.or.i.eq.1) r(i)=sr/small
109 1 call base(nmax,y,1,m,jh,jpntr,eps)
111 do 50 nn=1,iu ! find neighbours
113 call neigh(nmax,y,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
114 if(nfound.lt.kmin) then ! not enough neighbours found
115 iunp=iunp+1 ! mark for next sweep
117 else ! fine: enough neighbours
118 do 90 i=1,m ! centre of mass vector
121 100 s=s+y(nlist(np)-m+i)
123 if(istep.eq.1) then ! just store centre of mass
127 do 120 i=1,m ! corrected centre of mass vector
130 130 s=s+zcm(i,nlist(np))
131 120 cm(i)=2*cm(i)-s/nfound
132 do 140 i=1,m ! compute covariance matrix
137 150 s=s+(y(jm+i)-cm(i))*(y(jm+j)-cm(j))
138 c(i,j)=r(i)*r(j)*s/nfound
140 call eigen(c,m) ! find eigenvectors (decreasing)
145 170 s=s+(y(n-m+j)-cm(j))*c(i,iq)*c(j,iq)*r(j)
146 160 yc(n-m+i)=yc(n-m+i)-s/r(i)/r(i)
151 if(iv_uncorr(iverb).eq.1)
152 . write(istderr(),*) "With ", eps, iunp, " uncorrected"
154 30 if(iunp.ne.0) goto 1
157 c driver for diagonalisation routines
158 c see H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
159 c University Press (1997)
160 c Copyright (C) T. Schreiber (1997)
162 subroutine eigen(c,kk)
164 dimension c(md,md),d(md),w1(md),w2(md),z(md,md)
165 if(kk.gt.md) stop "eigen: make md larger."
167 call rs(md,kk,c,d,1,z,w1,w2,ierr)
170 10 c(i,j)=z(i,kk+1-j)