离散傅立叶学习之FFT

FFT，快速傅立叶变换是个蛮有意思的玩意，应用很广，尤其在信息处理方面，最近需要用到，就学学啦。

万事开头难，先做个简单的离散傅立叶变换。

什么是傅立叶变换，就不多解释了，就是把函数打成一堆波，不懂的请出门左转去交重修费。

一，什么是离散傅立叶变换（DFT）

详细内容见：http://zh.wikipedia.org/wiki/%E7%A6%BB%E6%95%A3%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2

虽然我觉得我用图拍它拍的更清楚。

至于离散傅立叶变换，则是应用于实践的对信号采样的傅立叶变换算法，无妨我们先用mathematica生成一个4096位的信号，借用mathematica来了解下离散傅立叶变换，折腾半天终于搞懂了离散傅立叶的工作原理。首先我们通过一个演示来说下.

首先用mathematica做了一个简单的函数发生器。

(*Mathematica Code*)
SignalMk[t_, Omega_] := (
Omega = 200;
f[x_] := Sin[x*Omega/4096] + Sin[2 x*Omega/4096];
Signal = Table[
f[i]
, {i, 1, 4096*t}];(*假设1秒采样4096个 *)
Return[Signal]
)

我们如果取

ListLinePlot[SignalMk[0.1, 200]]

生成一个时长0.1秒，omega为200的图像如下

一个二二的图像，于是我们把它扔到mathematica自带的离散傅立叶变化里面观察下。

下面是老长的代码

Table[
{
ts = 10;
sig = SignalMk[ts, Omega];
Print[“omega is”, Omega];
dr = 100;
P1 = ListLinePlot[
(Abs[Fourier[sig, FourierParameters -> {1, –1}]]
)[[1 ;; dr*ts]],
PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi)/ts,
PlotStyle -> Red
];
P2 = ListLinePlot[
(Re[Fourier[sig, FourierParameters -> {1, –1}]])
[[1 ;; dr*ts]],
PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi])/ts,
PlotStyle -> Blue
];
P3 = ListLinePlot[
(Im[Fourier[sig, FourierParameters -> {1, –1}]]
)[[1 ;; dr*ts]],
PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi)/ts,
PlotStyle -> Yellow
];
Show[P1, P2, P3]
}

生成一系列从omega为100-150-200的傅立叶后的图像，如下

可见在正确的地方出现了尖峰，也就是我们得到了一个漂亮频谱的分布。

那么来解释一下dft（离散傅立叶）吧，

所谓离散傅立叶，就是把离散化的数据点转化为一个离散的傅立叶级数，用来获得频谱分析等，在个领域都有很大用途，其本质就是一堆数据点，按照时序进行。转换为一个频率－振幅的分布，

其中

$$!\omega_k=k\cdot\frac{2\pi}{NT}$$

T为采样频率，N为采样数量，NT为采样总时长，k为编号,k为变换后的样点（在计算机中一般是振幅的标号）

也就是我上面生成频谱图的系数。对应的将是一个振幅，反之我们可以用逆变换把它变回去。

于是原信号就可以表示成

$$\sum_{k=0}^{n-1}Re[f(k)]\cdot Cos(\omega_{k}t)+Im[f(k)]\cdot Sin(\omega_{k}t)$$

无妨一试，又是mathematica

(* Mathematica code*)
sig = SignalMk[0.5, 200];(*生成一个0.5秒长，200omega的信号*)
fx = Fourier[sig];(*把这玩意傅立叶了*)
omgk[k_] := k (2 Pi)/0.5//k对应的omega
y = Table[
Sum[Re[fx[[k + 1]]]*Cos[omgk[k]*t/4096](*注意上面提到的k是0开始编号，mathematica默认1开始，故有＋1，此项无妨认为是示波器的那个“直流成分”*)
+ Im[fx[[k + 1]]]*Sin[omgk[k]*t/4096]
, {k, 1, 500}]
(*这个式子启示就是定义了*)
, {t, 1, 2048}];
ListLinePlot[sig]
ListLinePlot[y]

得到两幅图

分辨不出来那个是信号源的那个是重新逆傅立叶回去的？那就是验证成功了。

二，简单的一维离散傅立叶变换

Now，我们先撇开Mathematica和wiki，自己推倒一下离散傅立叶变换。

打开我们的微积分教材，连续傅立叶变换如是定义

$$!F(\omega)=\frac{1}{l}\intop_{0}^{+l}f(t)e^{-i\omega t}dt$$

这是对于连续函数的积分，当然我们不能直接实现，于是乎，我们得把它离散化，很开心的直接矩阵积个分好啦

我们让采样周期为最小分度dt也就是$$\frac{L}{N}(s)=T$$

那么得到式子

$$!F(\omega)=\frac{1}{l}\sum_{k=0}^{n-1}f(k\cdot T)e^{-i\omega k\cdot T}T$$

