FFT

作为一名做音频研究的人，如果对傅里叶变化不是很精通，那就不专业了。傅里叶变换实际上已经诞生了很长时间了，而且经过了实践的考验，现在几乎被烧到了所有的嵌入式数字信号处理设备上。由于楼主对硬件不是非常了解，但是猜测一下，估计傅里叶变换的快速算法FFT已经被烧到了现在个人电脑的主板上了。正是由于它的高效而广泛的应用，被评为20世纪最伟大的10个算法之一。

一、数学解释

1、多项式的表示——系数表示和点值表示

作为这么一个伟大的算法，后面的数学思想当然是非常有内涵了。

这要从多项式说起。如何来对一个多项式进行标示呢？一个非常直观的方法是使用多项式的系数，假设一个$N$阶多项式，可以使用下面的向量来表示，

$B=[b_0,b_1,...,b_{n-1}]$

还有一个不是非常直观的表示方式，是使用点值来进行表示。比如一个特定的N阶多项式就可以使用相异的N个点来进行表示。比如

$D=[{x_0,y_0},{x_1,y_1},...,{x_{n-1},y_{n-1}}]$

那么问题来了，两者的关系是什么呢？实际上两者是等价的。

 

假设我们知道N个点值对，那么我们就可以将其表示为上面的矩阵表示方式。从线性代数的角度我们知道上述？？是否有解取决于左侧的方阵的阶是否是满的，如果是也就是它是否是非奇异矩阵，判断方阵是否非奇异的方式是使用它的行列式。上述方阵的范德蒙矩阵为，

$|X|=\sum_{i<j}(x_{i}-x_{j})$

根据这样的推导定义我们知道，如果X的行列式不是零，则充分必要条件是N个点对互异。这一点在一个连续函数当时非常容易满足的，我们可以随意找到这样的N个点。这就证明了从上述的点值表示，一定可以找到所对应的系数表示。我们可以发现实际上这是一个对多的关系。本质上已知点值求系数是一个差值的过程，已知系数求点值是一个求值的过程，显然在这里，求值运算和差值运算时互逆的。

 

2、系数表示和点值表示的利弊

1) 系数表示是最直观的表示方式，而且从求值的角度看系数表示更加方便。

2) 但是对于某些多项式的运算使用点值表示来的更加方便，比如多项式的乘法运算，使用点成即可，复杂度是O（n）的。对于加法运算也是相似的道理，对于点值加法来说，复杂度也是O(N)的。

因此在涉及多项式有些运算时，一个解决方案是，先将两个多项式转化为点值表示，然后在对转化之后的点值表示进行点乘，然后在将点值表示的成绩反变换到系数表示。这里的复杂度取决于两种表示之间进行转换所需要的时间。

 当然对于多项式的加法运算来说，使用系数表示复杂度同样为$O(N)$,但是系数表示进行乘法运算复杂度必然是$O(N^2)$的，因为其中一个多项式的每一项要乘以另一个多项式的每一项，两层for循环。

 

3、单位复根

从上面的分析当中，我们知道，系数表示和点值表示之间是一对多的关系，所以在进行从系数表示想点值表示的转换时，我们有非常多的选择的余地。选择合适的求知点，可以将计算的复杂度降低。 在实际当中选择的求值点是单位复根。

所谓单位复根是满足下面的等式的复数。

$w^n=1$，这个方程的根总共有n个，并且等间距的排列在复平面上，分别为$e^{\frac{2\pi * i * k}{n}}, k=0,1,2,...,n-1$

为了方便表示，我们定义主n次单位复根为

$w_n=e^{\frac{2\pi * i }{n}}$

相应的n个单位复根表示为，

$w_{n}^{0},w_{n}^{1},w_{n}^{2},...,w_{n}^{n-1}$

关于单位复根有几个非常重要的引理

一、相消引理

$w_{nl}^{kl}=w_{n}^{k}$

证明：

根据$w_{n}^{k}$的定义，上面的等式可以写成，

$w_{nl}^{kl}=e^{\frac{2\pi kl i}{nl}}=e^{\frac{2\pi k i}{n}}=w_{n}^{k}$

问题得证

二、折半引理

根据引理一，我们可以得到一个结论

$w_{n}^{\frac{n}{2}}=-1$

