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

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

2、truct 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)函数功能：对输入的

3、复数组进行快速傅里叶变换（FFT）输入参数：*xin复数结构体组的首地址指针，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;i<nm1;i+) if(i<j) /如果i<j,即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求j的下一个倒位序 while(k<=j) /如果k<=j,表示j的最高位

4、为1 j=j-k; /把最高位变成0 k=k/2; /k/2，比较次高位，依次类推，逐个比较，直到某个位为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为蝶

5、形结运算系数，初始值为1 u.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(xinip,u); /蝶形运算，详见公式 xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag

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

7、;FFT_N;i+) /求变换后结果的模值，存入复数的实部部分 si.real=sqrt(si.real*si.real+si.imag*si.imag); while(1); %/二、%/%/%/MATLAB仿真信号源的源程序:Clear;Clc;t=O:O.01:3;yl=100*sin(pi/3*t);n=l;for t=-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;

8、clf;N=input('Node number')T=input('cai yang jian ge')f=input('frenquency')choise=input('add zero or not? 1/0 ')n=0:T:(N-1)*T; %采样点k=0:N-1;x=sin(2*pi*f*n);if choise=1 e=input('Input the number of zeros!') x=x zeros(1,e) N=N+e;elseend %加零k=0:N-1; %给k重新赋值，因为有可能出现

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

10、se p=X(i1)*W; X(i1)=X(i)-p; X(i)=X(i)+p; %蝶形运算 end end W=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); %调用系统的函数绘制结果，与自己的结果进行

11、比较/三、/FFT的C语言算法实现/程序如下： /*FFT*/ #include <stdio.h> #include <math.h> #include <stdlib.h> #define N 1000 typedef struct double real; double img; complex; void fft(); /*快速傅里叶变换*/ void ifft(); /*快速傅里叶逆变换*/ void initW(); void change(); void add(complex ,complex ,complex *); /*复数加法*/ vo

12、id mul(complex ,complex ,complex *); /*复数乘法*/ void sub(complex ,complex ,complex *); /*复数减法*/ void divi(complex ,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;

13、printf("Please input the size of x:n"); /*输入序列的长度*/ scanf("%d",&size_x); printf("Please input the data in xN:(such as:5 6)n"); /*输入序列对应的值*/ for(i=0;i<size_x;i+) scanf("%lf %lf",&xi.real,&xi.img); initW(); /*选择FFT或逆FFT运算*/ printf("Use FFT(0)

14、 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)/log(2) ;i+) l=1<<i; for(j=0;j<size_x;j+= 2*l ) for(k=0;k<l;k+) mul(xj

15、+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;j<size_x;j+= 2*l ) for(k=0;k<l;k+) add(xj+k,xj+k+l,&up

16、); 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;i<size_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-1*sin(2*PI/size_x*i); void

17、 change() complex temp; unsigned short i=0,j=0,k=0; double 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=temp; void output() /*输出结果*/ int i; printf("The result are as followsn"); for(i

18、=0;i<size_x;i+) printf("%.4f",xi.real); if(xi.img>=0.0001) printf("+%.4fjn",xi.img); else if(fabs(xi.img)<0.0001) printf("n"); else printf("%.4fjn",xi.img); void add(complex a,complex b,complex *c) c->real=a.real+b.real; c->img=a.img+b.img; void

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

20、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 （由于经过了移频，所以数据不是实数了）%频率调整 （将负半轴的频率成分移到正半轴）function f, y = zfft(x,

21、0;fi, fa, fs)  % x为采集的数据  % fi为分析的起始频率  % fa为分析的截止频率  % fs为采集数据的采样频率  % f为输出的频率序列  % y为输出的幅值序列(实数）    f0 = (fi + fa) / 2;     

22、;         %中心频率  N = length(x);                 %数据长度    r = 0:N-1;  b = 2*pi*f0.*r ./ fs;

23、                 x1 = x .* exp(-1j .* b);          %移频    bw = fa - fi;    

24、                                                  

25、           B = fir1(32, bw / fs);             %滤波 截止频率为0.5bw  x2 = filter(B, 1, x1);                   c = x2(1:floor(fs/bw):N);           %重新采样  N1 = length(c);&#