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

C語言、MATLAB實現(xiàn)FFT幾種方法總結(jié)前人經(jīng)驗，僅供參考///一、/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////c語言程序/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////#include<iom128.h>#include<intrinsics.h>#include<math.h>#definePI3.1415926535897932384626433832795028841971〃定義圓周率值#defineFFT_N128〃定義福利葉變換的點數(shù)structcompx{floatreal,imag;};〃定義一個復(fù)數(shù)結(jié)構(gòu)structcompxs[FFT_N];//FFT輸入和輸出：從S[1]開始存放，根據(jù)大小自己定義,*******************************************************************structcompxEE(structcompxb1,structcompxb2)對兩個復(fù)數(shù)進行乘法運算兩個以聯(lián)合體定義的復(fù)數(shù)a,ba和b的乘積，以聯(lián)合體的形式輸出structcompxEE(structcompxb1,structcompxb2)對兩個復(fù)數(shù)進行乘法運算兩個以聯(lián)合體定義的復(fù)數(shù)a,ba和b的乘積，以聯(lián)合體的形式輸出structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/*****************************************************************函數(shù)原型：voidFFT(structcompx*xin,intN)對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT)*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型函數(shù)功能輸入?yún)?shù)*****************************************************************/對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT)*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型intf,m,nv2,nm1,i,k,l,j=0;structcompxu,w,t;nv2=FFT_N/2;〃變址運算，即把自然順序變成倒位序，采用雷德算法nm1=FFT_N-1;for(i=0;i<nm1;i++){if(i<j)〃如果i<j,即進行變址{t=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2;〃求j的下一個倒位序while(k<=j)〃如果k<=j,表示j的最高位為1{j=j-k;〃把最高位變成0k=k/2;//k/2，比較次高位，依次類推，逐個比較，直到某個位為0}j=j+k;〃把0改為1}{//FFT運算核，使用蝶形運算完成FFT運算//FFT運算核，使用蝶形運算完成FFT運算//計算l的值，即計算蝶形級數(shù)f=FFT_N;;for(m=1;m<=l;m++){le=2<<(m-1);;for(m=1;m<=l;m++){le=2<<(m-1);lei=le/2;u.real=1.0;u.imag=0.0;//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為系數(shù)商，即當(dāng)前系數(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(xin[ip],u);〃蝶形運算，詳見公式xin[ip].real=xin[i].real-t.real;xin[ip].imag=xin[i].imag-t.imag;xin[i].real=xin[i].real+t.real;xin[i].imag=xin[i].imag+t.imag;}u=EE(u,w);〃改變系數(shù)，進行下一個蝶形運算}}}/************************************************************函數(shù)原型：voidmain()函數(shù)功能：測試FFT變換，演示函數(shù)使用方法輸入?yún)?shù)：無輸出參數(shù)：無************************************************************/voidmain(){inti;for(i=0;i<FFT_N;i++)〃給結(jié)構(gòu)體賦值{s[i].real=sin(2*3.141592653589793*i/FFT_N);〃實部為正弦波FFT_N點采樣，賦值為1s[i].imag=0;//虛部為0FFT(s);〃進行快速福利葉變換FFT(s);for(i=0;i<FFT_N;i++)〃求變換后結(jié)果的模值，存入復(fù)數(shù)的實部部分s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);while(1);}%////二、%/////////////////////////////////////////////////////////////////////////////////////////////////////////////%///////////////////////////////////////////////////////////////////////////////////////////////////////////%////////////////////////////////MATLAB仿真信號源的源程序::Clear;Clc;t=O:O.01:3;yl=100*sin(pi/3*t);n=l;fort=-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;clf;N=input('Nodenumber')T=input('caiyangjiange')