证明：

$w_{n}^{\frac{n}{2}}=e^{\frac{2 * \pi * i * \frac{n}{2}}{n}}$

即

$w_{n}^{\frac{n}{2}}=e^{\pi * i}$

转化为三角函数表示，即

$w_{n}^{\frac{n}{2}}=e^{\pi * i}=cos(\pi)+i*sin(\pi)=-1$

问题得证

那么问题来了，假设我们要对所有的n个n次单位复根进行平方，我们可以得到多少个不同的单位复根呢？

1）首先所得到的单位复根一定是$\frac{n}{2}$次的。

$(w_{n}^{k})^{2}=e^{\frac{4 * \pi * i * k}{n}}=e^{\frac{2 * \pi * i * k}{0.5*n}}$

也就是

$w_{\frac{n}{2}}^{k}$

2）那么问题来了这样的单位复根有多少个呢？显然这里的n个$\frac{n}{2}$次单位复根当中是有重复的，每一个都冲了一次。

那么对应关系是什么呢？

$(w^{k+\frac{n}{2}}_{n})^{2}= (w^{k}_{\frac{n}{2}})^{2} *(w^{\frac{n}{2}}_{n})^{2}$

由于

$w^{\frac{n}{2}}_{n}=-1$

所以上式转化为

$(w^{k+\frac{n}{2}}_{n})^{2}= (w^{k}_{\frac{n}{2}})^{2}$

4、DFT

Discrete Fourier Transform (DFT) 离散傅里叶变换的定义，实际上可以从多项式求值的角度来定义，假设我们知道多项式的定义，现在我们想求多项式在n个不同的点的值（X,Y），当然这里的X是单位复根（单位复根在数字信号当中有特殊的意义），这里的Y就是离散傅里叶变换的结果。所以本质上，离散傅里叶变换可以看做是多项式求值。

5、FFT

上述算法如果不进行完善，那么算法的复杂度将会是$O(N^2)$, 想一想吧，对于每一个单位复根都会带入多项式系数表示的多项式当中求出对应的Y，复杂度为$O(N)$,总共有N个单位复根，所以总体的复杂度为$O(N^{2})$

那么有没有更快的呢？FFT是一个会更快的算法。

主要的算法思想是分治的思想。(在上面已经提到了，这个算法的改进的一个关键是对求值点的选择，我们知道多项式的系数表示和点值表示之间的一对多的关系，我们在求值点的选择上有非常大的灵活性，合理的选择求值点可以大幅度的降低计算复杂度。这里的求值点就是上面介绍的N个N次单位复根，利用上面证明的两个性质来进行加速。)括号中的话，应该是有问题的，实际上单位复根是傅里叶变换定义的一部分，而不是FFT的特点，但是单位复根的某些特点有利于算法的实现，FFT的主要思想应该是分治算法。许多地方都在说，利用了旋转因子的对称性和周期性，不知所云！

 

废话不说了，上算法！

假设要计算的多项式为：

$F(x)=a_{0}+a_{1}*x+a_{2}*x^{2}+...+a_{n-1}*x^{n-1}$

我们将其按照每一项的下标的奇偶性，将其分解为两个多项式，即

$Odd(x)  = a_{1}+a_{3}*x+...+a_{n-1}*x^{\frac{n}{2}-1}$

$Even(x) = a_{0}+a_{2}*x+...+a_{n-2}*x^{\frac{n}{2}-1}$

那么$F(x)=Odd(x^{2})+x*Even(x^{2})$

显然$Odd(x)$和$Even(x)$的规模是原来的$\frac{1}{2}$,在原来的$F(x)$中是$N$次单位复根，那么这里是$\frac{N}{2}$次单位复根，根据上面的引理我们知道N个N次单位复根分别进行平方式实际上得到的就是$\frac{N}{2}$个$\frac{N}{2}$次单位复根，只是他们每一个重复了一次。也就是在小规模的多项式对$x$进行计算的结果实际上就是在大规模的多项式计算$x^{2}$的结果。这是一个非常关键的地方，这在直观上，我们发现，工作量减少了一半，其实不止是这些，还有哪些？？？这里需要思考一下。

