rfft.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <complex>
 
#include "NumCpp/Core/Constants.hpp"
#include "NumCpp/Core/Internal/StaticAsserts.hpp"
#include "NumCpp/Core/Internal/StlAlgorithms.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Functions/complex.hpp"
#include "NumCpp/NdArray.hpp"
 
namespace nc::fft
{
    namespace detail
    {
        //===========================================================================
        // Method Description:
        inline NdArray<std::complex<double>> rfft_internal(const NdArray<std::complex<double>>& x, uint32 n)
        {
            if (n == 0)
            {
                return {};
            }
 
            const auto realN  = n / 2 + 1;
            auto       result = NdArray<std::complex<double>>(1, realN);
 
            stl_algorithms::for_each(result.begin(),
                                     result.end(),
                                     [&](auto& resultElement)
                                     {
                                         const auto k = static_cast<double>(&resultElement - result.data());
                                         const auto minusTwoPiKOverN = -constants::twoPi * k / static_cast<double>(n);
                                         resultElement               = std::complex<double>{ 0., 0. };
                                         for (auto m = 0u; m < std::min(n, x.size()); ++m)
                                         {
                                             const auto angle = minusTwoPiKOverN * static_cast<double>(m);
                                             resultElement += (x[m] * std::polar(1., angle));
                                         }
                                     });
 
            return result;
        }
    } // namespace detail
 
    //===========================================================================
    // Method Description:
    template<typename dtype>
    NdArray<std::complex<double>> rfft(const NdArray<dtype>& inArray, uint32 inN, Axis inAxis = Axis::NONE)
    {
        STATIC_ASSERT_ARITHMETIC(dtype);
 
        switch (inAxis)
        {
            case Axis::NONE:
            {
                const auto data = nc::complex<dtype, double>(inArray);
                return detail::rfft_internal(data, inN);
            }
            case Axis::COL:
            {
                auto        data           = nc::complex<dtype, double>(inArray);
                const auto& shape          = inArray.shape();
                const auto  realN          = inN / 2 + 1;
                auto        result         = NdArray<std::complex<double>>(shape.rows, realN);
                const auto  dataColSlice   = data.cSlice();
                const auto  resultColSlice = result.cSlice();
 
                for (uint32 row = 0; row < data.numRows(); ++row)
                {
                    const auto rowData   = data(row, dataColSlice);
                    const auto rowResult = detail::rfft_internal(rowData, inN);
                    result.put(row, resultColSlice, rowResult);
                }
 
                return result;
            }
            case Axis::ROW:
            {
                return rfft(inArray.transpose(), inN, Axis::COL).transpose();
            }
            default:
            {
                THROW_INVALID_ARGUMENT_ERROR("Unimplemented axis type.");
                return {};
            }
        }
    }
 
    //===========================================================================
    // Method Description:
    template<typename dtype>
    NdArray<std::complex<double>> rfft(const NdArray<dtype>& inArray, Axis inAxis = Axis::NONE)
    {
        STATIC_ASSERT_ARITHMETIC(dtype);
 
        switch (inAxis)
        {
            case Axis::NONE:
            {
                return rfft(inArray, inArray.size(), inAxis);
            }
            case Axis::COL:
            {
                return rfft(inArray, inArray.numCols(), inAxis);
            }
            case Axis::ROW:
            {
                return rfft(inArray, inArray.numRows(), inAxis);
            }
            default:
            {
                THROW_INVALID_ARGUMENT_ERROR("Unimplemented axis type.");
                return {};
            }
        }
    }
} // namespace nc::fft