[音频处理]专科狗眼中的傅里叶变换人工智能工农村贴膜小哥的博客-

写在前面

不是科研狗，基础理论薄弱，写的比较匆忙，有理解有误的地方还请理解和指正。
网上大佬们写的傅里叶公式推导，证明已经很多了（瑟瑟发抖），我这里主要是讲傅里叶的应用，不涉及公式证明，而是直接拿起公式使用。
由于自己获取知识也是看大佬们博文理解学习得来的，所以图片中多少有一些是别人的图，不过我附上了别人的链接。
看完这篇你能收获到：
1 傅里叶变换的原理
2 傅里叶变换在音频的应用
3 离散傅里叶变换处理音频的C语言代码及讲解

背景

最近接触音视频处理比较多，就遇到了采集的音频数据有噪音的情况。可不可以用傅里叶变换去除音频噪音呢？在此之前，我对傅里叶变换的概念只有一句话，“时域转频域的一个玩意儿”。我将音频时域转为频域，再去除噪音的频段，不就达到去噪的效果吗？ 嗯，试试就逝世。
嗯，经过几天的折腾，效果（上图为原始音频信号，下图为在频域去掉高频正交基后的音频信号）
结论：
1，有点效果，但噪音只是变小了，没有去除。使用的是DFT，还有待改进。
2，琢磨了一段时间，还有诸多地方没有理解到，傅里叶变换是需要长期研究和推敲的学问。把自己理解到的知识点做下笔记，以后要走图像处理什么的路子，再来细细研究吧。

一 什么是傅里叶变换（分析）？

可以参考这篇
在回答这个问题的时候，我们需要先知道什么是变换?为什么要变换？

（1）什么是变换？

举个大家常用的例子,两个向量，a，b在坐标系中表示，可以变换为一组坐标。

（2）为什么变换？

比如想求c点坐标，在图中是不好处理的问题，但是变换后,c=a+b=(2+1,1+2)=(3,3) ，在数上就很好处理了。总之，一个问题不好处理的时候，换个思路（变换），就能迎刃而解。嗯嗯，有点像线性代数的味道。
傅里叶变换的思路就是时域和频域的相互变换。

二 时域和频域

时域和频域一一对应，通过频域可以变换为时域

傅里叶变换类型

傅里叶变换有许多类型，通常不加前缀说的是连续傅里叶变换，离散傅里叶变换常用在计算机处理上，在此之上进行算法改进，降低时间复杂度，叫快速傅里叶变换。

对于时域是周期性连续的函数进行变换叫傅里叶级数
对于时域是非周期连续的函数进行变换叫傅里叶变换

三傅里叶级数（Fourier Series）

所以啊，傅里叶就提出了这句话：任何周期函数都由多个不同的正弦波（sin和cos）叠加构成（满足狄里赫利条件情况下）
如图，频域由若干个正弦波构成，时域的矩形波可以由多个正弦波叠加构成（当正弦波个数趋向无穷大，逼近为矩形波，有点极限的赶脚）

周期函数其实说的就是时域信号，多个不同的正弦波说的是频域的信号。

当然，如上图所示，一个正弦波的确定，需要知道他的相位，振幅和角频率，所以时域转频域得到的是一系列的初相，振幅和角频率。
傅里叶说的这句话写成公式就是

又或者可以这样写

这个公式该这么理解呢？这就需要知道正弦波的正交性。

（1）正弦波的正交性

正弦波自己给自己内积为1，自己和别人内积为0，也就是sin（nwt），cos（nwt）都是标准正交基，再加上常量1， 任何的周期函数都可以用着三个表示 。所以呢，只需要将时域的f(t)分别和各种标准正交基做内积，等于0代表没有该正弦波，否则存在该分量。

什么是标准正交基

线性代数：在空间中找到一组向量，他们的内积为0，自身的模为1时，称这组向量为标准正交基，空间中的任意向量都可以由这些标准正交基构成。

什么是内积？

嗯，你需要学习到晚上2点[手动滑稽]

三傅里叶变换公式

（1） 欧拉公式

但上图的f(t)没有和各种正弦波做内积啊，而是和e^-iwt 做内积。其实这是欧拉公式将正弦波与指数之间进行了变换，本质上上图公式就是和各种正弦波做内积。现在就说说欧拉公式。

数系

回顾下初中知识
自然数：1，2，3，4。。。
整数：-1，1，-2，2.。。。
有理数：-1，2，2/3，4/5.。。
无理数：根号2等
这些一起就构成了实数，如图 实数轴。
但遇到x^2=-1这样的解，用实数就无法解决了，引入虚数。虚数+实数=复数，如图为复平面，任意点可以用a+bi 来表示

在复平面上画圆，由三角公式，可以求得a点为cosx+i *sinx。
复平面上乘法的意义：
乘-1，逆时针旋转180度。乘-i ，逆时针旋转90度，复数乘法如图

