FPGA上定点atan2的计算方法

我需要在具有连续输入/输出数据流的FPGA上计算atan2(x,y)。我设法使用展开的流水线CORDIC内核来实现它，但是要获得所需的精度，我必须执行32次迭代。这导致大量的LUT致力于这一任务。我尝试将流更改为使用部分展开的CORDIC内核，但随后我需要一个倍频的时钟频率以执行重复循环，同时仍保持连续的输入/输出流。有了这个，我无法满足时间要求。

所以现在我正在寻求替代的atan2(x,y)计算方法。

我考虑过使用带有插值的Block-RAM查找表，但是由于有2个变量，因此我需要2维查找表，这在Block-RAM的使用方面非常耗费资源。

然后我考虑使用象限调整将atan2(x,y)与atan(x/y)相关的事实。问题是x/y需要一个真正的除法，因为y不是常数，并且FPGA上的除法非常耗费资源。

是否还有更多新颖的方法可以在FPGA上实现atan2(x,y)在较低的LUT使用率下，但仍提供了较高的精度？

您的处理时钟速率和输入数据速率是多少？

您要求的精度是多少？我还假设您正在使用定点计算。您正在使用什么深度？具有象限调整的多项式逼近（或LUT）是实现atan2的常用方法。不过，不确定是否可以不分手而行。

输入时钟为150MHz，输入数据速率为150 MSamps / sec。基本上，每个时钟周期我都会得到一个新的输入。有延迟是可以的，但我也必须以150 MSamps / sec的速度产生输出。

我的模拟显示我可以忍受大约1 * 10 ^ -9。不确定绝对最小定点位数，但我一直在使用Q10.32定点格式进行仿真

本文介绍了atan2的定点实现。不过，您仍然需要一个部门。

#1 楼

您可以使用对数摆脱除法。对于第一象限中的$（x，y）$：

$$ z = \ log_2（y）-\ log_2（x）\\
\ text {atan2}（y， x）= \ text {atan}（y / x）= \ text {atan}（2 ^ z）$$

图1. $ \的图text {atan}（2 ^ z）$

您需要在$ -30
$$ b = \ text {floor}（\ log_2（a））\\
c = \ frac {a} {2 ^ b} \\
\ log_2（a）= b + \ log_2（c）$$

$ b $可以通过查找最高有效非零位的位置来计算。 $ c $可以通过位移来计算。您将需要在$ 1 \ le c <2 $范围内近似$ \ log_2（c）$。

图2. $ \ log_2（c）的图$

对于您的精度要求，线性插值和统一采样，$ 2 ^ {14} + 1 = 16385 $个$ \ log_2（c）$采样和$ 30 \ times 2 ^ {12} + 1 = $ 0
$ \ text {atan}（2 ^ z）$表可以拆分为多个子表，分别对应于$ 0 \ le z <1 $和$ z \不同的$ \ text {floor}（\ log_2（z））$。 ge 1 $，这很容易计算。表长度可以选择，如图3所示。子表内索引可以通过简单的位串操作来计算。为了您的准确性要求，如果您将$ z $的范围扩展到$ 0 \ le z <32 $，为简单起见，$ \ text {atan}（2 ^ z）$子表将总共有29217个样本。

为了便于以后参考，这是我用来计算近似误差的笨拙的Python脚本：