其中，k＝0, 此处与wikipedia略有不同，归一化得到

$$!F[l]=\sum_{k=0}^{k-1}f(k)e^{-i\frac{2\pi}{N}k\cdot l}$$

此时需要注意的是得到的系数并非我们最后使用的，中间还有一个归一化系数，不过这个暂且无关紧要，据此，我们可以写出第一个简单的o(n^2)的dft程序

首先，c等语言是不支持虚数的，浩浩又不想学fortran，于是先转为两个部分

$$!F[i].Re=\sum_{k=0}^{k-1}f(k)\cdot Cos(2\cdot\pi\cdot k\cdot i/N)$$

$$!F[i].Im=-\sum_{k=0}^{k-1}f(k)\cdot Sin(2\cdot\pi\cdot k\cdot i/N)$$

$$!F[i].Abs=Sqrt(F[i].Re^2+F[i].Im^2)$$

于是乎，rush一个小c的程序，就是我们的慢速傅立叶了，为了考虑向单片机迁移，代码风格是纯c的风格，大量指针毫无结构毫无美感，各位将就看哈。

#include“stdio.h”
#include“stdlib.h”
#include“math.h”
#define pi 3.1415926
float* dft(float *input,int Size,float *>re,float *>im)//return abs 慢速dft核心
{
float *abs=(float*)malloc(Size*sizeof(float));
re=(float *)malloc(Size*sizeof(float));
im=(float *)malloc(Size*sizeof(float));
for(int i=0;i<Size;i++)
{
float sumre=0;
float sumim=0;
for(int k=0;k<Size;k++)
{
sumre+=input[k]*cosf(2*pi*k*i/Size);//calculate re
sumim-=input[k]*sinf(2*pi*k*i/Size);
}
if(sumre>4000)
printf(“%d\n”,sumre);
re[i]=sumre/sqrtf(Size);
im[i]=sumim/sqrtf(Size);//取归一化系数都是sqrt(n)
abs[i]=sqrtf(sumre*sumre+sumim*sumim);
}
return abs;
}
float* SignalMk(float time,int omega)//信号发生器
{
int Size=(int)time*4096;
float *sig=(float *)malloc(Size*sizeof(float));
for(int i=0;i<Size;i++)
{
sig[i]=sinf( ((float)i)*omega/4096)+sinf( 2*((float)i)*omega/4096);
}
return sig;
}
int main()
{
float *re,*im;
float *sig=SignalMk(1,200);
printf(“Build Sucessful\n”);
float *abs=dft(sig,4096,re,im);
printf(“dft Sucessful\n”);
FILE *fp=fopen(“d.txt”,“w”);
for(int i=0;i<4096;i++)
fprintf(fp,“%f “,abs[i]);
}

于是，结果导入了mathematica得到了下面的图

So cool!

这就是那个简单的c程的结果。某种意义上，我们已经实现了多数的信号处理。

三，主题，FFT！效率，还是效率

对于数字信号处理，跑在12MHz单片机吃着无可忍耐的破烂环境还要干最重的活，我真提单片机程序员感到捉急（以前玩单片机的热情就是被硬件资源熄灭的，相比而言我更愿意在8 cores主机写非结构网格CFD）。于是，这个效率是$$O(n^2)$$简直不能忍有没有。。。设想我们要分析一个人的声音来确定他的身份，傅立叶变换是基础了，可是！mp3采样率是44,100 Hz ，也就是说，

某阴森的地下会场，男A道，“你好。”，耗时一秒

男B暗道，待我分析下它的声音是否正确。

经过了漫长的162秒（用我们的slowdft需要44100*44100次运算，51单片机需要162秒）。本拉登已经成功脱身。

或者

苦苦思考，柯南想到了证人是谁，打开变声器，说了一段60秒钟的推理，经过单片机自燃的温度的162小时运算，

犯人逃跑了。

这样不如找块豆腐撞死好了。于是乎，fft横空出世。

FFT是一种复杂度为NLog(N)的算法，一般看到这个效率，就会很明了这是一个二分构造的过程，整个算法必然存在一定的对称性，而且几乎可以认定是存在递归过程。剩下要做的就是找到那个简单的递归函数，然后得到我们的运算结果。不过在这里，算法对称性倒是由数学给出的。

参考：http://zh.wikipedia.org/wiki/%E5%BF%AB%E9%80%9F%E5%82%85%E9%87%8C%E5%8F%B6%E5%8F%98%E6%8D%A2

我们用$$W_N$$表示$$e^{-j\frac{2\pi}{N}}$$，于是乎

$$!W_{N}^{k+N}=W_{N}^{k}$$

$$!W_{N}^{k+\frac{N}{2}}=-W_{N}^{k}$$

$$!W_{N}^{m\cdot k_{n}}=W_{\frac{N}{m}}^{k_{n}}$$

