psdFilter.hpp

Go to the documentation of this file.

/** \file psdFilter.hpp
 * \brief Declares and defines a class for filtering with PSDs
 * \ingroup signal_processing_files
 * \author Jared R. Males (jaredmales@gmail.com)
 *
 */
 
//***********************************************************************//
// Copyright 2015, 2016, 2017, 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 psdFilter_hpp
#define psdFilter_hpp
 
#include <vector>
#include <complex>
#include <Eigen/Dense>
 
#include "../mxlib.hpp"
#include "../math/ft/fftT.hpp"
#include "../improc/eigenCube.hpp"
 
namespace mx
{
namespace sigproc
{
 
namespace psdFilterTypes
{
 
/// Types for different ranks in psdFilter
template <typename realT, size_t rank>
struct arrayT;
 
template <typename realT>
struct arrayT<realT, 1>
{
    typedef std::vector<realT> realArrayT;
    typedef std::vector<realT> *realArrayMapT;
    typedef std::vector<std::complex<realT>> complexArrayT;
 
    static void clear( realArrayT &arr )
    {
        arr.clear();
    }
 
    static void clear( complexArrayT &arr )
    {
        arr.clear();
    }
};
 
template <typename realT>
struct arrayT<realT, 2>
{
    typedef Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic> realArrayT;
    typedef Eigen::Map<Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic>> realArrayMapT;
 
    typedef Eigen::Array<std::complex<realT>, Eigen::Dynamic, Eigen::Dynamic> complexArrayT;
 
    static void clear( realArrayT &arr )
    {
        arr.resize( 0, 0 );
    }
 
    static void clear( complexArrayT &arr )
    {
        arr.resize( 0, 0 );
    }
};
 
template <typename realT>
struct arrayT<realT, 3>
{
    typedef improc::eigenCube<realT> realArrayT;
    typedef improc::eigenCube<realT> *realArrayMapT;
    typedef improc::eigenCube<std::complex<realT>> complexArrayT;
 
    static void clear( realArrayT &arr )
    {
        arr.resize( 0, 0, 0 );
    }
 
    static void clear( complexArrayT &arr )
    {
        arr.resize( 0, 0, 0 );
    }
};
 
} // namespace psdFilterTypes
 
// Forward declaration
template <typename _realT, size_t rank, int cuda = 0>
class psdFilter;
 
/** \defgroup psd_filter PSD Filter
 * \brief Filtering with a PSD to generate correlated noise.
 *
 * \ingroup psds
 */
 
/// A class for filtering noise with PSDs
/** 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 Requirements:
 * - the PSD must be in FFT storage order form.  That means including negative frequencies reversed from the end of the
 * array.
 * - the PSD used for this needs to be normalized properly, \ref psds "according to the mxlib standard", to produce
 * filtered noise with the correct statistics.
 *
 *
 * Array type varies based on rank.
 * - For rank==1, the array type is std::vector<realT>
 * - for rank==2, the array type is Eigen::Array<realT, -1, -1>
 * - for rank==3, the array type is mx::improc::eigenCube<realT>
 * and likewise for the complex types.
 *
 * \tparam _realT real floating type
 * \tparam _rank the rank, or dimension, of the PSD
 *
 * \ingroup psd_filter
 *
 * \todo once fftT has a plan interface with pointers for working memory, use it.
 *
 */
template <typename _realT, size_t _rank>
class psdFilter<_realT, _rank, 0>
{
  public:
    typedef _realT realT;                  ///< Real floating point type
    typedef std::complex<_realT> complexT; ///< Complex floating point type.
 
    static const size_t rank = _rank;
 
    typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayT
        realArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
    typedef typename psdFilterTypes::arrayT<realT, rank>::realArrayMapT
        realArrayMapT; ///< std::vector for rank==1, Eigen::Map for rank==2, eigenCube for rank==3.
    typedef typename psdFilterTypes::arrayT<realT, rank>::complexArrayT
        complexArrayT; ///< std::vector for rank==1, Eigen::Array for rank==2, eigenCube for rank==3.
 
  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.
    int m_planes{ 0 }; ///< Then number of planes in the filter.
 
    realT m_dFreq1{ 1.0 }; ///< The frequency scaling of the x-dimension.  Used to scale the output.
    realT m_dFreq2{ 1.0 }; ///< The frequency scaling of the y-dimension.  Used to scale the output.
    realT m_dFreq3{ 1.0 }; ///< The frequency scaling of the z-dimension.  Used to scale the output.
 
