/*! \file stfinvfourier.cc
 * \brief a base class for all engines which operate in the Fourier domain (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id$
 * \author Thomas Forbriger
 * \date 08/05/2011
 * 
 * a base class for all engines which operate in the Fourier domain (implementation)
 * 
 * Copyright (c) 2011 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 
 *  - 08/05/2011   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define STFINV_STFINVFOURIER_CC_VERSION \
  "STFINV_STFINVFOURIER_CC   V1.0   "
#define STFINV_STFINVFOURIER_CC_CVSID \
  "$Id$"

#include <sstream>
#include <stfinv/stfinvfourier.h>
#include <aff/subarray.h>
#include <aff/slice.h>

namespace stfinv {

  /*! 
   * - Refernces to signal storage are available
   * - Parameters are parsed and available through the base class
   * What has to be done here:
   * -# The FFT engine must be initialized appropriately
   *    -# Time domain padding has to be read from the parameters
   *    -# Number of sample modulo has to be read from parameters
   *    -# number of samples must be calculated from these values
   *    -# based on these values workspace arrays must be constructed
   *    -# The FFT engine must be constructed from the workspace arrays.
   *       This requires a copy constructor for the FFT engine, which
   *       currently is not present.
   *       Alternatively, we could only use a pointer to the engine.
   */
  STFFourierDomainEngine::STFFourierDomainEngine(const stfinv::Tvectoroftriples& triples,
                                                 const stfinv::Waveform& stf,
                                                 const std::string& parameters)
    :Tbase(triples, stf, parameters)
  { 
    this->initialize();
  } // STFFourierDomainEngine::STFFourierDomainEngine

  /*----------------------------------------------------------------------*/

  void STFFourierDomainEngine::help(std::ostream& os) const
  {
    STFFourierDomainEngine::classhelp(os);
  } // void STFFourierDomainEngine::help(std::ostream& os) const


  /*----------------------------------------------------------------------*/

  const char* STFFourierDomainEngine::name() const
  {
    return("STFFourierDomainEngine");
  } //  const char const* STFFourierDomainEngine::name() const
  /*----------------------------------------------------------------------*/

  /*! \brief online help text giving information on options
   *
   * This must be kept synchronous with the options used by
   * STFFourierDomainEngine::initialize()
   */
  void STFFourierDomainEngine::classhelp(std::ostream& os)
  {
    os << "Options and parameters in common for Fourier engines:\n"
      << "fpad=f       padding factor (default: 1.5)\n"
      << "fpow2        use power of two\n"
      << "fdiv=d       use integer multiple of d\n"
      << "These options define the number of samples N used for the FFT.\n"
      << "This should be larger than the number of samples M in the\n"
      << "original series to avoid wrap-around. N=M*f at least. If fpow2\n"
      << "is set, N will be the next power of 2 larger than M*f. Else if\n"
      << "fdiv is set, N will be the next integer multiple of d larger\n"
      << "than M*f." << std::endl;
    STFBaseEngine::classhelp(os);
  } // void STFFourierDomainEngine::classhelp(std::ostream& os)

  /*----------------------------------------------------------------------*/

  /*! \brief Create FFT engines
   *
   * \par Number of samples
   * For the FFT we usually use a number of samples larger than that of the
   * underlying time series. 
   * This way we at least partially avoid the undesired effect of the cyclic
   * discrete Fourier transform, which can lead to wrap-around.
   * Currently there are three ways to modify the number of samples actually
   * used:
   * -# Option \c fpad specifies an padding factor. 
   *    The number of samples simply is larger than the number of samples of
   *    the underlying time series by this factor.
   * -# Option \c fpow2:
   *    If set, the next power of two larger than the fpad*nsamples is used.
   * -# Option \c fdiv specifies a divisor.
   *    If set and if \c fpow2 is not set the number of samples will be the
   *    next integer multiple of the divisor larger than fpad*nsamples.
   *
   * \par Workspace
   * Two FFT engines will be created:
   * -# One engine (\c STFFourierDomainEngine::Mfftengineinput) being shared
   *  by recorded data and synthetic data, because both have to be transformed
   *  to Fourier domain at once.
   * -# One engine (\c STFFourierDomainEngine::Mfftengineoutput) being shared
   *  by the stf and the convolved synthetics, because both have to be
   *  transformed to the time domain at once.
   */
  void STFFourierDomainEngine::initialize()
  {
    // extract parameter values
    // ------------------------
    // padding factor
    double padfactor;
    {
      std::istringstream is(this->parameter("fpad","1.5"));
      is >> padfactor;
    }
    STFINV_assert(padfactor >= 1.,
                  "ERROR: parameter for option \"fpad\" not larger or equal 1");
    // flag: use power of two
    bool poweroftwo=(this->parameter("fpow2","false")=="true");
    // flag: use integer multiples of divisor
    bool divisorset=this->parameterisset("fdiv");
    // number of samples shall be integer power of divisor
    unsigned int divisor=1;
    if (divisorset)
    {
      std::istringstream is (this->parameter("fdiv","100")); 
      is >> divisor;
      STFINV_assert(divisor > 0,
                  "ERROR: parameter for option \"fdiv\" not larger than 0");
    }

    // define number of samples to be used by Fourier engine
    // -----------------------------------------------------
    // use at least the number of samples sepcified by the padfactor
    unsigned int nsamples=padfactor*this->nsamples();
    if (nsamples<this->nsamples()) { nsamples=this->nsamples(); }
    // power of two has precendence
    if (poweroftwo)
    {
      unsigned int n=2;
      while (n<nsamples) { n *= 2; }
      nsamples=n;
    }
    else if (divisorset)
    {
      unsigned int rest=nsamples % divisor;
      nsamples -= rest;
      nsamples += divisor;
    } // if (poweroftwo) or (divisorset)

    // allocate workspace by creating engines
    // --------------------------------------
    Mfftengineinput=fourier::fft::DRFFTWAFFArrayEngine(2*this->nreceivers(),
                                                       nsamples);
    Mfftengineoutput=fourier::fft::DRFFTWAFFArrayEngine(1+this->nreceivers(),
                                                        nsamples);
  } // void STFFourierDomainEngine::initialize() 

  /*----------------------------------------------------------------------*/

  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)
  } // void STFFourierDomainEngine::fftinput() 

  /*----------------------------------------------------------------------*/

  void STFFourierDomainEngine::fftoutput()
  {
    // 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::fftoutput() 

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::recordingspec() const
  {
    TAspectrum inspecarray=Mfftengineinput.spectrum();
    return(aff::subarray(inspecarray)()(0,this->nreceivers()-1));
  } // STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::recordingspec() const

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::syntheticspec() const
  {
    TAspectrum inspecarray=Mfftengineinput.spectrum();
    return(aff::subarray(inspecarray)()(this->nreceivers(),
                                       (2*this->nreceivers())-1));
  } // STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::syntheticspec() const

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::stfspec() const
  {
    return(Mfftengineoutput.spectrum(this->nreceivers()));
  } // STFFourierDomainEngine::Tspectrum STFFourierDomainEngine::stfspec() const

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum 
    STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const
  {
    STFFourierDomainEngine::TAspectrum rec=this->recordingspec();
    return(aff::slice(rec)(1,i));
  } // STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum 
    STFFourierDomainEngine::syntheticcoeff(const unsigned int& i) const
  {
    STFFourierDomainEngine::TAspectrum syn=this->syntheticspec();
    return(aff::slice(syn)(1,i));
  } // STFFourierDomainEngine::syntheticcoeff(const unsigned int& i) const

  /*----------------------------------------------------------------------*/

  STFFourierDomainEngine::TAspectrum::Tvalue& 
    STFFourierDomainEngine::stfcoeff(const unsigned int& i) const
  {
    STFFourierDomainEngine::TAspectrum stffft=this->stfspec();
    return(stffft(i));
  } // STFFourierDomainEngine::stfcoeff(const unsigned int& i) const

  /*----------------------------------------------------------------------*/

  double STFFourierDomainEngine::frequency(const unsigned int& i) const
  {
    return(static_cast<double>(i)/(this->dt()*Mfftengineinput.nsamples()));
  } // double STFFourierDomainEngine::frequency(const unsigned int& i) const

  /*----------------------------------------------------------------------*/

  unsigned int STFFourierDomainEngine::nfreq() const
  {
    return(Mfftengineinput.nfrequencies());
  } // unsigned int STFFourierDomainEngine::nfreq() const

} // namespace stfinv

/* ----- END OF stfinvfourier.cc ----- */