averagePeriodogram.hpp

Go to the documentation of this file.

/** \file averagePeriodogram.hpp
 * \brief A class to manage calculation of periodograms from time series data.
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup signal_processing_files
 *
 */
 
//***********************************************************************//
// Copyright 2020 Jared R. Males (jaredmales@gmail.com)
//
// This file is part of mxlib.
//
// mxlib 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.
//
// mxlib 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 mxlib.  If not, see <http://www.gnu.org/licenses/>.
//***********************************************************************//
 
#ifndef averagePeriodogram_hpp
#define averagePeriodogram_hpp
 
#include "psdUtils.hpp"
 
namespace mx
{
namespace sigproc
{
 
/// Calculate the average periodogram of a time-series.
/**
  * Implements the overlapped average periodogram, cf. pp 843 to 845 of \cite 10.5555.1795494. A.k.a. Welch's method
  * (https://en.wikipedia.org/wiki/Welch's_method).
  *
  * This produces
  * the variance normalized 1-sided PSD estimate of the periodogram.
  *
  * Optionally includes a window (by default no window is used, which is equivalent to the rectangle window).
  *
  * Optionally includes zero-padding, which is applied after windowing.  This can be used to increase frequency
  * resolution while still allowing overlapping.
  *
  * Can also be used to calculate the unaveraged non-overlapped periodogram of the entire time series, with or without
  * windowing and zero-padding.
  *
  * Example:
  \code
   typedef float realT;
 
   mx::fftwEnvironment<realT> fftwEnv; //for efficient fft planning
 
   std::vector<realT> ts = getTimeSeries(); //Some function to populate ts.
 
   math::vectorMeanSub(ts); //The time series should be overall mean subtracted.
 
   realT dt = 0.1; //The sampling of ts is 0.1 seconds.
 
   averagePeriodogram<realT> avgPgram(2.0/dt, 4.0/dt, dt); //This sets the averaging length to 2 seconds (20 samples),
  zero padding to 4 seconds, with a 1 second overlap (set automatically)
 
   avgPgram.win(window::hann); //Set the Hann window
 
   std::vector<realT> pgram = avgPgram(ts); //Calculate the half-overlapped periodogram estimate for ts.
 
   //Now output the periodogram estimate.  This will give two columns, the first frequency, the second power at that
  frequency. for(size_t n=0; n< avgPgram.size(); ++n)
   {
      std::cout << avgPgram[n] << " " << pgram[n] << "\n";
   }
 
  \endcode
  *
  * \ingroup psds
  */
template <typename realT>
class averagePeriodogram
{
 
  protected:
    /** \name Configuration Parameters
     *
     *@{
     */
 
    size_t m_avgLen{ 0 };  ///< The number of samples in each periodogram estimate.
    size_t m_overlap{ 0 }; ///< The number of samples by which to overlap.  This should almost always be 0.5*m_avgLen.
                           ///< Set 0 for the non-overlapped case.
 
    size_t m_padLen{ 0 }; ///< The size of the periodogram estimate with zero-padding
 
  protected:
    realT m_dt{ 1 }; ///< The time sampling.  Only used for normalization and calculation of the frequency scale.
 
    std::vector<realT> m_win; ///< The window function.  By default this is empty, which is equivalent to setting it to
                              ///< the rectangular window.
 
    ///@}
 
    size_t m_size{ 0 }; ///< The size of the periodogram vector, calculated as m_avgLen/2 + 1;
 
    int m_nOver{ 0 }; ///< The number of overlapping segments.  Calculated from m_avgLen and m_overlap;
 
    realT m_df{ 1 }; ///< The frequency sampling.  This is used only for normalization and frequency scale output.
 
    math::ft::fftT<realT, std::complex<realT>, 1, 0> m_fft;
 
    realT *m_tsWork{ nullptr };
 
    std::complex<realT> *m_fftWork{ nullptr };
    size_t m_fftWorkSize{ 0 };
 
  public:
    /// C'tor which sets up the optimum overlapped periodogram of the timeseries.
    /** Sets m_overlap = 0.5*m_avgLen.  If you desire the non-overlapped case use
      * the alternate constructor:
      * \code
        averagePeriogram p( avgLen, 0, dt);
        \endcode
      */
    explicit averagePeriodogram( size_t avgLen /**< [in] The length of averaging in samples.*/ );
 
    /// C'tor which sets up the optimum overlapped periodogram of the timeseries with zero-padding
    /** Sets m_overlap = 0.5*m_avgLen.  If you desire the non-overlapped case use
      * the alternate constructor:
      * \code
        averagePeriogram p( avgLen, 0, padLen, dt);
        \endcode
      */
    explicit averagePeriodogram(
        size_t avgLen, ///< [in] The length of averaging in samples.
        size_t padLen  ///< [in] [optional] the length of the padded sample.  If 0 or omitted then no padding used.
    );
 
