autocorrelation.hpp

Go to the documentation of this file.

/** \file autocorrelation.hpp
 * \brief Tools for working with autocorrelations
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup signal_processing_files
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017 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 autocorrelation_hpp
#define autocorrelation_hpp
 
#include "../math/ft/fftT.hpp"
 
namespace mx
{
 
namespace sigproc
{
 
/// Calculate the autocorrelation of a time-series
/**
 * \tparam T is the real type of the data and autocorrelation.
 *
 * \ingroup signal_processing
 */
template <typename T>
void autocorrelation(
    T *ac,      ///< [out] is the pre-allocated array of length Nac which will contain the autocorrelation on return
    size_t Nac, ///< [in] is the length of ac
    T *sig,     ///< [in] is the input time-series (signal)
    size_t Nsig ///< [in] is the length of the input time-series
)
{
 
    #pragma omp parallel for
    for( int i = 0; i < Nac; ++i )
    {
        ac[i] = 0;
        for( int j = 0; j < Nsig - i; ++j )
        {
            ac[i] += sig[j] * sig[j + i];
        }
    }
 
    T norm = ac[0];
    for( int i = 0; i < Nac; ++i )
    {
        ac[i] /= norm;
    }
}
 
/// Calculate the autocorrelation of a time-series
/**
 * \tparam T is the real type of the data and autocorrelation.
 *
 * \todo this should probably re-allocate if ac.size() != sig.size(), or at least <
 *
 * \ingroup signal_processing
 */
template <typename T>
void autocorrelation( std::vector<T> &ac, ///< [out] will contain the autocorrelation on return.  If ac.size()==0 then
                                          ///< it is resized to sig.size().
                      std::vector<T> &sig ///< [in] is the input time-series (signal)
)
{
    if( ac.size() == 0 )
    {
        ac.resize( sig.size() );
    }
 
    autocorrelation( ac.data(), ac.size(), sig.data(), sig.size() );
}
 
/// Functor for calculating the autocorrelation given a PSD
/** Stores the fftT object and related working memory so that
 * repeated calls do not re-allocate or re-plan the FFT.
 *
 * \tparam T is the real type of the PSD and resultig A.C.
 *
 * \ingroup signal_processing
 */
template <typename T>
struct autocorrelationFromPSD
{
    std::vector<std::complex<T>> fftOut;
    std::vector<std::complex<T>> fftIn;
 
    math::ft::fftT<std::complex<T>, std::complex<T>, 1, 0> fft;
 
    /// Calculate the A.C. as the inverse FFT of the PSD
    /** This calculates the circular autocorrelation from the PSD.
     *
     */
    void
    operator()( T *ac, ///<  [out] pre-allocated array, on output contains the first Nac points of the autocorrelation
                size_t Nac, ///< [in] the allocated size of ac.
                T *psd,     ///< [in] the 2-sided FFT storage order PSD
                size_t Npsd ///< [in] the number of points in the PSD
    )
    {
        fft.plan( Npsd, math::ft::dir::forward );
 
        fftOut.resize( Npsd );
        fftIn.resize( Npsd );
 
        for( size_t i = 0; i < Npsd; ++i )
            fftIn[i] = psd[i];
 
        fft( fftOut.data(), fftIn.data() );
 
        T norm = fftOut[0].real();
        for( size_t i = 0; i < Npsd && i < Nac; ++i )
            ac[i] = fftOut[i].real() / norm;
    }
 
    /// Calculate the A.C. as the inverse FFT of the PSD
    /** This calculates the circular autocorrelation from the PSD.
     *
     */
    void operator()( std::vector<T> &ac, ///< [out] On output contains the autocorrelation.  If ac.size()==0, it is
                                         ///< allocated to psd.size().
                     std::vector<T> &psd ///<  [in] the 2-sided FFT storage order PSD
    )
    {
        if( ac.size() == 0 )
            ac.resize( psd.size() );
        operator()( ac.data(), ac.size(), psd.data(), psd.size() );
    }
};
 
} // namespace sigproc
} // namespace mx
 
#endif // autocorrelation_hpp