引子

本文就二维离散傅里叶变换(2-D Discrete Fourier Transform)的计算机实现做一些分享，若有不正确的地方，还请多多指教。
傅里叶变换可以将图像从空间域(Spatial Domain)到频率域(Frequency Domain)，这对于受到周期性噪声(periodic noise)图像的修复有着很好的效果。

好啦，正文开始。
我将按照以下的顺序讲讲我的理解：

复数计算是如何体现在计算机上
初步的DFT算法实现 --- 时间复杂度： $T(n^4)$
更快的DFT算法 --- 时间复杂度： $T(n^3)$

复数计算

这是我初见傅里叶变换困扰了很久的东西（高中数学没学好 ಠ_ಠ）
其是也不是很难，只用先定义一个结构体：

typedef struct {
    double real;
    double imag;
} Complex;

即在DFT之前，IDFT(Inverse Discrete Fourier Transform)之后，图像只有实部。
而在这两者(DFT, IDFT)之间，图像有实部和虚部。

再贴一些复数计算公式，后面会用到：

加法 $(a+bi) +(c+di)=(a+c)+(b+d)i$
乘法 $(a+bi) \times(c+di)=(ac-bd)+(bc+ad)i$

初步的傅里叶变换

Discrete Fourier Transform

先来放个公式。
$F(u,v)=\frac{1}{MN}\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)e^{-j2\pi (ux/M+vy/N)}$

又因欧拉公式(Euler's formula)。
$e^{j\theta}=cos\theta+jsin\theta$

可得，
$F(u,v)=\frac{1}{MN}\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)(cos2\pi (ux/M+vy/N)-jsin2\pi (ux/M+vy/N))$

这里，前半部分的 $cos$ 放入complex.real中；后半部分的 $jsin$ 放入complex.imag中。

具体体现在代码里，是这样：

Complex *DFT(Image *image) {
    
    // Original Image
    int width = image->Width;
    int height = image->Height;
    unsigned char *imageData = image->data;
    
    // Malloc
    Complex *complexArr;
    complexArr = (Complex *)malloc(sizeof(Complex) * width * height);
    
    // Centralization
    int tempImage[width*height];
    for (int i = 0; i < width; i++) {
        for (int j = 0; j < height; j++) {
            tempImage[i * width + j] = int(round(imageData[i * width + j] * pow(-1, i + j)));
        }
    }
    
    // Loop the whole image
    for (int u = 0; u < width; u++) {
        for (int v = 0; v < height; v++) {
            // DFT Function
            for (int x = 0; x < width; x++) {
                for (int y = 0; y < height; y++) {
                    // theta
                    double theta = double(-2)*PI*(u*x/width + v*y/height);
                    // Real part
                    complexArr[u * width + v].real += tempImage[x * width + y] * cos(theta);
                    // Imaginary part
                    complexArr[u * width + v].imag += tempImage[x * width + y] * sin(theta);
                }
            }
            complexArr[u * width + v].real /= width * height;
            complexArr[u * width + v].imag /= width * height;
        }
    }
    
    return complexArr;
}

这里输入的图像格式为.pgm，输出的是Complex数组。
并且在DFT前，增加了中心化的处理，便于后面在频率域的观察。

好了，现在我们已经将图像从空间域转化到了频率域，可以在这里对图像进行很多操作，本文将只提及最基本的可视化操作。即，Fourier Spectrum (傅里叶谱)。

还需要提及的一点是，实部和虚部都包含着原图像的信息，在频率域做处理时，需要同时处理这两个部分（基本上都是同时乘以某个滤波器）。
后面做逆傅里叶变换时也是如此，而不是使用 Spectrum 后的结果。(这也是我曾经纠结过的问题 /_ \ )

Fourier Spectrum

现在来可视化一下，有两个可以将频率域可视化的方法：

Fourier Spectrum
Phase Angle

这里我们用第一个，因为第二个表达的是相角，其产生的视觉直观信息很少。
看看 Spectrum 的公式：
$|F(u,v))|=[R^2(u,v)+I^2(u,v)]^{1/2}$

然后是代码：

double *FourierSpectrum(Complex *imageArr, int width, int height) {
    
    double *outputArr = (double *)malloc(sizeof(double) * width * height);

    int rows = height;
    int cols = width;

    int enhance = 100;

    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            double real = imageArr[i * rows + j].real;
            double imag = imageArr[i * rows + j].imag;
            outputArr[i * rows + j] = sqrt(real * real + imag * imag) * enhance;
        }
    }

    return outputArr;   
}

这里将图像强度进行了提高，更便于对频率域的观察。
返回了浮点字符串，便可以保存为图片输出。

Inverse Discrete Fourier Transform

在这里，如果我们前面没有对频率域的做任何改变，那么经过逆傅里叶变换之后，我们将会得到和原图一模一样的图像。
再来看看公式，其实跟正傅里叶的几乎一样。
$f(x,y)=\sum^{M-1}_{u=0}\sum^{N-1}_{v=0}F(u,v)e^{j2\pi (ux/M+vy/N)}$

变换一下：
$f(x,y)=\sum^{M-1}_{u=0}\sum^{N-1}_{v=0}F(u,v)(cos2\pi (ux/M+vy/N)+jsin2\pi (ux/M+vy/N))$

只是一个正负号的差别。简单。

但需要注意的是，这里的 $F(x,y)$ 是复数，因此需要用到我们最开始提到的复数乘法（不记得的往上翻翻）。因为IDFT之后便是实数数据了，我们只需关注计算结果的实部部分即可。

再细看， $F(u,v)\to$ Complex.real, Complex.imag
即， $(real+imag)\times(cos2\pi (ux/M+vy/N)+jsin2\pi (ux/M+vy/N))$
$=real\times cos2\pi (ux/M+vy/N)-imag\times sin2\pi (ux/M+vy/N))$

