0
$\begingroup$

I am trying to implement a radix-4 DIT FFT. It looks like the forward transform is working correctly, but the backward transform output is not in the correct order. I'm pretty sure that the digitrevorder function is not the problem, as I also use it in radix-3 and radix-5 implementations. Am I missing something obvious here? This is my third attempt, using 2 books and a Python implementation as references. All have the same problem.

Code (C++)

#include <algorithm>
#include <cmath>
#include <complex>
#include <concepts>
#include <cstdio>
#include <numbers>
#include <span>

template<std::integral T>
[[nodiscard]] constexpr auto ipow(T base, T exponent) -> T
{
    T result = 1;
    for (T i = 0; i < exponent; i++) {
        result *= base;
    }
    return result;
}

/// Reorder input using base-radix digit reversal permutation.
template<std::size_t Radix>
auto digitrevorder(std::span<std::complex<double>> x) noexcept -> void
{
    auto const len = x.size();

    auto j = 0UL;
    for (auto i = 0UL; i < len - 1UL; i++) {
        if (i < j) {
            std::swap(x[i], x[j]);
        }
        auto k = (Radix - 1UL) * len / Radix;
        while (k <= j) {
            j -= k;
            k /= Radix;
        }
        j += k / (Radix - 1);
    }
}

// Chapter 2.4.3: The Cooley-Tukey Radix-4 Algorithm
// Computational frameworks for the fast fourier transform
// ISBN: 978-0-89871-285-8
auto c2c_dit4(std::span<std::complex<double>> x, std::size_t order, double sign) -> void
{
    static constexpr auto Radix = 4UL;
    static constexpr auto TwoPi = 3.14159265359 * 2.0;
    static constexpr auto I     = std::complex{0.0, 1.0};

    auto const t = order;
    auto const n = x.size();

    digitrevorder<Radix>(x);

    for (auto q{1UL}; q <= t; ++q) {
        auto const L  = ipow(Radix, q);
        auto const LR = L / Radix;
        auto const r  = n / L;

        for (auto j{0UL}; j < LR; ++j) {
            auto const angle = sign * TwoPi * static_cast<double>(j) / static_cast<double>(L);
            auto const w1    = std::polar(1.0, angle * 1.0);
            auto const w2    = std::polar(1.0, angle * 2.0);
            auto const w3    = std::polar(1.0, angle * 3.0);

            for (auto k{0UL}; k < r; ++k) {
                auto const a = x[k * L + LR * 0 + j];
                auto const b = x[k * L + LR * 1 + j] * w1;
                auto const c = x[k * L + LR * 2 + j] * w2;
                auto const d = x[k * L + LR * 3 + j] * w3;

                auto const t0 = a + c;
                auto const t1 = a - c;
                auto const t2 = b + d;
                auto const t3 = b - d;

                x[k * L + LR * 0 + j] = t0 + t2;
                x[k * L + LR * 1 + j] = t1 - t3 * I;
                x[k * L + LR * 2 + j] = t0 - t2;
                x[k * L + LR * 3 + j] = t1 + t3 * I;
            }
        }
    }
}

auto main() -> int
{
    auto print = [](auto const* msg, auto x, auto scale) {
        std::printf("%s\n", msg);
        for (auto z : x) {
            std::printf("%+.4f,%+.4f\n", z.real() * scale, z.imag() * scale);
        }
        std::puts("");
    };

    static constexpr auto Order = 1;
    static constexpr auto Size  = ipow(4, Order);

    auto x = std::array<std::complex<double>, Size>{};
    x[0]   = 1.0;
    x[1]   = 2.0;
    x[2]   = 3.0;
    print("input", x, 1.0);

    c2c_dit4(x, Order, -1.0);
    print("fwd", x, 1.0);

    c2c_dit4(x, Order, +1.0);
    print("bwd", x, 1.0 / double(Size));

    return 0;
}

Output

input
+1.0000,+0.0000
+2.0000,+0.0000
+3.0000,+0.0000
+0.0000,+0.0000

fwd
+6.0000,+0.0000
-2.0000,-2.0000
+2.0000,+0.0000
-2.0000,+2.0000

bwd
+1.0000,+0.0000
+0.0000,+0.0000
+3.0000,+0.0000
+2.0000,+0.0000
$\endgroup$

1 Answer 1

1
$\begingroup$

Notice that your result for ifft(fft([1,2,3,0])) = [1,0,3,2] coincides with fft(fft([1,2,3,0]))

Without going deep in your implementation details I am notice that you have

    static constexpr auto I     = std::complex{0.0, 1.0};

that is used as a precomputed value for std::polar(1, M_PI / 2) in the forward transform, so you probably should define it as

    static constexpr auto I     = std::complex{0.0, sign};

to reflect that the fact that the in the inverse fft every std::polar(1, theta) is replaced by std::polar(1, -theta).

$\endgroup$
2
  • 1
    $\begingroup$ Awesome. That solved it. Working now $\endgroup$ Jan 12 at 10:35
  • $\begingroup$ Thank you for the feedback, now take a look on dsp.stackexchange.com/tour. Welcome to the comunity :D. $\endgroup$
    – Bob
    Jan 12 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.