(完整word版)C語言、Matlab實現(xiàn)FFT幾種編程實例..

1、C語言、MATLAB實現(xiàn)FFT幾種方法 總結(jié)前人經(jīng)驗，僅供參考 III 一、 / IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIc語言程序 / / #i nclude #i nclude vintrin sics.h #in clude #defi ne PI 3.1415926535897932384626433832795028841971 / 定義圓周率值 #defi ne FFT_N 128/定義福利葉變換的點數(shù) struct compx float real,imag;/ 定義一個復(fù)數(shù)結(jié)構(gòu) struct compx s

2、FFT_N; /FFT輸入和輸出：從S1開始存放，根據(jù)大小自己定義 I* 函數(shù)原型：struct compx EE(struct compx b1,struct compx b2) 函數(shù)功能：對兩個復(fù)數(shù)進行乘法運算 輸入?yún)?shù)：兩個以聯(lián)合體定義的復(fù)數(shù)a,b 輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出 */ 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); I

3、* 函數(shù)原型：void FFT(struct compx *xi n,int N) 函數(shù)功能：對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT) 輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型 */ void FFT(struct compx *xi n) int f,m, nv2, nm1,i,k,l,j=0; struct compx u,w,t; n v2=FFT_N/2;變址運算，即把自然順序變成倒位序，采用雷德算法 nm仁 FFT_N-1; for(i=0;i n m1;i+) if(ivj) t=xi nj; xi nj=x in i; xi ni=t; k=nv2; whil

4、e(k=j) j=j-k; k=k/2; 某個位為0 j=j+k; /如果ij,即進行變址 /求j的下一個倒位序 /如果k=j,表示j的最咼位為1 /把最高位變成0 /k/2，比較次咼位，依次類推，逐個比較，直到 /把0改為1 int le,lei,ip; f=FFT_N; for(l=1;(f=f/2)!=1;l+) /FFT運算核，使用蝶形運算完成FFT運算 /計算1的值，即計算蝶形級數(shù) J for(m=1;m=l;m+) le=2(m-1); lei=le/2; u.real=1.0; u.imag=0.0; /控制蝶形結(jié)級數(shù) /m表示第m級蝶形，l為蝶形級總數(shù)l=log (2) N /

5、le蝶形結(jié)距離，即第m級蝶形的蝶形結(jié)相距l(xiāng)e點 /同一蝶形結(jié)中參加運算的兩點的距離 /u為蝶形結(jié)運算系數(shù)，初始值為1 w.real=cos(PI/lei);llw為系數(shù)商，即當(dāng)前系數(shù)與前一個系數(shù)的商 w.imag=-si n( PI/lei); for(j=0;jv=lei-1;j+)/控制計算不同種蝶形結(jié)，即計算系數(shù)不同的蝶形結(jié) for(i=j；iv=FFT_N-1;i=i+le) /控制同一蝶形結(jié)運算，即計算系數(shù)相同蝶形結(jié) ip=i+lei;/i，ip分別表示參加蝶形運算的兩個節(jié)點 t=EE(x in ip,u);/蝶形運算，詳見公式 xi nip.real=x in i.real-t.r

6、eal; xi nip.imag=x in i.imag-t.imag; xin i.real=x in i.real+t.real; xi ni.imag=x in i.imag+t.imag; u=EE(u,w);/改變系數(shù)，進行下一個蝶形運算 /* 函數(shù)原型：void mai n() 函數(shù)功能：測試FFT變換，演示函數(shù)使用方法 輸入?yún)?shù)：無 輸出參數(shù)：無 */ void mai n() int i; for(i=0;iFFT_N;i+)給結(jié)構(gòu)體賦值 si.real=sin(2*3.141592653589793*i/FFT_N); /實部為正弦波 FFT_N 點采 樣，賦值為1 si.i

