FAE讲堂：如何加快处理器的正弦计算

只看该作者 · 2013-3-15 21:31

有很多种算法可对单精度浮点数字的正弦值进行计算，但添加硬件加速器是功能最为强大的方法之一。之所以得出这一结论，是因为客户的应用要求使用此类正弦计算，而我们又针对能够提供良好、快速且高效的解决方案进行了多种方案的探索。

为了确定哪种实现方式最适合您的应用，首先需要对代码进行分析，以查找哪种功能需要改进；其次，由于修改软件比修改硬件更简便、迅速，因而请检查是否能通过修改软件来实现您所需的高速度（有时可以）。但是如果您还需要更高的性能，那么请考虑在硬件中实现部分算法。在硬件加速的支持下，您可以轻松胜过市场上任意微控制器或DSP。

为了解该流程，让我们以现实案例为例，探讨如何开发一个需要针对单精度浮点数字进行正弦计算的军事应用。出于对高性价比的原因考虑，客户已选择了一款采用嵌入式 MicroBlaze®的Spartan®-6 FPGA 作为主系统控制器。可处理正弦计算的软件算法应运行于MicroBlaze 之上。

客户的算法主要使用浮点运算。由于算法复杂，转而采用定点运算并不妥当。此外，客户还希望避免使用定点运算时可能出现的运行过度或运行不足的情况。

客户清楚 MicroBlaze IP 可提供两种类型的浮点单元 (FPU)，并已选用扩展版本（相对于基本版而言）来加速算法。但是，这样做就无法利用作为GNU工具链组成部分且随 EDK 一起交付的数学仿真库。数学库中的软件仿真例程程序运行速度非常慢，在任何情况下都应尽量避免将其用于算法中对性能起到关键作用的部分。

另外，客户还清楚 MicroBlaze FPU的两个版本都只能处理单精度数据，不能处理双精度数据。客户的算法可以明确地仅使用浮点精度数据 (float precision data)。但在开始使用数学函数时，有时也会进行隐式转换。这些转换会强制算法
在不知不觉中使用双精度数据。

步骤一：分析问题

我们的客户已经在运行他的算法，但发现该算法在MicroBlaze处理器上的运行速度偏慢。在对代码库进行特性描述后，客户发现引起速度慢的原因是正弦计算。下一步是找出其中原因并分析怎样做才能加快处理速度。

第一种方案是使用数学库提供的标准正弦函数，在客户将算法写入后，在不进行任何修改的情况下完整地运行它。主要的问题在于数学库函数仅针对双精度数据而创建，这就意味着正弦函数的原型应为如下所示：

double sin(double angle);

但客户希望以下列方式使用：
float sin_val;
float angle;
...
sin_val = sin(angle);

当然，这也是可能的，而且C编译器会自动从参数角添加所需的转换，进行“双精度化”，并将函数调用的结果转回浮点值。这样通常还是由数学库函数来执行两个额外的转换函数，甚至是正弦计算。

切记，MicroBlaze的FPU为单精度版本，只能完成如下执行指令：

sin_val = (float)sin((double)angle);

由于数学库的正弦函数是双精度的，因而FPU无法完成正弦计算，故需要纯软件的解决方案。但缺点在于速度太慢，无法满足客户的需求。

我们验证了使用双精度数据进行正弦值的计算是执行缓慢的原因。首先我们使用下列代码，从我们的执行文件中直接创建汇编代码：

mb-objdump.exe -D executable.elf
>dump.txt

检查汇编代码时，我们发现了如下代码行：

brlid r15,-15832 // 4400d300

其作用是调用数学库以进行双精度正弦计算。然后，我们测量了利用数学库函数完成单次正弦计算所需的时间，约为 38,700个CPU周期。
对于特定的任务，可以使用专用单精度函数，如计算平方根：

float sqrt_f( float h);

使用专用函数可以避免单、双精度函数之间的转换，而且还可充分利用MicroBlaze FPU。

但遗憾的是，在FPU上没有用于处理正弦计算的专用函数。此时，我们开始开发多个版本的算法来加速正弦值的计算，以实现更高的性能。

步骤二：创建更好的软件算法

创建硬件加速器通常需要一段时间而且也需要进行调试，因而我们试图避免在第一次运行中就采取这种方案。我们就性能问题与客户进行了沟通，获得了正弦计算的关键参数。

客户的算法要求正弦计算的参数角应具有1％的精度，而且计算出的正弦值精度应比数学库函数调用的结果高0.1%。
这些属于关键参数，而且客户告知我们，他有时必须按顺序计算多个正弦值（比如在处理之前先填入小表格）。

由于对表格的尺寸要求，使用填充了所有数值的查找表显然不太可能。条目的最小数量为360,000个浮点数值（每个值 4 个字节）。客户想找到高速解决方案，但在大小上也应该合适。我们建议的解决方案可使用下列等式：

sin(xi) with xi = x + d

