fft变换.docx

fft变换.docx

《fft变换.docx》由会员分享，可在线阅读，更多相关《fft变换.docx（14页珍藏版）》请在冰豆网上搜索。

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个点、第51个点、第76个点上出现峰值，其它各点应该接近0。

实际情况如何呢？

我们来看看FFT的结果的模值如图所示。

              图1FFT结果

  从图中我们可以看到，在第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结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。

总结：

假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：

Fn=（n-1）*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。

相位的计算可用函数atan2（b,a）计算。

atan2（b,a）是求坐标为（a,b）点的角度值，范围从-pi到pi。

要精确到xHz，则需要采样长度为1/x秒的信号，并做FFT。

要提高频率分辨率，就需要增加采样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。

解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0，使其长度达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。

具体的频率细分法可参考相关文献。

[附录：

本测试数据使用的matlab程序]

closeall;%先关闭所有图片

Adc=2;  %直流分量幅度

A1=3;  %频率F1信号的幅度

A2=1.5;%频率F2信号的幅度

F1=50;  %信号1频率（Hz）

F2=75;  %信号2频率（Hz）

Fs=256;%采样频率（Hz）

P1=-30;%信号1相位（度）

P2=90;  %信号相位（度）

N=256;  %采样点数

t=[0:

1/Fs:

N/Fs];%采样时刻

%信号

S=Adc+A1*cos（2*pi*F1*t+pi*P1/180）+A2*cos（2*pi*F2*t+pi*P2/180）;

%显示原始信号

plot（S）;

title（'原始信号'）;

figure;

Y=fft（S,N）;%做FFT变换

Ayy=（abs（Y））;%取模

plot（Ayy（1:

N））;%显示原始的FFT模值结果

title（'FFT模值'）;

figure;

Ayy=Ayy/（N/2）;  %换算成实际的幅度

Ayy

（1）=Ayy

（1）/2;

F=（[1:

N]-1）*Fs/N;%换算成实际的频率值

plot（F（1:

N/2）,Ayy（1:

N/2））;  %显示换算后的FFT模值结果

title（'幅度-频率曲线图'）;

figure;

Pyy=[1:

N/2];

fori="1:

N/2"

Pyy（i）=phase（Y（i））;%计算相位

Pyy（i）=Pyy（i）*180/pi;%换算为角度

end;

plot（F（1:

N/2）,Pyy（1:

N/2））;  %显示相位图

title（'相位-频率曲线图'）;

看完这个你就明白谐波分析了。

FFT算法原理和实现  

2011-05-1522:

16:

45|  分类：

 默认分类|举报|字号 订阅

在数字信号处理中常常需要用到离散傅立叶变换（DFT），以获取信号的频域特征。

尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。

因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散是傅立叶变换才在实际的工程中得到广泛应用。

需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。

本文就FFT的原理以及具体实现过程进行详尽讲解。

DFT计算公式

本文不加推导地直接给出DFT的计算公式：

其中x（n）表示输入的离散数字信号序列，WN为旋转因子，X（k）一组N点组成的频率成分的相对幅度。

一般情况下，假设x（n）来自于低通采样，采样频率为fs，那么X（k）表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。

因为DFT计算得到的一组离散频率幅度只实际上是在频率轴上从成周期变化的，即X（k+N）=X（k）。

因此任意取N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X（k）的周期特性，于是我们又可以认为X（k）表示了从零频率开始，频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。

N点DFT的计算量

根据

（1）式给出的DFT计算公式，我们可以知道每计算一个频率点X（k）均需要进行N次复数乘法和N-1次复数加法，计算N各点的X（k）共需要N^2次复数乘法和N*（N-1）次复数加法。

当x（n）为实数的情况下，计算N点的DFT需要2*N^2次实数乘法，2*N*（N-1）次实数加法。

旋转因子WN的特性

 1.WN的对称性

       2.WN的周期性

 3.WN的可约性

