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

resolved problem with non-unique number of samples

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.
fixed a bug un fftwaffar.h
added test to cxxfftwartest.cc


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 3960
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent d37a9625
......@@ -88,7 +88,7 @@ using std::endl;
// a struct to store values of command line arguments
struct Options {
// be verbose
bool verbose, verbsize;
bool verbose, verbsize, exchangeengines;
// file format to be used
std::string fileformat;
// overwrite output files in case they already exist
......@@ -130,7 +130,7 @@ int main(int iargc, char* argv[])
char usage_text[]=
{
CXXFFTWARTEST_VERSION "\n"
"usage: cxxfftwartest input output [-s] [-v] [-o] [-type f]" "\n"
"usage: cxxfftwartest input output [-s] [-v] [-o] [-x] [-type f]" "\n"
" or: cxxfftwartest --help|-h" "\n"
" or: cxxfftwartest --xhelp" "\n"
};
......@@ -146,7 +146,15 @@ int main(int iargc, char* argv[])
"-v verbose mode\n"
"-s report container sizes\n"
"-o overwrite output file\n"
"-x reverse order of DRFFTWAFF and DRFFTWAFFArrayEngine\n"
" (see below)\n"
"-type f file format\n"
"\n"
"Fourier coefficients for input time series are obtained from\n"
"DRFFTWAFFArrayEngine. Then time series are constructed from these\n"
"coefficients by application of DRFFTWAFF. If the -x option is\n"
"selected, both transformation tools are exchanged in the\n"
"application order.\n"
};
// define commandline options
......@@ -165,6 +173,8 @@ int main(int iargc, char* argv[])
{"type",arg_no,"sff"},
// 5: report container sizes
{"s",arg_no,"-"},
// 6: exchange engines
{"x",arg_no,"-"},
{NULL}
};
......@@ -202,6 +212,7 @@ int main(int iargc, char* argv[])
opt.overwrite=cmdline.optset(3);
opt.fileformat=cmdline.string_arg(4);
opt.verbsize=cmdline.optset(5);
opt.exchangeengines=cmdline.optset(6);
// report program version if in verbose mode
if (opt.verbose)
......@@ -291,50 +302,90 @@ int main(int iargc, char* argv[])
/*----------------------------------------------------------------------*/
// prepare array engine
if (opt.verbose) { cout << "create FFT array engine" << endl; }
fourier::fft::DRFFTWAFFArrayEngine engine(ntraces, nsamples);
if (opt.verbose) { cout << "fill FFT array engine" << endl; }
engine.series()=0;
// fill sample array
for (int i=0; i<engine.nseries(); ++i)
if (opt.exchangeengines)
{
engine.series(i).copyin(vecoftraces[i].series);
if (opt.verbsize)
// prepare array engine
if (opt.verbose) { cout << "create FFT array engine" << endl; }
fourier::fft::DRFFTWAFFArrayEngine engine(ntraces, nsamples);
if (opt.verbose)
{
cout << "engine series #" << i << " has " << engine.series(i).size()
<< " samples" << endl;
cout << "fill FFT array engine with Fourier coefficients" << endl;
}
}
engine.spectrum()=0;
if (opt.verbose) { cout << "execute FFT array engine" << endl; }
engine.r2c();
fourier::fft::DRFFTWAFF fft;
fft.size(nsamples);
if (opt.verbose)
{
cout << "read coefficients and transform to time domain" << endl;
for (int i=0; i<engine.nseries(); ++i)
{
fourier::fft::DRFFTWAFF::Tspectrum
spectrum=fft(vecoftraces[i].series,vecoftraces[i].wid2.dt);
engine.spectrum(i).copyin(spectrum);
}
if (opt.verbose) { cout << "execute FFT array engine" << endl; }
engine.c2r();
if (opt.verbose)
{
cout << "read time series samples from FFT array engine" << endl;
}
for (int i=0; i<engine.nseries(); ++i)
{
Trace& trace=vecoftraces[i];
Tseries& series=trace.series;
series.copyin(engine.series(i)*engine.scale_series(trace.wid2.dt));
}
}
fourier::fft::DRFFTWAFF fft;
for (int i=0; i<engine.nseries(); ++i)
else // if (opt.exchangeengines)
{
Trace& trace=vecoftraces[i];
Tseries& series=trace.series;
fourier::fft::DRFFTWAFFArrayEngine::TAspectrum coeff
=engine.spectrum(i)*engine.scale_spectrum(trace.wid2.dt);
fourier::fft::DRFFTWAFF::Tspectrum spectrum(coeff.size());
spectrum.copyin(coeff);
series=fft(spectrum,trace.wid2.dt);
if (opt.verbsize)
// prepare array engine
if (opt.verbose) { cout << "create FFT array engine" << endl; }
fourier::fft::DRFFTWAFFArrayEngine engine(ntraces, nsamples);
if (opt.verbose) { cout << "fill FFT array engine" << endl; }
engine.series()=0;
// fill sample array
for (int i=0; i<engine.nseries(); ++i)
{
cout << "after return from engine for trace #" << i << endl;
cout << " size of engine spectrum: " << engine.spectrum(i).size()
<< endl;
cout << " size of copy of spectrum: " << spectrum.size()
<< endl;
cout << " returned series size: " << series.size() << endl;
engine.series(i).copyin(vecoftraces[i].series);
if (opt.verbsize)
{
cout << "engine series #" << i << " has " << engine.series(i).size()
<< " samples" << endl;
}
}
}
if (opt.verbose) { cout << "execute FFT array engine" << endl; }
engine.r2c();
if (opt.verbose)
{
cout << "read coefficients and transform to time domain" << endl;
}
fourier::fft::DRFFTWAFF fft;
fft.size(nsamples);
for (int i=0; i<engine.nseries(); ++i)
{
Trace& trace=vecoftraces[i];
Tseries& series=trace.series;
fourier::fft::DRFFTWAFFArrayEngine::TAspectrum coeff
=engine.spectrum(i)*engine.scale_spectrum(trace.wid2.dt);
fourier::fft::DRFFTWAFF::Tspectrum spectrum(coeff.size());
spectrum.copyin(coeff);
series=fft(spectrum,trace.wid2.dt);
if (opt.verbsize)
{
cout << "after return from engine for trace #" << i << endl;
cout << " size of engine spectrum: " << engine.spectrum(i).size()
<< endl;
cout << " size of copy of spectrum: " << spectrum.size()
<< endl;
cout << " returned series size: " << series.size() << endl;
}
}
} // if (opt.exchangeengines)
/*----------------------------------------------------------------------*/
......
......@@ -26,6 +26,7 @@
#define TF_FFTWAFF_CC_CVSID \
"$Id$"
#include <iostream>
#include <fourier/fftwaff.h>
#include <fourier/error.h>
#include <fourier/error.h>
......@@ -95,7 +96,7 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! prepare FFT settings for size n.
void DRFFTWAFF::set_size(const int& n) const
void DRFFTWAFF::set_size(const unsigned int& n) const
{
if (n != Msize)
{
......@@ -244,12 +245,18 @@ namespace fourier {
DRFFTWAFF::Tseries DRFFTWAFF::operator()(const Tspectrum::Tcoc& s,
const bool& debug) const
{
// adjust FFT size
int seriessize=s.size()*2-1;
// is Nyquits coefficients real?
if (std::abs(s(s.size()).imag()) < 1.e-8*std::abs(s(s.size()).real()))
{ --seriessize; }
this->set_size(seriessize);
// check number of expected Fourier coefficients
if (this->ssize() != s.size())
{
// adjust FFT size
int seriessize=s.size()*2-1;
// is Nyquits coefficients real?
// this measure can only be effective if the Nyquist coefficient is
// finite, which will not be the case for most signals
if (std::abs(s(s.size()).imag()) < 1.e-8*std::abs(s(s.size()).real()))
{ --seriessize; }
this->set_size(seriessize);
}
#ifdef FFTWFALLBACK
Tseries retval(Msize);
......
......@@ -89,6 +89,25 @@ namespace fourier {
* transform (usual convention with \f$dt\f$ and \f$df\f$
* integrals - not \f$d \omega\f$).
*
* \note
* The appropriate number of samples for the time series obtained from a
* given set of Fourier coefficients is non-unique, if only the Fourier
* coefficients for positive frequencies are given, as is the case here.
* A time series with \f$2n\f$ samples and a time series with \f$2n+1\f$
* samples both will results in a set of \f$n+1\f$ Fourier coefficients.
* The method to determine the required number of time samples from the
* imaginary part of the Nyquist Fourier coefficients only works if the
* Nyquist coefficient (at least real part) is finite - which is not the
* case for most of our signals. This way we will obtain \f$n+1\f$ Fourier
* coefficients from \f$2n\f$ time series samples. When reconstructing the
* time series, we will obtain \f$2n+1\f$ samples with a sampling interval
* now being \f$2n/(2n+1)\f$ of the original interval. For this reason it
* is appropriate to set the size of the expected time series explicitely
* either by using the constructor
* DRFFTWAFF::DRFFTWAFF(const unsigned int& n) or by
* calling DRFFTWAFF::size(const unsigned int& s) prior to
* DRFFTWAFF::operator()(const Tspectrum::Tcoc& s).
*
* \sa \ref page_fftw3
*/
class DRFFTWAFF {
......@@ -98,12 +117,12 @@ namespace fourier {
typedef aff::Series<Tsample> Tseries;
typedef aff::Series<Tcoeff> Tspectrum;
#ifdef FFTWFALLBACK
DRFFTWAFF(const int& n):
DRFFTWAFF(const unsigned int& n):
Msize(n), Mplan_forward(0), Mplan_backward(0) { }
DRFFTWAFF():
Msize(0), Mplan_forward(0), Mplan_backward(0) { }
#else
DRFFTWAFF(const int& n, const bool& deletearrays=false):
DRFFTWAFF(const unsigned int& n, const bool& deletearrays=false):
Msize(n), Mplan_forward(0), Mplan_backward(0),
Mseriesarray(0), Mspectrumarray(0),
Mdeletearrays(deletearrays) { }
......@@ -125,7 +144,8 @@ namespace fourier {
const bool& debug=false) const;
Tsample scale_series(const Tsample& dt) const;
Tsample scale_spectrum(const Tsample& dt) const;
int size() const { return(Msize); }
unsigned int size() const { return(Msize); }
void size(const unsigned int& s) const { this->set_size(s); }
private:
void create_plan_forward() const;
void create_plan_backward() const;
......@@ -133,10 +153,10 @@ namespace fourier {
#ifndef FFTWFALLBACK
void create_arrays() const;
void delete_arrays() const;
int ssize() const { return(this->size()/2+1); }
unsigned int ssize() const { return(this->size()/2+1); }
#endif
void set_size(const int& n) const;
mutable int Msize;
void set_size(const unsigned int& n) const;
mutable unsigned int Msize;
#ifdef FFTWFALLBACK
mutable rfftw_plan Mplan_forward;
mutable rfftw_plan Mplan_backward;
......
......@@ -215,11 +215,11 @@ namespace fourier {
double* out=Cseriesarray.castedpointer<double>();
// create plan
Mplanr2c=fftw_plan_many_dft_c2r(rank, &n, howmany,
Mplanc2r=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");
FOURIER_assert(Mplanc2r!=0, "ERROR: creating c2r plan");
}
} // void DRFFTWAFFArrayEngine::createplanc2r()
......
......@@ -76,6 +76,11 @@ namespace fourier {
* the signal index.
* This is due to the column major arrangement used by aff::Array
* containers by default (this is the common Fortran layout).
*
* \note
* The c2r transform destroys its input data.
* See
* http://fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data
*/
class DRFFTWAFFArrayEngine {
public:
......
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