当前位置:首页 >> 计算机软件及应用 >>

C语言、Matlab实现FFT几种编程实例


C 语言、MATLAB 实现 FFT 几种方法 总结前人经验,仅供参考 ///一、 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////c 语言程序////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #include <iom128.h> #include <intrinsics.h> #include<math.h> #define PI 3.1415926535897932384626433832795028841971 #define FFT_N 128 struct compx {float real,imag;}; //定义圆周率值 //定义福利叶变换的点数 //定义一个复数结构

struct compx s[FFT_N]; //FFT 输入和输出:从 S[1]开始存放,根据大小自己定义 /******************************************************************* 函数原型:struct compx EE(struct compx b1,struct compx b2) 函数功能:对两个复数进行乘法运算 输入参数:两个以联合体定义的复数 a,b 输出参数:a 和 b 的乘积,以联合体的形式输出 *******************************************************************/ 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); } /***************************************************************** 函数原型:void FFT(struct compx *xin,int N) 函数功能:对输入的复数组进行快速傅里叶变换(FFT) 输入参数:*xin 复数结构体组的首地址指针,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) { t=xin[j]; xin[j]=xin[i]; xin[i]=t; } k=nv2; while(k<=j) { j=j-k; k=k/2; 某个位为 0 } j=j+k; } { int le,lei,ip; f=FFT_N; for(l=1;(f=f/2)!=1;l++) ; for(m=1;m<=l;m++) { le=2<<(m-1); lei=le/2; u.real=1.0; u.imag=0.0; // 控制蝶形结级数 //m 表示第 m 级蝶形,l 为蝶形级总数 l=log(2)N //le 蝶形结距离,即第 m 级蝶形的蝶形结相距 le 点 //同一蝶形结中参加运算的两点的距离 //u 为蝶形结运算系数,初始值为 1 //计算 l 的值,即计算蝶形级数 //FFT 运算核, 使用蝶形运算完成 FFT 运算 //把 0 改为 1 //把最高位变成 0 //k/2,比较次高位,依次类推,逐个比较,直到 //求 j 的下一个倒位序 //如果 k<=j,表示 j 的最高位为 1 //如果 i<j,即进行变址 //变址运算,即把自然顺序变成倒位序,采用雷德算法

w.real=cos(PI/lei); w.imag=-sin(PI/lei); for(j=0;j<=lei-1;j++) {

//w 为系数商,即当前系数与前一个系数的商 //控制计算不同种蝶形结,即计算系数不同的蝶形结

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算, 即计算系数相同蝶形结 { ip=i+lei; 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); } } } } /************************************************************ 函数原型:void main() 函数功能:测试 FFT 变换,演示函数使用方法 输入参数:无 输出参数:无 ************************************************************/ void main() { int i; for(i=0;i<FFT_N;i++) { s[i].real=sin(2*3.141592653589793*i/FFT_N); //实部为正弦波 FFT_N 点采 样,赋值为 1 s[i].imag=0; } // 虚部为 0 // 给结构体赋值 //改变系数,进行下一个蝶形运算 //i,ip 分别表示参加蝶形运算的两个节点 // 蝶形运算,详见公式

FFT(s); for(i=0;i<FFT_N;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; 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 yang jian ge')

f=input('frenquency') choise=input('add zero or not? n=0:T:(N-1)*T; 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; else end k=0:N-1; for l=1:N X(l)=x(bianzhi(l));%将采样后的值按照变址运算后的顺序放入 X 矩阵中 end d1=1; for m=1:log2(N) d2=d1; d1=d1*2; W=1; dw=exp(-j*pi/d2); for t=1:d2 for 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; end end W=W*dw; endend %蝶形运算系数的变化 % 蝶形运算 %判断是否超出范围 %做蝶形运算的两个数之间的距离 %同一级之下蝶形结之间的距离 %蝶形运算系数的初始值 %蝶形运算系数的增加量 % %加零 %给 k 重新赋值,因为有可能出现加零状况 1/0 ') % 采样点

bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1;%利用库函数进行变址运算

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))); subplot(2,2,4); stem(k,abs(fft(x))/max(abs(fft(x)))); %调用系统的函数绘制结果,与自己的结果进 行比较 %画自己的 fft 的运算结果 % 画采样后的离散信号 %画原曲线

//////三、 /////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////FFT 的 C 语言算法实现 ////////////程序如下: /************FFT***********/ #include #include #include #define typedef { double double }complex; void void void void void void void fft(); /*快速傅里叶变换*/ ifft(); /*快速傅里叶逆变换*/ initW(); change(); add(complex mul(complex sub(complex ,complex ,complex ,complex ,complex ,complex ,complex *); *); *); /*复数加法*/ /*复数乘法*/ /*复数减法*/ real; img; <stdio.h> <math.h> <stdlib.h> N 1000 struct

void void

divi(complex

,complex

,complex

*);/*复数除法*/

output(); /*输出结果*/ x[N], PI; *W;/*输出序列的值*/

complex int int { int double

size_x=0;/*输入序列的长度,只限 2 的 N 次方*/ main() i,method;

system("cls"); PI=atan(1)*4; printf("Please input the size of x:\n"); /*输入序列的长度*/ scanf("%d",&size_x); printf("Please input the data in x[N]:(such as:5 6)\n"); /*输入序列对应的值*/ for(i=0;i<size_x;i++) scanf("%lf %lf",&x[i].real,&x[i].img); initW(); /*选择 FFT 或逆 FFT 运算*/ printf("Use if(method==0) fft(); else ifft(); output(); return } /*进行基-2 FFT 运算*/ void { int i=0,j=0,k=0,l=0; up,down,product; complex fft() 0; FFT(0) or IFFT(1)?\n"); scanf("%d",&method);

change(); for(i=0;i< { l=1<<i; for(j=0;j<size_x;j+= { 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; } } } } void { int i=0,j=0,k=0,l=size_x; up,down; (int)( log(size_x)/log(2) );i++) /*蝶形运算*/ complex for(i=0;i< { l/=2; for(j=0;j<size_x;j+= { 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; 2*l ) ifft() 2*l ) log(size_x)/log(2) ;i++)

x[j+k+l]=down; } } } change(); } void { int i; *)malloc(sizeof(complex) * size_x); W=(complex { W[i].real=cos(2*PI/size_x*i); W[i].img=-1*sin(2*PI/size_x*i); } } void { complex unsigned double { k=i;j=0; t=(log(size_x)/log(2)); while( { j=j<<1; j|=(k } if(j>i) { temp=x[i]; & 1); k=k>>1; (t--)>0 ) t; temp; short i=0,j=0,k=0; change() initW()

for(i=0;i<size_x;i++)

for(i=0;i<size_x;i++)

x[i]=x[j]; x[j]=temp; } } } void { int i; result are as follows\n"); printf("The { printf("%.4f",x[i].real); if(x[i].img>=0.0001) printf("+%.4fj\n",x[i].img); else else printf("%.4fj\n",x[i].img); } } void { c->real=a.real+b.real; c->img=a.img+b.img; } void { c->real=a.real*b.real c->img=a.real*b.img } void { c->real=a.real-b.real; c->img=a.img-b.img; sub(complex a,complex b,complex *c) + a.img*b.img; a.img*b.real; mul(complex a,complex b,complex *c) add(complex a,complex b,complex *c) if(fabs(x[i].img)<0.0001) printf("\n"); output() /*输出结果*/

for(i=0;i<size_x;i++)

} void { c->real=( a.real*b.real+a.img*b.img )/( divi(complex a,complex b,complex *c)

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); 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*c os(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t); [f, y] = zfft(x, 109, 115, fs); plot(f, y); %将负半轴的幅值移过来 %重新采样 %滤波 截止频率为 0.5bw

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 图片效果:


赞助商链接
相关文章:
FFTmatlab实例
FFTmatlab实例_计算机软件及应用_IT/计算机_专业资料。FFTmatlab中几个实例 close all; %先关闭所有图片 Adc=2; %直流分量幅度 A1=3; %频率 F1 信号的幅度...
用Matlab编写fft
Matlab编写fft_计算机软件及应用_IT/计算机_专业资料。使用Matlab编写fftFFT算法编写 1、程序 function [A] = myfft(A,M) N=2^M; LH=N/2; J=LH; N1...
Matlab实现FFT变换
Matlab 实现 FFT 变换 Matlab 实现 FFT 变换(单边...同样输入 x 及变换的点数, 还有一个布尔控制, 是否...下面的例子,先进行 fourier transform,即双边谱 程序...
FFT算法(用matlab实现)
健康减肥10种吃不胖的食物 吃哪些食物不发胖 在家全套...FFT的matlab程序 5页 免费 FFT快速傅里叶变换Matla...MATLAB中FFT的使用方法 5页 1下载券 C语言实现FFT...
matlab编程实例100例
matlab编程实例100例_计算机软件及应用_IT/计算机_...(fft(eye(j+16))) set(h,'value',j) m(:,...(3,5); c=rand(3,5); fill3(x,y,z,c) ...
Matlab编程实现FFT变换及频谱分析的程序代码
Matlab编程实现FFT变换及频谱分析的程序代码_IT/计算机_专业资料。用matlab实现的fft变换和频谱分析,附代码 Matlab 编程实现 FFT 变换及频谱分析的程序代码 (喜欢进行...
基于matlab的FFT算法程序设计
基于matlabFFT算法程序设计_信息与通信_工程科技_专业...除了指数的符号相反、并多了一 1/n 的因子,...李永东.AT89C51 单片机在多层楼宇对讲系统中的应用[...
fft源代码(matlab)
FFT频谱测试法Matlab代码... 2页 免费 MATLAB中的FFT实例讲解 7页 免费 喜欢...FFT的matlab程序 5页 免费 C语言实现FFT(快速傅里叶... 12页 1下载券 ...
基于matlab的fft算法程序设计
课题名称 姓学学专名号院业 基于 matlabFFT 算法程序设计 指导教师 一、设计任务及要求: 设计任务: 设计利用 FFT 的算法程序画出对正弦信号进行频谱分析。...
C语言实现FFT变换
C语言实现FFT变换_计算机软件及应用_IT/计算机_专业资料。#include <iom128.h>...此程序包采用联合体的形式表示一复数,输入为自然顺序的复 数(输入实数是可令...
更多相关标签: