C 的离散傅里叶变换，实现问题？

Question

我正在尝试了解 DFT 的一些基础知识、一些数学方程式，并尝试用 C 来实现它。好吧，这是我从一本书（图像处理和计算机视觉算法）中使用的函数

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;
      }
    }
}

所以实际上，我卡在了常数因子部分。没看懂：

1-) 它来自什么（特别是 arctan(1)）和

2-) 它的目的是什么。

这是DFT的方程式：

这些是我使用的其他功能：

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

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

上面的代码注释为您提供了线索：

/* Constant factor -2 pi */

... 虽然实际上正在计算的是 -2 pi / n（在产生复数的更广泛背景下，将其作为其虚部的系数） .观察正弦和余弦相等的角的切线值为 1。具有属性且在 [0, pi) 范围内的角度是 pi / 4，因此 atan(1.0) * -8.0 是（一个很好的近似值）-2 pi.

2-) what its purpose of it.

它（或者实际上是它的加法逆）出现在您提供的 DFT 方程中，因此它很自然地出现在旨在实现该公式的函数中。

Answer 2

这是带有注释解释的代码。

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.
            */
    }
}

标准 C 支持复数运算，因此包含 <complex.h> 并以这种方式编写代码会更容易和更清晰：

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);
    }
}

C 的离散傅里叶变换，实现问题？

Discrete Fourier Transform With C, Implementation problems?

c

dct