基于FPGA的数字图像处理【2.3】

4.2.5 浮点计算

        有些算法的映射比较麻烦，这是由于在设计的过程中涉及大量的浮点运算。在映射为FPGA逻辑时，如果都把浮点运算转换为定点运算将会带来庞大的工作量。此外，在某些应用中，定点算法是不可行的，动态范围要求使用浮点算法的一个常见的例子是矩阵求逆运算。本节将介绍如何使用FPGA来实现浮点运算以便减少移植的工作量。
1.IEEE 754标准
        浮点数（float）是属于有理数中某特定子集的数的数字表示，在计算机中用以近似表示任意某个实数。具体来说，这个实数由一个整数或定点数（即尾数）乘以某个基数（计算机中通常是2）的整数次幂得到，这种表示方法类似于基数为10的科学记数法。
        IEEE二进制浮点数算术标准（IEEE 754）是20世纪80年代以来最广泛使用的浮点数运算标准，为许多CPU与浮点运算器所采用。这个标准定义了表示浮点数的格式（包括负零-0）与反常值（denormalnumber）），一些特殊数值（无穷（Inf）与非数值（NaN）），以及这些数值的“浮点数运算符”；它也指明了四种数值舍入规则和五种例外状况（包括例外发生的时机与处理方式）。
        IEEE 754规定了四种表示浮点数值的方式：单精确度（32位）、双精确度（64位）、延伸单精确度（43比特以上，很少使用）与延伸双精确度（79比特以上，通常以80位实现）。只有32位模式有强制要求，其他都是选择性的。大部分编程语言都有提供IEEE浮点数格式与算术，但有些将其列为非必需的。例如，IEEE 754在问世之前就有的C语言，现在有包括IEEE算术，但不算作强制要求（C语言的float通常是指IEEE单精确度，而double是指双精确度）。
1）浮点格式
        二进制浮点数是以符号数值表示法的格式存储——最高有效位被指定为符号位（sign bit）；“指数部分”，即次高有效的e个比特，存储指数部分；最后剩下的f个低有效位的比特，存储“有效数”（significand）的小数部分（在非规约形式下整数部分默认为0，其他情况下一律默认为1），如图4-13所示。

        符号位sign表示数的正负（0为正，1为负）。exponent表示科学计数法的指数部分，这里所填的指数并不是前面算出来的实际指数，而 是 等 于 实 际 指 数 加 上 指 数 偏 移 ， 偏 移 量 为 2e-1 -1 ， 其 中 e 是exponent的宽度（位数）。例如，对于32位单精度浮点数，exponent宽度为8，因此偏移量为127。之后的fraction表示尾数，即科学计数法中的小数部分。
2）单精度浮点数
        单精度浮点数用32位二进制表示，其中最高位Bit[31],MSB为符号位，即sign域。Bit[30：23]为exponent域，这8位数据表示指数；最低的23位Bit[22：0]为fraction域，用于表示浮点数的小数部分。
3）双精度浮点数
        单精度浮点数用64位二进制表示，其中最高位Bit[63],MSB为符号位，即sign域。
Bit[62：52]为exponent域，这11位数据表示指数；最低的52位
Bit[51：0]为fraction域，用于表示浮点数的小数部分
4）单精度扩展浮点数
        对于单精度扩展浮点数用43～64位二进制表示，其中最高MSB为符号位，即sign域。Bit[62：52]为exponent域，这11位数据表示指数；最低的52位Bit[51：0]为fraction域，用于表示浮点数的小数部分。

