此贴已删除

kexuexyz

哥花了一周的时间彻底搞通了FFT，特地写了篇日志分享一下研究成果。 我日，图片贴不上来

miaoguoqiang

粘贴的
-----------------------------------------

快速傅里叶变换FFT的C语言算法彻底研究fficeffice" />

LED音乐频谱显示的核心算法就是快速傅里叶变换，FFT的理解和编程还是比较难的，特地撰写此文分享一下研究成果。

 一、彻底理解傅里叶变换

快速傅里叶变换（Fast Fourier Transform)是离散傅里叶变换的一种快速算法，简称FFT，通过FFT可以将一个信号从时域变换到频域。

模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。

假设采样频率为fs，采样点数为N，那么FFT结果就是一个N点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/N。

举例说明：用1kHz的采样频率采样128点，则FFT结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 Hz。

这个频率点的幅值为：该点复数的模值除以N/2（n=1时是直流分量，其幅值是该点的模值除以N）。

二、傅里叶变换的C语言编程

1、对于快速傅里叶变换FFT，第一个要解决的问题就是码位倒序。

    假设一个N点的输入序列，那么它的序号二进制数位数就是t=log2N.

    码位倒序要解决两个问题：①将t位二进制数倒序；②将倒序后的两个存储单元进行交换。

①将t=3位二进制数倒序

②将倒序后的两个存储单元进行交换

 如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢？利用C语言的移位功能！

程序如下，我不多说，看不懂者智商一定在180以下！

复数类型定义及其运算

#define N 64      //64点

#define log2N 6   //log2N=6

/*复数类型*/

typedef struct

{

 float real;

 float img;

}complex;

complex xdata x[N]; //输入序列

/*复数加法*/

complex add(complex a,complex b)

{

 complex c;

 c.real=a.real+b.real;

 c.img=a.img+b.img;

 return c;

}

/*复数减法*/

complex sub(complex a,complex b)

{

 complex c;

 c.real=a.real-b.real;

 c.img=a.img-b.img;

 return c;

}

/*复数乘法*/

complex mul(complex a,complex b)

{

 complex c;

 c.real=a.real*b.real - a.img*b.img;

 c.img=a.real*b.img + a.img*b.real;

 return c;

}

/***码位倒序函数***/

void Reverse(void)

{

 unsigned int i,j,k;

 unsigned int t;

 complex temp;//临时交换变量

 for(i=0;i<N;i++)//从第0个序号到第N-1个序号

 {

  k=i;//当前第i个序号

  j=0;//存储倒序后的序号，先初始化为0

  for(t=0;t<log2N;t++)//共移位t次，其中log2N是事先宏定义算好的

  {

   j<<=1;

   j|=(k&1);//j左移一位然后加上k的最低位

   k>>=1;//k右移一位，次低位变为最低位

  }

  if(j>i)//如果倒序后大于原序数，就将两个存储单元进行交换（判断j>i是为了防止重复交换）

  {

   temp=x;

   x=x[j];

   x[j]=temp;

  }

 }

}

2、第二个要解决的问题就是蝶形运算

①第1级（第1列）每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，

 第m级蝶形运算，每个蝶形的两节点“距离”L=2m-1。

②对于16点的FFT，第1级有16组蝶形，每组有1个蝶形；第2级有4组蝶形，每组有2个蝶形；第3级有2组蝶形，每组有4个蝶形；第4级有1组蝶形，每组有8个蝶形。由此可推出，

对于N点的FFT，第m级有N/2L组蝶形，每组有L=2m-1个蝶形。

③旋转因子的确定

以16点FFT为例，第m级第k个旋转因子为，其中k=0～2m-1-1，即第m级共有2m-1个旋转因子，根据旋转因子的可约性，，所以第m级第k个旋转因子为，其中k=0～2m-1-1。

为提高FFT的运算速度，我们可以事先建立一个旋转因子数组，然后通过查表法来实现。