得到：

sin(x+d) = sin(x)*cos(d) +cos(x)*sin(d)

在这里，d是一个始终小于 x最小可能值（大于0）的值。这种解决方案有什么优势呢？我们需要缩小表格的大小，但会带来计算量的增加。表格从开始就划分为四个表格：

cos(x)
sin(x)
cos(d)
sin(d)

图1和图2显示了所有4个表格所需的分辨率以及这些值通常情况下的表现。这些表格仅显示了16个值的条目，用于说明需要填入我们的查找表中的值。我们在我们最终的解决方案中所使用的值要多得多。

图 1 — x 值的正弦与余弦表，范围介于0到360度之间

图 2 — d 值的正弦与余弦表，范围介于0到360/16度之间

实际上，我们在每个表格中都使用了1 , 0 2 4 个值。X的最小值为360/1024=0.3515625 度。d 的所有值都将小于等于该值。该方法可以减少存储的占用，因为完整的查找表需要 4,096个条目（每条目 4 个字节）。

使用这种方案，我们能够实现的自变量总体精度为：

360/(1024*1024) = 0,000343 degree

而且这个精度非常好。计算充分利用了MicroBlaze FPU。

真正的计算会占用一些时钟周期，具体来说，需要进行两次fmul运算和一次fadd运算。不过，我们还需要进行一些其它计算。首先，我们必须把自变量 xi拆分成两个值，对应x和d；然后，我们将这两个值从表格中读出；最后，我们必须使用新的算法才能计算结果。

我们在软件中实现算法并对其进行测试时，我们耗用的时钟周期总数为6,520个。

为了进一步提高分辨率，我们可以使用下列的象限关系：

第一象限

sin(x) = sin(x)

第二象限

sin(x) = sin(π - x)

第三象限

sin(x) = -sin(π + x)

第四象限：

sin(x) = -sin(2* π - x)

这在保持表格大小不变的同时还可将总体分辨率提高4倍。另一方面，我们需要进行更多的计算才能找出我们必须进行计算的象限是哪一个。仍然需要改进算法或缩小表格的大小（缩小四分之几）。我们还没有进行到这一步。

步骤三：优化算法

由于我们的解决方案到目前为止，速度还不能满足我们客户的需要，因而我们需要稍做算法优化，不过仍然完全采用运行在 MicroBlaze 处理器上的软件。这是一种简单的优化方案，不过会降低部分精度。因此，我们创建了软件模型（在PC上运行以提升运行速度）以运行所有可能的值，同时使用 sin()计算出的原始双精度值与使用我们的软件算法计算出的正弦值进行比较。我们决定在标准的PC上运行算法，因为在MicroBlaze上进行比较和计算需要花较长的时间（注意，我们的MicroBlaze运行速度远低于PC）。

现在我们开始优化计算以获得正弦值：

sin(x+d) = sin(x)*cos(d) +cos(x)*sin(d)

由于在每个表格中我们都使用了1,024个值，这意味着d始终小于360度／1,024个步进，即：

cos(2* π /1024) = 0.99998

而且该值约等于1.0。对较小的d值，适用下列等式：

cos(d) = ~1.0

这样可以将我们的公式简化为如下等式：

sin(x+d) = sin(x) + cos(x)*sin(d)

在我们在MicroBlaze上实现新等式之前，我们使用PC模式对新等式的精度进行了检验，发现最大误差仍然低于我们客户的目标。

现在我们将该算法当作软件算法在MicroBlaze上实现，仍然使用每张带有1,024个条目的表。新的算法只需要三个表，比之前的实现方案少一个。这样既节省了存储空间，也为更多的计算留出了时间。

我们在我们的硬件上测量了算法。一次正弦计算需要6,180个周期。

步骤四：进一步优化

另一种看似可行的优化方式是转换正弦计算的浮点值，并在此使用整数自变量。我们使用的算法使我们能够创建~1E6 个不同的值 (1,024*1,024)。整数自变量足以处理这个数量的值。

这种优化方式使我们能够使用简单得多的计算来将 xi 值拆分为 x 和 d。拆分只是一种简单的“与”运算加上部分10 位的移位。我们参数角的上10位是xi，下10位是 d。

我们再次在PC上创建了一个软件模型，并对其进行检验，然后在MicroBlaze处理器系统上实现模型，这需要5,460个周期才能完成一次正弦计算。

步骤五：考虑硬件实现

虽然与数学库的原始计算相比，算法的速度有了明显的改善，但客户需要的是速度快得多的实现。不过前文所述的最后一步给我们提供了一种能够轻松转向硬件实现的方法。

这种实现方法需要某些用于拆分 xi值的运算。要在硬件中做到这一点，只需将所需的位进行连接即可。然后我们需要三个表；我们使用以我们的PC模型计算出的预定义值推导出ROM，然后将其转入IP的VHDL代码中。该IP能够一次读取所有三个表，从而能够再度节省时间。最后，我们需要进行一次浮点MUL和一次浮点ADD运算。

