离散傅里叶变换的实现
最近在朋友的点拨下，以及通过网上查阅的一些资料来看，实现了一维离散傅里叶变换到二维离散傅里叶变换，以至于到FFT的实现及相应的逆变换。对傅里叶变换这个很长时间以来都没有理解的东西，有了一个深刻的认识。所以就想总结一下其中的原理以及具体的实现过程。
本文将会从一维离散傅里叶变换开始，逐步讲解到FFT的实现及相应逆变换的实现方法。以及FFT实现时使用的蝶形变换的具体操作方法。其中，还会罗列出所参考的资料。
一维傅里叶变换
在实际操作之前，一定要有一个重要的认知，那就是傅里叶变换本质上是把一个函数转换成用另一种形式来表示。这就是经常说的，时域信号转成频域信号。但这都不重要，数学公式上来看，傅里叶公式如下，即等号左边的f(x)可以表示成等号右边的一个和式，这个公式的通俗解释就是，一个函数可以表示成正弦函数和余弦函数的和。而其中an和bn表示的就是一个权值。
an和bn也是可以通过公式计算得到的，公式如下，可以看到是两个积分式。
总结这两个公式，可以获得一个认知：
- 一个函数可以表示成不同频率的正弦函数和余弦函数的加权和。权值可以通过公式计算。
- 只要确定了权值an和bn，结合第一个公式，通过对不同频率的正弦函数和余弦函数进行加权求和，我们可以算出来对应f(x)的值。
这里面有一个容易被忽视的点，那就是不同频率的正弦函数和余弦函数的组合，是固定的。所以最后只要确定权值an和bn即可，至此就可以知道，对一个函数进行傅里叶变换，要做的事情就是确定an和bn的值。通过手动计算当然很难，但是借助计算机却是可以的，对于原函数很难确定的积分公式，计算机也只能进行离散地计算以获得一个趋近于正确答案的结果，只要误差允许范围内，这没什么大问题。
DFT做的事情，主要是在长度为N的离散信号中，针对k=(0,1,2…)，分别找出在长度N内振动k个周期的三角波分量的权值。举个例子，针对某个余弦信号，在两个周期内采样40次：
然后通过DFT可以知道它在40次采样周期内，震动了几个周期。算法的处理很暴力：
首先，选取40个长度为40个点的基信号，它们分别长这样：
第一个，40次采样内振动0个周期：$cos(2\pi\frac{0n}{40})$，即常值:
第二个，40次采样内振动1个周期：$cos(2\pi\frac{n}{40})$
以此类推，一直到40个采样内振动39个周期。
接下来，对于上述每个基信号，判断它们跟原信号的相关程度，就是用他们在同一点的函数值相乘，并对结果求和（向量的内积），即如下的公式：
这个值越大，则x[k]与y[k]越像。于是DFT把$cos(2\pi\frac{0n}{40})$到$cos(2\pi\frac{39n}{40})$这40个基函数与当前函数$cos(2\pi\frac{2n}{40})$比较了一下，发现$cos(2\pi\frac{2n}{40})$和$\cos(2\pi\frac{38n}{40})$长得最像。
这也和那显然，因为$cos(2\pi\frac{38n}{40}) = cos(2\pi n-2\pi \frac{38n}{40}) = cos(2\pi \frac{2n}{40})$
下面，如果我们把这40次每次比较的correlation值记下来，就得到了原信号在每个频率上的分量大小。就得到了所谓原信号的频域X：
如果用频域信号替换原信号，则：
问题来了，虽然貌似联系很紧密，但这怎么跟DFT的公式长得不一样。。。DFT的公式应该是这样的：
用欧拉公式展开，我们得到的时：
这又是为什么呢？这是因为，对于一个信号，如果只跟余弦函数比较，会损失一些信息，比如相位。如果原信号有一些相位偏移，$x = cos(2\pi \frac{2n}{40} + \frac{\pi}{4})$
对这个函数同样按照上面的方法计算频域，结果会有些不一样:
如果再找一个信号y, 没有相位偏移，而是把幅值砍到$\frac{\sqrt2}{2}$，即:
那么这个信号的DFT结果：
跟x信号的记过一模一样，这样就由于损失信息，无法通过频域恢复信号了。
解决方法是另选一组以正弦函数（实际上选了负正弦）为基准的“基信号”，即$-sin(2\pi\frac{0n}{40})$到 $-sin(2\pi\frac{39n}{40})$，计算另一组原信号与正弦基的相关系数，这两组系数一起作为DFT的最终结果。而复数只是一个工具，用来方便地同时存储两组计算结果。当然它还有一个好处就是能够比较直观地表现出模和相位。
选负正弦还是正弦做基信号其实无所谓，只是最后的结果算出来的相位反一下而已，幅值是一样的。如果一个信号跟某频率余弦和负正弦的相关系数分别为a,b，那么代表这个信号差不多型如：
根据高中数学可以求得其模为$\sqrt{a^2+b^2}$，相对余弦的相位为$arctan(\frac{b}{a})$，这与复数$a+bi$的模和相位是相同的，因此DFT的公式相当于同时把x[n]做了跟N个余弦基、N个负正弦基的比对，结果用N个复数存储。如果想要看某个频率的相位和模，就看对应复数的相位和模。
我们再来看看上面有相位偏移的那个例子：
原信号：$x[n] = cos(2\pi\frac{2n}{40} + \frac{\pi}{4})$
与余弦比对： $X_2 = \sum_{n=0}^{39}cos(2\pi\frac{2n}{40}+\frac{\pi}{4})cos(2\pi\frac{2n}{40})=10\sqrt{2}$
与负正弦比对： $X_2=\sum_{n=0}{39}cos(2\pi\frac{2n}{40} + \frac{\pi}{4})sin(-2\pi\frac{2n}{40})=10\sqrt{2}$
在40个点内振动两个周期这个频率上，其DFT的结果为$10\sqrt{2}+10\sqrt{2}j$
其模为20，与其相位偏移前相同
其相位也为$\frac{\pi}{4}$，也没有问题
如果将原信号变为$x[n]=sin(2\pi\frac{2n}{40}+\frac{\pi}{4})$,会求得该频率DFT结果为$10\sqrt{2}-10\sqrt{2}j$，求得其相位为$-\frac{\pi}{4}$。因此，根据DFT结果求得的相位是相对余弦信号的相位。
一维离散傅里叶变换的逆变换
一维离散傅里叶逆变换的公式如下：
具体的公式推导我不太理解，所以就不讲了。其中，要还原的目标，即原函数的值。则是由经过DFT变换得到的结果，即一组各个频率正交基的系数。当初在大学期间，看这个公式时，有一个很疑惑的点是左边的值是一个实数，右边的公式中是DFT的结果，看着也像是实数，然后又是一个e的指数形式，很困惑。这个疑惑不知道大家有没有，这里就特别提一下，在这个公式中，其实是复数。是DFT的结果，也是复数，而，由欧拉公式可以知道，也是复数。所以这个公式中右边部分是复数相乘并求和的结果，左边自然是也是复数。最后的得到的复数，实部就是我们想要的结果，虚部会计算变为0
C++实现一维离散傅里叶变换
二维离散傅里叶变换及其逆变换
二维离散傅里叶变换的公式有一维离散傅里叶变换相比，就是多了一层求和，公式如下：
逆变换的公式如下：
在代码的实现上，就是从两层循环，变成了4层循环。我用一张图片作为例子，来实现二维离散傅里叶变换及其逆变换。程序加载图片但只对它的R通道进行处理，将傅里叶变换的结果输出成图像。并将逆变换的结果输出成图像，和原图像进行对比。
程序需要依赖stb_image.h及stb_image_write.h图形库，具体的使用方法可以直接参考代码的示例。图片用的是缩小版本的lena的照片，不过只取了R通道。
首先实现一个Image的类，用于对图形数据的提取，下面的图给出类的声明、构造和析构方法。实现加载图片的逻辑。加载出来的图片，用一维浮点数组保存，且拆分成RGB三个通道，索引时根据x和y计算整数索引来取数图片中的值。在这篇示例中，只需要关注R通道即可。
在main函数里，使用刚刚定义的类实现加载和保存的逻辑。
然后我们需要在加载图片之后，实现DFT方法，并将结果保存为dft_result.png。
如果将DFT的模输出成图像，是下面这张图
为了让其显像，做一个取对数的操作，是下面这张图。
这也和我们常见频谱图像有区别，而如何从上面这张图变换得到下面这张图，下文再进行介绍。
接着实现IDFT方法，并将结果保存为output.png。
可以看到，输出的结果与原图一致。如下图是转换过程，最左边是原图，最右边是复原图，中间是DFT的结果。肉眼无法分辨出差别，而实际上，两张图的差值是0。
而对于下面这张频谱图像的获得，实际上是对DFT的结果做了一些操作。并且对图像进行处理，使其能够显像，这在冈萨雷斯的《数字图像处理》一书中也有体现。代码就不贴了，可以看Image类的GenerateSpectrum方法。
所以，对于DFT结果和频谱图像之间的关系，就讲清楚了。其实如果只是写程序实现这个过程，频谱图像是不必要关心的，但是频谱图像对于。更应该关注的是DFT的结果是如何存储的，以及如何复原原图像。下面要介绍的FFT利用一些定理对DFT加速，得到的结果和DFT一致。二维DFT的代码也上传到仓库里。
参考资料
如何通俗地解释什么是离散傅里叶变换？ - 知乎 (zhihu.com)
