simple FFT (template, std::complex, C++11, header-only)

#ifndef FFT_H_INCLUDED
#define FFT_H_INCLUDED


//   fft.h - template functions for FFT transform
//   The code is derived from LIBROW FFT code (www.librow.com)

#include <complex>


namespace fft {


    namespace _ {


        template<typename T> struct ComplexType { typedef std::complex< typename std::remove_reference<typename std::remove_const<T>::type>::type> type;};
        template<typename T> struct ComplexType<std::complex<T>> {typedef typename std::complex<T> type;};

        inline void updateTarget(size_t& target, const size_t size)
        {
            size_t Mask = size;
            while (target & (Mask >>= 1))
                target &= ~Mask;
            target |= Mask;
        }



        //! @return destination's end iterator
        template<typename IterSrc, typename IterDst>
        IterDst rearrange(IterSrc srcFirst, IterSrc srcLast, IterDst dstFirst)
        {
            const size_t size = srcLast - srcFirst;
            size_t target = 0;
            for (; srcFirst != srcLast; ++srcFirst)
            {
                *(dstFirst + target) = *srcFirst;
                updateTarget(target, size);
            }
            return dstFirst + size;
        }



        //! Inplace version of rearrange function
        template<typename IterDst>
        void rearrange(IterDst dstFirst, IterDst dstLast)
        {
            size_t target = 0;
            const size_t size = dstLast - dstFirst;
            for (auto i = dstFirst; i != dstLast; ++i)
            {
                if (target > size_t(i - dstFirst))
                    std::swap(*i, *(dstFirst + target));
                updateTarget(target, size);
            }
        }


        //! Actual direct or inverse FFT on rearranged data
        template<typename IterDst>
        void perform(IterDst dstFirst, IterDst dstLast, const bool inverse = false)
        {
            typedef typename std::remove_reference<decltype(*dstFirst)>::type Complex;
            typedef typename Complex::value_type Number;
            const Number pi = inverse ? Number(3.14159265358979323846) : Number(-3.14159265358979323846);
            const size_t size = dstLast - dstFirst;
            for (size_t step = 1; step < size; step <<= 1)
            {
                const size_t jump = step << 1;
                const Number delta = pi / Number(step);
                const Number sine = std::sin(delta * .5);
                const Complex multiplier(-2. * sine * sine, std::sin(delta));
                Complex factor(1.);
                for (size_t group = 0; group < step; ++group)
                {
                    for (size_t pair = group; pair < size; pair += jump)
                    {
                        const size_t match = pair + step;
                        const Complex product(factor * dstFirst[match]);
                        dstFirst[match] = dstFirst[pair] - product;
                        dstFirst[pair] += product;
                    }
                    factor = multiplier * factor + factor;
                }
            }
        }

        //! Scaling of inverse FFT result
        template<typename IterDst>
        void scale(IterDst dstFirst, IterDst dstLast)
        {
            typedef typename std::remove_reference<decltype(*dstFirst)>::type Complex;
            typedef typename Complex::value_type Number;
            const Number factor = 1. / Number(dstLast - dstFirst);
            for (; dstFirst != dstLast; ++dstFirst)
                *dstFirst *= factor;
        }



    }

    //! Plain FFT
    //! @tparam IterSrc source iterator (read-only, random access), dereferences to numeric type Number1
    //! @tparam IterDst destination iterator (writeable, random access), dereferences to complex<Number2>, Number2 must be constructible from Number1
    //! @pre dstFirst has capacity at least srcLast - srcFirst, and srcLast - srcFirst must be power of 2
    //! @return destination's end iterator
    template<typename IterSrc, typename IterDst>
    IterDst transform (IterSrc srcFirst, IterSrc srcLast, IterDst dstFirst)
    {
        const auto dstLast = ::fft::_::rearrange(srcFirst, srcLast, dstFirst);
        ::fft::_::perform(dstFirst, dstLast);
        return dstLast;
    }

    //! In-place FFT
    //! @tparam IterDst destination iterator (writeable, random access), dereferences to complex<Number>
    //! @pre srcLast - srcFirst must be power of 2
    template<typename IterDst>
    void transform (IterDst dstFirst, IterDst dstLast)
    {
        ::fft::_::rearrange(dstFirst, dstLast);
        ::fft::_::perform(dstFirst, dstLast);
    }

    //! Inverse FFT
    //! @tparam IterDst destination iterator (writeable, random access), dereferences to complex<Number>
    //! @pre dstFirst has capacity at least srcLast - srcFirst, and srcLast - srcFirst must be power of 2
    template<typename IterDst>
    void inverse (IterDst srcFirst, IterDst srcLast, IterDst dstFirst, const bool doScale = true)
    {
        const auto dstLast = ::fft::_::rearrange(srcFirst, srcLast, dstFirst);
        ::fft::_::perform(dstFirst, dstLast, true);
        if(doScale)
            ::fft::_::scale(dstFirst, dstLast);
    }

    //! In-place inverse FFT
    //! @tparam IterDst destination iterator (writeable, random access), dereferences to complex<Number>
    //! @pre srcLast - srcFirst must be power of 2
    template<typename IterDst>
    void inverse (IterDst dstFirst, IterDst dstLast, const bool doScale = true)
    {
        ::fft::_::rearrange(dstFirst, dstLast);
        ::fft::_::perform(dstFirst, dstLast, true);
        if(doScale)
            ::fft::_::scale(dstFirst, dstLast);
    }


}

#endif