根据以上这些性质，我们可以得到式（5）的一些列有用结果

基-2FFT算法推导

假设采样序列点数为N=2^L，L为整数，如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。

首先我们将序列x（n）按照奇偶分为两组如下：

于是根据DFT计算公式

（1）有：

至此，我们将一个N点的DFT转化为了式（7）的形式，此时k的取值为0到N-1,现在分为两段来讨论，当k为0～N/2-1的时候，因为x1（r），x2（r）为N/2点的序列，因此，此时式（7）可以写为：

而当k取值为N/2~N-1时，k用k’+N/2取代，k’取值为0~N/2-1。

对式（7）化简可得：

综合以上推导我们可以得到如下结论：

一个N点的DFT变换过程可以用两个N/2点的DFT变换过程来表示，其具体公式如式（10）所示DFT快速算法的迭代公式：

上式中X'（k’）为偶数项分支的离散傅立叶变换，X''（k’’）为奇数项分支的离散傅立叶变换。

式（10）的计算过程可以用图1的蝶形算法流图直观地表示出来。

图1时间抽取法蝶形运算流图

在图1中，输入为两个N/2点的DFT输出为一个N点的DFT结果，输入输出点数一致。

运用这种表示方法，8点的DFT可以用图2来表示：

图28点DFT的4点分解

根据公式（10），一个N点的DFT可以由两个N/2点的DFT运算构成，再结合图1的蝶形信号流图可以得到图2的8点DFT的第一次分解。

该分解可以用以下几个步骤来描述：

 1.将N点的输入序列按奇偶分为2组分别为N/2点的序列

 2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2

 3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果

根据式（10）我们可以对图2中的4点DFT进一步分解，得到图3的结果，分解步骤和前面一致。

图38点DFT的2点分解

最后对2点DFT进一步分解得到最终的8点FFT信号计算流图：

图48点DFT的全分解

从图2到图4的过程中关于旋转系数的变化规律需要说明一下。

看起来似乎向前推一级，在奇数分组部分的旋转系数因子增量似乎就要变大，其实不是这样。

事实上奇数分组部分的旋转因子指数每次增量固定为1，只是因为每向前推进一次，该分组序列的数据个数变少了，为了统一使用以原数据N为基的旋转因子就进行了变换导致的。

每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。

以上边的8点DFT为例，第一次分组N=8,第二次分组N为4，为了统一根据式（4）进行了变换将N变为了8，但指数相应的需要乘以2。

N点基-2FFT算法的计算量

从图4可以看到N点DFT的FFT变换可以转为log2（N）级级联的蝶形运算，每一级均包含有N/2次蝶形计算。

而每一个蝶形运算包含了1次复数乘法，2次复数加法。

因此N点FFT计算的总计算量为：

复数乘法——N/2×log2（N）   复数加法——N×log2（N）。

假设被采样的序列为实数序列，那么也只有第一级的计算为实数与复数的混合计算，经过一次迭代后来的计算均变为复数计算，在这一点上和直接的DFT计算不一致。

因此对于输入序列是复数还是实数对FFT算法的效率影响较小。

一次复数乘法包含了4次实数乘法，2次实数加法，一次复数加法包含了2次复数加法。

因此对于N点的FFT计算需要总共的实数乘法数量为：

2×N×log2（N）；总的复数加法次数为：

2xNxlog2（N）。

N点基-2FFT算法的实现方法

从图4我们可以总结出对于点数为N=2^L的DFT快速计算方法的流程：

  1.对于输入数据序列进行倒位序变换。

该变换的目的是使输出能够得到X（0）~X（N-1）的顺序序列，同样以8点DFT为例，该变换将顺序输入序列x（0）~x（7）变为如图4的x（0）,x（4）,x

（2）,x（6）,x

（1）,x（5）,x（3）,x（7）序列。

其实现方法是：

