psdVarMean.hpp

Go to the documentation of this file.

/** \file psdVarMean.hpp
 * \brief Tools for calculating the variance of the mean of a PSD
 *
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 * \ingroup signal_processing_files
 *
 */
 
//***********************************************************************//
// Copyright 2018,2019,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 psdVarMean_hpp
#define psdVarMean_hpp
 
#include <vector>
 
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
 
#include "../math/func/bessel.hpp"
#include "../math/constants.hpp"
#include "../math/gslInterpolator.hpp"
 
#include "psdFilter.hpp"
 
namespace mx
{
namespace sigproc
{
 
/// Parameters for calculating the variance of the mean from a numerical PSD.
/** Interpolates on the PSD for the integration.
 *
 *  \tparam realT the real floating point type of the frequency and PSD vectors. [Must be double due to GSL]
 */
template <typename _realT>
struct psdVarMeanParams
{
    typedef _realT realT;
 
    std::vector<realT> *freq{ nullptr }; ///< Pointer to the frequency vector
    std::vector<realT> *psd{ nullptr };  ///< Pointer to the psd vector
    realT T{ 0 };                        ///< The sample length (units appropriate for freq)
 
    realT minf{ 0 }; ///< The minimum frequency
    realT maxf{ 0 }; ///< The maximum frequency
 
    realT extrapAlpha{ 0 };
 
    bool initialized{ false }; ///< Flag controlling whether the interpolator is initialized
 
#ifdef MX_OLD_GSL
    math::gslInterpolator<math::gsl_interp_linear<realT>> terp; ///< Interpolator for the PSD
#else
    math::gslInterpolator<math::gsl_interp_steffen<realT>> terp; ///< Interpolator for the PSD
#endif
 
    /// Get the value of the PSD at k using interpolation
    /** Extrapolates outside the range of frequencies in freq.
     * Right now this is just constant.
     *
     * \todo implement better extrapolation past max-freq.
     */
    realT psdVal( realT k /**< [in] The transformed frequency coordinate */ )
    {
        if( !initialized )
        {
            if( freq == nullptr || psd == nullptr )
            {
                std::cerr << "Vectors not set.\n";
                exit( -1 );
            }
            if( T == 0 )
            {
                std::cerr << "T not set.\n";
                exit( -1 );
            }
 
            terp.setup( *freq, *psd );
            // #ifdef MX_OLD_GSL
            // terp.setup(gsl_interp_linear, *freq, *psd);
            // #else
            // terp.setup(gsl_interp_steffen, *freq, *psd);
            // #endif
 
            minf = ( *freq )[0];
            maxf = ( *freq )[freq->size() - 1];
 
            initialized = true;
        }
 
        realT f = ( 2.0 / T ) * k;
 
        if( f < minf )
            return ( *psd )[0];
 
        if( f > maxf )
        {
            // Enter extrapolation here
            // return 0;//(*psd)[psd->size()-1];
            if( extrapAlpha == 0 )
                return 0;
 
            return ( *psd )[psd->size() - 1] * pow( f / maxf, -1 * extrapAlpha );
        }
 
        return terp( f );
    }
    realT psdVal( realT k /**< [in] The transformed frequency coordinate */ ) {…}
};
struct psdVarMeanParams {…};
 
/// Integration worker function for psdVarMean::operator()
/**
 * \tparam paramsT the type of the parameter structure. See psdVarMean.
 *
 * \returns the value of the integrand at k
 *
 *
 */
template <typename paramsT>
typename paramsT::realT psdVarMeanFunc( typename paramsT::realT k, ///< [in] the scaled frequency coordinate
                                        void *params               ///< [in] pointer to the parameter structure
)
{
    typedef typename paramsT::realT realT;
 
    paramsT *pvar = (paramsT *)params;
 
    realT J = math::func::bessel_j<realT, realT>( 0.5, math::two_pi<realT>() * k );
 
    return J * J / k * pvar->psdVal( k );
}
typename paramsT::realT psdVarMeanFunc( typename paramsT::realT k, ///< [in] the scaled frequency coordinate {…}
 
/// Calculate the variance of the mean for a process given its PSD
/** Functor which manages the GSL integration workspace memory.
 *
 * \todo document the paramsT interface
 *
 * \tparam paramsT is a parameter structure with the correct interface.
 *
 *  \ingroup psds
 */
template <typename paramsT>
struct psdVarMean
{
    typedef typename paramsT::realT realT;
 
    realT extrapAlpha{ 0 }; ///< The extrapolation exponent.  No extrapolation if 0.
 
  protected:
    gsl_integration_workspace *w{ nullptr }; ///< The GSL integration workspace memory
 
    size_t wSize{ 10000 }; ///< The size of the GSL workspace
 
    /// Initializes the gsl library.
    void init()
    {
        ///\todo We'll get the occasional failure to reach tolerance error, just ignore them all for now.
        gsl_set_error_handler_off();
    }
    void init() {…}
 
  public:
    /// Default c'tor
    psdVarMean()
    {
        init();
    }
    psdVarMean() {…}
 
    /// Constructor which initializes workspace size
    psdVarMean( size_t wm /**< [in] Size of the GSL integration workspace*/ )
    {
        init();
        wSize = wm;
    }
    psdVarMean( size_t wm /**< [in] Size of the GSL integration workspace*/ ) {…}
 
    /// d'tor cleans up the workspace memory
    ~psdVarMean()
    {
        if( w != nullptr )
        {
            gsl_integration_workspace_free( w );
        }
    }
    ~psdVarMean() {…}
 
    /// Calculate the variance of the mean for a psd over a sample length
    /**
     * \returns the variance of the mean for the process governed by the PSD.
     */
    realT operator()( realT &error,             ///< [out] The error estimate returned by gsl_integration_qagiu
                      std::vector<realT> &freq, ///< [in] The frequency vector
                      std::vector<realT> &psd,  ///< [in] The psd vector
                      realT T                   ///< [in] The sample length (units approprate for freq)
    )
    {
        if( w == nullptr )
        {
            w = gsl_integration_workspace_alloc( wSize );
        }
 
        // Setup the integration parameters structure
        paramsT pvar;
 
        pvar.freq = &freq;
        pvar.psd = &psd;
        pvar.T = T;
 
        pvar.extrapAlpha = extrapAlpha;
 
        // Setup the function
        gsl_function func;
 
        func.function = &psdVarMeanFunc<paramsT>;
        func.params = &pvar;
 
        realT result;
 
        int ec = gsl_integration_qagiu( &func, 0.0, 1e-14, 1e-8, wSize, w, &result, &error );
        static_cast<void>( ec );
 
        return result / ( 2 * T );
    }
    realT operator()( realT &error,             ///< [out] The error estimate returned by gsl_integration_qagiu {…}
};
struct psdVarMean {…};
 
} // namespace sigproc
} // namespace mx
 
#endif // psdVarMean_hpp