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

1、C語言、MATLAB 實現(xiàn)FFT幾種方法總結(jié)前人經(jīng)驗，僅供參考/一、/c語言程序/#include <iom128.h>#include <intrinsics.h>#include<math.h>#define PI 3.1415926535897932384626433832795028841971/定義圓周率值#define FFT_N 128/定義福利葉變換的點數(shù)struct compx float real,imag;定義一個復數(shù)結(jié)構(gòu)struct compx sFFT_N; FFT輸入和輸出：從S1開始存放，根據(jù)大小自己定義I*函數(shù)原型：struc

2、t compx EE(struct compx b1,struct compx b2)函數(shù)功能：對兩個復數(shù)進行乘法運算輸入?yún)?shù)：兩個以聯(lián)合體定義的復數(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);/*函數(shù)原型：void FFT(struct compx *xin,int N)函數(shù)功能：對輸入的復數(shù)組進行快速傅里

3、葉變換(FFT)輸入?yún)?shù)：*xin復數(shù)結(jié)構(gòu)體組的首地址指針，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的最高位為1j=j-k;/把最高位變成0k=k/2;/k/

4、2,比較次高位，依次類推，逐個比較，直到某個位為0j=j+k;把0改為1/FFT運算核，使用蝶形運算完成FFT運算/計算l的值，即計算蝶形級數(shù)int le,lei,ip;f=FFT_N;Jfor(m=1;m<=l;m+)le=2<<(m-1);lei=le/2;u.real=1.0;u.imag=0.0;for(l=1；(f=f/2)!=l；l+)/控制蝶形結(jié)級數(shù)/m表示第m級蝶形，l為蝶形級總數(shù)l=log (2) N /le蝶形結(jié)距離，即第m級蝶形的蝶形結(jié)相距l(xiāng)e點 /同一蝶形結(jié)中參加運算的兩點的距離/u為蝶形結(jié)運算系數(shù)，初始值為1w.real=cos(PI/lei);/w

5、為系數(shù)商，即當前系數(shù)與前一個系數(shù)的商w.imag=-sin(PI/lei);for(j=0;j<=lei-1;j+)/控制計算不同種蝶形結(jié)，即計算系數(shù)不同的蝶形結(jié)for(i=j；i<=FFT_N-1；i=i+le)/控制同一蝶形結(jié)運算，即計算系數(shù)相同蝶形結(jié) ip=i+lei;/i, ip分別表示參加蝶形運算的兩個節(jié)點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.ima

6、g; u=EE(u,w);/改變系數(shù)，進行下一個蝶形運算 I*函數(shù)原型：void main()函數(shù)功能：測試FFT變換，演示函數(shù)使用方法輸入?yún)?shù)：無輸出參數(shù)：無*/void main()int i;for(i=0;i<FFT_N;i+)/給結(jié)構(gòu)體賦值si.real=sin(2*3.141592653589793*i/FFT_N); /實部為正弦波 FFT_N 點采樣，賦值為1si.imag=0;/虛部為 0FFT(s)；/進行快速福利葉變換for(i=0;i<FFT_N;i+)/求變換后結(jié)果的模值，存入復數(shù)的實部部分si.real=sqrt(si.real*si.real+si.i

7、mag*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+; end min(y2) y=yl, y2; figure(1); plot(y); %funboxing(O.001+1/3) %/ %/快速傅里葉變換 matlab程序: %/clc;clear;clf;N=input('Node number')T=input('cai yan

8、g 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=1e=input('Input the number of zeros!')x=x zeros(1,e) N=N+e;elseend%加零k=0:N-1;%給卜重新賦值，因為有可能出現(xiàn)加零狀況bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)+1;% 利用庫

9、函數(shù)進行變址運算 for l=1:NX(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:d2for i=t:d1:N i1=i+d2; if(i1>N)break; else p=X(i1)*W; X(i1)=X(i)-p; X(i)=X(i)+p; endendW=W*dw;%做蝶形運算的兩個數(shù)之間的距離%同一級之下蝶形結(jié)之間的距離%蝶形運算系數(shù)的初始值%蝶形運算系數(shù)的增加量%判斷是否超出范圍%蝶形運算%蝶形運算系數(shù)的變化en

10、dendsubplot(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);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)/程序如下：/*FFT*/#include <stdio.h> #include <math.h&

11、gt; #include <stdlib.h> #define N 1000 typedef struct double real; double img; complex; voidfft(); /*快速傅里葉變換*/voidifft(); /*快速傅里葉逆變換*/void initW(); void change();voidadd(complex,complex,complex*);/*復數(shù)加法*/voidmul(complex,complex,complex*);/*復數(shù)乘法*/voidsub(complex,complex,complex*);/*復數(shù)減法*/void d

12、ivi(complex ,complex ,complex *);/* 復數(shù)除法 */ void output(); /* 輸出結(jié)果 */complex xN,*W;/* 輸出序列的值*/intsize_x=0;/*輸入序歹！J的長度，只限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("

13、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) or IFFT(1)?n");scanf("%d”,&method);if(method=0)fft();elseifft();output();return 0;/*進行基-2 FFT運算*/void

14、 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+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=s

15、ize_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);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 ini

16、tW()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 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>&

17、gt;1;if(j>i)temp=xi;xi=xj;xj=temp;void output()/* 輸出結(jié)果*/int i;printf("The result are as followsn");for(i=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");elseprintf("%.4fjn

18、",xi.img);void add(complexc->real=a.real+b.real; c->img=a.img+b.img; void mul(complexc->real=a.real*b.real c->img=a.real*b.img void sub(complexc->real=a.real-b.real; c->img=a.img-b.img;a,complex b,complexa,complex b,complex- a.img*b.img;+ a.img*b.real;a,complex b,complex*c)*c)

19、*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%移頻(將選帶的中心頻率移動到零頻)%數(shù)字低通濾波器(防止頻率混疊)%重新采樣(將采樣的數(shù)據(jù)再次間隔采樣，間隔的數(shù)據(jù)取決于分析的帶寬，就是放大倍數(shù))%復55丁 (由于經(jīng)過了移頻，所以數(shù)據(jù)不是實數(shù)了)%頻率調(diào)整(將負半軸的頻率成分移到正半軸)fu