离散傅里叶变换学习——从一维到海洋模拟（一）最近在朋友的点拨下，以及通过网上查阅的一些资料来看，实现了一维离散傅里叶变换

离散傅里叶变换的实现

最近在朋友的点拨下，以及通过网上查阅的一些资料来看，实现了一维离散傅里叶变换到二维离散傅里叶变换，以至于到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次：

x[n] = cos(2\pi\frac{n}2),n = 0,1,...,39

然后通过DFT可以知道它在40次采样周期内，震动了几个周期。算法的处理很暴力：

首先，选取40个长度为40个点的基信号，它们分别长这样：

cos(2\pi\frac{0n}{40}),cos(2\pi\frac{n}{40}),cos(2\pi\frac{2n}{40})...cos(2\pi\frac{39n}{40})

第一个，40次采样内振动0个周期： $cos(2\pi\frac{0n}{40})$ ，即常值:

第二个，40次采样内振动1个周期： $cos(2\pi\frac{n}{40})$

以此类推，一直到40个采样内振动39个周期。

接下来，对于上述每个基信号，判断它们跟原信号的相关程度，就是用他们在同一点的函数值相乘，并对结果求和（向量的内积），即如下的公式：

correlation(x,y) = \Sigma_k x[k]y[k]

这个值越大，则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：

如果用频域信号替换原信号，则：

X_k = \sum_{n=0}^{N-1}x[n]cos(2\pi \frac{kn}{N}), (N=40)

问题来了，虽然貌似联系很紧密，但这怎么跟DFT的公式长得不一样。。。DFT的公式应该是这样的：

X_k = \sum_{n=0}^{N-1} x[n]e^{-2\pi j \frac{kn}{N}}

用欧拉公式展开，我们得到的时：

X_k = \sum_{n-0}^{N-1} x[n]cos(2\pi \frac{kn}{N}) - j\sum_{n=0}^{N-1}x[n]sin(2\pi \frac{kn}{N})

这又是为什么呢？这是因为，对于一个信号，如果只跟余弦函数比较，会损失一些信息，比如相位。如果原信号有一些相位偏移， $x = cos(2\pi \frac{2n}{40} + \frac{\pi}{4})$

对这个函数同样按照上面的方法计算频域，结果会有些不一样:

X_2 = X_38 = \sum_{n=0}^{39}cos(2\pi\frac{2n}{40} + \frac{\pi}{4})cos(2\pi \frac{2n}{40}) = 10\sqrt2

如果再找一个信号y, 没有相位偏移，而是把幅值砍到 $\frac{\sqrt2}{2}$ ，即:

y = \frac{\sqrt 2}{2} cos(2\pi \frac{2n}{40})

那么这个信号的DFT结果：

Y_2 = Y_38 = \sum_{n=0}^{39}\frac{\sqrt 2}{2}cos(2\pi\frac{2n}{40})cos(2\pi\frac{2n}{40})=10\sqrt2

跟x信号的记过一模一样，这样就由于损失信息，无法通过频域恢复信号了。

解决方法是另选一组以正弦函数（实际上选了负正弦）为基准的“基信号”，即 $-sin(2\pi\frac{0n}{40})$ 到 $-sin(2\pi\frac{39n}{40})$ ，计算另一组原信号与正弦基的相关系数，这两组系数一起作为DFT的最终结果。而复数只是一个工具，用来方便地同时存储两组计算结果。当然它还有一个好处就是能够比较直观地表现出模和相位。

选负正弦还是正弦做基信号其实无所谓，只是最后的结果算出来的相位反一下而已，幅值是一样的。如果一个信号跟某频率余弦和负正弦的相关系数分别为a,b，那么代表这个信号差不多型如：

acos(\frac{kn}{N}) - bsin(\frac{kn}{N})

根据高中数学可以求得其模为 $\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