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

stfinv debugging: added debugging output and array index shifts

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: 4005
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 43cd55e1
......@@ -41,6 +41,7 @@
#include <sstream>
#include <cmath>
#include <stfinv/stfinvbase.h>
#include <stfinv/debug.h>
namespace stfinv {
......@@ -49,7 +50,6 @@ namespace stfinv {
const std::string& parameters)
: Mtriples(triples), Mstf(stf), Mweights(0,triples.size()-1)
{
this->checkconsistency();
this->parseparameters(parameters);
// extract parameters
......@@ -57,6 +57,8 @@ namespace stfinv {
std::istringstream is(this->parameter("DEBUG","0"));
is >> Mdebug;
}
STFINV_debug(Mdebug>0, "STFBaseEngine::STFBaseEngine",
"initializing base class");
Mverbose=(this->parameter("verbose","false")=="true");
if (this->parameterisset("exp"))
{
......@@ -74,6 +76,8 @@ namespace stfinv {
Mweights=1.;
}
this->checkconsistency();
// align index ranges; first index will be zero
if (Mstf.series.first()!=0) { Mstf.series.shift(-Mstf.series.first()); }
stfinv::Tvectoroftriples::iterator I=Mtriples.begin();
......@@ -111,6 +115,7 @@ namespace stfinv {
"inconsistent number of samples");
STFINV_assert(std::abs(1.-header.sampling.dt/dt)<tolerance,
"inconsistent values of sampling interval");
++I;
} // while (I!=Mtriples.end())
}
......
......@@ -41,6 +41,7 @@
#include <iostream>
#include <aff/functions/sqrsum.h>
#include <stfinv/stfinvblinddeconv.h>
#include <stfinv/debug.h>
namespace stfinv {
......@@ -70,6 +71,8 @@ namespace stfinv {
// scale energy
std::istringstream is (this->parameter("waterlevel","1.e-3"));
is >> Mwaterlevel;
STFINV_debug(Mdebug&1, "STFEngineBlindDeconvolution::initialize()",
"Mwaterlevel=" << Mwaterlevel);
STFINV_assert(Mwaterlevel > 0,
"ERROR: parameter for option \"waterlevel\" not larger than 0");
} // void STFEngineBlindDeconvolution::initialize()
......@@ -95,8 +98,6 @@ namespace stfinv {
void STFEngineBlindDeconvolution::exec()
{
STFINV_abort("STFEngineBlindDeconvolution::run not yet implemented");
// read signals and calculate FFT
this->fftinput();
......@@ -109,20 +110,32 @@ namespace stfinv {
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
waterlevel *= 2.*Mwaterlevel;
STFINV_debug(Mdebug&2, "STFEngineBlindDeconvolution::exec()",
"waterlevel=" << waterlevel);
STFINV_debug(Mdebug&2, "STFEngineBlindDeconvolution::exec()",
"nfreq=" << this->nfreq());
// cycle through frequencies
for (unsigned int i=0; i<this->nfreq(); ++i)
{
TAspectrum synthetic=this->syntheticcoeff(i);
TAspectrum recording=this->recordingcoeff(i);
STFINV_debug(Mdebug&4, "STFEngineBlindDeconvolution::exec()",
"frequency=" << i);
TAspectrum syntheticref=this->syntheticcoeff(i);
TAspectrum recordingref=this->recordingcoeff(i);
STFINV_debug(Mdebug&4, "STFEngineBlindDeconvolution::exec()",
"syntheticref.first(0)=" << syntheticref.first(0) << " " <<
"syntheticref.last(0)=" << syntheticref.last(0) << " "
"syntheticref.first(1)=" << syntheticref.first(1) << " " <<
"syntheticref.last(1)=" << syntheticref.last(1) << " ");
TAspectrum::Tvalue numerator, denominator;
denominator=waterlevel;
numerator=0.;
for (unsigned int j=0; j<this->nreceivers(); ++j)
{
numerator += this->weight(j)*this->weight(j)
*conj(synthetic(j))*recording(j);
*conj(syntheticref(j))*recordingref(j);
denominator += this->weight(j)*this->weight(j)
*conj(synthetic(j))*synthetic(j);
*conj(syntheticref(j))*syntheticref(j);
}
this->stfcoeff(i)=numerator/denominator;
} // for (unsigned int i=0; i<this->nfreq(); ++i)
......
......@@ -40,6 +40,7 @@
#include <sstream>
#include <stfinv/stfinvfourier.h>
#include <stfinv/debug.h>
#include <aff/subarray.h>
#include <aff/slice.h>
......@@ -183,6 +184,12 @@ namespace stfinv {
nsamples);
Mfftengineoutput=fourier::fft::DRFFTWAFFArrayEngine(1+this->nreceivers(),
nsamples);
STFINV_debug(Mdebug&1, "STFEngineBlindDeconvolution::initialize()",
"this->nsamples()=" << this->nsamples() << " " <<
"padfactor=" << padfactor << " " <<
"divisor=" << divisor << " " <<
"nsamples=" << nsamples << " ");
} // void STFFourierDomainEngine::initialize()
/*----------------------------------------------------------------------*/
......@@ -207,7 +214,10 @@ namespace stfinv {
STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::recordingspec() const
{
TAspectrum inspecarray=Mfftengineinput.spectrum();
return(aff::subarray(inspecarray)()(0,this->nreceivers()-1));
TAspectrum subarray=aff::subarray(inspecarray)()(0,this->nreceivers()-1);
subarray.shape().setfirst(0,0);
subarray.shape().setfirst(1,0);
return(subarray);
} // STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::recordingspec() const
/*----------------------------------------------------------------------*/
......@@ -215,8 +225,11 @@ namespace stfinv {
STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::syntheticspec() const
{
TAspectrum inspecarray=Mfftengineinput.spectrum();
return(aff::subarray(inspecarray)()(this->nreceivers(),
(2*this->nreceivers())-1));
TAspectrum subarray=aff::subarray(inspecarray)()(this->nreceivers(),
(2*this->nreceivers())-1);
subarray.shape().setfirst(0,0);
subarray.shape().setfirst(1,0);
return(subarray);
} // STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::syntheticspec() const
/*----------------------------------------------------------------------*/
......@@ -232,7 +245,10 @@ namespace stfinv {
STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const
{
STFFourierDomainEngine::TAspectrum rec=this->recordingspec();
return(aff::slice(rec)(1,i));
STFFourierDomainEngine::TAspectrum slice=aff::slice(rec)(0,i);
slice.shape().setfirst(0,0);
slice.shape().setfirst(1,0);
return(slice);
} // STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const
/*----------------------------------------------------------------------*/
......@@ -240,8 +256,18 @@ namespace stfinv {
STFFourierDomainEngine::TAspectrum
STFFourierDomainEngine::syntheticcoeff(const unsigned int& i) const
{
STFINV_debug(Mdebug&8, "STFFourierDomainEngine::syntheticcoeff",
"i=" << i);
STFFourierDomainEngine::TAspectrum syn=this->syntheticspec();
return(aff::slice(syn)(1,i));
STFFourierDomainEngine::TAspectrum slice=aff::slice(syn)(0,i);
slice.shape().setfirst(0,0);
slice.shape().setfirst(1,0);
STFINV_debug(Mdebug&8, "STFFourierDomainEngine::syntheticcoeff",
"slice.first(0)=" << slice.first(0) << " " <<
"slice.last(0)=" << slice.last(0) << " "
"slice.first(1)=" << slice.first(1) << " " <<
"slice.last(1)=" << slice.last(1) << " ");
return(slice);
} // STFFourierDomainEngine::syntheticcoeff(const unsigned int& i) const
/*----------------------------------------------------------------------*/
......@@ -271,7 +297,7 @@ namespace stfinv {
void STFFourierDomainEngine::getinput()
{
// clear workspace
// clear workspace through reference
TAseries sarray=Mfftengineinput.series();
sarray=0.;
......@@ -279,11 +305,11 @@ namespace stfinv {
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
// get references to time series in workspace
TAseries recording=Mfftengineinput.series(i);
TAseries synthetic=Mfftengineinput.series(i+this->nreceivers());
TAseries recordingref=Mfftengineinput.series(i);
TAseries syntheticref=Mfftengineinput.series(i+this->nreceivers());
// copyin function copies as many elements as possible
recording.copyin(this->recording(i));
synthetic.copyin(this->synthetic(i));
recordingref.copyin(this->recording(i));
syntheticref.copyin(this->synthetic(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
} // void STFFourierDomainEngine::getinput()
......@@ -311,7 +337,9 @@ namespace stfinv {
{
TAspectrum synthetic=this->syntheticspec();
TAspectrum convolvedsynthetic
=aff::subarray(Mfftengineinput.spectrum())()(0,this->nreceivers()-1);
=aff::subarray(Mfftengineoutput.spectrum())()(0,this->nreceivers()-1);
convolvedsynthetic.shape().setfirst(0,0);
convolvedsynthetic.shape().setfirst(1,0);
TAspectrum stfcoeff=this->stfspec();
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
......
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