2022年C语言、Matlab实现FFT几种编程实例...docx

《2022年C语言、Matlab实现FFT几种编程实例...docx》由会员分享，可在线阅读，更多相关《2022年C语言、Matlab实现FFT几种编程实例...docx（22页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -C 语言、 MATLAB 实现 FFT 几种方法总结前人体会,仅供参考/一、/ /c语言程序 / / #include #include #include #define PI 3.1415926535897932384626433832795028841971 /定义圆周率值#define FFT_N 128 struct compx float real,imag; /定义福利叶变换的点数/定义一个复数结构struct compx sFFT_N; /FFT 输入和输出：从 S1开头存放,依据大小自己定义

2、/* 函数原型： struct compx EEstruct compx b1,struct compx b2 函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数 a,b 输出参数： a 和 b 的乘积,以联合体的形式输出 */ struct compx EEstruct 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; returnc; /* 函数原型： void FFTstruct compx *xin,i

3、nt N 函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数： *xin 复数结构体组的首地址指针,struct 型*/ void FFTstruct compx *xin 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 1 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -int f,m,nv2,nm1,i,k,l,j=0; struct compx u,w,t; nv2=FFT_N/2; /变址运算,即把自然次序变成倒位序,采纳雷德算法 nm1=F

4、FT_N-1; fori=0;inm1;i+ ifij /假如 ij,即进行变址 t=xinj; xinj=xini; xini=t; k=nv2; /求 j 的下一个倒位序whilek=j /假如 k=j, 表示 j 的最高位为 1 j=j-k; /把最高位变成 0 k=k/2; /k/2,比较次高位,依次类推,逐个比较,直到 某个位为 0 j=j+k; /把 0 改为 1 int le,lei,ip; /FFT 运算核,使用蝶形运算完成FFT 运算f=FFT_N; forl=1;f=f/2.=1;l+ ; /运算 l 的值,即运算蝶形级数form=1;m=l;m+ / 掌握蝶形结级数 /m

5、 表示第 m 级蝶形, l 为蝶形级总数 l=log（2）N le=2m-1; /le 蝶形结距离,即第 m 级蝶形的蝶形结相距 le 点lei=le/2; /同一蝶形结中参与运算的两点的距离u.real=1.0; /u 为蝶形结运算系数,初始值为 1 u.imag=0.0; 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 2 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -w.real=cosPI/lei; /w 为系数商,即当前系数与前一个系数的商w.im

6、ag=-sinPI/lei; forj=0;j=lei-1;j+ /掌握运算不同种蝶形结,即运算系数不同的蝶形结fori=j;i=FFT_N-1;i=i+le / 掌握同一蝶形结运算, 即运算系数相同蝶形结 ip=i+lei; /i,ip 分别表示参与蝶形运算的两个节点 t=EExinip,u; /蝶形运算,详见公式 xinip.real=xini.real-t.real; xinip.imag=xini.imag-t.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag; u=EEu,w; /转变系数,进行下一个蝶形运算 /

7、* 函数原型： void main 函数功能：测试 FFT 变换,演示函数使用方法 输入参数：无输出参数：无*/ void main int i; fori=0;iFFT_N;i+ /给结构体赋值 si.real=sin2*3.141592653589793*i/FFT_N; /实部为正弦波 FFT_N 点采 样,赋值为 1 si.imag=0; /虚部为 0 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 3 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -

8、FFTs; /进行快速福利叶变换fori=0;iNbreak; %判定是否超出范畴 else p=Xi1*W; Xi1=Xi-p; Xi=Xi+p; %蝶形运算 end end W=W*dw; %蝶形运算系数的变化 endend 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 5 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -subplot2,2,1; t=0:0.0000001:N*T; plott,sin2*pi*f*t; %画原曲线 subplot2,

9、2,2; stemk,x; %画采样后的离散信号 subplot2,2,3; stemk,absX/maxabsX; subplot2,2,4; %画自己的 fft 的运算结果stemk,absfftx/maxabsfftx; % 调用系统的函数绘制结果,与自己的结果进 行比较/ 三、/ / /FFT 的 C 语言算法实现 / 程序如下：/*FFT*/ #include #include #include #define N 1000 typedef struct double real; double img; complex; 细心整理归纳 精选学习资料 void fft; /* 快速傅里

10、叶变换 */ *; /* 复数加法 */ 第 6 页,共 12 页 void ifft; /* 快速傅里叶逆变换 */ void initW; void change; void addcomplex ,complex ,complex void mulcomplex ,complex ,complex *; /* 复数乘法 */ void subcomplex ,complex ,complex *; /* 复数减法 */ - - - - - - - - - - - - - - - - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - -

11、 - - -void divicomplex ,complex ,complex *;/* 复数除法 */ void output; /* 输出结果 */ complex xN, *W;/* 输出序列的值 */ int size_x=0;/*输入序列的长度,只限 2 的 N 次方 */ double PI; int main int i,method; systemcls; PI=atan1*4; printfPlease input the size of x:n; /* 输入序列的长度 */ scanf%d,&size_x; printfPlease input the data in x

12、N:such as:5 6n; /* 输入序列对应的值 */ fori=0;isize_x;i+ scanf%lf %lf,&xi.real,&xi.img; initW; /* 挑选 FFT 或逆 FFT 运算*/ printfUse FFT0 or IFFT1.n; scanf%d,&method; ifmethod=0 fft; else ifft; output; return 0; /* 进行基 -2 FFT 运算 */ void fft 细心整理归纳 精选学习资料 int i=0,j=0,k=0,l=0; 第 7 页,共 12 页 complex up,down,product;

13、- - - - - - - - - - - - - - - - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -change; fori=0;i logsize_x/log2 ;i+ l=1i; forj=0;jsize_x;j+= 2*l fork=0;kl;k+ mulxj+k+l,Wsize_x*k/2/l,&product; addxj+k,product,&up; subxj+k,product,&down; xj+k=up; xj+k+l=down; void ifft int i=0,j=0,k=0,l=siz

14、e_x; complex up,down; fori=0;i int logsize_x/log2 ;i+ /* 蝶形运算 */ l/=2; forj=0;jsize_x;j+= 2*l fork=0;kl;k+ addxj+k,xj+k+l,&up; up.real/=2;up.img/=2; subxj+k,xj+k+l,&down; down.real/=2;down.img/=2; dividown,Wsize_x*k/2/l,&down; xj+k=up; 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 8 页,共 12 页 - - - -

15、 - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -xj+k+l=down; change; void initW int i; *mallocsizeofcomplex * size_x; W=complex fori=0;isize_x;i+ Wi.real=cos2*PI/size_x*i; Wi.img=-1*sin2*PI/size_x*i; void change complex temp; i=0,j=0,k=0; unsigned short double t; fori=0;i0 j=j1; ifji temp=xi;

16、细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 9 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -xi=xj; xj=temp; void output /* 输出结果 */ int i; result are as followsn; printfThe fori=0;i=0.0001 printf+%.4fjn,xi.img; else iffabsxi.imgreal=a.real+b.real; c-img=a.img+b.img; void mul

17、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 subcomplex a,complex b,complex *c c-real=a.real-b.real; c-img=a.img-b.img; 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 10 页,共 12 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - -

18、- void divicomplex 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 %移频 （将选带的中心频率移动到零频）%数字低通滤波器（防止频率混叠）%重新采样（将采样的数据再次间隔采样,间隔的数据取决于分析的带宽,就 是放大倍数）%复 FFT （由于经过了移频,所以数据不是实数了

19、）%频率调整（将负半轴的频率成分移到正半轴）function f, y = zfftx, fi, fa, fs % x 为采集的数据 % fi 为分析的起始频率 % fa 为分析的截止频率 % fs 为采集数据的采样频率 % f 为输出的频率序列 % y 为输出的幅值序列 实数）f0 = fi + fa / 2; % 中心频率 N = lengthx; %数据长度r = 0:N-1; b = 2*pi*f0.*r ./ fs; x1 = x .* exp-1j .* b; % 移频细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 11 页,共 12 页

20、 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -bw = fa - fi; B = fir132, bw / fs; %滤波 截止频率为 0.5bw x2 = filterB, 1, x1; c = x21:floorfs/bw:N; % 重新采样 N1 = lengthc; f = linspacefi, fa, N1; y = absfftc ./ N1 * 2; y = circshifty, 0, floorN1/2; % 将负半轴的幅值移过来 end / / %上边 四程序 测试应用实例：fs = 2048; T = 100; t = 0:1/fs:T; x = 30 * cos2*pi*110.*t + 30 * cos2*pi*111.45.*t + 25*cos2*pi*112.3*t + 48*c os2*pi*113.8.*t+50*cos2*pi*114.5.*t; f, y = zfftx, 109, 115, fs; plotf, y; / 图片成效：细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 12 页,共 12 页 - - - - - - - - -