C语言、Matlab实现FFT几种编程实例..

1、C语言、MATLAB实现FFT几种方法总结前人经验，仅供参考/ 一、/c语言程序 / / #in elude <iom128.h> #i nclude vintrin sics.h> #in clude<math.h>#defi ne PI 3.1415926535897932384626433832795028841971/ 定义圆周率值#define FFT N 128/定义福利叶变换的点数struct compx float real,imag;/ 定义一个复数结构struct compx sFFT_N; /FFT输入和输出：从S1开始存放，根据大小自己定义

2、 /*m*mmm*mmm*mmstruct compx EE(struct compx b1,struct compx b2)对两个复数进行乘法运算两个以联合体定义的复数a,ba和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 *xi n,i

3、 nt N)对输入的复数组进行快速傅里叶变换(FFT)*xin复数结构体组的首地址指针，struct型函数原型函数功能输入参数*/ void FFT(struct compx *xi n)int f,m, nv2, nm1,i,k,l,j=0;struct compx u,w,t;nv2=FFT_N/2; n m仁FFT_N-1; for(i=0;i vn m1;i+) if(i<j)t=xi nj;xi nj=x in i; xin i=t;k=nv2;/变址运算，即把自然顺序变成倒位序，采用雷德算法/如果ivj,即进行变址while(k<=j) j=j-k; k=k/2;某个位

4、为0/求j的下一个倒位序/如果k<=j,表示j的最高位为1/把最高位变成0/k/2，比较次高位，依次类推，逐个比较，直到j=j+k；/把0改为1int le,lei,i p;f=FFT_N;for(l=1;(f=f/2)!=1;l+)/FFT运算核，使用蝶形运算完成FFT运算/计算I的值，即计算蝶形级数for(m=1;mv=l;m+) le=2<<(m-1);lei=le/2;u.real=1.0;u.imag=0.0;/控制蝶形结级数/m表示第m级蝶形，l为蝶形级总数l=log（ 2）N /le蝶形结距离，即第m级蝶形的蝶形结相距le点/同一蝶形结中参加运算的两点的距离/u

5、为蝶形结运算系数，初始值为1/w为系数商，即当前系数与前一个系数的商w.real=cos (P l/lei);w.imag=-si n(P l/lei);for(j=0;j<=lei-1;j+)for(i=j；iv=FFT_N-1;i=i+le) /控制同一蝶形结运算，即计算系数相同蝶形结ip=i+lei;t=EE(xi nip ,u);/控制计算不同种蝶形结，即计算系数不同的蝶形结/i，ip分别表示参加蝶形运算的两个节点/蝶形运算，详见公式xi ni p.real=x in i.real-t.real;xi ni p.imag=x in i.imag-t.imag;xin i.real

6、=x in i.real+t.real;xi ni.imag=xi n i.imag+t.imag;u=EE(u,w);/改变系数，进行下一个蝶形运算/*函数原型：void mai n()函数功能：测试FFT变换，演示函数使用方法 输入参数：无输出参数：无*/void mai n()int i;for(i=0;i<FFT_N;i+)/给结构体赋值样,si.real=sin(2*3.141592653589793*i/FFT_N); /实部为正弦波 FFT_N 点采赋值为1si.imag=O;虚部为0FFT(s);/进行快速福利叶变换for(i=0;i<FFT_N;i+)/求变换后结

7、果的模值，存入复数的实部部分si.real=sqrt(si.real*si.real+si.imag*si.imag);while(1);%/ 二%/ %/%/MATLAB仿真信号源的源程序:Clear;CIc;t=0:0.01:3; yl=100*si n(p i/3*t); n=l;for t=-O:O.01:10;y2(1, n)=-61.1614*ex p(-0.9*t);n=n+;endmin( y2)y=yl, y2;figure(1);plot(y); %fun box in g(O.001+1/3)%/%/快速傅里叶变换 matlab程序: %/clc;clear;clf;N=

8、inpu t('Node nu mber')T=inpu t('cai yang jia n ge')f=inpu t('fre nquen cy')choise=inpu t('add zero or not? n=0:T:(N-1)*T;k=0:N-1;1/0 ')%采样点x=si n(2*pi*f*n);if choise=1e=inpu t('I nput the nu mber of zeros!')x=x zeros(1,e)N=N+e;else%加零 k=0:N-1;%给k重新赋值，因为有可能出现加零

