用Xilinx Vivado HLS实现浮点复数QRD矩阵分解

作者：
Xilinx公司资深DSP专家 王宏强
Xilinx公司高级FAE 徐坚 邓涛
Xilinx公司工具与方法学高级专家 徐天容

在数字信号处理领域，如自适应滤波、DPD系数计算、MIMO Decoder 等，常常需要矩阵解方程运算以获得其系数，因此需对矩阵进行求逆运算。然而，由于直接对矩阵求逆会导致庞大的运算量，所以在实际工程中往往需要先将矩阵分解成几个特殊矩阵（正规正交矩阵或上、下三角矩阵以求其逆矩阵需要更小的运算量）的乘积。目前，QRD矩阵分解法是求一般矩阵全部特征值的最有效且广泛应用的方法之一。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R，称为QRD矩阵分解。

由于浮点具有更大的数据动态范围，所以在众多多算法中具有只需要一种数据类型的优势，所以很多QRD矩阵分解是基于浮点数据类型的。不过在通信应用中，更多的场景还是复数类型。因此，随着通信技术的日益发展，算法的复杂度越来越高，QRD矩阵的维度也越来越大。如果是用传统的手写RTL，浮点复数超大维度QRD矩阵分解的FPGA实现将变得非常复杂，需要很长的时间来编写RTL代码、仿真和进行验证等工作，使得开发效率不是很高。

本文将介绍如何使用Xilinx Vivado HLS（Vivado 高层次综合）工具实现浮点复数QRD矩阵分解并提升开发效率。使用Vivado HLS可以快速、高效地基于FPGA实现各种矩阵分解算法，降低开发者对算法FPGA的实现难度。其中包括：

• 使用Vivado HLS开发效率比手写RTL实现快5-10倍，而实现的FPGA资源效率与手写RTL接近
• 由于C/C++仿真验证比传统FPGA RTL要快100倍，Vivado HLS实现可以大大缩短用户的代码开发时间和仿真验证时间，从而大幅提高生产效率。
  

Vivado HLS工具介绍

Vivado HLS 是Xilinx公司2010 年收购 AutoESL 以后重新打造的高层次综合工具，它可以让用户通过添加适当的directives（制导语句） 和 constrains（约束）, 将其C/C++/System C代码直接转换成 FPGA RTL（ Verilog, VHDL, System C ）代码。让用户可以在算法开发环境而非通常的硬件开发环境中只需专注于算法规格和算法的 C 实现，Vivado HLS 工具会自动考虑 FPGA 微观实现架构，并可生成可综合的FPGA RTL代码。如图1所示。


图1 Vivado HLS设计介绍

Vivado HLS FPGA设计流程：

    首先用 C/C++/System C 将算法实现，并编写 C testbench 验证 C 的功能，确保其功能正确；
    然后就可以通过 Vivado HLS 工具进行 C 综合,将 C 转换成 RTL；
    接着做 C/RTL 的协同仿真（Co-simulation），以保证生成的RTL代码功能与C的功能完全一致。

最后 Vivado HLS 生成的 RTL 代码可直接用于 Xilinx 设计开发环境下做系统集成、仿真和生成bit文件。
如图2所示。


图2 Vivado HLS设计流程

QRD矩阵分解算法

QRD矩阵分解是将一个矩阵A分解成Q与R相乘,
A=QR
    其中R是个上三角矩阵， Q是个正交矩阵，
Q^TQ = 1
Q^T 是Q的转置共轭矩阵

QRD矩阵分解 Vivado HLS 实现C++代码构架

QRD矩阵分解C++实现代码的顶层模块是qrd_engine.cpp， 它调用cal_core.cpp（核心计算函数）, coef_cal.cpp（系数计算函数），以及浮点加、减、乘、除法等子函数来实现。

struct cf_t {
    float re;
    float im;
};
void qrd_engine
    (
    cf_t    in_u[(R_DIM+Y_DIM)/div_NUM][div_NUM],
    cf_t    pd_err_in,
    float   lamda,
    float   lamda_sqrt,
    float   diag[R_DIM],
    cf_t    r[R_DIM][X_DIM],
    cf_t    p[R_DIM]
    )
    {
#pragma HLS ARRAY_PARTITION variable=in_u complete dim=2
#pragma HLS ARRAY_PARTITION variable=r complete dim=2
//注：这里对数组加入完全分割directive，目的是提高数据的并行带宽，从而获得并行计算。
#pragma HLS RESOURCE variable=in_u core=RAM3S
//注：这里对数组加入RAM3S directive，目的是控制生成的BRAM的延迟。
for (i=0;i<R_DIM;i++)   
{
#pragma HLS PIPELINE II =124
#pragma HLS latency max = 123
//注：这里加入II（interval是指下一个新的可以输入QRD模块的数据与前一个数据之间所间隔的时钟周期数）及latency directives (延迟制导语句)，目的是控制这个模块在指定的延迟节拍范围内完成所有的计算。
coef_calc(lamda_sqrt,lamda,
,pre_in_u,pd_err_in,&s_o,&s_conj_o,&lamda_sqrtxs_o,&c_o,&lamda_sqrtxc_o,&diag_out,&p_o,&pd_err);
   
cal_core(u_tmp, r_tmp, s_n, i, j, k, c_o, lamda_sqrtxc_o, lamda_sqrtxs_o, s_conj_o, &in_u_w2, &r【i][r_addr]);
}
}
void coef_calc
(
    float    lamda_sqrt,
    float    lamda,
    float    r_diag,
    cf_t    u_diag,
    cf_t    pd_err_in,
    cf_t    *s,
    cf_t    *s_conj,
    cf_t    *lamda_sqrtxs,
    float    *c,
    float    *lamda_sqrtxc,
    float    *diag,
    cf_t    *p_o,
    cf_t    *pd_err
)

void calc_core
(
cf_t in_u,
cf_t r,
int s_n,
int i,
unsigned char j,
unsigned char k,
float c_o,
float lamda_sqrtxc_o,
cf_t lamda_sqrtxs_o,
cf_t s_conj_o,
cf_t* u_ret,
cf_t* r_ret
)