基于 FPGA 实现高速 FFT 处理的方法 【技术领域】
本发明公开了一种应用 FPGA 芯片实现高速 FFT 处理的方法， 它涉及数字信号处理领域。 背景技术 DFT(Discrete Fourier Transfer) 及其快速算法 FFT 是数字信号处理领域的核 心组成部分， FFT 算法多种多样， 按数据抽取方式不同又可分为基 2、 基 4 等。各算法优缺点 视不同制约因素而不同。FFT 实现方法也多种多样， 可以用软件实现， 也可以用硬件实现。 用软件在 PC 机或工作站上实现则计算速度很慢， 一般多结合具体系统用硬件实现， 例如用 单片机或 DSP(DigitalSignal Process) 处理器实现， 但是速度仍然不能满足某些应用领 域的需求， 难以与快速的 A/D 器件匹配。在雷达、 电子对抗等应用领域追求运算速度， 对实 时性要求非常高。FPGA 是英文 Field Programmable Gate Array 的缩写， 即现场可编程门 阵列， 它是在 PAL(Programmable Array Logic)、 GAL(GateArray Logic)、 EPLD(Embedded Programmable Logic Device) 等可编程器件的基础上进一步发展的产物， 具有更高的集成 度、 更强的逻辑实现能力和更好的设计灵活性。FPGA 由许多独立的可编程逻辑模块组成， 用户可以通过编程将这些模块连接起来实现不同的设计， 它是作为专用集成电路 (ASIC， Application Specific Integrated Chip) 领域中的一种半定制电路而出现的， 既解决了 定制电路的不足， 又克服了原有可编程器件门电路数有限的缺点。近些年来， FPGA 技术发 展迅速。一方面， 各种大容量、 高性能、 低功耗的可编程逻辑器件不断推出。另一方面， 出现 了许多 FPGA 设计辅助工具， 这些工具大大提高了以 FPGA 为基础的新型集成电路的设计效 率， 使更低成本、 更短周期的复杂数字系统开发成为可能。
数字信号处理中， 设输入信号 X(n) 是长为 N 的复序列， 其 DFT(DiscreteFourier Transfer) 定义为 ：
其 IDFT(Inverse Discrete Fourier Transfer) 定义为 ：其中，式 (1) 称为离散傅里叶正变换， 式 (2) 称为离散傅里叶反变换，x(n) 与 X(k) 构成了离散傅里叶变换对。
根据上述公式， 计算一个 X(k)， 需要 N 次复数乘法和 N-1 次复数加法， 而计算全部 2 X(k)(0 ≤ k ≤ N-1)， 共需要 N 次复数乘法和 N(N-1) 次复数加法， 直接计算全部 X(k) 共需 2 2 要 4N 次实数乘法和 2N(2N-1) 次实数加法。可见工作量与 N 成正比， N 越大， 运算量将显 著增加。
参考文献 ：
[1]Cooley J W， Tukey.An algorithm for the machine calculationof complexFourier series[J].Math.Compute， 1965， Vol.19， No.4， p297 ～ 301
[2] 任淑艳 . 应用 VHDL 语言的 FFT 算法实现 [J]， 哈尔滨理工大学学报， Vol.8， No.6， 2003.12， P.24-26
[3] 刘 朝 晖， 韩 月 秋 . 用 FPGA 实 现 FFT 的 研 究 [J]. 北 京 理 工 大 学 学 报 [J]， Vol.19， No.2， 1999.4， P.234-238
[4] 王旭东， 刘渝 . 全并行结构 FFT 及其 FPGA 实现 [J]. 南京航空航天大学学报， Vol.38 No.1 Feb.2006， p96-100
[5]Lo Sing Cheng， Ali Miri， Tet Hin Yeap.Efficient FPGAimplementation of FFT based Multipliers[C]， Proc.IEEE CCECE/CCGEI， Saskatoon， May 2005， pp.1300-1303
[6]Ayan Banerjee， Anindya Sundar Dhar， Swapna Banerjee.FPGArealization of a CORDIC based FFT processor for biomedical signalprocessing[J] ， Microprocessor and Microsystems 2001， 25(1)， pp.131-142
[7]PEREZ-PASCUAL ，A. ，SANSALONI ，T. ，and VALLS ，J.On-lineradix-2 butterflies on FPGA[C].Proc.WSES International Conferenceon Speech， signal and image processing， 2002， pp.2201-2205 [8]PEREZ-PASCUAL ， A. ， SANSALONI ， T. ， and VALLS ， J.FPGA basedradix-4 butterflies for HiperLAN2[C].Proc.IEEE InternationalSymposium on Circuits and Systems(ISCAS 2002)， 2002， pp.III.277-III.280
[9]Z.Szadkowski.16-point Discrete Fourier transform based onthe Radix-2 FFT algorithm implemented into Cyclone FPGA as the UHERCtrigger for horizontal air showers[C]， Proc.of the 29th ICRC， Pune， 2005.
[10]Wen-Chang Yeh， Chein-Wei Jen.High-Speed and Low-PowerSplit-Radix FFT[J].IEEE Tans.Signal Processing， 2003， Vol.51， No.3， p864-874.
[11]Yutai Ma， Lars Wanhammar.A Hardware Efficient Control ofMemory Addressing for High-Performance FFT Proces sors[J].IEEE Tans.Signal Processing， 2000， Vol.48， No.3， p917-920.
发明内容 技术问题 ： 本发明目的是为减少运算量， 提高运算速度， 针对背景技术存在的缺陷 提供一种基于 FPGA 实现高速 FFT 处理的方法。
技术方案 ： 本发明为实现上述目的， 采用如下技术方案 ：
本发明基于 FPGA 实现高速 FFT 处理的方法， 其特征在于离散傅里叶变换 X(k) 与 输入信号 x(n) 构成了离散傅里叶变换对， 都是长为 N 的复序列， 序列长度 N 是 2 的整数幂 M 次方， 即N＝2， 其中 M 为正整数， 表示幂次数， 其中 k、 n 分别表示离散傅里叶变换 X(k) 与 输入信号 x(n) 的序列 ；
首先将离散傅里叶变换 x(n) 分解为二组， 偶数项为一组， 奇数项为一组， 得到两 个 N/2 点的子序列， 即：
x1(r) ＝ x(2r)， x2(r) ＝ x(2r+1)， (0 ≤ r ≤ N/2-1) (3)
利用式 (4) 可以写成 ：其中 X1(k) 和 X2(k) 分别为 x1(r) 和 x2(r) 的 DFT。
有益效果 ：
Altera 公司采用传统结构的 FFT 算法其 32 点 FFT IP 核的运算时间大于 1.0us。 用 DSP 做的 32 点 FFT 时间也要 1.0us 以上。本发明利用 FPGA 器件丰富的逻辑、 RAM(Read Access Memory)、 ROM(Read Only Memory) 和 DSP 等资源， 采用一种新结构 FFT 算法， 使得 FFT 运算速度较传统方法有了很大提高。 随着芯片集成度的不断提高， 这种高速 FFT 处理方 法其优越性将越来越明显。而且该方法易扩展。
附图说明
图1： 蝶形运算图 ； 图2： 按时间抽取 ： 将一个 8 点 DFT 分解为 4 个 2 点 DFT ； 图3： 基 2 时间抽取全并行 FFT 算法结构图 ； 图4： 一步交换后 FFT 结构图 ； 图5： 高速并行 FFT 处理结构 ； 图6： 高速 FFT 的 FPGA 实现框图 ； 图7： 仿真时序图 ； 图8： FPGA 运算结果 ； 图9： 理论输出结果。具体实施方式
本发明采用一种基 -2 新型结构 FFT 算法， 应用 Altera 公司先进 FPGA 器件 Stratix 系列 EP1S25F780-5 来做硬件仿真。 Stratix 器件是一款采用高性能结构体系的 FPGA 器件。 它结合了强大内核性能， 大存储带宽， 数字信号处理 (DSP) 功能， 高速 I/O 性能和模块化设 计与一体的 FPGA。其内嵌的 DSP 模块具有很高的乘法运算速度， 在用 VHDL 编程时可以用 MegaWizard 方法指定用器件内部 DSP 模块生成乘法器， 用这种乘法器来做蝶形， 用多个蝶 形来构成 FFT 运算级， 通过流水操作可实现 FFT 运算并行化。本发明还将采用 Altera 公司 Quartus II 软件做逻辑设计和时序分析， 并利用 Quartus II 软件强大的硬件仿真和逻辑分 析功能， 对 VHDL 硬件描述语言进行设计编译、 综合、 布局、 布线， 最后下载到 FPGA 芯片中生 成了专用 FFT 处理芯片。
计算 DFT 过程中需要完成的运算的系数里， 存在相当多的对称性。通过研究这种 对称性， 可以简化计算过程中的运算， 从而减少计算 DFT 所需的时间。
利用的周期性、 对称性、 可约性三大特性， 可将 x(n) 或 X(k) 序列按一定规律分解成短序列进行运算， 这样可以避免大量的重复运算， 提高计算 DFT 的运算速度。算法形式可以分为两大类， 即按时间抽取 FFT 算法和按频率抽取 FFT 算法。
本专利采用基 -2 时间抽取 FFT 算法。 设序列长度 N 是 2 的整数幂次方， 即 N ＝ 2M， 其中 M 为正整数， 表示幂次数。
首先将序列 x(n) 分解为二组， 偶数项为一组， 奇数项为一组， 得到两个 N/2 点的子 序列， 即：
x1(r) ＝ x(2r)， x2(r) ＝ x(2r+1)， (0 ≤ r ≤ N/2-1) (3)
其中 X1(k) 和 X2(k) 分别为 x1(r) 和 x2(r) 的 DFT。利用式 (4) 可以写成 ：用图形化的方式表示如图 1 所示。每个蝶形运算需要一次复数相乘和两次复数加 法运算。采用这种表示方法， 上述分解运算的过程流图如图 2 所示。通过分解后， 每个 N/2 2 2 点 DFT 需要 N /4 次复数相乘， 两个 N/2 点 DFT 共需 N /2 次复数相乘， 组合运算共需 N/2 个 2 2 蝶形运算， 需 N/2 次复数相乘， 因此共需 N /2+N/2 ≈ N /2 次复数相乘， 与直接运算相比节省 近一半的运算量。
由于 N/2 ＝ 2M-1 依然为整数， 可将该序列一直分解下去， 直到最后是 2 点的 DFT 为 止。对于一个 8 点的 FFT， 根据上述算法可以得到一个完整的 N ＝ 8 的基 -2FFT 的运算流 图， 如图 3 所示。
根据上述 FFT 算法原理及图 3， 可以归纳出基 -2FFT 算法的一些规律和特点 ：
1、 整个 FFT 流图全由蝶形运算组成， 因此蝶形运算是 FFT 运算的核心， 该算法 的具体实现， 就是如何按一定顺序依次计算完成全部蝶形运算。每个蝶形运算需要一 次复数的相乘和两次复数的相加。对于一个 N ＝ 2M 的序列， 可以逐步分解到最后全为
2 点的 DFT 运算。全部 N 点的 FFT 运算共有个蝶形运算， 共需复数乘法次， 复数加法 af ＝ NM ＝ Nlog2N 次， 与直接计算相比运算量显著减少。 2、 流图中各蝶形的输入输出量是互相不相重的， 任何一个蝶形的两个输入量经蝶 形运算后， 便失去了价值， 不需要保存， 因此可以实现同址运算， 即经过一级运算后的结果 可以存放在原来贮存输入数据同一地址的单元中。因此， 只需要 N 个复数的存储单元， 既可 存放输入的原始数据， 又可存放中间结果， 而且还可以存放最后的运算结果。 这种原位运算 方式节省了大量的存储单元， 是 FFT 算法的一大优点， 从而有效降低了设备成本。
3、 对于同址运算结构， 运算完毕后输出结果 X(k) 仍按自然顺序存放。而输入序列 x(n)， 由于逐次按偶、 奇时间顺序抽取的分解结果， 重新排列了序列数据的存放顺序， 因此 它是按码位颠倒的顺序排放的。所谓码位颠倒， 就是将二进制的最高有效位到最低有效位 的位序颠倒放置而得的二进制数。例如， N ＝ 8 时， n ＝ 1 的二进制码位 001， 码位颠倒后为 100， 即相应的二进制数为 4。 实际运算中， 不直接将 x(n) 按照倒位序输出， 通常先将其按照 自然顺序输入存储单元， 并根据 FFT 算法， 然后进行变址运算变换得到倒位序的排列， 如表 1。
4、 蝶形运算所需系数各级有所不同。 每级自上而下观察， 均是以开始， 按等比级数依次递增， 周期重复。例如， 第 m 级运算， 系数为1 ＝ 0， 1，…， 2m-1-1， 共 2m-1 个系数， 指数 1 逐次增 1， 周期重复 2M-m 次。计算时所需系数可以事前计算好后存在一个数表 中， 这样运算速度快， 但需要开销内存。 亦可以在需要时依次递推计算， 这样可节省内存， 但 要增加一定的运算工作量。
表 1 码位排序 (N ＝ 8)
图 3 所示这种全并行结构特点， 如果采用同步时序电路， 则可以实现每个时钟节 拍输出一组 FFT 计算结果， 从而充分发挥了并行加流水结构的快速处理特点， 使得 FFT 运算 速度得到极大提高。
【高速 FFT 处理算法】
由前面的介绍知， 传统基 2、 8 点 FFT 算法如图 3 所示。箭头上数字代表旋转因子
中的 k。输入按码位颠倒的顺序排放， 输出是自然顺序。这种结构的特点是每个蝶形输 出数据仍然放在原来输入数据存储单元内， 这样只需要 2N 个存储单元 (FFT 中数据是复数 形式， 每点需要两个单元存储， N 是 FFT 运算点数 )。其缺点是不同级同一位置蝶形输入数 据寻址不固定， 难以实现循环控制。用 FPGA 编程时难以并行实现， 数据处理速度慢。当 FFT 点数增加时更是如此。
通过观察传统结构 FFT 算法可以发现， 如果将第一级中间两个蝶形交换， 则可以 得到如图 4 所示结构。对此结构进一步变换， 将第二级输出不送回原处而是将其存储起来 并按顺序存放， 则第三级中间两个蝶形跟着调换， 并把输入按顺序排列， 则可构成如图 5 所 示高速并行 FFT 处理结构。在蝶形变换同时， 其旋转因子亦跟随变换。
这种结构输入是顺序的， 而输出是位码倒序的， 每级旋转因子都被存放在 FPGA 片 内 ROM 里， 其中旋转因子的寻址是关键。由以上 FFT 结构图变换过程可推导每级旋转因子 寻址规律， 对于 N 点的 FFT(N ＝ 2k， k 为级数 )， 旋转因子有 个， 将它们按位码倒序形式排成一个含有 N/2 个元素的数组计为 ：
i共 N/2则第 i 级 (i ＝ 0， 1， 2…， k-1) 的旋转因子排列顺序是将 w(0)， w(1)， w(2)，… .， w(2 ) 重复 2k-i-1 次得到的。此结构特点是每 级输入、 输出数据顺序不变， 因此每级几何结构固定， 易于用 FPGA 实现并行 FFT 硬件结构， 从而明显加快 FFT 运算速度。
【FPGA 硬件实现】
FPGA 器件的特点是可用硬件描述语言对其进行灵活编程。利用 FPGA 厂商提供的软件可仿真硬件的功能， 使硬件设计如同软件设计一样灵活方便。 缩短了系统研发周期。 利 用 JTAG(Joint Test Action Group) 接口可对其进行 ISP(In System Programmable， 在系 统编程 ) 提高了系统的灵活性。随着芯片集成度的提高， 单片 FPGA 内不仅拥有大量的逻辑 单元而且还能集成 RAM、 ROM、 I/O 及 DSP 块等， 从而使 SOC(System On a Chip， 片上系统 ) 成为现实。本发明采用的是 Altera 公司的 Stratix 系列芯片 EP1S25F780-5， 用 QuartusII 软件做硬件仿真和逻辑分析， 并将输出结果与 Matlab 理论运算结果进行了比较。系统结构 框图如图 6 所示。本系统的结构特点是 ：
1. 为 提 高 数 据 精 度， 系 统 全 部 用 16 位 宽。 用 data-array， write-array 和 fly-array 三个数组实现了内核的并行处理， 可在 10 个时钟周期内算完 32 点复 FFT。时钟 周期为 25 纳秒， 因此 32 点 FFT 只需 250 纳秒。
2. 实现了数据的流水输入输出。在计算第 i 组数据同时， 第 i-1 组数据 FFT 结果 正在串行输出， 第 i+1 组数据则正在串行输入。因为内核计算是并行的、 速度快， 输入带宽 大。 本系统数据输入速率可达 200MHz。 为验证设计的正确性， 下面对其进行仿真， 输入信号 如下 ：
x(n) ＝ sin(2*n*pi/4.7)+sin(2*n*pi/16.3)+w(n) (n ＝ 0 ， 1， 2… 31) (6) 仿真数据由两个正弦波信号叠加噪声构成， 信噪比为 3dB， 正弦信号频率分别为 0.213fs 与 0.061fs， fs ＝ 1 为归一化采样频率， w(n) 是均值为 0 方差为 1 的高斯噪声。系 统仿真波形如图 7 所示。
将 FPGA 输出的 FFT 运算结果和用 Matlab 计算的 FFT 理论结果进行比较， 如图 8、 9 所示。 仿真信号采用的是实信号， 其频谱具有对称性， 图中只画了 32 点仿真结果的一半即 前 16 点输出结果。