Commit b3919285 authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

new feature

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 2133
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent f62d112d
c this is <gabor.f>
c------------------------------------------------------------------------------
c ($Id: gabor.f,v 1.10 2005-01-10 19:31:36 tforb Exp $)
c ($Id: gabor.f,v 1.11 2006-06-13 11:21:10 tforb Exp $)
c
c 22/02/2001 by Thomas Forbriger (IMGF Frankfurt)
c
......@@ -10,15 +10,17 @@ c REVISIONS and CHANGES
c 22/02/2001 V1.0 Thomas Forbriger
c 02/08/2001 V1.1 allow trace order rather than offset order
c 09/09/2004 V1.2 set lower limit of time window
c 13/06/2006 V1.3 provide Hanning taper
c
c==============================================================================
c
program gabor
c
character*(*) version
parameter(version='GABOR V1.2 calculate gabor matrix')
parameter(version='GABOR V1.3 calculate gabor matrix')
character*(*) GABOR_CVS_ID
parameter(GABOR_CVS_ID='$Id: gabor.f,v 1.10 2005-01-10 19:31:36 tforb Exp $')
parameter(GABOR_CVS_ID=
& '$Id: gabor.f,v 1.11 2006-06-13 11:21:10 tforb Exp $')
c
c input dataset
character*80 filename
......@@ -26,9 +28,9 @@ c input dataset
parameter(maxtraces=40, totmaxsamples=4000000)
integer lu
parameter(lu=12)
real data(totmaxsamples)
real fdata(totmaxsamples)
integer idata(totmaxsamples)
equivalence(data,idata)
equivalence(fdata,idata)
real toffset(maxtraces), tracedt(maxtraces), roffset(maxtraces)
integer innsamples(maxtraces), firstsample(maxtraces)
integer ntraces
......@@ -52,31 +54,32 @@ c greens function
real tref(maxtime), om(maxom)
c
c processing
real fmax,tau,tmax,winfac,stdfac,tmin
real fmax,tau,tmax,winfac,stdfac,tmin, factor, now, wlen
integer ntime, nptraces, nom
character*80 outbase, outfile
logical overwrite, traceorder
logical overwrite, traceorder, hanning
parameter(stdfac=8.)
c
integer i,j,k,itrace,itime
c constants
complex ime
real pi2,hin
real pi2,hin,pi
parameter(hin=1.,pi2=2.*3.14159265358979,ime=(0.,1.))
parameter(pi=3.14159265358979)
c debugging
logical debug, verbose
c commandline
integer maxopt, lastarg, iargc
character*80 argument
parameter(maxopt=10)
parameter(maxopt=11)
character*2 optid(maxopt)
character*40 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
data optid/2h-d, 2h-v, '-f', '-n', '-o', '-t', '-T', '-w', '-S',
& 2h-m/
data opthasarg/2*.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE./
data optarg/2*1h-,'100.','100','-','1.','1','1.','1','0.'/
& 2h-m, '-h'/
data opthasarg/2*.FALSE.,2*.TRUE.,.FALSE.,5*.TRUE.,.false./
data optarg/2*1h-,'100.','100','-','1.','1','1.','1','0.','-'/
c
c------------------------------------------------------------------------------
c basic information
......@@ -87,10 +90,12 @@ c
if ((argument(1:5).eq.'-help').or.(iargc().lt.2)) then
print *,version
print *,'Usage: gabor datafile outfile [-v] [-d]'
print *,' [-f freq] [-n N] [-t tmax] [-o] [-T N | -S N]'
print *,' [-f freq] [-n N] [-t tmax] [-o]'
print *,' [-T N | -S N] [-h]'
print *,' [-w fac] [-m tmin]'
print *,' or: gabor -help'
if (argument(1:5).ne.'-help') stop 'ERROR: wrong number of arguments'
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
print *,'calculate gabor matrix'
print *,' '
......@@ -103,11 +108,13 @@ c
print *,'-f freq analyse frequencies up to ''freq'' '
print *,'-n N use ''N'' time steps'
print *,'-t tmax use windows up to time ''tmax'' '
print *,'-o overwrite mode when writing results to file'
print *,'-o '
print *, 'overwrite mode when writing results to file'
print *,'-T N use first ''N'' traces in offset order'
print *,'-S N use first ''N'' traces in trace order'
print *,'-w fac broaden time windoe by factor ''fac'' '
print *,'-m tmin use windows starting at tmin'
print *,'-h use Hanning taper rather than Gaussian'
print *,' '
print *,GABOR_CVS_ID
stop
......@@ -129,6 +136,7 @@ c
traceorder=optset(9)
if (optset(9)) read(optarg(9), *) nptraces
read(optarg(10), *) tmin
hanning=optset(11)
overwrite=optset(5)
......@@ -138,8 +146,9 @@ c
call getarg(2, outbase)
c------------------------------------------------------------------------------
c go
call sffu_simpleread(lu, filename, maxtraces, totmaxsamples, data, idata,
& toffset, tracedt, roffset, innsamples, firstsample, ntraces, verbose)
call sffu_simpleread(lu, filename, maxtraces, totmaxsamples,
& fdata, idata, toffset, tracedt, roffset, innsamples,
& firstsample, ntraces, verbose)
c
call tf_rchain(roffset, chain, ntraces, first, 1)
c
......@@ -170,7 +179,6 @@ c
nsamples=nsamples*2
enddo
if (nsamples.gt.maxsamples) stop 'ERROR: too many samples'
if (verbose) print *,'going to use ',nsamples,' samples'
df=1./(float(nsamples)*dt)
nfny=nsamples/2
c
......@@ -178,7 +186,10 @@ c
c----------------------------------------------------------------------
tmax=min(tmax,float(nsamples)*dt)
k=first
c effective window length
tau=stdfac*winfac*tmax/float(ntime)
wlen=2.*tau*sqrt(-log(0.5))
if (hanning) tau=wlen
if (debug) print *,'tau: ',tau, ' tmax: ',tmax, ' ntime: ',ntime
nom=fmax/df
do i=1,nom
......@@ -187,16 +198,52 @@ c----------------------------------------------------------------------
do itime=1,ntime
tref(itime)=tmin+float(itime-1)*(tmax-tmin)/float(ntime-1)
enddo
if (verbose) then
print *,'analysis parameters:'
print *,' number of samples used per trace: ',nsamples
print *,' number of timesteps to output: ',ntime
print *,' number of frequencies to output: ',nom
print *,' Nyquist frequency: ',nfny*df,' Hz'
print *,' smallest frequency to output: ',df,' Hz'
print *,' maximum frequency to output: ',fmax,' Hz'
print *,' first sample taken at: ',tmin,' s'
print *,' last sample taken at: ',tmax,' s'
print *,' window halfwidth: ',wlen,' s'
if (hanning) then
print *,' use Hanning taper'
else
print *,' use Gaussian taper'
print *,' Gaussian time constant: ',tau,' s'
endif
endif
do itrace=1,nptraces
if (traceorder) k=itrace
c----------------------------------------------------------------------
c inner loop, calculates Gabor matrix
do itime=1,ntime
do i=1,nsamples
process(i)=(0.,0.)
enddo
c copy tapered time series
do i=1,innsamples(k)
process(i)=cmplx(data(i+firstsample(k)-1))*
& exp(-1.*((dt*(i-1)-tref(itime))/tau)**2)
now=dt*(i-1)
if (hanning) then
if (now.lt.(tref(itime)-tau)) then
factor=0.
c elseif (now.lt.tref(itime)) then
c factor=0.5*(1.+cos(pi*(tref(itime)-now)/tau))
elseif (now.lt.(tref(itime)+tau)) then
factor=0.5*(1.+cos(pi*(tref(itime)-now)/tau))
else
factor=0.
endif
if (debug.and.(itime.eq.20).and.(itrace.eq.1)) then
print *,'DEBUG ',i,' factor= ',factor
endif
else
factor=exp(-1.*((now-tref(itime))/tau)**2)
endif
process(i)=cmplx(fdata(i+firstsample(k)-1))*factor
enddo
c calculate Fourier coefficients
call tf_fork(nsamples,process,hin)
......@@ -204,6 +251,8 @@ c store them
do i=1,nom
gabormat(itime,i)=process(i)
enddo
c end of inner loop, this trace's Gabor matrix is ready
c----------------------------------------------------------------------
enddo
k=chain(k)
c
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment