C語(yǔ)言、Matlab實(shí)現(xiàn)FFT幾種編程實(shí)例

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

2、truct compx EE(struct compx b1,struct compx b2) 函數(shù)功能：對(duì)兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算輸入?yún)?shù)：兩個(gè)以聯(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);/*函數(shù)原型：void FFT(struct compx *xin,int N)函數(shù)功能：對(duì)輸入的

3、復(fù)數(shù)組進(jìn)行快速傅里葉變換（FFT）輸入?yún)?shù)：*xin復(fù)數(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; /變址運(yùn)算，即把自然順序變成倒位序，采用雷德算法 nm1=FFT_N-1; for(i=0;i<nm1;i+) if(i<j) /如果i<j,即進(jìn)行變址 t=xinj; xinj=xini; xini=t; k=nv2; /求j的下一個(gè)倒位序 while(k<=j) /如果k<=j,表示j的最高位

4、為1 j=j-k; /把最高位變成0 k=k/2; /k/2，比較次高位，依次類(lèi)推，逐個(gè)比較，直到某個(gè)位為0 j=j+k; /把0改為1 int le,lei,ip; /FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算 f=FFT_N; for(l=1;(f=f/2)!=1;l+) /計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù) ; for(m=1;m<=l;m+) / 控制蝶形結(jié)級(jí)數(shù) /m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)l=log（2）N le=2<<(m-1); /le蝶形結(jié)距離，即第m級(jí)蝶形的蝶形結(jié)相距l(xiāng)e點(diǎn) lei=le/2; /同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離 u.real=1.0; /u為蝶

5、形結(jié)運(yùn)算系數(shù)，初始值為1 u.imag=0.0; w.real=cos(PI/lei); /w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(shù)的商 w.imag=-sin(PI/lei); for(j=0;j<=lei-1;j+) /控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié) for(i=j;i<=FFT_N-1;i=i+le) /控制同一蝶形結(jié)運(yùn)算，即計(jì)算系數(shù)相同蝶形結(jié) ip=i+lei; /i，ip分別表示參加蝶形運(yùn)算的兩個(gè)節(jié)點(diǎn) t=EE(xinip,u); /蝶形運(yùn)算，詳見(jiàn)公式 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); /改變系數(shù)，進(jìn)行下一個(gè)蝶形運(yùn)算 /*函數(shù)原型：void main() 函數(shù)功能：測(cè)試FFT變換，演示函數(shù)使用方法輸入?yún)?shù)：無(wú)輸出參數(shù)：無(wú)*/void main() int i; for(i=0;i<FFT_N;i+) /給結(jié)構(gòu)體賦值 si.real=sin(2*3.141592653589793*i/FFT_N); /實(shí)部為正弦波FFT_N點(diǎn)采樣，賦值為1 si.imag=0; /虛部為0 FFT(s); /進(jìn)行快速福利葉變換 for(i=0;i<

7、;FFT_N;i+) /求變換后結(jié)果的模值，存入復(fù)數(shù)的實(shí)部部分 si.real=sqrt(si.real*si.real+si.imag*si.imag); while(1); %/二、%/%/%/MATLAB仿真信號(hào)源的源程序: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; %采樣點(diǎn)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重新賦值，因?yàn)橛锌赡艹霈F(xiàn)

9、加零狀況 bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)+1;%利用庫(kù)函數(shù)進(jìn)行變址運(yùn)算for l=1:N X(l)=x(bianzhi(l);%將采樣后的值按照變址運(yùn)算后的順序放入X矩陣中endd1=1;for m=1:log2(N) d2=d1; %做蝶形運(yùn)算的兩個(gè)數(shù)之間的距離 d1=d1*2; %同一級(jí)之下蝶形結(jié)之間的距離 W=1; %蝶形運(yùn)算系數(shù)的初始值 dw=exp(-j*pi/d2); %蝶形運(yùn)算系數(shù)的增加量 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; %蝶形運(yùn)算 end end W=W*dw; %蝶形運(yùn)算系數(shù)的變化 endendsubplot(2,2,1);t=0:0.0000001:N*T;plot(t,sin(2*pi*f*t); %畫(huà)原曲線subplot(2,2,2);stem(k,x); %畫(huà)采樣后的離散信號(hào)subplot(2,2,3);stem(k,abs(X)/max(abs(X); %畫(huà)自己的fft的運(yùn)算結(jié)果subplot(2,2,4);stem(k,abs(fft(x)/max(abs(fft(x); %調(diào)用系統(tǒng)的函數(shù)繪制結(jié)果，與自己的結(jié)果進(jìn)行

11、比較/三、/FFT的C語(yǔ)言算法實(shí)現(xiàn)/程序如下： /*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 *); /*復(fù)數(shù)加法*/ vo

12、id 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 size_x=0;/*輸入序列的長(zhǎng)度，只限2的N次方*/ double PI; int main() int i,method; system("cls"); PI=atan(1)*4;

13、printf("Please input the size of x:n"); /*輸入序列的長(zhǎng)度*/ scanf("%d",&size_x); printf("Please input the data in xN:(such as:5 6)n"); /*輸入序列對(duì)應(yīng)的值*/ for(i=0;i<size_x;i+) scanf("%lf %lf",&xi.real,&xi.img); initW(); /*選擇FFT或逆FFT運(yùn)算*/ printf("Use FFT(0)

14、 or IFFT(1)?n"); scanf("%d",&method); if(method=0) fft(); else ifft(); output(); return 0; /*進(jìn)行基-2 FFT運(yùn)算*/ 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+) /*蝶形運(yùn)算*/ 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() /*輸出結(jié)果*/ 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%移頻 （將選帶的中心頻率移動(dòng)到零頻）%數(shù)字低通濾波器  （防止頻率混疊）%重新采樣  （將采樣的數(shù)據(jù)再次間隔采樣，間隔的數(shù)據(jù)取決于分析的帶寬，就是放大倍數(shù)）%復(fù)FFT （由于經(jīng)過(guò)了移頻，所以數(shù)據(jù)不是實(shí)數(shù)了）%頻率調(diào)整 （將負(fù)半軸的頻率成分移到正半軸）function f, y = zfft(x,

21、0;fi, fa, fs)  % x為采集的數(shù)據(jù)  % fi為分析的起始頻率  % fa為分析的截止頻率  % fs為采集數(shù)據(jù)的采樣頻率  % f為輸出的頻率序列  % y為輸出的幅值序列(實(shí)數(shù)）    f0 = (fi + fa) / 2;     

22、;         %中心頻率  N = length(x);                 %數(shù)據(jù)長(zhǎng)度    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);&#