离散傅里叶变换（DFT）的C语言实现疑问：常数因子的来源与作用解析

我正在学习离散傅里叶变换（DFT）的基础原理与相关数学公式，并且尝试用C语言实现DFT。以下是我从《Algorithms for Image Processing And Computer Vision》一书中引用的DFT实现函数：

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. 该常数因子的来源是什么，尤其是其中的atan(1)的含义？
2. 这个常数因子的作用是什么？

问题解答 #

1. 常数因子的来源与atan(1)的含义

先回忆DFT的核心定义公式：
$$Y[m] = \sum_{k=0}^{n-1} X[k] \cdot e^{-j\frac{2\pi}{n}km}$$
公式里的复指数项$e^{-j\frac{2\pi}{n}km}$是DFT的核心旋转因子，而代码里的常数因子正是为了生成这个旋转因子。

拆解代码里的计算式：atan(1.0)的数学结果是$\frac{\pi}{4}$（因为$\tan(\frac{\pi}{4})=1$），所以$\frac{\pi}{4} \times -8 = -2\pi$，最终整个常数因子就是$\frac{-2\pi}{n}$。代码把这个值赋值给了复数tmp的虚部（cmplx(0.0, ..., &tmp)），所以tmp是纯虚数：$0 + j\frac{-2\pi}{n}$。

作者用atan(1.0)*8来表示$2\pi$，本质是绕开直接使用标准库的M_PI宏——有些早期编译环境或者特定平台可能没有定义这个圆周率常量，用三角函数计算的方式可以保证代码的兼容性。

2. 常数因子的作用

这个$\frac{-2\pi}{n}$是DFT旋转因子的核心组成部分，代码里用它来生成DFT公式要求的复指数项：

• 先计算tmp * k * m，也就是$j\frac{-2\pi}{n} \times k \times m$
• 再调用cexp函数计算这个纯虚数的指数，根据欧拉公式$e^{j\theta} = \cos\theta + j\sin\theta$，这里的$\theta = \frac{-2\pi}{n}km$，所以计算结果就是：
  $$e^{j\frac{-2\pi}{n}km} = \cos\left(\frac{2\pi}{n}km\right) - j\sin\left(\frac{2\pi}{n}km\right)$$
  这完全匹配DFT公式中的复指数项。

简单来说，这个常数因子的作用就是生成DFT所需的旋转因子，用来给时域信号的每个样本按不同频率分量（对应不同的m）加权，最终通过求和得到频域结果Y[m]。

我使用的其他辅助函数 #

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

内容的提问来源于stack exchange，提问作者badcode