from numpy import *
from math import *
N = 10
M = 20
x = array(range(N + 1))/double(N) + 1
y = empty(N + 1, double)
for i in range(N + 1):
    y[i] = log(x[i], 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y[i] + (y[i + 1] - y[i])*j/M
        if N*M < 1000: 
            print str((i*M + j)/double(N*M) + 1) + ' ' + str(a)
        b = log((i*M + j)/double(N*M) + 1, 2)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2 = empty(N + 1, double)
for i in range(1, N):
    y2[i] = -1.0/16.0*y[i-1] + 9.0/8.0*y[i] - 1.0/16.0*y[i+1]


y2[0] = -1.0/16.0*log(-1.0/N + 1, 2) + 9.0/8.0*y[0] - 1.0/16.0*y[1]
y2[N] = -1.0/16.0*y[N-1] + 9.0/8.0*y[N] - 1.0/16.0*log((N+1.0)/N + 1, 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000: 
            print a
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2[0] = 15.0/16.0*y[0] + 1.0/8.0*y[1] - 1.0/16.0*y[2]
y2[N] = -1.0/16.0*y[N - 2] + 1.0/8.0*y[N - 1] + 15.0/16.0*y[N]

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000: 
            print str(a) + ' ' + str(b)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

P = 32
NN = 13
M = 8
for k in range(NN):
    N = 2**k
    x = array(range(N*P + 1))/double(N)
    y = empty((N*P + 1, NN), double)
    maxErr = zeros(P)
    for i in range(N*P + 1):
        y[i] = atan(2**x[i])

    for i in range(N*P):
        for j in range(M):
            a = y[i] + (y[i + 1] - y[i])*j/M
            b = atan(2**((i*M + j)/double(N*M)))
            err = abs(a - b)
            if (i*M + j > 0 and err > maxErr[int(i/N)]):
                maxErr[int(i/N)] = err

    print N
    for i in range(P):
        print str(i) + " " + str(maxErr[i])

近似函数$ f（通过线性插值$ f（x）$的样本中的$ \ hat {f}（x）$（采用采样间隔为$ \ Delta x $的均匀采样所获得的值），可以通过以下分析近似得出：

$$ \ widehat {f}（x）-f（x）\ approx（\ Delta x）^ 2 \ lim _ {\ Delta x \ rightarrow 0} \ frac {\ frac {f（x）+ f（ x + \ Delta x）} {2}-f（x + \ frac {\ Delta x} {2}）} {（\ Delta x）^ 2} = \ frac {（\ Delta x）^ 2 f'' （x）} {8}，$$

其中$ f''（x）$是$ f（x）$的二阶导数，$ x $是绝对值的局部最大值错误。通过以上方法我们可以得出近似值：

$$ \ widehat {\ text {atan}}（2 ^ z）-\ text {atan}（2 ^ z）\ approx \ frac {（\ Delta z）^ 2 2 ^ z（1-4 ^ z）\ ln（2）^ 2} {8（4 ^ z + 1）^ 2}，\\
\ widehat {\ log_2}（a ）-\ log_2（a）\ approx \ frac {-（\ Delta a）^ 2} {8 a ^ 2 \ ln（2）}。$$

因为函数是凹的，所以样本与函数匹配，误差始终是一个方向。如果使误差的符号在每个采样间隔内来回交替一次，则可以将局部最大绝对误差减半。使用线性插值，可以通过以下方式对每个表进行预过滤来获得接近最佳的结果：

$$ y [k] = \ cases {\ begin {array} {rrrrrl} && b_0x [k]＆\ negthickspace \ negthickspace \ negthickspace + b_1x [k + 1]＆\ negthickspace \ negthickspace \ negthickspace + b_2x [k + 2]＆\ text {if} k = 0，\\
＆c_1x [k-1]＆\ negthickspace \ negthickspace \ negthickspace + c_0x [k]＆\ negthickspace \ negthickspace \ negthickspace + c_1x [k + 1] && \ text {if} 0 b_2x [ k-2]＆\ negthickspace \ negthickspace \ negthickspace + b_1x [k-1]＆\ negthickspace \ negthickspace \ negthickspace + b_0x [k] &&& \ text {if} k = N，
\ end {array}} $$

其中$ x $和$ y $是原始表，过滤后的表都跨越$ 0 \ le k \ le N $，权重为$ c_0 = \ frac {9} {8}，c_1 =-\ frac {1} {16}，b_0 = \ frac {15} {16}，b_1 = \ frac {1} {8}，b_2 =-\ frac {1} {16} $。与使用表外函数的样本相比，结束条件（上式中的第一行和最后一行）减少了表末尾的错误，因为不需要调整第一个和最后一个样本以减少内插误差它和桌子外的样品之间。具有不同采样间隔的子表应单独进行预过滤。权重$ c_0，c_1 $的值是通过按顺序最小化以增加指数$ N $的近似误差的最大绝对值而得出的：

$$（\ Delta x）^ N \ lim_ { \ Delta x \ rightarrow 0} \ frac {\ left（c_1f（x-\ Delta x）+ c_0f（x）+ c_1f（x + \ Delta x）\ right）（1-a）+ \ left（c_1f（x ）+ c_0f（x + \ Delta x）+ c_1f（x + 2 \ Delta x）\ right）a-f（x + a \ Delta x）} {（\ Delta x）^ N}
= \ left \ {\ begin {array} {ll}（c_0 + 2c_1-1）f（x）＆\ text {if} N = 0，\ bigg | c_1 = \ frac {1-c_0} {2} \\
0＆\文本{如果} N = 1，\\
\ frac {1 + aa ^ 2-c_0} {2}（\ Delta x）^ 2 f''（x）＆\ text {if} N = 2，\ bigg | c_0 = \ frac {9} {8} \ end {array} \ right。$$

对于样本间插值位置$ 0 \ le a <1 $，具有凹函数或凸函数$ f（x）$（例如$ f（x）= e ^ x $）。解决了这些权重后，通过类似地最小化以下项的最大绝对值来找到最终条件权重$ b_0，b_1，b_2 $的值：

$$（\ Delta x）^ N \ lim _ {\ Delta x \ rightarrow 0} \ frac {\ left（b_0f（x）+ b_1f（x + \ Delta x）+ b_2f（x + 2 \ Delta x）\ right ）（1-a）+ \ left（c_1f（x）+ c_0f（x + \ Delta x）+ c_1f（x + 2 \ Delta x）\ right）a-f（x + a \ Delta x）} {（ \ Delta x）^ N}
= \ left \ {\ begin {array} {ll} \ left（b_0 + b_1 + b_2-1 + a（1-b_0-b_1-b_2）\ right）f（ x）＆\ text {if} N = 0，\ bigg | b_2 = 1-b_0-b_1 \\
（a-1）（2b_0 + b_1-2）\ Delta x f'（x）＆\ text {if} N = 1，\ bigg | b_1 = 2- 2b_0 \\
\ left（-\ frac {1} {2} a ^ 2 + \ left（\ frac {23} {16}-b_0 \ right）a + b_0-1 \ right）（\ Delta x）^ 2f''（x）＆\ text {if} N = 2，\ bigg | b_0 = \ frac {15} {16} \ end {array} \ right。$$

$ 0 \ le a <1 $。使用预过滤器可使近似误差减半，并且比对表进行全面优化要容易。

图4. $ \ log_2（a）的近似误差来自11个样本的$，带和不带预过滤器，带和不带末端调节器。在没有结束条件的情况下，预过滤器可以访问表外的函数值。

本文可能提出一种非常相似的算法：R. Gutierrez，V。Torres和J. Valls，“ FPGA对数转换和基于LUT的技术对atan（Y / X）的实现”，《系统结构学报》，第1卷。 2010年5月56日。摘要表示，它们的实现在速度方面击败了以前的基于CORDIC的算法，而在占用空间方面又击败了基于LUT的算法。

$ \ begingroup $
我和Matthew Gambrell对1985 Yamaha YM3812声音芯片（通过显微镜）进行了反向工程，并在其中发现了类似的日志/扩展只读存储器（ROM）表。 Yamaha使用了另一个技巧，即将每个表中的第二个条目替换为与前一个条目的不同之处。对于平稳的功能，差异代表的位数和芯片面积少于函数。他们已经在芯片上有一个加法器，可以用来将差异添加到上一个条目。
$ \ endgroup $
–奥利·尼米塔洛（Olli Niemitalo）
16-2-16在6:51

$ \ begingroup $
非常感谢！我喜欢数学性质的这些利用。我一定会对此进行一些MATLAB模拟，如果一切看起来不错，请转到HDL。完成所有步骤后，我将报告我的LUT节省。
$ \ endgroup $
–user2913869
16年2月16日在17:32

$ \ begingroup $
我以您的描述为指导，我很高兴能将LUT减少了近60％。我确实需要减少BRAM，因此我发现我可以通过进行非均匀采样来在ATAN表上获得一致的最大错误：我有多个LUT BRAM（都具有相同数量的地址位），越接近零，采样越快。我将表范围选择为2的幂，因此我可以轻松地检测到自己所在的范围，并通过位操作进行自动表索引。我也应用了atan对称性，所以我只存储了一半的波形。
$ \ endgroup $
–user2913869
16-2-19的2:55

$ \ begingroup $
另外，我可能已经错过了您的一些修改，但是我设法通过将其拆分为2 ^ {if} = 2 ^ i * 2 ^ {0.f}来实现2 ^ z，其中i是整数f是小数部分。 2 ^ i很简单，只需进行位操作即可，而2 ^ {0.f}的范围有限，因此它很适合插值使用的LUT。我还处理了否定的情况：2 ^ {-if} = 2 ^ {-i} * 1 /（2 ^ {0.f}。因此，又有一张表格用于1/2 ^ {0.f}。我的下一步可能是对log2（y）LUT施加2个测距/非均匀采样的幂，因为这似乎将是这种事情的理想候选波形。
$ \ endgroup $
–user2913869
16-2-19的3:01

$ \ begingroup $
大声笑，我完全错过了这一步。我现在要尝试。应该为我节省更多的LUT和更多的BRAM
$ \ endgroup $
–user2913869
16-2-19在14:46

#2 楼

MATLAB中没有乘法算法：
（x，y）的32 / Pi X反正切使用（x或y）的4个最高有效位

%16x16 quadrant look-up table:

ATAN2(1,1:16)= fi([ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 ],0,4,0);

ATAN2(2,1:16)= fi([ 0  8 11 13 13 14 14 14 15 15 15 15 15 15 15 15 ],0,4,0);

ATAN2(3,1:16)= fi([ 0  5  8 10 11 12 13 13 13 14 14 14 14 14 14 14 ],0,4,0);

ATAN2(4,1:16)= fi([ 0  3  6  8  9 10 11 12 12 13 13 13 13 13 14 14 ],0,4,0);

ATAN2(5,1:16)= fi([ 0  2  5  6  8  9 10 11 11 12 12 12 13 13 13 13 ],0,4,0);

ATAN2(6,1:16)= fi([ 0  2  4  5  7  8  9 10 10 11 11 11 12 12 12 13 ],0,4,0);

ATAN2(7,1:16)= fi([ 0  2  3  5  6  7  8  9  9 10 10 11 11 11 12 12 ],0,4,0);

ATAN2(8,1:16)= fi([ 0  1  3  4  5  6  7  8  9  9 10 10 10 11 11 11 ],0,4,0);

ATAN2(9,1:16)= fi([ 0  1  2  4  5  6  6  7  8  8  9  9 10 10 11 11 ],0,4,0);

ATAN2(10,1:16)=fi([ 0  1  2  3  4  5  6  7  7  8  8  9  9 10 10 10 ],0,4,0);

ATAN2(11,1:16)=fi([ 0  1  2  3  4  5  5  6  7  7  8  8  9  9 10 10 ],0,4,0);

ATAN2(12,1:16)=fi([ 0  1  2  3  3  4  5  6  6  7  7  8  8  9  9  9 ],0,4,0);

ATAN2(13,1:16)=fi([ 0  1  2  2  3  4  5  5  6  6  7  7  8  8  9  9 ],0,4,0);

ATAN2(14,1:16)=fi([ 0  1  2  2  3  4  4  5  6  6  7  7  7  8  8  9 ],0,4,0);

ATAN2(15,1:16)=fi([ 0  1  1  2  3  3  4  5  5  6  6  7  7  8  8  8 ],0,4,0);

ATAN2(16,1:16)=fi([ 0  1  1  2  3  3  4  4  5  5  6  6  7  7  8  8 ],0,4,0);

CMP7=fi([0 127 0 0],0,7,0); % using look-up instead if if-end 

CMPA=fi([31 0 31 0],0,5,0); % using look-up instead if if-end
 

X=int32(x);

Y=int32(y);

XsN=bitget(X,32);

YsN=bitget(Y,32);

XSM=uint16(bitshift(XsN+bitxor(CMP(1+XsN),X),-16));

YSM=uint16(bitshift(YsN+bitxor(CMP(1+YsN),Y),-16));

XYSM=bitor(bitor(XSM,YSM),uint16(15));

XYSM=bitset(XYSM,15,or(bitget(XYSM,16),bitget(XYSM,15)));

XYSM=bitset(XYSM,14,or(bitget(XYSM,15),bitget(XYSM,14)));

XYSM=bitset(XYSM,13,or(bitget(XYSM,14),bitget(XYSM,13)));

XYSM=bitset(XYSM,12,or(bitget(XYSM,13),bitget(XYSM,12)));

XYSM=bitset(XYSM,11,or(bitget(XYSM,12),bitget(XYSM,11)));

XYSM=bitset(XYSM,10,or(bitget(XYSM,11),bitget(XYSM,10)));

XYSM=bitset(XYSM,9,or(bitget(XYSM,10),bitget(XYSM,9)));

XYSM=bitset(XYSM,8,or(bitget(XYSM,9),bitget(XYSM,8)));

XYSM=bitset(XYSM,7,or(bitget(XYSM,8),bitget(XYSM,7)));

XYSM=bitset(XYSM,6,or(bitget(XYSM,7),bitget(XYSM,6)));

XYSM=bitset(XYSM,5,or(bitget(XYSM,6),bitget(XYSM,5)));

XYSM=bitset(XYSM,4,or(bitget(XYSM,5),bitget(XYSM,4)));

XYSM=bitset(XYSM,3,or(bitget(XYSM,4),bitget(XYSM,3)));

XYSM=bitset(XYSM,2,or(bitget(XYSM,3),bitget(XYSM,2)));

XYSM=bitset(XYSM,1,or(bitget(XYSM,2),bitget(XYSM,1)));

XYSh=bitget(XYSM,1)+bitget(XYSM,2)+bitget(XYSM,3)+bitget(XYSM,4)+bitget(XYSM,5)+bitget(XYSM,6)+bitget(XYSM,7)+bitget(XYSM,8)+bitget(XYSM,9)+bitget(XYSM,10)+bitget(XYSM,11)+bitget(XYSM,12)+bitget(XYSM,13)+bitget(XYSM,14)+bitget(XYSM,15)+bitget(XYSM,16);

XS = fi(0,0,5,0);

YS = fi(0,0,5,0);

XS = bitset(XS,4,bitget(XSM,XYSh));

XS = bitset(XS,3,bitget(XSM,XYSh-1));

XS = bitset(XS,2,bitget(XSM,XYSh-2));

XS = bitset(XS,1,bitget(XSM,XYSh-3));



YS = bitset(YS,4,bitget(YSM,XYSh));

YS = bitset(YS,3,bitget(YSM,XYSh-1));

YS = bitset(YS,2,bitget(YSM,XYSh-2));

YS = bitset(YS,1,bitget(YSM,XYSh-3));



XsZ = XS==0;
XsP = ~or(XsZ,XsN);
XsN = ~or(XsZ,XsP);

YsZ = YS==0;
YsP = ~or(YsZ,YsN);
YsN = ~or(YsZ,YsP);

Quad2 = or(and(XsZ,YsP),and(XsN,YsP));

Quad3 = or(and(XsN,YsZ),and(XsN,YsN));

Quad4 = or(and(XsZ,YsN),and(XsP,YsN));

Quadf = or(Quad4,Quad2)+1;

Quad=bitand(CMP7(Quad2+1),fi(16,0,7,0))+bitand(CMP7(Quad3+1),fi(32,0,7,0))+bitand(CMP7(Quad4+1),fi(48,0,7,0));

AnswerX32overPi=fi(ATAN2(1+bitand(CMPA(Quadf),XS)+bitand(CMPA(Quadf+1),YS)+bitshift(fi(bitand(CMPA(Quadf),YS)+bitand(CMPA(Quadf+1),XS),0,8,0),4))+Quad,1,7,0);

编程黑洞网

FPGA上定点atan2的计算方法

评论

#1 楼

评论

#2 楼