GPLIB++
MultiRecCalc.cpp
Go to the documentation of this file.
1 #include "MultiRecCalc.h"
2 #include "TimeFrequency.h"
3 #include "FatalException.h"
4 #include "CRFSpecElement.h"
5 #include "TsSpectrum.h"
6 #include "WFunc.h"
7 #include "VecMat.h"
8 #include "Adaptors.h"
9 
10 namespace ublas = boost::numeric::ublas;
11 using namespace std;
12 
13 namespace gplib
14  {
15  MultiRecCalc::MultiRecCalc(const int myshift, const double mysigma,
16  const double myc) :
17  c(myc), sigma(mysigma), shift(myshift)
18  {
19  }
20 
22  {
23  }
24 
28  SeismicDataComp &Receiver)
29  {
30  if (RadComp.size() != VerComp.size())
31  throw FatalException(
32  "Not the same number of stations for vertical and radial components !");
33  const unsigned int compsize = RadComp.front()->GetData().size();
34  const unsigned int ncomp = RadComp.size();
35  const double omegastep = 2. * PI / (RadComp.front()->GetDt()
36  * RadComp.front()->GetData().size());
37  gplib::rvec ConcatRad(compsize * ncomp), ConcatVer(compsize * ncomp);
38  for (unsigned int i = 0; i < ncomp; ++i)
39  {
40  if (RadComp.at(i)->GetData().size() != compsize
41  || VerComp.at(i)->GetData().size() != compsize)
42  throw FatalException(
43  "Incompatible sizes for Vertical or Radial component.");
44  copy(RadComp.at(i)->GetData().begin(),
45  RadComp.at(i)->GetData().end(), ConcatRad.begin() + i
46  * compsize);
47  copy(VerComp.at(i)->GetData().begin(),
48  VerComp.at(i)->GetData().end(), ConcatVer.begin() + i
49  * compsize);
50  }
51  gplib::cmat radtfspectrum = TimeFrequency(ConcatRad.begin(),
52  ConcatRad.end(), compsize, Steep());
53  gplib::cmat vertfspectrum = TimeFrequency(ConcatVer.begin(),
54  ConcatVer.end(), compsize, Steep());
55  gplib::cvec RecSpec(compsize / 2 + 1);
56  RecSpec *= 0.0;
57  for (unsigned int i = 0; i < radtfspectrum.size1(); ++i)
58  {
59  complex<double> delta = *max_element(
60  column(vertfspectrum, i).begin(),
61  column(vertfspectrum, i).end(), gplib::absLess<tcomp, tcomp>())
62  * c;
63  delta *= conj(delta);
64  //delta /= radtfspectrum.size1();
65  for (unsigned int j = 0; j < radtfspectrum.size2(); ++j)
66  {
67  RecSpec(j) += CalcSpectralElement(delta)(radtfspectrum(i, j),
68  vertfspectrum(i, j));
69  }
70  }
71  const double denom = 1. / (4.0 * sigma * sigma);
72  for (unsigned int i = 0; i < RecSpec.size(); ++i)
73  {
74  double freq = omegastep * i;
75  RecSpec(i) *= exp(-(freq * freq) * denom);
76  }
77  Receiver.SetB(-shift);
78  Receiver.SetDt(RadComp.front()->GetDt());
79  Receiver.GetData().assign(compsize, 0.0);
80  TsSpectrum Spectrum(false);
81  Spectrum.CalcTimeSeries(RecSpec.begin(), RecSpec.end(),
82  Receiver.GetData().begin(), Receiver.GetData().end());
83  const int shiftpoints = boost::numeric_cast<size_t>(shift
84  / RadComp.front()->GetDt());
85  rotate(Receiver.GetData().begin(), Receiver.GetData().end()
86  - shiftpoints, Receiver.GetData().end());
87 /* const double maxamp = 1. / *max_element(Receiver.GetData().begin(),
88  Receiver.GetData().end());
89  transform(Receiver.GetData().begin(), Receiver.GetData().end(),
90  Receiver.GetData().begin(), boost::bind(multiplies<double> (), _1,
91  maxamp));*/
92  }
93  }
std::vector< double > & GetData()
Access for data vector, for ease of use and efficiency we return a reference.
std::vector< boost::shared_ptr< SeismicDataComp > > tseiscompvector
This functor rises steeply at the edges and then has a wide range where it is unity.
Definition: WFunc.h:51
gplib::cmat TimeFrequency(InputIterator tsbegin, InputIterator tsend, const size_t seglength, WindowFunctype WFunc)
Calculate a sliding windowed fourier transform for a time series and store the results for each segme...
Definition: TimeFrequency.h:23
void SetB(const double theb)
The class CTsSpectrum is used to calculate spectra from (real) time series data.
Definition: TsSpectrum.h:21
void SetDt(const double dt)
Set delta t in s.
void CalcTimeSeries(_InputIterator freqbegin, _InputIterator freqend, _OutputIterator tsbegin, _OutputIterator tsend)
The member function to calculate a time series from (complex) spectra.
Definition: TsSpectrum.h:106
This class calculates one spectral element of the receiver function from the two input spectral eleme...
void CalcRecData(const SeismicStationList::tseiscompvector &RadComp, const SeismicStationList::tseiscompvector &VerComp, SeismicDataComp &Receiver)
Calculate the receiver function by spectral division.
The basic exception class for all errors that arise in gplib.