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

proceeding; next: read out shape explicitely

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/branches/libfourier
SVN Revision: 3937
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 1cac2e6f
......@@ -161,7 +161,7 @@ HEADERS=$(shell find . -name \*.h)
# in the binary version of the library
# (see below for the configuration of a preinstantiated version of template
# code)
SRC=fcommand.cc filters.cc polesnzeroes.cc fftwaff.cc error.cc
SRC=fcommand.cc filters.cc polesnzeroes.cc fftwaff.cc error.cc fftwaffar.cc
# test programs are placed in a subdirectory
TESTS=$(wildcard tests/*.cc)
# whereever we find a README, we will use it
......
......@@ -38,7 +38,8 @@
#define TF_FFTWAFFAR_CC_CVSID \
"$Id$"
#include <fourier::fftwaffar.h>
#include <fourier/fftwaffar.h>
#include <fourier/error.h>
namespace fourier {
......@@ -46,18 +47,18 @@ namespace fourier {
DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const int& nseis,
const int& nsamp)
: Mseries(TAseries(nseis,nsamp)),
Mspectrum(TAspectrum(nseis,DRFFTWAFFArrayEngine::ncoeff(nsamp))),
: Mseriesarray(TAseries(nseis,nsamp)),
Mspectrumarray(TAspectrum(nseis,DRFFTWAFFArrayEngine::ncoeff(nsamp))),
Mplanr2c(0), Mplanc2r(0)
{
this->checkconsistency;
this->checkconsistency();
}
/*----------------------------------------------------------------------*/
DRFFTWAFFArrayEngine::DRFFTWAFFArrayEngine(const TAseries& series,
const TAspectrum& spec)
: Mseriesarray(series), Mspectrumarrays(spec),
: Mseriesarray(series), Mspectrumarray(spec),
Mplanr2c(0), Mplanc2r(0)
{ this->checkconsistency(); }
......@@ -72,7 +73,7 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! delete plans.
void DRFFTWAFFArrayEngine::delete_plans() const
void DRFFTWAFFArrayEngine::delete_plans()
{
if (Mplanr2c != 0) { fftw_destroy_plan(Mplanr2c); }
if (Mplanc2r != 0) { fftw_destroy_plan(Mplanc2r); }
......@@ -94,69 +95,138 @@ namespace fourier {
/*----------------------------------------------------------------------*/
void DRFFTWAFFArrayEnginge::checkconsistency()
void DRFFTWAFFArrayEngine::checkconsistency()
{
FOURIER_assert(Mseries.size(0)==
DRFFTWAFFArrayEnginge::ncoeff(Mspectrum.size(0)),
FOURIER_assert(Mseriesarray.size(0)==
DRFFTWAFFArrayEngine::ncoeff(Mspectrumarray.size(0)),
"ERROR: sample size inconsistent");
FOURIER_assert(Mseries.size(1)==Mspectrum.size(1),
FOURIER_assert(Mseriesarray.size(1)==Mspectrumarray.size(1),
"ERROR: inconsistent number of signals");
FOURIER_assert(((Mseries.size(2)==1) && (Mspectrum.size(2)==1)),
FOURIER_assert(((Mseriesarray.size(2)==1) && (Mspectrumarray.size(2)==1)),
"ERROR: two-dimensional arrays only");
FOURIER_assert(((Mseries.size(3)==1) && (Mspectrum.size(3)==1)),
FOURIER_assert(((Mseriesarray.size(3)==1) && (Mspectrumarray.size(3)==1)),
"ERROR: two-dimensional arrays only");
} // void DRFFTWAFFArrayEnginge::checkconsistency()
} // void DRFFTWAFFArrayEngine::checkconsistency()
/*----------------------------------------------------------------------*/
void DRFFTWAFFArrayEnginge::createplanr2c()
void DRFFTWAFFArrayEngine::createplanr2c()
{
if (Mplanr2c==0)
{
// overall parameters
// ------------------
// tranformation only along one dimension
const int rank=1;
const int n[]={1};
// size of each series
int n=Mseriesarray.size(0);
// number of signals to be transformed
const int howmany=Mseries.size(1);
const int howmany=Mseriesarray.size(1);
// quick design
const unsigned flags=FFTW_ESTIMATE;
Mplanr2c=fftw_plan_many_dft_r2c(int rank,
const int *n,
int howmany,
double *in,
const int *inembed,
int istride,
int idist,
fftw_complex *out,
const int *onembed,
int ostride,
int odist,
unsigned flags);
// 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 DRFFTWAFFArrayEnginge::createplanr2c()
} // void DRFFTWAFFArrayEngine::createplanr2c()
/*----------------------------------------------------------------------*/
void DRFFTWAFFArrayEnginge::createplanc2r()
void DRFFTWAFFArrayEngine::createplanc2r()
{
if (Mplanc2r==0)
{
Mplanc2r=fftw_plan_many_dft_c2r(int rank,
const int *n,
int howmany,
fftw_complex *in,
const int *inembed,
int istride,
int idist,
double *out,
const int *onembed,
int ostride,
int odist,
unsigned flags);
// 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 DRFFTWAFFArrayEnginge::createplanc2r()
} // void DRFFTWAFFArrayEngine::createplanc2r()
} // namespace fft
......
......@@ -45,6 +45,7 @@
#include<complex>
#include<fftw3.h>
#include<aff/array.h>
#include<aff/series.h>
namespace fourier {
......@@ -58,6 +59,10 @@ namespace fourier {
* persistent workspace and makes full use of the reference semantics
* implemented in aff::Array.
*
* The class is designed for strided arrays. aff::Array has to be used
* together with aff::Strided. Currently (5/2011) there is no alternative
* in the aff::Array class template.
*
* When initializing DRFFTWAFFArrayEngine either the dimensions of the
* time series array must be passed or a readily prepared workspace is
* passed to the engine. In both cases the engine assumes that the user
......@@ -66,6 +71,9 @@ namespace fourier {
* actual transform. New values to be transformed are received through the
* common reference with the user program. A reference to the arrays is
* provided through member functions.
*
* The first index to the array is the sample index, the second index is
* the signal index.
*/
class DRFFTWAFFArrayEngine {
public:
......@@ -73,6 +81,8 @@ namespace fourier {
typedef std::complex<Tsample> Tcoeff;
typedef aff::Array<Tsample> TAseries;
typedef aff::Array<Tcoeff> TAspectrum;
typedef aff::Series<Tsample> Tseries;
typedef aff::Series<Tcoeff> Tspectrum;
DRFFTWAFFArrayEngine(const int& nseis,
const int& nsamp);
DRFFTWAFFArrayEngine(const TAseries& series,
......@@ -82,14 +92,22 @@ namespace fourier {
void c2r();
Tsample scale_series(const Tsample& dt) const;
Tsample scale_spectrum(const Tsample& dt) const;
//! \brief return a reference to the time series arrays
TAseries series() const { return Mseriesarray; }
//! \brief return a reference to the Fourier coefficient arrays
TAspectrum spectrum() const { return Mspectrumarray; }
//! \brief return the number of series in the arrays
int nseries() const;
//! \brief return a reference to the time series i
Tseries series(const int& i) const;
//! \brief return a reference to the Fourier coefficients of series i
Tspectrum specturm(const int& i) const;
private:
void checkconsistency();
void createplanr2c();
void createplanc2r();
void delete_plans() const;
static int ncoeff(const int& nsamples) const { return(nsamples/2+1); }
void delete_plans();
static int ncoeff(const int& nsamples) { return(nsamples/2+1); }
TAseries Mseriesarray;
TAspectrum Mspectrumarray;
fftw_plan Mplanr2c;
......
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