假设顺序输入序列一次村在A（0）~A（N-1）的数组元素中，首先我们将数组下标进行二进制化（例：

对于点数为8的序列只需要LOG2（8）=3位二进制序列表示，序号6就表示为110）。

二进制化以后就是将二进制序列进行倒位，倒位的过程就是将原序列从右到左书写一次构成新的序列，例如序号为6的二进制表示为110，倒位后变为了011，即使十进制的3。

第三步就是将倒位前和倒位后的序号对应的数据互换。

依然以序号6为例，其互换过程如下：

temp=A（6）; A（6）=A（3）; A（3）=A（6）;

实际上考虑到执行效率，如果对于每一次输入的数据都需要这个处理过程是非常浪费时间的。

我们可以采用指向指针的指针来实现该过程，或者是采用指针数组来实现该过程。

2.蝶形运算的循环结构。

 从图4中我们可以看到对于点数为N=2^L的fft运算，可以分解为L阶蝶形图级联，每一阶蝶形图内又分为M个蝶形组，每个蝶形组内包含K个蝶形。

根据这一点我们就可以构造三阶循环来实现蝶形运算。

编程过程需要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。

3.浮点到定点转换需要注意的关键问题 

上边的分析都是基于浮点运算来得到的结论，事实上大多数嵌入式系统对浮点运算支持甚微，因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点方式。

对于简单的DFT运算从浮点到定点显得非常容易。

根据式

（1），假设输入x（n）是经过AD采样的数字序列，AD位数为12位，则输入信号范围为0~4096。

为了进行定点运算我们将旋转因子实部虚部同时扩大2^12倍，取整数部分代表旋转因子。

之后，我们可以按照

（1）式计算，得到的结果与原结果成比例关系，新的结果比原结果的2^12倍。

但是，对于使用蝶形运算的fft我们不能采用这种简单的放大旋转因子转为整数计算的方式。

因为fft是一个非对称迭代过程，假设我们对旋转因子进行了放大，根据蝶形流图我们可以发现其最终的结果是，不同的输入被放大了不同的倍数，对于第一个输入x（0）永远也不会放大。

举一个更加形象的例子，还是以图4为例。

从图中可以看出右侧的X（0）可以直接用下式表示：

从上式我们可以看到不同输入项所乘的旋转因子个数（注意这里是个数，就算是wn^0，也被考虑进去了，因为在没有放大时wn^0等于1，放大后所有旋转因子指数模均不为1，因此需要考虑）。

这就导致输入不平衡，运算结果不正确。

经查阅相关资料，比较妥善的做法是，首先对所有旋转因子都放大2^Q倍，Q必须要大于等于L，以保证不同旋转因子的差异化。

旋转因子放大，为了保证其模为1，在每一次蝶形运算的乘积运算中我们需要将结果右移Q位来抵消这个放大，从而得到正确的结果。

之所以采用放大倍数必须是2的整数次幂的原因也在于此，我们之后可以通过简单的右移位运算将之前的放大抵消，而右移位又代替了除法运算，大大节省了时间。

   4.计算过程中的溢出问题 

最后需要注意的一个问题就是计算过程中的溢出问题。

在实际应用中，AD虽然有12位的位宽，但是采样得到的信号可能较小，例如可能在0~8之间波动，也就是说实际可能只有3位的情况。

这种情况下为了在计算过程中不丢失信息，一般都需要先将输入数据左移P位进行放大处理，数据放大可能会导致溢出，从而使计算错误，而溢出的极限情况是这样：

假设我们数据位宽为D位（不包括符号位），AD采样位数B位，数字放大倍数P位，旋转因此放大倍数Q位，FFT级联运算带来的最大累加倍数L位。

我们得到:

假设AD位宽12，数据位宽32，符号位1位，因此有效位宽31位，采样点数N，那么我们可以得到log2（N）+P+Q<=19,假设点数128,又Q>=L可以得到放大倍数P<=5。

FFT代码示例 

