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);&#