9、状况bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)+1;% 利用库函数进行变址运算 for l=1:NendX(l)=x(bianzhi(l);%将采样后的值按照变址运算后的顺序放入X矩阵中endd1=1;for m=1:log2(N) d2=d1; d1=d1*2;W=1;dw=ex p(-j*p i/d2);for t=1:d2%做蝶形运算的两个数之间的距离%同一级之下蝶形结之间的距离%蝶形运算系数的初始值%蝶形运算系数的增加量for i=t:d1:Ni1=i+d2;if(i1>N)break; else%判断是否超出范围p=X(i1)*

10、W;X(i1)=X(i)-p; X(i)=X(i)+p; endendW=W*dw;%蝶形运算%蝶形运算系数的变化endendsub plot(2,2,1);t=0:0.0000001:N*T;%画原曲线p lot(t,si n(2* pi*f*t);sub plot(2,2,2);%画采样后的离散信号stem(k,x);sub plot(2,2,3);%画自己的fft的运算结果stem(k,abs(X)/max(abs(X);sub plot(2,2,4);stem(k,abs(fft(x)/max(abs(fft(x); %调用系统的函数绘制结果，与自己的结果进 行比较/ 三/ /FFT

11、的 C 语言算法实现/ 程序如下：/*Fft*/#in eludevstdio.h>#in elude<math.h>#in elude<stdlib.h>#defi neN 1000structtyp edefreal;doubledoubleimg;eo mpl ex;voidfft(); /*快速傅里叶变换*/ifft(); /*快速傅里叶逆变换voidin itW();voideha nge();voidadd(eo mp lex,eo mplex,eo mp lex*);/*复数加法*/voidmul(eo mp lex,eo mp lex,eo mpl

12、ex *);/*复数乘法*/voidsub(eo mp lex,eo mp lex,eom plex *);/*复数减法*/void*/voiddivi(complex,complex,complex*);/* 复数除法 */void output(); /* 输出结果 */comp lex xN, *W;/* 输出序列的值 */ int size_x=0;/*输入序列的长度，只限2 double PI;int mai n()int i,method;的N次方*/system("cls");P l=ata n( 1)*4;prin tf(" PI easeinpu

13、t the size of/*输入序列的长度*/x:n");scan f("%d"，&size_x);prin tf(" PI easeinput the data in/*输入序列对应的值*/xN:(such as:5 6)n");for(i=0;ivsize_x;i+)scan f("%lf %lf", &xi.real, &xi.img);ini two；/*选择FFT或逆FFT运算*/prin tf("UseFFT(0) or IFFT(1)?n");scan f(&quo

14、t;%d",&method);if(method=0)fft();elseifft();out put();return 0;/*进行基-2 FFT运算*/void fft()int i=0,j=0,k=0,l=0; complexup, dow n,p roduct;cha nge();for(i=0;i< log(size_x)/log(2) ;i+)l=1<<i;2*1for(j=0;jvsize_x;j+=for(k=0;kv|;k+)mul(x|j+k+l,Wsize_x*k/2/l,&p roduct); add(xj+k ,p rodu

15、ct，&up);sub(xj+k, product，&dow n);xj+k=up; xj+k+l=dow n;voidifftO int i=0,j=0,k=0,l=size_x;complex up, dow n;for(i=0;i< (int)( log(size_x)/log(2) );i+) /* 蝶形运算 */ l/=2;2*l)for(j=0;j<size_x;j+=for(k=0;k<l;k+)add(xj+k , xj+k+l,&up);up. real/=2; up .img/=2;sub(xj+k,xj+k+l,&dow

16、n);dow n.real/=2;dow n.img/=2;divi(dow n,Wsize_x*k/2/l, &dow n);xj+k=up;x|j+k+l=dow n;cha nge();voidin itW() int i;size_x);W=(co mp lex*)malloc(sizeof(co mp lex)for(i=0;i<size_x;i+)Wi.real=cos(2* Pl/size_x*i);Wi.img=-1*si n(2* Pl/size_x*i); complex un sig ned doublevoidcha nge()temp;short i=O

17、,j=O,k=O; 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=xi;xi=xj;xj=te mp;voidoutput()/* 输出结果 */ int i;follows' n");prin tf("The result are asfor(i=0;i<size_x;i+)prin tf("%.4f",xi.real);if(xi.

18、img>=0.0001)prin tf("+%.4fjn",xi.img);else if(fabs(xi.img)<0.0001)prin tf("n");elseprin tf("%.4fjn",xi.img);voidadd(co mp lexa,co mp lexb,co mp lex*c) c->real=a.real+b.real; c->img=a.img+b.img;voidmul(co mp lexa,co mplexb,co mplex*c)c->real=a.real*b.realc

19、->img=a.real*b.imgvoidsub(co mp lex- a.img*b.img;+ a.img*b.real;a,co mp lexb,co mp lex*c)c->real=a.real-b.real;c->img=a.img-b.img;voiddivi(co mplexa,co mplexb,co mplex*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);/ 四、/ / %/ 选带傅里叶变换(zoom-fft) MATLAB%移频（将选带的中心频率移动到零频）%数字低通滤波器（防止频率混叠）%重新采样（将采样的数据再次间隔采样，间隔的数据取决于分析的带宽，就是放大倍数）%复FFT （由于经过了移频，所以数据不是实数了）%频率调整（将负半轴的频率成分移