整理的可用FFT原理及程序.docx
- 文档编号:30075220
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:24
- 大小:21.62KB
整理的可用FFT原理及程序.docx
《整理的可用FFT原理及程序.docx》由会员分享,可在线阅读,更多相关《整理的可用FFT原理及程序.docx(24页珍藏版)》请在冰豆网上搜索。
整理的可用FFT原理及程序
FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。
有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。
这就是很多信号分析采用FFT变换的原因。
另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。
虽然很多人都知道FFT是什么,可以用来做什么,怎么去做,但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。
现在就根据实际经验来说说FFT结果的具体物理意义。
一个模拟信号,经过ADC采样之后,就变成了数字信号。
采样定理告诉我们,采样频率要大于信号频率的两倍,这些我就不在此啰嗦了。
采样得到的数字信号,就可以做FFT变换了。
N个采样点,经过FFT之后,就可以得到N个点的FFT结果。
为了方便进行FFT运算,通常N取2的整数次方。
假设采样频率为Fs,信号频率F,采样点数为N。
那么FFT之后结果就是一个为N点的复数。
每一个点就对应着一个频率点。
这个点的模值,就是该频率值下的幅度特性。
具体跟原始信号的幅度有什么关系呢?
假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。
而第一个点就是直流分量,它的模值就是直流分量的N倍。
而每个点的相位呢,就是在该频率下的信号的相位。
第一个点表示直流分量(即0Hz),而最后一个点N的再下一个点(实际上这个点是不存在的,这里是假设的第N+1个点,也可以看做是将第一个点分做两半分,另一半移到最后)则表示采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率依次增加。
例如某点n所表示的频率为:
Fn=(n-1)*Fs/N。
由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。
1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。
如果要提高频率分辨力,则必须增加采样点数,也即采样时间。
频率分辨率和采样时间是倒数关系。
假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是An=根号a*a+b*b,相位就是Pn=atan2(b,a)。
根据以上的结果,就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:
An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。
对于n=1点的信号,是直流分量,幅度即为A1/N。
由于FFT结果的对称性,通常我们只使用前半部分的结果,即小于采样频率一半的结果。
下面以一个实际的信号来做说明。
假设我们有一个信号,它含有2V的直流分量,频率为50Hz、相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。
用数学表达式就是如下:
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。
式中cos参数为弧度,所以-30度和90度要分别换算成弧度。
我们以256Hz的采样率对这个信号进行采样,总共采样256点。
按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。
我们的信号有3个频率:
0Hz、50Hz、75Hz,应该分别在第1个点、第50个点、第76个点上出现峰值,其它各点应该接近0。
实际情况如何呢?
我们来看看FFT的结果的模值如图所示。
从图中我们可以看到,在第1点、第51点、和第76点附近有比较大的值。
我们分别将这三个点附近的数据拿上来细看:
1点:
512+0i
2点:
-2.6195E-14-1.4162E-13i
3点:
-2.8586E-14-1.1898E-13i
50点:
-6.2076E-13-2.1713E-12i
51点:
332.55-192i
52点:
-1.6707E-12-1.5241E-12i
75点:
-2.2199E-13-1.0076E-12i
76点:
3.4315E-12+192i
77点:
-3.0263E-14+7.5609E-13i
很明显,1点、51点、76点的值都比较大,它附近的点值都很小,可以认为是0,即在那些频率点上的信号幅度为0。
接着,我们来计算各点的幅度值。
分别计算这三个点的模值,结果如下:
1点:
512
51点:
384
76点:
192
按照公式,可以计算出直流分量为:
512/N=512/256=2;50Hz信号的幅度为:
384/(N/2)=384/(256/2)=3;75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。
可见,从频谱分析出来的幅度是正确的。
然后再来计算相位信息。
直流信号没有相位可言,不用管它。
先计算50Hz信号的相位,atan2(-192,332.55)=-0.5236,结果是弧度,换算为角度就是180*(-0.5236)/pi=-30.0001。
再计算75Hz信号的相位,atan2(192,3.4315E-12)=1.5708弧度,换算成角度就是180*1.5708/pi=90.0002。
可见,相位也是对的。
根据FFT结果以及上面的分析计算,我们就可以写出信号的表达式了,它就是我们开始提供的信号。
Ver1.0
---------------------------------------------------------------------------------------------------------------------------------
#include
#include
/*********************************************************************
快速福利叶变换C函数
函数简介:
此函数是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依
赖硬件。
此函数采用联合体的形式表示一个复数,输入为自然顺序的复
数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的
复数
使用说明:
使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的
应该为2的N次方,不满足此条件时应在后面补0
函数调用:
FFT(s);
作者:
吉帅虎
时间:
2010-2-20
版本:
Ver1.0
参考文献:
**********************************************************************/
#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); } Ver1.1 #include #include /********************************************************************* 快速福利叶变换C程序包 函数简介: 此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依 赖硬件。 此程序包采用联合体的形式表示一个复数,输入为自然顺序的复 数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表, 以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度 使用说明: 使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的 应该为2的N次方,不满足此条件时应在后面补0。 若使用查表法计算sin值和 cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表 函数调用: FFT(s); 作者: 吉帅虎 时间: 2010-2-20 版本: Ver1.1 参考文献: **********************************************************************/ #include #defineFFT_N128//定义福利叶变换的点数 #definePI3.1415926535897932384626433832795028841971//定义圆周率值 structcompx{floatreal,imag;};//定义一个复数结构 structcompxs[FFT_N];//FFT输入和输出: 从S[0]开始存放,根据大小自己定义 floatSIN_TAB[FFT_N/2];//定义正弦表的存放空间 /******************************************************************* 函数原型: 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); } /****************************************************************** 函数原型: voidcreate_sin_tab(float*sin_t) 函数功能: 创建一个正弦采样表,采样点数与福利叶变换点数相同 输入参数: *sin_t存放正弦表的数组指针 输出参数: 无 ******************************************************************/ voidcreate_sin_tab(float*sin_t) { inti; for(i=0;i sin_t[i]=sin(2*PI*i/FFT_N); } /****************************************************************** 函数原型: voidsin_tab(floatpi) 函数功能: 采用查表的方法计算一个数的正弦值 输入参数: pi所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换 输出参数: 输入值pi的正弦值 ******************************************************************/ floatsin_tab(floatpi) { intn; floata; n=(int)(pi*FFT_N/2/PI); if(n>=0&&n a=SIN_TAB[n]; elseif(n>=FFT_N/2&&n a=-SIN_TAB[n-FFT_N/2]; returna; } /****************************************************************** 函数原型: voidcos_tab(floatpi) 函数功能: 采用查表的方法计算一个数的余弦值 输入参数: pi所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换 输出参数: 输入值pi的余弦值 ******************************************************************/ floatcos_tab(floatpi) { floata,pi2; pi2=pi+PI/2; if(pi2>2*PI) pi2-=2*PI; a=sin_tab(pi2); returna; } /***************************************************************** 函数原型: 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);//不适用查表法计算sin值和cos值 //w.imag=-sin(PI/lei); w.real=cos_tab(PI/lei);//w为系数商,即当前系数与前一个系数的商 w.imag=-sin_tab(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; create_sin_tab(SIN_TAB); 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); } Ver1.2 #include #include /********************************************************************* 快速福利叶变换C程序包 函数简介: 此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依 赖硬件。 此程序包采用联合体的形式表示一个复数,输入为自然顺序的复 数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 整理 可用 FFT 原理 程序