diff --git a/00-Preface.md b/00-Preface.md index 287acc60..64dc852a 100755 --- a/00-Preface.md +++ b/00-Preface.md @@ -4,7 +4,7 @@ 本书将着重介绍如何对算法进行高层次综合（HLS），并以此完成一些比较具体、细分的FPGA应用。我们的目的是让读者认识到用HLS创造并优化硬件设计的好处。当然，FPGA的并行编程肯定是有别于在多核处理器、GPU上实行的并行编程，但是一些最关键的概念是相似的。例如，设计者必须充分理解内存层级和带宽、空间局部性与时间局部性、并行结构和计算与存储之间的取舍与平衡。 -本书将更多的作为一个实际应用的向导，为那些对于研发FPGA系统有兴趣的读者提供帮助。对于大学教育来说，这本书将更适用于高阶的本科课程或研究生课程，同时也对应用系统设计师和嵌入式程序员有所帮助。我们不会对C/C++方面的知识做过多的阐述，而会以提供很多的代码的方式作为示范。另外，读者需要对基本的计算机架构有所熟悉，例如流水线（pipeline），加速，阿姆达尔定律（Amdahl's Law）。以寄存器传输级（RTL)为基础FPGA设计知识并不是必需的，但会对理解本书有所帮助。 +本书将更多的作为一个实际应用的向导，为那些对于研发FPGA系统有兴趣的读者提供帮助。对于大学教育来说，这本书将更适用于高阶的本科课程或研究生课程，同时也对应用系统设计师和嵌入式程序员有所帮助。我们不会对C/C++方面的知识做过多的阐述，而会以提供很多的代码的方式作为示范。另外，读者需要对基本的计算机架构有所熟悉，例如流水线（pipeline），加速，阿姆达尔定律（Amdahl's Law）。以寄存器传输级（RTL）为基础FPGA设计知识并不是必需的，但会对理解本书有所帮助。 本书囊括了很多对教学很有帮助的内容。每个章节均有一些小问题留给读者，这些问题将有助于加深对于材料的理解。在加州大学圣地亚哥分校（UCSD）的CSE 237C这门课里也有很多用HLS开发的项目，如果有出于教育目的的需要，我们可以对提出申请的读者分享这些课程项目的文件。这些HLS项目主要是与数字信号分析相关，并重点研究无线交流系统的开发。每个单独的项目都或多或少与书中的某一章节有所关联。这些项目以赛灵思大学计划（Xilinx University Program）使用的FPGA开发板为基础而开发，设计基础参考 。赛灵思也同时提供这些开发板的商业订单。同时我们鼓励读者在 申请Vivado HLS的试用许可。 diff --git a/03-CORDIC.md b/03-CORDIC.md index 2db5a38c..30472477 100755 --- a/03-CORDIC.md +++ b/03-CORDIC.md @@ -23,7 +23,7 @@ CORDIC算法是1950年由Jack Volder发明，它最开始是作为数字解决 ​ CORDIC核心思想是在一个二维平面上高效地执行一组矢量旋转。在这些旋转运算的基础上增加一些简单控制，我们就可以实现各种基础操作，例如，三角函数，双曲函数，对数函数，实乘和复乘，以及矩阵分解和因式分解。CORDIC已经广泛应用于信号处理、机器人技术、通信和许多科学计算领域。由于CORDIC占用资源少，所以常用在FPGA设计中。 ​ 在下文中，我们将介绍CORDIC如何执行给定输入角θ的正弦和余弦的过程。这是使用一系列矢量旋转来完成的，这些简单的操作在硬件中使用非常有效。在高层次，算法使用一系列旋转来工作，目标是达到目标输入角θ。实现这种效率的关键创新是旋转可以以需要最少计算的方式完成。尤其我们使用乘以2的常数幂来执行旋转。这意味着简单地在硬件中移动位是非常有效的，因为它不需要任何逻辑。 - + ​ 图3.1提供了用于计算cosφ和sinφ的CORDIC程序的高级概述。在这种情况下，我们在x轴上开始初始旋转矢量，即0°角。然后，我们执行一系列迭代旋转;在这个例子中，我们只执行四次旋转，但通常这是40次旋转。每个后续旋转使用越来越小的角度，这意味着每次迭代都会为输出值增加更多精度。在每次迭代中，我们决定以较小的角度进行正向或负向旋转。我们旋转的角度值是先验固定的;因此，我们可以轻松地将它们的值存储在一个小内存中，并保持我们到目前为止已经旋转的累积角度的运行总和。如果该累积角度大于我们的目标角度φ，则我们执行负旋转。如果它更小，那么旋转就是正的。一旦我们完成了足够数量的旋转，我们就可以通过直接读取最终旋转矢量的x和y值来确定cosφ和sinφ。如果我们的最终向量的幅度为1，则x =cosφ且y =sinφ。 我们从一些术语开始，目的是重新定义你的一些基本的三角函数和矢量概念。 如果熟悉的话就不必看了。 但请记住，创建高效硬件设计最重要的一个方面是真正理解应用程序; 只有这样，设计师才能有效地利用优化指令并执行代码重构，这是获得最有效设计所必需的。 @@ -130,7 +130,6 @@ $$ 计算旋转操作数值，我们可以得到 - $$ x_i = x_{i-1} \cos 45^\circ - y_{i-1} \sin 45^\circ = x_{i-1} \cdot \sqrt 2/2 - y_{i-1} \cdot \sqrt 2/2 \quad(3.10) $$ @@ -348,7 +347,6 @@ cordic(THETA_TYPE theta, COS_SIN_TYPE &s, COS_SIN_TYPE &c) 图3.3:CORDIC代码实现计算给定角度的正弦和余弦值。 - ![图3.4:图中显示了一个二维平面和一个同时用笛卡尔形式(x,y)和极坐标形式(r,θ)表示的矢量,通过该矢量说明这两个坐标系之间的关系。](images/rotation.jpg) ​ 此代码接近于“软件”版本。它有多种方式来提高其性能并减小资源面积。我们将在本章后面讨论如何优化这段代码。 @@ -402,7 +400,6 @@ $$ ​ cordic函数使用当前常用的变量类型。例如，变量sigma被定义为int类型，其他变量使用自定义数据类型(例如，THETA_TYPE和COS_SIN_TYPE)。在许多情况下，HLS工具能够进一步优化这些值，来实现硬件资源优化。举个例子,在图3.3中,变量sigma的值为1或−1，即使已经将变量声明为至少包含32位位宽的int类型，在不改变程序功能的情况下，可以用更少的位数用来实现。在其它情况下，特别是变量经常出现在函数输入、存储和递归中，HLS工具不能针对变量进行自动优化。在这些情况下，修改代码使用较小的数据类型是避免资源浪费的重要优化手段。 - ​ 虽然减少变量位宽通常是一个好的优化方法，但是这种优化方法可能改变程序的行为。一个小位宽数据类型不能像大位宽数据类型那样表征大量信息，而且无限的二进制数可以表示所有无限精度的实数。幸运的是，作为设计人员，我们可以选择符合特定需求的精度，并且可以在精度、资源占用和性能之间进行权衡。 ​ 在进一步讨论cordic函数数值表示方法之前，我们首先介绍数值表示的背景知识。我们提供基础知识，因为这在理解Vivado HLS所提供的数据类型进行数值表示时非常重要。下一节将从数值表示的基本背景开始讲解，然后讨论Vivado HLS中可用的任意精度变量。 @@ -556,12 +553,14 @@ float p2 = 0xB4p-4;//Initialize p2 to "11.25" | ------ | ----- | ----- | -------- | -------- | -------- | -------- | --------------- | | 0 | 1 | 1 | 0 | 1 | 0 | 0 | =3.25 | -​ 同样，当一个数值不能通过给定小数位数来精确地表示时，需要采用*舍入*来解决这个问题。同样，这有几种常见的方法来对数值进行舍入。最简单的方法就是去掉额外不能表示的小数位数，但这样就会使数值变的更负。这种舍入的方法通常被称为向下舍入或向负无穷舍入。当向下舍入到最近的整数时，它就能与floot()函数对应起来，尽管它也可能舍入到其他小数位。 +​ 同样，当一个数值不能通过给定小数位数来精确地表示时，需要采用*舍入*来解决这个问题。同样，这有几种常见的方法来对数值进行舍入。最简单的方法就是去掉额外不能表示的小数位数，但这样就会使数值变的更负。这种舍入的方法通常被称为向下舍入或向负无穷舍入。当向下舍入到最近的整数时，它就能与float()函数对应起来，尽管它也可能舍入到其他小数位。 ![](images/3.5.3_1.png) ​ 同样，也可以使用类似的方法使舍入后数值变得更大(这种称为向上舍入或向正无穷舍入与ceil()函数相对应),向绝对值较小的方向舍入(称为向零舍入与trunc()函数相对应),或向绝对值更大的值舍入(称为远离零舍入或向无穷舍入和round()函数相对应)。然而，这些操作中没有一个总可以最小化舍入误差。一种更好的方法叫做向最接近偶数舍入、收敛舍入或四舍六入五成双，它在lrint()函数中实现。正如你所期望的，这种舍入方式总是选择最近可表示数字。另外，如果有两个误差相等的数，则取偶数。如果最后一个数值是零，那么任意精度数是一个偶数。这种方法是处理IEEE浮点数舍入的默认方法，因为它不仅可以最小化舍入误差，而且还可以确保计算随机数舍入误差之和趋于零。 + ![](images/3.5.3_2.png) + ### 3.5.4 二进制运算 ​ 二进制加法与十进制加法相似，可以简单地对齐每位二进制数并进行加法运算，注意正确处理从一列到下一列的进位操作。注意，两个N位数值加或减的结果通常需要N+1位的数值来正确表示才不能溢出。对于小数，总是在符号位进行增加位宽 @@ -662,7 +661,7 @@ uint32_t p = ((uint32_t) a*(uint32_t) b) >> 12; {% endhint %} -​ 由于这些原因，最好使用c++和Vivado HLS模板类ap_int<>, ap_uint<>, ap_ fixed<>, ap_ufixed<>来表示任意精度数据。ap_int<>和ap_uint<>模板类需要整数参数来定义他们的位宽。计算函数通常产生结果的位宽足够宽，可以表示正确结果，参照3.5.4节中的规则。只有当结果分配给较小的位宽时，才会发生溢出或下溢。 +​ 由于这些原因，最好使用c++和Vivado HLS模板类ap_int<>, ap_uint<>, ap_fixed<>, ap_ufixed<>来表示任意精度数据。ap_int<>和ap_uint<>模板类需要整数参数来定义他们的位宽。计算函数通常产生结果的位宽足够宽，可以表示正确结果，参照3.5.4节中的规则。只有当结果分配给较小的位宽时，才会发生溢出或下溢。 ```C++ #include "ap_int.h" @@ -707,10 +706,10 @@ ap_ufixed<18,12> p = a*b; ​ 在本节中，我们将对优化CORDIC函数的方法提出一些简要的想法和建议。在权衡吞吐量、精度和资源的同时，我们关注不同优化方式如何对处理结果精度产生影响。 在实现CORDIC功能时，选择合适的数据类型来表示角度和结果十分重要。尽管源代码在浮点数和定点数下都能运行，CORDIC通常都会使用定点数数据类型使得乘法器数量减少。而当高效的乘法实现可行时，其它计算三角函数的方法就会被采用。回看图3.3中的源代码，包含了几个跟sigma和factor变量有关的乘法器。通过限制代码只工作在定点数下，我们可以用移位和加法操作来替代。代码如图3.6所示。 - + ```C++ // The file cordic.h holds definitions for the data types and constant values -#include ”cordic.h” +#include "cordic.h" // The cordic phase array holds the angle for the current rotation // cordic_phase[0] =˜ 0.785 @@ -752,6 +751,7 @@ void cordic(THETA_TYPE theta, COS_SIN_TYPE &s, COS_SIN_TYPE &c) c = current_cos; } ``` + Figure 3.6 定点数CORDIC代码实现给定角度sin和cos值 ​ 最后，CORDIC运算得到一个近似值。随着迭代次数的增加，近似值的误差通常会减小。这对应于我们在cordic函数中执行**for**循环的次数，该循环次数由NUM_ITERATION进行设置。即使我们执行了大量迭代，我们仍然可能得到一个近似值。这样做的一个原因是我们可以接近，但从来没有得到完全匹配的目标角度。然而，我们可以通过选择执行更大或更小的迭代角度来调整精度。上述算法中需要更改的地方都可以通过修改NUM_ITERATIONS的值来实现。NUM_ITERATIONS数值依赖于使用这个CORDIC算法应用程序所需精度。 diff --git a/04-Discrete-Fourier-Transform.md b/04-Discrete-Fourier-Transform.md index 378bc433..f40e9632 100755 --- a/04-Discrete-Fourier-Transform.md +++ b/04-Discrete-Fourier-Transform.md @@ -244,12 +244,12 @@ data_loop: 如果说数据的吞吐量我们需要考虑的头号问题，则所有数据都将存储在FF中。这将允许任何元素在每个时钟周期内被访问尽可能多的次数。但是，随着矩阵阵列变大，这种方案也将变得不可行。在矩阵向量乘法的情况下，存储1024位乘以1024位矩阵的32位整数将需要大约4兆字节的存储器。即使使用BRAM来存储，由于每个BRAM存储大约4KBytes，也需要大约1024个BRAM块。另一方面，使用单个大型基于BRAM的内存意味着我们一次只能访问两个元素。这明显降低了性能，如图4.7所示，它需要在每个时钟周期访问多个数组元素（V_In[]的所有8个元素以及M[][]的8个元素）。在实际工程中，大多数设计需要更大的阵列分布存放在更小的BRAM存储器中，这种方法称为阵列分区。较小的数组（通常用于索引较大的数组）可以完全划分为单独的标准变量并映射到FF。匹配流水线选择和数组分区以最大限度地提高运算符使用率和内存使用率是HLS设计探索的一个重要方面。 {% hint style='tip' %} -Vivado HLS 将自动执行一些阵列分区，但由于阵列分区倾向于某些特定设计，因此通常需要我们利用好工具以获得最佳结果。阵列分区的全局配置可在config_array_partiton选项中找到。单个数组可以使用array_patition指令来显示分区，并将指令数组分区完成将数组的每个元素分解到它自己的寄存器中，最终形成基于FF内存。与许多其他基于指令的优化一样，通过手动重写代码也可以实现相同的效果。一般情况下，最好使用工具指令，因为它避免了引入错误并易于代码维护。 +Vivado HLS 将自动执行一些阵列分区，但由于阵列分区倾向于某些特定设计，因此通常需要我们利用好工具以获得最佳结果。阵列分区的全局配置可在config_array_partiton选项中找到。单个数组可以使用array_partition指令来显示分区，并将指令数组分区完成将数组的每个元素分解到它自己的寄存器中，最终形成基于FF内存。与许多其他基于指令的优化一样，通过手动重写代码也可以实现相同的效果。一般情况下，最好使用工具指令，因为它避免了引入错误并易于代码维护。 {% endhint %} 回到图4.4中的矩阵向量乘法代码，我们可以通过添加几个指令来实现高度并行，如图4.11所示，最终的体系结构如图4.12所示。请注意，内部j循环由Vivado HLS 自动展开，因此j在每次使用的时候都被替换为常量。此设计演示了阵列分区的最常见用法，其中分区的数组维度（在本例中为V_In[]和第二维的M[][]）都被索引为常量（在本例中为循环索引j来展开循环）。这使得多路复用器可以无需访问分区阵列架构。 -我们还可以用更少的乘法器降低性能以实现其他设计。例如，在图4.10中，这些设计只使用三个乘法器，因此我们只需要在每个时钟周期读取三个矩阵M[][]和矢量V_in[]的元素。完全分割这些数组会导致额外的多路复用，如图4.13所示。实际上，阵列只需要分成三个物理存储器。同样，这种分区可以通过重写代码手动实现，也可以使用array_patition循环指令在Vivado HLS中实现。 +我们还可以用更少的乘法器降低性能以实现其他设计。例如，在图4.10中，这些设计只使用三个乘法器，因此我们只需要在每个时钟周期读取三个矩阵M[][]和矢量V_in[]的元素。完全分割这些数组会导致额外的多路复用，如图4.13所示。实际上，阵列只需要分成三个物理存储器。同样，这种分区可以通过重写代码手动实现，也可以使用array_partition循环指令在Vivado HLS中实现。 {% hint style='tip' %} 让我们从包含数据的矩阵X来分析： @@ -258,7 +258,7 @@ $$ 1 &2 &3 &4 &5 &6 &7 &8 &9 \end{bmatrix} $$ -数组中指令　array_patition　variable=x factor=2 cycle 将矩阵分为如下两个 +数组中指令　array_partition　variable=x factor=2 cycle 将矩阵分为如下两个 $$ \begin{bmatrix} 1 &3 &5 &7 &9 @@ -296,12 +296,11 @@ data_loop: ![图4.13 II = 3时具有特定的阵列分区选择的矩阵向量乘法体系结构。如左图所示，阵列已被分割得超过需要以至于使用了多路复用器。如右图所示，数组用factor=3分区，这样减少了多路复用器的使用，但j循环索引成为地址计算的一部分。](images/matrix_vector_poor_scaling.jpg) {% hint style='tip' %} -相似地，如果我们使用 **array_patition** variable=x factor=2 block 指令可以将它分为两个矩阵向量： +相似地，如果我们使用 **array_partition** variable=x factor=2 block 指令可以将它分为两个矩阵向量： [1 2 3 4 5] and[ 6 7 8 9] {% endhint %} - -{% hint style=‘info’ %} +{% hint style='info' %} 让我们来研究一下变化的流水线II和阵列分区对性能和面积的影响，比较根据每秒的矩阵向量乘法运算（吞吐量）和根据展开阵列分区因子的数量的性能高低，并绘制相同的区域趋势图（如显示LUT，FF，DSP模块，BRAM）。再思考这两种情况的总趋势是什么？你会选择哪种设计？为什么？ {% endhint %} @@ -330,17 +329,17 @@ data_loop: ``` 图4.14 内部循环展开2倍的矩阵向量乘法代码。 -{% hint style=‘tip’ %} +{% hint style='tip' %} HLS工具可以使用 unroll 指令自动展开循环。该指令采用一个因子作为参数，它是一个正整数，表示循环体应该展开的次数。 {% endhint %} -{% hint style=‘info’ %} +{% hint style='info' %} 使用 **array_partition** cyclic factor = 2 指令和将 M[][]和向量 V_in[]手动划分为单独的数组有着相同的效果。考虑一下应该如何修改代码才能更改访问模式呢？ 现在我们手动将循环展开为两倍。原始代码（没有数组分区和不展开），只执行数组分区的代码，同时执行数组分区和循环展开代码之间的性能结果有哪些不同呢？ 最后，使用指令执行数组分割和循环展开的结果与手动执行的结果相比有哪些不同呢？ {% endhint %} 在这段代码中，我们看到数组分区通常与流水线操作并行执行。通过2倍的数组分割可以使性能提高2倍，我们可以用将内环部分展开2倍或将外环的II减少2倍来实现。提升性能需要相应数量的阵列分区。在矩阵向量乘法中，这种关系相对简单，因为对内部循环中的每个变量只有一次访问权限。在其他代码中，关系可能更复杂。无论如何，设计者的目标应该是确保例化的FPGA资源得到有效利用。一般情况下，将性能提高2倍将使用大约两倍的资源，相反将性能降低2倍可以节约一半的资源。 -{% hint style=‘info’ %} +{% hint style='info' %} 让我们来研究一下循环展开和阵列分区对性能和面积的影响，比较根据每秒的矩阵向量乘法运算（吞吐量）和根据展开阵列分区因子的数量的性能高低，并绘制相同的区域趋势图（如显示LUT，FF，DSP模块，BRAM）。再思考这两种情况的总趋势是什么？你会选择哪种设计？为什么？ {% endhint %} @@ -400,11 +399,11 @@ void dft(IN_TYPE sample_real[N], IN_TYPE sample_imag[N]) { ``` 图4.15：DFT的baseline code -{% hint style=‘info’ %} +{% hint style='info' %} 如果你要使用你设计的CORDIC（如从第3章开始），那么此代码需要做什么修改？ 改变CORDIC核心的准确性会使DFT硬件资源使用情况发生变化吗？ 它会如何影响性能？ {% endhint %} -{% hint style=‘info’ %} +{% hint style='info' %} 请使用HLS实现DFT的基线代码，实现后查看报告，与乘法和加法相比，实现三角函数的相对成本是多少？ 对哪些操作尝试优化更有意义？通过流水线操作内循环可以实现什么性能？ {% endhint %} @@ -428,11 +427,11 @@ void dft(IN_TYPE sample_real[N], IN_TYPE sample_imag[N]) { ![图4.17：图4.16中的行为的流水线版本。在这个设计案例下，由于每个浮点的添加需要4个时钟周期才能完成，并且在下一个循环开始之前需要上一个循环的结果加入计算（以红色显示相关性），所以设置循环的启动间隔为4个间隔。所有迭代的依赖关系汇总在右图中。.](images/dft_pipelined_behavior.jpg) -{% hint style=‘info’ %} +{% hint style='info' %} 重新排列图4.15中代码的循环，并显示您可以使用1的II来管理内部循环。 {% endhint %} -根据DFT中S矩阵的结构，我们可以应用其他优化方式来完全消除三角运算。回想一下，S矩阵的每个元素的复矢量是通过单位圆的固定旋转角度来计算的。S矩阵的行S[0][]对应于单位圆周的零旋转，行S[1][]对应于单次旋转，并且随后的行对应于单位圆周更大角度的旋转。我们可以发现，第二行S[1][]相对应的向量覆盖了来自其他行的所有向量，因为8个列向量每个绕单位圆旋转45°，一共则围绕单位圆旋转了一圈。我们可以通过研究图9.11来直观地确认这个现象。这样我们可以只存储第二行这一次旋转中的正弦和余弦值，然后索引到这个存储器中的值以计算相应行的必要值。这只需要2×N=O（N）个存储单元可以有效减少存储器的O（N）的大小。对于1024个点的DFT，存储器的存储需求将减少到1024×2个条目。假设有一个32位的固定值或浮点值，则只需要8KB大小的片上存储器。与明确存储整个S矩阵相比，明显减少了存储容量。我们将矩阵S的这个一维存储表示为$$s^\prime$$ +根据DFT中S矩阵的结构，我们可以应用其他优化方式来完全消除三角运算。回想一下，S矩阵的每个元素的复矢量是通过单位圆的固定旋转角度来计算的。S矩阵的行$$S[0][]$$对应于单位圆周的零旋转，行$$S[1][]$$对应于单次旋转，并且随后的行对应于单位圆周更大角度的旋转。我们可以发现，第二行$$S[1][]$$相对应的向量覆盖了来自其他行的所有向量，因为8个列向量每个绕单位圆旋转45°，一共则围绕单位圆旋转了一圈。我们可以通过研究图9.11来直观地确认这个现象。这样我们可以只存储第二行这一次旋转中的正弦和余弦值，然后索引到这个存储器中的值以计算相应行的必要值。这只需要2×N=O（N）个存储单元可以有效减少存储器的O（N）的大小。对于1024个点的DFT，存储器的存储需求将减少到1024×2个条目。假设有一个32位的固定值或浮点值，则只需要8KB大小的片上存储器。与明确存储整个S矩阵相比，明显减少了存储容量。我们将矩阵S的这个一维存储表示为$$s^\prime$$ $$ S^\prime=S[1][.]=\begin{bmatrix} @@ -440,19 +439,19 @@ S^\prime=S[1][.]=\begin{bmatrix} \end{bmatrix} $$ -{% hint style=‘info’ %} +{% hint style='info' %} 导出一维数组$$s^\prime$$对应于数组S的行号i和列元素j的输入的访问模式的公式。也就是说，我们如何才能索引到一维S数组以访问元素S（i; j）呢？ {% endhint %} 为了进一步提高性能，我们可以应用一种与矩阵向量乘法非常相似的技术。之前我们发现了提高矩阵向量乘法的性能需要对M[][]数组进行分区。但是如果用$$s^\prime$$ 表示S矩阵则意味着不再有一种有效的方法来划分$$s^\prime$$ 以增加我们在每个时钟上读取的数据量。S的每一个奇数行和列都包括$$s^\prime$$ 的每个元素。因此，我们无法像对S一样对$$s^\prime$$ 的值进行分区。这样增加我们从存储$$s^\prime$$ 的内存中读取数据端口的数量的唯一方法是复制存储。幸运的是，不像与必须读取和写入的内存，复制只读的数组的存储是相对容易的。事实上，Vivado HLS将只对在初始化且从未例化的只读存储器（ROM）自动执行此优化。这种功能的一个优点是我们可以简单地将sin（）和cos（）调用移动到数组初始化中。在大多数情况下，如果此代码位于函数的开头并仅初始化阵列，则Vivado HLS能够完全优化三角函数计算并自动计算ROM的内容。 -{% hint style=‘info’ %} +{% hint style='info' %} 设计一个利用$$s^\prime$$——S矩阵的一维版本的体系结构。这种体系结构是如何影响到所需的存储空间的？与使用二维的S矩阵的实现相比，这会改变逻辑利用率吗？ {% endhint %} 为了有效地优化设计，我们必须考虑代码的每个部分。往往是最薄弱的环节决定了设计的整体性能，这意味着，如果设计有一个瓶颈，将显著影响设计的性能。当前版本的DFT可以对输入和输出数据进行就地操作，即它存储结果相同的数组作为输入数据，输入数组 sample_real 和 sample_imag 都是有效的存储器端口，也就是说，你可以把这些参数的数组存储在相同的存储位置。这样，在任意给定的周期，我们只能获取其中一个阵列的一个数据，这可能会在函数中并行的乘法和加法运算方面产生瓶颈。这就是我们为什么必须将所有的输出结果存储在一个临时数组的原因，然后将所有这些结果复制到函数结尾处的“sample”数组中。如果我们没有执行就地操作，则不需要这样做。 -{% hint style=‘info’ %} +{% hint style='info' %} 修改DFT函数接口，使输入和输出存储在单独的数组中。这会如何影响你的可以执行优化？它如何改变性能？ 区域结果如何？ {% endhint %} diff --git a/06-Sparse-Matrix-Vector-Multiplication.md b/06-Sparse-Matrix-Vector-Multiplication.md index 66ede97a..c12f3423 100755 --- a/06-Sparse-Matrix-Vector-Multiplication.md +++ b/06-Sparse-Matrix-Vector-Multiplication.md @@ -4,13 +4,10 @@ 本章会引入几个与HLS相关的新概念，并进一步深入之前讨论过的优化。本章的目标之一是引入一种更复杂的数据结构。我们用压缩行存储（CRS）来保存稀疏矩阵。另一个目标是演示如何进行性能测试。我们编写了简单的激励用来检验设计是否正确。这在硬件设计中十分重要，**Vivado®HLS** 工具采用HLS C编写激励，并能轻松的对工具生成的RTL代码进行多方面的验证。这是基于HLS设计比基于RTL设计的巨大优势之一。章节中也会讲解如何采用Vivado®HLS工具进行C/RTL联合仿真。不同**SpMV**设计会带来性能上差异，因为执行时间和稀疏矩阵是密切相关的，所以我们必须通过输入数据来确定任务执行之间的间隔以及任务延迟。 - -## 6.1 背景 ## +## 6.1 背景 图6.1显示了一个4x4的矩阵M表示的2种方式。其中图6.1-a采用通用的二维方式16个元素来表示矩阵，每个元素存储在自己对应的位置上。图**6.1**-b采用**CRS**的方式表示相同的矩阵。**CRS** 作为一种数据结构，由3个数组组成。值(**values**)数组保存矩阵中非零元素的值。列索引(**columnIndex**)数组和行指针（**rowPtr**）数组对非零元素的位置信息进行编码。列索引存储每个元素的列数，行指针包含每一行第一个元素在**values**中的索引。**CRS** 结构避免存储矩阵中的0值，确实在数值数组中确实没有存储0。但是在这个例子中，虽然数值数组不保存0，但是列索引数组和行指针数组作为标记信息，表示了矩阵的形态。**CRS** 广泛用于大型的矩阵但是仅仅有少量的非零元素（少于10%或者更低），这样可以简化这类矩阵的存储以及相关的运算。 - - ![图 6.1: M是一个4x4矩阵，用两种方式表示：同"密集"矩阵一样存在二维数组之中；作为稀疏矩阵，以行压缩存储的形式保存，行压缩存储是一种由3个数组组成的数据结构。](images/crs.jpg) 但是，CRS对矩阵的稀疏性没有要求，可以适用于任何矩阵。作为一种针对矩阵的通用方法，但不见得是最高效的。CRS结构也不见得是表示稀疏矩阵最高效的方式，其他稀疏矩阵表示方法也在被使用。 @@ -34,20 +31,16 @@ L1: for (int i = 0; i < NUM_ROWS; i++) { y[i] = y0; } } - ``` ![图6.2: 主体代码演示了系数矩阵向量乘（SpMV）y=M*x的计算。采用CRS的方式，通过rowPtr、columnIndex 和 value 保存矩阵M。第一个for循环通过迭代访问每一行，第二个for循环访问每一列，实现矩阵M中非0元素和向量中对应的元素相乘并保存值在向量y中。](images/placeholder.png) - {% hint style='info' %} 给定一个二维数组表示一个矩阵，通过C代码实现矩阵CRS格式。编写对应的C代码实现将矩阵从CRS格式转化为二维数组的形式。 {% endhint %} 结果表明，通过采用**CRS**的方式，我们能高效的实现稀疏矩阵乘法，不需要将矩阵转化为二维形式。实际上， 对于大型的矩阵仅仅只有一小部分非0元素，稀疏矩阵向量乘法会比第四章中讨论的密集矩阵向量乘高效很多。因为我们直接找到非0元素，并执行非0元素对应的运算。 - - ## 6.2 基本实现 图6.2 提供了基本代码对系数矩阵乘法的实现。函数**spmv**函数有5个参数，分别是**rowPtr**、**columnIndex** ，以及 **values** 对应矩阵 **M** 的 **CRS** 格式中包含的3个参数，这和图6.1中描述的数据结构等价。参数 $$y$$ 用于保存输出的结果，参数x表示输入的被乘向量$$x$$。变量**NUM_ROWS**表示矩阵**M**中行号。变量**NNZ**表示矩阵中非0元素的个数。最后，变量SIZE表示数组x和数组y中元素的个数。 @@ -67,20 +60,17 @@ void spmv(int rowPtr[NUM_ROWS+1], int columnIndex[NNZ], #endif // __MATRIXMUL_H__ not defined ``` -![图6.3： spmv函数和激励的头文件](images/placeholder.png) - +![图6.3： spmv函数和激励的头文件](images/placeholder.png) ## 6.3 测试平台 图6.4 展示了一个针对**spmv**函数测试平台。测试平台通过定义**matrixvector**函数，直接实现矩阵向量乘法，它不考虑矩阵是否为稀疏矩阵以及矩阵是否采用**CRS**方式表示。我们比较**matrixvector**函数输出和**spmv**函数的输出。 - {% hint style='info' %} 在通常的测试平台中，需要实现的函数都会有个“黄金"参考(gold reference)，作为用户期望综合的结果。测试平台会比较黄金参考的输出和通过**Vivado®HLS**综合的代码执行结果。最好的实践方式是，测试平台既可以用于黄金参考，也可用于被综合的代码。这样就保证了两者实现的正确性。 {% endhint %} - 测试平台在主函数**main**中执行。这里我们通过设置**fail**变量初始化为0，当**spmv**函数的输出成结果与**matrixvector**函数输出结果不相同是时，变量置1。定义与矩阵**M**相关的变量、被乘向量$$x$$ 和输出结果$$y$$。对于矩阵**M**，即有普通模式，也有**CSR**模式（保存为**values**、**columnIndex**、**rowPtr**）。矩阵**M** 的**value**如图6.1中所示，输出向量$$y$$有两种，其中y_sw数组保存**matrixvector**函数输出的结果，y数组保存**spmv**函数输出的结果。 在定义好所有的输入变量和输出变量之后，分别调用**spmv**函数和**matrixvector**函数并输入合适的数据。 接下来的**for**循环用于比较y_sw和y中的每一个对应的结果。如果其中一个不相同，则将**fail** 标志置1。最后，程序会打印测试的结果并返回**fail**变量。 @@ -124,18 +114,15 @@ int main(){ return fail; } ``` + ![图6.4 ： 一个简单spmv函数的简单测试平台。测试平台生成了一个用例，并且计算矩阵的向量乘法通过稀疏矩阵乘法（spmv）和非系数矩阵乘法(matrixvector)。](images/placeholder.png) 这个测试平台相对简单并且可能无法充分验证所有的输入都能正常输出。最主要的原因是，它仅仅只用了一个矩阵作为例子，相反，一个好的激励会测试许多矩阵。通常，会通过随机的方式产生输入的测试用例，并且重点测试边界用例。在这个例子中，我们不仅要保证值正确计算，同时保证通过加速器正确的被执行了，而且编译相关的参数也会改变会，不同策略会在实现不同加速器。最重要的是在同一段程序中，我们能通过随机产生很多输入数据来进行测试。编译相关的参数每次发生变化，都需要我们重新编译代码。 - {% hint style='tip' %} 创建一个复杂testbench，通过随机数方式生成许多组测试数据。稀疏矩阵编译相关参数应该是可以修改的（例如，**SIZE**，**NNZ** 等）。创建一个HLS综合脚本，在编译时间参数合理范围改变时，能执行代码很多次。 {% endhint %} - - - ## 6.4 指定循环的属性 如果直接将上述代码进行综合，我们可以得到函数运行的时钟周期及资源占用率。但是，我们不能得到模块执行所需的时钟周期、任务执行的延迟和任务执行之间的间隔。因为这些都取依赖于输入数据，由**spmv**函数外部因素决定。最主要的因素是，内层循环执行的次数是由矩阵**M**中非0元素个数决定的。非0元素的个数在代码中是由常量**NNZ**决定的，虽然可以调用函数计算不同大小的矩阵，但是实际迭代次数是和输入数据相关的。另外，性能也会因为非0元素的分布、综合优化的约束产生不同。更复杂的是，迭代的次数由输入决定，许多可能的输入并没有被遍历。所以，对于工具而言，不通过复杂的分析和额外的信息，工具是不能知道spmv函数执行需要多少时钟周期。**Vivado®HLS** 工具也不能进行上述的分析。 @@ -162,8 +149,6 @@ spmv函数能正常工作的前提条件是什么？证明给定的前提条件 **loop_tripcount** 引导能帮助设计者对函数的性能有个原始的估计。这样,相同的函数通过使用不同的directives优化或者重构代码后的各种方案可以进行比较。但是，这不能确定**min**、**max**和 **avg** 参数。这也很难确定边界条件min和max的值。如果有测试平台，就有一种更准确的方式用于计算**spmv**函数执行的时钟周期数，那就是**C/RTL** 协同仿真。 - - ## 6.5 C/RTL 协同仿真 **C/RTL** 协同仿真能自动化测试Vivado®HLS工具生成的RTL代码，通过已执行综合的代码和testbench实现。每次执行综合以后的代码和提供的测试平台，记录输入和输出结果。输入的值按照时钟转换成**输入向量**。这里的输入向量用于生成RTL设计的RTL级仿真，同时记录**输出向量**。更新综合后的代码， 再次运行测试平台并保存输入和输出数据。测试平台如果返回值是0，则表示成功；若激励返回非0值，则表示失败。 @@ -226,7 +211,6 @@ C/RTL协同仿真能提供可变循环边界函数的延迟。它反馈函数运 | case10 | - | pipeline,unroll=8,cyclic=8 | | case11 | - | pipeline,unroll=8,block=8 | - 如果你完成了上述练习，你会发现盲目的使用优化directives，可能不会得到你期望的结果。通常在设计时， 在思考下考虑应用的特性，选择针对设计的特定优化方式。当然，这也需要一些直觉能力和一些专用工具投入使用。虽然，搞清楚像**Vivado®HLS**这样复杂工具中每一个细节是困难乃至不可能的，但是我们能基于关键的方面建立思考模型。 上面我们在用例3和4中考虑对外层循环**L1**进行流水化操作而不是对内层循环。这种变化针对一个任务，可以提高潜在的并行程度。为了完成优化，**Vivado®HLS** 工具必须展开代码中所有的内层循环**L2** 。如果循环能全部展开，这样能减少计算循环边界的时间，同时也能消除递归（recurrences）。但是代码中的内层循环Vivado HLS是无法完全展开的，因为循环边界不是常量。 @@ -242,6 +226,7 @@ C/RTL协同仿真能提供可变循环边界函数的延迟。它反馈函数运 注意图**6.6**中的代码和自动循环展开的代码是由一点点区别的。自动循环展开复制计算，但是保留每次计算先后顺序（除了当前的例子）。这就导致了计算顺序由内层循环决定，如图6.7左所示。对计算顺序进行调整后，操作上的依赖关系如图6.7 左边所示。在当前的代码中，最后累加求和是一个递归（recurrence ）。当使用浮点数据类型时，这种调整计算顺序的操作可能对程序产生改变，所以Vivado HLS对这种类型的代码不进行操作顺序自动调整。 这个设计可能会被综合、实现如图**6.8**所示的结果。在这个例子中，$$S=3$$与$$II$$最匹配，乘法器的延迟正好是3。所有的运算过程都是在一个乘法器和加法器上执行。比较这个例子与图**6.5**中的例子，我们可以发现一些缺点。最明显的是，内层循环的流水线长度很长，实现的时候需要多个更多的周期才能执行下一次外层**L1**循环。处理一行中非零元素和执行块**S** 相同。一行有个3个元素和一行有一个元素计算的时间是相同的。剩下的运算也需要在循环流水线中执行，即使他们的结果没有用。为了严格的比较两个设计的特性，我们需要了解设计对矩阵每行非零元素个数的预期。 + ```c #include "spmv.h" @@ -266,14 +251,13 @@ void spmv(int rowPtr[NUM_ROWS+1], int columnIndex[NNZ], } } ``` + ![图6.6 局部展开图6.2中**smpv**函数](images/placeholder.png) ![图6.7 针对累加的两种不同方式的局部展开。左边的版本有3个加法器进行递归操作，相反右边的版本只有1个加法器进行递归累加](images/spmv_partial_unroll.jpg) - ![图6.8 图6.6中 spmv函数基于部分展开和内部流水线处理后执行过程](images/spmv_unrolled_behavior.jpg) - 如果矩阵每行非零元素很少，则采用第一种实现方式较优；如果矩阵中每行非零元素较多，则第二种实现方式更好。 需要注意，这里存在一个关于先有鸡还是先有蛋的问题。我们需要知道目标器件和时钟周期，这样才能确定流水线中加法器能不能满足时序要求。只有在我们知道流水线的级数之后（也许S=1时，Vivado HLS才能识别到加法递归），我们才能选择合适版本的参数S，来满足$$II=1$$。一旦我们确定了**S**，我们能通过**C/RTL**协同仿真来，通过一组测试数据，确定是不是达到了性能上的要求。因为循环边界是可变的，所以得到的性能参数是依赖于数据的，所以我们需要设定不同的**S**，来找到性能的最大值。改变器件的类型和工作频率会影响之前所有的设计！尽管看来去高层次综合（**HLS**）对解决问题提供的帮助不多，相比于RTL开发新版本然后进行验证，它开发起来快（代码编写方便）。 diff --git a/07-Matrix-Multiplication.md b/07-Matrix-Multiplication.md index 2a743c8e..39dc9a44 100755 --- a/07-Matrix-Multiplication.md +++ b/07-Matrix-Multiplication.md @@ -79,7 +79,6 @@ $$ \qquad(7.4) $$ - ```c void matrixmul(int A[N][M], int B[M][P], int AB[N][P]) { #pragma HLS ARRAY_RESHAPE variable=A complete dim=2 @@ -97,7 +96,6 @@ void matrixmul(int A[N][M], int B[M][P], int AB[N][P]) { } } } - ``` ![图7.1 一个通用的3层for循环实现了矩阵乘法。](images/placeholder.png) @@ -138,7 +136,7 @@ void matrixmul(int A[N][M], int B[M][P], int AB[N][P]) { 去掉**array_reshape** directive。这对性能有什么影响？这对资源使用率有什么影响？通过改变**array_reshape** 的参数，对这些数组有没有影响？这种情况下，和只使用**array_reshape** 有什么区别？ {% endhint %} -## 7.3 块矩阵乘法 +## 7.3 块矩阵乘法 **块矩阵** 是被分割出来子矩阵。直观的理解是，可以通过在水平方向和垂直方向画线对矩阵元素进行分割。得到的块作为原始矩阵的子矩阵。反过来，我们可以认为原始矩阵是块矩阵的集合。这就自然的派生出许多算法在线性代数运算中，例如矩阵的乘法，通过把大型的矩阵切分成许多小的矩阵，然后基于块进行运算。 @@ -190,8 +188,8 @@ typedef struct { DTYPE out[BLOCK_SIZE][BLOCK_SIZE]; } blockmat; void blockmatmul(hls::stream &Arows, hls::stream &Bcols, blockmat & ABpartial, DTYPE iteration); #endif - ``` + ![图7.4：分块矩阵乘法的头文件。文件定义了函数内部使用的数据结构，最关键的是，函数blockmatmul的接口。](images/placeholder.png) {% hint style='info' %} @@ -253,12 +251,13 @@ void blockmatmul(hls::stream &Arows, hls::stream &Bcols, } } ``` + ![图7.5： 函数blockmatmul从矩阵A输入一系列大小为BLOCK_SIZE的行，从矩阵 B 输入一系列大小为BLOCK_SIZE 的列，创建一个BLOCK_SIZWxBLOCK_SIZE 分布结果，作为矩阵AB 的一部分。代码的第一部分（标记为loadA ）保存 A 的行到局部存储，第二部分是嵌套的partialsum for循环执行部分结果的计算，最后的一部分（标记为 writeoutput ）将之前计算返回的分布结果放到合适格式中。](images/placeholder.png) 需要对**stream**类进行一些解释说明才能完全掌握这段代码，并且使用它。**stream** 类型变量**Arows**由许多**blockvec**类型的元素构成。**blockvec** 是一个**BLOCK_SIZE**x**BLOCK_SIZE**的矩阵。每个**Arows**中的元素都是一个数组，用于保存矩阵**A**中**BLOCK_SIZE**行的数据。所以在每次调用**blockmatmul**函数时，**Arow** 流会有**SIZE**个元素，每个用于保存**BLOCK_SIZE**行。语句**tempA = Arows.read()** 从**Arows** 流中取一个元素。然后，将这些对应的元素赋值给局部矩阵**A** 中对应的位置。 {% hint style='info' %} -**stream** 类对操作符 << 进行了重载，它等价于函数**read()** 。所以，语句** tempA = Arows.read()** 和 ** tempA << Arows** 执行时等价的。 +**stream** 类对操作符 << 进行了重载，它等价于函数**read()** 。所以，语句 **tempA = Arows.read()** 和 **tempA << Arows** 执行时等价的。 {% endhint %} 剩余部分的计算式是部分求和。大部分都是在函数**blockmatmul** 中进行。 @@ -312,7 +311,7 @@ $$ 外部的两层**for**循环每一步都以块的形式实现输入矩阵。你可以发现每次迭代都是以**BLOCK_SIZE**作为步进。接下来的两个**for**循环把矩阵**A**的行写到**strm_matrix1_element**，把矩阵**B**的列写到**strm_martix2_element**中。在执行上述步骤时都是一个元素一个元素进行的，通过使用变量**k**访问单个的值从行（列）然后将它们写入到一维数组中。注意，**strm_matrix1_element** 和**strm_martix2_element**都是**blockvec**数据类型，即长度为**BLOCK_SIZE**的一维数组。它们用于保存行或列的**BLOCK_SIZE**个元素。内层循环迭代**BLOCK_SIZE**次。流类型变量**strm_matrix1**和**strm_matrix2**被写入**SIZE**次。换句话说，就是对整行（列）进行缓存，缓存中每个元素保存**BLOCK_SIZE**个值。 {% hint style='info' %} -类**stream**通过重载**>>** 操作符，它等价于函数 **write（data)** 。这和对 **read()** 函数重载，使其等价于操作符**<<** 。因此，表达式 **strm_matrix1.write(strm_matrix1_element)** 和 **strm_matrix1_element >> strm_matrix1** 执行是相同的。 +类**stream**通过重载**>>** 操作符，它等价于函数 **write(data)** 。这和对 **read()** 函数重载，使其等价于操作符**<<** 。因此，表达式 **strm_matrix1.write(strm_matrix1_element)** 和 **strm_matrix1_element >> strm_matrix1** 执行是相同的。 {% endhint %} ```c @@ -354,6 +353,7 @@ int main() { } //the remainder of this testbench is displayed in the next figure ``` + ![图7.6： 分块矩阵乘法测试平台的第一部分。函数被分为两个部分， 因为太长不能在单页上一次显示。测试平台剩下的部分在图7.7中。其中有一个“软件”版本的矩阵乘法、变量的声明和初始化。](images/placeholder.png) ```c @@ -393,8 +393,8 @@ int main() { return 0; } - ``` + ![图7.7：分块矩阵乘法的第二部分。第一部分在图7.6中。这部分演示如何将流类型数据送入函数blockmatmul，代码测试这个函数结果和更简单的三重for循环得出的结果是否匹配。](images/placeholder.png) 这段代码最后一部分从高亮的**if**语句开始。这和矩阵**A**的值相关。本质上，这里是为了我们不用经常向**strm_matrix1**中写入相同的值。矩阵**A**的值可以对应多次**blockmatmul**函数的调用。图7.3是关于这的讨论。这里的**if**语句高亮是为了强调，不能在这里一次次写入相同值。因为函数**blockmatmul**只在有必要的时候读这些数据。所以当我们连续写数据，代码执行不正常，因为流写入的数据是多于读出的数据。 @@ -419,6 +419,6 @@ int main() { 比较分块矩阵乘法和矩阵乘法之间的性能。当增大矩阵的大小，性能如何变化？分块大小是否在性能中起到了重要作用？选择两种有相同资源占用率的架构，这两种架构之间性能差异有多大？ {% endhint %} -## 7.4 总结 ## +## 7.4 总结 矩阵乘法提供了一种不同的方式来实现矩阵乘法。函数通过流的方式输入矩阵的部分，然后计算矩阵的部分结果。函数通过多次计算，完成对整个矩阵乘法的计算。 diff --git a/08-Prefix-Sum-and-Histogram.md b/08-Prefix-Sum-and-Histogram.md index ce84c438..efbea5f2 100755 --- a/08-Prefix-Sum-and-Histogram.md +++ b/08-Prefix-Sum-and-Histogram.md @@ -28,7 +28,6 @@ $$ ​这段代码的编写方式是将每个输出数值都写入输出寄存器out[]中，然后在下一次迭代中再次从寄存器中读出上一次输出的数值。由于读取寄存器的延迟是1个时钟周期，因此从寄存器中读取的数据只有在下一个时钟周期才能被处理。结果是，这样的代码设计只能实现II (Initiation interval)=2 的循环设计。在这种情况下，有一种简单的改良此代码的方法：我们可以使用一个单独的局部变量来进行累加操作，而不是像以前一样从数组中读回前一次累加的数值。在CPU处理器代码设计中，避免使用额外的外部存储器来替代寄存器作为数据访问的方式更有优势，同样在HLS设计中这样的数据访问方式更重要，因为其他的处理操作很少能成为系统的性能瓶颈。该操作方式代码如图8.2所示。 - ![图8.1 前缀和实现代码及行为描述](images/8_1.png) ![图8.2 优化后的前缀和代码及行为描述](images/8_2.png) @@ -73,7 +72,6 @@ void histogram(int in[INPUT_SIZE], int hist[VALUE_SIZE]) { ## 8.3 直方图优化和错误依赖 - ​让我们更深入的观察上面的处理过程。在循环的第一次迭代中，我们读取位置x0处的hist数组值并将其写回到相同的位置x0处hist数组中。由于读操作需要一个时钟延迟，所以写入数组操作必须在下一个时钟周期发生。然后在下一次迭代循环中，我们读取另一个hist数组中另一个位置x1中数组。x0和x1都取决于输入值，并可以取任何值。因此我们考虑到综合成电路时最坏的情况，如果x0 = x1,则在前一个写入完成前，位置x1处的读取无法开始。因此，我们必须在读写之间进行切换。 ​事实证明，只要x0和x1是独立的，我们就必须在读写之间进行切换。如果它们实际上不是独立的呢？例如，我们可能知道数据源实际不会连续产生两个完全相同的数据。那我们现在该怎么做呢？如果我们可以将x0和x1是不同的地址额外信息提供给HLS工具，那么它就能够同时在位置x1处读取，而在位置x0处写入数据了。在Vivado@HLS中，可以通过设置**dependence**指令来完成。 @@ -82,7 +80,7 @@ void histogram(int in[INPUT_SIZE], int hist[VALUE_SIZE]) { ```c #include -#include ”histogram.h” +#include "histogram.h" // Precondition: hist[] is initialized with zeros. // Precondition: for all x, in[x] != in[x+1] void histogram(int in[INPUT SIZE], int hist[VALUE SIZE]) { @@ -116,7 +114,7 @@ void histogram(int in[INPUT SIZE], int hist[VALUE SIZE]) { {% hint style='tip' %}对于图8.10中的代码，你可能会问为什么像Vivado@HLS这样的工具无法确定性能。事实上，虽然在这样的一些简单情况下，更好的代码分析可以将if条件属性传播到每个分支，但我们必须接受存在一些代码段，其中内存访问的性能是不可判定的。在这种情况下最高的性能只能通过添加用户信息的静态进程来实现。最近的一些研究工作通过引入一些研究来寻求改进设计中的动态控制逻辑[[60](./BIBLIOGRAPHY.md#60), [40](./BIBLIOGRAPHY.md#40),[19](./BIBLIOGRAPHY.md#19)]。 {% endhint %} ```c -#include ”histogram.h” +#include "histogram.h" void histogram(int in[INPUT_SIZE], int hist[VALUE_SIZE]) { int acc = 0; int i, val; @@ -134,7 +132,7 @@ void histogram(int in[INPUT SIZE], int hist[VALUE SIZE]) { old = val; } hist[old] = acc; -} +} ``` ![图8.10：从for循环中删除写入依赖后的读取。 这需要一个if /else结构，看起来似乎会增加设计的不必要的复杂性。尽管数据路径更复杂，但它允许更有效的流水线操作。](images/placeholder.png) @@ -149,7 +147,7 @@ void histogram(int in[INPUT SIZE], int hist[VALUE SIZE]) { ​不过，并不是没有指望了。因为我们可以通过将直方图计算分解为两个阶段来达到更多的并行性。在第一阶段，我们将输入数据分成若干独立的分块。每个分块的直方图可以使用我们之前的直方图解决方案独立计算。第二阶段，将各个直方图组合起来生成完整数据集的直方图。这种分块（或映射）和合并（或还原）过程与MapReduce框架[[20](./BIBLIOGRAPHY.md#20)]采用的过程非常相似，并且是并行计算的常见模式。Map-reduce模式适用于包含交换和关联操作的循环，例如这种情况下的加法。完整方案如图8.12所示。 -​图8.13中式实现这种架构的代码。直方图函数实现了Map-reduce模式的’map’部分，并且会多次实例化。该代码与图8.10中的代码非常相似。主要的区别在于我们添加了额外的代码来初始化hist数组。**Histogram_map** 函数输入数组是hist数组中一个分区数据。**Histogram_reduce** 函数实现了模式中的“还原”部分。它将分区数据的直方图作为输入，并通过将每个直方图的计数相加，将它们组合成完整的直方图。在我们的图8.13的代码示例中，我们只有两个处理对象，因此将两个输入数组hist1和hist2合并。这可以很容易的扩展以处理更多的元素。 +​图8.13中式实现这种架构的代码。直方图函数实现了Map-reduce模式的’map’部分，并且会多次实例化。该代码与图8.10中的代码非常相似。主要的区别在于我们添加了额外的代码来初始化hist数组。**histogram_map** 函数输入数组是hist数组中一个分区数据。**histogram_reduce** 函数实现了模式中的“还原”部分。它将分区数据的直方图作为输入，并通过将每个直方图的计数相加，将它们组合成完整的直方图。在我们的图8.13的代码示例中，我们只有两个处理对象，因此将两个输入数组hist1和hist2合并。这可以很容易的扩展以处理更多的元素。 ![图8.12：使用map-reduce模式实现的直方图计算。 a）部分中的处理单元（PE）与图8.11所示的结构相同。 in数组是分块后的，每个分块由一个单独的PE处理。合并块组合各个直方图以创建最终直方图。](images/architectures_histogram_parallel.jpg) @@ -157,10 +155,10 @@ void histogram(int in[INPUT SIZE], int hist[VALUE SIZE]) { {% hint style='info' %}修改图8.13中的代码，使得其支持可参数化改变的**PE**数目（变量名记为**NUM_PE**）？提示：你需要根据数组分块数量**NUM_PE**以及循环将数组合并成一个数组。当你改变PE的数量时，吞吐量和任务间隔会发生什么变化？{% endhint %} -​我们在**histogram**函数中使用**dataflow**指令来达到任务级流水线设计。在这种情况下有三个模块：两个**partial_histogram**函数实例和一个**histogram_reduce**函数实例。在一个任务中，因为两个**partial_histogram**处理的数据相互独立，所以可以同时执行。**Histogram_reduce** 函数处理过程必须在 **partial_histogram** 处理完成后才开始。因此，**dataflow** 指令本质上是创建了一个两阶段的任务管道。第一阶段执行**partial_histogram**函数，而第二阶段执行**histogram_reduce**函数。与任何数据流设计一样，整个histogram函数的间隔取决于两个阶段的最大启动间隔。第一个阶段的两个**partial_histogram**函数时相同的，并且具有相同的间隔（$$I{I_{histogram\_map}}$$）。**Histogram_reduce** 函数则是另一个间隔（$$I{I_{histogram\_reduce}}$$）。顶层 **histogram** 函数的启动间隔$$I{I_{histogram}}$$是max($$I{I_{histogram\_map}}$$ ,$$I{I_{histogram\_reduce}}$$ ). +​我们在**histogram**函数中使用**dataflow**指令来达到任务级流水线设计。在这种情况下有三个模块：两个**partial_histogram**函数实例和一个**histogram_reduce**函数实例。在一个任务中，因为两个**partial_histogram**处理的数据相互独立，所以可以同时执行。**histogram_reduce** 函数处理过程必须在 **partial_histogram** 处理完成后才开始。因此，**dataflow** 指令本质上是创建了一个两阶段的任务管道。第一阶段执行**partial_histogram**函数，而第二阶段执行**histogram_reduce**函数。与任何数据流设计一样，整个histogram函数的间隔取决于两个阶段的最大启动间隔。第一个阶段的两个**partial_histogram**函数时相同的，并且具有相同的间隔（$$I{I_{histogram\_map}}$$）。**histogram_reduce** 函数则是另一个间隔（$$I{I_{histogram\_reduce}}$$）。顶层 **histogram** 函数的启动间隔$$I{I_{histogram}}$$是max($$I{I_{histogram\_map}}$$ ,$$I{I_{histogram\_reduce}}$$ ). ```c -#include ”histogram parallel.h” +#include "histogram_parallel.h" void histogram_map(int in[INPUT_SIZE/2], int hist[VALUE_SIZE]) { #pragma HLS DEPENDENCE variable=hist intra RAW false for(int i = 0; i < VALUE_SIZE; i++) { diff --git a/09-Video-Systems.md b/09-Video-Systems.md index 490f0d7a..7e2395b7 100755 --- a/09-Video-Systems.md +++ b/09-Video-Systems.md @@ -18,7 +18,9 @@ ​第三个方面是由于眼睛中的眼杆和感光部分对绿光比对红光更加敏感，并且大脑主要从绿光中获取亮度信息。因此，图像传感器和显示器采用马赛克模式，比如采样2个绿色像素和1个红像素与1个蓝像素比例的**Bayer**模式。于是在相同数量焦平面采集单元或显示器显示单元条件下系统可以采集或显示更高分辨率的图像，从而降低图像采集传感器或显示器的制造成本。 -{% hint style='info' %}视频系统多年来一直围绕人类视觉系统进行设计。最早的黑白相机对蓝绿光敏感，以匹配眼睛对该颜色范围内亮度的敏感度。然而，不幸的是它们对红光不太敏感，因此红色（如红妆）在相机上看起来不正常。当时针对这个问题的解决方案绝对是低技术含量的：拍照时演员穿着华丽的绿色和蓝色化的妆。 {% endhint %} +{% hint style='info' %} +视频系统多年来一直围绕人类视觉系统进行设计。最早的黑白相机对蓝绿光敏感，以匹配眼睛对该颜色范围内亮度的敏感度。然而，不幸的是它们对红光不太敏感，因此红色（如红妆）在相机上看起来不正常。当时针对这个问题的解决方案绝对是低技术含量的：拍照时演员穿着华丽的绿色和蓝色化的妆。 +{% endhint %} ### 9.1.2 数字视频格式 @@ -28,12 +30,14 @@ ![图9.2：1080P@60Hz高清视频信号中的典型的同步信号](images/video_syncs.jpg) -{% hint style='info' %}数字视频信号均是对原来电视标准的模拟信号（美国的NTSC和许多欧洲国家的PAL）通过抽样、量化和编码后形成。由于用于模拟阴极射线管扫描的硬件具有有限的电平转换速率，所以水平和垂直同步间隔的时间要满足从硬件恢复扫描到下一行的开始。同步信号由若干个时钟周期的低电平组成。另外，由于电视无法有效地显示靠近同步信号的像素，因此引入前肩信号和后肩信号，这样就可以显示更多的图片像素。即使如此，由于许多电视机都是基于**电子束扫描**原理设计的，所以其中的图像数据有20％的边框处的像素是不可见。{% endhint %} +{% hint style='info' %} +数字视频信号均是对原来电视标准的模拟信号（美国的NTSC和许多欧洲国家的PAL）通过抽样、量化和编码后形成。由于用于模拟阴极射线管扫描的硬件具有有限的电平转换速率，所以水平和垂直同步间隔的时间要满足从硬件恢复扫描到下一行的开始。同步信号由若干个时钟周期的低电平组成。另外，由于电视无法有效地显示靠近同步信号的像素，因此引入前肩信号和后肩信号，这样就可以显示更多的图片像素。即使如此，由于许多电视机都是基于**电子束扫描**原理设计的，所以其中的图像数据有20％的边框处的像素是不可见。 +{% endhint %} ​如图9.2所示，典型1080P（60Hz）视频帧共包含2200 * 1125个数据采样像素点。在每秒60帧的情况下，这相当于每秒总共148.5百万的采样像素点。这比1080P（60Hz）帧中的有效视频像素的平均速率（1920 * 1080 * 60 =每秒124.4百万像素）高很多。现在的大多数FPGA都可以满足以这个时钟速率进行高效的处理，并且是在一个时钟周期进行一次采样的方式下运行的。对于更高分辨率的处理系统需求，如4K或2K的数字影院，或者每秒120帧甚至240帧的处理需求，就需要一个时钟周期采样更多像素点，通过增加数据处理的吞吐量方式来运行。注意，我们通常可以通过在HLS中展开循环的方式来生成此类结构（请参见1.4.2小节）。同样，当处理低分辨率和低帧率需求时，优选的方案是采用操作共享的方式，多个时钟周期去处理一次采样。这样的处理结构是通过指定循环的处理间隔方式实现。 ```c -#include ”video_common.h” +#include "video_common.h" unsigned char rescale(unsigned char val, unsigned char offset, unsigned char scale) { return ((val − offset) ∗ scale) >> 4; } @@ -78,7 +82,9 @@ rgb_pixel pixel_out[MAX_HEIGHT][MAX_WIDTH]) { } ``` -{% hint style='info' %}高速计算机视觉应用是需要满足每秒处理10000帧的200*180像素的视频帧的速度要求。这样的应用通常是使用一个高速图像传感器直接和FPGA直接连接，并且中间不需要有同步信号。在这种情况下，你会每个时钟周期内进行多少次采样处理呢？FPGA可以很好的完成处理么？答案是你可以使用HLS编写嵌套循环的结构来完成这样的设计。{% endhint %} +{% hint style='info' %} +高速计算机视觉应用是需要满足每秒处理10000帧的200*180像素的视频帧的速度要求。这样的应用通常是使用一个高速图像传感器直接和FPGA直接连接，并且中间不需要有同步信号。在这种情况下，你会每个时钟周期内进行多少次采样处理呢？FPGA可以很好的完成处理么？答案是你可以使用HLS编写嵌套循环的结构来完成这样的设计。 +{% endhint %} ### 9.1.3 视频处理系统架构 @@ -86,7 +92,9 @@ rgb_pixel pixel_out[MAX_HEIGHT][MAX_WIDTH]) { ​默认情况下，Vivado@HLS软件会将代码中的函数数组接口综合成简单的存储器接口。其中，存储器写数据接口由地址总线、数据总线和写使能信号线组成。存储器每次读取和写入数据均有相应地址变化，并且数据读取和写入时间固定。如图9.4所示，将这种接口与片上存储资源Block RAM资源集成在一起使用很简单方便。但是即使在大容量的芯片，片上Block RAM资源也是很稀缺的，因此如果存储大容量的视频资源会很快消耗完片上Block RAM资源。 -{% hint style='info' %}针对每像素24位的1920x1080视频帧，芯片需要消耗多少BlockRAM资源才能完成存储存一帧视频帧？片上Block RAM资源最多又能存储多少帧呢？{% endhint %} +{% hint style='info' %} +针对每像素24位的1920x1080视频帧，芯片需要消耗多少BlockRAM资源才能完成存储存一帧视频帧？片上Block RAM资源最多又能存储多少帧呢？ +{% endhint %} ![图9.5 将视频处理设计与外部存储接口集成](images/video_DDR_interface.jpg) @@ -109,6 +117,7 @@ void video_filter(pixel_t pixel_in[MAX_HEIGHT][MAX_WIDTH], pixel_t pixel_out[MAX_HEIGHT][MAX_WIDTH]) { #pragma HLS interface s_axi port = pixel_out #pragma HLS interface s_axi port = pixel_in +} ``` ```c @@ -144,7 +153,9 @@ void video_filter(pixel_t pixel_in[MAX_HEIGHT][MAX_WIDTH], ​在视频处理算法中，通常计算一个输出像素的值时需要输入像素及其周围的像素的值做参考，我们把储存计算所需要的输入像素值的区域叫**窗口**。从概念上讲，按照窗口大小去Z型扫描图片，然后根据扫描到的像素值计算出输出结果就可以。如图9.9所示示例代码中，示例代码是针对视频帧进行二维滤波。示例代码中在计算每一个输出像素结果前需要从输入视频帧中读取相应的窗口像素数据（存储在数组中）。 -{% hint style='info' %} 在图9.9中，代码中包含的`int wi =row + i – 1;int wj = col +j -1;`, 解释这些表达式为什么包含”-1”这一项。提示：如果滤波器核换成7×7，而不是3×3，”-1”这个数字项会改变么？ {% endhint %} +{% hint style='info' %} +在图9.9中，代码中包含的`int wi =row + i – 1;int wj = col +j -1;`, 解释这些表达式为什么包含”-1”这一项。提示：如果滤波器核换成7×7，而不是3×3，”-1”这个数字项会改变么？ +{% endhint %} ​注意，在此代码中，每计算一个输出像素值，必须多次读入像素并填充到窗口缓冲区中。如果每个时钟周期只能执行一次读操作，示例代码的执行性能会受限于读入像素速率大小。因此示例代码基本上是图1.8中所示一维滤波器的二维版本。另外，因为输入像素不是以正常扫描线顺序读取的，所以接口约束也只有有限的几个可选。（本主题将在9.1.3小节中更详细讨论）。 @@ -197,20 +208,25 @@ row_loop: for (int row = 0; row < MAX_HEIGHT; row++) { ``` ![图9.9: 没有使用线性缓冲区的2D滤波代码](images/placeholder.png) - ​如图9.10所示代码，是使用线性缓冲区和窗口缓冲区方式实现的，其中代码实现过程如图9.11所示。代码每次执行一次循环时，窗口会移动一次，使用来自于1个输入像素和两个来自于行缓冲区缓存的像素来填充窗口缓冲区。另外，输入的新像素被移入线性缓冲区，准备被下一行的窗口运算过程所使用。请注意，由于为了每个时钟周期处理一个像素并输出结果，所以系统必须在每个时钟周期内完成对窗口缓冲区数据的读取和写入。另外，当展开”i”循环后，对于窗口缓冲区数组而言，每个数组索引都是一个常量。在这种情况下，Vivado@HLS将转换数组中的每个元素成一个标量变量（一个成为**标量化**的过程）。窗口数组中的元素随后使用触发器(FF)编译实现。同样，线性缓冲区的每一行数据都会被访问两次（被读取一次和写入一次）。示例代码中明确约束行缓冲区中每一行数组元素被分割成到一块单独存储区中。根据MAX WIDTH的可能取值情况，最后决定使用一个或者多个Block RAM实现。注意，每个Block RAM可在每个时钟周期支持两次独立访问。 {% hint style='info' %}线性缓冲区是应用于模块式计算的**重用缓冲区**的一种特殊例子。如图9.9中所示，重用缓冲区和线性缓冲区的高级综合是目前一个热点研究领域。参考[[8](./BIBLIOGRAPHY.md#8), [31](./BIBLIOGRAPHY.md#31)]。{% endhint %} -{% hint style='info' %}Vivado@HLS包含`hls::linebuffer<>`和`hls::window buffer<>`类，它们可以简化窗口缓冲区和行缓冲区的管理。{% endhint %} +{% hint style='info' %} +Vivado@HLS包含`hls::linebuffer<>`和`hls::window buffer<>`类，它们可以简化窗口缓冲区和行缓冲区的管理。 +{% endhint %} -{% hint style='info' %}对于3×3的图像滤波器核，存储每个像素占用4字节的1920×1080图像的一行数据需要使用多少个FPGA片上Block RAM?{% endhint %} +{% hint style='info' %} +对于3×3的图像滤波器核，存储每个像素占用4字节的1920×1080图像的一行数据需要使用多少个FPGA片上Block RAM? +{% endhint %} ### 9.2.2 因果滤波器 ​图9.10中实现的滤波器是每个时钟周期读取一个输入像素，计算出一个输出像素。该滤波器实现原理与图9.9中所示不完全相同。输出结果是从先前读取像素的窗口缓冲区中数据计算出来的。窗口缓冲区中取数顺序是“向上和向左”。因此，输出图像相对于输入图像是“向下和向右”移动的。这种情况类似于信号处理中的因果滤波器和非因果滤波器的原理。大多数信号处理理论的重点是因果滤波器。因为只有因果滤波器对于时间采样信号（例如，其中$$x[n] = x(n*T)$$和​$$y[n] = y(n*T)$$）才有实用价值。 -{% hint style='info' %}在**因果**滤波器$$h[n]$$中，$$\forall k < 0,h[k] = 0$$。虽然，有限滤波器$$h[n]$$不是因果滤波器，但是可以通过延时方法使得而被$$\wedge h[n] = h[n - D]$$转换为因果滤波器$$h[n]$$。新滤波器$$\wedge y = x \otimes \wedge h$$的输出与旧滤波器$$y = x \otimes h$$的延时输出相同。具体而言，$$ \wedge y = y[n - D]$$。{% endhint %} +{% hint style='info' %} +在**因果**滤波器$$h[n]$$中，$$\forall k < 0,h[k] = 0$$。虽然，有限滤波器$$h[n]$$不是因果滤波器，但是可以通过延时方法使得而被$$\wedge h[n] = h[n - D]$$转换为因果滤波器$$h[n]$$。新滤波器$$\wedge y = x \otimes \wedge h$$的输出与旧滤波器$$y = x \otimes h$$的延时输出相同。具体而言，$$ \wedge y = y[n - D]$$。 +{% endhint %} ```c void video_2dfilter_linebuffer(rgb_pixel pixel_in[MAX_HEIGHT][MAX_WIDTH], @@ -248,8 +264,8 @@ row_loop: for (int row = 0; row < MAX_HEIGHT; row++) { ![图9.11:图中所示为图9.10中的代码实现的存储器。 存储器在特定的循环周期时将从线性缓冲区中读取出数据存储在窗口缓冲区的最右端区域。黑色的像素存储在行缓冲区中，红色的像素存储在窗口缓冲区中。](images/video_buffers.jpg) - -{% hint style='info' %}通过卷积的定义证明前面原理： +{% hint style='info' %} +通过卷积的定义证明前面原理： $$ y = x \otimes h:y[n] = \sum\limits_{k = - \infty }^\infty {x[k] * h[n - k]} $$ diff --git a/10-Sorting-Algorithms.md b/10-Sorting-Algorithms.md index 588359e1..4cd262fc 100755 --- a/10-Sorting-Algorithms.md +++ b/10-Sorting-Algorithms.md @@ -22,8 +22,6 @@ $$ 最终确认A(5)确实包含在序列中。 - - 在上面的例子中，序列A可以代表各种不同的概念。A可以代表一个数学集合，可以搜索元素是否在集合内。在集合中，A也可以表示数据中的一部分，通常称为关键字，用于索引其余信息。关键字可以是一个人的名字，基于关键字的搜索，可以获取该人的其余数据信息，比如他们的出生日期等等。在有些场合下，关键字还可能是更抽象的对象，比如一些数据或密钥的加密散列。在这种情况下，数据存储顺序可能是随机的，但只要获得正确的密码散列，仍然可以搜索它。不管何种情况，排序和搜索所需的基本操作类似，都要需要比较两个不同的值，在这一章中将忽略这些问题。 在通用处理器系统中，有各种各样的排序算法应用，参考[[36](./BIBLIOGRAPHY.md#36)]。它们在空间复杂度和时间复杂度方面有所不同，但绝大数都需要$$O(NlogN)$$次比较才能对$$N$$个元素进行排序。对一个排序序列使用二叉树插入1个新元素时，搜索$$O(logN)$$次即能找到新值的正确位置；当需要插入$$N$$个元素时，这个过程需要重复$$N$$次。 @@ -69,11 +67,11 @@ void insertion_sort(DTYPE A[SIZE]) { } } ``` + ![图10.2 : 插入排序的参考代码。外部循环for用来遍历数组元素，内部循环while用来遍历已经排序的序列，将当前元素移动到对应位置。](images/placeholder.png) 图10.2是插入排序的C参考代码。外部循环标记为L1，由于单个元素$$A[0]$$已经排序，这里仅需要从元素$$A[1]$$ 迭代到元素$$A[SIZE-1]$$，其中$$SIZE$$ 表示数组元素的长度。L1的每次循环首先拷贝当前要插入有序序列中的元素$$A[i]$$，然后再执行内部L2循环寻找该值索引的合适位置。内部循环的判断条件是还没有循环到数组末尾(即条件为$$j>0$$)，并且数组元素是大于即将要插入的元素(即条件$$A[j-1]>item$$)。只要该判断条件满足，排序子数组的元素右移(操作$$A[j]=A[j-1]$$)；当判断条件不满足时，即已经查询到要插入元素的坐标位置；此时退出循环，找到索引的正确位置直接插入元素。当第​$$i$$次迭代完成后，从​$$A[0]$$到​$$A[i]$$的元素按排序顺序排列。 - 图10.2中的插入排序代码是一个简单的实现，没有做任何优化。这里可以使用Vivado HLS不同的优化指令 ，如**pipeline**, **unroll** 和**array_paration**等进行优化。最简单的优化策略即是使用**pipeline**指令，使得内部循环L2支持流水功能。对于插入排序算法，内部循环访问不同的数组元素时，没有数据相关性，因此设置期望的流水线启动间隔为1（即II=1）。生成的加速器平均需要 $$N^2/4$$次数据比较，由于每个时钟周期都需要比较一次，因此流水线的处理延迟大概需要 $$N^2/4$$个时钟周期, 参考[[58](./BIBLIOGRAPHY.md#58)]。实际上，外部循环的顺序执行的计算延迟稍高。为了获取更好的性能，可以尝试将**pipeline**指定应用到外部循环L或者函数本身，或者还可以使用部分循环展开组合。插入排序的一些优化选项可以参考表10.1。 表10.1：图10.2中插入排序的C代码实现可以尝试的优化选项 @@ -116,6 +114,7 @@ void insertion_sort(DTYPE A[SIZE]) { } } ``` + ![图10.3 ： 根据表10.1中的优化选项1重构插入排序](images/placeholder.png) 表中的优化选项2是以factor因数为2展开内部循环，即在每一个时钟周期内完成两个移位操作，可能会降低插入排序的延迟。由于每次存储访问不能分配到不同的存储区域，因此即使添加数组分区指令优化，Vivado HLS也不能达到启动间隔为1的循环。 @@ -177,8 +176,6 @@ void insertion_sort_parallel(DTYPE A[SIZE], DTYPE B[SIZE]) { ![图10.5：插入单元的体系结构。每个单元只存放一个本地值，在每次执行中，获取输入值，将其与本地值进行比较，并把较小的值写入输出端口。排序N个元素需要N个单元。](images/insertion_cell.jpg) - - 一个插入单元的核心代码如图10.6所示。输入和输出变量都声明为HLS的数据流类型。DTYPE是一个类型参数，允许对不同类型进行操作。本地变量存放阵列中的某一个元素，添加了**static**关键字是为了在多个函数调用中保存它的值。这里需要注意的是使用相同的函数，复制单元功能$$N$$次；每个单元必须有一个单独的单元静态变量，一个静态变量不能跨$$N$$个函数共享。 ```c @@ -305,7 +302,6 @@ int main () { else std::cout << "FAIL\n"; return fail; } - ``` ![图10.8：insertion_cell_sort函数的测试平台](images/placeholder.png) @@ -330,8 +326,6 @@ int main () { ![图10.10：合并两个有序数组的操作示例。顶部是初始状态，算法的每一步都需要考虑下划线的元素，并将其中的一个元素排序在输出有序数组中。](images/10.10.jpg) - - 图10.10的第一步是初始化状态，输入的两个数组各有四个元素，输出数组是空的。第二步，比较两个数组中的第一个元素，并且把最小的元素写入输出数组中（即比较3和1中小元素1并排序到数组out中）；此时数组in2的指针指向下一个元素。第三步，继续比较两个数组中当前指针所指向元素的大小，并把小的元素排序到输出数组中（即比较3和2中的小元素2并排序到数组out中）。继续类似的操作，直到其中一个数组元素移空；最后只需要拷贝未处理完的数组中剩下的元素到输出数组out中。 归并排序算法是递归函数调用的典型示例 。大多数的高层综合语言不支持递归，或者以有边界的方式支持递归。这里我们关注归并排序算法的非递归实现方式，代码看起来与通常CPU架构下的完全不同，但是算法的核心思想是完全相同的。 diff --git a/11-Huffman-Encoding.md b/11-Huffman-Encoding.md index bbb7f3b9..67d1211f 100755 --- a/11-Huffman-Encoding.md +++ b/11-Huffman-Encoding.md @@ -1,4 +1,5 @@ # 第十一章 霍夫曼编码 + ## 11.1 背景 无损数据压缩是高效数据存储中的一个关键要素，而霍夫曼编码则是其中最流行的可变长度编码算法[[33](./BIBLIOGRAPHY.md#33)]。 给定一组数据符号以及它们出现的频率，霍夫曼编码将以更短的代码分配给更频繁（出现的）符号这种方式生成码字来最小化平均代码长度。由于它保证了最优性，霍夫曼编码已被各种应用广泛采用[[25](./BIBLIOGRAPHY.md#25)]。 在现代多级压缩设计中，它经常用作系统的后端，来提高特定领域前端的压缩性能，如GZIP[[23](./BIBLIOGRAPHY.md#23)]，JPEG[[57](./BIBLIOGRAPHY.md#57)]和MP3[[59](./BIBLIOGRAPHY.md#59)]。 尽管算术编码[[61](./BIBLIOGRAPHY.md#61)]（霍夫曼编码的一个广义版本，它将整个消息转换为单个数字）可以对大多数场景实现更好的压缩，但是由于算术编码专利问题[[38](./BIBLIOGRAPHY.md#38)]，霍夫曼编码通常成为许多系统的首选算法。 @@ -7,7 +8,6 @@ Canonical霍夫曼编码与传统的霍夫曼编码相比有两大优势。 在 在基本的霍夫曼编码中，解码器通过遍历霍夫曼树，从根开始直到它到达叶节点来解压缩数据。这种方式有两个主要缺点：第一、它需要存储整个霍夫曼树，从而增加了内存使用量。此外，为每个符号遍历树的计算成本是很高的。Canonical霍夫曼编码通过使用标准的规范格式创建代码来解决这两个问题。使用Canonical编码的好处在于：我们只需要传输每个霍夫曼码字的长度。一个Canonical霍夫曼代码具有两个额外的属性：首先，较长代码比相同长度前缀的较短代码具有更高的数值。其次，具有相同长度的代码随着符号值的增加而增加。这意味着如果我们知道每个代码长度的起始符号，我们可以很容易重建Canonical霍夫曼编码。霍夫曼树本质上等同于“排序”版本的原始霍夫曼树，以便更长的码字在树的最右边的分支上，树的同一级别上的所有节点按照符号的顺序排序。 - ![图11.1：Canonical霍夫曼编码过程。 符号被过滤、排序，并且被用来建造霍夫曼树，而不是将整个树传递给解码器（就像“基本的” 霍夫曼编码所完成的那样），这样编码使得只有树的符号长度是解码器所需要的。 请注意，最终Canonical树与初始树在流程开始创建时是不同的。](images/canonical_huffman_flow.jpg) 图11.1展示了创建Canonical霍夫曼编码的过程。filter模块仅传递非零频率的符号。sort模块根据它们的频率按升序重新排列符号。接下来，create tree模块使用以下三个步骤构建霍夫曼树： @@ -97,6 +97,7 @@ dataflow指令限制函数中的信息流动。许多限制强化了子函数之 {% endhint %} ## 11.2 实现 + Canonical霍夫曼编码过程本质上是被分成子函数。因此，我们可以逐一处理这些子函数。 在我们这样做之前，我们应该考虑这些函数的每个接口。 图11.4显示了函数及其输入和输出的数据。 为了简单起见，它只显示与数组的接口，由于它们很大，我们可以假设它们存储在RAM块中（BRAMs）。 在我们描述这些函数以及他们的输入和输出之前，我们需要讨论huffman.h中定义的常量，自定义数据类型和函数接口。 图11.3显示了这个文件的内容。 @@ -144,6 +145,7 @@ void huffman_encoding ( int *num_nonzero_symbols ); ``` + 如图11.3：顶层函数huffman_encoding的参数、自定义数据类型和函数接口。 ![图11.4：Canonical霍夫曼编码的硬件实现框图。灰色块表示由不同子函数生成和消耗的重要输入和输出数据。白色块对应于函数（计算核心）。注意，数组初始位长度出现了两次（图中initial bit length模块），以使图形更加清晰。](images/che_dataflow.jpg) @@ -173,9 +175,11 @@ void filter( * n = j; } ``` + ![图11.5：filter函数在输入数组in中迭代，并将具有非零频率字段的任何符号条目添加到输出数组out。此外，它记录非零频率元素的数量并将其传递到输出参数n。](images/placeholder.png) ### 11.2.1 过滤 + 霍夫曼编码过程的第一个函数是filter，如图11.5所示。 这个函数的输入是Symbol数组。 输出是另一个Symbol数组，它是输入数组in的一个子集。Filter函数删除任何频率等于0的条目。函数本身只是简单的迭代in数组，如果其frequency域不为零，则将每个元素存储到out数组中。 另外，该函数计算输出的非零条目数，作为输出参数n传递，使其他函数只能处理“有用的”数据。 {% hint style='tip' %} @@ -183,6 +187,7 @@ Vivado HLS可以自动内联函数以便生成更高效的架构。 大多数情 {% endhint %} ### 11.2.2 分类 + 如图11.6所示，sort函数根据对输入符号的frequency值进行排序。该函数由两个for循环组成，标记为copy_in_to_sorting和radix_sort。copy_in_to_sorting循环把输入的数据从in数组搬移到sorting数组中，这确保了in数组是只读的，以满足在顶层使用的dataflow指令的要求。 sorting函数在整个执行过程中读取和写入sorting数组。即使对于这样的简单循环，使用pipeline指令生成最有效的结果和最准确的性能估计是很重要的。 radix_sort 循环实现核心基数排序算法。 通常，基数排序算法通过一次处理一个数字或一组比特来对数据进行排序。 每个数字的大小决定了排序的基数。 我们的算法在32比特的Symbol.frequency变量时处理4比特。因此我们使用基数r = 2^4 = 16排序。 对于32位数字中的每个4比特数字，我们进行计数排序。radix_sort循环执行这8个计数排序操作，以4为单位迭代到32。基数排序算法也可以从左到右（首先是最低有效位）或从右到左（首先是最高有效位）操作。算法从最低有效位运行到最高有效位。 在代码中，基数可以通过设置RADIX和BITS_PER_LOOP参数来配置。 @@ -252,6 +257,7 @@ void sort( } } ``` + ![图11.6：sort函数根据输入符号的频率值对其进行基数排序。](images/placeholder.png) {% hint style='tip' %} @@ -283,6 +289,7 @@ dataflow指令对于执行任务级别流水线优化有几个要求。其中之 {% endhint %} ### 11.2.3 创建树 + 霍夫曼编码过程中的下一个功能形成了代表霍夫曼编码的二叉树。这在图11.8所示的create_tree函数中实现。in[]包含num_symbols符号元素，按最低频率到最高频率排序。该函数创建一个这些符号的二叉树,存储在三个名为parent，left和right的输出数组中。left和right数组表示树中每个中间节点的左侧和右侧子结点。 如果子结点是一个叶结点，那么左边或右边数组的相应元素将包含孩子的符号值，否则包含特殊符号INTERNAL_NODE。同样，父数组保存每个中间结点的父结点的索引。树的根结点的父结点被定义为索引零点。树也是有序的，就像父母一样总是比孩子有更高的指数。因此，我们可以很好地实现树的自下而上和自上而下遍历。 图11.7显示了这些数据结构的一个例子。 六个符号按其频率排序并存储在数组中。 霍夫曼树的结果存储在parent，left和right三个数组中。 另外，每个中间节点的频率都存储在frequency数组中。为了便于说明，我们直接表示左右数组的节点号（例如，n0，n1等）。这些将在现实中保存一个特殊的内部结点值。 @@ -363,6 +370,7 @@ void create_tree ( parent[tree_count] = 0; //Set parent of last node (root) to 0 } ``` + ![图11.8：霍夫曼树的创建完整的代码。 该代码将已排序的符号数组in作为输入，该数组中的元素数量为n，并输出left, right和 parent三个数组的霍夫曼树。](images/placeholder.png) 主循环包含两个相似的代码块。每个代码块将[in_count].frequency中下一个可用符号的频率与下一个[tree_count]可用的的中间结点的频率相比较，然后选择两个最低频率合并为一个新的中间结点的叶子。第一个代码块对新节点的左侧子节点执行这个过程，并将结果存储在left[i]中。第二个代码块选择新节点的右侧子节点执行这个过程，并将结果存储在right[i]中。这两种情况我们都需要小心确保这个比较是有意义的。 在第一次循环迭代中，tree_count==0和i==0，所以没有有效的中间节点需要考虑，我们必须始终选择一个输入符号。在循环的最后迭代期间，所有的输入符号可能都会被使用，因此in_count==num symbols，我们必须始终使用一个中间节点。 @@ -377,6 +385,7 @@ void create_tree ( {% endhint %} ### 11.2.4 计算比特长度 + Compute_bit_length函数为每个符号确定了树的深度。深度很重要，因为它决定了用于编码的每个符号的位数。计算树中每个节点的深度是使用下面的循环完成的： $$ @@ -450,6 +459,7 @@ traverse_tree: } } ``` + ![图11.9：用于确定每个比特长度上符号数量的完整代码。](images/placeholder.png) {% hint style='tip' %} @@ -465,7 +475,7 @@ traverse_tree: 霍夫曼编码过程的下一部分是重组具有大于MAX_CODEWORD_LENGTH所指定的深度的节点。这是通过寻找任何具有更大深度的符号，并将它们移动到比所指定的目标更小的级别来完成的。有趣的是，这完全可以通过操作符号深度直方图来完成，只要直方图以与原始树上相同的模式一致的方式被修改。 输入直方图包含在input_length_histogram中，它是前面章节中compute_bit_length()函数导出的。有两个相同的输出数组truncated_length_histogram1和truncated_length histogram2。它们在稍后的过程中传递给两个单独的函数(canonize_tree and create_codewords)，因此我们必须使两个数组遵守单个生产者，单个消费者的dataow指令的约束。 -代码如图11.10所示。 copy_input循环从输入数组input_length_ histogram复制数据。 move_nodes循环包含修改直方图的大部分处理。 最后，input_length_histogram函数将内部结果复制到函数结尾的其他输出。 +代码如图11.10所示。 copy_input循环从输入数组input_length_histogram复制数据。 move_nodes循环包含修改直方图的大部分处理。 最后，input_length_histogram函数将内部结果复制到函数结尾的其他输出。 {% hint style='tip' %} copy_in for循环并未被优化。如果我们在这个循环中使用pipeline或unroll指令,延迟和启动间隔会发生什么？设计的总延时和启动间隔会发生什么变化？ @@ -525,6 +535,7 @@ void truncate_tree( } } ``` + ![图11.10：重新组织霍夫曼树的完整代码，以便任何节点的深度低于参数MAX_CODEWORD_LENGTH指定的目标值。](images/placeholder.png) 重新排序 while循环在每次迭代中移动一个节点。第一个 if语句用于找到最大深度的叶节点。然后，我们将这个结点设为中间结点，并将其添加到比目标的深度大的叶节点作为子节点。这个if从句有一个do/while循环，从目标向下迭代，寻找truncated_length_histogram数组中非零的记录。它的工作方式与开头的节点移动for循环相似。当它发现最深的叶节点小于目标深度时就停止，该节点的深度存储在j中。 @@ -534,6 +545,7 @@ void truncate_tree( 该函数是通过创建新的位长度的另一副本来完成的。这是通过存储数组truncated_length_histogram1的更新位长度到truncated_length_histogram2中完成的。我们将这两个数组传递给huffman_encoding顶层函数的最后的两个函数，我们需要这两个数组来确保满足dataflow指令的约束。 ### 11.2.6 Canonize树 + 编码过程的下一步是确定每个符号的位数，我们在图11.11所示的canoniz_tree函数中执行此操作。该函数按排序顺序的符号数组、符号总数（num_symbols）和霍夫曼树直方图的长度作为输入。symbol_bits[]包含每个符号使用的编码位数。 因此，如果具有值0x0A的符号以4位编码，则symbol_bits[10] = 4。 canonization过程由两个循环组成，分别是init_bits 和 process_symbols。 init_bits循环首先执行，初始化symbol_bits[]数组为0。然后process_symbols循环按照从最小频率到最大频率的排序顺序处理符号。通常，频率最低的符号被分配最长的码，而频率最高的符号被分配最短的码。 每次通过process_symbols循环时，我们都分配到一个符号的长度。符号的长度由内部do / while循环决定，它遍历长度直方图。该循环查找尚未分配码字的最大位长，并将该长度中的码字数存储在count中。 每次通过外循环，count递减，直到我们用完码字。 当count变为零时，内部do/while循环再次执行去寻找要分配的代码字的长度。 @@ -575,6 +587,7 @@ void canonize_tree( } } ``` + ![图11.11：完整的canonizing霍夫曼树代码，它定义了每个符号的位数](images/placeholder.png) ```c @@ -587,7 +600,8 @@ for(length = TREE_DEPTH; length >= 0; length--) { symbol_bits[sorted[k++].value] = length; } } - ``` +``` + ![图11.12：图11.11中process_symbols循环的替代循环结构。](images/placeholder.png) {% hint style='info' %} @@ -599,6 +613,7 @@ process_symbols循环无法进行流水线操作，因为内部do / while循环 {% endhint %} ### 11.2.7 创建码字 + 编码过程中的最后一步是为每个符号创建码字。 该流程根据Canonical霍夫曼代码的特性，按顺序简单地分配每个符号。第一个特性是具有相同长度前缀的较长的代码比短码有更高的数值。第二个特性是相同长度的代码随着符号值的增加而增加。为了在保持代码简单的同时实现这些属性，确定每个长度的第一个码字是有用的。如果我们知道由码字长度直方图给出的每个长度的码字的数量，则可以使用以下循环来得出结果： $$ \begin{aligned} @@ -666,6 +681,7 @@ void create_codeword( } } ``` + ![图11.13:为每个符号生成canonical霍夫曼码字的完整代码。可以在知道每个符号的比特数的情况下计算码字（存储在输入数组symbol_bits[]中）。另外，我们有另一个输入数组codeword_length_histogram[]，其在每个条目中存储具有该位长度的码字的符号数，输出是存储在encoding[]数组中的每个符号的代码字。](images/placeholder.png) 现在让我们来看看我们的运行示例，并说明如何使用它来导出最初的码字。 在这个例子中，符号A，D和E对于它们的编码有两比特; 符号C有三比特; 而符号B和F有四比特。 因此，我们有： @@ -677,7 +693,7 @@ $$ \text{bit\_length}(3) &= 1 \\ \text{bit\_length}(4) &= 2 \end{aligned} -\qquad(11.4) +\qquad(11.4) $$ 使用公式11.3来计算第一个码字的值，我们定义： @@ -706,6 +722,7 @@ $$ $$ ### 11.2.8 测试平台 + 代码的最后部分是测试平台，如图11.14所示。 总体结构是从文件中读取输入的频率值，使用huffman_encoding函数处理它们，并将结果码字与存储在文件中的现有黄金参考进行比较。 main()函数首先创建从文件中读取频率所需的变量（在这里，文件是huffman.random256.txt）并将它们放入in[]中，这是在file_to_array函数中完成的，该函数将输入数据的文件名和数据长度（数组长度）作为输入，并将该文件中的输入存储到array[]变量中。 这个文件包含了每个符号的频率， 频率值按照符号顺序存储，因此文件第一个值表示符号'0'的频率，依此类推。 @@ -779,7 +796,9 @@ int main() { return return_val; } ``` + ![图11.14：完整的canonical霍夫曼编码测试平台代码。该代码使用输入文件中的数据初始化in数组，并将它传递给顶层的huffman_encoding函数，然后它将结果码字存储到一个文件中，并将其与另一个黄金参考文件进行比较，最后打印出比较结果，并返回适当的值](images/placeholder.png) ## 11.3 结论 + 霍夫曼编码是许多应用中常用的数据压缩类型。 虽然使用霍夫曼码进行编码和解码是相对简单的操作，但生成霍夫曼码本身可能是计算上具有挑战性的问题。 在许多系统中，拥有相对较小的数据块是有利的，这意味着必须经常创建新的霍夫曼码,使之值得被加速。与我们在本书中研究过的其他算法相比，创建霍夫曼码包含许多具有完全不同代码结构的步骤。 有些相对容易实现并行化，而另一些实现起来则更具有挑战性。 算法的某些部分自然具有较高的O(n)复杂度，这意味着它们必须更加并行化以实现平衡的流水线。但是，使用Vivado HLS中的数据流指令，这些不同的代码结构可以相对容易地链接在一起。