tucao

先要吐槽一下简书，居然用不了MathJax，导致所有的公式要写成Latex，真是反人类
||Φ|(|T|Д|T|)|Φ|||

void 吐槽()
{
最近刚看完FFT，顿时心中真是成百上千的mmp。倒不是方法本身的问题，而是现在的书籍和资料很多真是对萌新极其不友好，尤其是国内一些所谓的“入门到入坑”系列，总是把一些简单的问题说复杂，而且描述同一个公式还版本多种多样，并且把一些应该放到后面的知识排到了前面，要么让你望而生畏，要么间接up你查阅资料的能力。
好的，吐槽完毕。
return
}

这篇简书只是关于DFT笔记的汇总，大部分是个人主观的理解，对DFT理解不深，文中一些不严谨的地方欢迎指出纠正。

开篇之前要先感谢一位师兄，将BASIC版本的Bit Reversal Sorting 翻译成了C版本，本文FFT的程序才能完整实现。
你部长还是你部长 <)。(>

目录

• DFT基本理论
• DFT四种算法
• 使用DFT去噪声与解调

DFT基本理论

DFT全名Discrete Fourier Transform (离散傅里叶变换) ，是傅氏变换家族中用于处理在横纵轴上数值离散且具有周期性的信号。相比于另外三个，这种变换方法处理有限并且周期延拓的信号，所以适用于计算机处理。

FT主要的目的是把时域信号中的各频率成分的振幅和相位分别提取出来，从而将时域信号转换位频域信号中的幅频谱和相位谱上。
对于一维的时域信号，要得到频谱只需进行一次DFT；而对于二维（图片）以及二维以上的信号，需要进行多次的不同方向上的DFT，本文只涉及一维信号的DFT，对于变换图片的频谱请参考相关资料并类推即可。

DFT四种算法

对时域信号进行DFT的方法有很多，简wo单zhi的hui有以下几种：

一、联立方程组求解 Σ( ° △ °|||)︴

根据傅里叶分解的原理，任何信号都可以多组正交的正弦曲线叠加而成
每组正弦曲线的频率ω不同，相位Ф不同，振幅A不同
根据香农采样定理，采样率应为信号最高频率的两倍以上，取最低要求情况，若刚好位两倍的情况，ω的范围是0（直流）到0.5f（采样频率的一半）。
列出方程组求解得到多组的解A和Ф，求得幅频数组Mag[]和相位数组Phase[]。

线性方程组

数学专业的同学了解一下……
嗯，1024个点为1024个线性方程组，解吧:-)

二、相关性求解法（传统方法）

方法一是从数学的角度来理解DFT，可以作为理论上的用法，但是由于计算量过大，不能用于计算机处理。

若从求相关性的方法来计算，则可用于计算机处理DFT。
相关性求解的方法属于一种类卷积的运算，通过将原始信号与目标频率信号求相关，可以得到二者的相关度，再求代数和，即可得到该频率下的信号相关度。

求相关

这里改方法根据原始信号的种类不同又可以分为两种：

• 实信号处理模式（更快，但是公式混乱）
  对于N点的实数信号，根据以上的方法，可以分解为N/2点的频域数组Re[]和Im[]。所以求相关性时应该使用对应频率的正弦函数sin和相同频率的余弦函数cos，求得的相关值分别放入Im[]和Re[]。

• 复信号处理模式（慢一些，但是统一度好）
  对于N点的信号，将其当作复数处理，相当于有N点的实部和N点的虚部（对于纯实数信号，N点虚部全为0），再将其与旋转因子（有的版本称相位因子）W_N

进行相关度运算，得到实部Re和虚部Im。
值得注意的是，此时Re数组和Im数组长度均为N而不是N/2。

DFT&IDFT

（0 ~ N/2）部分的数值和“实数模式中得到的结果一摸一样，称之为正频率，而（N/2 ~ N）的数值为负频率。
负频率由于本身没有物理意义，且和正频率水平镜像对称，为冗余信息，故可以舍弃不运算加快速度（但是负频率在FFT中有着其存在意义）。

得到相关性数组Re[]和Im[]，是直角坐标系下的频谱，在运算上有作用，但是不直观。

极坐标形式的幅频图

根据高中学过的Asin(ωt+Ф)=Xsin（ωt）+Ycos（ωt）可知，Re[]和Im[]的信息可以等价保存为幅频数组Mag[]和相位数组Phase[]，将直角坐标系转化为极坐标系，使得数据更加直观（尤其是频谱部分）。

三、对称性求解法

对称性方法比较容易理解，但是计算量却十分的大（仅次于方程组求解法）

我们知道，频域的一个δ(n)冲激响应（一个点）相当于时域的一个响应频率的正弦曲线，根据时域与频域的对称性，时域的一个δ(n)也相当于频域的一条正弦曲线
（算是是自然界中最神奇的对称之一）
所以，可以从时域逐点求得对应的正弦曲线，在频域进行累加运算。
时域上每个δ冲击函数都有位移，相当于δ(n)与Δ(n+1)进行了一次卷积，所以频域的每次累加要先乘于一个对应频率的正弦曲线。（因为Δ(n)的傅里叶变换对是sin）。
说白了就是左边一个个数幅度，然后右边慢慢加

四、基2法（“饥饿法”！？）

还有一种叫基4法（这篇比较全面）：
http://www.doc88.com/p-7979920042436.html

