psdFilterCuda.hpp

Go to the documentation of this file.

/** \file psdFilterCuda.hpp
 * \brief Declares and defines a class for filtering with PSDs on CUDA GPUs
 * \ingroup signal_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2018 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 psdFilterCuda_hpp
#define psdFilterCuda_hpp
 
#include "../cuda/templateCudaPtr.hpp"
#include "../cuda/templateCuda.hpp"
 
#include "psdFilter.hpp"
 
#include <cufft.h>
#include <cufftXt.h>
// #include <helper_functions.h>
#include <helper_cuda.h>
 
namespace mx
{
namespace sigproc
{
 
/// A class for filtering noise with PSDs using CUDA cuFFT
/** The square-root of the PSD is maintained by this class, either as a pointer to an external array or using internally
  allocated memory (which will be
  * de-allocated on destruction.
  *
  * PSD Normalization: the PSD used for this needs to be normalized properly to produce filtered noise with the correct
  statistics.  Given an
  * array of size N which contains the PSD as "power per unit frequency" vs frequency (i.e. the PSD) on some frequency
  scale with uniform spacing \f$ \Delta f \f$,
  * the normalization is
  \f[
       PSD_{norm} = PSD * N * \Delta f
  \f]
  * for a 1-dimensional PSD and
  \f[
       PSD_{norm} = PSD * N_0 * \Delta f_0 * N_1 * \Delta f_1
  \f]
  * for a 2-dimensional PSD. Remember that these are applied before taking the square root.
  *
  *
  * \ingroup psd_filter
  *
  * \todo once fftT has a plan interface with pointers for working memory, use it.
  */
template <typename _realT>
class psdFilter<_realT, 1>
{
  public:
    typedef _realT realT;                                                   ///< Real floating point type
    typedef std::complex<_realT> complexT;                                  ///< Complex floating point type.
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> realArrayT; ///< Eigen array type with Scalar==realT
    typedef Eigen::Array<complexT, Eigen::Dynamic, Eigen::Dynamic>
        complexArrayT; ///< Eigen array type with Scalar==complexT
 
    typedef realT deviceRealPtrT;
    typedef complexT deviceComplexPtrT;
 
  protected:
    int m_rows{ 0 }; ///< The number of rows in the filter, and the required number of rows in the noise array.
    int m_cols{ 0 }; ///< The number of columns in the filter, and the required number of columns in the noise array.
 
    mx::cuda::cudaPtr<complexT> m_psdSqrt; ///< Pointer to the real array containing the square root of the PSD.
 
    mx::cuda::cudaPtr<complexT> m_scale; ///< the scale factor.
 
    /// Cuda FFT plan.  We only need one since the forward/inverse is part of execution.
    cufftHandle m_fftPlan{ 0 };
 
  public:
    /// C'tor.
    psdFilter();
 
    /// Destructor
    ~psdFilter();
 
  protected:
    /// Set the size of the filter.
    /** Handles allocation of the _ftWork array and fftw planning.
     *
     * Requires _psdSqrt to be set first.  This is called by the psdSqrt() and psd() methods.
     *
     */
    int setSize();
 
  public:
    /// Get the number of rows in the filter
    /**
     * \returns the current value of m_rows.
     */
    int rows();
 
    /// Get the number of columns in the filter
    /**
     * \returns the current value of m_cols.
     */
    int cols();
 
    /// Set the sqaure-root of the PSD.
    /** This allocates _npsdSqrt and fills it with th evalues in the array.
     *
     * See the discussion of PSD normalization above.
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int psdSqrt( const realArrayT &npsdSqrt /**< [in] an array containing the square root of the PSD.*/ );
 
    int psdSqrt(
        const cuda::cudaPtr<std::complex<realT>> &npsdSqrt, /**< [in] an array containing the square root of the PSD.*/
        size_t rows,
        size_t cols );
 
    /// De-allocate all working memory and reset to initial state.
    void clear();
 
    /// Apply the filter.
    /**
     *
     * \returns 0 on success
     * \returns -1 on error
     */
    int
    filter( deviceComplexPtrT *noise ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
    );
 
    /// Apply the filter.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int operator()(
        realArrayT &noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ );
 
    /// Apply the filter.
    /**
     * \returns 0 on success
     * \returns -1 on error
     */
    int
    operator()( realArrayT &noise,  ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
                realArrayT &noiseIm ///< [out] [optional] an array to fill with the imaginary output of the filter,
                                    ///< allowing 2-for-1 calculation.
    );
};
 
