本系列的2-7篇分别介绍了FIR和IIR滤波器的FPGA实现。除了数字滤波器外,快速傅里叶变换(FFT)也是DSP系统常用的运算单元,用于对信号进行频域分析。FFT算法的实现很复杂,但Altera和Xilinx都提供了可快速上手使用的IP核。本文将先介绍如何使用Quartus的FFT IP核进行频谱分析。
IP核概述
由于版本的关系,Quartus提供的IP核有两种,老版本集成在“MegaWizard Plug-In Manager”中;新版本集成在 “IP Catalog”和qsys中:
两种Quartus版本下的IP核,从使用者的角度来看仅仅是配置界面不同,在参数设置和使用方法上基本一致。前面在讲述FIR IP核的使用时以MegaWizard Plug-In Manager为例;本文以新版本的qsys为例。
Quartus的FFT IP核属于收费IP,如果是个人学习使用需要对IP核单独破解。破解方法与当时破解FIR Compiler IP核时相同,只不过FFT IP核对应的编号为0034(参考https://blog.csdn.net/fpgadesigner/article/details/80605869 中的破解步骤)。保存后重新打开Quartus,在License Setup页面看到FFT/IFFT的破解信息,则可以正常使用。
IP核参数设置
Quartus中FFT IP核的参数配置很简单(相比之下Vivado的FFT IP核就更复杂、更专业),配置界面如下:
1.Transform:“Length”设置FFT的长度即点数。“Direction”设置该IP核工作于FFT模式(Forward)或者IFFT(Reverse)模式,也可以设置为Bi-directional,IP核接口会多出一个名为reverse的信号,用于在代码中改变工作模式。
2.I/O:“Data Flow”其实设置的是FFT IP核的实现结构,包括4种:缓存(Buffered Burst)、突发(Burst)、流(Streaming)和可变流(Variable Streaming),计算速度和消耗资源依次增加。“Input Order”和“Output Order”表示输入数据顺序和输出数据顺序,有顺序nature和位逆序bit reverse。
3.Data and Twiddle:“Representation”设置FFT计算过程以定点(Fixed Point)、块浮点(Block Floating Point)或单精度浮点(Single Floating Point)表示,计算精度和消耗资源依次增加。块浮点是一种不错的折中方案,每个值的尾数各自表示,同时共用一个指数,有着不错的精度且能大大节省资源。
“Data Input Width”和“Twiddle Width”设置输入数据位宽和旋转因子位宽,系统会根据前面的设置计算出FFT输出数据的位宽“Data Output Width”。
4.Latency Estimates:系统会根据设置得到FFT计算所需的两个延时时间。
把标签卡切换到Advanced,配置界面如下:
这里是设置FFT引擎,结构“Architecture”可以设置为单输出(Single Output)或四输出(Quad Output),这里的“单”和“四”指的是内部FFT蝶形处理器的吞吐量,前者时分复用一个复数乘法器,单个时钟内得到1个输出;后者同时使用4个复数乘法器,单个时钟内得到4个输出。显然单输出消耗资源少,四输出计算速度快。“Number of Parallel Engines”是设置并行工作的FFT引擎数。
注意:集成在qsys中的IP核在配置完成后不会自动加入到工程中,需要在File管理界面手动添加qip文件到工程。
IP核接口说明
Quartus的很多IP核采用的是Avalon-ST接口,主要有数据(data)、准备好(ready)、有效(valid)和错误(error)等信号。这里主要介绍按照上述IP核配置(Burst模式、块浮点运算)所需的接口信号:
注意FFT IP核的数据总线皆以带符号数二进制补码表示。在设置为其它模式时,还会用到其它的接口。其它接口在后面的设计中使用到FFT的其它模式时,再做介绍。
FPGA设计
IP核的接口在Verilog HDL中进行设计时,一定要参考官方文档中给出的时序图。在IP核的配置界面点击“details”,可以找到IP核的user guide和release notes。
FFT IP核的相关文档如上图所示。注意,在altera卖身给intel后,官网进行了大改,目前老版本的quartus中的IP核点击“documentation”时已经成了error,即官网已经没有了相关文档。
在User Guide中找到Burst模式的时序图如下(注意信号名和波形之间有一定错位,不过这确实是原图):
Burst模式的接口时序相对比较简单,Verilog HDL示例代码如下所示:
`timescale 1ns / 1ps
//-------------------------------------------------------------------------
// 调用Quartus FFT IP核进行傅里叶变换
//-------------------------------------------------------------------------
module altera_fft(
input clk, //时钟
input rst_n, //低电平有效复位
input signed [11:0] data_in, //AD采集数据输入
output signed [24:0] amp //FFT计算结果
);
/********** 定义FFT IP核使用端口 **********/
wire inverse; // 输入,为1时进行IFFT,为0时进行FFT
wire sink_ready; // 输出,FFT引擎准备好接收数据时该信号置位
wire source_ready; // 输入,下传流模块在可以接收数据时将该信号置位
reg sink_valid; // 输入,有效标记信号,sink_valid和sink_ready都置位时开始数据传输
reg sink_sop; // 输入,高电平表示1帧数据载入开始
reg sink_eop; // 输入,高电平表示1帧数据载入结束
wire signed [11:0]sink_imag; // 输入,输入数据的虚部,二进制补码数据
wire [1:0] sink_error; // 输入,表示载入数据状态,一般置0
wire [1:0] source_error; // 输出,表示FFT转换出现的错误
wire source_sop; // 输出,高电平表示一帧数据转换开始
wire source_eop; // 输出,高电平表示一帧数据转换结束
wire [5:0]source_exp;
wire source_valid;
wire signed [11:0] xkre; // 输出,输出数据的实部,二进制补码数据
wire signed [11:0] xkim; // 输出,输出数据的虚部
assign sink_error = 2'b00;
assign source_ready = 1'b1; // 该信号置1表示永远准备好接收FFT数据
assign inverse = 1'b0; // 进行FFT正变换
assign sink_imag = 12'd0; // 输入数据虚部接地
/********** 控制FFT数据的载入 **********/
//在sink_valid为高电平期间,通过sink_sop、sink_eop控制载入数据
//设置FFT变换起始脉冲,sink_eop/sink_sop高电平后开始载入数据
//由于Burst模式下,FFT变换时延不超过2048个时钟周期,因此每2048个周期进行一次FFT变换
reg [10:0] count;
always @ (posedge clk or negedge rst_n)
if (!rst_n) begin
sink_eop <= 'b0;
sink_sop <= 'b0;
sink_valid <= 'b0;
count <= 'b0;
end
else begin
count <= count + 1'd1;
if (count == 1) sink_sop <= 1'b1; //计数1,置位sop,开始载入AD数据
else sink_sop <= 1'b0;
if (count == 512) sink_eop <= 1'b1; //计数512,置位eop,结束载入AD数据,进行512点FFT
else sink_eop <= 1'b0;
if (count>=1 & count<=512) sink_valid <= 1'b1; //载入数据期间,置位sink_valid
else sink_valid <= 1'b0;
end
/********** 调用IP核进行FFT变换,Burst模式 **********/
//FFT核,实现512点的FFT正变换,在sink_valid\sink_sop\sink_eop的控制下,每2048个数据进行一次正变换
fft3 u0 (
.clk (clk),
.reset_n (rst_n),
.sink_valid (sink_valid),
.sink_ready (sink_ready),
.sink_error (sink_error),
.sink_sop (sink_sop),
.sink_eop (sink_eop),
.sink_real (data_in),
.sink_imag (sink_imag),
.inverse (inverse),
.source_valid (source_valid),
.source_ready (source_ready),
.source_error (source_error),
.source_sop (source_sop),
.source_eop (source_eop),
.source_exp (source_exp),
.source_real (xkre),
.source_imag (xkim)
);
/********** 计算频谱的幅值信号 **********/
wire signed [23:0] xkre_square, xkim_square;
assign xkre_square = xkre * xkre;
assign xkim_square = xkim * xkim;
assign amp = xkre_square + xkim_square;
endmodule
程序中通过sink_valid、sink_sop、sink_eop的组合来控制ADC采样数据的载入,选取模2048计数器的512个时钟周期输入数据,其它时间等待计算结果。计数器的模值一定要大于IP核配置界面估计的“throughput latency”。输入数据只有ADC的采样数据作为实部,虚部设置为0。
Burst模式是一种计算比较慢的结构,必须在前一帧的FFT计算结果完全输出后,才能加载下一帧需要计算FFT的数据。上述代码中没有对sink_ready信号做判断是因为计数器模值较大,可以保证下一帧数据载入前上一帧的FFT结果已经全部输出。
FFT计算结果处理
波形数据经过FFT之后得到的是频域中的一串复数,我们需要从这串数据中提取出想要的信息。
1. 如果对结果不做处理,直接取复数的模,得到幅度谱(Magnitude);
2. 把第一个和最后一个对应直流分量的点的模值除以N,其余点的模值除N/2,得到幅值谱(Amplitude);
3. 工程中通常取对数坐标,转化为dB单位。
计算公式如下:
不过对于FPGA而言,开平方根sqrt和对数log都是极其消耗资源的运算(尤其是后者),通常可以将平方和数据发送到其它平台(PC、单片机)处理。
如果需要在FPGA中计算,Quartus也提供了相应的IP核。开平方可以使用ALTSQRT IP核实现;对数需要将其转换为浮点数,使用ALTFT_LOG IP核实现。这两种运算也可以自己使用查找表、拉格朗日逼近等方法设计,牺牲一定的精度来节省资源,相关设计将在本系列后续的文章中给出。
根据理论知识可知,N点FFT对应的频谱范围为0~Fs,Fs为采样频率,即每一个点对应的频率分辨率Δf=Fs/N。根据采样定理可知,FFT计算结果只有0~Fs/2是有效的,另一半是前者的对称。据此我们可以计算出每一个点在频谱中所对应的频率。
测试
本文对FFT IP核的使用测试以实例的方式给出。ADC采集正弦波信号,FFT计算之后将幅度谱发送至HMI串口屏上显示:
可以看到单频信号在幅度谱中对应一个频率点。50MHz ADC采集方波信号,FFT计算之后将幅度谱发送至PC上,用MATLAB显示:
有效频谱为0~25MHz,清楚可见方波信号为一系列谐波的叠加。具体工程示例及下载会在“FPGA综合系统设计(五):频谱分析系统“中给出。