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

blind deconvolution is implemented and compiles

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: 3989
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 6fda7fb5
......@@ -39,6 +39,7 @@
"$Id$"
#include <iostream>
#include <aff/functions/sqrsum.h>
#include <stfinv/stfinvblinddeconv.h>
namespace stfinv {
......@@ -95,6 +96,39 @@ namespace stfinv {
void STFEngineBlindDeconvolution::exec()
{
STFINV_abort("STFEngineBlindDeconvolution::run not yet implemented");
// read signals and calculate FFT
this->fftinput();
// calculate waterlevel
double waterlevel=0.;
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
waterlevel += this->weight(i)*this->weight(i)
*aff::func::sqrsum(this->synthetic(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
waterlevel *= 2.*Mwaterlevel;
// cycle through frequencies
for (unsigned int i=0; i<this->nfreq(); ++i)
{
TAspectrum synthetic=this->syntheticcoeff(i);
TAspectrum recording=this->recordingcoeff(i);
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);
denominator += this->weight(j)*this->weight(j)
*conj(synthetic(j))*synthetic(j);
}
this->stfcoeff(i)=numerator/denominator;
} // for (unsigned int i=0; i<this->nfreq(); ++i)
// provide results to user
this->fftoutput();
} // void STFEngineBlindDeconvolution::exec()
} // namespace stfinv
......
......@@ -111,13 +111,13 @@ namespace stfinv {
* \f$d_{lk}-s_{lk}q_l\f$ equals the scaled energy of the synthetics
* \f$s_{lk}\f$ and
* \f[
* \eta^2=\frac{1}{\sum\limits_k f_k\sum\limits_l \left|s_{lk}\right|^2}
* \eta^2=\frac{1}{\sum\limits_k f_k^2\sum\limits_l \left|s_{lk}\right|^2}
* \f]
* is the reciprocal of the scaled energy of the synthetics.
* If we then choose
* \f[
* \frac{\lambda^2}{\eta^2}=\frac{\epsilon^2}{N\eta^2}=
* \frac{\epsilon^2}{N}\sum\limits_k f_k\sum\limits_{l=0}^{N-1}
* \frac{\epsilon^2}{N}\sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}
* \left|s_{lk}\right|^2
* \f]
* where \f$N\f$ is the number of frequencies, then \f$\epsilon^2\f$
......@@ -142,10 +142,10 @@ namespace stfinv {
* In the above calculation the energy sum only uses the positive
* frequencies
* \f[
* \sum\limits_k f_k\sum\limits_{l=0}^{N-1}\left|\tilde{s}_{lk}\right|^2
* \sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}\left|\tilde{s}_{lk}\right|^2
* =
* 2N\,(\Delta t)^2\,
* \sum\limits_k f_k
* \sum\limits_k f_k^2
* \sum\limits_{k=0}^{2N-1}\left|S_{lk}\right|^2.
* \f]
* Fourier coefficients \f$s_{lk}\f$ calculated by the
......@@ -156,10 +156,10 @@ namespace stfinv {
* \f]
* Consequently
* \f[
* \sum\limits_k f_k\sum\limits_{l=0}^{N-1}\left|s_{lk}\right|^2
* \sum\limits_k f_k^2\sum\limits_{l=0}^{N-1}\left|s_{lk}\right|^2
* =
* 2N\,
* \sum\limits_k f_k
* \sum\limits_k f_k^2
* \sum\limits_{k=0}^{2N-1}\left|S_{lk}\right|^2.
* \f]
*
......@@ -169,7 +169,7 @@ namespace stfinv {
* q_l=\frac{
* \sum\limits_{k}f_k^2\,s_{kl}^\ast\,d_{kl}
* }{
* 2\epsilon^2\,\sum\limits_k f_k
* 2\epsilon^2\,\sum\limits_k f_k^2
* \sum\limits_{k=0}^{2N-1}\left|S_{lk}\right|^2
* +\sum\limits_{k}f_k^2\,s_{kl}^\ast\,s_{kl}
* }
......
......@@ -189,39 +189,17 @@ namespace stfinv {
void STFFourierDomainEngine::fftinput()
{
// clear workspace
TAseries sarray=Mfftengineinput.series();
sarray=0.;
// cycle through receivers
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());
// copyin function copies as many elements as possible
recording.copyin(this->recording(i));
synthetic.copyin(this->synthetic(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
this->getinput();
Mfftengineinput.r2c();
} // void STFFourierDomainEngine::fftinput()
/*----------------------------------------------------------------------*/
void STFFourierDomainEngine::fftoutput()
{
STFINV_abort("convolution is not yet implemented");
// cycle through receivers
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
// get references to time series
stfinv::Tseries convolvedsynthetics=this->convolvedsynthetic(i);
// copyin function copies as many elements as possible
convolvedsynthetics.copyin(Mfftengineoutput.series(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
// copy stf too
stfinv::Tseries stf=this->stf();
stf.copyin(Mfftengineoutput.series(this->nreceivers()));
this->convolve();
Mfftengineoutput.c2r();
this->putoutput();
} // void STFFourierDomainEngine::fftoutput()
/*----------------------------------------------------------------------*/
......@@ -289,6 +267,61 @@ namespace stfinv {
return(Mfftengineinput.nfrequencies());
} // unsigned int STFFourierDomainEngine::nfreq() const
/*----------------------------------------------------------------------*/
void STFFourierDomainEngine::getinput()
{
// clear workspace
TAseries sarray=Mfftengineinput.series();
sarray=0.;
// cycle through receivers
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());
// copyin function copies as many elements as possible
recording.copyin(this->recording(i));
synthetic.copyin(this->synthetic(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
} // void STFFourierDomainEngine::getinput()
/*----------------------------------------------------------------------*/
void STFFourierDomainEngine::putoutput()
{
// cycle through receivers
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
// get references to time series
stfinv::Tseries convolvedsynthetics=this->convolvedsynthetic(i);
// copyin function copies as many elements as possible
convolvedsynthetics.copyin(Mfftengineoutput.series(i));
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
// copy stf too
stfinv::Tseries stf=this->stf();
stf.copyin(Mfftengineoutput.series(this->nreceivers()));
} // void STFFourierDomainEngine::putoutput()
/*----------------------------------------------------------------------*/
void STFFourierDomainEngine::convolve()
{
TAspectrum synthetic=this->syntheticspec();
TAspectrum convolvedsynthetic
=aff::subarray(Mfftengineinput.spectrum())()(0,this->nreceivers()-1);
TAspectrum stfcoeff=this->stfspec();
for (unsigned int i=0; i<this->nreceivers(); ++i)
{
for (unsigned int j=0; j<this->nfreq(); ++j)
{
convolvedsynthetic(j,i)=synthetic(j,i)*stfcoeff(j);
} // for (unsigned int j=0; j<this->nfreq(); ++j)
} // for (unsigned int i=0; i<this->nreceivers(); ++i)
} // void STFFourierDomainEngine::convolve()
} // namespace stfinv
/* ----- END OF stfinvfourier.cc ----- */
......@@ -136,6 +136,9 @@ namespace stfinv {
* to user memory
*/
void putoutput();
/*! \brief convolve synthetics with stf
*/
void convolve();
// member data
// -----------
......
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