    realArrayT *m_psdSqrt{ nullptr }; ///< Pointer to the real array containing the square root of the PSD.
    bool m_owner{ false }; ///< Flag indicates whether or not m_psdSqrt was allocated by this instance, and so must be
                           ///< deallocated.
 
    mutable complexArrayT
        m_ftWork; ///< Working memory for the FFT.  Declared mutable so it can be accessed in the const filter method.
 
    math::ft::fftT<complexT, complexT, rank, 0> m_fft_fwd;  ///< FFT object for the forward transform.
    math::ft::fftT<complexT, complexT, rank, 0> m_fft_bwd; ///< FFT object for the backward transfsorm.
 
  public:
    /// C'tor.
    /**
     */
    psdFilter();
 
    /// Destructor
    ~psdFilter();
 
  protected:
    /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
    /// PSD.
    /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
     * deallocation, etc.
     *
     * See the discussion of PSD normalization above.
     *
     * This private version handles the actual setting of m_psdSqrt, which is rank independent.
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    int psdSqrt( realArrayT *npsdSqrt /**< [in] a pointer to an array containing the square root of the PSD. */ );
 
    /// Set the sqaure-root of the PSD.
    /** This allocates _npsdSqrt and fills it with the values in the array.
     *
     * See the discussion of PSD normalization above.
     *
     * This private version handles the actual setting of m_psdSqrt, which is rank independent.
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    int psdSqrt( const realArrayT &npsdSqrt /**< [in] an array containing the square root of the PSD. */ );
 
    /// Set the size of the filter.
    /** Handles allocation of the m_ftWork array and fftw planning.
     *
     * Requires m_psdSqrt to be set first.  This is called by the psdSqrt() and psd() methods.
     *
     * This version compiles when rank==1
     *
     */
    template <size_t crank = rank>
    int setSize( typename std::enable_if<crank == 1>::type * = 0 );
 
    /// Set the size of the filter.
    /** Handles allocation of the m_ftWork array and fftw planning.
     *
     * Requires m_psdSqrt to be set first.  This is called by the psdSqrt() and psd() methods.
     *
     * This version compiles when rank==2
     *
     */
    template <size_t crank = rank>
    int setSize( typename std::enable_if<crank == 2>::type * = 0 );
 
    /// Set the size of the filter.
    /** Handles allocation of the m_ftWork array and fftw planning.
     *
     * Requires m_psdSqrt to be set first.  This is called by the psdSqrt() and psd() methods.
     *
     * This version compiles when rank==3
     *
     */
    template <size_t crank = rank>
    int setSize( typename std::enable_if<crank == 3>::type * = 0 );
 
  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();
 
    /// Get the number of planes in the filter
    /**
     * \returns the current value of m_planes.
     *
     */
    int planes();
 
    /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
    /// PSD.
    /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
     * deallocation, etc.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==1
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
                 realT df,             ///< [in] the frequency spacing
                 typename std::enable_if<crank == 1>::type * = 0 );
 
