FPGA数字信号处理（八）Quartus FFT IP核实现

本系列的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综合系统设计（五）：频谱分析系统“中给出。