搜索
C语言、MATLAB实现FFT几种方法总结前人经验，仅供参考///一、/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////c语言程序/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////#includeiom128.h#includeintrinsics.h#includemath.h#definePI3.1415926535897932384626433832795028841971//定义圆周率值#defineFFT_N128//定义福利叶变换的点数structcompx{floatreal,imag;};//定义一个复数结构structcompxs[FFT_N];//FFT输入和输出：从S[1]开始存放，根据大小自己定义/*******************************************************************函数原型：structcompxEE(structcompxb1,structcompxb2)函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b输出参数：a和b的乘积，以联合体的形式输出*******************************************************************/structcompxEE(structcompxa,structcompxb){structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/*****************************************************************函数原型：voidFFT(structcompx*xin,intN)函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数：*xin复数结构体组的首地址指针，struct型*****************************************************************/voidFFT(structcompx*xin){intf,m,nv2,nm1,i,k,l,j=0;structcompxu,w,t;nv2=FFT_N/2;//变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1;for(i=0;inm1;i++){if(ij)//如果ij,即进行变址{t=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2;//求j的下一个倒位序while(k=j)//如果k=j,表示j的最高位为1{j=j-k;//把最高位变成0k=k/2;//k/2，比较次高位，依次类推，逐个比较，直到某个位为0}j=j+k;//把0改为1}{intle,lei,ip;//FFT运算核，使用蝶形运算完成FFT运算f=FFT_N;for(l=1;(f=f/2)!=1;l++)//计算l的值，即计算蝶形级数;for(m=1;m=l;m++)//控制蝶形结级数{//m表示第m级蝶形，l为蝶形级总数l=log（2）Nle=2(m-1);//le蝶形结距离，即第m级蝶形的蝶形结相距le点lei=le/2;//同一蝶形结中参加运算的两点的距离u.real=1.0;//u为蝶形结运算系数，初始值为1u.imag=0.0;w.real=cos(PI/lei);//w为系数商，即当前系数与前一个系数的商w.imag=-sin(PI/lei);for(j=0;j=lei-1;j++)//控制计算不同种蝶形结，即计算系数不同的蝶形结{for(i=j;i=FFT_N-1;i=i+le)//控制同一蝶形结运算，即计算系数相同蝶形结{ip=i+lei;//i，ip分别表示参加蝶形运算的两个节点t=EE(xin[ip],u);//蝶形运算，详见公式xin[ip].real=xin[i].real-t.real;xin[ip].imag=xin[i].imag-t.imag;xin[i].real=xin[i].real+t.real;xin[i].imag=xin[i].imag+t.imag;}u=EE(u,w);//改变系数，进行下一个蝶形运算}}}}/************************************************************函数原型：voidmain()函数功能：测试FFT变换，演示函数使用方法输入参数：无输出参数：无************************************************************/voidmain(){inti;for(i=0;iFFT_N;i++)//给结构体赋值{s[i].real=sin(2*3.141592653589793*i/FFT_N);//实部为正弦波FFT_N点采样，赋值为1s[i].imag=0;//虚部为0}FFT(s);//进行快速福利叶变换for(i=0;iFFT_N;i++)//求变换后结果的模值，存入复数的实部部分s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);while(1);}%////二、%/////////////////////////////////////////////////////////////////////////////////////////////////////////////%///////////////////////////////////////////////////////////////////////////////////////////////////////////%////////////////////////////////MATLAB仿真信号源的源程序::Clear;Clc;t=O:O.01:3;yl=100*sin(pi/3*t);n=l;fort=-O:O.01:10;y2(1,n)=-61.1614*exp(-0.9*t);n=n+;endmin(y2)y=[yl，y2];figure(1);plot(y);%funboxing(O.001+1/3)%////////////////////////%/////////快速傅里叶变换matlab程序:%////////////////////////clc;clear;clf;N=input('Nodenumber')T=input('caiyangjiange')f=input('frenquency')choise=input('addzeroornot?1/0')n=0:T:(N-1)*T;%采样点k=0:N-1;x=sin(2*pi*f*n);ifchoise==1e=input('Inputthenumberofzeros!')x=[xzeros(1,e)]N=N+e;elseend%加零k=0:N-1;%给k重新赋值，因为有可能出现加零状况bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;%利用库函数进行变址运算forl=1:NX(l)=x(bianzhi(l));%将采样后的值按照变址运算后的顺序放入X矩阵中endd1=1;form=1:log2(N)d2=d1;%做蝶形运算的两个数之间的距离d1=d1*2;%同一级之下蝶形结之间的距离W=1;%蝶形运算系数的初始值dw=exp(-j*pi/d2);%蝶形运算系数的增加量fort=1:d2%fori=t:d1:Ni1=i+d2;if(i1N)break;%判断是否超出范围elsep=X(i1)*W;X(i1)=X(i)-p;X(i)=X(i)+p;%蝶形运算endendW=W*dw;%蝶形运算系数的变化endendsubplot(2,2,1);t=0:0.0000001:N*T;plot(t,sin(2*pi*f*t));%画原曲线subplot(2,2,2);stem(k,x);%画采样后的离散信号subplot(2,2,3);stem(k,abs(X)/max(abs(X)));%画自己的fft的运算结果subplot(2,2,4);stem(k,abs(fft(x))/max(abs(fft(x))));%调用系统的函数绘制结果，与自己的结果进行比较//////三、/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////FFT的C语言算法实现////////////程序如下：/************FFT***********/#includestdio.h#includemath.h#includestdlib.h#defineN1000typedefstruct{doublereal;doubleimg;}complex;voidfft();/*快速傅里叶变换*/voidifft();/*快速傅里叶逆变换*/voidinitW();voidchange();voidadd(complex,complex,complex*);/*复数加法*/voidmul(complex,complex,complex*);/*复数乘法*/voidsub(complex,complex,complex*);/*复数减法*/voiddivi(complex,complex,complex*);/*复数除法*/voidoutput();/*输出结果*/complexx[N],*W;/*输出序列的值*/intsize_x=0;/*输入序列的长度，只限2的N次方*/doublePI;intmain(){inti,method;system(cls);PI=atan(1)*4;printf(Pleaseinputthesizeofx:\n);/*输入序列的长度*/scanf(%d,&size_x);printf(Pleaseinputthedatainx[N]:(suchas:56)\n);/*输入序列对应的值*/for(i=0;isize_x;i++)scanf(%lf%lf,&x[i].real,&x[i].img);initW();/*选择FFT或逆FFT运算*/printf(UseFFT(0)orIFFT(1)?\n);scanf(%d,&method);if(method==0)fft();elseifft();output();return0;}/*进行基-2FFT运算*/voidfft(){inti=0,j=0,k=0,l=0;complexup,down,produ