Conjugate split-radix reference FFT.

gistfile1.cpp

typedef std::complex<float> complexf;

static size_t const kMaxN = 2048; // largest supported FFT size.

// FFT twiddles. Rather than access strided for smaller sub-FFTs, just
// have a "mip chain" of pow2 sizes up to kMaxN.
//
// The N-element subarray starts at position N.
static complexf s_twiddles[(kMaxN / 4) * 2];

static void init_fft()
{
    // init twiddle tables
    for (size_t N = 1; N <= kMaxN / 4; N *= 2)
    {
        double const kPi = 3.1415926535897932384626433832795;
        double step = -2.0 * kPi / (double)(N * 4);
        complexf *twiddle = s_twiddles + N;

        for (size_t k = 0; k < N; k++)
        {
            double phase = step * (double)k;
            twiddle[k] = complexf((float)cos(phase), (float)sin(phase));
        }
    }
}

// Reference conjugate split-radix FFT (DIT) - recursion step
static void conj_fft_rec(complexf out[], complexf const in[], size_t offs, size_t mask, size_t stride, size_t N)
{
    if (N == 1)
        out[0] = in[offs & mask];
    else if (N == 2)
    {
        complexf a = in[offs & mask];
        complexf b = in[(offs + stride) & mask];
        out[0] = a + b;
        out[1] = a - b;
    }
    else
    {
        conj_fft_rec(out, in, offs, mask, stride * 2, N / 2);
        conj_fft_rec(out + N/2, in, offs + stride, mask, stride * 4, N / 4);
        conj_fft_rec(out + 3*N/4, in, offs - stride, mask, stride * 4, N / 4);

        complexf const *twiddle = s_twiddles + N/4;
        for (size_t k = 0; k < N/4; k++)
        {
            complexf Uk  = out[k];
            complexf Zk  = out[k + N/2];
            complexf Uk2 = out[k + N/4];
            complexf Zdk = out[k + 3*N/4];
            complexf const &w = twiddle[k];

            // Twiddle
            Zk *= w;
            Zdk *= std::conj(w);

            // Z butterflies
            complexf Zsum = Zk + Zdk;
            complexf Zdif = Zk - Zdk;

            out[k]       = Uk  + Zsum;
            out[k+N/2]   = Uk  - Zsum;
            out[k+N/4]   = Uk2 - complexf(0.0f, 1.0f)*Zdif;
            out[k+3*N/4] = Uk2 + complexf(0.0f, 1.0f)*Zdif;
        }
    }
}

// Call this
static void conj_fft(complexf out[], complexf const in[], size_t N)
{
    conj_fft_rec(out, in, 0, N - 1, 1, N);
}