FFT运算C语言基2蝶形运算程序

#include"reg51.h" #include"math.h" #include"fft.h" volatile struct FFT_DATA FFTDATA[Sample_Num]; volatile struct Result result[Sample_Num/2]; volatile struct WN Wn[Sample_Num]; /********************************************************************** *函数：void FFT_N(void)* *功能：FFT运算，基2蝶形运算 * *参数描述：xx是倒序操作后的新的数组的编号，Num是根据Sample_Num计算得到的数据位数， **********比如说Samle_Num=8，则Num=3* * **********************************************************************/ void FFT_N(void) { Uint16 xx,Num; Uint16 i,j,k,b,p,L; float TR,TI,temp; Num = (Uint16)(log(Sample_Num)/log(2)+0.0001); ///=======following code invert sequence 倒序+原位*********** for ( i=0;i<Sample_Num;i++ ) { xx = 0; for(j=0;j<Num;j++) { if(i&(1<<j)) xx|=(1<<(Num-1-j)); } FFTDATA[xx].dataI=FFTDATA[i].dataR;//先把倒序后的数据存放到虚部 } for ( i=0;i<Sample_Num;i++ ) { FFTDATA[i].dataR=FFTDATA[i].dataI; //原位操作，放到原来存储位置 FFTDATA[i].dataI=0; } //=************** following code FFT ******************* //L是第L级蝶形运算，与Num的值有关，总共是Num级的FFT运算 //b是第L级做蝶形运算的两个点的距离，p是蝶形运算的旋转因子Wnp,根据旋转因子的对称性，知道旋转因子只有N/2个值 //由于第K个点和k+b个点进行蝶形运算结果还得保存在第K个点的存储位置，所以要定义中间变量TR,TI,temp防止上次结果被覆盖 //FFT算法流图旋转因子规律 //第一级的蝶形系数p均为0，蝶形节点的距离为1 //第二级的蝶形系数p为0,N/4，蝶形节点的距离为2。 //第三级的蝶形系数p为0,N/8,2N/8,3N/8，蝶形节点的距离为4。 //第M级的蝶形系数p为0,1,...N/2-1，蝶形节点的距离为N/2。 for ( L=1;L<=Num;L++ ) { b = pow(2,L-1);//计算每一级蝶形运算的距离 for ( j=0;j<=b-1;j++ ) //0到b-1是p的取值范围，根据b的值，求蝶形系数p的范围 { p = pow(2,Num-L)*j;//计算每一级蝶形运算的旋转因子 for ( k=j;k<Sample_Num;k=k+2*b )//由k和k+b距离的两个点进行蝶形运算，要注意原位计算 { TR=FFTDATA[k].dataR; TI=FFTDATA[k].dataI; temp=FFTDATA[k+b].dataR; FFTDATA[k].dataR=FFTDATA[k].dataR+FFTDATA[k+b].dataR*Wn[p].cos+FFTDATA[k+b].dataI*Wn[p].sin; FFTDATA[k].dataI=FFTDATA[k].dataI-FFTDATA[k+b].dataR*Wn[p].sin+FFTDATA[k+b].dataI*Wn[p].cos; FFTDATA[k+b].dataR=TR-FFTDATA[k+b].dataR*Wn[p].cos-FFTDATA[k+b].dataI*Wn[p].sin; FFTDATA[k+b].dataI=TI+temp*Wn[p].sin-FFTDATA[k+b].dataI*Wn[p].cos; } } } for ( i=0;i<Sample_Num/2;i++ ) { result[i].freq = 10 * i; if(i==0) result[0].U = (sqrt(FFTDATA[0].dataR*FFTDATA[0].dataR+FFTDATA[0].dataI*FFTDATA[0].dataI)/Sample_Num+0.001);//直流分量结果 else result[i].U = 1.414*(sqrt(FFTDATA[i].dataR*FFTDATA[i].dataR+FFTDATA[i].dataI*FFTDATA[i].dataI)/Sample_Num+0.001); //各次谐波有效值 result[i].angle = atan(FFTDATA[i].dataI/FFTDATA[i].dataR); result[i].P = result[i].U*result[i].U/50; //U*U/R; //各次谐波的功率 } } //初始化正余弦数表，旋转因子的计算通过查表计算 void InitForFFT() { int i; for ( i=0;i<Sample_Num;i++ ) { Wn[i].sin=sin(PI*2*i/Sample_Num); Wn[i].cos=cos(PI*2*i/Sample_Num); } } void MakeWave() { int i; for ( i=0;i<Sample_Num;i++ ) { FFTDATA[i].dataR =20+30*1.414*cos(2*PI*i/Sample_Num-PI*30/180)+10*1.414*cos(2*PI*3*i/Sample_Num+PI*90/180); } } void main(void) { InitForFFT(); MakeWave(); FFT_N(); while(1); }