基于FFT谱分析测频算法的FPGA实现

田西柱

摘 要：在信号频谱分析试验中，通过FPGA实现FFT。在MAX+plusⅡ系统环境下，介绍了流水线结构FFT的蝶形单元设计，详解了旋转因子的生成，通过地址产生单元和块浮点单元实现了运算结果的输出，并将其输出结果与Matlab结果进行比较。

关键词：FPGA；QuartusⅡ；FFT处理器；旋转因子

中图分类号：TP391 文献标志码：A 文章编号：2095-1302（2014）10-00-03

0 引 言

在FPGA实验中，主要是用FPGA来实现FFT，使其完成对信号的频谱分析。流水线结构FFT的设计主要是蝶形单元的设计，通过旋转参数的生成，将运算结果写入地址并完成输出。

1 实验原理及步骤

1.1 QuartusⅡ开发环境

QuartusⅡ是Altera公司提供的FPGA/CPLD集成开发软件，在QuartusⅡ上可以完成设计输入、HDL综合、布新布局（适配）、仿真和选择以及硬件测试等流程，它提供了一种与结构无关的设计环境，使设计者能方便地进行设计输入、开始处理和器件编程。QuartusⅡ具备仿真功能，同时支持第三方的仿真工具（如ModelSin）。此外，QuartusⅡ与Matlab和DSP Builder结合，可用进行基于FPAG的DSP系统开发，是DSP硬件系统实现的工具EDA工具。

QuartusⅡ设计与开发的流程如图1所示。

图1 QuartusII设计与开发的流程

1.2 快速傅里叶变换（FFT）

FFT是DFT的快速算法。设离散的有限长时间序列为 x（n），0≤n≤N-1，则其离散的傅里叶变换为：

（1）

k=0，1，…，N-1 WN=e-j（2nN） （2）

其中：x（n）为时域点；X（k）为频域点；WN为旋转因子；x（n）、X（k）、WN都是复数。完成整个DFT运算共需要N2次复数乘法以及N（N-1）次复数加法运算。当N很大时，运算量很大，对于实时信号处理，要求CPU运算速度很高，难以工程实现。因此，出现了快速傅里叶变换（FFT）算法。FFT算法的基本思想是利用旋转因子WN的周期性、对称性、特殊性以及周期N的可互换性，将长度为N点序列的DFT运算逐次分级为较短序列的DFT运算，并将相同项合并，因为DFT的运算量与N2成正比，当N减小时，就大大减少了运算量，提高了运算效率。N=2n个点的DFT复数乘法量由N2次降为（N/2）log2N次，复数加法由N（N-1）次降为（N/2）log2N次。

FFT算法种类很多，基本上可分为两大类：一类是针对N等于2的整数次幂的算法，如基2算法、基4算法和分裂基算法等；另一类是针对N不等于2的整数次幂的算法，以Winograd为代表，它们有重要的理论价值，但是不适于硬件实现。基2算法结构简单，但运算量大。基4算法相对于基2算法更为复杂，但是计算量减少了。FFT算法按分解方式的不同又可以分为时域抽取算法（decimation in time，DIT）和频域抽取算法（decimation in frequency，DIF）两种。这两种算法在本质上都是一种基于标号分解的算法，在运算量和复杂性等方面完全一样。考虑到本设计FFT运算的点数不是太多，故选用了时域抽取基2算法（DIT）。

1.3 按时间抽取的基2-FFT算法（DIT-基2-FFT）原理

FFT运算的基本单元是蝶形运算单元，基2蝶形运算符号如图2所示。设蝶形运算输入数据为A=Ap+Aqj，B=Bp+Bqj，旋转因子W rN=Wp+Wqj。则蝶型运算输出为：

（3）

图2 基2蝶形运算

FFT算法由多级蝶形运算构成，具体运算流图也有多种形式。本设计选用了输入倒序、输出顺序的运算流图，图3所示为N=8点时的DIT-FFT运算流图。这种运算流图是同址运算，其优点是：在同一级运算中， 每个蝶形的两个输入数据只对计算本蝶形有用，而且蝶形的输入输出数据节点又同在一条水平线上，这就意味着计算完一个蝶形运算后，所得输出数据可以立即存入元出入数据所占用的存储器。因此，在硬件实现时可以节省存储单元。

图3 N=8点时的DIT-FFT运算流图

一个长度为N的序列x（n），满足N=2M，M为整数。那么此序列x（n）的FFT运算流图由M级蝶形运算构成，每一级有N/2个蝶形运算，第L级蝶形运算中使用旋转因子的个数为2L，L=0，1，2，…，M-1。64点FFT运算，分6级，每级有32个蝶形运算。

1.4 FFT处理器结构设计

FFT算法的FPGA硬件实现在Altera公司的MAX+plusⅡ系统环境下开发完成，选用基于查找表结构内嵌存储器的APEX20系列FPGA器件。图4为FFT处理器的结构图。本设计采用单元结构设计思路，整个处理器由数据接收单元、运算单元、旋转因子存储单元、地址产生单元和中央控制单元5个单元组成，各单元在中央控制单元的控制下协调工作。其中，内部接收单元采用乒乓RAM结构，扩大了数据吞吐量，计算单元采用流水与并行结合的结构，加快了运算速。

图4 FFT处理器的结构图

1.5 中央控制单元

中央控制单元是整个系统的控制核心，其主要功能是控制数据流向，协调各单元之间的运行。中央控制单元根据系统时钟确定当前蝶型运算所处的级数m和个数n，并把m、n传送给地址产生单元。地址产生单元产生蝶型运算两个输入数据和旋转因子的地址，并把地址传送给运算RAM和旋转因子存储器。在中央控制单元读使能信号控制下两个输入数据和旋转因子被读出。读出的数据进行必要的延迟和定标处理后，送给运算单元。经过蝶型运算后，运算结果按原址写入RAM。

1.6 数据接收单元

数据接收单元主要功能是按帧接收外部数据，并将每帧数据按码位倒置的顺序乒乓存入接收RAM1或接收RAM2。中央控制单元交替的对接收RAM中的数据进行处理，当中央控制单元将接收RAM1中的数据取出，经过蝶型运算，结果存入运算RAM1的同时上一帧数据的FFT运算结果从运算RAM2取出。接收RAM用FPGA的片上双口RAM实现，接收单元控制写端口，中心控制单元控制读端口。

1.7 运算单元

运算单元由蝶型运算器和运算RAM组成。蝶型运算器完成对输入数据的蝶型运算，运算RAM作为FFT的中间数据缓存。蝶型运算器输入数据为A=Ap+Aqj，B=Bp+Bqj，旋转因子W rN=Wp+Wqj，蝶型运算输出如式（3）所示。根据式（3），蝶型运算器可由一个复数乘法和两个复数加（减）法器组成。为了提高运算速度采用并行运算，用四个实数乘法器、三个实数加法器、三个实数减法器组成。蝶型运算器实现框图如图5所示。蝶型运算各个模块利用MAX+plusⅡ开发软件中所提供的宏单元生成。

图5 蝶形运算器

运算RAM1和运算RAM2作为FFT的中间数据缓存。两块RAM交替作为数据读出和运算结果写入单元，直到第6级蝶型运算完成。

1.8 旋转因子存储单元

旋转因子存储单元，存储FFT运算所需要的旋转因子WrN，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。旋转因子先在Matlab中分实部和虚部产生，转化成16位定点数，并将结果保存成hex文件格式。利用MAX+plusII软件提供的ROM宏模块 “lpm_rom”产生两个（N/2）×16 b的ROM，并分别用旋转因子实部和虚部对应的hex文件对两个ROM初始化，这样旋转因子的值就固化在了FPGA中。对应不同级的蝶型运算，地址产生器产生相应的地址送给ROM将旋转因子读出。

1.9 地址产生单元

地址产生单元产生蝶型运算两个输入数据和旋转因子的地址。实现的方法是根据地址产生的算法，通过逻辑运算产生。前面已介绍本设计FFT实现结构为同址运算，即蝶型运算的结果仍然写回输入数据读出单元。因此，将读数据地址延迟若干时钟周期，就可作为运算结果写入地址。对于N点FFT运算，用m∈（0∧log2N-1），n∈（0∧N/2-1）表示第m级的第n个蝶型运算。addr_A，addr_B（addr_A

（4）

式中：（n/2m）×2m+1描述为将n的低m位清零，再左移一位。n%2m描述为取n的低m位。蝶型运算所对应的旋转因子存储器入口地址设为addr_w，对于N点FFT共需要N/2个旋转因子，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。根据第m级蝶型运算所需旋转因子的排列规律，旋转因子存储器入口地址应为：

addr_w=（n×2（log2N-1-m））（N/2） （5）