2.用FPGA实现浮点运算
        当在FPGA实现图像处理计算时，尾数和指数的大小可以根据计算精度的需求进行调整，以便减少对硬件的需求。但是，这样做会使计算结果与使用IEEE标准在软件中执行同样的算法得到的结果不同。与定点数相同，这需要对近似误差进行仔细的分析，确保结果是有意义的。
        对于浮点数来说，乘法和除法是相对比较简单的操作。对于乘法运算、位数相乘及指数相加，若位数结果大于2，则需要重新规范化，将其右移一位并且增加指数。由于成绩的位数可能会比表示位多，因此，需要进行舍位处理。两个指数的和包括两个偏移量，因此必须减掉一个。输出值的符号位是输入符号位的异或。需要附加的逻辑来检测下溢出、上移除及处理其他的错误情况，例如，处理无穷大和非数。除法与惩罚类似，只是尾数相除，指数相减并且在重新规范化时可能包含左移。
        加法和减法的实现要更加复杂。与原码表示相同，实际操作的执行取决于输入的符号。指数必须相同，因为数必须要对齐。根据指数位的差值，将较小的指数右移相应的位数。在FPGA上实现时，一个这样的移位或者较慢（用一些列较小的移位器实现）或者价格比较昂贵（用许多宽的多路复用器实现）。有必要为移位的数保留一个额外的位来减少操作引起的误差。然后根据操作对尾数进行加法或减法运算并对齐进行重新规范化。如果两个非常相近的数相减，许多高有效位可能会消失，因此需要确定最左边的1的位置，将其移动到最高位来重新规范化，并对指数进行相应的调整。

        浮点操作相对于定点操作无疑要消耗更多的资源，这不是因为浮点操作有多复杂，而是因为处理异常时需要很多逻辑，尤其是在需要符合IEEE标准时。大部分FPGA在进行浮点运算时，为符合IEEE 754标准，每次运算都需要去归一化和归一化步骤，导致了极大的性能瓶颈。因为这些归一化和去归一化步骤一般通过FPGA中的大规模桶形移位寄存器实现，需要大量的逻辑和布线资源。通常一个单精度浮点加法器需要500个查找表(LUT)，单精度浮点要占用30%的LUT，指数和自
然对数等更复杂的数学函数需要大约1000个LUT。因此，随着DSP算法越来越复杂，FPGA性能会明显劣化，对占用80%～90%逻辑资源的FPGA会造成严重的布线拥塞，阻碍FPGA的快速互连，最终会影响时序收敛。图4-14为传统的计算浮点计算方法。
3.Altera的浮点IP核
        针对浮点计算的这个缺点，FPGA厂家也提出了各种解决方案。Altera的DSP Builder高级模块库引入了融合数据通路设计。它将基本算子组合在一个函数或者数据通路中，通过分析数据通路的位增长，选择最优归一化输入，为数据通路分配足够的精度，尽可能消除归一化和去归一化步骤。这一优化平台将定点DSP模块与可编程软核逻辑相结合，避免了大量使用这类桶形移位寄存器。与使用几种基本IEEE754算子构成的等价数据通路相比，减少了50%的逻辑，延时减小了50%。并且这一方法总的数据精度一般高于使用基本IEEE 754浮点算子库的方法。

        Altera在其Arria 10和Stratix 10系列器件中增加了硬核浮点DSP模块。浮点模式的一个DSP模块提供了IEEE 754单精度浮点乘法器及IEEE 754单精度加法器，可以实现单精度的乘法、乘加、乘减、累加、支持浮点矢量运算、卷积、点乘和其他线性算术函数，以及使用快速傅里叶变换的复数乘法等。硬浮点模块的出现大大降低了浮点运算的资源消耗，可以大大减轻设计人员对复杂耗时的算法的转换工作。
        Altera现在提供业界最全面的单精度和双精度浮点IP内核，其性能非常高。目前提供的浮点IP内核包括以下几项：
（1）加减乘除运算。
（2）倒数。
（3）指数与指数。
（4）三角函数。
（5）平方根与逆平方根。
（6）矩阵乘法与求逆。
（7）快速傅里叶变换（FFT）。
（8）定点与浮点转换。
        altera的浮点IP输入和输出均为浮点数，并且支持32位单精度浮点、64位双精度浮点及自定义扩展性单精度浮点类型。在使用时首先需将定点数或整形数转换为浮点类型，调用浮点IP计算后，有必要时再将结果转换为整形，计算步骤如图4-15所示。

        这里例化以下几个浮点运算IP核：sin，cos，sqrt，exp分别求取正弦函数，余弦函数，求根号函数，指数函数。输入格式为定点小数 ： 整 数 位 8 位 ， 小 数 位 8 位 。 所 调 用 的 IP 核 分 别 为ALTFP_CONVERT,ALTFP_CONVERT,ALTFP_SINCOS,ALTFP_SQRT,ALTFP_EXP。
设计Testbench如下：