    /// Set the square-root of the PSD to be a pointer to an array containing the square root of the properly normalized
    /// PSD.
    /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
     * deallocation, etc.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==2
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
                 realT dk1,            ///< [in] the frequency spacing along dimension 1
                 realT dk2,            ///< [in] the frequency spacing along dimension 2
                 typename std::enable_if<crank == 2>::type * = 0 );
 
    /// Set the sqaure-root of the PSD to be a pointer to an array containing the square root of the properly normalized
    /// PSD.
    /** This does not allocate _npsdSqrt, it merely points to the specified array, which remains your responsibility for
     * deallocation, etc.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==3
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( realArrayT *npsdSqrt, ///< [in] a pointer to an array containing the square root of the PSD.
                 realT dk1,            ///< [in] the frequency spacing along dimension 1
                 realT dk2,            ///< [in] the frequency spacing along dimension 2
                 realT df,             ///< [in] the frequency spacing along dimension 3
                 typename std::enable_if<crank == 3>::type * = 0 );
 
    /// 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.
     *
     * This version compiles when rank==1
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
                 realT df,                   ///< [in] the frequency spacing
                 typename std::enable_if<crank == 1>::type * = 0 );
 
    /// 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.
     *
     * This version compiles when rank==2
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
                 realT dk1,                  ///< [in] the frequency spacing along dimension 1
                 realT dk2,                  ///< [in] the frequency spacing along dimension 2
                 typename std::enable_if<crank == 2>::type * = 0 );
 
    /// 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.
     *
     * This version compiles when rank==3
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psdSqrt( const realArrayT &npsdSqrt, ///< [in] an array containing the square root of the PSD.
                 realT dk1,                  ///< [in] the frequency spacing along dimension 1
                 realT dk2,                  ///< [in] the frequency spacing along dimension 2
                 realT df,                   ///< [in] the frequency spacing along dimension 3
                 typename std::enable_if<crank == 3>::type * = 0 );
 
    /// Set the sqaure-root of the PSD from the PSD.
    /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==1
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
             const realT df,         ///< [in] the frequency spacing
             typename std::enable_if<crank == 1>::type * = 0 );
 
    /// Set the sqaure-root of the PSD from the PSD.
    /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==2
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
             const realT dk1,        ///< [in] the frequency spacing along dimension 1
             const realT dk2,        ///< [in] the frequency spacing along dimension 2
             typename std::enable_if<crank == 2>::type * = 0 );
 
    /// Set the sqaure-root of the PSD from the PSD.
    /** This allocates _npsdSqrt and fills it with the square root of the values in the array.
     *
     * See the discussion of PSD normalization above.
     *
     * This version compiles when rank==3
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int psd( const realArrayT &npsd, ///< [in] an array containing the PSD
             const realT dk1,        ///< [in] the frequency spacing along dimension 1
             const realT dk2,        ///< [in] the frequency spacing along dimension 2
             const realT df,         ///< [in] the frequency spacing along dimension 3
             typename std::enable_if<crank == 3>::type * = 0 );
 
    /// De-allocate all working memory and reset to initial state.
    /**
     *
     */
    void clear();
 
    /// Apply the filter.
    /**
     * This version compiles when rank==1
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
                realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
                                               ///< filter, allowing 2-for-1 calculation.
                typename std::enable_if<crank == 1>::type * = 0 ) const;
 
    /// Apply the filter.
    /**
     * This version compiles when rank==2
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
                realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
                                               ///< filter, allowing 2-for-1 calculation.
                typename std::enable_if<crank == 2>::type * = 0 ) const;
 
    /// Apply the filter.
    /**
     * This version compiles when rank==2
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int filter( realArrayMapT noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
                realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
                                               ///< filter, allowing 2-for-1 calculation.
                typename std::enable_if<crank == 2>::type * = 0 ) const;
 
    /// Apply the filter.
    /**
     * This version compiles when rank==3
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    template <size_t crank = rank>
    int filter( realArrayT &noise, ///< [in.out] the noise field of size rows() X cols(), which is filtered in-place.
                realArrayT *noiseIm = nullptr, ///< [out] [optional] an array to fill with the imaginary output of the
                                               ///< filter, allowing 2-for-1 calculation.
                typename std::enable_if<crank == 3>::type * = 0 ) const;
 
    /// 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. */ ) const;
 
    /// Apply the filter.
    /**
     * \overload
     *
     * \returns 0 on success
     * \returns -1 on error
     *
     */
    int operator()(
        realArrayMapT noise /**< [in.out] the noise field of size rows() X cols(), which is filtered in-place. */ )
        const;
 
    /// 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.
    ) const;
};
 
template <typename realT, size_t rank>
psdFilter<realT, rank>::psdFilter()
{
}
 
