portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2024-11-12 17:35:35 +00:00
Code Releases Activity

FFT with GSL

Browse Source
This commit is contained in:
Antonin Portelli 2017-10-10 16:56:23 +01:00
parent 839159831f
commit 6b624186cd
6 changed files with 301 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
lib/FFT.hpp Normal file
View File

@ -0,0 +1,53 @@
/*
* FFT.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
*
* LatAnalyze 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 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 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 LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_FFT_hpp_
#define Latan_FFT_hpp_
#include <LatAnalyze/Global.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* FFT abstract class *
******************************************************************************/
class FFT
{
public:
enum
{
Forward = 0,
Backward = 1
};
public:
// constructor
FFT(void) = default;
FFT(const Index size);
// destructor
virtual ~FFT(void) = default;
// size
virtual void resize(const Index size) = 0;
// FFT
virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward) = 0;
};
END_LATAN_NAMESPACE
#endif // Latan_FFT_hpp_

89
lib/GslFFT.cpp Normal file
View File

@ -0,0 +1,89 @@
/*
* GslFFT.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
*
* LatAnalyze 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 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 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 LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/GslFFT.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* GslFFT implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
GslFFT::GslFFT(const Index size)
{
resize(size);
}
// destructor //////////////////////////////////////////////////////////////////
GslFFT::~GslFFT(void)
{
clear();
}
// size ////////////////////////////////////////////////////////////////////////
void GslFFT::resize(const Index size)
{
if (size_ != size)
{
size_ = size;
wavetable_ = gsl_fft_complex_wavetable_alloc(size_);
workspace_ = gsl_fft_complex_workspace_alloc(size_);
}
}
// fft /////////////////////////////////////////////////////////////////////////
void GslFFT::operator()(CMat &x, const unsigned int dir)
{
if (x.size() != size_)
{
LATAN_ERROR(Size, "wrong input vector size");
}
else
{
switch (dir)
{
case FFT::Forward:
gsl_fft_complex_forward((double *)x.data(), 1, size_,
wavetable_, workspace_);
break;
case FFT::Backward:
gsl_fft_complex_backward((double *)x.data(), 1, size_,
wavetable_, workspace_);
break;
default:
LATAN_ERROR(Argument, "invalid FT direction");
break;
}
}
}
// destroy GSL objects /////////////////////////////////////////////////////////
void GslFFT::clear(void)
{
if (!wavetable_)
{
gsl_fft_complex_wavetable_free(wavetable_);
}
if (!workspace_)
{
gsl_fft_complex_workspace_free(workspace_);
}
}

57
lib/GslFFT.hpp Normal file
View File

@ -0,0 +1,57 @@
/*
* GslFFT.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
*
* LatAnalyze 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 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 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 LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_GslFFT_hpp_
#define Latan_GslFFT_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/FFT.hpp>
#include <gsl/gsl_fft_complex.h>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* GSL FFT *
******************************************************************************/
class GslFFT: public FFT
{
public:
// constructors
GslFFT(void) = default;
GslFFT(const Index size);
// destructor
virtual ~GslFFT(void);
// size
void resize(const Index size);
// fft
virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward);
private:
// destroy GSL objects
void clear(void);
private:
Index size_{0};
gsl_fft_complex_wavetable *wavetable_{nullptr};
gsl_fft_complex_workspace *workspace_{nullptr};
};
END_LATAN_NAMESPACE
#endif // Latan_GslFFT_hpp_

3
lib/Makefile.am
View File

@ -31,6 +31,7 @@ libLatAnalyze_la_SOURCES = \
FitInterface.cpp \ FitInterface.cpp \
Function.cpp \ Function.cpp \
Global.cpp \ Global.cpp \
GslFFT.cpp \
GslHybridRootFinder.cpp\ GslHybridRootFinder.cpp\
GslMinimizer.cpp \ GslMinimizer.cpp \
GslQagsIntegrator.cpp \ GslQagsIntegrator.cpp \
@ -59,10 +60,12 @@ libLatAnalyze_la_HEADERS = \
Dataset.hpp \ Dataset.hpp \
Derivative.hpp \ Derivative.hpp \
Exceptions.hpp \ Exceptions.hpp \
FFT.hpp \
Function.hpp \ Function.hpp \
File.hpp \ File.hpp \
FitInterface.hpp \ FitInterface.hpp \
Global.hpp \ Global.hpp \
GslFFT.hpp \
GslHybridRootFinder.hpp\ GslHybridRootFinder.hpp\
GslMinimizer.hpp \ GslMinimizer.hpp \
GslQagsIntegrator.hpp \ GslQagsIntegrator.hpp \