根据以上的分析，我整理了一份128点的FFT代码如下，该代码在keil中仿真运行，未发现问题。

#defineN128

#definePOWER6//该值代表了输入数据首先被放大的倍数，放大倍数为2^POWER

#defineP_COEF8//该值代表了旋转因子被放大的倍数，放大倍数为2^POWER

#if（N==4）

#defineL2//L的定义满足L=log2（N）

#elif（N==8）

#defineL3//L的定义满足L=log2（N）

#elif（N==16）

#defineL4//L的定义满足L=log2（N）

#elif（N==32）

#defineL5//L的定义满足L=log2（N）

#elif（N==64）

#defineL6//L的定义满足L=log2（N）

#elif（N==128）

#defineL7//L的定义满足L=log2（N）

#endif

//AD采样位数为12位，本可以采用s16x[N]定义数据空间，但是为了节省存储空间，fft结果也将存储在该变量空间。

由于受

//fft影响变量需要重新定义，定义的方式及具体原因如下：

//fft过程会乘以较大旋转因子，因此需要32位定义

//fft过程会产生负数，因此需要有符号定义

//fft过程会出现复数，因此需要定义二维数组，[][0]存放实部，[][1]存放虚部

s32x[N][2]={};

//定义*p[N]是为了在第一次指针初始化以后，该数组指针按照位倒序指向数据变量x

//p[i][0]代表了指向数据的实部，p[i][1]代表了指向数据的虚部

s32*p[N];

//定义旋转因子矩阵

//旋转因子矩阵存储了wn^0,wn^1,wn^2...wn^（N/2-1）这N/2个复数旋转因子

s16w[N>>1][2]={256,0,256,-13,255,-25,253,-38,251,-50,248,-62,245,-74,241,-86,237,-98,231,-109,226,

                 -121,220,-132,213,-142,206,-152,198,-162,190,-172,181,-181,172,-190,162,-198,152,

                 -206,142,-213,132,-220,121,-226,109,-231,98,-237,86,-241,74,-245,62,-248,50,-251,38,

                 -253,25,-255,13,-256,0,-256,-13,-256,-25,-255,-38,-253,-50,-251,-62,-248,-74,-245,-86,

                 -241,-98,-237,-109,-231,-121,-226,-132,-220,-142,-213,-152,-206,-162,-198,-172,-190,-181,

                 -181,-190,-172,-198,-162,-206,-152,-213,-142,-220,-132,-226,-121,-231,-109,-237,-98,-241,

                 -86,-245,-74,-248,-62,-251,-50,-253,-38,-255,-25,-256,-13};

voidinit_pointer（void）

{

unsignedchari=0;

unsignedcharj=0;

unsignedchark=0;

for（i=0;i

{

j=0;

for（k=0;k

{j|=（（（i>>k）&0x01）<<（L-k-1））;

}

p[i]=&x[j][0];

}

}

/**description:

基2fft主函数,该函数将借助指针数组p对全局变量数组x进行fft计算

*并将结果存储在数组x中

*globalvar:

rw->x;r->p,w。

（r表示读，w表示写，rw表示读写）*/

voidfft2（void）

{

u8i;//i用于表示蝶形图级联的阶数

u8j;//表示蝶形分组起始点序列，蝶形分组跨度为2^i

u8k;//k表示蝶形组内偶数分支序列点号

u8gp_distance=1;//蝶形分组跨度

u8temp;//temp用于记录当前组间距离的一半,同时也表示了同一碟形两分支间的距离

u8gp_hf=0;//记录当前组的中间点序号

u8delta=N;//旋转因子下标增量，本来下标初始值应该为N/2，但是。

。

s16*pw=&（w[0][0]）;

inttp_result[2];//用于临时存放旋转因子和奇数分组的乘积

//输入信号序列放大

for（i=0;i

{x[i][0]<<=POWER;

x[i][1]<<=POWER;

}for（i=0;i