Simple implementation of the Complex FFT in C99

fft.c

double complex *
fft(double complex iq[], int bins)
{
    int i, j, k;
    int q;
    double t;
    complex double alpha, uri, ur;
    complex double *foo;

    /* and they must be a power of 2 */
    t = log(bins) / log(2);
    if (modf(t, &t) > 0) {
        return NULL;
    }

    foo = calloc(bins, sizeof(complex double));

    q = (int) t;

    /* This first bit is a reflection sort,
     * Most people do a 'sort in place' of
     * the source data, but I'm trying to preserve
     * that original data for other use, so I
     * 'sort into place' from the source into
     * my temporary array foo.
     *
     * The end result is each entry is 2^n away
     * from its sibling.
     */
    for (i = 0; i < bins; i++) {
        /* compute reflected index */
        for (k = 0, j = 0; j < q; j++) {
            k = (k << 1) | ((i >> j) & 1);
        }
        *(foo + i)  = iq[k];
    }

    /*
     * now synthesize the frequency domain, 1 thru n
     * frequency domain 0 is easy, its just the value
     * in the bin because a 1 bin DFT is the spectrum
     * of that DFT.
     */

    for (i = 1; i <= q; i++) {
        int bfly_len = 1 << i;          /* Butterfly elements */
        int half_bfly = bfly_len / 2;       /* Half-the butterfly */

        /* unity root value (complex) */
        ur = 1.0;

        /* This is the unity root increment (complex)
         * If you multiply ur by this value 'n' times then
         * ur will return to [1 + 0i]. 'n' in this case is
         * bfly_len times. Technically the argument to the
         * trancendentals would be 2 * pi / butterfly-span
         * but since we calculate butterfly-span / 2 (half_bfly)
         * we use algebra to simplify math to pi / half-bfly.
         */
        uri = cos(M_PI / half_bfly) - sin(M_PI / half_bfly) * I;

        /*
         * Combine two 2^i DFTs into a single
         * 2^(i+1) DFT. So two 1 bin DFTs to
         * a 2 bin DFT, two 2 bin DFTs to a 4
         * bin DFT, etc. The number of times
         * we convert is a function of how many
         * 2^(i+1) bin DFTs are in the total
         * number of bins. So for 512 bins (example)
         * there are two hundred and fifty six  2-bin DFTs,
         * one hundred and twenty eight 4-bin DFTs, all
         * the way up to exactly one 512-bin DFT.
         */
        for (j = 0; j < half_bfly ; j++) {
            for (k = j; k < bins; k += bfly_len) {
                /*
                 * Apply the FFT butterfly function to
                 * P[n], P[n + half_bfly]
                 *
                 * Alpha is the 180 degrees out point multiplied by
                 * the current unit root.
                 */

                alpha = *(foo + k + half_bfly) * ur;
                /*
                 * Mathematically P[n + half_bfly] is 180 degrees
                 * 'further' than P[n]. So changing the sign on
                 * alpha (1 + 0i) => (-1 - 0i) is the equivalent
                 * value of alpha for that point in the transform.
                 *
                 * step 2, P[n]             = P[n] + alpha
                 *         P[n + half_bfly] = P[n] - alpha
                 *
                 * Sequence is important here, set P[n + h]
                 * before you change the value of P[n].
                 */
                *(foo + k + half_bfly) = *(foo + k) - alpha;
                *(foo + k) += alpha;
              
                /*
                 * and that is it, except for scaling perhaps, if you want
                 * to (or need to) keep the bins within the precision
                 * of their numeric representation.
                 */
            }
            /*
             * Now my multiplying UR by URI we save ourselves from
             * recomputing the sin and cos, Euler tells us we can just
             * keep multiplying and the values will go through a
             * sequence from 1 + 0, to 1 - i, to -1 + 0, to 1 + i and
             * then back to 1 + 0.
             */
            ur = ur * uri;
        }
    }
    return foo;
}