64点FFT旋转因子有32个，共需要5位表示地址。第m级、第n个蝶型运算的旋转因子地址可由逻辑移位的方法快速得到。将n（5位）左移位（5-m），低位补0，得到了在（0～31）中的addr_w。

1.10 块浮点单元

块浮点单元的实现思路是每级蝶型运算结果动态扩展但最大扩展2位。块浮点单元对蝶型运算结果的高3位进行检测，判断当前结果动态范围扩展位数，记录当前级的最大扩展位数。下一级蝶型运算时，根据前一级的最大扩展位数，对读出的数据进行定标，选取数据送入蝶型运算器。块浮点单元将每一级运算结果动态范围扩展位数进行累加，和FFT运算结果一同输出。

1.11 FFT处理器功能仿真与设计验证

仿真结果如图6所示：

由仿真结果可以看出，该FFT处理器为串行流水线结构，各级运算模块没有实现并行运行。

图6 时序仿真图

2 实验情况记录

为验证仿真结果的正确性，采用上述方法实现256点FFT处理器，同时为提高精度将输入数据的实部和虚部采用16位二进制数。对函数：

x=cos（200πt）+cos（400π）

以120 MHz的频率进行抽样，取256点作为FFT处理器输入数据进行快速傅里叶变换，并将其输出结果与Matlab结果进行比较，结果如图7～图9所示。

图7 抽样函数及Matlab计算频谱

图8 FPGA计算频谱

图9 各采样点相对误差

3 结 语

由于处理器采用定点运算，在进行乘法和加法运算时不可避免地造成一定误差，尤其是在功率谱接近零值的这些点上，相对误差较大，但是在我们更为关心的功率谱幅值点上，相对误差仅为1%上下，完全可以满足大多数应用对于运算精度的要求。

参考文献

[1]张辉. 张记龙.基-2 FFT处理器的FPGA实现[J].计算机与现代化.2009（9）：154-157.

[2]程佩青.数字信号处理教程[M].3版.北京：清华大学出版社，2007.

[3]袁俊泉.孙敏琪.曹瑞.Verilog HDL数字系统设计及其应用[M].西安：西安电子科技大学出版社，2002.

[4]谢彦林.可变点流水线结构FFT处理器的设计及其FPGA实现[D].西安：西安电子科技大学，2007.

[5]云霄.可配置FFT/IFFT处理器的设计及其FPGA构造[D].西安：西安电子科技大学，2009.

[6]蔡可红.基于FPGA的FFT设计与实现[D].南京：南京理工大学，2006.

1.6 数据接收单元

数据接收单元主要功能是按帧接收外部数据，并将每帧数据按码位倒置的顺序乒乓存入接收RAM1或接收RAM2。中央控制单元交替的对接收RAM中的数据进行处理，当中央控制单元将接收RAM1中的数据取出，经过蝶型运算，结果存入运算RAM1的同时上一帧数据的FFT运算结果从运算RAM2取出。接收RAM用FPGA的片上双口RAM实现，接收单元控制写端口，中心控制单元控制读端口。

1.7 运算单元

运算单元由蝶型运算器和运算RAM组成。蝶型运算器完成对输入数据的蝶型运算，运算RAM作为FFT的中间数据缓存。蝶型运算器输入数据为A=Ap+Aqj，B=Bp+Bqj，旋转因子W rN=Wp+Wqj，蝶型运算输出如式（3）所示。根据式（3），蝶型运算器可由一个复数乘法和两个复数加（减）法器组成。为了提高运算速度采用并行运算，用四个实数乘法器、三个实数加法器、三个实数减法器组成。蝶型运算器实现框图如图5所示。蝶型运算各个模块利用MAX+plusⅡ开发软件中所提供的宏单元生成。

图5 蝶形运算器

运算RAM1和运算RAM2作为FFT的中间数据缓存。两块RAM交替作为数据读出和运算结果写入单元，直到第6级蝶型运算完成。

1.8 旋转因子存储单元

旋转因子存储单元，存储FFT运算所需要的旋转因子WrN，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。旋转因子先在Matlab中分实部和虚部产生，转化成16位定点数，并将结果保存成hex文件格式。利用MAX+plusII软件提供的ROM宏模块 “lpm_rom”产生两个（N/2）×16 b的ROM，并分别用旋转因子实部和虚部对应的hex文件对两个ROM初始化，这样旋转因子的值就固化在了FPGA中。对应不同级的蝶型运算，地址产生器产生相应的地址送给ROM将旋转因子读出。

1.9 地址产生单元

地址产生单元产生蝶型运算两个输入数据和旋转因子的地址。实现的方法是根据地址产生的算法，通过逻辑运算产生。前面已介绍本设计FFT实现结构为同址运算，即蝶型运算的结果仍然写回输入数据读出单元。因此，将读数据地址延迟若干时钟周期，就可作为运算结果写入地址。对于N点FFT运算，用m∈（0∧log2N-1），n∈（0∧N/2-1）表示第m级的第n个蝶型运算。addr_A，addr_B（addr_A

（4）

式中：（n/2m）×2m+1描述为将n的低m位清零，再左移一位。n%2m描述为取n的低m位。蝶型运算所对应的旋转因子存储器入口地址设为addr_w，对于N点FFT共需要N/2个旋转因子，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。根据第m级蝶型运算所需旋转因子的排列规律，旋转因子存储器入口地址应为：

addr_w=（n×2（log2N-1-m））（N/2） （5）

64点FFT旋转因子有32个，共需要5位表示地址。第m级、第n个蝶型运算的旋转因子地址可由逻辑移位的方法快速得到。将n（5位）左移位（5-m），低位补0，得到了在（0～31）中的addr_w。

1.10 块浮点单元

块浮点单元的实现思路是每级蝶型运算结果动态扩展但最大扩展2位。块浮点单元对蝶型运算结果的高3位进行检测，判断当前结果动态范围扩展位数，记录当前级的最大扩展位数。下一级蝶型运算时，根据前一级的最大扩展位数，对读出的数据进行定标，选取数据送入蝶型运算器。块浮点单元将每一级运算结果动态范围扩展位数进行累加，和FFT运算结果一同输出。

1.11 FFT处理器功能仿真与设计验证

仿真结果如图6所示：

由仿真结果可以看出，该FFT处理器为串行流水线结构，各级运算模块没有实现并行运行。

图6 时序仿真图

2 实验情况记录

为验证仿真结果的正确性，采用上述方法实现256点FFT处理器，同时为提高精度将输入数据的实部和虚部采用16位二进制数。对函数：

x=cos（200πt）+cos（400π）

以120 MHz的频率进行抽样，取256点作为FFT处理器输入数据进行快速傅里叶变换，并将其输出结果与Matlab结果进行比较，结果如图7～图9所示。

图7 抽样函数及Matlab计算频谱

图8 FPGA计算频谱

图9 各采样点相对误差

3 结 语

由于处理器采用定点运算，在进行乘法和加法运算时不可避免地造成一定误差，尤其是在功率谱接近零值的这些点上，相对误差较大，但是在我们更为关心的功率谱幅值点上，相对误差仅为1%上下，完全可以满足大多数应用对于运算精度的要求。

参考文献

[1]张辉. 张记龙.基-2 FFT处理器的FPGA实现[J].计算机与现代化.2009（9）：154-157.

[2]程佩青.数字信号处理教程[M].3版.北京：清华大学出版社，2007.

[3]袁俊泉.孙敏琪.曹瑞.Verilog HDL数字系统设计及其应用[M].西安：西安电子科技大学出版社，2002.

[4]谢彦林.可变点流水线结构FFT处理器的设计及其FPGA实现[D].西安：西安电子科技大学，2007.

[5]云霄.可配置FFT/IFFT处理器的设计及其FPGA构造[D].西安：西安电子科技大学，2009.

[6]蔡可红.基于FPGA的FFT设计与实现[D].南京：南京理工大学，2006.

1.6 数据接收单元

数据接收单元主要功能是按帧接收外部数据，并将每帧数据按码位倒置的顺序乒乓存入接收RAM1或接收RAM2。中央控制单元交替的对接收RAM中的数据进行处理，当中央控制单元将接收RAM1中的数据取出，经过蝶型运算，结果存入运算RAM1的同时上一帧数据的FFT运算结果从运算RAM2取出。接收RAM用FPGA的片上双口RAM实现，接收单元控制写端口，中心控制单元控制读端口。

1.7 运算单元

运算单元由蝶型运算器和运算RAM组成。蝶型运算器完成对输入数据的蝶型运算，运算RAM作为FFT的中间数据缓存。蝶型运算器输入数据为A=Ap+Aqj，B=Bp+Bqj，旋转因子W rN=Wp+Wqj，蝶型运算输出如式（3）所示。根据式（3），蝶型运算器可由一个复数乘法和两个复数加（减）法器组成。为了提高运算速度采用并行运算，用四个实数乘法器、三个实数加法器、三个实数减法器组成。蝶型运算器实现框图如图5所示。蝶型运算各个模块利用MAX+plusⅡ开发软件中所提供的宏单元生成。

