CN107102841A - 一种坐标变换并行计算方法和装置 - Google Patents

Abstract

本发明公开一种新型的低z通路判断，低延迟，高精度的并行CORDIC计算方法和装置，方法步骤为：(1)对输入的XIN，YIN，ZIN进行变换，使输入角度值调整到区间内，对调整后的弧度拆分成三个部分：高位，中位，低位；(2)对的高m位进行预测得出前m次旋转方向值，对前m次旋转的目标角度进行量化，产生新的旋转矩阵，将前m次旋转误差累计到低位，得到新的Z值，对剩余的角度Z分段采用不同的预测方法，可同时得出后面所有的旋转方向值；(3)对m到N/2次旋转采用传统矩阵，对N/2到N次旋转利用乘法直接求得，最后进行象限恢复输出计算得到的三角函数值。

Description

一种坐标变换并行计算方法和装置

技术领域

本发明涉及数字处理技术领域，尤其涉及基于CORDIC算法的三角函数计算与矢量变换的并行计算。

背景技术

坐标旋转数字计算机CORDIC(Coordinate Rotation Digital Computer，简称CORDIC)由J.Volder于1959年提出的，它是一种用于计算常用超越函数的循环迭代算法，其基本思想是是通过一系列的只与运算基数有关的固定小角度的不断偏摆来逼近所要旋转的角度。为了扩展可计算函数的个数，1971年J.Wahher提出了统一CORDIC算法，即把圆周坐标、线性坐标和双曲坐标统一到同一个CORDIC迭代方程中。

CORDIC算法因为能够将多种难以用硬件电路直接实现的复杂运算分解为简单的加减法和移位操作，所以很适合用数字电路来实现，因而其应用范围也就极其广泛，比如数控振荡器、正余弦函数发生器、数字频率合成器等。

下面介绍旋转模式圆周坐标下的CORDIC算法：

在平面坐标系下，将向量(x₁，y₁)旋转到如附图1的向量(x₂，y₂)，两个向量之间满足如式1所示的运算关系。

对所旋转θ角度进行分解，将其分成N个递减的小旋转角θ_i，即其中θ_i，δi为判决算子，用于确定旋转的方向。当顺时针旋转θ_i角度时，δ_i为-1；当逆时针旋转θ_i角度时，δ_i为1.每次小的角度旋转都有

再令θ_i＝tan^-12^-i，即tanθ_i＝2^-i。

这样，平面坐标系中的向量(x₁，y₁)经过N次圆周旋转之后达到同一坐标平面中的向量(x2，y2)，该旋转过程可表示为

其中K为伸缩因子

其中，K值可以提前计算出，K＝0.607529350088。

此时若令x₁＝K，y₁＝0，可得出x₂＝cosθ，y₂＝sinθ。

通过上述的理论推导，平面向量(x₁，y₁)的旋转计算问题就可以由如下基本计算式迭代N次实现。此时第i次的迭代公式就转变为下式：

由上式可知道，传统的CORDIC算法中，每一步的旋转方向都要根据上一步的结果进行判断，即δ_i＝sign(z_i)。要想得到N位的精度，需要至少迭代N次，需要N个时钟周期才能完成一次计算。为了提高运算速度，最有效的办法就是减少z通路的计算，提高运算的并行度。

Tso-Bing Juang，Shen-Fu Hsiao等人于2004年提出并行CORDIC算法。其思想是根据输入角度的n位二进制表示，将其分为高m位和低n-m位，先对高位实施Binary-to-bipolar(BBR)转换，一次产生高位的旋转次序，然后将高位产生的角度误差叠加在低位，然后在对低位采用BBR转换，一次产生低位的旋转次序。为了维持统一的矫正因子K，对于一次高位旋转2^-i的弧度，采用若干次统一方向的微旋转进行补偿(此方法称为microrotationangle recoding)，使其在第i次旋转产生的误差可以在低位进行补偿而不会产生高位的精度损失。该方法的虽然实现了CORDIC的并行处理，但增加了很多补偿的微旋转(在输入角度位数为32位时，需要额外8次微旋转)，造成X/Y路径的延迟加大。公布号为CN 102073471 A，已授权生效的发明专利“一种处理器Cordic迭代运算方法及电路”将三次迭代运算并行展开，目的是想一次计算三次迭代，但是这三次迭代的方向无法在开始迭代前确定，还是需要在计算通路中确定，实际上并没有使运算加速多少。目前大部分CORDIC算法都是采用流水线计算或者迭代结构。得到一个计算的值往往需要很多时钟周期，非常不利于CORDIC算法在FPGA上的快速实现。

发明内容