template <typename realT>
psdFilter<realT, 1>::psdFilter()
{
}
 
template <typename realT>
psdFilter<realT, 1>::~psdFilter()
{
    if( m_fftPlan )
    {
        checkCudaErrors( cufftDestroy( m_fftPlan ) );
    }
}
 
template <typename realT>
int psdFilter<realT, 1>::setSize()
{
 
    m_scale.resize( 1 );
    std::complex<realT> scale( 1. / ( m_rows * m_cols ), 0 );
    m_scale.upload( &scale, 1 );
 
    checkCudaErrors( cufftPlan2d( &m_fftPlan, m_rows, m_cols, CUFFT_C2C ) );
 
    return 0;
}
 
template <typename realT>
int psdFilter<realT, 1>::rows()
{
    return m_rows;
}
 
template <typename realT>
int psdFilter<realT, 1>::cols()
{
    return m_cols;
}
 
// template<typename realT>
// int psdFilter<realT,1>::psdSqrt( realArrayT * npsdSqrt )
// {
//    if(_psdSqrt && _owner)
//    {
//       delete _psdSqrt;
//    }
//
//    _psdSqrt = npsdSqrt;
//    _owner = false;
//
//    setSize();
//
//    return 0;
// }
 
template <typename realT>
int psdFilter<realT, 1>::psdSqrt( const realArrayT &npsdSqrt )
{
    m_rows = npsdSqrt.rows();
    m_cols = npsdSqrt.cols();
 
    complexArrayT tmp( m_rows, m_cols );
    tmp.setZero();
    tmp.real() = npsdSqrt;
 
    m_psdSqrt.upload( tmp.data(), m_rows * m_cols );
 
    setSize();
 
    return 0;
}
 
template <typename realT>
int psdFilter<realT, 1>::psdSqrt( const cuda::cudaPtr<std::complex<realT>> &npsdSqrt, size_t rows, size_t cols )
{
    m_rows = rows;
    m_cols = cols;
 
    m_psdSqrt.resize( m_rows * m_cols );
 
    ///\todo move this into cudaPtr
    cudaMemcpy( m_psdSqrt.m_devicePtr,
                npsdSqrt.m_devicePtr,
                m_rows * m_cols * sizeof( std::complex<realT> ),
                cudaMemcpyDeviceToDevice );
 
    setSize();
 
    return 0;
}
 
template <typename realT>
void psdFilter<realT, 1>::clear()
{
    m_scale.resize( 0 );
    m_rows = 0;
    m_cols = 0;
 
    m_psdSqrt.resize( 0 );
}
 
template <typename realT>
int psdFilter<realT, 1>::filter( deviceComplexPtrT *noise )
{
 
    // Transform complex noise to Fourier domain.
 
    cufftExecC2C( m_fftPlan, (cuComplex *)noise, (cuComplex *)noise, CUFFT_FORWARD );
 
    // Apply the filter.
    mx::cuda::pointwiseMul<cuComplex>
        <<<32, 256>>>( (cuComplex *)noise, (cufftComplex *)m_psdSqrt.m_devicePtr, m_rows * m_cols );
 
    cufftExecC2C( m_fftPlan, (cuComplex *)noise, (cuComplex *)noise, CUFFT_INVERSE );
 
    mx::cuda::scalarMul<cuComplex>
        <<<32, 256>>>( (cuComplex *)noise, (cuComplex *)m_scale.m_devicePtr, m_rows * m_cols );
    // Now take the real part, and normalize.
 
    // noise = noise/(noise.rows()*noise.cols());
 
    return 0;
}
 
template <typename realT>
int psdFilter<realT, 1>::operator()( realArrayT &noise )
{
    return filter( noise );
}
 
template <typename realT>
int psdFilter<realT, 1>::operator()( realArrayT &noise, realArrayT &noiseIm )
{
    return filter( noise, &noiseIm );
}
 
} // namespace sigproc
} // namespace mx
 
#endif // psdFilterCuda_hpp