'timescale 1 ps / 1 psmodule float_tb;
parameter DW = 16;
reg [DW-1:0]din;
reg clock;
reg rst;
wire [DW-1:0]sin_dout;
wire [DW-1:0]cos_dout;
wire [DW-1:0]sqrt_dout;
wire [DW-1:0]exp_dout;
wire [32-1:0]float_temp;
wire [32-1:0]sin_temp;
wire [32-1:0]cos_temp;
wire [32-1:0]sqrt_temp;
wire [32-1:0]exp_temp;
initial
begin
clock ＜= 1'b0;
din ＜= {DW{1'b0}};
rst ＜= 1'b1;
#12003;
rst ＜= 1'b0;
end
//将输入数据转换为浮点
int_to_float u0(
.aclr(rst),
.clock(clock),
.dataa(din),.result(float_temp)
);
//求正弦
sin sin_ins(
.aclr(rst),
.clock(clock),
.data(float_temp),
.result(sin_temp)
);
//求余弦
cos cos_ins(
.aclr(rst),
.clock(clock),
.data(float_temp),
.result(cos_temp)
);
//求根号
sqrt sqrt_ins(
.aclr(rst),
.clock(clock),
.data(float_temp),
.result(sqrt_temp)
);
//求指数，以e为底
exp_fp exp_ins(
.aclr(rst),
.clock(clock),.data(float_temp),
.result(exp_temp)
);
//将结果转换为整形
float_to_int u2(
.aclr(rst),
.clock(clock),
.dataa(sin_temp),
.result(sin_dout)
);
float_to_int u3(
.aclr(rst),
.clock(clock),
.dataa(cos_temp),
.result(cos_dout)
);
float_to_int u4(
.aclr(rst),
.clock(clock),
.dataa(sqrt_temp),
.result(sqrt_dout)
);
float_to_int u5(
.aclr(rst),
.clock(clock),
.dataa(exp_temp),
.result(exp_dout));
always @(clock)
clock ＜= #50 ～clock;
always @(posedge clock)
din＜=din+1'b1;
endmodule

        仿真结果如图4-16所示。仿真结果与实际的函数波形是一致的，需要注意的是越界问题。
        如果选择内部没有浮点硬核的FPGA，那么系统会采用内部资源来实现浮点运算，单精度正弦和余弦函数的资源消耗情况分别如图4-17和图4-18所示。

        由资源消耗情况可见（左下角），浮点运算所消耗的逻辑资源是相当惊人的。这也是为什么在一般情况下宁愿采用定点计算的原因。需要注意的是，一个细节是余弦函数明显比正弦函数消耗的资源小，如果必须使用浮点三角运算，那么可以考虑将正弦函数转换成余弦进行计算，这样可以节省一定的资源。

4.2.6 Cordic技术

        Cordic （ Coordinate Rotation Digital Computer；Volder ，1959）是用于计算初等函数的一种迭代方法。最初被设计用于计算初等三角函数（Volder，1959），但后来被推广到包括双曲函数、乘法和除法操作上（Walther，1971）。Cordic计算的基本原理是基于向量旋转，考虑如图4-19所示的向量变换。

        对于图4-19所示的旋转，可以通过多次旋转来达到目的，这个多次旋转的过程也被称为迭代旋转。我们设定初始化坐标为（X0 ,Y0），经过k次迭代旋转到最终的目标向量，i为第i次旋转，每次的旋转角度为θi 。则每次迭代的数学模型为

算法的技巧在于选择合适的旋转角度，并重新整理上式使矩阵内的因数是2的负次幂，于是乘法操作可以通过简单的移位来实现：

其中，di 是旋转的方向，+1代表逆时针旋转，-1代表顺时针旋转，而角度由以下式给出：

角度综合取决于每次迭代旋转的方向，因此，增加一个寄存器（标记为z）用于累加角度是十分有用的。

有两种模式可以用于选择每部的旋转方向：旋转模式和向量化模式。
旋转模式始于Z寄存器中的一个角度，选择与Z符号相同的di ，随
着向量的旋转，Z逐渐减小直至为0：
di ＝sign(Zi )
k次迭代后的结果为

换句话说，如果当前角度为0，经过一个旋转移动后，那么剩余的旋转必须能再次将该角度移回到0。在这种情况下，该级数收敛，因为

        该模式与旋转模式有同样的收敛域，因此，若X0 ﹤0，则有必要使该范围内的向量旋转180°。向量化模式相当于将向量从笛卡尔矩形坐标系转换为极坐标系。如果X0 和Y0 都非常小，那么舍入误差就会非常大，尤其是在计算反正切时（Kota和Cavallaro），如有必要，可以通过使用2的幂来缩放Xi 和Yi 使其规范化来解决这个问题。
        对于旋转模式和引导模式来说，每部迭代都把角度的精度提高了接近一个二进制位。Andraka（1998）提出了另外一种操作模式用于计算反正弦和反余弦。开始时，X0 ＝Y0 ＝ 1，按照如下方式进行旋转：

这里假设角度用弧度表示。要注意这种方法并不适用于向量化模式，因为这种方法要求角度趋向于0。