C语言Matlab实现FFT几种编程实例.docx
- 文档编号:5601435
- 上传时间:2022-12-28
- 格式:DOCX
- 页数:14
- 大小:30.16KB
C语言Matlab实现FFT几种编程实例.docx
《C语言Matlab实现FFT几种编程实例.docx》由会员分享,可在线阅读,更多相关《C语言Matlab实现FFT几种编程实例.docx(14页珍藏版)》请在冰豆网上搜索。
C语言Matlab实现FFT几种编程实例
C语言、MATLAB实现FFT几种方法
总结前人经验,仅供参考
///一、
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////c语言程序//////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#include
#include
#include
#definePI3.1415926535897932384626433832795028841971//定义圆周率值
#defineFFT_N128//定义福利叶变换的点数
structcompx{floatreal,imag;};//定义一个复数结构
structcompxs[FFT_N];//FFT输入和输出:
从S[1]开始存放,根据大小自己定义
/*******************************************************************
函数原型:
structcompxEE(structcompxb1,structcompxb2)
函数功能:
对两个复数进行乘法运算
输入参数:
两个以联合体定义的复数a,b
输出参数:
a和b的乘积,以联合体的形式输出
*******************************************************************/
structcompxEE(structcompxa,structcompxb)
{
structcompxc;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
/*****************************************************************
函数原型:
voidFFT(structcompx*xin,intN)
函数功能:
对输入的复数组进行快速傅里叶变换(FFT)
输入参数:
*xin复数结构体组的首地址指针,struct型
*****************************************************************/
voidFFT(structcompx*xin)
{
intf,m,nv2,nm1,i,k,l,j=0;
structcompxu,w,t;
nv2=FFT_N/2;//变址运算,即把自然顺序变成倒位序,采用雷德算法
nm1=FFT_N-1;
for(i=0;i { if(i { t=xin[j]; xin[j]=xin[i]; xin[i]=t; } k=nv2;//求j的下一个倒位序 while(k<=j)//如果k<=j,表示j的最高位为1 { j=j-k;//把最高位变成0 k=k/2;//k/2,比较次高位,依次类推,逐个比较,直到某个位为0 } j=j+k;//把0改为1 } { intle,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为蝶形结运算系数,初始值为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(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);//改变系数,进行下一个蝶形运算 } } } } /************************************************************ 函数原型: voidmain() 函数功能: 测试FFT变换,演示函数使用方法 输入参数: 无 输出参数: 无 ************************************************************/ voidmain() { inti; for(i=0;i { s[i].real=sin(2*3.141592653589793*i/FFT_N);//实部为正弦波FFT_N点采样,赋值为1 s[i].imag=0;//虚部为0 } FFT(s);//进行快速福利叶变换 for(i=0;i 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+; end min(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==1 e=input('Inputthenumberofzeros! ') x=[xzeros(1,e)] N=N+e; else end%加零 k=0: N-1;%给k重新赋值,因为有可能出现加零状况 bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;%利用库函数进行变址运算 forl=1: N X(l)=x(bianzhi(l));%将采样后的值按照变址运算后的顺序放入X矩阵中 end d1=1; form=1: log2(N) d2=d1;%做蝶形运算的两个数之间的距离 d1=d1*2;%同一级之下蝶形结之间的距离 W=1;%蝶形运算系数的初始值 dw=exp(-j*pi/d2);%蝶形运算系数的增加量 fort=1: d2% fori=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;%蝶形运算 end end W=W*dw;%蝶形运算系数的变化 endend subplot(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))));%调用系统的函数绘制结果,与自己的结果进行比较 //////三、 /////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////FFT的C语言算法实现 ////////////程序如下: /************FFT***********/ #include #include #include #defineN1000 typedefstruct { doublereal; doubleimg; }complex; voidfft();/*快速傅里叶变换*/ voidifft();/*快速傅里叶逆变换*/ voidinitW(); voidchange(); voidadd(complex,complex,complex*);/*复数加法*/ voidmul(complex,complex,complex*);/*复数乘法*/ voidsub(complex,complex,complex*);/*复数减法*/ voiddivi(complex,complex,complex*);/*复数除法*/ voidoutput();/*输出结果*/ 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"); /*输入序列对应的值*/ for(i=0;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(); else ifft(); output(); return0; } /*进行基-2FFT运算*/ voidfft() { inti=0,j=0,k=0,l=0; complexup,down,product; change(); for(i=0;i (2);i++) { l=1< for(j=0;j { for(k=0;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 { for(k=0;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 { 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 { 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()/*输出结果*/ { inti; printf("Theresultareasfollows\n"); for(i=0;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"); else printf("%.4fj\n",x[i].img); } } voidadd(complexa,complexb,complex*c) { c->real=a.real+b.real; c->img=a.img+b.img; } voidmul(complexa,complexb,complex*c) { c->real=a.real*b.real-a.img*b.img; c->img=a.real*b.img+a.img*b.real; } voidsub(complexa,complexb,complex*c) { c->real=a.real-b.real; c->img=a.img-b.img; } 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 %移频(将选带的中心频率移动到零频) %数字低通滤波器 (防止频率混叠) %重新采样 (将采样的数据再次间隔采样,间隔的数据取决于分析的带宽,就是放大倍数) %复FFT(由于经过了移频,所以数据不是实数了) %频率调整(将负半轴的频率成分移到正半轴) function [f, y] = zfft(x, fi, fa, fs) % x为采集的数据 % fi为分析的起始频率 % fa为分析的截止频率 % fs为采集数据的采样频率 % f为输出的频率序列 % y为输出的幅值序列(实数) f0 = (fi + fa) / 2; %中心频率 N = length(x); %数据长度 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 = length(c); f = linspace(fi, fa, N1); y = abs(fft(c)) ./ N1 * 2; y = circshift(y, [0, floor(N1/2)]); %将负半轴的幅值移过来 end ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////// %上边四程序测试应用实例: fs = 2048; T = 100; t = 0: 1/fs: T; x = 30 * cos(2*pi*110.*t) + 30 * cos(2*pi*111.45.*t) + 25*cos(2*pi*112.3*t) + 48*cos(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t); [f, y] = zfft(x, 109, 115, fs); plot(f, y); ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 图片效果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 Matlab 实现 FFT 编程 实例
![提示](https://static.bdocx.com/images/bang_tan.gif)