fftwaffar.cc

/*! \file fftwaffar.cc
 * \brief engine to transfrom several signals at once (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * $Id$
 * \author Thomas Forbriger
 * \date 13/05/2011
 * 
 * engine to transfrom several signals at once (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 
 *  - 13/05/2011   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define TF_FFTWAFFAR_CC_VERSION \
  "TF_FFTWAFFAR_CC   V1.0   "
#define TF_FFTWAFFAR_CC_CVSID \
  "$Id$"

#include <fourier/fftwaffar.h>
#include <fourier/error.h>

namespace fourier {

  namespace fft {

    DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const int& nseis,
                                               const int& nsamp)
      : Mseriesarray(TAseries(nseis,nsamp)),
      Mspectrumarray(TAspectrum(nseis,DRFFTWAFFArrayEngine::ncoeff(nsamp))),
      Mplanr2c(0), Mplanc2r(0) 
    {
      this->checkconsistency();
    }

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

    DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const TAseries& series,
                                               const TAspectrum& spec)
      : Mseriesarray(series), Mspectrumarray(spec),
      Mplanr2c(0), Mplanc2r(0) 
    { this->checkconsistency(); }

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

    //! delete plan.
    DRFFTWAFFArrayEngine::~DRFFTWAFFArrayEngine()
    {
      this->delete_plans();
    } // DRFFTWAFFArrayEngine::~DRFFTWAFFArrayEngine()

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

    //! delete plans.
    void DRFFTWAFFArrayEngine::delete_plans() 
    {
      if (Mplanr2c != 0) { fftw_destroy_plan(Mplanr2c); }
      if (Mplanc2r != 0) { fftw_destroy_plan(Mplanc2r); }
      Mplanr2c=0;
      Mplanc2r=0;
    } // DRFFTWAFFArrayEngine::delete_plans()

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

    //! execute r2c plan
    void DRFFTWAFFArrayEngine::r2c() 
    {
      if (Mplanr2c == 0)
      {

      }

    } // DRFFTWAFFArrayEngine::r2c()

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

    void DRFFTWAFFArrayEngine::checkconsistency()
    {
      FOURIER_assert(Mseriesarray.size(0)==
                     DRFFTWAFFArrayEngine::ncoeff(Mspectrumarray.size(0)),
                     "ERROR: sample size inconsistent");
      FOURIER_assert(Mseriesarray.size(1)==Mspectrumarray.size(1),
                     "ERROR: inconsistent number of signals");
      FOURIER_assert(((Mseriesarray.size(2)==1) && (Mspectrumarray.size(2)==1)),
                     "ERROR: two-dimensional arrays only");
      FOURIER_assert(((Mseriesarray.size(3)==1) && (Mspectrumarray.size(3)==1)),
                     "ERROR: two-dimensional arrays only");
    } // void DRFFTWAFFArrayEngine::checkconsistency()

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

    void DRFFTWAFFArrayEngine::createplanr2c()
    {
      if (Mplanr2c==0)
      {
        // overall parameters
        // ------------------
        // tranformation only along one dimension
        const int rank=1;
        // size of each series
        int n=Mseriesarray.size(0);
        // number of signals to be transformed
        const int howmany=Mseriesarray.size(1);
        // quick design
        const unsigned flags=FFTW_ESTIMATE;

        // input array
        // -----------
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
        const int istride=Mseriesarray.stride(0);
        // distance to next signal
        const int idist=Mseriesarray.stride(1);
        // reference to memory
        TAseries::Trepresentation irep=Mseriesarray.representation();
        // offset in array representation
        aff::Tsubscript ioffset=Mseriesarray.first_offset();
        // pointer to first array element in memory
        TAseries::Trepresentation::Tpointer ipointer=&irep[ioffset];
        // casted pointer
        double* in=aff::SizeCheckedCast<TAseries::Tvalue,double>::cast(ipointer);

        // output array
        // ------------
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
        const int ostride=Mspectrumarray.stride(0);
        // distance to next signal
        const int odist=Mspectrumarray.stride(1);
        // reference to memory
        TAspectrum::Trepresentation orep=Mspectrumarray.representation();
        // offset in array representation
        aff::Tsubscript ooffset=Mspectrumarray.first_offset();
        // pointer to first array element in memory
        TAspectrum::Trepresentation::Tpointer opointer=&orep[ooffset];
        // casted pointer
        fftw_complex* out=aff::SizeCheckedCast<TAspectrum::Tvalue,fftw_complex>::cast(opointer);

        // create plan
        Mplanr2c=fftw_plan_many_dft_r2c(rank, &n, howmany,
                                        in, inembed, istride, idist,
                                        out, onembed, ostride, odist,
                                        flags);
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating r2c plan");
      }
    } // void DRFFTWAFFArrayEngine::createplanr2c()

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

    void DRFFTWAFFArrayEngine::createplanc2r()
    {
      if (Mplanc2r==0)
      {
        // overall parameters
        // ------------------
        // tranformation only along one dimension
        const int rank=1;
        // size of each series
        int n=Mseriesarray.size(0);
        // number of signals to be transformed
        const int howmany=Mseriesarray.size(1);
        // quick design
        const unsigned flags=FFTW_ESTIMATE;

        // input array
        // -----------
        // one-dimensional transform: use default
        int* inembed=0;
        // distance to next sample
        const int istride=Mspectrumarray.stride(0);
        // distance to next signal
        const int idist=Mspectrumarray.stride(1);
        // reference to memory
        TAspectrum::Trepresentation irep=Mspectrumarray.representation();
        // offset in array representation
        aff::Tsubscript ioffset=Mspectrumarray.first_offset();
        // pointer to first array element in memory
        TAspectrum::Trepresentation::Tpointer ipointer=&irep[ioffset];
        // casted pointer
        fftw_complex* in=aff::SizeCheckedCast<TAspectrum::Tvalue,fftw_complex>::cast(ipointer);

        // output array
        // ------------
        // one-dimensional transform: use default
        int* onembed=0;
        // distance to next sample
        const int ostride=Mseries.stride(0);
        // distance to next signal
        const int odist=Mseries.stride(1);
        // reference to memory
        TAseries::Trepresentation orep=Mseries.representation();
        // offset in array representation
        aff::Tsubscript ooffset=Mseries.first_offset();
        // pointer to first array element in memory
        TAseries::Trepresentation::Tpointer opointer=&orep[ooffset];
        // casted pointer
        double* out=aff::SizeCheckedCast<TAseries::Tvalue,double>::cast(opointer);

        // create plan
        Mplanr2c=fftw_plan_many_dft_c2r(rank, &n, howmany,
                                        in, inembed, istride, idist,
                                        out, onembed, ostride, odist,
                                        flags);
        FOURIER_assert(Mplanr2c!=0, "ERROR: creating c2r plan");
      }
    } // void DRFFTWAFFArrayEngine::createplanc2r()

  } // namespace fft

} // namespace fourier

/* ----- END OF fftwaffar.cc ----- */