整理的可用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;is[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;isin_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&&na=SIN_TAB[n];
elseif(n>=FFT_N/2&&na=-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;is[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()函数创建正弦函数表,