我们可以举一个非常明显的例子，将$Odd(x_{i})$简记为$Odd_{i}$假设我们现在已经计算得到了$[Odd_{0},Odd_{1},Odd_{2},...,Odd_{\frac{n}{2}-1}]$和$[Even_{0},Even_{1},Even_{2},...,Even_{\frac{n}{2}-1}]$

下面我们通过上面的组装公式将其进行组装为$F(x)$

由于$F(x)$实际上是针对的N个N次单位复根，而Odd和Even是针对的$\frac{N}{2}$个$\frac{N}{2}$次单位复根，但是我们知道对前者进行平方之后有下面的对应关系

也就是说前半部分和后半部分是相等的，显然我们在组装加大规模的多项式的时候，只需要考虑前半部就可以了。那么前半部的N次单位复根和后半部分的N次单位复根是啥关系呢？上实际上有个引理的证明中就已经提到了，这里在重复一遍。

$w_{n}^{k+\frac{n}{2}}=w_{n}^{k}$,因为

$w_{n}^{\frac{n}{2}}=-1$

那么实际上，上面的单位复根前后是呈现相反数的关系。那么问题就好解决了。

我们现在只关于一次组装过程的伪代码

assemble-FFT(Odd,Even)

　　n=length(Odd);

　　res=zeros(1,2*n);

　　$w_{n}=e^{\frac{2* \pi * i}{2*n}}$

　　$w=1$

　　for k=0:n

  　　res(k)=Even(k)+w*Odd(k)

  　　res(k+n)=Even(k)-w*Odd(k)

  　　$w=w*w_{n}$

　　end

　　return res

二、数字信号解释

 $Y=FFT(X)$

     通过$magnitude=abs(Y)$我们可以得到每种频率分量所对应的能量值，由于存在对称性，我们只会考虑折半之后的频带能量分布。

　  ONE：　　

     下面我们来证明一下频谱分布的共轭性质，我们知道。

　　$Y_{k}=\sum_{n=0}^{N-1}x_{n}e^{i*n*k*\frac{2* \pi}{N}}$

　　将$Y_{k}$的共轭记为$\bar{Y_{k}}$,我们需要证明$Y_{N-k}=\bar{Y_{k}}$，下面我们证明之。

　　$Y_{N-k}=\sum_{n=0}^{N-1}x_{n}e^{i*n*(N-k)*\frac{2 * \pi}{N})}$，可以化简为

　　$Y_{N-k}=\sum_{n=0}^{N-1}x_{n}e^{-i*n*k*\frac{2 * \pi}{N}}e^{i*n*2*\pi}$

　　又由于$e^{i*n*2* \pi}=cos(2*n*\pi)-i*sin(2*n*\pi)=1$,所以上式可以化简为

　　$Y_{N-k}=\sum_{n=0}^{N-1}x_{n}e^{-i*n*k*\frac{2 * \pi}{N}}=\bar{Y_{k}}$

　　问题得证！

　　如果对上述频谱进行求取绝对值，从而得到相对应的频谱能量，由于存在前后对称的共轭性质，因此能量的分布也是对称的。所以我们一般会去掉上半部分的高频部分，值考虑低频部分，因为两者是对称的，对于高频的考虑可以使用低频对称的转换过去即可。

　　TWO：

　　对于频谱求取能量又应该怎么解释呢？

　　实际上最根本的原因是在于，我们刚开始对于原始信号的转化，这要从数字信号本身谈起，我们试图使用离散的三角函数来表示现在的信号，而且从数学的角度我们知道，任何函数都可以使用三角级数来加以逼近。大体的意思可以形式化为如下所示。

　　$X=b_{0}+\sum_{n=0}^{N-1}(a_{n}sin\frac{2*\pi*n}{N}+b_{n}cos\frac{2*\pi *n}{N})$

通过欧拉公式，我们可以将上市转化为复数形式

　　$X=c_{0}+\sum_{n=1-N}^{N-1}c_{n}e^{\frac{n*2*\pi}{N}}$

　　这里的$c_{n}$实际上就是上面的傅里叶变换的结果，具体值得是什么呢？

　　$c_{n}=\frac{1}{2}(b_{n}-i*a_{n})$

　　因此$c_{n}$的平方实际上就是前面的表示的同一种频率的两个三角函数的系数之和，这可以看做这种频率分量的能量。