然后我们就能看懂代码了。

double *IDFT(Complex *complexArr, int width, int height) {

    // Malloc
    doube *outpurArr = (double *)malloc(sizeof(double) * width * height);

    // Loop the whole image
    for (int x = 0; x < width; x++) {
        for (int y = 0; y < height; y++) {
            // IDFT Function
            for (int u = 0; u < width; u++) {
                for (int v = 0; v < height; v++) {
                    double theta= 2 * PI * ((double)x * u / width + (double)y * v / height);
                    outpurArr[x * width + y] += complexArr[u * width + v].real * cos(theta) - complexArr[u * width + v].imag * sin(theta);
                }
            }
        }
    }
    
    return outpurArr;
}

至此，便是正向、反向傅里叶变换的完整过程了，只需要加上一些基本图像读取、保存的操作，就可以愉快的运行了。

但我们看代码，或者运行的时候就能很明显的感觉到，这太慢了，跑一个稍大一些的图像就需要很长的时间。因为这 4 个嵌套的 for 循环，使时间复杂度达到了 $T(n^4)$ ，我们能不能再快一点呢？

这里将一个不需要快速傅里叶变换的方式，便可以将时间复杂度降到 $T(n^3)$ ，这是一个很大的进步，大概能加速100倍。
接着...

更快的傅里叶变换

这里需要引入核(Kernel)的概念。书中是这样说的：

A particularly important class of 2-D linear transforms, denoted $T(u,v)$ , can be expressed in the general form:
$T(u,v)=\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)r(x,y,u,v)$
where $f(x,y)$ is the input image, $r(x,y,u,v)$ is called the forward transformation kernel.

对于傅里叶变换来说：

The 2-D Fourier transform discussed in Example 2.11 has the following forward and inverse kernels:
$r(x,y,u,v)=e^{-j2\pi (ux/M+vy/N)}$ and $s(x,y,u,v)=e^{j2\pi (ux/M+vy/N)}$

这里是不是就有些眼熟了？并且我们可以证明，Fourier kernels 是可分离的(separable)且对称的(symmetric)。这就很重要了，意味着我们可以把一个 $T(n^4)$ 的计算拆分成两个 $T(n^3)$ 的子计算，体现在二维图像就是，先算行(列)的DFT，再计算列(行)的DFT。

公式是这样的：
$T(u,v)=\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}f(x,y)r(x,y,u,v)$
$=\sum^{M-1}_{x=0}r_1(x,u)\sum^{N-1}_{y=0}f(x,y)r_1(y,v)$
$=\sum^{M-1}_{x=0}T(x,v)r_1(x,u)$
where, $T(x,v)=\sum^{N-1}_{y=0}f(x,y)r_1(y,v)$

在傅里叶变换中， $T(u,v)$ 就是 $F(u,v)$ 。并且在逆傅里叶变换中也是类似的。

DFT 和 IDFT 一起贴上吧 :)

Complex *DFT(double *imageArr, int width, int height) {

    Complex *t = (Complex *)malloc(sizeof(Complex) * width * height);
    Complex *g = (Complex *)malloc(sizeof(Complex) * width * height);

    int rows = height;
    int cols = width;

    printf("Starting DFT ...\n");

    // DFT for cols
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            for (int m = 0; m < rows; m++) {
                double theta = double(-2) * PI * (i * m) / rows;
                t[i * rows + j].real += imageArr[m * rows + j] * cos(theta);
                t[i * rows + j].imag += imageArr[m * rows + j] * sin(theta);
            }
        }
    }

    // DFT for rows
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            for (int n = 0; n < cols; n++) {
                double theta = double(-2) * PI * (j * n) / cols;
                g[i * rows + j].real += t[i * rows + n].real * cos(theta) - t[i * rows + n].imag * sin(theta);
                g[i * rows + j].imag += t[i * rows + n].real * sin(theta) + t[i * rows + n].imag * cos(theta);
            }
        }
    }

    // Divide by image size
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            g[i * rows + j].real /= rows * cols;
            g[i * rows + j].imag /= rows * cols;
        }
    }

    printf("DFT finished\n");

    return g;
}

And...

double *IDFT(Complex *imageArr, int width, int height) {

    Complex *t = (Complex *)malloc(sizeof(Complex) * width * height);
    double *g = (double *)malloc(sizeof(double) * width * height);

    int rows = height;
    int cols = width;

    printf("Starting IDFT ...\n");

    // IDFT for cols
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            for (int m = 0; m < rows; m++) {
                double theta = double(2) * PI * (i * m) / rows;
                t[i * rows + j].real += imageArr[m * rows + j].real * cos(theta) - imageArr[m * rows + j].imag * sin(theta);
                t[i * rows + j].imag += imageArr[m * rows + j].real * sin(theta) + imageArr[m * rows + j].imag * cos(theta);
            }
        }
    }

    // IDFT for rows
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            for (int n = 0; n < cols; n++) {
                double theta = double(2) * PI * (j * n) / cols;
                g[i * rows + j] += t[i * rows + n].real * cos(theta) - t[i * rows + n].imag * sin(theta);
            }
        }
    }

    printf("IDFT finished\n");

    return g;
}

傅里叶变换的计算方式和第一种是一样的，如果有不清楚的地方，可以看看上面那部分写的，相对比较详细一些。

尾巴

内容有点多，Markdown也没有拼写纠错，可能会有失误的地方，还请指出。

后面我可能还会更新一些，图像处理的学习笔记，挖坑ing...

也很感谢你们能看到这里，祝你们在频率域玩的愉快！¬_¬

彻底理解傅里叶变换 2-D Discrete Fourier Transform -- C/C++ 实现