傅里叶变换，合数的FFT实现和源码
题外话：
这是很久之前写的代码，当时是因为想回顾一下以前的小霸王游戏，所以在去了解一个NES游戏的过程中，遇到了声音处理上的一些问题，不得以才看了一下傅里叶。虽然很多大神已经写了一些通俗易懂的傅里叶，但是我还是想从代码层面来说一些我理解的傅里叶。不想看介绍请直接跳到本文末c++关于快速傅里叶的实现代码，文中会有实现介绍。希望文本对你有所帮助。
一些必要的概念：
学习概念呢，是为了对后续更深的知识进行一些自我挖掘，从而不断完善知识系统，让世界趋于完美。
傅里叶变换是一种线性积分变换，用于信号在时域（或空域）和频域之间的变换，在物理学和工程学中有许多应用。
时域，通俗上说就是我们的所见所闻，也就是我们最容易理解的世界随时间的变化而发生的最直观的改变。严格上说，是信号与时间的对应关系。
频域 ，相对于时域而言，这个就是我们不容易通过直接观察得到的，但是却是一种最简单的方式可以让我们用来描述这个世界的运动规律。严格上说，频域是时间与状态的对应关系。
如果不想看下面我的一些我自己对傅里叶的理解步骤（自我陶醉），请从标题实现介绍看起或请直接拉到最后看源码，或点击直接下载源码，有问题的小伙伴请直接留言或发邮件@1354647573@qq.com联系我。
时域主要用来直观的观察变换，而频域则是描述我们观察变换的一个趋势，我这里暂时把频域作为一个矢量来看，类似于用速度描述一个人的状态，而路程和时间的关系则体现速度的快慢，但是在一个很微小的时间，比如低于人体反应的时间0.02s，我们是感觉不到变化的，所以同样的，我们就很难描述这个时间里面发生的事。但是频域这个就能时刻展现这个变换的趋势，通过计算我们就可以很直观的看到整个变换了。这两个概念我一直不知道要怎么做一个通俗易懂深入人心的解释，用如下一个图简单举个例子。我们直观看到的是小球在t1时刻的位置和t2时刻的位置，但是我们看不到两个时刻的速度是多少，只知道速度的快慢，也就是小球在下一刻移动的距离远近，而这个速度就是我们在数学上严谨的一个描述性的矢量，就是这里说的频域。
image
说实话，上面的那段话就是为了绕晕一些新手的，毕竟本文的目的不在于介绍傅里叶的一些概念，而是根据概念来实现。
基础公式
下面摆出两个公式，为什么我会手写一次呢，大概是我觉得学习还是得自己动手吧，毕竟没有自己写过，或许提及某一个重要概念的时候，甚至都不知道它的具体内容，举个例子，比如你们看一下是否还记得几个排序算法的各种复杂度呢。
image
image
看到这两个公式，一般来说有部分人就头大。先不说连续傅里叶，就离散傅里叶变换公式，不妨自己亲手写一下，然后慢慢展开这两个离散的公式，看一下能不能发现什么，如果什么都没发现，那一定是偷懒了。
在说我发现了什么之前，先来说一个数学家，欧拉。非要说为什么要说他的话，就得从他的发现说起，欧拉公式
image
一个能把指数转化为三角函数的公式，方便了我们的计算方式。然后我们试一下把欧拉公式带入上述写好的展开式里，就会有一些比较清晰的思路了。三角函数转化，合并一些项...
到这里，基础的一些东西就算给完了，下面就给一些具体的做法。
具体实现
离散傅里叶变换，上面的公式，也叫DFT，很容易用编程语言实现，计算时间，和我们的输入数量成正比，简单的说一个问题，按照现在1920 1080个像素点，约210的6次方,DFT呢，每个点计算210的6次方，也就是总共需要410的12次方才可以计算完一张图片，我们简单算一下时间，假设cpu是4GHz，计算的时候假设是单线程满cpu负载，约等于1s处理4*10的9次方，我们做一个简单的除法，一张图片需要4000s计算出其频域，也就是比1小时还大，这个速度你们能忍吗，好吧，我知道了。这还是最直接的计算，涉及到的复数和实数的计算，加法和乘法的计算，如果是ARM架构的CPU我们做一下cpu指令集的优化计算，但是最后的速度也是很难让人忍受的。（假设每个像素点包含透明像素，则每个像素点4 byte。）
下面是一个java的dft的实现
/**
*这个是一维的DFT的简单实现，就是套公式
*arr 输入的一维的数值，因为现在的应用领域，所以输入一般都是整数
*flag 值为-1，考虑到逆变换也可以用同一个公式，所以就不多写一个函数了
*/
public static List<Complex> dft(int[] arr,int flag){
lengthC = arr.length;
if(lengthC == 0 ) {
System.err.println("---wrong length of array");
return null;
}
List<Complex> lc =new ArrayList<Complex>();
double vi,ri;
for(int j=0;j<lengthN;j++){
Complex cm = new Complex();
for(int i=0;i<lengthN;i++){
vi = Math.sin((flag*2.0*pi * i*j) / lengthN);
ri = Math.cos((flag*2.0*pi * i*j) / lengthN);
cm.virtual = cm.virtual + vi * arr[i];
cm.real = cm.real + ri * arr[i];
}
lc.add(cm);
cm = null;
}
return lc;
}
上面的代码有两个for循环，也就是O(n2)的时间复杂度。虽然计算机的处理能力在不断增强，但是这个时间复杂度不能忍。
所以在数学层面就有先驱者做了一些优化，当输入个数是合数时可以采用快速傅里叶变换，当n为大素数的时候，可以采用线性调频Z变换(Chirp- Z 变换)算法。具体概念请参照wiki百科。
我们先从简单的学起吧，毕竟我只是一个认真不起来的学渣。
当输入的个数为2的幂次方时，这个就给人一种熟悉的感觉，是不是好像很多东西都会要求是2的幂次方，在这个时候，公式可以做一些简单的化简，之前有展开的公式的同学可以把奇数项和偶数项分开，然后奇数项的特征呢，2n +1,偶数项的特征呢，2n, 这个时候奇数和偶数项的个数分别都是原来的一半，这时候，在思考的同学大概想到了，每一项都是从n=0开始到N/2-1，然后根据公式的规律，变成了两个小的傅里叶变换。一直重复上述操作，最终就可以划分为 N/2 个 只有两个值的傅里叶变换。这样计算的方式的复杂度就是2*2 *N/2 * logN， logN是分出每两项的次数。
图片来自wiki
见到那说一下这幅图应该怎么看吧，如下图，这样子看起来就很直观了，实现起来的时间也就是对编程语言熟悉程度的问题了。
简单分析
授人与渔呢，还是得授人如何渔。
我们在这里采用DIT的实现，所以我们一般会对数据进行一些处理，这里由于本人计算机基础拙劣，没有采用一个好的数据结构，就使用了最基础的数据交换方式
void bitReverse(Complexf *cm,u_int N){
u_int temp,idx;
Complexf *res = (Complexf *)malloc(sizeof(Complexf) * N);
u_int L = log2(N);
Complexf tpcm = Complexf();
for(u_int i=0;i<N;i++)
{
temp = 0;
idx= i;
for(int j = 0;j<L;j++)
{
temp = temp << 1;
temp += idx%2;
idx = idx >> 1;
}
res[i].setData(cm[temp]);
}
for(int i=0;i<N;i++){
cm[i].setData(res[i]);
}
free(res);
};
进行上面处理之后再进行FFT
/**
*这里只是我用C++实现的FFT, Complexf是自己写的一个复数实现。
*wnk是e^(2*PI*i*n/N) n=0到N-1的计算出来的一张表
*cm是做了预处理之后传入的数值
*N传入的个数
*/
void fft_main_base2(Complexf *cm,int N,Complexf *wnk){
int m = log2(N);
//辅助变量
Complexf temp = Complexf();
Complexf temp2 = Complexf();
/*
*ni 当前深度
*cn当前深度的数量，也就是当前有互不干扰的计算，比如上图深度最大时计算的x[0],x[4]，计算x[2],x[6]互不影响，
*pn暂时没用
*nd 辅助变量
*/
u_int nd,cn,pn,ni,wnid1,wnid2;
//级数，logN,cn means current is in which level.
for(int i=0;i<m;i++){
cn = 2 << i;
ni = cn / 2;
pn = N / ni;
//每一级对应点
for(int j=0;j<ni;j++){
for(int k=j;k<N;k+=cn){
nd = ni + k;
temp2.setData(cm[k]);
temp.setData(cm[nd]);
wnid1=((k%cn)* (N/cn))%N;
temp *= wnk[wnid1];
cm[k] += temp;
temp.setData( cm[nd] );
wnid2 = ((nd%cn) *(N/cn)) %N;
temp *= wnk[wnid2];
cm[nd].setData( temp2);
cm[nd] += temp;
}
}
}//理论总时间 = M * pn * 2 * N / cn = N *logN
/*if(wnk!=nullptr)
delete(wnk);*/
}
上面的两段就是输入个数为2的幂次方的时候FFT的基本实现思路了。
合数下的FFT
这个数为2的幂次方的数确有点特殊了，适合用来做FFT的入门。在数学上，我们追求的是捕获事物的本质，但是我们可以从特殊下手，找寻规律，再去捕获本质。
当N为普通合数的时候，这时我们之前的方法是不是行不通了呢。简单一想的确是这样，转念一下，假设 2 * 2 * 2 变成 2 * 2 * 3呢，规律始终是有的，就看我们能不能发现它。
就这样，假设N是2 * 2 * 3， 我们把之前的展开式，先分成2 堆，在分为4堆，这4堆里面每堆都有三个数，那么我们做一次个数为3的DFT不就行了吗？这样考虑起来是不是感觉合数的FFT就是把2基底FFT(上面介绍的FFT)的2换成了其他素数，其实也不是那么难，所谓的困难只是你那颗不敢尝试的心。
ok，这样子一换，之前的运算参考图就不管用了，所以又是一波新的图。但是原理，只能说没变化，下面是我的一个思路和代码。
首先我们直接拿到输入个数的所有因子。这里呢，暂时设置素数因子小于19，获取因子这里可以自由发挥。(这些个因子的作用主要是为了分解合数嘛)
/*
get number of prime the N can divide
*/
int CN_prime[] = {2,3,5,7,11,13,17,19};
int CN_num = 8;
#define ffor(a) for(int i=0;i<a;i++)
u_int* getPrime(int N){
int j = 0;
int k = 0;
u_int tepN = N;
u_int *temp = (u_int *) malloc(sizeof(temp) * CN_num);
ffor(CN_num){
temp[i] = 0;
}
while(j<CN_num){
if(tepN%CN_prime[j] == 0) {
tepN = tepN /CN_prime[j];
temp[j]++;
k++;
}
else{
j++;
}
}
u_int* res = new u_int[k+1];
res[0] = k;
k = 1;
ffor(CN_num){
while(temp[i]>0) {
res[k++] = CN_prime[i];
--temp[i];
}
}
//delete temp
free(temp);
//remember to delete
return res;
}
接着呢，我们通过拿到的素数，开始进行翻转数据了，主要是为了处理好数据可以直接使用。这里没什么特别的，和上面的2的幂次方的处理也是差不多，2的幂次方的数所有因子都是2，这里的因子不是2而已，简单思考一下是很容易实现的。
/**
*合数bit reverse,O(N)
*pN是素数数组，pN[0]表示数组长度，N复数长度
*/
void CN_reverse(Complexf *cm,u_int *pN,u_int N){
u_int *temp = new u_int[N];
Complexf *tmp =new Complexf[N];
for(u_int i=0;i<N;i++){
temp[i] = 0;
tmp[i].setData(cm[i]);
}
u_int idx,md,t;
idx = 1;
t = 1;
for(u_int i= 1;i<pN[0]+1;i++){
idx *= pN[i];
md = N/idx;
for(u_int k=t;k<idx;k++){
temp[k] = temp[k-t]+md;
}
t = idx;
}
for(u_int i=0;i<N;i++){
cm[i].setData(tmp[temp[i]]);
}
std::cout<<" reverse end\n";
//删除临时复数数组
delete[] tmp;
//删除临时下标数组
delete[] temp;
}
数据处理之后，就可以开始我们的合数傅里叶变换了
/**
*这里用来快速计算的，主要是wnk的叠加，就类似2的幂次方的那种快速计算。
*这个算法是缩减运算dft时间的一个算法，也就是合数fft的理解精髓。
*/
void fast_compute(Complexf *cm,Complexf* res,Complexf* wnk,u_int n,u_int *nk){
Complexf temp = Complexf();
Complexf temp1 = Complexf();
for(int i=0;i<n;i++){
temp1.setReal(0);
temp1.setIm(0);
for(int j=0;j<n;j++)
{
temp.setData(cm[j]);
temp *= wnk[(j *(nk[(i)%n]))%nk[n]];
temp1 += temp;
}
res[i].setData(temp1);
}
}
/*
*prim是所有的因子，prim[0]为因子的数量
*wnk这里不做赘述
*/
void fft_main_CN(Complexf *cm,Complexf *wnk,int N,u_int *prim){
int nL = prim[0];
//长度为1，直接用dft计算。
if(nL==1)
{dft(cm,wnk,N);return;}
int M = nL+1;
u_int *nk = new u_int[N];
Complexf *res = new Complexf[N];
Complexf *tempf = new Complexf[N];
Complexf tempmm = Complexf();
u_int cn,md,ni,mp,nidx;
cn = 1;
ni = 1;
//当前的因子
for(int i=1;i<M;i++){
cn *= prim[i];
md = N /cn;
mp = prim[i];
//快速计算所需要
nk[mp] = N;
for(int m=0;m<N;m+=cn){
for(int j=0;j<ni;j++){
for(int k=0;k<mp;k++){
nidx = k*ni+j+m;
tempf[k].setData(cm[nidx]);
nk[k]=((nidx%cn)* md)%N;
}
//首先要做一个dft，然后才去做快速计算
if(i==1){ dft(tempf,res,wnk,mp,nk); }
else{ fast_compute(tempf,res,wnk,mp,nk); }
//把计算好的设置到原输入去
for(int k=0;k<mp;k++){
nidx = k*ni+j+m;
cm[nidx].setData(res[k]);
}
}
}
ni = cn;
}
delete [] res;
delete [] tempf;
delete [] nk;
}
以上就是这次说的一个快速傅里叶计算的方式，实现起来就和描述一样，不过要理解这个描述可能还是需要一些时间的。有疑问的同学可以评论，也可以发邮件到1354647573@qq.com，很开心你能看完我的拙劣理解与实现。