5
utils/Makefile.am
View File

@ -10,6 +10,7 @@ bin_PROGRAMS = \
latan-sample-combine \ latan-sample-combine \
latan-sample-element \ latan-sample-element \
latan-sample-fake \ latan-sample-fake \
latan-sample-ft \
latan-sample-plot \ latan-sample-plot \
latan-sample-plot-corr\ latan-sample-plot-corr\
latan-sample-read \ latan-sample-read \
@ -27,6 +28,10 @@ latan_sample_fake_SOURCES = sample-fake.cpp
latan_sample_fake_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_fake_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_fake_LDFLAGS = -L../lib/.libs -lLatAnalyze latan_sample_fake_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_sample_ft_SOURCES = sample-ft.cpp
latan_sample_ft_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_ft_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_sample_plot_corr_SOURCES = sample-plot-corr.cpp latan_sample_plot_corr_SOURCES = sample-plot-corr.cpp
latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_plot_corr_LDFLAGS = -L../lib/.libs -lLatAnalyze latan_sample_plot_corr_LDFLAGS = -L../lib/.libs -lLatAnalyze

94
utils/sample-ft.cpp Normal file
View File

@ -0,0 +1,94 @@
/*
* sample-ft.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
*
* LatAnalyze 3 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 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 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 LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatCore/OptParser.hpp>
#include <LatAnalyze/GslFFT.hpp>
#include <LatAnalyze/Io.hpp>
using namespace std;
using namespace Latan;
int main(int argc, char *argv[])
{
// argument parsing ////////////////////////////////////////////////////////
OptParser opt;
bool parsed;
string inFilename, outFilename;
unsigned int dir = FFT::Forward;
opt.addOption("o", "output", OptParser::OptType::value , true,
"output file name (default: result not saved)", "");
opt.addOption("b", "backward", OptParser::OptType::trigger, true,
"backward Fourier transform (forward by default)");
opt.addOption("" , "help" , OptParser::OptType::trigger, true,
"show this help message and exit");
parsed = opt.parse(argc, argv);
if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
{
cerr << "usage: " << argv[0];
cerr << " <options> <input file>" << endl;
cerr << endl << "Possible options:" << endl << opt << endl;
return EXIT_FAILURE;
}
inFilename = opt.getArgs()[0];
outFilename = opt.optionValue("o");
if (opt.gotOption("b"))
{
dir = FFT::Backward;
}
// Fourier transform ///////////////////////////////////////////////////////
DMatSample in = Io::load<DMatSample>(inFilename);
Index nSample = in.size(), l = in[central].rows();
bool isInComplex = (in[central].cols() > 1);
CMatSample res(nSample, l, 1);
DMatSample out(nSample, l, 2);
GslFFT ft(l);
cout << "-- computing Fourier transform..." << endl;
FOR_STAT_ARRAY(in, s)
{
res[s].real() = in[s].col(0);
if (isInComplex)
{
res[s].imag() = in[s].col(1);
}
else
{
res[s].imag() = DVec::Constant(l, 0.);
}
ft(res[s], dir);
out[s].col(0) = res[s].real();
out[s].col(1) = res[s].imag();
}
// output /////////////////////////////////////////////////////////////////
cout << scientific;
cout << "central value:\n" << out[central];
cout << endl;
cout << "standard deviation:\n" << out.variance().cwiseSqrt();
cout << endl;
if (!outFilename.empty())
{
Io::save<DMatSample>(out, outFilename);
}
return EXIT_SUCCESS;
}