kissfft

How fast is Kiss FFT?


Topic question: how fast is Kiss FFT in big O notation? I can't find this stated anywhere, it's probably somewhere obvious and I'm just blind.


Solution

  • I'm pretty sure the point of a Fast Fourier Transform is that it's θ(n log n).
    If it's anything slower then it wouldn't be called a "Fast" Fourier Transform.
    (And it's probably not faster.)


    Speaking of "simple" FFTs, here's one (powers of 2 only):

    // To get the discrete Fourier _series_, divide by the period
    template<class RanInIt, class RanOutIt, class RanScratchIt>
    void fft(RanInIt begin, RanInIt end,
        RanOutIt outputs, RanScratchIt scratch /* at least twice as big as output */,
        bool const inverse = false, size_t const stride = 1)
    {
        typedef std::iterator_traits<RanInIt>::value_type complex;
        static long double const two_pi = 6.2831853071795864769252867665590057683943L;
        size_t const n = 1 + static_cast<size_t>(std::distance(begin, end) - 1) / stride;
        if ((n & (n - 1)) != 0) { throw std::logic_error("FFTs must be in powers of 2!"); }
        complex const divisor(inverse ? n : 1);
        if (n > 2)
        {
            fft(begin + 0 * stride, end, scratch + 0 * stride, outputs + 0 * stride, false, 2 * stride);
            fft(begin + 1 * stride, end, scratch + 1 * stride, outputs + 1 * stride, false, 2 * stride);
            size_t const n_2 = n / 2;
            for (size_t i = 0; i < n_2; i++)
            {
                using std::polar;
                complex const
                    t(polar(1.0L, -two_pi * i / n)),
                    s  (*(scratch + (2 * i + 0) * stride)),
                    psi(*(scratch + (2 * i + 1) * stride) * t);
                *(outputs + ((i +   0) * stride)) = (s + psi) / divisor;
                *(outputs + ((i + n_2) * stride)) = (s - psi) / divisor;
            }
        }
        else if (n > 0)
        {
            *(outputs + (0 * stride)) = (*begin + *(begin + 1 * stride)) / divisor;
            *(outputs + (1 * stride)) = (*begin - *(begin + 1 * stride)) / divisor;
        }
        if (inverse)
        {
            std::swap_ranges(outputs + 1, outputs + n / 2, std::reverse_iterator<RanOutIt>(outputs + n));
        }
    }
    
    template<class RanInIt>
    std::vector<typename std::iterator_traits<RanInIt>::value_type>
        fft(RanInIt begin, RanInIt end, bool const inverse = false)
    {
        size_t const n = static_cast<size_t>(std::distance(begin, end));
        size_t log2ceil_n = 0;
        for (size_t i = 1; i < n; i *= 2) { log2ceil_n++; }
        std::vector<typename std::iterator_traits<RanInIt>::value_type>
            outputs(static_cast<size_t>(1) << log2ceil_n),
            scratch(2 * outputs.size());
        fft(begin, end, outputs.begin(), scratch.begin(), inverse);
        outputs.resize(n);
        return outputs;
    }
    
    template<class Range>
    std::vector<typename Range::value_type> fft(Range range, bool const inverse = false)
    { return fft(range.begin(), range.end(), inverse); }
    
    template<class T, size_t N>
    std::vector<T> fft(T const (&range)[N], bool const inverse = false)
    { return fft(&range[0], &range[N], inverse); }