针对现有技术的缺点，提出一种新型的低z通路判断，低延迟，高精度的并行CORDIC计算方法和装置。

旋转弧度可用二进制表示为其中i∈{0,1}，θ≤π/4。

将输入的弧度拆分成三个部分：高位，中位，低位：

m是满足2^-m-tan^-12^-m<2^-N的最小值，可得出m＝[(N-log₂3)/3]。当i≥m时，tan2^-i＝2^-i。l为N/2，此时模矫正因子

根据现有技术对θ^H进行BBR转换，可得到：

根据(9)式，因为b_i-1∈{0,1}，所以r_i∈{-1,1}，由此可根据输入弧度的高m-1位直接得到1到m次旋转方向及r₁,r₂,…,r_m的值。

此时，理想的情况下1到m次的旋转的角度都为2^-i，i＝1,2,…,m。但此时 i＝1,2,…,m-1.硬件上难以实现tan 2^-i。

对弧度2^-i进行量化处理可由下式表示：

为第i次旋转的角度值。如图2所示，实线为第i次旋转目标角度，虚线为实际旋转角度，若旋转超过目标角度，角度误差-ε_i，若旋转未达到目标角度，角度误差为ε_i。

此时第i次变换矩阵变为：

σ_j的取值满足式12的前提下，取

1到m-1次旋转由式(11)替代，由式(10)可得θ^H产生的剩余角度为：

结合式(12)和式(16)可知所以σ_j满足式(12)就可以将高位产生的误差传递到低位而不会丢失精度。变换矩阵(11)硬件上可由简单的移位和相加实现。

进一步的经过前m次转换，剩余弧度为:

对θ^M′继续采用与θ^H同样的转换，得到:

其中r_m＝1-2b′_m-1,r_i＝(2b_i-1-1)，i＝m+1，m+2,…，l.

根据式(18)可得r_m,r_m+1,…,r_l次旋转的值，因为i在这个区间内，满足tan^-12^-i＝2^-i，转换矩阵变为传统矩阵如下：

经过第二阶段l-m次旋转，剩余弧度为:

在i∈[l,N]，由于模矫正因子cosθ_i＝1，可以得出若后续的某些旋转没有进行，也能保持统一的模校正因子，所以无需进行变换，此时可根据θ^L″进行简单的单向旋转，对于b_i″＝0的不进行旋转，b_i″＝1的进行单向旋转。直接根据二进制所在位的值，直接确定剩余旋转的次序。更进一步的，当i>N/2，任意的两次相邻的旋转由下式表示：

式(21)中2^-(m+n)<2^-N，已超出了最小精度，所以式(21)可简化为：

由式(22)可得N/2到N次旋转可压缩成下面矩阵表示：

式(23)中，若剩余弧度为负，则s为-1，否则为1。式(23)可由两次无符号乘法，两次加法实现。

上述中旋转方向的预测可分为两步，第一步根据θ^H得出r₁,r₂,…,r_m的值，第二步根据θ^M′+θ^L′，对θ^M′和θ^L′采用不同的预测方法，可同时得出r_m,r_m+1,…,r_l,r_l,…,r_N的值。所以整体的并行CORDIC运算只要两拍。

根据上述推导，伪旋转产生的总校正因子可以提前计算出来，不会因为输入角度不同而不同，保留了传统CORDIC算法的优点，且不需要额外的微旋转进行补偿。

附图说明：

图1是CORDIC算法旋转示意图。

图2是第i次旋转误差产生示意图。

图3是式(11)的旋转矩阵结构示意图。

图4是本实施例32位系统第一旋转矩阵硬件结构示意图。

图5是本实施例总体结构示意图。

具体实施方式：

本发明的具体实施步骤为：

(1)对输入的XIN，YIN，ZIN进行变换，使输入角度值调整到(0，π/4)区间内。对调整后的弧度拆分成三个部分：高位，中位，低位：

(2)对θ的前m-1次进行预测得出r₁,r₂,…,r_m的值，对前m-1旋转的目标角度进行量化，产生新的旋转矩阵。将前m次旋转误差累计到低位，得到新的θ^M′+θ^L′，根据θ^M′+θ^L′，对θ^M′和θ^L′采用不同的预测方法，可同时得出r_m,r_m+1,…,r_l,r_l,…,r_N的值。

(3)对m到l次旋转采用传统矩阵，对l到N次旋转利用乘法直接求得。最后进行象限恢复输出计算得到的三角函数值。

本实施例以32位的系统为例进行说明，对于一个32位系统，最大的弧度为2π，所以采用3位二进制表示整数和29位二进制表示小数。对于输入弧度不在(0，π/4)区间中的，根据现有的区间折叠技术，很容易将输入的弧度调整到(0，π/4)中，这里不在赘述。