complex code WN[N]=//旋转因子数组
{ //为节省CPU计算时间，旋转因子采用查表处理
  //★根据实际FFT的点数N，该表数据需自行修改
  //以下结果通过Excel自动生成
  //  WN[k].real=cos(2*PI/N*k);
  //  WN[k].img=-sin(2*PI/N*k);

  {1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
  {0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
  {0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
  {0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
  {0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
  {-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
  {-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
  {-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
  {-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
  {-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
  {-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
  {-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
  {0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
  {0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
  {0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
  {0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}    
};

3、算法实现

我们已经知道，N点FFT从左到右共有log2N级蝶形，每级有N/2L组，每组有L个。所以FFT的C语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算（整个FFT共log2N级），中间层循环完成每一组的蝶形运算（每一级有N/2L组），最内层循环完成单独1个蝶形运算（每一组有L个）。

  /***【快速傅里叶变换】***/

  void FFT(void)

  {

   unsigned int i,j,k,l;

   complex top,bottom,xW;

   Reverse(); //码位倒序

   for(i=0;i<log2N;i++)   /*共log2N级*/

   {    //一级蝶形运算

    l=1<<i;//l等于2的i次方

    for(j=0;j<N;j+=2*l)  /*每L个蝶形是一组，每级有N/2L组*/

    {   //一组蝶形运算

     for(k=0;k<l;k++)   /*每组有L个*/

     {  //一个蝶形运算

       xW=mul(x[j+k+l],WN[N/(2*l)*k]); //碟间距为l

       top=add(x[j+k],xW); //每组的第k个蝶形

       bottom=sub(x[j+k],xW);

       x[j+k]=top;

       x[j+k+l]=bottom;

     }

    }

   }

  }

三、FFT计算结果验证

随便输入一个64点序列，例如
x[N]={{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0}};
在Keil中Debug查看数组变量x的FFT计算结果并与MATLAB计算结果进行比对，可以发现非常准确，说明程序编写正确！

     FFT验证：C语言与MATLAB计算结果对比

正点原子

直接发论坛吧...链来链去的,不好.

kexuexyz

回复【2楼】正点原子:
直接发论坛吧...链来链去的,不好.
---------------------------------
有图片，粘贴后不显示

正点原子

回复【3楼】kexuexyz:
---------------------------------
复制一下,上传就可以了.

正点原子

qq空间的图片拒绝外链,狗日的TX...

Pony279

回复【6楼】正点原子:
---------------------------------
TX估计是考虑到安全性吧, 毕竟QQ空间和新浪微博不一样, 新浪微博是谁都可以看的, 而QQ空间里面的东西就有限制了

SprayFly

顶一个！谢谢楼主

the_ether

回复【楼主位】kexuexyz:
---------------------------------
额。图片不能看，能否把这篇资料的WORD发我邮箱527717930@QQ.COM

tian123chi

上传个VC++ 6.0程序 
/*时间抽选基2 FFT及IFFT算法C语言实现*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000
/*定义复数类型*/
typedef struct{
 double real;
 double img;
}complex;

void fft();    /*快速傅里叶变换*/
void ifft();
void initW();  /*初始化变换核*/
void change(); /*变址*/
void add(complex ,complex ,complex *); /*复数加法*/
void mul(complex ,complex ,complex *); /*复数乘法*/
void sub(complex ,complex ,complex *); /*复数减法*/
void divi(complex ,complex ,complex *);/*复数除法*/
void output();                            /*输出结果*/

complex x[N], *W; /*输入序列,变换核*/
int size_x=0;     /*输入序列的大小，在本程序中仅限2的次幂*/
double I;        /*圆周率*/

int main(){
 int i,method;
 void fft();    /*快速傅里叶变换*/
 void ifft();
 void initW();  /*初始化变换核*/
 void change(); /*变址*/
 void add(complex ,complex ,complex *); /*复数加法*/
 void mul(complex ,complex ,complex *); /*复数乘法*/
 void sub(complex ,complex ,complex *); /*复数减法*/
 void divi(complex ,complex ,complex *);/*复数除法*/
 void output();                            /*输出结果*/
 system("cls");
 I=atan(1)*4;
 printf("Please input the size of x:\n");
 scanf("%d",&size_x);
 printf("Please input the data in x[N]:\n");
 for(i=0;i<size_x;i++)
  scanf("%lf%lf",&x.real,&x.img);
 initW();
 printf("Use FFT(0) or IFFT(1)?\n");
 scanf("%d",&method);
 if(method==0)
  fft();
 else
  ifft();
 output();
 return 0;
}

/*快速傅里叶变换*/
void fft(){
 int i=0,j=0,k=0,l=0;
 complex up,down,product;
 change();
 for(i=0;i< log(size_x)/log(2) ;i++){  /*一级蝶形运算*/
  l=1<<i;
  for(j=0;j<size_x;j+= 2*l ){            /*一组蝶形运算*/
   for(k=0;k<l;k++){       /*一个蝶形运算*/
     mul(x[j+k+l],W[size_x*k/2/l],&product);
     add(x[j+k],product,&up);
     sub(x[j+k],product,&down);
     x[j+k]=up;
     x[j+k+l]=down;
   }
  }
 }
}



/*初始化变换核*/
void initW(){
 int i;
 W=(complex *)malloc(sizeof(complex) * size_x);
 for(i=0;i<size_x;i++){
  W.real=cos(2*PI/size_x*i);
  W.img=-1*sin(2*PI/size_x*i);
 }
}

/*变址计算，将x(n)码位倒置*/
void change(){
 complex temp;
 unsigned short i=0,j=0,k=0;
 double t;
 for(i=0;i<size_x;i++){
  k=i;j=0;
  t=(log(size_x)/log(2));
  while( (t--)>0 ){
   j=j<<1;
   j|=(k & 1);
   k=k>>1;
  }
  if(j>i){
   temp=x;
   x=x[j];
   x[j]=temp;
  }
 }
}

/*输出傅里叶变换的结果*/
void output(){
 int i;
 printf("The result are as follows\n");
 for(i=0;i<size_x;i++){
  printf("%.4f",x.real);
  if(x.img>=0.0001)printf("+%.4fj\n",x.img);
  else if(fabs(x.img)<0.0001)printf("\n");
  else printf("%.4fj\n",x.img);
 }
}
void add(complex a,complex b,complex *c){
 c->real=a.real+b.real;
 c->img=a.img+b.img;
}

void mul(complex a,complex b,complex *c){
 c->real=a.real*b.real - a.img*b.img;
 c->img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c){
 c->real=a.real-b.real;
 c->img=a.img-b.img;
}
void divi(complex a,complex b,complex *c){
 c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img);
 c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);
}


/*************************************
比如求fft([1,2j,3,4j])
运行过程：
Please input the size of x:
4
Please input the data in x[N]:
1 0
0 2
3 0
0 4
Use FFT(0) or IFFT(1)?
0
The result are as follows
4.0000+6.0000j
-4.0000
4.0000-6.0000j
0.0000
****************************************/

kexuexyz

回复【11楼】tian123chi:
---------------------------------
你这个程序计算结果不对，虚部相差一个符号！不知哪里出问题了？

Tardis

其实我很想问，神马是傅立叶变换？
我大学的时候因为这个东西数字信号处理一科目清考的说。

tian123chi

回复【14楼】Tardis:
---------------------------------
傅里叶变换是从时域转换成频域的变换，跟波形打交道的基本都用！

tian123chi

回复【12楼】kexuexyz:
---------------------------------
我跟matlab对比了几次是对的，当然用的是8点的。

Tardis

回复【15楼】tian123chi:
---------------------------------
哦，原来如此

tian123chi

回复【17楼】Tardis:
---------------------------------
玩电源的应该懂得把!

Tardis

回复【18楼】tian123chi:
---------------------------------
问题是，我业余不玩电源的

liyanchao153

回复【楼主位】kexuexyz:
---------------------------------
lz  你的这篇日志地址是什么   我想看你发表这表文章的图片

liklon

感觉好复杂的样子

hqr1021

只要涉及信号的处理，离不开傅里叶。FFT程序貌似《数字信号处理》已经有了。实时数据处理最好用dsp吧~

小小一书生

最近正烦扰这个，谢谢分享哈！！！

倒拔萝卜

回复【5楼】miaoguoqiang:
---------------------------------
图片呢？麻花藤又不听话了

lsj9383

回复【19楼】Tardis:
---------------------------------
不仅仅是电源吧，很多领域，都会涉及到信号处理。

Tardis

回复【26楼】lsj9383:
---------------------------------
问题是，傅立叶变换有专用的DSP芯片，不用人手来弄，可以硬解

Tardis

回复【26楼】lsj9383:
---------------------------------
就像HDMI一样，哪有人手动解码的？S5PV210芯片，三星官方就有现成驱动，直接调用就得，人家的驱动和STM32的库函数一样，你不用管内部是怎么写的，直接调用，调用了，功能就出了

lsj9383

回复【28楼】Tardis:
---------------------------------
无法抓住核心理论。
顶多就是将各个函数和模块衔接起来的体力活。
这叫搞起耍。

小小一书生

回复【11楼】tian123chi:
---------------------------------
你这个程序编译什么的都行，就是运行结果出不来，这是哪里出错了？

shendade

感觉好详细的样子，我现在是初级，还不太懂fft的过程，所以困惑有待解决。

kuangnan

正在捣鼓傅立叶，通常学会单片机可能很简单，其实难的是你在处理实际的工程数据，也就是数据处理，才发现其实我们学得远远不够呀。我在坛子里共享了“漫画学习傅立叶解析”。
今天看到贵贴，也来学习学习，希望坛子里能有更深，更广的学习交流。

qingkongzlr

回复【9楼】qq371833846:
---------------------------------
想问问有DSP库吗？我在官网没有找到，急需，可以发到我的邮箱吗？791141003@qq.com,十分感谢。

临海观雪

求链接，这图看不到

闪客小子

回复【6楼】正点原子:
---------------------------------
原子哥，我遇到一个问题，没搞懂，我在做FFT时，移植的人家的程序，最后用串口打印出来发现数据是对称的，数组是用C语言生成的，生成1024个点，然后相当于1024HZ的AD在采集，当我设置1HZ为1时，他的数组的最后一个也就是相当于1023HZ的幅值也是1，这就是对称吧，这个问题不知道怎么办，如果我本来就有1HZ和1023HZ的信号，那我怎么知道哪一个才是真的，呜呜呜呜，谢谢原子哥

正点原子

回复【35楼】闪客小子:
---------------------------------
这个我也不太懂。

2593

做音乐频谱分析的时候根据采样定理对声音信号进行采样的采样频率fs>=40KHZ，那么根据楼主说的如果采样128点的话，那么FFT之后序列是关于中间点对称的，按照楼主的说法从第64点之后代表的频率分量是大于20KHZ？既然是对称的那么说明后面的64点是没有用的，真正有用的是前面的64点吧？后面的64点应该只是原频率点的周期延拓吧？

狼烟四起

不错，做数据和波形分析用得上。

adaphoto

MARK，好贴要顶

文刀言身寸

在程佩青的《数字信号处理教程》中有一个对FFT用C实现的程序，而且还有流程图。不过里面处理的数据好像都是比较简单的，没涉及到复数。

伤剑使者

真心是 搞不懂傅立叶，没办法

ALIEN_WALI

回复【3楼】kexuexyz:
---------------------------------
在哪能看完整版的？链接在哪里？

合肥-文盲

顶一个  楼主搞个带图的完整版本的呀

涛涛

我现在用fft 通过stm32 F103 AD 采集出来了，送到FFT运算，输入分别是 市电三相 电压和三相电流，经过FFT 计算出各相的 电压 电流，基波 和各次谐波。实部 和虚部，。怎么计算电压，电流 的 有功功率 和无功功率 ，和功率因数怎么算啊。能给个例程吗

awkwardeli

回复【5楼】miaoguoqiang:
---------------------------------
不知道为什么，我使用这个程序做了8点FFT,求出来的结果和matlab不一样，，求指教

op417583795

这是错的，采样率完全没用到啊

1860xq

感谢分享！！

budy

请问有没有图啊  想结合图更直观点

budy

miaoguoqiang 发表于 2013-3-19 15:18
粘贴的
-----------------------------------------

没有图 不是看得很懂 求图

sunshine

正点原子 发表于 2013-3-19 00:16
回复【3楼】kexuexyz:
---------------------------------
复制一下,上传就可以了.

挖坟。。。想不到现在出来工作了还有机会来同学院师兄原子哥这里......最近工作在做FFT相关的东西

heart小十一

好的  看看 学习