    /// C'tor which sets up the optimum overlapped periodogram of the timeseries with zero-padding and sets the
    /// sampling.
    /** Sets m_overlap = 0.5*m_avgLen.  If you desire the non-overlapped case use
      * the alternate constructor:
      * \code
        averagePeriogram p( avgLen, 0, dt);
        \endcode
      */
    averagePeriodogram( size_t avgLen, ///< [in] The length of averaging in samples.
                        size_t padLen, ///< [in] the length of the padded sample.  Set to 0 for no padding.
                        realT dt       ///< [in] the sampling interval of the time-series
    );
 
    /// C'tor setting up an arbitrary overlap.
    /** Set olap to 0 for the unoverlapped case.
     */
    averagePeriodogram( size_t avgLen, ///< [in] The number of samples in each periodgoram estimate.
                        size_t padLen, ///< [in] the length of the padded sample.  Set to 0 for no padding.
                        size_t olap, ///< [in] The number of samples by which to overlap.  This should almost always be
                                     ///< 0.5*m_avgLen.  Set 0 for the non-overlapped case.
                        realT dt     ///< [in] the sampling interval of the time-series
    );
 
    /// D'tor, frees all working memory.
    ~averagePeriodogram();
 
    /// Set the sampling interval of the time-series
    /** This also sets the frequency scale of the output.
     */
    void dt( realT ndt );
 
    /// Get the sampling interval of the time-series
    /** \returns m_dt
     */
    realT dt();
 
    /// Get the frequency interval of the periodogram
    /** \returns m_df
     */
    realT df();
 
    /// Resize the periodogram working memory, setting up the 1/2-overlapped optimum case.
    /** This sets the overlap to 0.5*avgLen.
     *
     * Also performs fft planning.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int resize( size_t avgLen /**< [in] The number of samples in each periodgoram estimate. */ );
 
    /// Resize the periodogram working memory.
    /** Also performs fft planning.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int resize( size_t avgLen, ///< [in] The number of samples in each periodgoram estimate.
                size_t padLen, ///< [in] the length of the sample after zero-padding.
                size_t olap    ///< [in] The number of samples by which to overlap.  This should almost always be
                               ///< 0.5*m_avgLen.  Set 0 for the non-overlapped case.
    );
 
    /// Get a reference to the window vector.
    /** This allows population of the window with any arbitrary function.  You should call
     * resizeWin() before using this to make sure that the vector has the correct length.
     *
     * \returns a reference to m_win.
     */
    std::vector<realT> &win();
 
    /// Set the window vector using a function.
    /** This will resize m_win, and then call the function with m_win as the argument.
      * For example, to use the Hann window from the \ref signal_windows1D "1-D window functions":
      * \code
        avgPgram.win( window::hann );
        \endcode
      * which will set the Hann window.
      */
    void win( void(*winFunc)(std::vector<realT> &) /**<[in] pointer to a function which takes a pre-allocated vector and fills it in with a window function*/);
 
    /// Resize the window and, unless not desired, initialize to the rectangular window by setting all 1s.
    /** If `setRect=false` then the window is not initialized.
     */
    void resizeWin( bool setRect = true /**< [in] [optional] if false, then the window is not initialized to 1 */ );
 
    /// Calculate the periodogram for a time-series.
    /**
     *
     */
    void operator()( realT *pgram,    ///< [out] a pre-allocated to size() array which will be populated with the
                                      ///< periodogram of the time-series
                     const realT *ts, ///< [in] the time-series
                     size_t sz        ///< [in] the length of the time-series array
    );
 
    /// Calculate the periodogram for a time-series.
    /** \overload
     *
     */
    void operator()( std::vector<realT> &pgram,   ///< [out] a vector which will be allocated and populated with the
                                                  ///< periodogram of the time-series
                     const std::vector<realT> &ts ///< [in] the time-series
    );
 
    /// Calculate the periodogram for a time-series.
    /** \overload
     *
     * \returns the periodogram as a vector.
     */
    std::vector<realT> operator()( std::vector<realT> &ts /**< [in] the time-series*/ );
 
    /// Return the size of the periodogram.
    /**
     *
     * \returns the size of the periodogram estimate for the current setup.
     */
    size_t size();
 
    /// Get the frequency at a given point in the periodogram.
    /**
     * \returns the frequency at point n in the periodogram.
     */
    realT operator[]( size_t n /**<[in] the point in the periodogram at which the frequency is desired*/ );
};
 
template <typename realT>
averagePeriodogram<realT>::averagePeriodogram( size_t avgLen )
{
    resize( avgLen, 0, 0.5 * avgLen );
    dt( 1 );
}
 
template <typename realT>
averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen )
{
    resize( avgLen, padLen, 0.5 * avgLen );
    dt( 1 );
}
 
template <typename realT>
averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen, realT ndt )
{
    resize( avgLen, padLen, 0.5 * avgLen );
    dt( ndt );
}
 