现在就可以 拿出欧拉公式了

公式推导就是使用的泰勒公式，e，sin，cos都是十大必背泰勒展开式，两边展开是相等的。
到这里就打住吧，总之我们到这里能知道傅里叶变换公式e^iwt其实就是与各个sin，cos做内积。具体欧拉公式讲解可以看这篇
欧拉公式讲解
欧拉恒等式
也就是角度到达180度时候，又回到实数轴

离散傅里叶变换（DFT）

离散傅里叶变换可以参考这篇
离散傅里叶变换讲解
离散和连续的区别在于，连续是积分而离散是求和

原理，如下图，a,b图是待检测信号，c，d是3个周期的正弦信号，很显然a图含有正弦波，e=ac，将e图的各点相加，很显然值是正的，这就说明a图含频率为3的正弦波，f=bd，显然将f图中各点相加结果约等于0了，说明b图不含有周期为3的正弦波

在计算机中用这个公式更好处理一点

n和N是在一个正弦周期内采样N个点，采样间隔为2piN,n用来步进，一次步进2piN,最后进行累加求和，就得出了X(k)

最后 离散傅里叶变换完整代码

1，从文件读取8000个音频数据，由于现实中的音频没有虚部，所以只设置实部。
2，离散傅里叶变换关键处
temp的re就是对应上图公式的cos，同理im就是对应上图的sin，每个X[k]进行累加求和

 for (int k = 0; k < N; k++) {   X[k].re = 0;   X[k].im = 0; for (int n = 0; n < N; n++) {    temp.re = (float)cos(2 * pi*k*n / N);    temp.im = -(float)sin(2 * pi*k*n / N);    X[k] = complexadd(X[k], complexMult(x[n], temp)); } } 

3，离散傅里叶逆变换就是X（k）乘上e^(j2tkn/N),也就是实部虚部都为正

temp.re = (float)cos(2 * pi*k*n / N);    temp.im = (float)sin(2 * pi*k*n / N);    x[k] = complexadd(x[k], complexMult(X[n], temp)); 

4，下图公式求幅度

在代码中表示

printf("%d ", int(2.0/N*sqrt(X[k].re*X[k].re + X[k].im*X[k].im))); //打印的为幅度 

完整代码

 #include <stdio.h> #include <math.h> #define pi 3.1415926 struct complex { float re; float im; };  complex complexadd(complex a, complex b) { //复数加  complex rt;  rt.re = a.re + b.re;  rt.im = a.im + b.im; return rt; }  complex complexMult(complex a, complex b) { //复数乘  complex rt;  rt.re = a.re*b.re - a.im*b.im;  rt.im = a.im*b.re + a.re*b.im; return rt; } complex complexSet(complex *a, short *b, int N)//复数设置 { for (int i = 0; i < N; i++) {   a[i].re = b[i];   a[i].im = 0; } } //离散傅里叶变换 //X[]标识变换后频域，x[]为时域采样信号 void dft(complex X[], complex x[], int N) {   complex temp; int k, n; for (int k = 0; k < N; k++) {   X[k].re = 0;   X[k].im = 0; for (int n = 0; n < N; n++) {    temp.re = (float)cos(2 * pi*k*n / N);    temp.im = -(float)sin(2 * pi*k*n / N);    X[k] = complexadd(X[k], complexMult(x[n], temp)); } printf("%d ", int(2.0/N*sqrt(X[k].re*X[k].re + X[k].im*X[k].im))); //打印的为幅度 if (k >= 6000)//去除高频信号 {    X[k].re = 0.0;    X[k].im = 0.0; } } } //离散傅里叶逆变换 //X[]标识变换后频域，x[]为时域采样信号 void idft(complex X[], complex x[], int N) {  complex temp; //int k, n; for (int k = 0; k < N; k++) {   x[k].re = 0;   x[k].im = 0; for (int n = 0; n < N; n++) {    temp.re = (float)cos(2 * pi*k*n / N);    temp.im = (float)sin(2 * pi*k*n / N);    x[k] = complexadd(x[k], complexMult(X[n], temp)); }   x[k].re /= N;   x[k].im /= N; } } #define N 8000//采样率为8000 int main() {  complex samples[N], X[N], x[N]; //原始数据，变换后的频域数据，逆变换后的时域数据  FILE *fp;  FILE *fp2; short buf[N]; int re=0;  fp=fopen("./ttt.pcm", "rb");  fp2 = fopen("./trans.pcm", "wb"); if (!fp ||!fp2) { printf("I could not open file.n"); return 1; } else { while (fread(buf, sizeof(short)*N,1 , fp) > 0)//末尾数据忽略 { complexSet(samples, buf, N); dft(X, samples, N); idft(X, x, N); for (int i = 0; i < N; i++) {     buf[i] = x[i].re; } fwrite(buf, sizeof(short)*N,1 , fp2); } } fclose(fp); fclose(fp2); return 0; }