目录
傅里叶变换原理
讲解傅里叶变换的文章有许多，而且洞见深刻，我作为一个学习者就不再从根源开始探讨了，这里给出几个大佬的文章链接。
好了，下面默认你对DFT思想有一定了解，接下来开始推导公式吧。
首先我们要明确傅里叶变换的目的，那就是某些信号在时域很难看出特征并区分，但是转换到频率域，就使得某些特征被明显区分。
傅里叶变换的本质是内积，傅里叶变换就是将满足一定条件的某个函数表示成三角函数（正弦/余弦）或其积分的线性组合，将一个连续的信号（不方便处理）转换成一个个小信号的叠加（好处理）。
傅里叶原理表明：对于任何连续测量的数字信号，都可以用不同频率的正弦波信号的无限叠加来表示。
理论推导
（DFT）序列x(n)(n=0,1,⋯,N−1)的离散傅里叶变换定义为：
（IDFT）序列X(k)的离散傅里叶反变换定义为：
为了把难搞的指数函数替换掉，请出欧拉公式：
将欧拉公式带入DFT公式，式中，k=0,1,⋯,N−1，这就是接下来我们做正变换的公式了：
同样地，将欧拉公式带入IDFT公式，式中，n=0,1,⋯,N−1：
接下来很重要的一步，将正变换的公式带进反变换公式，x(n)可以进一步展开为四项：
现在我们来分析一下，①和④中不包含虚数单位（④中两个j等价为-1），也就是说，它俩个为实数部分，很好，很符合时间域，毕竟x(n)表示时间域的离散序列数组嘛。
但是②和③呢，烦！都有虚数单位，怎么办？
要不直接忽视这两项吧，时间域不应该存在虚数，那这两项势必会消掉。
由三角函数我们知道，sin(-x)=-sinx;cos(-x)=cosx，那么②和③中的展开式的每一项都正好互相抵消！
舒服了！
现在我们把上式整理一下，就是接下来做反变换的公式了：
代码实现
正变换代码实现
那我们从hello world开始吧。
先创建一个文本：touch DFT.c
首先我们要做DFT，那必须要有用于变换的数据，为了适用于大多数看这篇文章的csdner，使用雷克子波做为演示，同时为了直观，把结果写到磁盘。
（为了下一盘大棋，使用1024数据长度的雷克子波。如果未来我有时间的话，会继续讲解FFT、二维FFT、褶积等等，使用2的幂长度的数据更方便理解）
int it,fm;
float t1,t0,N,dt;
FILE *fp;
N=1024;dt=0.0005;t0=0.1;fm=25;
float *ricker=(float *)malloc(sizeof(float)*N);
for(it=0;it<N;it++){
t1=it*dt-t0;
ricker[it]=(1.0-2.0*PI*PI*fm*fm*t1*t1)*exp(-(PI*PI*fm*fm*t1*t1));
}
fp=fopen("ricker.bin","wb");
fwrite(ricker,sizeof(float),N,fp);
fclose(fp);
现在数据有了，那就开始做正变换吧！
我们再把刚才推出的公式摘出来，照着敲代码吧，冲呀！
这里我们虽然使用的一维的ricker，但是为了匹配以后的二维数据，我们使用二维数据的DFT写法，做三层for循环，那雷克子波就相当于Nx==1。
傅里叶变换之后会产生实数部分和虚数部分两部分数据，同时，我们还可以使用两部分数据产生振幅数据amp，最后我们把三块数据写到磁盘。
最后，为了看看DFT用了多久时间，我们include<time.h>，调用时间函数。
void DFT(float * input,int Nx,int Nz,float * real,float * imag)
{
float *amp=(float *)malloc(sizeof(float)*Nx*Nz);
float R,I;
float re,im;
int ix,in,ik;
for(ix=0;ix<Nx;ix++)
{
for(ik=0;ik<Nz;ik++)
{
real[ix*Nz+ik]=0;
imag[ix*Nz+ik]=0;
R=0;
I=0;
for(in=0;in<Nz;in++)
{
R=cos(((-2)*PI*ik*in)/Nz);
I=sin(((-2)*PI*ik*in)/Nz);
real[ix*Nz+ik] += input[ix*Nz+in]*R;
imag[ix*Nz+ik] += input[ix*Nz+in]*I;
}
re=real[ix*Nz+ik];
im=imag[ix*Nz+ik];
amp[ix*Nz+ik]=sqrt(pow(re,2)+pow(im,2));
}
}
FILE *fp=fopen("ricker_dft_real.bin","wb");
FILE *fm=fopen("ricker_dft_imag.bin","wb");
FILE *fr=fopen("ricker_dft_amplitude.bin","wb");
fwrite(real,sizeof(float),Nx*Nz,fp);
fwrite(imag,sizeof(float),Nx*Nz,fm);
fwrite(amp,sizeof(float),Nx*Nz,fr);
fclose(fp);
fclose(fm);
fclose(fr);
free(amp);
}现在看看实部、虚部、振幅谱数据吧！
这就是实部、虚部、振幅数据，很直观地，我们可以发现最突出的特点，他们在频率域是对称的。
这启发我们为了减小运算量，是否可以只计算一半呢，这部分内容会在以后讲解。
同时可以发现，振幅谱有一个波峰，那实际上就是我们给的主频25HZ。恰好验证了单一的频率在频率域显得十分空旷，数据中间有大量零值。（我们可以试想一下，如果时间域一个复杂的信号在频率域会被分解成什么样子呢，下一章节会介绍，歪七扭八、平平无奇的时域数据，在频域可谓是散如满天星！）
正变换讲完了，那反变换是如何开展的呢，我们继续！
反变换代码实现
通过前面的学习，我们有了正变换产生的实部和虚数数据，这两块数据就是用来做反变换的输入数据，等再一次反变换回到时间域，我们拿反变换的信号与最开始做的雷克子波做个残差，看看反变换的结果是否合理。
同样的，我们把反变换的公式摘下来，照着敲就好。
公式中
表示之前算出来的实部数据，对应的
表示之前正变换算出来的虚部数据
这里要务必注意！！！！
因为在正变换的时候，虚部数据的j我们是无法在表达式中体现的，所以j就被无视了，现在做反变换的时候，反变换表达式x(n)中的虚数单位可不能视而不见，因为后跟的表达式 j*sin(2*PI*k*n/N) 中又出现了一个虚数单位，j*j=-1；
也就是说，在实际计算中，反变换的表达式x(n)中间不再是加号了，而是减号。
void IDFT(float *real_part,float *imag_part,int Nx,int Nz,float *idft_part)
{
float *real=(float *)malloc(sizeof(float)*Nx*Nz);
float *imag=(float *)malloc(sizeof(float)*Nx*Nz);
float R,I;
int ix,in,ik;
for(ix=0;ix<Nx;ix++)
{
for(in=0;in<Nz;in++)
{
real[ix*Nz+in]=0;
imag[ix*Nz+in]=0;
R=0;
I=0;
for(ik=0;ik<Nz;ik++)
{
R=cos((2*PI*ik*in)/Nz);
I=sin((2*PI*ik*in)/Nz);
real[ix*Nz+in] += real_part[ix*Nz+ik]*R;
imag[ix*Nz+in] += imag_part[ix*Nz+ik]*I;
}
idft_part[ix*Nz+in]=((float)(real[ix*Nz+in]/Nz)-(float)(imag[ix*Nz+in]/Nz));//a-b !!!!
}
}
FILE *fp=fopen("ricker_idft.bin","wb");
fwrite(idft_part,sizeof(float),Nx*Nz,fp);
fclose(fp);
free(real);
free(imag);
}
残差
现在我们有两份时间域数据了，一份是根据雷克子波公式写的，一份是反变换出的，目测两个波形一致，为了更严谨一点，做个残差。
残差做了不算完，为了定量分析两者是否匹配，拿残差数据求取其L2范数，作为评判。
L2=0;
for(it=0;it<N;it++){
residual[it]=ricker[it]-Idft_part[it];
L2+=residual[it]*residual[it];
}
fprintf(stderr,"Residual between DFT with IDFT is %f\n",sqrt(L2));
fp=fopen("ricker_residual.bin","wb");
fwrite(residual,sizeof(float),N,fp);
fclose(fp);残差如下图，实际上，这一点也不奇怪，因为受限于数据类型的精度，两部分数据一定会有差异的，如果将数据类型由float改为double误差就会稍微小点，但是一定不能消除，毕竟计算机底层计算原理是离散的。
程序执行结束打印输出结果：
DFT time taken: 0.000000
IDFT time taken: 0.000000
Residual between DFT with IDFT is 0.000002
时间耗费分析：只能说数据量太小，咻的一下就算完了，之后搞的大的数据来测试。
残差分析：DFT和IDFT做的残差求取的L2范数为0.000002（2e-6），那就是说在十的负五次方精度下，两者一致。验证了DFT和IDFT算法的正确性。
完美！
全部源代码
编译方法：gcc DFT.c -o DFT -lm
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define PI 3.1415926
void DFT(float * input,int Nx,int Nz,float * real,float * imag);
void IDFT(float *real_part,float *imag_part,int Nx,int Nz,float *idft_part);
int main()
{
int it,fm;
float t1,t0,N,dt,L2;
clock_t start,end;
FILE *fp;
N=1024;dt=0.0005;t0=0.1;fm=25;
float *ricker=(float *)malloc(sizeof(float)*N);
for(it=0;it<N;it++){
t1=it*dt-t0;
ricker[it]=(1.0-2.0*PI*PI*fm*fm*t1*t1)*exp(-(PI*PI*fm*fm*t1*t1));
}
fp=fopen("ricker.bin","wb");
fwrite(ricker,sizeof(float),N,fp);
fclose(fp);
float *Real_part=(float *)malloc(sizeof(float)*N);
float *Imag_part=(float *)malloc(sizeof(float)*N);
float *Idft_part=(float *)malloc(sizeof(float)*N);
float *residual =(float *)malloc(sizeof(float)*N);
start=clock();
DFT(ricker,1,N,Real_part,Imag_part);
end=clock();
fprintf(stderr,"DFT time taken: %lf\n",(double)((end-start)/CLOCKS_PER_SEC));
start=clock();
IDFT(Real_part,Imag_part,1,N,Idft_part);
end=clock();
fprintf(stderr,"IDFT time taken: %lf\n",(double)((end-start)/CLOCKS_PER_SEC));
L2=0;
for(it=0;it<N;it++){
residual[it]=ricker[it]-Idft_part[it];
L2+=residual[it]*residual[it];
}
fprintf(stderr,"Residual between DFT with IDFT is %f\n",sqrt(L2));
fp=fopen("ricker_residual.bin","wb");
fwrite(residual,sizeof(float),N,fp);
fclose(fp);
free(ricker);
free(Real_part);
free(Imag_part);
free(Idft_part);
free(residual);
return 0;
}
void DFT(float * input,int Nx,int Nz,float * real,float * imag)
{
float *amp=(float *)malloc(sizeof(float)*Nx*Nz);
float R,I;
float re,im;
int ix,in,ik;
for(ix=0;ix<Nx;ix++)
{
for(ik=0;ik<Nz;ik++)
{
real[ix*Nz+ik]=0;
imag[ix*Nz+ik]=0;
R=0;
I=0;
for(in=0;in<Nz;in++)
{
R=cos(((-2)*PI*ik*in)/Nz);
I=sin(((-2)*PI*ik*in)/Nz);
real[ix*Nz+ik] += input[ix*Nz+in]*R;
imag[ix*Nz+ik] += input[ix*Nz+in]*I;
}
re=real[ix*Nz+ik];
im=imag[ix*Nz+ik];
amp[ix*Nz+ik]=sqrt(pow(re,2)+pow(im,2));
}
}
FILE *fp=fopen("ricker_dft_real.bin","wb");
FILE *fm=fopen("ricker_dft_imag.bin","wb");
FILE *fr=fopen("ricker_dft_amplitude.bin","wb");
fwrite(real,sizeof(float),Nx*Nz,fp);
fwrite(imag,sizeof(float),Nx*Nz,fm);
fwrite(amp,sizeof(float),Nx*Nz,fr);
fclose(fp);
fclose(fm);
fclose(fr);
free(amp);
}
void IDFT(float *real_part,float *imag_part,int Nx,int Nz,float *idft_part)
{
float *real=(float *)malloc(sizeof(float)*Nx*Nz);
float *imag=(float *)malloc(sizeof(float)*Nx*Nz);
float R,I;
int ix,in,ik;
for(ix=0;ix<Nx;ix++)
{
for(in=0;in<Nz;in++)
{
real[ix*Nz+in]=0;
imag[ix*Nz+in]=0;
R=0;
I=0;
for(ik=0;ik<Nz;ik++)
{
R=cos((2*PI*ik*in)/Nz);
I=sin((2*PI*ik*in)/Nz);
real[ix*Nz+in] += real_part[ix*Nz+ik]*R;
imag[ix*Nz+in] += imag_part[ix*Nz+ik]*I;
}
idft_part[ix*Nz+in]=((float)(real[ix*Nz+in]/Nz)-(float)(imag[ix*Nz+in]/Nz));//a-b !!!!
}
}
FILE *fp=fopen("ricker_idft.bin","wb");
fwrite(idft_part,sizeof(float),Nx*Nz,fp);
fclose(fp);
free(real);
free(imag);
}一维离散傅里叶变换就到这里，但是呢，这个代码写的太死了，只能用于测试算法正确性，接下来为了让它灵活，为了可以计算任意给定数据的DFT和IDFT，下一篇文章我会用argc、argv、atoi等工具对它稍做升级，实现任意数据的离散傅里叶变换。
知识分享是我的初衷，在我自己学习过程中，鲜有案例讲解，踩了许多坑，如今窥探到了一点傅里叶的深邃思想，所以才有这篇文章。
如果你觉得这篇文章解决了你某个疑惑，不妨点个赞，让更多人看到，谢谢！