调整后的弧度可表示为其中b_i∈{0,1}，θ≤π/4。

将输入的弧度拆分成三个部分：高位，中位，低位：

对于θ^H进行BBR转换，得到：

根据式(27)可一次得到r₁,r₂,…,r₁₀的值。

前10次的微旋转角度和矩阵为：

θ₁＝tan^-1(2^-1+2^-5+2^-6) (28)

θ₂＝tan^-1(2^-2+2^-8+2^-10) (30)

θ₃＝tan^-1(2^-3+2^-11+2^-13) (32)

θ₄＝tan^-12^-4 (34)

θ₅＝tan^-12^-5 (35)

θ₆＝tan^-12^-6 (37)

θ₇＝tan^-12^-7 (38)

θ₈＝tan^-12^-8 (39)

θ₉＝tan^-12^-9 (40)

θ₁₀＝tan^-12^-10＝2^-10 (41)

4到10次采用传统的旋转矩阵

由上述各个微旋转产生的误差为：

ε₁＝|2^-1-θ₁|＝0.00044＝o(2^-12) (42)

ε₂＝|2^-2-θ₂|＝0.00043＝o(2^-12) (43)

ε₃＝|2^-3-θ₃|＝0.00008＝o(2^-14) (44)

ε₄＝|2^-4-θ₄|＝0.000081＝o(2^-14) (45)

ε₅＝|2^-5-θ₅|＝0.00001＝o(2^-17) (46)

ε₆＝|2^-6-θ₆|＝0.000001＝o(2^-20) (47)

ε₇＝|2^-7-θ₇|＝o(2^-23) (48)

ε₈＝|2^-8-θ₈|<o(2^-23) (49)

ε₉＝|2^-9-θ₉|<o(2^-23) (50)

由θ^H产生的误差最大弧度为

预先将每一个微旋转的误差值制表保存。每一步旋转的误差值用保留进位加法器叠加产生新的

对中位θ^M′进行数学转换可得到

r₁₀＝1-2b′₉,r_i＝(2b_i-1-1)，i＝11,12,…,15 (51)

可直接得到r₁₀,r₁₁,r₁₂,r₁₃,r₁₄,r₁₅的值。对i∈[10,15]采用传统的旋转矩阵。

后15到29次旋转矩阵可由下式表示：

根据式的原码表示，若为负，s＝-1，反之，s＝1。r_i与位权2^-i所在位的值相等，r_i∈{0,1}。15到29次旋转直接由两次乘法，两次加法实现。上述中旋转方向的预测可分为两步，第一步根据θ^H得出r₁,r₂,…,r₉的值，第二步根据θ^M′+θ^L′，对θ^M′和θ^L′采用不同的预测方法，可同时得出r₁₀,…,r₁₅,r₁₅,…,r₂₉的值。

下面结合附图5进行说明：

附图5中1为前处理模块，采用区间折叠技术，对输入的XIN，YIN，ZIN进行变换，使输入角度值调整到(0，π/4)区间内，并输出象限信息。

附图5中3为第一角度预测模块，对调整后的弧度拆分成三个部分：高位，中位，低位：

对θ^H进行数学转换，重编码后一次产生r₁,r₂,…,r₁₀。各个值的确定方法如下：由于调整后的弧度一直为正，所以r₁＝1，r_i＝1-2b_i-1。确定r₁,r₂,…,r₁₀值后，结合α₁,α₂,…,α₉的值，每次旋转后产生的角度误差与θ^M+θ^L相加形成θ^M′+θ^L′。

上面提到的α₁,α₂,…,α₉由每次旋转后的误差值确定，若实际转过的角度大于目标角度，则当前次α_i＝-1，若实际转过的角度小于目标角度，则当前次α_i＝1。

附图5中2为第一旋转模块包含10次旋转。第一级旋转的角度为θ₁＝tan^-1(2^-1+2^-5+2^-6)，转换矩阵为：第二级旋转的角度为θ₂＝tan^-1(2^-2+2^-8+2^-10)，转换矩阵为第三级旋转的角度为θ₃＝tan^-1(2^-3+2^-11+2^-13)，转换矩阵为后面四级到十级旋转角度依次为θ₄＝tan^-12^-4，θ₅＝tan^-12^-5，θ₆＝tan^-12^-6，θ₇＝tan^-12^-7，θ₈＝tan^-12^-8，θ₉＝tan^-12^-9，θ₁₀＝tan^-12^-10＝2^-10，转换矩阵为传统的矩阵。r₁,r₂,…,r₁₀由第一角度预测模块得到后，开始由第一级到第十级依次计算，X/Y通路的延迟严格保持一样。在X/Y加法的计算中，根据现有技术，可选择CSA加法器合成的4-2压缩加法器加速计算，计算完的X/Y值输入第二选择模块。