7、mag=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; %做蝶形運算的兩個數(shù)之間的距離 %同一級之下蝶形結(jié)之間的距離 %蝶形運算系數(shù)的初始值 %蝶形運算系數(shù)的增加量 % %判斷是否超出范圍 %蝶形運算 %蝶形運算系數(shù)的變化 endend subplot(2,2,1); t=0:0.0000001:N*T; plot(t,si n(2*pi*f*t); subplot(2,2,2); stem(k,x); subplot(2,2,3);

8、stem(k,abs(X)/max(abs(X); subplot(2,2,4); %畫原曲線 %畫采樣后的離散信號 %畫自己的fft的運算結(jié)果 stem(k,abs(fft(x)/max(abs(fft(x); %調(diào)用系統(tǒng)的函數(shù)繪制結(jié)果，與自己的結(jié)果進 行比較 / 三、 / / /FFT 的 C 語言算法實現(xiàn) / 程序如下： /* CCT* */ #i nclude #in elude #in elude #defi neN 1000 typedef struct double real; double img; complex; void fft(); /*快速傅里葉變換*/ void i

9、fft(); /*快速傅里葉逆變換 */ void in itW(); void cha nge(); void add(complex ,complex ,complex *) ;/*復(fù)數(shù)加法*/ void mul(complex ,complex ,complex * );/*復(fù)數(shù)乘法*/ void sub(complex ,complex ,complex *) ;/*復(fù)數(shù)減法*/ void divi(complex complex complex *);/* 復(fù)數(shù)除法 */ void output(); /* 輸出結(jié)果 */ complex xN,*W;/* 輸出序列的值 */ int

10、size_x=0;/*輸入序列的長度，只限2的N次方*/ double PI; int mai n() int i,method; system(cls); PI=ata n( 1)*4; prin tf(Pleasein put the size of x:n); /*輸入序列的長度*/ scan f(%d， prin tf(Pleasein put the data in xN:(such as:5 6)n); /*輸入序列對應(yīng)的值*/ for(i=0;isize_x;i+) scan f(%lf %lf, ini tW(); /*選擇FFT或逆FFT運算*/ prin tf(UseFFT

11、(0) or IFFT(1)?n); scan f(%d, if(method=0) fft(); else ifft(); output(); return 0; /*進行基-2 FFT運算*/ void fft() int i=0,j=0,k=0,l=0; complex up,dow n,product; cha nge(); for(i=0;i log(size_x)/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, add(xj+k,product， sub(xj+k

12、,product， xj+k=up; xj+k+l=dow n; void ifft() int i=0,j=0,k=0,l=size_x; );i+) /*蝶形運算*/ complexup,dow n; for(i=0;i (int)(log(size_x)/log (2) l/=2; for(j=0;jsize_x;j+= 2*l) for(k=0;kl;k+) add(xj+k , xj+k+l, up.real/=2;up.img/=2; sub(xj+k,xj+k+l, dow n.real/=2;dow n.img/=2; divi(dow n,Wsize_x*k/2/l, xj+

13、k=up; xj+k+l=dow n; cha nge(); void ini tW() size_x); int i; W=(complex *)malloc(sizeof(complex) for(i=0;isize_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-1*si n(2*PI/size_x*i); voidcha nge() complex temp; un sig nedshort i=0,j=0,k=0; double t; for(i=0;i0) j=j1; if(ji) temp=xi; xi=xj; xj=temp; void ou

14、tput()/* 輸出結(jié)果 */ int i; prin tf(The result are as followsn); for(i=0;i=0.0001) prin tf(+%.4fjn,xi.img); else if(fabs(xi.img)real=a.real+b.real; c-img=a.img+b.img; void mul(complex c-real=a.real*b.real c-img=a.real*b.img void sub(complex c-real=a.real-b.real; c-img=a.img-b.img; a,complex b,complex a,

15、complex b,complex - a.img*b.img; + a.img*b.real; a,complex b,complex *c) *c) *c) void divi(complex a,complex b,complex *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 %移頻(將選帶的中

16、心頻率移動到零頻) %數(shù)字低通濾波器(防止頻率混疊) %重新采樣(將采樣的數(shù)據(jù)再次間隔采樣，間隔的數(shù)據(jù)取決于分析的帶寬，就 是放大倍數(shù)) %復(fù)FFT (由于經(jīng)過了移頻，所以數(shù)據(jù)不是實數(shù)了) %頻率調(diào)整(將負(fù)半軸的頻率成分移到正半軸) fun cti on f, y = zfft(x, fi, fa, fs) %x為采集的數(shù)據(jù) % fi為分析的起始頻率 % fa為分析的截止頻率 % fs為采集數(shù)據(jù)的采樣頻率 % f為輸出的頻率序列 % y為輸出的幅值序列(實數(shù)) fO = (fi + fa) / 2;% 中心頻率 N = len gth(x);%數(shù)據(jù)長度 r = 0:N-1; b = 2*pi*f0.*r ./ fs; x1 = x * exp(-1j * b);% 移頻 bw = fa - fi; B = fir1(32, bw / fs);% 濾波 截止頻率為 0.5bw x2 = filter(B, 1, x1); c = x2(1:floor(fs/bw):N); % 重新采樣 N1 = len gth(c); f = lin space(fi, fa, N1); y = abs(fft(c) ./ N1 * 2; y = circs