对于该任务，我们发现用于浮点运算的CORE GeneratorTM模块非常适合。

图 3 — 无流水线功能的加速器 IP

我们使用一些Slice和乘法器，对这些硬件模块中的两个进行例化。两个内核都要求4到5个周期的延迟，以匹配我们设计的时序要求。延迟在此不是什么问题，我们将在下面的步骤中进行讨论。

我们将最终的IP以MicroBlaze的快速单工链路 (FSL) IP 的形式进行实现。对时序的第一次估算结果表明：

• 将数据从MicroBlaze传输到FSL总线需用一个时钟周期
• 将数据从FSL总线传输至FSL IP（当正弦计算的自变量从FSL总线读出时，将立即从BRAM读取数据，因而无需时钟周期）需用一个时钟周期
• 完成MUL运算 (cos(x)*sin(d)) 需用四个时钟周期
• 将方程的结果存储到寄存器中需用一个时钟周期
• 完成ADD运算需用四个时钟周期
• 将数据发送回FSL总线需用一个时钟周期
• MicroBlaze从FSL IP读取数据需用一个时钟周期。

请注意，在没有使用任何额外流水线（我们将在下一步骤中讨论这一点）的情况下，自变量数据在整个过程中必须保持稳定。这就意味着MicroBlaze仅能请求一次正弦计算，且必须读取该值，然后至少要等上13个时钟周期，才能请求下一次计算。

因此，我们估计进行该实现需要13个时钟周期。当然，要处理软件上的函数调用以及某些其他运算，还需要更多的时钟周期。

我们简单地把一些标准时钟组合在一起，不到一天就实现了该IP，随即在硬件中对该算法进行测量。整个算法（软硬件混合）耗用了360个时钟周期（包括所有的函数调用）。虽然这已是显著的进步，但是仍不足以充分满足客户的需求。

在我们的加速器IP处理所有数据之前，我们使用一个SRL16来延迟信号的写入。

虽然该算法现在可与我们的MicroBlaze并行运行，但它每次只能计算一个值。

步骤六：添加流水线和适配客户代码

设计到了这一步，我们就可以开始向我们的内核添加流水线。浮点ADD和浮点MUL的CORE Generator模块已采用流水线实现，因而我们在此无需再做什么。第一个版本的算法要求自变量保持恒定，直至计算完成。在开始新计算之前（自变量数据到达FSL IP内部），立刻读取两个BRAM并执行浮点MUL。运算的结果在数个时钟周期后生效。

我们的 sin(xi) 的自变量 xi 是一个20位宽的整数，它分为 x 和 d 两个部分。因此，我们必须对自变量 xi的MSB部分 x 进行几个时钟周期的延迟，以读取 BRAM 的内容，存储自变量xi，并将其与MUL运算的结果相匹配。

我们为我们的10位宽数值使用了少量SRL16元件（总共 10 个），共占用了10个LUT（但由于Spartan-6具有LUT组合功能，如果采用该器件较宽的LUT6结构，则仅需 5 个 LUT 即可）。

最后的工作量相当小。在图4中已对增加的SRL16x10位用红圈进行了标注。

图 4：带流水线的加速器内核

然后我们使用EDK向导来修改我们的FSL总线FIFO，以便存储多个值（我们确定能够存储8个值就足以达到我们的目的，但可根据需要轻松增加更多）。

这就意味着我们的客户甚至在请求第一个结果之前即能获得多达8个值。这足以满足我们客户当前的需求，但如果想请求更多正弦值的话，则可以轻松将FIFO缓冲参数扩展为较大的值。

我们在与客户讨论这种新的方案时，发现可将正弦计算进一步划分为两个部分：

1. 请求正弦计算（fslput 运算）
2. 请求正弦计算的结果（fslget运算）

由于我们在运算中有一个固定时延，所以如果这两个运算依次衔接、紧密地按顺序执行，那么MicroBlaze将停顿，并等待FSL IP完成对请求的处理。如果能够将这两组运算分开（这在客户的算法中是可以的），那么我们即可进一步提
升运算的总体速度。通过增加流水线，在MicroBlaze上执行的最终代码如下：

putfsl(arg1,fsl1_id);
putfsl(arg2,fsl1_id);
putfsl(arg3,fsl1_id);
putfsl(arg4,fsl1_id);
putfsl(arg5,fsl1_id);
putfsl(arg6,fsl1_id);
putfsl(arg7,fsl1_id);
putfsl(arg8,fsl1_id);
...
getfsl(result1,fsl1_id);
getfsl(result2,fsl1_id);
getfsl(result3,fsl1_id);
getfsl(result4,fsl1_id);
getfsl(result5,fsl1_id);
getfsl(result6,fsl1_id);
getfsl(result7,fsl1_id);
getfsl(result8,fsl1_id);