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 power spectrum of spike trains
23 c author T. Schreiber (1999)
24 c===========================================================================
26 dimension x(nx), lx(nx), sp(nx)
27 character*72 file, fout
30 call whatido("power spectrum of spike trains",iverb)
35 nfreq=min(ican("#",0),nx)
38 isout=igetout(fout,iverb)
40 do 10 ifi=1,nstrings()
41 call nthstring(ifi,file)
43 call readfile(nmax,x,nexcl,jcol,file,iverb)
47 1 call sort(nmax,x,lx)
48 if(fmax.le.0.) fmax=2*nmax/(x(nmax)-x(1))
49 if(nfreq.le.0) nfreq=fmax*(x(nmax)-x(1))/2
50 write(istderr(),*) "spikespec: total time covered: ",
52 write(istderr(),*) "spikespec: computing ", nfreq,
53 . " frequencies up to ", fmax
56 30 sp(n)=sspect(nmax,x,f)
58 if(ibin.gt.0) write(istderr(),*)
59 . "spikespec: binning", 2*ibin+1," frequencies"
60 if(file.eq."-") file="stdin"
61 if(isout.eq.1) call addsuff(fout,file,"_ss")
62 call outfile(fout,iunit,iverb)
63 do 40 n=1+ibin,nfreq-ibin,2*ibin+1
66 do 50 ib=n-ibin,n+ibin
68 40 write(iunit,*) f, p
69 if(iunit.eq.istdout()) write(iunit,*)
70 if(iunit.eq.istdout()) write(iunit,*)
71 10 if(iunit.ne.istdout()) close(iunit)
74 function sspect(nmax,x,f)
83 10 si=si+sin(omega*x(n))
91 . "[-F# -## -w# -i -o outfile -l# -x# -c# -V# -h] file(s)")
92 call popt("F","maximal frequency [2*l / total time]")
93 call popt("#","number of frequencies [F* total time /2]")
94 call popt("w","frequency resolution [0]")
95 call popt("i","input data: intervals rather than times")
96 call popt("l","number of values to be read [all]")
97 call popt("x","number of values to be skipped [0]")
98 call popt("c","column to be read [1 or file,#]")