用FPGA实现FFT算法_fpga 32点fft运算时间

(原文地址:http://www.ednchina.com/blog/chwb/3709/message.aspx)

引言

　　DFT(Discrete Fourier Transformation)是数字信号分析与处理如图形、语音及图像等领域的重要变换工具，直接计算DFT的计算量与变换区间长度N的平方成正比。当N较大时，因计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。快速傅立叶变换(Fast Fourier Transformation，简称FFT)使DFT运算效率提高1～2个数量级。其原因是当N较大时，对DFT进行了基4和基2分解运算。FFT算法除了必需的数据存储器ram和旋转因子rom外，仍需较复杂的运算和控制电路单元，即使现在,实现长点数的FFT仍然是很困难。本文提出的FFT实现算法是基于FPGA之上的，算法完成对一个序列的FFT计算，完全由脉冲触发，外部只输入一脉冲头和输入数据，便可以得到该脉冲头作为起始标志的N点FFT输出结果。由于使用了双ram，该算法是流型(Pipelined)的，可以连续计算N点复数输入FFT，即输入可以是分段N点连续复数数据流。采用DIF(Decimation In Frequency)-FFT和DIT(Decimation In Time)-FFT对于算法本身来说是无关紧要的，因为两种情况下只是存储器的读写地址有所变动而已，不影响算法的结构和流程，也不会对算法复杂度有何影响。算法实现的可以是基2/4混合基FFT，也可以是纯基4FFT和纯基2FFT运算。

傅立叶变换和逆变换

对于变换长度为N的序列x(n)其傅立叶变换可以表示如下：

　　　　　　　　　　　　　　　　　　　　　式（１）

其中，W="exp"(-2π/N)。 当点数N较大时，必须对式(1)进行基4/基2分解，以短点数实现长点数的变换。而IDFT的实现在DFT的基础上就显得较为简单了：

　　式（２）

由式(2)可以看出，在FFT运算模块的基础上，只需将输入序列进行取共轭后再进行FFT运算，输出结果再取一次共轭便实现了对输入序列的IDFT运算，因子1/N对于不同的数据表示格式具体实现时的处理方式是不一样的。IDFT在FFT的基础上输入和输出均有一次共轭操作，但它们共用一个内核，仍然是十分方便的。

基4和基2
基4和基2运算流图及信号之间的运算关系如图1所示：

　
以基4为例，令A="r0"+j×i0；B="r1"+j×i1；C="r2"+j×i2；D="r3"+j×i3；Wk0=c0+j×s0：Wk1=c1+j×s1；Wk2=c2+j×s2；Wk3=c3+j×s3。分别代入图1中的基4运算的四个等式中有：
A'=[r0+(r1×c1-i1×s1)+(r2×c2-i2×s2)+(r3×c3-i3×s3)]+j[i0+(i1×c1+r1×s1)+(i2×c2+r2×s2)+(i3×c3+r3×s3)] 式(3)
B'=[r0+(i1×c1+r1×s1)-(r2×c2-i2×s2)-(i3×c3+r3×s3)]+j[i0-(r1×c1-i1×s1)-(i2×c2+r2×s2)+(r3×c3-i3×s3)] 式(4)
C'=[r0-(r1×c1-i1×s1)+(r2×c2-i2×s2)-(r3×c3-i3×s3)]+j[i0-(i1×c1+r1×s1)+(i2×c2+r2×s2)-(i3×c3+r3×s3)] 式(5)
D'=[r0-(i1×c1+r1×s1)-(r2×c2-i2×s2)+(i3×c3+r3×s3)]+j[i0+(r1×c1-i1×s1)-(i2×c2+r2×s2)-(r3×c3-i3×s3)] 式(6)

　　可以看出，式(3)至式(6)有多个公共项和类似项，这一点得到充分利用之后可以大大缩减基4和基2运算模块中的乘法器的个数，如上面A'至D'的四个等式中的这三对类似项：(r1×c1-i1×s1)与(i1×c1+r1×s1)、(r2×c2-i2×s2)与(i2×c2+r2×s2)、(r3×c3-i3×s3)与(i3×c3+r3×s3)以高于输入数据率的时钟进行时分复用，最终可以做到只需要3个甚至1个复数乘法器便可以实现。基2运算之所以采用图1-(b)中的形式进行基2运算，是为了将基本模块做成基4/2复用模块，它对于N有着更大的适用性和可借鉴性。在基4、基2和基4/2模块的基础上，构建基16、基8和基16/8模块有着非常大的意义。

算法实现

　　傅立叶变换实现时首先进行基2、基4分解，一般来说，如果算法使用基4实现，虽然使用的资源多了一些，但速度上的好处足以弥补。如果资源充足，使用基16、基8或基16/8复用模块，速度可以大大提高。一般FFT实现简单框图如图2所示。

 

　　在图2中，运算模块即为基 2/4/8 /16模块或它们的复用模块，Rom表中存储的是N点旋转因子表。控制模块产生所有的控制信号，存储器1和2的读写地址、写使能、运算模块的启动信号及因子表的读地址等信号。当然对于运算模块为基16/8复用模块时，控制模块就需要产生模式选择信号，如对于运算模块是基4/2模块时，该信号就决定了内部运算模块是进行基4运算还是基2运算。存储器1作为当前输入标志对应输入N点数据的缓冲器，存储器2作为中间结果存储器，用于存储运算模块计算出的各Pass的结果。在图中的各种地址、使能和数据的紧密配合下，经过一定延时后输出计算结果及其对应指示标志。图2只是一定点或浮点的FFT实现模块，如果是块浮点运算，则必须加入一个数据因子控制器，控制每遍运算过程中的数据大小，并根据各个Pass的乘性因子之和的大小，对最终输出进行大小控制，以保证每段FFT运算输出增益一致。

　　外部输入为N点数据段流和启动信号(N点之间如无间隔，则每N数据点输入一脉冲信号)，一方面，外部数据存入存储器1中，同时通过控制模块的控制，读出存储器1中的前段N点数据和Rom表中的因子及相关控制信号送入运算核心模块进行各个Pass的运算，每个Pass的输出都存入存储器2中，最后一个Pass的计算结果存入存储器2中，并在下一个启动头到来后，输出计算结果。对图2的实现，除去运算模块，关键是各个Pass数据因子读写地址及控制信号的配合。

速度、资源和精度

　　假定输入数据的速率为fin，则每数据的持续时间T="1"/fin，运算模块的计算时钟频率为fa，对于N(N="2p"，p即为Pass数目)点FFT计算时延与Pass数目直接相关。如果使用基2运算不考虑控制开销，纯粹的计算时延为td="p"×N×T×fin/fa。显然在fa>p× fin时，在N点内可完成FFT运算。否则不能完成，即不能实现流型的变换。这在N很大且输入数据速率较高时以FPGA实现几乎是不可能的，而且内部计算时钟过高容易导致电路的工作不稳定。设基2时的最小可流型工作运算频率为fa0，则使用基4实现流型的变换，计算时钟fa= fa0就可以。而使用基8时计算时钟fa= fa0便可完成，基16时为fa0的1/4。上面所讨论的是纯基运算，当N不为4的幂次方时(如N="2048"=16×16×8，运算模块为基16/8复用模块)，而又希望使用较低倍的时钟完成运算时，图2中的运算模块必然包括基4/2复用模块(即基16/8复用模块)，这也就是前面提到复用模块的主要用意。由上面的分析可以得出结论，如果计算使用的基越大，完成速度越快。　　但是，使用基16/8模块所使用的逻辑资源要比基4/2模块多将近一倍，这是因为基16/8复用模块是以基4模块和基4/2复用模块构建而成。当然，可以直接实现基16/8复用模块，但用FPGA很难解决复杂度和成本问题。另外，如果流型运算间隔比N点数据长度长一倍以上，可以考虑在较低的计算时钟下使用基2运算模块实现流型FFT。

　　运算结果的精度直接与计算过程中数据和因子位数(浮点算法)相关，如果中间计算的位数、存储数据位数和Rom表中的位数越大，输出精度就越大。当然，位数增大后逻辑运算资源和存储资源都会直线上升。

浮点、块浮点和定点FFT
　　根据运算过程中对数据位数取位和表示形式的不同，可以将FFT分为浮点FFT、块浮点FFT和定点FFT。它们在实现时对于系统资源的要求是不同的，而且有着不同的适用范围。　浮点FFT是基于数据表示为浮点的基础之上的，即数据是由一纯小数和一因子组成，输入要转成纯小数和因子的浮点表示形式，所有计算过程中保存应得计算结果大小，而输出要变成所需大小的定点表示形式。只要因子位数足够大，浮点FFT计算是不会溢出的。而定点则是所有计算过程中都是定点运算，如果各个Pass的截位规则不适当，很容易出现溢出，必须要有溢出控制。块浮点是介于它们之间的一种运算机制，它是根据本Pass的输入数据的大小，在计算之前进行控制(数据上移一比特或下移一比特或乘以一特定因子)，可以保证不溢出，但一般也需要溢出控制。　　浮点运算没有溢出，信号平均信噪比高，但由于因子的运算必然导致电路复杂，实现困难。定点运算实现简单，难以保证不溢出，需要统计得出合适的截位规则，否则溢出严重导致输出结果错误。块浮点由于每个Pass(包括最后输出前)结束后有一统计控制过程，延时较大，但是可以保证不溢出而且电路又相对浮点来说简单得多。　　应根据具体应用的具体要求，选择合适的FFT。如果要求精度，并且要解决频域很高的单频干扰，就必须使用浮点的FFT，使用数据位数很大的定点和块浮点也能解决这个问题，但位数的确定十分困难。如果不要求高精度，逻辑资源和Rom比较紧张，可考虑定点运算。如果输入在频域集中于几个点上或者对精度要求一般，可以慢速处理，可以采用块浮点运算，就能够保证这几点的信噪比，而忽略其他点处的信噪比。

基于FPGA的快速傅立叶变换

 

摘要：在对FFT（快速傅立叶变换）算法进行研究的基础上，描述了用FPGA实现FFT的方法，并对其中的整体结构、蝶形单元及性能等进行了分析。

关键词：FPGA FFT

傅立叶变换是数字信号处理中的基本操作，广泛应用于表述及分析离散时域信号领域。但由于其运算量与变换点数Ｎ的平方成正比关系，因此，在Ｎ较大时，直接应用ＤＦＴ算法进行谱变换是不切合实际的。然而，快速傅立叶变换技术的出现使情况发生了根本性的变化。本文主要描述了采用ＦＰＧＡ来实现２ｋ／４ｋ／８ｋ点ＦＦＴ的设计方法。

１　整体结构

一般情况下，Ｎ点的傅立叶变换对为：

其中，ＷＮ＝ｅｘｐ(－２ ｐｉ／Ｎ)。Ｘ(ｋ)和ｘ(ｎ)都为复数。与之相对的快速傅立叶变换有很多种,如ＤＩＴ(时域抽取法)、ＤＩＦ（频域抽取法）、Ｃｏｏｌｅｙ－Ｔｕｋｅｙ和Ｗｉｎｏｇｒａｄ等。对于２ｎ傅立叶变换，Ｃｏｏｌｅｙ－Ｔｕｋｅｙ算法可导出ＤＩＴ和ＤＩＦ算法。本文运用的基本思想是Ｃｏｏｌｅｙ－Ｔｕｋｅｙ算法，即将高点数的傅立叶变换通过多重低点数傅立叶变换来实现。虽然ＤＩＴ与ＤＩＦ有差别，但由于它们在本质上都是一种基于标号分解的算法，故在运算量和算法复杂性等方面完全一样，而没有性能上的优劣之分，所以可以根据需要任取其中一种，本文主要以ＤＩＴ方法为对象来讨论。

Ｎ＝８１９２点ＤＦＴ的运算表达式为：

式中，ｍ＝(４ｎ１＋ｎ２)(２０４８ｋ１＋ｋ２)(ｎ＝４ｎ１＋ｎ２，ｋ＝２０４８ｋ１＋ｋ２)其中ｎ１和ｋ２可取０,１,．．．,２０４７,ｋ１和ｎ２可取０,１,２,３。

由式（３）可知，８ｋ傅立叶变换可由４×２ｋ的傅立叶变换构成。同理，４ｋ傅立叶变换可由２×２ｋ的傅立叶变换构成。而２ｋ傅立叶变换可由１２８×１６的傅立叶变换构成。１２８的傅立叶变换可进一步由１６×８的傅立叶变换构成，归根结底，整个傅立叶变换可由基２、基４的傅立叶变换构成。２ｋ的ＦＦＴ可以通过５个基４和１个基２变换来实现；４ｋ的ＦＦＴ变换可通过６个基４变换来实现；８ｋ的ＦＦＴ可以通过６个基４和１个基２变换来实现。也就是说：ＦＦＴ的基本结构可由基２／４模块、复数乘法器、存储单元和存储器控制模块构成，其整体结构如图１所示。

图１中，ＲＡＭ用来存储输入数据、运算过程中的中间结果以及运算完成后的数据，ＲＯＭ用来存储旋转因子表。蝶形运算单元即为基２／４模块，控制模块可用于产生控制时序及地址信号，以控制中间运算过程及最后输出结果。

２　蝶形运算器的实现

 

基４和基２的信号流如图２所示。图中，若Ａ＝ｒ０＋ｊ＊ｉ０，Ｂ＝ｒ１＋ｊ＊ｉ１，Ｃ＝ｒ２＋ｊ＊ｉ２，Ｄ＝ｒ３＋ｊ＊ｉ３是要进行变换的信号，Ｗｋ０＝ｃ０＋ｊ＊ｓ０＝１，Ｗｋ１＝ｃ１＋ｊ＊ｓ１，Ｗｋ２＝ｃ２＋ｊ＊ｓ２，Ｗｋ３＝ｃ３＋ｊ＊ｓ３为旋转因子，将其分别代入图２中的基４蝶形运算单元，则有：

Ａ′＝[ｒ０＋(ｒ１×ｃ１－ｉ１×ｓ１)＋(ｒ２×ｃ２－ｉ２×ｓ２)＋(ｒ３×ｃ３－ｉ３×ｓ３)]＋ｊ[ｉ０＋(ｉ１×ｃ１＋ｒ１×ｓ１)＋(ｉ２×ｃ２＋ｒ２×ｓ２)＋(ｉ３×ｃ３＋ｒ３×ｓ３)] 　（４）

Ｂ′＝[ｒ０＋(ｉ１×ｃ１＋ｒ１×ｓ１)－(ｒ２×ｃ２－ｉ２×ｓ２)－(ｉ３×ｃ３＋ｒ３×ｓ３)]＋ｊ[ｉ０－(ｒ１×ｃ１－ｉ１×ｓ１)－(ｉ２×ｃ２＋ｒ２×ｓ２)＋(ｒ３×ｃ３－ｉ３×ｓ３)] 　(５）

Ｃ′＝[ｒ０－(ｒ１×ｃ１－ｉ１×ｓ１)＋(ｒ２×ｃ２－ｉ２×ｓ２)－(ｒ３×ｃ３－ｉ３×ｓ３)]＋ｊ[ｉ０－(ｉ１×ｃ１＋ｒ１×ｓ１)＋(ｉ２×ｃ２＋ｒ２×ｓ２)－(ｉ３×ｃ３＋ｒ３×ｓ３)] （６）

Ｄ′＝[ｒ０－(ｉ１×ｃ１＋ｒ１×ｓ１)－(ｒ２×ｃ２－ｉ２×ｓ２)＋(ｉ３×ｃ３＋ｒ３×ｓ３)]＋ｊ[ｉ０＋(ｒ１×ｃ１－ｉ１×ｓ１)－(ｉ２×ｃ２＋ｒ２×ｓ２)－(ｒ３×ｃ３－ｉ３×ｓ３)] （７）

而在基２蝶形中，Ｗｋ０和Ｗｋ２的值均为１，这样，将Ａ，Ｂ，Ｃ和Ｄ的表达式代入图２中的基２运算的四个等式中，则有：

Ａ′＝ｒ０＋(ｒ１×ｃ１－ｉ１×ｓ１)＋ｊ[ｉ０＋(ｉ１×ｃ１＋ｒ１×ｓ１)] （８）

Ｂ′＝ｒ０－ (ｒ１×ｃ１－ｉ１×ｓ１)＋ｊ[ｉ０－(ｉ１×ｃ１＋ｒ１×ｓ１)] 　（９）

Ｃ′＝ｒ２＋(ｒ３×ｃ３－ｉ３×ｓ３)＋ｊ[ｉ０＋(ｉ３×ｃ３＋ｒ３×ｓ３)] （１０）

Ｄ′＝ｒ２－(ｒ３×ｃ３－ｉ３×ｓ３)＋ｊ[ｉ０－(ｉ３×ｃ３＋ｒ３×ｓ３)] （１１）

在上述式（４）～（１１）中有很多类同项，如ｉ１×ｃ１＋ｒ１×ｓ１和ｒ１×ｃ１－ｉ１×ｓ１等，它们仅仅是加减号的不同，其结构和运算均类似，这就为简化电路提供了可能。同时，在蝶形运算中，复数乘法可以由实数乘法以一定的格式来表示，这也为设计复数乘法器提供了一种实现的途径。

以基４为例，在其运算单元中，实际上只需做三个复数乘法运算，即只须计算ＢＷｋ１、ＣＷｋ２和ＤＷｋ３的值即可，这样在一个基４蝶形单元里面，最多只需要３个复数乘法器就可以了。在实际过程中，在不提高时钟频率下，只要将时序控制好便可利用流水线（Ｐｉｐｅｌｉｎｅ）技术并只用一个复数乘法器就可完成这三个复数乘法，大大节省了硬件资源。

图2 基2和基4蝶形算法的信号流图

３　ＦＦＴ的地址

ＦＦＴ变换后输出的结果通常为一特定的倒序,因此，几级变换后对地址的控制必须准确无误。

倒序的规律是和分解的方式密切相关的，以基８为例，其基本倒序规则如下：

基８可以用２×２×２三级基２变换来表示，则其输入顺序则可用二进制序列（ｎ１ｎ２ ｎ３）来表示，变换结束后，其顺序将变为（ｎ３ ｎ２ ｎ１），如：Ｘ０１１→ ｘ１１０，即输入顺序为３，输出时顺序变为６。

更进一步，对于基１６的变换，可由２×２×２×２，４×４，４×２×２等形式来构成，相对于不同的分解形式，往往会有不同的倒序方式。以４×４为例，其输入顺序可以用二进制序列（ｎ１ｎ２ ｎ３ ｎ４）来表示变换结束后，其顺序可变为（（ｎ３ ｎ４）（ｎ１ ｎ２）），如： Ｘ０１１１→ ｘ１１０１。即输入顺序为７，输出时顺序变为１３。

在２ｋ／４ｋ／８ｋ的傅立叶变换中，由于要经过多次的基４和基２运算，因此，从每次运算完成后到进入下一次运算前，应对运算的结果进行倒序，以保证运算的正确性。

４　旋转因子

Ｎ点傅立叶变换的旋转因子有着明显的周期性和对称性。其周期性表现为：

FFT之所以可使运算效率得到提高，就是利用

ＦＦＴ之所以可使运算效率得到提高，就是利用了对称性和周期性把长序列的ＤＦＴ逐级分解成几个序列的ＤＦＴ，并最终以短点数变换来实现长点数变换。

根据旋转因子的对称性和周期性，在利用ＲＯＭ存储旋转因子时，可以只存储旋转因子表的一部分，而在读出时增加读出地址及符号的控制，这样可以正确实现ＦＦＴ。因此,充分利用旋转因子的性质，可节省７０％以上存储单元。

实际上，由于旋转因子可分解为正、余弦函数的组合，故ＲＯＭ中存的值为正、余弦函数值的组合。对２ｋ／４ｋ／８ｋ的傅立叶变换来说，只是对一个周期进行不同的分割。由于８ｋ变换的旋转因子包括了２ｋ／４ｋ的所有因子，因此，实现时只要对读ＲＯＭ的地址进行控制，即可实现２ｋ／４ｋ／８ｋ变换的通用。

５　存储器的控制

因ＦＦＴ是为时序电路而设计的，因此，控制信号要包括时序的控制信号及存储器的读写地址，并产生各种辅助的指示信号。同时在计算模块的内部，为保证高速，所有的乘法器都须始终保持较高的利用率。这意味着在每一个时钟来临时都要向这些单元输入新的操作数，而这一切都需要控制信号的紧密配合。

为了实现ＦＦＴ的流形运算，在运算的同时，存储器也要接收数据。这可以采用乒乓ＲＡＭ的方法来完成。这种方式决定了实现ＦＦＴ运算的最大时间。对于４ｋ操作，其接收时间为４０９６个数据周期，这样ＦＦＴ的最大运算时间就是４０９６个数据周期。另外，由于输入数据是以一定的时钟为周期依次输入的，故在进行内部运算时，可以用较高的内部时钟进行运算，然后再存入ＲＡＭ依次输出。

为节省资源，可对存储数据ＲＡＭ采用原址读出原址写入的方法，即在进行下一级变换的同时，首先应将结果回写到读出数据的ＲＡＭ存贮器中；而对于ＲＯＭ，则应采用与运算的数据相对应的方法来读出存储器中旋转因子的值。

在２ｋ／４ｋ／８ｋ傅立叶变换中，要实现通用性，控制器是最主要的模块。２ｋ、４ｋ、８ｋ变换具有不同的内部运算时间和存储器地址，在设计中，针对不同的点数应设计不同的存储器存取地址，同时，在完成变换后，还要对开始输出有用信号的时刻进行指示。

６　硬件的选择

本设计的硬件实现选用的是现场可编程门阵列(ＦＰＧＡ)来满足较高速度的需要。本系统在设计时选用的是ＡＬＴＥＲＡ公司的ＳＴＲＡＴＩＸ芯片，该芯片中包含有ＤＳＰ单元，可以完成较为耗费资源的乘法器单元。同时，该器件也包含有大量存储单元，从而可保证旋转因子的精度。

除了一些专用引脚外，ＦＰＧＡ上几乎所有的引脚均可供用户使用，这使得ＦＰＧＡ信号处理方案具有非常好的Ｉ／Ｏ带宽。大量的Ｉ／Ｏ引脚和多块存储器可使设计获得优越的并行处理性能。其独立的存储块可作为输入／工作存储区和结果的缓存区，这使得Ｉ／Ｏ可与ＦＦＴ计算同时进行。在实现的时间方面，该设计能在４０９６个时钟周期内完成一个４０９６点的ＦＦＴ。若采用１０ＭＨｚ的输入时钟，其变换时间在２００μｓ左右。而由于最新的ＦＰＧＡ使用了ＭｕｌｔｉＴｒａｃｋ互连技术，故可在２５０ＭＨｚ以下频率稳定地工作，同时，ＦＦＴ的实现时间也可以大大缩小。

ＦＦＴ运算结果的精度与输入数据的位数及运算过程中的位数有关，同时和数据的表示形式也有很大关系。一般来说，浮点方式比定点方式精度高。而在定点计算中，存储器数据的位数越大，运算精度越高，使用的存储单元和逻辑单元也越多。在实际应用中，应根据实际情况折衷选择精度和资源。本设计通过ＭＡＴＬＡＢ进行仿真证明：其实现的变换结果与ＭＡＴＬＡＢ工具箱中的ＦＦＴ函数相比，信噪比可以达到６５ｄｂ以上，完全可以满足一般工程的实际应用要求。