C语言、Matlab实现FFT几种编程实例(10页).doc

《C语言、Matlab实现FFT几种编程实例(10页).doc》由会员分享，可在线阅读，更多相关《C语言、Matlab实现FFT几种编程实例(10页).doc（10页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、-C语言、Matlab实现FFT几种编程实例-第 10 页C语言、MATLAB实现FFT几种方法总结前人经验，仅供参考/一、/c语言程序/#include #include #include#define PI 3.1415926535897932384626433832795028841971 /定义圆周率值#define FFT_N 128 /定义福利叶变换的点数struct compx float real,imag; /定义一个复数结构struct compx sFFT_N; /FFT输入和输出：从S1开始存放，根据大小自己定义函数原型：struct compx EE(struct c

2、ompx 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)函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数：*xin复

3、数结构体组的首地址指针，struct型void FFT(struct compx *xin) int f,m,nv2,nm1,i,k,l,j=0; struct compx u,w,t; nv2=FFT_N/2; /变址运算，即把自然顺序变成倒位序，采用雷德算法 nm1=FFT_N-1; for(i=0;inm1;i+) if(ij) /如果ij,即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求j的下一个倒位序 while(k=j) /如果k=j,表示j的最高位为1 j=j-k; /把最高位变成0 k=k/2; /k/2，比较次高位，依次类推，逐个比较，直到

4、某个位为0 j=j+k; /把0改为1 int le,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）N le=2(m-1); /le蝶形结距离，即第m级蝶形的蝶形结相距le点 lei=le/2; /同一蝶形结中参加运算的两点的距离 u.real=1.0; /u为蝶形结运算系数，初始值为1 u.imag=0.0; w.real=cos(PI/lei); /w为系数商，即当前系数与前一个系数

5、的商 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(xinip,u); /蝶形运算，详见公式 xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag; u=EE(u,w); /改

6、变系数，进行下一个蝶形运算函数原型：void main() 函数功能：测试FFT变换，演示函数使用方法输入参数：无输出参数：无void main() int i; for(i=0;iFFT_N;i+) /给结构体赋值 si.real=sin(2*3.141592653589793*i/FFT_N); /实部为正弦波FFT_N点采样，赋值为1 si.imag=0; /虚部为0 FFT(s); /进行快速福利叶变换 for(i=0;iN)break; %判断是否超出范围 else p=X(i1)*W; X(i1)=X(i)-p; X(i)=X(i)+p; %蝶形运算 end end W=W*dw;

7、 %蝶形运算系数的变化 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*/ #include #include #include #de

8、fine 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(comple

9、x ,complex ,complex *);/*复数除法*/ void output(); /*输出结果*/ complex xN, *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度，只限2的N次方*/ double PI; int main() int i,method; system(cls); PI=atan(1)*4; printf(Please input the size of x:n); /*输入序列的长度*/ scanf(%d,&size_x); printf(Please input the data in xN:(such as:5 6)n); /

10、*输入序列对应的值*/ for(i=0;isize_x;i+) scanf(%lf %lf,&xi.real,&xi.img); initW(); /*选择FFT或逆FFT运算*/ printf(Use FFT(0) or IFFT(1)?n); scanf(%d,&method); if(method=0) fft(); else ifft(); output(); return 0; /*进行基-2 FFT运算*/ void fft() int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i log(size_x)

11、/log(2) ;i+) l=1i; for(j=0;jsize_x;j+= 2*l ) for(k=0;kl;k+) mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; void ifft() int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i (int)( log(size_x)/log(2) );i+) /*蝶形运算*/ l/=2; for(j=0;jsize_x;j+= 2*

12、l ) for(k=0;kl;k+) add(xj+k,xj+k+l,&up); up.real/=2;up.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2; divi(down,Wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change(); void initW() int i; W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;isize_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-

13、1*sin(2*PI/size_x*i); void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i0 ) j=j1; if(ji) temp=xi; xi=xj; xj=temp; void output() /*输出结果*/ int i; printf(The result are as followsn); for(i=0;i=0.0001) printf(+%.4fjn,xi.img); else if(fabs(xi.img)real=a.real+b.real; c-img=a.img+b

14、.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

15、*b.img); c-img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); /四、%/选带傅里叶变换(zoom-fft) MATLAB%移频 （将选带的中心频率移动到零频）%数字低通滤波器 （防止频率混叠）%重新采样 （将采样的数据再次间隔采样，间隔的数据取决于分析的带宽，就是放大倍数）%复FFT （由于经过了移频，所以数据不是实数了）%频率调整 （将负半轴的频率成分移到正半轴）functionf,y=zfft(x,fi,fa,fs)%x为采集的数据%fi为分析的起始频率%fa为分析的截止频率%fs为采集数据的采样频率%f为

16、输出的频率序列%y为输出的幅值序列(实数）f0=(fi+fa)/2;%中心频率N=length(x);%数据长度r=0:N-1;b=2*pi*f0.*r./fs;x1=x.*exp(-1j.*b);%移频bw=fa-fi;B=fir1(32,bw/fs);%滤波截止频率为0.5bwx2=filter(B,1,x1);c=x2(1:floor(fs/bw):N);%重新采样N1=length(c);f=linspace(fi,fa,N1);y=abs(fft(c)./N1*2;y=circshift(y,0,floor(N1/2);%将负半轴的幅值移过来end%上边 四程序 测试应用实例：fs=2048;T=100;t=0:1/fs:T;x=30*cos(2*pi*110.*t)+30*cos(2*pi*111.45.*t)+25*cos(2*pi*112.3*t)+48*cos(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t);f,y=zfft(x,109,115,fs);plot(f,y);图片效果：