图5 蝶形运算器

运算RAM1和运算RAM2作为FFT的中间数据缓存。两块RAM交替作为数据读出和运算结果写入单元，直到第6级蝶型运算完成。

1.8 旋转因子存储单元

旋转因子存储单元，存储FFT运算所需要的旋转因子WrN，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。旋转因子先在Matlab中分实部和虚部产生，转化成16位定点数，并将结果保存成hex文件格式。利用MAX+plusII软件提供的ROM宏模块 “lpm_rom”产生两个（N/2）×16 b的ROM，并分别用旋转因子实部和虚部对应的hex文件对两个ROM初始化，这样旋转因子的值就固化在了FPGA中。对应不同级的蝶型运算，地址产生器产生相应的地址送给ROM将旋转因子读出。

1.9 地址产生单元

地址产生单元产生蝶型运算两个输入数据和旋转因子的地址。实现的方法是根据地址产生的算法，通过逻辑运算产生。前面已介绍本设计FFT实现结构为同址运算，即蝶型运算的结果仍然写回输入数据读出单元。因此，将读数据地址延迟若干时钟周期，就可作为运算结果写入地址。对于N点FFT运算，用m∈（0∧log2N-1），n∈（0∧N/2-1）表示第m级的第n个蝶型运算。addr_A，addr_B（addr_A

（4）

式中：（n/2m）×2m+1描述为将n的低m位清零，再左移一位。n%2m描述为取n的低m位。蝶型运算所对应的旋转因子存储器入口地址设为addr_w，对于N点FFT共需要N/2个旋转因子，W rN=e（-j2π/N）r （r=0，1，…，N/2-1）。根据第m级蝶型运算所需旋转因子的排列规律，旋转因子存储器入口地址应为：

addr_w=（n×2（log2N-1-m））（N/2） （5）

64点FFT旋转因子有32个，共需要5位表示地址。第m级、第n个蝶型运算的旋转因子地址可由逻辑移位的方法快速得到。将n（5位）左移位（5-m），低位补0，得到了在（0～31）中的addr_w。

1.10 块浮点单元

块浮点单元的实现思路是每级蝶型运算结果动态扩展但最大扩展2位。块浮点单元对蝶型运算结果的高3位进行检测，判断当前结果动态范围扩展位数，记录当前级的最大扩展位数。下一级蝶型运算时，根据前一级的最大扩展位数，对读出的数据进行定标，选取数据送入蝶型运算器。块浮点单元将每一级运算结果动态范围扩展位数进行累加，和FFT运算结果一同输出。

1.11 FFT处理器功能仿真与设计验证

仿真结果如图6所示：

由仿真结果可以看出，该FFT处理器为串行流水线结构，各级运算模块没有实现并行运行。

图6 时序仿真图

2 实验情况记录

为验证仿真结果的正确性，采用上述方法实现256点FFT处理器，同时为提高精度将输入数据的实部和虚部采用16位二进制数。对函数：

x=cos（200πt）+cos（400π）

以120 MHz的频率进行抽样，取256点作为FFT处理器输入数据进行快速傅里叶变换，并将其输出结果与Matlab结果进行比较，结果如图7～图9所示。

图7 抽样函数及Matlab计算频谱

图8 FPGA计算频谱

图9 各采样点相对误差

3 结 语

由于处理器采用定点运算，在进行乘法和加法运算时不可避免地造成一定误差，尤其是在功率谱接近零值的这些点上，相对误差较大，但是在我们更为关心的功率谱幅值点上，相对误差仅为1%上下，完全可以满足大多数应用对于运算精度的要求。

参考文献

[1]张辉. 张记龙.基-2 FFT处理器的FPGA实现[J].计算机与现代化.2009（9）：154-157.

[2]程佩青.数字信号处理教程[M].3版.北京：清华大学出版社，2007.

[3]袁俊泉.孙敏琪.曹瑞.Verilog HDL数字系统设计及其应用[M].西安：西安电子科技大学出版社，2002.

[4]谢彦林.可变点流水线结构FFT处理器的设计及其FPGA实现[D].西安：西安电子科技大学，2007.

[5]云霄.可配置FFT/IFFT处理器的设计及其FPGA构造[D].西安：西安电子科技大学，2009.

[6]蔡可红.基于FPGA的FFT设计与实现[D].南京：南京理工大学，2006.