Go to the documentation of this file.00001 #include "ForIA/AnalysisTools/FFT.hh"
00002 #include "ForIA/Utils.hh"
00003
00004 #include <gsl/gsl_fft_real.h>
00005 #include <gsl/gsl_errno.h>
00006 #include <gsl/gsl_fft_halfcomplex.h>
00007
00008 #include "boost/lexical_cast.hpp"
00009
00010 #include <string>
00011 #include <stdexcept>
00012 #include <iostream>
00013
00014 namespace ForIA{
00015
00016 FFT::FFT(): m_freshMagnitudes(false),
00017 m_freshReals(false), m_freshImgs(false){
00018
00019 }
00020
00021 bool FFT::setGrid(const vector<double> &grid){
00022
00023 unsigned int nSegments = grid.size();
00024
00025 if( ((nSegments-1) & nSegments)){
00026 throw std::runtime_error("Number of segments must be divisible by 2. Number given is " + boost::lexical_cast<std::string>(nSegments));
00027 return false;
00028 }
00029
00030 vector<double> coeffs(grid.begin(), grid.end());
00031
00032 int fftStatus = gsl_fft_real_radix2_transform(&(coeffs[0]), 1, nSegments);
00033
00034 if(fftStatus != GSL_SUCCESS){
00035 throw std::runtime_error("Could not run gsl_fft_real_radix2_transform - probably because the input number of data bins is not a power of 2");
00036 return false;
00037 }
00038
00039 m_coefficients.clear();
00040 m_coefficients.push_back(Coefficient(coeffs[0], 0.0));
00041
00042 for(unsigned int n=1; n != nSegments / 2; ++n){
00043 m_coefficients.push_back(Coefficient(coeffs[n], coeffs[nSegments-n]));
00044 }
00045
00046 m_coefficients.push_back(Coefficient(coeffs[nSegments/2], 0.));
00047
00048 m_freshMagnitudes = true;
00049 m_freshReals = true;
00050 m_freshImgs = true;
00051
00052 return true;
00053 }
00054
00055 const FTCoefficients &FFT::coefficients() const{
00056 return m_coefficients;
00057 }
00058
00059 const vector<double> &FFT::magnitudes() const{
00060 if(m_freshMagnitudes){
00061 m_magnitudes.clear();
00062 for(FTCoefficients::const_iterator coeff = m_coefficients.begin();
00063 coeff != m_coefficients.end(); ++coeff){
00064 m_magnitudes.push_back(abs(*coeff));
00065 }
00066 m_freshMagnitudes = false;
00067 }
00068 return m_magnitudes;
00069 }
00070
00071 const vector<double> &FFT::reals() const{
00072 if(m_freshReals){
00073 m_reals.clear();
00074 for(FTCoefficients::const_iterator coeff = m_coefficients.begin();
00075 coeff != m_coefficients.end(); ++coeff){
00076 m_reals.push_back(coeff->real());
00077 }
00078 m_freshReals = false;
00079 }
00080 return m_reals;
00081 }
00082
00083 const vector<double> &FFT::imaginaries() const{
00084 if(m_freshImgs){
00085 m_imgs.clear();
00086 for(FTCoefficients::const_iterator coeff = m_coefficients.begin();
00087 coeff != m_coefficients.end(); ++coeff){
00088 m_imgs.push_back(coeff->imag());
00089 }
00090 m_freshImgs = false;
00091 }
00092 return m_imgs;
00093 }
00094
00095 }