template <typename realT>
averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen, size_t olap, realT ndt )
{
    resize( avgLen, padLen, olap );
    dt( ndt );
}
 
template <typename realT>
averagePeriodogram<realT>::~averagePeriodogram()
{
    if( m_tsWork )
    {
        fftw_free( m_tsWork );
    }
 
    if( m_fftWork )
    {
        fftw_free( m_fftWork );
    }
}
 
template <typename realT>
void averagePeriodogram<realT>::dt( realT ndt )
{
    m_dt = ndt;
    m_df = 1.0 / ( m_padLen * m_dt );
}
 
template <typename realT>
realT averagePeriodogram<realT>::dt()
{
    return m_dt;
}
 
template <typename realT>
realT averagePeriodogram<realT>::df()
{
    return m_df;
}
 
template <typename realT>
int averagePeriodogram<realT>::resize( size_t avgLen )
{
    return resize( avgLen, 0, 0.5 * avgLen );
}
 
template <typename realT>
int averagePeriodogram<realT>::resize( size_t avgLen, size_t padLen, size_t olap )
{
    m_avgLen = avgLen;
    m_padLen = padLen;
 
    if( m_padLen < m_avgLen )
    {
        m_padLen = m_avgLen;
    }
 
    m_size = m_padLen / 2 + 1;
 
    m_df = 1.0 / ( m_padLen * m_dt );
 
    m_overlap = olap;
 
    m_nOver = ( m_avgLen - m_overlap );
 
    m_fft.plan( m_padLen, math::ft::dir::forward, false );
 
    if( m_tsWork )
    {
        fftw_free( m_tsWork );
    }
 
    m_tsWork = math::ft::fftw_malloc<realT>( m_padLen );
 
    if( m_fftWork )
    {
        fftw_free( m_fftWork );
    }
 
    m_fftWork = math::ft::fftw_malloc<std::complex<realT>>( ( m_padLen / 2 + 1 ) );
 
    return 0;
}
 
template <typename realT>
std::vector<realT> &averagePeriodogram<realT>::win()
{
    return m_win;
}
 
template <typename realT>
void averagePeriodogram<realT>::win( void ( *winFunc )( std::vector<realT> & ) )
{
    resizeWin( false );
    winFunc( m_win );
}
 
template <typename realT>
void averagePeriodogram<realT>::resizeWin( bool setRect )
{
    if( setRect )
    {
        m_win.resize( m_avgLen, static_cast<realT>( 1 ) );
    }
    else
    {
        m_win.resize( m_avgLen );
    }
}
 
template <typename realT>
void averagePeriodogram<realT>::operator()( realT *pgram, const realT *ts, size_t sz )
{
    if( m_win.size() > 0 && m_win.size() != m_avgLen )
    {
        std::cerr << "averagePeriodogram: Window size not correct.\n";
    }
 
    int Navg = sz / m_nOver;
 
    while( Navg * m_nOver + m_avgLen > sz )
        --Navg;
 
    if( Navg < 1 )
        Navg = 1; // Always do at least 1!
 
    for( int i = 0; i < Navg; ++i )
    {
        if( m_win.size() == m_avgLen ) // no if inside the for loop
        {
            for( size_t j = 0; j < m_avgLen; ++j )
            {
                m_tsWork[j] = ts[i * m_nOver + j] * m_win[j];
            }
        }
        else
        {
            for( size_t j = 0; j < m_avgLen; ++j )
            {
                m_tsWork[j] = ts[i * m_nOver + j];
            }
        }
 
        // Zero-pad
        for( size_t j = m_avgLen; j < m_padLen; ++j )
        {
            m_tsWork[j] = 0;
        }
 
        m_fft( m_fftWork, m_tsWork );
 
        for( size_t j = 0; j < m_size; ++j )
            pgram[j] += norm( m_fftWork[j] );
    }
 
    realT pgramVar = psdVar1sided( m_df, pgram, m_size );
 
    realT tsVar = mx::math::vectorVariance( ts, sz );
 
    for( size_t j = 0; j < m_size; ++j )
        pgram[j] *= tsVar / pgramVar;
}
 
template <typename realT>
void averagePeriodogram<realT>::operator()( std::vector<realT> &pgram, const std::vector<realT> &ts )
{
    pgram.resize( m_size );
    operator()( pgram.data(), ts.data(), ts.size() );
}
 
template <typename realT>
std::vector<realT> averagePeriodogram<realT>::operator()( std::vector<realT> &ts )
{
    std::vector<realT> pgram;
    operator()( pgram, ts );
    return pgram;
}
 
template <typename realT>
size_t averagePeriodogram<realT>::size()
{
    return m_size;
}
 
template <typename realT>
realT averagePeriodogram<realT>::operator[]( size_t n )
{
    return n * m_df;
}
 
} // namespace sigproc
} // namespace mx
 
#endif // averagePeriodogram_hpp