Discrete Fourier Transform With C, Implementation problems?

Question

I'm trying to understand some basics of DFT, some math equations, and try to implement it with C. Well, this is the function i used from a book (Algorithms for Image Processing And Computer Vision)

void slowft (float *x, COMPLEX *y, int n)
{
    COMPLEX tmp, z1, z2, z3, z4;
    int m, k;

/* Constant factor -2 pi */
    cmplx (0.0, (float)(atan (1.0)/n * -8.0), &tmp);
    
    printf (" constant factor -2 pi %f ", (float)(atan (1.0)/n * -8.0));
    for (m = 0; m<=n; m++)
    {
      NEXT();
      cmplx (x[0], 0.0, &(y[m]));
      for (k=1; k<=n-1; k++)
      {
/* Exp (tmp*k*m) */
        cmplx ((float)k, 0.0, &z2);
        cmult (tmp, z2, &z3);
        cmplx ((float)m, 0.0, &z2);
        cmult (z2, z3, &z4);
        cexp (z4, &z2);
/* *x[k] */
        cmplx (x[k], 0.0, &z3);
        cmult (z2, z3, &z4);
/* + y[m] */
        csum (y[m], z4, &z2);
        y[m].real = z2.real; y[m].imag = z2.imag;
      }
    }
}

So actually, I'm stuck on the Constant Factor part. I didn't understand:

1-) what it came from(especially arctan(1)) and

2-) what its purpose of it.

This is the equation of DFT:

And these are other functions that i used:

void cexp (COMPLEX z1, COMPLEX *res)
{
    COMPLEX x, y;

    x.real = exp((double)z1.real);
    x.imag = 0.0;
    y.real = (float)cos((double)z1.imag);
    y.imag = (float)sin((double)z1.imag);
    cmult (x, y, res);
}

void cmult (COMPLEX z1, COMPLEX z2, COMPLEX *res)
{
    res->real = z1.real*z2.real - z1.imag*z2.imag;
    res->imag = z1.real*z2.imag + z1.imag*z2.real;
}

void csum (COMPLEX z1, COMPLEX z2, COMPLEX *res)
{
    res->real = z1.real + z2.real;
    res->imag = z1.imag + z2.imag;
}

void cmplx (float rp, float ip, COMPLEX *z)
{
    z->real = rp;
    z->imag = ip;
}

float cnorm (COMPLEX z)
{
    return z.real*z.real + z.imag*z.imag;
}

Answer 1

Here is the code with comments explaining it.

void slowft (float *x, COMPLEX *y, int n)
{
    COMPLEX tmp, z1, z2, z3, z4;
    int m, k;

/* Constant factor -2 pi */
    cmplx (0.0, (float)(atan (1.0)/n * -8.0), &tmp);
        /*  atan(1) is π/4, so this sets tmp to -2πi/n.  Note that the i
            factor, the imaginary unit, comes from putting the expression in
            the second argument, which gives the imaginary portion of the
            complex number being assigned.  (It is written as "j" in the
            equation displayed in the question.  That is because engineers use
            "j" for i, having historically already used "i" for other purposes.)
        */
    
    printf (" constant factor -2 pi %f ", (float)(atan (1.0)/n * -8.0));
    for (m = 0; m<=n; m++)
    {
      NEXT();
        // Well, that is a frightening thing to see in code.  It is cryptic.

      cmplx (x[0], 0.0, &(y[m]));
        /*  This starts to calculate a sum that will be accumulated in y[m].
            The sum will be over k from 0 to n-1.  For the first term, k is 0,
            so -2πiwk/n will be 0.  The coefficient is e to the power of that,
            and e**0 is 1, so the first term is x[0] * 1, so we just put x[0]
            diretly in y[m] with no multiplication.
        */

      for (k=1; k<=n-1; k++)
        //  This adds the rest of the terms.
      {

/* Exp (tmp*k*m) */
        cmplx ((float)k, 0.0, &z2);
            //  This sets z2 to k.

        cmult (tmp, z2, &z3);
            /*  This multiplies the -2πi/n from above with k, so it puts
                -2πi/n from above, and This computes -2πik/n it in z3.
             */

        cmplx ((float)m, 0.0, &z2);
            //  This sets z2 to m.  m corresponds to the ω in the equation.

        cmult (z2, z3, &z4);
            //  This multiplies m by -2πik/n, putting -2πiwk/n in z4.

        cexp (z4, &z2);
            /*  This raises e to the power of -2πiwk/n, finishing the
                coefficient of the term in the sum.
            */

/* *x[k] */
        cmplx (x[k], 0.0, &z3);
            //  This sets z3 to x[k].

        cmult (z2, z3, &z4);
            //  This multiplies x[k] by the coefficient, e**(-2πiwk/n).

/* + y[m] */
        csum (y[m], z4, &z2);
            /*  This adds the term (z4) to the sum being accumulated (y[m])
                and puts the updated sum in z2.
            */

        y[m].real = z2.real; y[m].imag = z2.imag;
            /*  This moves the updated sum to y[m].  This is not necessary
                because csum is passed its operands as values, so they are
                copied when calling the function, and it is safe to update its
                output.  csum(y[m], z4, &y[m]) above would have worked.  But
                this works too.
            */
    }
}

Standard C has support for complex arithmetic, so it would be easier and clearer to include <complex.h> and write code this way:

void slowft(float *x, complex float *y, int n)
{
    static const float TwoPi = 0x3.243f6a8885a308d313198a2e03707344ap1f;

    float t0 = -TwoPi/n;

    for (int m = 0; m <=n; m++)
    {
        float t1 = t0*m;

        y[m] = x[0];
        for (int k = 1; k < n; k++)
            y[m] += x[k] * cexpf(t1 * k * I);
    }
}