Tag: FFT原理

FFT原理与实现

在 数字信号处理中常常需要用到离散傅立叶变换(DFT),以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征,但是算法计算量大,耗时长,不 利于计算机实时对信号进行处理。因此至DFT被发现以来,在很长的一段时间内都不能被应用到实际的工程项目中,直到一种快速的离散傅立叶计算方法—— FFT,被发现,离散是傅立叶变换才在实际的工程中得到广泛应用。需要强调的是,FFT并不是一种新的频域特征获取方式,而是DFT的一种快速实现算法。 本文就FFT的原理以及具体实现过程进行详尽讲解。

其 中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个频率点的相对幅度。

从 图2到图4的过程中关于旋转系数的变化规律需要说明一下。看起来似乎向前推一级,在奇数分组部分的旋转系数因子增量似乎就要变大,其实不是这样。事实上奇 数分组部分的旋转因子指数每次增量固定为1,只是因为每向前推进一次,该分组序列的数据个数变少了,为了统一使用以原数据N为基的旋转因子就进行了变换导 致的。每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。以上边的8点DFT为例,第一次分组N=8,第二次分组N为4,为了统一根据式 (4)进行了变换将N变为了8,但指数相应的需要乘以2。

从 图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)。

该 变换的目的是使输出能够得到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为例,其互换过程如下:

从图4中我们可以看到对于点数为N = 2^L的fft运算,可以分解为L阶蝶形图级联,每一阶蝶形图内又分为M个蝶形组,每个蝶形组内包含K个蝶形。根据这一点我们就可以构造三阶循环来实现蝶 形运算。编程过程需要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。

上边的分析都是基于浮点运算来得到的结论,事实上大多数嵌入式系统对浮点运算支持甚微,因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点 方式。对于简单的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的整数次 幂的原因也在于此,我们之后可以通过简单的右移位运算将之前的放大抵消,而右移位又代替了除法运算,大大节省了时间。

最后需要注意的一个问题就是计算过程中的溢出问题。在实际应用中,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快速傅立叶变换的工作原理

FFT是计算DFT的快速算法,但是它是基于复数的,所以计算实数DFT的时候需要将其转换为复数的格式,下图展示了实数DFT和虚数DFT的情况,实数DFT将时域中N点信号转换成2个(N/2+1)点的频域信号,其中1个(N/2+1)点的信号称之为实部,另一个(N/2+1)点的信号称之为虚部,实部和虚部分别是正弦和余弦信号的幅度。

相比较而言,复数DFT将2个N点的时域信号转换为2个N点的频域信号。时域和频域中,1个N点信号是实部,另1个N点信号是虚部。

如果要计算N点实数DFT,则将这个N个点作为时域中的实部,另取N个0点作为时域的虚部,用FFT计算这样一个复数信号的DFT得到2个N点的频域信号,1个N点是实部另1个N点是虚部,在这两个N点的信号中,从0到N/2个点就是须计算的N点实数的DFT频域。

对于实数DFT来说,就像前几章讲的那样,它的频域也是离散周期信号,其周期为N点,从0到N/2点和1-N到-1点具有对称性,这个你可以从下面一张图看出。图中坐标不是用N表示是用采样频率的分数表示,如果你看不懂,请看前面几章。

FFT的计算可以分为三步:首先将1个N点的时域信号分成N个1点的时域信号,然后计算这N个1点时域信号的频域,得到N个频域的点,然后将这个N个频域的点按照一定的顺序加起来,就得到了我们需要的频谱。这里每个点的意思是复数,都有实部和虚部。

这一步很简单,因为一个时域的点的频谱的数值就是它自己,所以这一步什么也不需做,但需明白这时候N个点不是时域信号了,而是频域信号。

这一步是最麻烦的一步。就是和前面时域分解的顺序相反,将2个1点的频域信号变成1个2点的频域信号,再将2个2点的频域信号变成1个4点的频域信号,一直到结束。这里看下如何将2个4点的频域信号变成1个8点的频域信号。

首先对1个4点的频域信号进行复制,这样能稀释时域信号,也对另1个4点的频域信号进行复制不过复制之前需要乘上正弦函数,这样得到的稀释时域信号时经过了平移的,然后将这两个频域信号加起来,如下图所示。之所以这么做的目的是在时域分解的时候就是用这种交织的分解方式的。

比如使用基4或者基8,这样不是2点一计算,而是4点或者8点一计算,可以提高速度。

FFT对DSP来说就像是晶体管对电子学来说,都是领域的基础,每个人都知道怎么使用它们,但是只有很少一部分真正了解它们的原理。