附图5中5为第二角度预测模块，第一角度预测模块累加的剩余角度传入第二角度预测模块。对剩余角度的第10到第14位进行与第一预测模块高位值相同的等式变换，产生r₁₀,r₁₁,r₁₂,r₁₃,r₁₄,r₁₅的值，同时产生低N/2位的剩余角度θ^L″，根据θ^L″的值产生r₁₅,…,r₂₉的值。r₁₀＝1-2b′₉,r₁₁＝2b′₁₀-1,r₁₂＝2b′₁₁-1,r₁₃＝2b′₁₂-1,r₁₄＝2b′₁₃-1,r₁₅＝2b′₁₄-1。r₁₅,…,r₂₉与θ^L″的原码对应，r₁₅＝b″₁₅,…,，r₂₉＝b″₂₉。

附图5中4为第二旋转计算模块包括6次传统旋转计算和一次乘法计算。第二旋转模块接收第一旋转模块的X/Y的输出，接收第二角度预测模块输出的r₉,r₁₀,…,r₁₅,r₁₅,…,r₂₉的值。r₉,r₁₀,…,r₁₅用于确定第10级到第15级传统旋转的方向。r₁₅,…,r₂₉的值确定15位的乘数。

附图5中6为后处理模块，后处理模块对第二旋转计算模块的X/Y根据1中输出的象限信息进行象限的恢复。最后输出X/Y的值。

Claims (4)

1.一种坐标变换并行计算方法和装置，其特征在于包括步骤：

(1)输入角度值，对输入的XIN,YIN,ZIN进行预处理，保存象限信息，该系统有效精度为N；

(2)将处理后的弧度值根据二进制表示值，划分为三段，先对第一段进行旋转预测，确定第一段的m次旋转方向值，将第一段的弧度残余误差累加到后两段，然后根据剩余弧度预测剩余旋转方向；

(3)根据前处理模块保存的象限信息，进行象限恢复，输出三角函数值。

2.根据权利要求1所述的一种坐标变换并行计算方法和装置，其特征在于，所述步骤(2)中的第一段旋转预测和计算方法具体步骤为：

(2.1)第一段的旋转索引为1到m次，所述的m为小于[(N-log₂3)/3]的最大整数，根据弧度位权为2^-1到2^-m的值，直接得到1到m次旋转方向值；

(2.2)对于第i次目标旋转弧度2^-i，所述的i属于1到m之间，对所述弧度2^-i进行量化处理，σ_j∈{-1,0,1}，α_i∈{-1,1}，σ_j的取值满足的前提下，取

(2.3)第i次旋转矩阵由传统旋转矩阵变换为

<mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mn>1</mn> </mtd> <mtd> <mrow> <mo>-</mo> <msub> <mi>r</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>i</mi> </mrow> </msup> <mo>+</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>A</mi> </msubsup> <msub> <mi>&amp;sigma;</mi> <mi>j</mi> </msub> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>r</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>i</mi> </mrow> </msup> <mo>+</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>A</mi> </msubsup> <msub> <mi>&amp;sigma;</mi> <mi>j</mi> </msub> <msup> <mn>2</mn> <mrow> <mo>-</mo> <mi>j</mi> </mrow> </msup> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mn>1</mn> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> </mrow> σj∈{-1,0,1}，ri∈{-1,1}；

(2.4)将前m次旋转的残余误差叠加到低位弧度，形成新的弧度值。

3.根据权利要求1所述的一种坐标变换并行计算方法和装置，其特征在于，所述步骤(2)中的第二，三段旋转预测和计算方法具体步骤为：

(3.1)第二段的计算方法为根据权利要求2，2.4中所述形成的新的弧度值，直接预测m到N/2次旋转方向值，对m到N/2次旋转采用传统旋转矩阵；

(3.2)对第二段旋转产生的残余弧度误差叠加到第三段，产生新的弧度值，然后直接预测得到剩余弧度索引为N/2到N次旋转方向值。

4.根据权利要求1所述的一种坐标变换并行计算方法和装置，其特征在于，所述步骤(2)中的第三段的计算方法为：

根据权利要求3，3.2中得到的N/2到N次旋转方向值，将N/2到N次旋转矩阵压缩为：

r_i∈{0,1}，若剩余弧度为负，则s为-1，否则为1。