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

lisousi no longer lives here

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: 5132
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent e22cc09f
......@@ -26,13 +26,14 @@
# siggen and siggenx (being linked against libsff and
# libfapidxx, respectively)
# 11/01/2013 V1.8 added sigscale to program targets
# 17/04/2013 V1.9 lisousi has its own place now (moved source code)
#
# ============================================================================
PROGRAMS=tsfilt stufi rotate coro xyz2uvw susei evelo tesiff teswf \
phasedsignals hamres siggen smoos dise foutra teseco resaseda gatherdiff \
autocorr cross tidofi fredofi sigfit noisymize sigval fidasexx soutifu \
deconv geophone siggenx lisousi sigscale
deconv geophone siggenx sigscale
.PHONY: all
all: install doc
......@@ -142,18 +143,6 @@ siggenx: %x: %.o
-ldatrwxx -lsffxx \
-lgsexx -ltime++ -laff -lgsl -lm -lgslcblas
lisousi.cc: lisousi_description_text.h
lisousi.cc: lisousi_help_text.h
lisousi: \
%: %_description_text.o %_help_text.o %.o \
%_filterresponse.o \
%_applytaper.o
$(CXX) -o $@ $^ -I$(LOCINCLUDEDIR) -lfourierxx -lfftw3 -lm \
-lsffxx -ldatrwxx -llinearxx -lgsl -lgslcblas \
-ltsxx -ltfxx -lsffxx -lgsexx -ltime++ -laff \
-lsffxx $(FOTRANLIB) -lm\
-L$(LOCLIBDIR) $(CXXFLAGS) $(FLAGS) $(LDFLAGS)
deconv foutra sigval: \
%: %.o
$(CXX) -o $@ $^ -I$(LOCINCLUDEDIR) -lfourierxx -lfftw3 -lm \
......
This diff is collapsed.
/*! \file lisousi.h
* \brief prototypes and structs for lisousi (prototypes)
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 17/04/2013
*
* prototypes and structs for lisousi (prototypes)
*
* Copyright (c) 2013 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ----
*
*
* REVISIONS and CHANGES
* - 17/04/2013 V1.0 Thomas Forbriger
*
* ============================================================================
*/
// include guard
#ifndef TF_LISOUSI_H_VERSION
#define TF_LISOUSI_H_VERSION \
"TF_LISOUSI_H V1.0 "
#define TF_LISOUSI_H_CVSID \
"$Id$"
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <gsl/gsl_sf_bessel.h>
#include <tfxx/commandline.h>
#include <tfxx/xcmdline.h>
#include <tfxx/error.h>
#include <tfxx/stringfunc.h>
#include <tfxx/rangestring.h>
#include <tfxx/rangelist.h>
#include <tfxx/misc.h>
#include <tsxx/anyfilter.h>
#include <tsxx/convolve.h>
#include <datrwxx/readany.h>
#include <datrwxx/writeany.h>
#include <sffxx.h>
#include <sffostream.h>
#include <aff/dump.h>
#include <aff/seriesoperators.h>
#include <aff/subarray.h>
#include <tsxx/ovtaper.h>
#include <fourier/fftwaff.h>
using std::cout;
using std::cerr;
using std::endl;
typedef ts::filter::Ttimeseries Ttimeseries;
typedef Ttimeseries::Tseries Tseries;
/*----------------------------------------------------------------------*/
// type of Fourier domain solution
enum Efdtype {
Ffdfarfield,
Ffdexplosion,
Ffdzforce,
Ffdlamb
};
struct Options {
bool verbose, debug;
bool overwrite, taperfirst, sqrttaper, limitlength;
bool fredomain, fdfilter, tdfilter, nointeg, transition, tapsloset;
double tfac, tshift, tlim, integshift, tapdel;
double tapslo, transition1, transition2;
int npad;
std::string inputformat, outputformat;
// single-velocity switches:
Efdtype fdtype;
// propagation parameters:
double velocity, vpvsratio, pquality, squality;
// trapezoid integration:
double tzedgefactor, tzedgetaper;
int tznsteps;
}; // struct Options
/*======================================================================*/
/* functions */
#endif // TF_LISOUSI_H_VERSION (includeguard)
Tseries applytaper(const Tseries::Tcoc& input, const Tseries::Tcoc& taper,
const double& factor, const double& dt,
const double& offset, const Options& opt);
Tseries filterresponse(const int& size, const double& dt, const Options& opt);
/* ----- END OF lisousi.h ----- */
/*! \file lisousi_applytaper.cc
* \brief apply a taper to the time series (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 17/04/2013
*
* apply a taper to the time series (implementation)
*
* Copyright (c) 2013 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ----
*
*
* REVISIONS and CHANGES
* - 17/04/2013 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define TF_LISOUSI_APPLYTAPER_CC_VERSION \
"TF_LISOUSI_APPLYTAPER_CC V1.0 "
#define TF_LISOUSI_APPLYTAPER_CC_CVSID \
"$Id$"
#include "lisousi.h"
Tseries applytaper(const Tseries::Tcoc& input, const Tseries::Tcoc& taper,
const double& factor, const double& dt,
const double& offset, const Options& opt)
{
if (opt.verbose)
{
cout << " apply taper function"
<< " with delay: " << dt*int(opt.tapdel/dt) << "s"
<< endl;
}
// calculate taper delay
double taperdelay=opt.tapdel;
if (opt.tapsloset)
{
taperdelay=taperdelay<opt.tapslo*offset ? taperdelay : opt.tapslo*offset;
}
// prepare return value
Tseries retval(0,input.size()-1);
retval=0.;
// prepare iterators
aff::Browser<Tseries::Tcoc> IT(taper);
aff::Browser<Tseries::Tcoc> II(input);
aff::Iterator<Tseries> IS(retval);
unsigned int i=0;
// apply taper
while (IT.valid() && IS.valid() && II.valid())
{
double t=i*dt;
(*IS) = (*II) * (*IT) * factor;
++IS;
++II;
if (t >= taperdelay) { ++IT; }
++i;
}
return(retval);
}
/* ----- END OF lisousi_applytaper.cc ----- */
# this is <lisousi_description_text.txt>
# ============================================================================
# description of lisousi processing concepts
# ------------------------------------------
# $Id$
# ============================================================================
#
$Id$
This program takes seismograms from the input file which are
expected to be vertical or radial components of the surface
displacement (or particle velocity) excited by a point source
(vertical single force or explosion) at the surface. The program
is able to apply transformation procedures to the data which
simulate equivalent waveforms as being excited by a line source
perpendicular to the direction of wave propagation. Different
approaches are offered. All of them operate on single traces
and just require the offset of the trace as additional parameter
(in contrast to a Fourier-Bessel transformation, which requires
a complete common source gather).
All approaches are based on the ratio of the 2D and the 3D Greens
function in the frequency domain. This ratio can be understood as
a convolution with 1/sqrt(t) in the time domain and appropriate
scaling (tapering) of the samples in the time domain. The appropriate
scaling factor depends of wave travel time and wave travel distance.
In the case of this program, we assume direct waves, such that the
offset equals the travel distance. The sample time is taken as
travel time. This way the waveforms must be tapered by 1/sqrt(t)
and scaled by the offset times sqrt(2).
This differs from approaches used in reflection seismic analysis,
where waveforms are tapered by sqrt(t). The latter is offered as an
option.
The program offers options to select different ways of constructing
the 1/sqrt(t) filter and the way this filter is applied. Application
of the filter (convolution with 1/sqrt(t)) and the application of
the scaling (tapering with 1/sqrt(t)) are not commutative in a
mathematical sense and the results may change for a different
order of the operations.
Available approaches:
---------------------
Default:
Convolution with 1/sqrt(t) and tapering with 1/sqrt(t).
So-called "direct wave transformation"
Application of filter:
1) Frequency domain application:
recommended approach
default; fastest approach
Convolution with the 1/sqrt(t) filter response is done in the
Fourier domain by multiplication with the analytically derived
Fourier coefficients of the filter response function.
2) Frequency domain application of dicrete Fourier transform:
select option: -fdfilter
The 1/sqrt(t) filter response is constructed in the time
domain. Convolution with the time series is done in the
frequency domain after FFT.
3) Time domain application:
slowest approach
select option: -tdfilter
The 1/sqrt(t) filter response is constructed in the time
domain. Convolution with the time series is done in the
time domain by discrete convolution.
Options:
For the frequency domain applications (1 and 2) the following
options are available:
-pad n padding of time series
The 1/sqrt(t) filter response becomes singular for t=0. This
requires special measures for the approaches which construct
the filter response in the time domain (2 and 3) in order to
create a response which behaves equivalent to its analytical
counterpart. The following options are available:
-integshift t
-notinteg
If -nointeg is selected, the following options are available:
-tshift f
-tlim f
-tfac f
Application of the taper:
By default a 1/sqrt(t) taper is used.
Option -sqrttaper selects sqrt(t) as taper function and scales
seismograms additionally by a factor sqrt(2)*velocity.
Tapers are applied prior to filtering by default. Option -taperlast
selects to taper the time series after filtering.
Tapers are always constructed in the time domain. The default
1/sqrt(t) taper is singular at t=0. This requires special
measures to set up the taper appropriately, such that it behaves
equivalent to its analytical counterpart. The following options
are available:
-integshift t
-notinteg
If -nointeg is selected, the following options are available:
-tshift f
-tlim f
-tfac f
They have a lesser (or even unnoticeable) impact when compared
to their impact whne constructing the filter resonse. This is
because the time series to be tapered usually are small or
vanish close to t=0.
Sampling time is taken as travel time in order to determine the
appropriate scaling (taper) factor. Since a finite bandwidth
wavelet appears slightly after the ray-theoretical arrival
in the seismogram, scaling factors are typically too small
for near offset traces, where pulse duration is close to
pulse travel time. The program offers the option -tapdel to
apply a correction for this effect. -tapslo can be used to
set an upper limit for tapdel at small offsets.
Single velocity transformation:
select by option: -fredomain
additional options:
-velocity v
-pad n
Operation:
A Fourier transform of the waveform is divided by the Fourier
transform of the 3D Greens function and multiplied by the 2D
Greens function. This approach works only for waves of a single
given wave velocity. It is hence only appropriate for the
impulse response in homogeneous full space, but also performs
surprisingly well when applied to dispersive wavefields.
# ----- END OF lisousi_description_text.txt -----
/*! \file lisousi_filterresponse.cc
* \brief create a 1/sqrt(t) filter response for given parameters (implementation)
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 17/04/2013
*
* create a 1/sqrt(t) filter response for given parameters (implementation)
*
* Copyright (c) 2013 by Thomas Forbriger (BFO Schiltach)
*
* ----
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
* ----
*
*
* REVISIONS and CHANGES
* - 17/04/2013 V1.0 Thomas Forbriger
*
* ============================================================================
*/
#define TF_LISOUSI_FILTERRESPONSE_CC_VERSION \
"TF_LISOUSI_FILTERRESPONSE_CC V1.0 "
#define TF_LISOUSI_FILTERRESPONSE_CC_CVSID \
"$Id$"
#include "lisousi.h"
/* create a 1/sqrt(t) filter response */
Tseries filterresponse(const int& size, const double& dt, const Options& opt)
{
if (opt.verbose)
{
cout << " construct 1/sqrt(t); ";
if (opt.nointeg)
{
cout << "tshift: " << opt.tshift << "; "
<< "tlim: " << opt.tlim << "; "
<< "tfac: " << opt.tfac << endl;
}
else
{
cout << "integshift: " << opt.integshift << " means "
<< opt.integshift*dt << "s time shift" << endl;
}
}
Tseries retval=Tseries(0,size-1);
aff::Iterator<Tseries> I(retval);
for (unsigned int i=0; i<retval.size(); ++i)
{
if (opt.nointeg)
{
double t=dt*(double(i)+opt.tshift);
if (t < (opt.tlim*dt)) { t=opt.tfac*dt; }
*I=sqrt(1./t);
}
else
{
if (i == 0)
{
double tref=1.-opt.integshift;
*I=2.*sqrt(tref/dt);
}
else
{
double tref=(double(i)-opt.integshift);
*I=2.*(sqrt(tref+1)-sqrt(tref))/sqrt(dt);
}
}
++I;
}
return(retval);
}
/* ----- END OF lisousi_filterresponse.cc ----- */
# this is <lisousi_help_text.txt>
# ============================================================================
# description of lisousi program options
# --------------------------------------
# $Id$
# ============================================================================
#
$Id$
outfile name of output file to contain results
infile input file
t:T select traces T, where T may be any range
specification like '3-4' or '5,6,7-12,20'
f:F specifies an input file format differing from
the format selected by "-type"
informational output:
-xhelp print detailed information regarding file formats
-v be verbose
-DEBUG produce debug output
data file input and output:
-o overwrite output
-type type choose input file format (default: sff)
-Type type choose output file format (default: sff)
filter alternatives:
default: Fourier domain application of analytic filter coefficents
-fdfilter apply numerically derived Fourier coefficients of 1/sqrt(t)
-tdfilter apply 1/sqrt(t) in the time domain by discrete convolution
-pad n in the case of Fourier domain processing, the time series
is padded with zeroes to factor n of the number of samples
of the input series prior to the Fourier transformation
taper alternatives:
-taperlast tapers are applied prior to filtering by default;this
option selects to taper the time series after filtering
If this option is selected, the order of the processing
steps (which are not commutative) is reversed.
-sqrttaper By default we use sqrt(2)*r/sqrt(t) as the taper and
scaling function. This is appropriate for direct
waves. In the case of reflected waves commonly the
propagation velocity is assumed to be known (see
option -velocity). Then the scaling function is
c*sqrt(2*t), where c is the assumed phase velocity.
The latter is selected by this option.
-tapdel t delay taper function by t seconds; this accounts for
the finite duration of a band limited source signal
-tapslo s delay taper by not more than s*offset, where s is a
slowness in s/m; serves as an upper limit for tapdel
at small offsets
1/sqrt(t) construction:
the default operation to handle the singularity at t=0 is:
choose samples of 1/sqrt(t) such that they produce the
same mean values as the continuous function when being
integrated over one sampling interval
-integshift f shift time axis by a fraction f of the sampling interval
in order to adjust centroid of integration interval
the following options are for experimental purposes only:
-nointeg do not derive samples by integration; sample 1/sqrt(t)
directly and optionally use the following parameters to
control the bevahiour close to t=0
deprecated method
if -nointeg is selected:
-tshift f to avoid the singularity of 1/sqrt(t) at t=0, sampling time
is shifted by a fraction f of the sampling interval
-tlim f to avoid the singularity of 1/sqrt(t) at t=0, the sample
time for samples at times smaller than f*dt, where dt is
the sampling interval, is set to a given fraction of the
sampling interval (see option -tfac)
-tfac f to avoid the singularity of 1/sqrt(t) at t=0, the sample
time at small times (see options -tlim) is set to a
fraction f of the sampling interval
Single velocity transformation:
-fredomain apply transformation in the frequency domain for a
single direct wave velocity.
-exact apply exact transformation for single velocity approach
(default is: far-field approximation)
-velocity v set assumed wave velocity in km/s
-transition r1,r2 mix both approaches; up to offset r1 the "single
velocity transformation" is used; for offsets larger than
r2 the "direct wave transformation" is used; in between
both are mixed with continously varying percentage
from one to the other
-limitlength some convolution schemes return a larger number of samples
than present in the input series (due to padding for
example); this option limits to numer of output samples
to the number of input samples
# ----- END OF lisousi_help_text.txt -----
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