快速傅里叶变换FFT结果的物理意义附已验证C51AVR单片机程序

《快速傅里叶变换FFT结果的物理意义附已验证C51AVR单片机程序》由会员分享，可在线阅读，更多相关《快速傅里叶变换FFT结果的物理意义附已验证C51AVR单片机程序（8页珍藏版）》请在装配图网上搜索。

1、FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有 些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看 出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信 号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT 之后的结果是什意思、如何决定要使用多少点来做FFT。现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC 采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的 两倍，这些我就不在此啰嗦了。采样得到的数字信号，就可以做FFT变

2、换了。N个采样点，经过FFT之后，就可 以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为 N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的 幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那 么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。 而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢， 就是在该频率下的信号的相位。第一个点表示直流分量(即0Hz)，而最后一个 点N的再下一个点(实际上这个点是不存在的，这里

3、是假设的第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。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分 辨率和采样时间是倒

4、数关系。假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是入根号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

5、、相 位为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 的

6、结果的模值如图所示。从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们 分别将这三个点附近的数据拿上来细看：1 点：512+0i2 点：-2.6195E-14 - 1.4162E-13i3 点：-2.8586E-14 - 1.1898E-13i50 点：-6.2076E-13 - 2.1713E-12i51 点：332.55 - 192i52 点：-1.6707E-12 - 1.5241E-12i75 点：-2.2199E-13 -1.0076E-12i76 点：3.4315E-12 + 192i77 点：-3.0263E-14 +7.5609E-13i很明显，1点、51

7、点、76点的值都比较大，它附近的点值都很小，可以认 为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。 分别计算这三个点的模值，结果如下：1 点：51251 点：38476 点：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,

8、结果是弧度，换算为角度就是 180*(-0.5236)/pi=-30.0001。再计算 75Hz 信号的相位，atan2(192, 3.4315E-12)=1.5708 弧度，换算成角度就是 180*1.5708/pi=90.0002。可见， 相位也是对的。根据FFT结果以及上面的分析计算，我们就可以写出信号的表达 式了，它就是我们开始提供的信号。#include /AT89C52/#include /atmeg128#include /* *快速福利叶变换C函数函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不 依赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序

9、的 复数(输入实数是可令复数虚部为0)，输出为经过FFT变换的自然顺序的复数使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N 的应该为2的N次方，不满足此条件时应在后面补0函数调用：FFT(s);作者：吉帅虎时间：2010-2-20版 本：Ver1.0参考文献：*/#include#define PI 3.1415926535897932384626433832795028841971/定义圆周率值#define FFT_N 128/定义福利叶变换的点数struct compx (float real,imag;/定义一个复数结构struct compxsFFT_

10、N;/FFT 输入和输出：从S1开始存放，根据大小自己定义/* 函数原型：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

11、(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;inm1;i+)(if(ij)如果ij,即进行变址(t=xinj;xinj=xini;xini=t;k=nv2;求j的下一个倒位序while(k=j)如果k=j,表示j的最高位为

12、1(j=j-k;/把最高位变成0k=k/2;/k/2，比较次高位，依次类推，逐个比较，直到某个位为0j=j+k;把0改为1(int le,lei,ip;算完成FFT运算f=FFT_N;for(l=1;(f=f/2)! = 1;l+) 数/FFT运算核，使用蝶形运计算l的值，即计算蝶形级7for(m=1;m=l;m+)(形级总数l=log (2) Nle=2(m-1);蝶形的蝶形结相距le点lei=le/2;两点的距离u.real=1.0;值为1u.imag=0.0;w.real=cos(PI/lei);前一个系数的商w.imag=-sin(PI/lei);for(j=0;j=lei-1;j+)

13、计算系数不同的蝶形结(for(i=j;i=FFT_N-1;i=i+le)算系数相同蝶形结(ip=i+lei;算的两个节点t=EE(xinip,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=EE(u,w);运算/控制蝶形结级数/m表示第m级蝶形，l为蝶/le蝶形结距离，即第m级 /同一蝶形结中参加运算的/u为蝶形结运算系数，初始/w为系数商，即当前系数与控制计算不同种蝶形结，即控制同一蝶形结运算，即计/i，ip分别

14、表示参加蝶形运 蝶形运算，详见公式改变系数，进行下一个蝶形/* 函数原型：void main()函数功能：测试FFT变换，演示函数使用方法输入参数：无输出参数：无*/void main()(int i;for(i=0;iFFT_N;i+)/给结构体赋值(si.real=1+2*sin(2*3.141592653589793*i/FFT_N); /实部为正弦波FFT_N点采样，赋值为1si.imag=0;/虚部为 0FFT(s);/进行快速福利叶变换for(i=0;iFFT_N;i+)/求变换后结果的模值，存入复数的实部部分si.real二sqrt(si.real*si.real+si.imag*si.imag);while(1);