GPLIB++
StackedSpectrum.h
Go to the documentation of this file.
1 #ifndef CSTACKEDSPECTRUM_H
2 #define CSTACKEDSPECTRUM_H
3 #include "TsSpectrum.h"
4 #include "FatalException.h"
5 #include <functional>
6 #include <iterator>
7 #include "statutils.h"
8 #include <boost/bind.hpp>
9 
10 namespace gplib
11  {
12 
13  /** \addtogroup tstools Time series analysis methods */
14  /* @{ */
15 
16  //! This template is used to calculate stacked spectra for example for power spectrum estimation
17  /*! The template StackedSpectrum calculates stacked spectra from an input time series. Two input iterators
18  * define the input range, one output iterator the start of the output range, we have to provide the length
19  * of individual segments and a windowing function. Make sure the container the output is stored in can hold at least seglength/2+1.
20  * The Windowing function gets the relative position within the segment and should therefore be defined on the interval 0..1
21  */
22  template<typename InputIterator, typename OutputIterator,
23  typename WindowFunctype>
24  void StackedSpectrum(InputIterator tsbegin, InputIterator tsend,
25  OutputIterator freqbegin, const size_t seglength, WindowFunctype WFunc)
26  {
27  //some convenience type definitions for the vectors that hold the input and output values
28  typedef std::vector<
29  typename std::iterator_traits<OutputIterator>::value_type>
30  toutvector;
31  typedef std::vector<
32  typename std::iterator_traits<InputIterator>::value_type> tinvector;
33  const size_t insize = distance(tsbegin, tsend);
34  if (insize < seglength)
35  throw FatalException(
36  "Time Series shorter than requested segment length");
37  const unsigned int nsegments = insize / seglength; //we do integer division, so the last (not full) segment will be discarded
38  TsSpectrum SpecEst(false); //we do not want to perform multiple calculations
39  std::fill_n(freqbegin, seglength / 2 + 1,
40  typename std::iterator_traits<OutputIterator>::value_type()); // fill output with 0
41  toutvector segspec(seglength / 2 + 1); //the spectrum of the current segment
42 
43  InputIterator currstart = tsbegin; //the start of the current segment
44  InputIterator currend = tsbegin;
45  advance(currend, seglength);
46  for (unsigned int i = 0; i < nsegments; ++i)
47  {
48  tinvector currseg(currstart, currend); //temporary copy of time series
49 
50  SubMean(currseg.begin(), currseg.end());
51 
52  for (typename tinvector::iterator it = currseg.begin(); it
53  != currseg.end(); ++it)
54  {
55  //multiply the current segment with the windowing function
56  *it *= WFunc(
57  static_cast<double> (distance(currseg.begin(), it))
58  / seglength); //by passing the relative position to the windowing function
59  }
60 
61  //calculate the fourier transform for the current segment
62  SpecEst.CalcSpectrum(currseg.begin(), currseg.end(),
63  segspec.begin(), segspec.end());
64  //add the spectral values to the output vector
65  std::transform(segspec.begin(), segspec.end(), freqbegin,
66  freqbegin, std::plus<typename std::iterator_traits<
67  OutputIterator>::value_type>());
68  advance(currstart, seglength); //shift the boundaries of current segment
69  advance(currend, seglength);
70  }
71  }
72  /* @} */
73  }
74 #endif // CSTACKEDSPECTRUM_H
void CalcSpectrum(_InputIterator tsbegin, _InputIterator tsend, _OutputIterator freqbegin, _OutputIterator freqend)
The member function to calculate a spectrum from a time series.
Definition: TsSpectrum.h:80
void SubMean(InputIterator begin, InputIterator end, typename std::iterator_traits< InputIterator >::value_type mean)
Substract the mean from each element in the container, mean is passed as a parameter.
Definition: statutils.h:85
The class CTsSpectrum is used to calculate spectra from (real) time series data.
Definition: TsSpectrum.h:21
void StackedSpectrum(InputIterator tsbegin, InputIterator tsend, OutputIterator freqbegin, const size_t seglength, WindowFunctype WFunc)
This template is used to calculate stacked spectra for example for power spectrum estimation...
The basic exception class for all errors that arise in gplib.