FFT最详细的源代码和解释.docx

我自己的一些详细标注，有利于深入了解FFT，后面附加几位网友对FFT的理解及源代码，让广大朋友更迅速的掌握FFT#include DSP281x_Device.h / DSP281x Headerfile Include File，添加所有头文件#include DSP281x_Examples.h / DSP281x Examples Include File，条件编译而已#include f2812a.h /一些变量的宏定义而已#includemath.h#define PI 3.1415926 /前变后常#define SAMPLENUMBER 128/#define SAMPLENUMBER 512void InitForFFT();void MakeWave();/void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER);int INPUTSAMPLENUMBER,DATASAMPLENUMBER;float fWaveRSAMPLENUMBER,fWaveISAMPLENUMBER,wSAMPLENUMBER;float sin_tabSAMPLENUMBER,cos_tabSAMPLENUMBER;/逐级计算FFT，一级一级递推void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER)int x0,x1,x2,x3,x4,x5,x6,xx;int i,j,k,b,p,L;float TR,TI,temp;/* following code invert sequence */倒序for ( i=0;iSAMPLENUMBER;i+ )/就是码位倒置嘛，二进制各个位独立出来再反向 /128七位二进制表示，/号代表右移嘛x0=x1=x2=x3=x4=x5=x6=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;/最低位，最高位反过来dataIxx=dataRi;for ( i=0;iSAMPLENUMBER;i+ )dataRi=dataIi; dataIi=0; /对应过来/* following code FFT */for ( L=1;L0 ) b=b*2; i-; /* b= 2(L-1) */for ( j=0;j0 ) /* p=pow(2,7-L)*j; */p=p*2; i-;p=p*j;for ( k=j;k128;k=k+2*b ) /* for (3) */TR=dataRk; TI=dataIk; temp=dataRk+b;dataRk=dataRk+dataRk+b*cos_tabp+dataIk+b*sin_tabp;dataIk=dataIk-dataRk+b*sin_tabp+dataIk+b*cos_tabp;dataRk+b=TR-dataRk+b*cos_tabp-dataIk+b*sin_tabp;dataIk+b=TI+temp*sin_tabp-dataIk+b*cos_tabp; /递推嘛，防止立马调用结果， /* END for (3) */ /引入一个中间变量存原始值， /* END for (2) */ /防止上一步对下一步的影响 /* END for (1) */for ( i=0;iSAMPLENUMBER/2;i+ ) /对称性，前半部分即可 wi=sqrt(dataRi*dataRi+dataIi*dataIi); /* END FFT */main()int i;InitForFFT();MakeWave();for ( i=0;iSAMPLENUMBER;i+ )fWaveRi=INPUTi;fWaveIi=0.0f;wi=0.0f;FFT(fWaveR,fWaveI);/输入波形进行FFT变换,此处引入起始实参即可递推下去for ( i=0;iSAMPLENUMBER;i+ )DATAi=wi;/变换后的波形转换到输出接口while ( 1 );/ break point/旋转因子事先初始化好，方便调用void InitForFFT()int i;for ( i=0;iSAMPLENUMBER;i+ )sin_tabi=sin(PI*2*i/SAMPLENUMBER);cos_tabi=cos(PI*2*i/SAMPLENUMBER);/旋转因子事先初始化好，方便调用 /利用这个，确实能产生各种各样的谐波void MakeWave() /利用这个，确实能产生各种各样的谐波int i;for ( i=0;iSAMPLENUMBER;i+ )/1024是相应的幅值嘛，只要弄出个基波，以之为标准即可INPUTi=sin(PI*2*i/SAMPLENUMBER)*1024 +sin(PI*2*i/SAMPLENUMBER*3)*1024/3 +sin(PI*2*i/SAMPLENUMBER*5)*1024/5 +sin(PI*2*i/SAMPLENUMBER*7)*1024/7 +sin(PI*2*i/SAMPLENUMBER*9)*1024/9; /INPUTi=sin(PI*2*i/SAMPLENUMBER*3)*1024;FFT原理及实现(Radix-2)哈!经过连续几个晚上的奋战,终于弄懂了FFT推导过程及实现!HappyJ基2 FFT总的思想是将输入信号对半分割,再对半分割,再再对半分割(以下省略10000个再再.J)直至分割到2点.两点DFT简化假设输入为x0,x1;输出为X0,X1.伪代码如下:/ -#defineN2#definePI3.1415926/ -inti, jfor(i=0, Xi=0.0; iN; i+)for(j=0; jN; j+)Xi += xj * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );注意到(我想Audio编解码很多时候都是对cos,sin进行优化!)j=0j=1i=0cos1sin0tw1cos1sin0tw1i=1cos1Sin0tw1cos-1sin0tw-1X0=x0*(1-0) + x1*(1-0)=x0 + 1*x1;X1=x0*(1-0) + x1*(-1-0) =x0 - 1*x1;这就是单个2点蝶形算法.FFT实现流程图分析(N=8,以8点信号为例)FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-point DFTs8点FFT流程图(Layer表示层, gr表示当前层的颗粒)下面以LayerI为例.LayerI部分,具有4个颗粒,每个颗粒2个输入(注意2个输入的来源,由时域信号友情提供,感谢感谢J)我们将输入xk分为两部分x_rk, x_ik.具有实部和虚部,时域信号本没有虚部的,因此可以让x_ik为0.那么为什么还要画蛇添足分为实部和虚部呢?这是因为LayerII, LayerIII的输入是复数,为了编码统一而强行分的.当然你编码时可以判断当前层是否为1来决定是否分.但是我想每个人最后都会倾向分的.旋转因子tw = cos(2*PI*k/N)-j*sin(2*PI*k/N);也可以分为实部和虚部,令其为tw_r, tw_i;则tw = tw_r - j*tw_i;Xk = (x_rk + j*x_ik) + (tw_rj*tw_i) * (x_rk+N/2+j*x_ik+N/2)则X_Rk = x_rk + tw_r*x_rk+N/2 + tw_i*x_ik+N/2;X_Ik = x_ik - tw_i*x_rk+N/2 + tw_r*x_ik+N/2;LayerII部分,具有2个颗粒,每个颗粒4个输入(注意4个输入的来源,由LayerI友情提供,感谢感谢J)LayerIII部分,具有1个颗粒,每个颗粒8个输入(注意8个输入的来源,由LayerII友情提供,感谢感谢J)LayerI, LayerII, LayerIII从左往右,蝶形信号运算流非常明显!假令输入为xk, xk+N/2,输出为Xk, Xk+N/2. xk分解为x_rk, x_ik部分则该蝶形运算为Xk= (x_rk-j*x_ik) + (x_rk+N/2-j*x_ik+N/2)*(cos(2*PI*k/N)-j*sin(2*PI*k/N);再令cos(2*PI*k/N)为tw1, sin(2*PI*k/N)为tw2则Xk = (x_rk-j*x_ik) + (x_rk+N/2-j*x_ik+N/2)*(tw1-j*tw2);X_Rk = x_rk + x_rk+N/2*tw1 -x_ik+N/2*tw2;X_IK = x_ikx_rk = x_rk + x_rk+b*tw1 + x_ik+b*tw2;x_ik = x_ik - x_rk+b*tw2 + x_ik+b*tw1;譬如8点输入x81.先分割成2部分:x0, x2, x4, x6和x1, x3, x5, x72.信号x0, x2, x4, x6再分割成x0, x4和x2, x6信号x1, x3, x5, x7再分割成x1, x5和x3, x73.无法分割了,已经分割成2点了J.如上图:在LayerI的时候,我们是对2点进行DFT.(一共4次DFT )输入为x0&x4;x2&x6;x1&x5;x3&x7输出为y0,y1;Y2,y3;Y4,y5;Y6,y7;流程:I.希望将输入直接转换为x0, x4, x2, x6, x1, x5, x3, x7的顺序II.对转换顺序后的信号进行4次DFT步骤I代码实现/*反转算法.这个算法效率比较低!先用起来在说,之后需要进行优化.*/staticvoidbitrev(void)intp=1, q, i;intbit_rev N ;floatxx_r N ;bit_rev 0 = 0;while( p N )for(q=0; qp; q+)bit_rev q = bit_rev q * 2;bit_rev q + p = bit_rev q + 1;p *= 2;for(i=0; iN; i+)xx_r i = x_r i ;for(i=0; iN; i+)x_ri = xx_r bit_revi ;/ -此刻序列x重排完毕-步骤II代码实现int j;float TR;/临时变量float tw1;/旋转因子/*两点DFT */for(k=0; kN; k+=2)/两点DFT简化告诉我们tw1=1TR = x_rk;/ TR就是A, x_rk+b就是B.x_rk= TR + tw1*x_rk+b;x_rk+b = TR - tw1*x_rk+b;在LayerII的时候,我们希望得到z,就需要对y进行DFT.y0,y2;y1,y3;y4,y6;y5,y7;z0,z1;z2,z3;z4,z5;z6,z7;在LayerIII的时候,我们希望得到v,就需要对z进行DFT.z0,z4;z1,z5;z2,z6;z3,z7;v0,v1;v2,v3;v4,v5;v6,v7;准备令输入为xs, xs+N/2,输出为ys, ys+N/2这个N绝对不是上面的8,这个N是当前颗粒的输入样本总量对于LayerI而言N是2;对于LayerII而言N是4;对于LayerIII而言N是8复数乘法:(a+j*b) * (c+j*d)实部= a*c bd;虚部= ad + bc;旋转因子:实现(C描述)#include#include#include/#include complex.h/ -#defineN8/64#defineM3/6/2m=N#definePI3.1415926/ -floattwiddleN/2 = 1.0, 0.707, 0.0, -0.707;floatx_rN = 1, 1, 1, 1, 0, 0, 0, 0;floatx_iN;/N=8/*floattwiddleN/2 = 1,0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,0.7075, 0.6349,0.5561,0.4721,0.3835,0.2912,0.1961,0.0991,0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951;/N=64floatx_rN=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,;floatx_iN;*/FILE *fp;/ - func -/*初始化输出虚部*/staticvoidfft_init(void)inti;for(i=0; iN; i+)x_ii = 0.0;/*反转算法.将时域信号重新排序.*这个算法有改进的空间*/staticvoidbitrev(void)intp=1, q, i;intbit_rev N ;/floatxx_r N ;/bit_rev 0 = 0;while( p N )for(q=0; qp; q+)bit_rev q = bit_rev q * 2;bit_rev q + p = bit_rev q + 1;p *= 2;for(i=0; iN; i+)xx_r i = x_r i ;for(i=0; iN; i+)x_ri = xx_r bit_revi ;/* - add by sshc625 - */staticvoidbitrev2(void)return;/*/voiddisplay(void)printf(/n/n);inti;for(i=0; iN; i+)printf(%f/t%f/n, x_ri, x_ii);/*/voidfft1(void)fp = fopen(log1.txt,a+);intL, i, b, j, p, k, tx1, tx2;floatTR, TI, temp;/临时变量floattw1, tw2;/*深M.对层进行循环. L为当前层,总层数为M. */for(L=1; L 0)b *= 2;i-;/ -是否外层对颗粒循环,内层对样本点循环逻辑性更强一些呢! -/* outter对参与DFT的样本点进行循环* L=1,循环了1次(4个颗粒,每个颗粒2个样本点)* L=2,循环了2次(2个颗粒,每个颗粒4个样本点)* L=3,循环了4次(1个颗粒,每个颗粒8个样本点)*/for(j=0; j 0)p = p*2;i-;p= p * j;tx1 = p % N;tx2 = tx1 + 3*N/4;tx2 = tx2 % N;/ tw1是cos部分,实部; tw2是sin部分,虚数部分.tw1 = ( tx1=N/2)? -twiddletx1-N/2: twiddle tx1 ;tw2 = ( tx2=N/2)? -twiddletx2-(N/2) : twiddletx2;/* inner对颗粒进行循环* L=1,循环了4次(4个颗粒,每个颗粒2个输入)* L=2,循环了2次(2个颗粒,每个颗粒4个输入)* L=3,循环了1次(1个颗粒,每个颗粒8个输入)*/for(k=j; kN; k=k+2*b)TR = x_rk;/ TR就是A, x_rk+b就是B.TI = x_ik;temp = x_rk+b;/*如果复习一下(a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么*就能理解为什么会如下运算了,只有在L=1时候输入才是实数,之后层的*输入都是复数,为了让所有的层的输入都是复数,我们只好让L=1时候的*输入虚部为0* x_ik+b*tw2是两个虚数相乘*/fprintf(fp,tw1=%f, tw2=%f/n, tw1, tw2);x_rk= TR + x_rk+b*tw1 + x_ik+b*tw2;x_ik= TI - x_rk+b*tw2 + x_ik+b*tw1;x_rk+b = TR - x_rk+b*tw1 - x_ik+b*tw2;x_ik+b = TI + temp*tw2- x_ik+b*tw1;fprintf(fp,k=%d, x_rk=%f, x_ik=%f/n, k, x_rk, x_ik);fprintf(fp,k=%d, x_rk=%f, x_ik=%f/n, k+b, x_rk+b, x_ik+b);/* - add by sshc625 -*该实现的流程为* for( Layer )*for( Granule )*for( Sample )*/voidfft2(void)fp = fopen(log2.txt,a+);intcur_layer, gr_num, i, k, p;floattmp_real, tmp_imag, temp;/临时变量,记录实部floattw1, tw2;/旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.intstep;/步进intsample_num;/颗粒的样本总数(各层不同,因为各层颗粒的输入不同)/*对层循环*/for(cur_layer=1; cur_layer 0)i-;gr_num *= 2;/*每个颗粒的输入样本数N */sample_num= (int)pow(2, cur_layer);/*步进.步进是N/2 */step= sample_num/2;/*/k = 0;/*对颗粒进行循环*/for(i=0; igr_num; i+)/*对样本点进行循环,注意上限和步进*/for(p=0; psample_num/2; p+)/旋转因子,需要优化.tw1 = cos(2*PI*p/pow(2, cur_layer);tw2 = -sin(2*PI*p/pow(2, cur_layer);tmp_real = x_rk+p;tmp_imag = x_ik+p;temp = x_rk+p+step;/*(tw1+jtw2)(x_rk+jx_ik)* real : tw1*x_rk - tw2*x_ik* imag : tw1*x_ik + tw2*x_rk*我想不抽象出一个* typedef struct *double real;/实部*double imag;/虚部* complex;以及针对complex的操作*来简化复数运算是否是因为效率上的考虑!*/*蝶形算法*/x_rk+p= tmp_real + ( tw1*x_rk+p+step - tw2*x_ik+p+step );x_ik+p= tmp_imag + ( tw2*x_rk+p+step + tw1*x_ik+p+step );/* Xk = A(k)+WB(k)* Xk+N/2 = A(k)-WB(k)的性质可以优化这里*/旋转因子,需要优化.tw1 = cos(2*PI*(p+step)/pow(2, cur_layer);tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer);x_rk+p+step= tmp_real + ( tw1*temp - tw2*x_ik+p+step );x_ik+p+step= tmp_imag + ( tw2*temp + tw1*x_ik+p+step );printf(k=%d, x_rk=%f, x_ik=%f/n, k+p, x_rk+p, x_ik+p);printf(k=%d, x_rk=%f, x_ik=%f/n, k+p+step, x_rk+p+step, x_ik+p+step);/*开跳!:) */k += 2*step;/*后记:*究竟是颗粒在外层循环还是样本输入在外层,好象也差不多,复杂度完全一样.*但以我资质愚钝花费了不少时间才弄明白这数十行代码.*从中我发现一个于我非常有帮助的教训,很久以前我写过一部分算法,其中绝大多数都是递归.*将数据量减少,减少再减少,用归纳的方式来找出数据量加大代码的规律*比如FFT* 1.先写死LayerI的代码;然后再把LayerI的输出作为LayerII的输入,又写死代码; .*大约3层就可以统计出规律来.这和递归也是一样,先写死一两层,自然就出来了!* 2.有的功能可以写伪代码,不急于求出结果,降低复杂性,把逻辑结果定出来后再添加.*比如旋转因子就可以写死,就写1.0.流程出来后再写旋转因子.*寥寥数语,我可真是流了不少汗! Happy!*/voiddft(void)inti, n, k, tx1, tx2;floattw1,tw2;floatxx_rN,xx_iN;/* clear any data in Real and Imaginary result arrays prior to DFT*/for(k=0; k=N-1; k+)xx_rk = xx_ik = x_ik = 0.0;/ caculate the DFTfor(k=0; k=(N-1); k+)for(n=0; n= (N/2)tw1 = -twiddletx1-(N/2);elsetw1 = twiddletx1;if(tx2 = (N/2)tw2 = -twiddletx2-(N/2);elsetw2 = twiddletx2;xx_rk = xx_rk+x_rn*tw1;xx_ik = xx_ik+x_rn*tw2;xx_ik = -xx_ik;/ displayfor(i=0; iN; i+)printf(%f/t%f/n, xx_ri, xx_ii);/ -intmain(void)fft_init( );bitrev( );/ bitrev2( );/fft1( );fft2( );display( );system(pause);/ dft();return1;#include #include /* 快速福利叶变换C函数函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依 赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复 数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的 复数使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的 应该为2的N次方，不满足此条件时应在后面补0函数调用：FFT(s);时 间：2010-2-20版 本：Ver1.0参考文献： */#include#define PI 3.1415926535897932384626433832795028841971 /定义圆周率值#define FFT_N 128 /定义福利叶变换的点数struct compx float real,imag; /定义一个复数结构struct compx sFFT_N; /FFT输入和输出：从S1开始存放，根据大小自己定义/*函数原型：struct compx EE(struct compx b1,struct compx b2) 函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b输出参数：a和b的乘积，以联合体的形式输出*/struct compx EE(struct compx a,struct compx b) struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c);/*函数原型：void FFT(struct compx *xin,int N)函数功能：对输入的复数组进行快速傅里叶变