/* dftmerge.c - DFT merging for radix 2 decimation-in-time FFT */

#include <cmplx.h>

void dftmerge(N, XF)
complex  *XF;
int N;
{
       double pi = 4. * atan(1.0);
       int k, i, p, q, M;
       complex  A, B, V, W;

       M = 2;
       while (M <= N) {                          /* two \((M/2)\)-DFTs into one \(M\)-DFT  */
            W = cexp(cmplx(0.0, -2 * pi / M));   /* order-\(M\) twiddle factor */
            V = cmplx(1., 0.);                   /* successive powers of \(W\) */
            for (k = 0; k < M/2; k++) {          /* index for an \((M/2)\)-DFT */
                 for (i = 0; i < N; i += M) {    /* \(i\)th butterfly; increment by \(M\) */
                      p = k + i;                 /* absolute indices for */
                      q = p + M / 2;             /* \(i\)th butterfly */
                      A = XF[p];
                      B = cmul(XF[q], V);        /* \(V = W\sp{k}\) */
                      XF[p] = cadd(A, B);        /* butterfly operations */
                      XF[q] = csub(A, B);
                      }
                 V = cmul(V, W);                 /* \(V = VW = W\sp{k+1}\) */
                 }
            M = 2 * M;                           /* next stage */
            }
}