template <typename realT, size_t rank>
psdFilter<realT, rank>::~psdFilter()
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt )
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
 
    m_psdSqrt = npsdSqrt;
    m_owner = false;
 
    setSize();
 
    return 0;
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt )
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
 
    m_psdSqrt = new realArrayT;
 
    ( *m_psdSqrt ) = npsdSqrt;
    m_owner = true;
 
    setSize();
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 1>::type * )
{
    if( m_psdSqrt == 0 )
    {
        internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
        return -1;
    }
 
    if( m_rows == m_psdSqrt->size() )
    {
        return 0;
    }
 
    m_rows = m_psdSqrt->size();
    m_cols = 1;
    m_planes = 1;
 
    m_ftWork.resize( m_rows );
 
    m_fft_fwd.plan( m_rows, math::ft::dir::forward, true );
 
    m_fft_bwd.plan( m_rows, math::ft::dir::backward, true );
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 2>::type * )
{
    if( m_psdSqrt == 0 )
    {
        internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
        return -1;
    }
 
    if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() )
    {
        return 0;
    }
 
    m_rows = m_psdSqrt->rows();
    m_cols = m_psdSqrt->cols();
    m_planes = 1;
 
    m_ftWork.resize( m_rows, m_cols );
 
    m_fft_fwd.plan( m_rows, m_cols, math::ft::dir::forward, true );
 
    m_fft_bwd.plan( m_rows, m_cols, math::ft::dir::backward, true );
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::setSize( typename std::enable_if<crank == 3>::type * )
{
    if( m_psdSqrt == 0 )
    {
        internal::mxlib_error_report(error_t::paramnotset, "m_psdSqrt has not been set yet, is still NULL." );
        return -1;
    }
 
    if( m_rows == m_psdSqrt->rows() && m_cols == m_psdSqrt->cols() && m_planes == m_psdSqrt->planes() )
    {
        return 0;
    }
 
    m_rows = m_psdSqrt->rows();
    m_cols = m_psdSqrt->cols();
    m_planes = m_psdSqrt->planes();
 
    m_ftWork.resize( m_rows, m_cols, m_planes );
 
    m_fft_fwd.plan( m_planes, m_rows, m_cols, math::ft::dir::forward, true );
 
    m_fft_bwd.plan( m_planes, m_rows, m_cols, math::ft::dir::backward, true );
 
    return 0;
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::rows()
{
    return m_rows;
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::cols()
{
    return m_cols;
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::planes()
{
    return m_planes;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
{
    m_dFreq1 = df;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt( realArrayT *npsdSqrt,
                                     realT dk1,
                                     realT dk2,
                                     typename std::enable_if<crank == 2>::type * )
{
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt(
    realArrayT *npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
{
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
    m_dFreq3 = df;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt, realT df, typename std::enable_if<crank == 1>::type * )
{
    m_dFreq1 = df;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt( const realArrayT &npsdSqrt,
                                     realT dk1,
                                     realT dk2,
                                     typename std::enable_if<crank == 2>::type * )
{
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psdSqrt(
    const realArrayT &npsdSqrt, realT dk1, realT dk2, realT df, typename std::enable_if<crank == 3>::type * )
{
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
    m_dFreq3 = df;
    return psdSqrt( npsdSqrt );
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psd( const realArrayT &npsd, const realT df1, typename std::enable_if<crank == 1>::type * )
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
 
    m_psdSqrt = new realArrayT;
 
    // Vector
    m_psdSqrt->resize( npsd.size() );
    for( size_t n = 0; n < npsd.size(); ++n )
        ( *m_psdSqrt )[n] = sqrt( npsd[n] );
 
    m_owner = true;
 
    m_dFreq1 = df1;
 
    setSize();
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psd( const realArrayT &npsd,
                                 const realT dk1,
                                 const realT dk2,
                                 typename std::enable_if<crank == 2>::type * )
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
 
    m_psdSqrt = new realArrayT;
 
    ( *m_psdSqrt ) = npsd.sqrt();
    m_owner = true;
 
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
 
    setSize();
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::psd( const realArrayT &npsd,
                                 const realT dk1,
                                 const realT dk2,
                                 const realT df,
                                 typename std::enable_if<crank == 3>::type * )
{
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
    }
 
    m_psdSqrt = new realArrayT;
 
    // Cube
    m_psdSqrt->resize( npsd.rows(), npsd.cols(), npsd.planes() );
    for( int pp = 0; pp < npsd.planes(); ++pp )
    {
        for( int cc = 0; cc < npsd.cols(); ++cc )
        {
            for( int rr = 0; rr < npsd.rows(); ++rr )
            {
                m_psdSqrt->image( pp )( rr, cc ) = sqrt( npsd.image( pp )( rr, cc ) );
            }
        }
    }
 
    m_owner = true;
 
    m_dFreq1 = dk1;
    m_dFreq2 = dk2;
    m_dFreq3 = df;
 
    setSize();
 
    return 0;
}
 
template <typename realT, size_t rank>
void psdFilter<realT, rank>::clear()
{
    // m_ftWork.resize(0,0);
    psdFilterTypes::arrayT<realT, rank>::clear( m_ftWork );
 
    m_rows = 0;
    m_cols = 0;
    m_planes = 0;
 
    if( m_psdSqrt && m_owner )
    {
        delete m_psdSqrt;
        m_psdSqrt = 0;
    }
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::filter( realArrayT &noise,
                                    realArrayT *noiseIm,
                                    typename std::enable_if<crank == 1>::type * ) const
{
    for( int nn = 0; nn < noise.size(); ++nn )
        m_ftWork[nn] = complexT( noise[nn], 0 );
 
    // Transform complex noise to Fourier domain.
    m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
 
    // Apply the filter.
    for( int nn = 0; nn < m_ftWork.size(); ++nn )
        m_ftWork[nn] *= ( *m_psdSqrt )[nn];
 
    m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
 
    // Now take the real part, and normalize.
    realT norm = sqrt( noise.size() / m_dFreq1 );
    for( int nn = 0; nn < m_ftWork.size(); ++nn )
        noise[nn] = m_ftWork[nn].real() / norm;
 
    if( noiseIm != nullptr )
    {
        for( int nn = 0; nn < m_ftWork.size(); ++nn )
            ( *noiseIm )[nn] = m_ftWork[nn].imag() / norm;
    }
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::filter( realArrayT &noise,
                                    realArrayT *noiseIm,
                                    typename std::enable_if<crank == 2>::type * ) const
{
    // Make noise a complex number
    for( int ii = 0; ii < noise.rows(); ++ii )
    {
        for( int jj = 0; jj < noise.cols(); ++jj )
        {
            m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
        }
    }
 
    // Transform complex noise to Fourier domain.
    m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
 
    // Apply the filter.
    m_ftWork *= *m_psdSqrt;
 
    m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
 
    realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
 
    // Now take the real part, and normalize.
    noise = m_ftWork.real() / norm;
 
    if( noiseIm != nullptr )
    {
        *noiseIm = m_ftWork.imag() / norm;
    }
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::filter( realArrayMapT noise,
                                    realArrayT *noiseIm,
                                    typename std::enable_if<crank == 2>::type * ) const
{
    // Make noise a complex number
    for( int ii = 0; ii < noise.rows(); ++ii )
    {
        for( int jj = 0; jj < noise.cols(); ++jj )
        {
            m_ftWork( ii, jj ) = complexT( noise( ii, jj ), 0 );
        }
    }
 
    // Transform complex noise to Fourier domain.
    m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
 
    // Apply the filter.
    m_ftWork *= *m_psdSqrt;
 
    m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
 
    realT norm = sqrt( noise.rows() * noise.cols() / ( m_dFreq1 * m_dFreq2 ) );
 
    // Now take the real part, and normalize.
    noise = m_ftWork.real() / norm;
 
    if( noiseIm != nullptr )
    {
        *noiseIm = m_ftWork.imag() / norm;
    }
 
    return 0;
}
 
template <typename realT, size_t rank>
template <size_t crank>
int psdFilter<realT, rank>::filter( realArrayT &noise,
                                    realArrayT *noiseIm,
                                    typename std::enable_if<crank == 3>::type * ) const
{
    // Make noise a complex number
    for( int pp = 0; pp < noise.planes(); ++pp )
    {
        for( int ii = 0; ii < noise.rows(); ++ii )
        {
            for( int jj = 0; jj < noise.cols(); ++jj )
            {
                m_ftWork.image( pp )( ii, jj ) = complexT( noise.image( pp )( ii, jj ), 0 );
            }
        }
    }
 
    // Transform complex noise to Fourier domain.
    m_fft_fwd( m_ftWork.data(), m_ftWork.data() );
 
    // Apply the filter.
    for( int pp = 0; pp < noise.planes(); ++pp )
        m_ftWork.image( pp ) *= m_psdSqrt->image( pp );
 
    m_fft_bwd( m_ftWork.data(), m_ftWork.data() );
 
    // Now take the real part, and normalize.
 
    realT norm = sqrt( m_rows * m_cols * m_planes / ( m_dFreq1 * m_dFreq2 * m_dFreq3 ) );
    for( int pp = 0; pp < noise.planes(); ++pp )
        noise.image( pp ) = m_ftWork.image( pp ).real() / norm;
 
    if( noiseIm != nullptr )
    {
        for( int pp = 0; pp < noise.planes(); ++pp )
            noiseIm->image( pp ) = m_ftWork.image( pp ).imag() / norm;
    }
 
    return 0;
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::operator()( realArrayT &noise ) const
{
    return filter( noise );
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::operator()( realArrayMapT noise ) const
{
    return filter( noise );
}
 
template <typename realT, size_t rank>
int psdFilter<realT, rank>::operator()( realArrayT &noise, realArrayT &noiseIm ) const
{
    return filter( noise, &noiseIm );
}
 
} // namespace sigproc
} // namespace mx
 
#endif // psdFilter_hpp