离散傅立叶学习之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做了一个简单的函数发生器。

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

我们如果取

  1. ListLinePlot[SignalMk[0.1, 200]]

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

anal

 

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

下面是老长的代码

  1. Table[
  2.  {
  3.   ts = 10;
  4.   sig = SignalMk[ts, Omega];
  5.   Print["omega is", Omega];
  6.   dr = 100;
  7.   P1 = ListLinePlot[
  8.     (Abs[Fourier[sig, FourierParameters -> {1, -1}]]
  9.       )[[1 ;; dr*ts]],
  10.     PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi)/ts, 
  11.     PlotStyle -> Red
  12.     ];
  13.   
  14.   P2 = ListLinePlot[
  15.     (Re[Fourier[sig, FourierParameters -> {1, -1}]])
  16.      [[1 ;; dr*ts]],
  17.     PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi])/ts, 
  18.     PlotStyle -> Blue
  19.     ];
  20.   P3 = ListLinePlot[
  21.     (Im[Fourier[sig, FourierParameters -> {1, -1}]]
  22.       )[[1 ;; dr*ts]],
  23.     PlotRange -> All, DataRange -> {1, dr*ts}*(2 Pi)/ts, 
  24.     PlotStyle -> Yellow
  25.     ];
  26.   
  27.   Show[P1, P2, P3]
  28.   
  29.   }

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

50fly 100fly 200fly

 

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

那么来解释一下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

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

得到两幅图

redft0 redft1

 

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

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

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的风格,大量指针毫无结构毫无美感,各位将就看哈。

  1. #include"stdio.h"
  2. #include"stdlib.h"
  3. #include"math.h"
  4. #define pi 3.1415926
  5. float* dft(float *input,int Size,float *>re,float *>im)//return abs 慢速dft核心
  6. {
  7.     float *abs=(float*)malloc(Size*sizeof(float));
  8.     re=(float *)malloc(Size*sizeof(float));
  9.     im=(float *)malloc(Size*sizeof(float));
  10.     for(int i=0;i<Size;i++)
  11.     {
  12.         float sumre=0;
  13.         float sumim=0;
  14.         for(int k=0;k<Size;k++)
  15.         {
  16.             sumre+=input[k]*cosf(2*pi*k*i/Size);//calculate re
  17.             sumim-=input[k]*sinf(2*pi*k*i/Size);
  18.         }
  19.         if(sumre>4000)
  20.             printf("%d\n",sumre);
  21.         re[i]=sumre/sqrtf(Size);
  22.         im[i]=sumim/sqrtf(Size);//取归一化系数都是sqrt(n)
  23.         abs[i]=sqrtf(sumre*sumre+sumim*sumim);
  24.     }
  25.     return abs;
  26. }
  27. float* SignalMk(float time,int omega)//信号发生器
  28. {
  29.     int Size=(int)time*4096;
  30.     float *sig=(float *)malloc(Size*sizeof(float));
  31.     for(int i=0;i<Size;i++)
  32.     {
  33.         sig[i]=sinf( ((float)i)*omega/4096)+sinf( 2*((float)i)*omega/4096);
  34.     }
  35.     return sig;
  36. }
  37. int main()
  38. {
  39.     float *re,*im;
  40.     float *sig=SignalMk(1,200);
  41.     printf("Build Sucessful\n");
  42.     float *abs=dft(sig,4096,re,im);
  43.     printf("dft Sucessful\n");
  44.     FILE *fp=fopen("d.txt","w");
  45.     for(int i=0;i<4096;i++)
  46.         fprintf(fp,"%f ",abs[i]);
  47. }

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

slowdft

 

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)

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

  1. Complex dfftk(float *input ,int Size,int k)
  2. {
  3.     float*even=(float *)malloc(sizeof(float)*Size/2);
  4.     float*odd=(float *)malloc(sizeof(float)*Size/2);
  5.     for(int i=0;i<Size;i++)//筛选器
  6.         if(i%2==0)
  7.             even[i/2]=input[i];
  8.         else
  9.             odd[i/2]=input[i];
  10.     if(Size==1)
  11.         return input[0]*Complex(1,0);
  12.     if(k<=Size/2)
  13.         return (dfftk(even,Size/2,k)+W(Size,k)*dfftk(odd,Size/2,k));
  14.     else
  15.     {
  16.         return (dfftk(even,Size/2,k-Size/2)-W(Size,k-Size/2)*dfftk(odd,Size/2,k-Size/2));
  17.     }
  18.     return Complex(0,0);
  19. }

然后主函数这么搞

  1. Complex *fft(float*input,int Size)
  2. {
  3.     Complex*res=(Complex *)malloc(sizeof(Complex)*Size);
  4.     int i;
  5.     for(i=0;i<Size;i++)
  6.         res[i]=dfftk(input,Size,i);
  7.     return res;
  8. }

于是得到了漂亮的结果

fft0

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

四:雕虫小技

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

1.点集合:input,包含Size信息。

2.k的值。

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

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

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

于是乎:

  1. bool Travel[4096][13][4096]={0};//纪录是否访问过
  2. Complex val[4096][13][4096];//纪录访问过的值
  3. int numheight[15]={0};
  4.  
  5. Complex dfftk(float *input,int*nums ,int Size,int k)
  6. {
  7.     int h=log(Size)/log(2);
  8.     int top;
  9.     top=nums[0];
  10.     if(Travel[top][h][k])//访问过就直接用
  11.         return Complex(val[top][h][k]);
  12.     Travel[top][h][k]=1;//标定为访问过
  13.     float*even=(float *)malloc(sizeof(float)*Size/2);
  14.     int*evenums=(int *)malloc(sizeof(int)*Size/2);
  15.     float*odd=(float *)malloc(sizeof(float)*Size/2);
  16.     int*oddnums=(int *)malloc(sizeof(int)*Size/2);
  17.     for(int i=0;i<Size;i++)//奇偶分流
  18.         if(i%2==0)
  19.         {
  20.             even[i/2]=input[i];
  21.             evenums[i/2]=nums[i];
  22.         }
  23.         else
  24.         {
  25.             odd[i/2]=input[i];
  26.             oddnums[i/2]=nums[i];
  27.         }
  28.     if(Size==1)
  29.     {
  30.         val[nums[0]][h][k]=(input[0]*Complex(1,0));
  31.         return Complex(val[top][h][k]);
  32.     }
  33.     if(k<=Size/2)
  34.     {
  35.         Complex everes,oddres;
  36.         everes=dfftk(even,evenums,Size/2,k);
  37.         oddres=dfftk(odd,oddnums,Size/2,k);
  38.         val[top][h][k]=(everes+W(Size,k)*oddres);
  39.         Travel[top][h][k+Size/2]=1;
  40.         val[top][h][k+Size/2]=(everes-W(Size,k)*oddres);
  41.         return Complex(val[top][h][k]);
  42.     }
  43.     printf("You are Wrong\n");
  44.     return Complex(0,0);
  45. }
  46. Complex *fft(float*input,int Size)
  47. {
  48.     Complex*res=(Complex *)malloc(sizeof(Complex)*Size);
  49.     int i;
  50.     int *nums=(int *)malloc(sizeof(int)*Size);
  51.     for(i=0;i<Size;i++)
  52.         nums[i]=i;
  53.     for(i=0;i<Size;i++)
  54.         res[i]=dfftk(input,nums,Size,i);
  55.     return res;
  56. }

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

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

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

五,丧心病狂

状态压缩:

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

首先,当然是毫无节操的直接stl::Map了(本来准备自己写个hash来着。。困了。。。去他的纯c平台)。。。于是。