基2法包括DIT时域抽取和DIF频域抽取，这里只简单介绍以下DIT。
基2法全名是 DIT-FFT-基2 算法。
此处FFT是Fast Fourier Transform即快速傅里叶算法，FFT有很多种，需要注意的是FFT不是新的傅里叶分解方法，它只是DFT的一种方法改进，就是公式还是原来的，只是优化了。
相比于传统的DFT，FFT大约能快个300倍左右，并且数据量越大，FFT越有优势。

• 传统的DFT为：

DFT

对于N点的时域信号，算一次DFT，需要进行 _N2 次复数乘法和累加运算，这就是传统DFT的局限性之一。
在100Mz的处理器上，计算 ₂10 个采样点，即便预先算好旋转因子制成一个数组表，估计也要十几秒，那么对于更大的点的速度将会是不可接受的（Time out啦兄dei）。

之所以这么慢，是因为传统DFT的计算有着大量的冗余。利用旋转因子W_N^nk本身的对称性、周期性、可约性以及特殊值，简化DFT的计算步骤得到FFT算法。

• 分而治之
  FFT算法利用的是“分而治之”的算法思想，对于DIT基2型的FFT，将原始时域信号X[k]一分为二，一边是偶次点X₀[k]，另一边是奇次点X₁[k]。

则DFT(X[K])=DFT(X₀[k])+DFT(X₁[k])

根据旋转因子的特性，将公式变形后可以得到

X[k]=F₀[k]+(W_N^nk)F₁[k] k=0,1,2,……,N-1

二点运算

而F[k]又可以同样地往下拆，直到只剩下一个点。
然而一个点的时域数值就是频域数值，由此回溯，最终得到Re[]和Im[]。

蝶形算法

使用DFT去噪声与解调

DFT的应用很多
不如说，DSP就是围绕着傅里叶变换而发展起来的

例如，语音通话去背景噪音、音乐人声去除、语音识别的声学特征提取、雷达索敌、声纳探测、神经网络……

多倒是挺多的，不过要做成实用的还有很多边沿技术和交叉学科要学

看了FFT之后兴冲冲地想去亲手写一次变换程序，结果还是太年轻了，卡在了第一步的位反转分类算法，写了几个版本就是得不到正确的结果。
BASIC版本的又太混乱看不下去，无奈向部长求助，结果是一直把1(0001)当成F(1111)了......
(;´༎ຶД༎ຶ`)

这次就不用Matlab（MatLab功能太强，一些细节问题都帮我们处理了，就不知道问题出在哪了），然后 .NET 又不会，Qt又不会，GTK+也不会
（那我到底会什么(ノへ￣、)。。。。。。）
好吧，会一点Java，那就用swing吧

code

利用一个int数组（double也行，不过Panel的绘图是像素位单位，还要再转回去比较麻烦），保存1024个抽样点的时域值，再创建ReX[]和ImX[]保存直角坐标系下的频域值，MagX[]和PhaseX[]保存极坐标系下的频域值。

需要注意的几个点（坑）：

• IDFT逆离散傅里叶变换之后要乘一个-1，不然图是反的，至今不知道为什么
• 可以定义一个PI常量为3.1415926
• 若ReX[k]为0，则不能求arctan（ReX[k]/ImX[k]），要根据ImX[k]的正负判断相位
• 当MagX[k]小到可以忽略时，相位Phase[k]无意义

image.png

窗口左边的按钮是产生时域波形，使用sin函数生成一个1024长度的数组，右边按钮是FFT，将波形转化为频域信号。

--------------------------分割线----------------------------

先来一个sin波形，频率为（1/512）Hz

F=(1/512)

FFT后：

FFT

虽然不是很明显，但是最左边有一个高高立起的（flag）频率为一的δ冲激信号，其他地方都是0，说明波形中只有一种频率。
（右边一坨的是相位，用于显示突变的点）

换成f=(10/512)后：

F=(10/512)

FFT

频谱往右挪动了一丢丢，但还不是很明显

F=350/512：

F=(350/512)

FFT

很好，已经确定DFT分析成功，但是现在有什么用呢？
对于一种波形DFT确实没什么好做的，但是现实中的波形肯定不止一种，而是多种波形混合而成

比如：
f=3/512混合上f=50/512混合上f=150/512的情况是：

混合信号

假如f=3/512是语音信号（当然不可能这么低，只是举例说明）
f=50/512是载波信号（语音信号的载体）
f=150/512是噪音信号

现在三者一起传过去，怎么样从上面混乱的信号里得到我们想要的声音呢？
不妨先用傅里叶分析看一下频谱：

FFT

for(int i=140;i<MagX.length;i++)
{
MagX[i]=0;
}

纯信号

现在是纯净信号，比之前好看很多

要提取我们的声音信号只需把f=50/512的载波去除就好了，方法和上面一样：

for(int i=40;i<MagX.length;i++)
{
MagX[i]=0;
}

语音信号

显然有一些失真，这是因为在FFT的时候没有给原时域信号加窗函数平滑处理，导致时域分帧后能量泄露，引起边缘混叠，再变回来就变成这样。
可以选择sinc窗函数或者汉明窗函数：
再变回来就好一些了：

平滑处理

总结

程序员的优点除了懒惰外，最大的就是动手的机会多，DFT看似容易，其实很多问题要在实际编程中才能发现，解决细节问题。

傅里叶变换是DSP的核心，一切从这里出发，一切又回到这里。

之前在某乎上看到一句话写的很好：

时域上的变化无常的错综复杂的曲线，只是频域上的一个点
你眼中看似落叶纷飞变化无常的世界，实际只是躺在上帝怀中一份早已谱好的乐章。

最后上一张傅里叶的微笑：

B in your heart