f=input('frenquency')choise=input('addzeroornot?1/0')n=0:T:(N-1)*T;%采樣點k=0:N-1;x=sin(2*pi*f*n);ifchoise==1e=input('Inputthenumberofzeros!')x=[xzeros(1,e)]N=N+e;elseend%加零k=0:N-1;%給k重新賦值，因為有可能出現(xiàn)加零狀況bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;%利用庫函數(shù)進行變址運算forl=1:NX(l)=x(bianzhi(l));%將采樣后的值按照變址運算后的順序放入X矩陣中endd1=1;form=1:log2(N)%做蝶形運算的兩個數(shù)之間的距離%同一級之下蝶形結(jié)之間的距離%蝶形運算系數(shù)的初始值%蝶形運算系數(shù)的增加量%%判斷是否超出范圍d2=d1;d1=d1*2;W=1;dw=exp(-j*pi/d2);fort=1:d2fori=t:d1:Ni1=i+d2;if(i1>N)break;elsep=X(i1)*W;X(i1)=X(i)-p;%做蝶形運算的兩個數(shù)之間的距離%同一級之下蝶形結(jié)之間的距離%蝶形運算系數(shù)的初始值%蝶形運算系數(shù)的增加量%%判斷是否超出范圍%蝶形運算系數(shù)的變化subplot(2,2,1);%畫原曲線%畫原曲線plot(t,sin(2*pi*f*t));subplot(2,2,2);%畫采樣后的離散信號%畫采樣后的離散信號%畫自己的%畫自己的fft的運算結(jié)果stem(k,abs(X)/max(abs(X)));subplot(2,2,4);stem(k,abs(fft(x))/max(abs(fft(x))));%調(diào)用系統(tǒng)的函數(shù)繪制結(jié)果，與自己的結(jié)果進行比較//////三、IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////FFT的C語言算法實現(xiàn)////////////程序如下：3“““““““““““I'Ivii“““““““““““//************FFT***********/#include<stdio.h>#include<math.h>#include<stdlib.h>#defineN1000structtypedefstruct{real;doublereal;img;}complex;voiddoublefft();/*快速傅里葉變換*/ifft();/*快速傅里葉逆變換*/voidinitW();voidchange();voidadd(complex,complex,complex*);/*復(fù)數(shù)加法*/voidmul(complex,complex,complex*);/*復(fù)數(shù)乘法*/voidsub(complex,complex,complex*);/*復(fù)數(shù)減法*/voidvoiddivi(complex,complex,complex*);/*復(fù)數(shù)除法*/voidoutput。；/*輸出結(jié)果*/complexx[N],*W;/*輸出序列的值*/intsize_x=0;/*輸入序列的長度，只限2的N次方*/doublePI;intmain(){inti,method;system("cls");PI=atan(1)*4;printf("Pleaseinputthesizeofx:\n");/*輸入序列的長度*/scanf("%d",&size_x);printf("Pleaseinputthedatainx[N]:(suchas:56)\n");/*輸入序列對應(yīng)的值*/for(i=0;i<size_x;i++)scanf("%lf%lf”,&x[i].real,&x[i].img);initW();/*選擇FFT或逆FFT運算*/printf("UseFFT(0)orIFFT(1)?\n");scanf("%d",&method);if(method==0)fft();elseifft();output();return0;}/*進行基-2FFT運算*/voidfft(){inti=0,j=0,k=0,l=0;complexup,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(x[j+k+l],W[size_x*k/2/l],&product);add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}voidifft(){inti=0,j=0,k=0,l=size_x;complexup,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(x[j+k],x[j+k+l],&up);up.real/=2;up.img/=2;sub(x[j+k],x[j+k+l],&down);down.real/=2;down.img/=2;divi(down,W[size_x*k/2/l],&down);x[j+k]=up;x[j+k+l]=down;}}}change();}voidinitW(){inti;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;i<size_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}voidchange(){complextemp;unsignedshorti=0,j=0,k=0;doublet;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=x[i];

x[i]=x[j];x[j]=temp;}}}voidoutput()/*輸出結(jié)果*/{inti;printf("Theresultareasfollows\n");for(i=0;i<size_x;i++){printf("%.4f",x[i].real);if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img);elseif(fabs(x[i].img)<0.0001)printf("\n");elseprintf("%.4fj\n",x[i].img);}}voidadd(complex{c->real=a.real+b.real;c->img=a.img+b.img;}voidmul(complex{c->real=a.real*b.realc->img=a.real*b.img}voidsub(complex{c->real=a.real-b.real;c->img=a.img-b.img;a,complexb,complexa,complexb,complex-a.img*b.img;+a.img*b.real;a,complexb,complex*c)*c)*c)}voiddivi(complexa,complexb,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ù))%復(fù)FFT(由于經(jīng)過了移頻，所以數(shù)據(jù)不是實數(shù)了)%頻率調(diào)整(將負半軸的頻率成分移到正半軸)function[f,y]=zfft(x,fi,fa,fs)%x為采集的數(shù)據(jù)%fi為分析的起始頻率%fa為分析的截止頻率%fs為采集數(shù)據(jù)的采樣頻率%f為輸出的頻率序列%y為輸出的幅值序列(實數(shù))f0=(fi+fa)/2;%