根据这个，我们可以得到这样一个递归构造过程（在此我倒更宁愿用计算机上的方式而非纯粹数学语言表述）。

一般，对于此类问题都是先分解为朴素递归，然后找到重复使用的调用（比如完全相同的参数输入），得到结果，是为计划搜索。

F(N,k)=

If(k<=$$\frac{N}{2}$$)

$$F_{even}(k)+W_N^kF_{odd}(k)$$

If(k>$$\frac{N}{2}$$)

$$F_{even}(k-\frac{N}{2})-W_N^{k-\frac{N}{2}}F_{odd}(k-\frac{N}{2})$$

由此可以递归构造出F(k)

于是，我写了这样一个递归

Complex dfftk(float *input ,int Size,int k)
{
float*even=(float *)malloc(sizeof(float)*Size/2);
float*odd=(float *)malloc(sizeof(float)*Size/2);
for(int i=0;i<Size;i++)//筛选器
if(i%2==0)
even[i/2]=input[i];
else
odd[i/2]=input[i];
if(Size==1)
return input[0]*Complex(1,0);
if(k<=Size/2)
return (dfftk(even,Size/2,k)+W(Size,k)*dfftk(odd,Size/2,k));
else
{
return (dfftk(even,Size/2,k–Size/2)-W(Size,k–Size/2)*dfftk(odd,Size/2,k–Size/2));
}
return Complex(0,0);
}

然后主函数这么搞

Complex *fft(float*input,int Size)
{
Complex*res=(Complex *)malloc(sizeof(Complex)*Size);
int i;
for(i=0;i<Size;i++)
res[i]=dfftk(input,Size,i);
return res;
}

于是得到了漂亮的结果

嗯！这就是我们想要的，现在万事俱备，只欠东风：如何让它快起来。

四：雕虫小技

让我们在回顾下我们的函数输入，假定计算的都是2的指数，函数的结果依赖于

1.点集合：input，包含Size信息。

2.k的值。

so，点集合是一个颇为虚无飘渺的概念，那我们怎么确定它呢，注意到我们假定数量一直是二的指数，而且一直二分，那么，只要说一个点集合的Size，我们就能确定二分了多少次，给出第一个点的全局编号，我们就可以在这棵查找树上确定他的位置。于是，对于每次搜索，我们可以用以下几个点标名其意义

深度，首位编号，以及k。

在计算中，我们要做对半优化，其实只需在k<Size/2时候把+Size/2加上。

于是乎：

bool Travel[4096][13][4096]={0};//纪录是否访问过
Complex val[4096][13][4096];//纪录访问过的值
int numheight[15]={0};
Complex dfftk(float *input,int*nums ,int Size,int k)
{
int h=log(Size)/log(2);
int top;
top=nums[0];
if(Travel[top][h][k])//访问过就直接用
return Complex(val[top][h][k]);
Travel[top][h][k]=1;//标定为访问过
float*even=(float *)malloc(sizeof(float)*Size/2);
int*evenums=(int *)malloc(sizeof(int)*Size/2);
float*odd=(float *)malloc(sizeof(float)*Size/2);
int*oddnums=(int *)malloc(sizeof(int)*Size/2);
for(int i=0;i<Size;i++)//奇偶分流
if(i%2==0)
{
even[i/2]=input[i];
evenums[i/2]=nums[i];
}
else
{
odd[i/2]=input[i];
oddnums[i/2]=nums[i];
}
if(Size==1)
{
val[nums[0]][h][k]=(input[0]*Complex(1,0));
return Complex(val[top][h][k]);
}
if(k<=Size/2)
{
Complex everes,oddres;
everes=dfftk(even,evenums,Size/2,k);
oddres=dfftk(odd,oddnums,Size/2,k);
val[top][h][k]=(everes+W(Size,k)*oddres);
Travel[top][h][k+Size/2]=1;
val[top][h][k+Size/2]=(everes–W(Size,k)*oddres);
return Complex(val[top][h][k]);
}
printf(“You are Wrong\n”);
return Complex(0,0);
}
Complex *fft(float*input,int Size)
{
Complex*res=(Complex *)malloc(sizeof(Complex)*Size);
int i;
int *nums=(int *)malloc(sizeof(int)*Size);
for(i=0;i<Size;i++)
nums[i]=i;
for(i=0;i<Size;i++)
res[i]=dfftk(input,nums,Size,i);
return res;
}

跑一下，结果未变，时间从24646变成了189毫秒。。。。。这差距

于是，时间得到了优化，再回头看我们的数组，有什么不对么？

$$N^2*Log_2N$$的空间？你在逗我嘛？

五，丧心病狂

状态压缩：

首先看这个记录，我们此时需要计算一下我们真正需要记录的数量，从一到k次次二分，共有$$NLog_2N$$个状态，其中有一半是作为k<Size/2return给出的，一半是作为记录的。于是我们需要精确构造一个新的记录，使得这个存储量大幅度减少。