CN107102841A - A kind of coordinate transform parallel calculating method and device - Google Patents

Abstract

The present invention discloses a kind of new low z paths and judged, low latency, high-precision Parallel CORDIC computational methods and device, and method and step is：(1) to the XIN of input, YIN, ZIN enters line translation, input angle angle value is adjusted in interval, and three parts are split into the radian after adjustment：A high position, middle position, low level；(2) to high m be predicted draw before m rotational direction value, the angle on target of first m times rotation is quantified, produce new spin matrix, preceding m rotation error is accumulated to low level, obtain new Z values, different Forecasting Methodologies is used to remaining angle Z segmentations, rotational direction value all below can be drawn simultaneously；(3) classical matrix is used to m to N/2 rotations, N/2 is directly tried to achieve to n times rotation using multiplication, finally carried out quadrant recovery output and calculate obtained trigonometric function value.

Description

A kind of coordinate transform parallel calculating method and device

Technical field

The present invention relates to digital processing technology field, more particularly to the trigonometric function calculating based on cordic algorithm and vector The parallel computation of conversion.

Background technology

CORDIC CORDIC (Coordinate Rotation Digital Computer, referred to as CORDIC) proposed by J.Volder in nineteen fifty-nine, it be it is a kind of be used to calculate the conventional loop iteration algorithm surmounted function, its Basic thought is to be rotated by the continuous beat of a series of fixation low-angle relevant with computing radix to approach Angle.In order to extend the number of calculable functions, J.Wahher in 1971 proposes unified cordic algorithm, i.e., circumference sat Mark, linear coordinate and hyperbolic coordinates are unified into same CORDIC iterative equations.

A variety of complex calculations for being difficult to be directly realized by with hardware circuit because can be decomposed into simply by cordic algorithm Addition and subtraction and shifting function, so be well suited for being realized with digital circuit, thus its application is also just extremely extensive, than in full Controlled oscillator, sin cos functionses generator, digital frequency synthesizer etc..

The cordic algorithm under rotary mode circumferential coordinates is described below：

Under plane coordinate system, by vector (x₁, y₁) rotate to vector (x such as accompanying drawing 1₂, y₂), it is full between two vectors Foot operation relation as shown in Equation 1.

Rotated θ angles are decomposed, N number of small rotation angle θ successively decreased is divided into_i, i.e.,Wherein θ_i, δ i are judgement operator, the direction for determining rotation.As the θ that turns clockwise_iDuring angle, δ_iFor -1；As rotate counterclockwise θ_iAngle When, δ_iHave for 1. each small angle rotations

θ is made again_i=tan^-12^-i, i.e. tan θ_i=2^-i。

So, the vector (x in plane coordinate system₁, y₁) reached after the rotation of n times circumference in same coordinate plane Vectorial (x2, y2), the rotary course is represented by

Wherein K is contraction-expansion factor

Wherein, K values can be calculated in advance, K=0.607529350088.

If now making x₁=K, y₁=0, x can be drawn₂=cos θ, y₂=sin θ.

Pass through above-mentioned theory deduction, plane vector (x₁, y₁) rotation computational problem just can be by following basic calculating formula Iteration n times are realized.Now the iterative formula of ith is just changed into following formula：

It may know that by above formula, in traditional cordic algorithm, the direction of rotation of each step will be according to the result of previous step Judged, i.e. δ_i=sign (z_i).Expect the precision of N, it is necessary at least iteration n times, it is necessary to N number of clock cycle could Completion is once calculated.In order to improve arithmetic speed, maximally effective method is exactly to reduce the calculating of z paths, improves the parallel of computing Degree.

Tso-Bing Juang, Shen-Fu Hsiao et al. proposed Parallel CORDIC algorithm in 2004.Its thought is root According to the n positions binary representation of input angle, high m is classified as and low n-m, Binary-to- first is implemented to a high position Bipolar (BBR) is changed, and once produces the rotation sequence of a high position, and the angular error for then producing a high position is superimposed upon low level, so Changed afterwards to low level using BBR, once produce the rotation sequence of low level.In order to maintain unified improvement factor K, for once High position rotation 2^-iRadian, compensate that (the method is referred to as microrotation using the micro- rotation for unifying direction several times Angle recoding), it is compensated in the error that ith rotates generation in low level without producing high-order essence Degree loss.Although this method realizes CORDIC parallel processing, the micro- rotation much compensated is added (in input angle Digit be 32 when, it is necessary to extra 8 micro- rotations), cause the delay in X/Y paths to increase.Publication No. is the A of CN 102073471, The patent of invention " a kind of processor Cordic interative computations method and circuit " come into force has been authorized to open up three interative computations parallel Open, it is therefore an objective to think once to calculate three iteration, but the direction of this three iteration can not be determined before iteration is started, it is desired nonetheless to Determined in path is calculated how much do not make computing acceleration actually.Major part cordic algorithm is all to use flowing water at present Line computation or iteration structure.The value for obtaining a calculating generally requires many clock cycle, is highly detrimental to cordic algorithm Quick realization on FPGA.

The content of the invention

In view of the shortcomings of the prior art, propose that a kind of new low z paths judge, low latency is high-precision parallel CORDIC computational methods and device.

Rotation radian is available to be represented in binary asWherein i ∈ { 0,1 }, θ≤π/4.

The radian of input is split into three parts：A high position, middle position, low level：

M is satisfaction 2^-m-tan^-12^-m<2^-NMinimum value, m=[(N-log can be drawn₂3)/3].As i >=m, tan2^-i=2^-i。l For N/2, now mould improvement factor

According to prior art to θ^HBBR conversions are carried out, be can obtain：

According to (9) formula, because b_i-1∈ { 0,1 }, so r_i∈ { -1,1 }, thus can be straight according to high m-1 that input radian Connect and obtain 1 to m direction of rotation and r₁,r₂,…,r_mValue.

Now, ideally the angle of 1 Dao m times rotation is all 2^-i, i=1,2 ..., m.But now I=1,2 ..., tan 2 is difficult on m-1. hardware^-i。

To radian 2^-iCarrying out quantification treatment can be expressed from the next：

The angle value rotated for ith.As shown in Fig. 2 solid line, which is ith, rotates mesh Angle is marked, dotted line is the actual anglec of rotation, if rotation exceedes angle on target, angular error-ε_iIf, rotation miss the mark angle Degree, angular error is ε_i。

Now ith transformation matrix is changed into：

σ_jValue meet formula 12 on the premise of, take

1 to m-1 times rotation is substituted by formula (11), and θ can be obtained by formula (10)^HThe remaining angle of generation is：

Convolution (12) and formula (16) are understoodSo σ_jMeet formula (12) just The error propagation that a high position can be produced is to low level without losing precision.Can be by simply shifting on transformation matrix (11) hardware Be added realization.

Further by first m times conversion, remaining radian is:

To θ^M′Continue to use and θ^HSame conversion, is obtained:

Wherein r_m=1-2b '_m-1,r_i=(2b_i-1- 1), i=m+1, m+2 ..., l.

R can be obtained according to formula (18)_m,r_m+1,…,r_lThe value of secondary rotation, because i meets tan in this interval^-12^-i=2^-i, it is as follows that transition matrix is changed into classical matrix：

By l-m rotation of second stage, remaining radian is:

At i ∈ [l, N], due to mould improvement factor cos θ_i=1, it can be deduced that if follow-up some rotations are not carried out, The unified Mod correction factor can be kept, so line translation need not be entered, now can be according to θ^L" simple single direction rotation is carried out, for b_i"=0 without rotation, b_i"=1 carry out single direction rotation.Directly according to binary system value in place, directly determine remaining The order of rotation.Further, i is worked as>N/2, arbitrarily rotation adjacent twice is expressed from the next：

2 in formula (21)^-(m+n)<2^-N, beyond minimum precision, so formula (21) can be reduced to：

N/2 can be obtained by formula (22) to represent into following matrix to n times rotation is compressible：

In formula (23), if remaining radian is negative, s is -1, is otherwise 1.Formula (23) can be by twice without sign multiplication, twice Addition is realized.

The prediction of above-mentioned middle direction of rotation can be divided into two steps, and the first step is according to θ^HDraw r₁,r₂,…,r_mValue, second step According to θ^M′+θ^L′, to θ^M′And θ^L′Using different Forecasting Methodologies, r can be drawn simultaneously_m,r_m+1,…,r_l,r_l,…,r_NValue.So As long as overall Parallel CORDIC computing two is clapped.

According to above-mentioned derivation, the overall correction factor that puppet rotation is produced can be calculated in advance, will not be because of input angle It is different and different, the advantage of traditional cordic algorithm is remained, and micro- rotation that need not be additionally is compensated.

Brief description of the drawings：

Fig. 1 is cordic algorithm rotation schematic diagram.

Fig. 2 is that ith rotation error produces schematic diagram.

Fig. 3 is the spin matrix structural representation of formula (11).

Fig. 4 is 32 system the first spin matrix hardware architecture diagrams of the present embodiment.

Fig. 5 is the present embodiment general structure schematic diagram.

Embodiment：

The present invention specific implementation step be：

(1) to the XIN of input, YIN, ZIN enters line translation, make input angle angle value be adjusted to (0, π/it is 4) interval in.To adjustment Radian afterwards splits into three parts：A high position, middle position, low level：

(2) first m-1 times of θ is predicted and draws r₁,r₂,…,r_mValue, to preceding m-1 rotate angle on target measure Change, produce new spin matrix.Preceding m rotation error is accumulated to low level, new θ is obtained^M′+θ^L′, according to θ^M′+θ^L′, to θ^M′ And θ^L′Using different Forecasting Methodologies, r can be drawn simultaneously_m,r_m+1,…,r_l,r_l,…,r_NValue.

(3) classical matrix is used to m to l rotations, l is directly tried to achieve to n times rotation using multiplication.Finally carry out quadrant Recover output and calculate obtained trigonometric function value.

The present embodiment is illustrated by taking the system of 32 as an example, for 32 systems, and maximum radian is 2 π, so Using 3 binary representation integers and 29 binary representation decimals.For input radian not (0, π/it is 4) interval in, root According to existing interval folding, it is easy to the radian of input is adjusted to (0, π/4) in, do not repeating here.

Radian after adjustment is represented byWherein b_i∈ { 0,1 }, θ≤π/4.

The radian of input is split into three parts：A high position, middle position, low level：

For θ^HBBR conversions are carried out, are obtained：

R can once be obtained according to formula (27)₁,r₂,…,r₁₀Value.

Micro- anglec of rotation and matrix of first 10 times be：

θ₁=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)

The traditional spin matrix of 4 to 10 uses

It is by each above-mentioned micro- error produced that rotates：

ε₁=| 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)

By θ^HThe maximum radian of the error of generation is

The error amount of each micro- rotation is tabulated in advance and preserved.The error amount carry save adder (CSA) of each step rotation Superposition produces new

To middle position θ^M′Carry out mathematics conversion available

r₁₀=1-2b '₉,r_i=(2b_i-1- 1), i=11,12 ..., 15 (51)

R can be directly obtained₁₀,r₁₁,r₁₂,r₁₃,r₁₄,r₁₅Value.Traditional spin matrix is used to i ∈ [10,15].

15 to 29 spin matrixs can be expressed from the next afterwards：

According to formulaTrue form represent, if negative, s=-1, conversely, s=1.r_iWith position power 2^-iPlace The value of position is equal, r_i∈{0,1}.15 to 29 rotations are directly by multiplication twice, and two sub-additions are realized.Above-mentioned middle direction of rotation Prediction can be divided into two steps, and the first step is according to θ^HDraw r₁,r₂,…,r₉Value, second step is according to θ^M′+θ^L′, to θ^M′And θ^L′Using not Same Forecasting Methodology, r can be drawn simultaneously₁₀,…,r₁₅,r₁₅,…,r₂₉Value.

5 illustrate below in conjunction with the accompanying drawings：

1 is pre-processing module in accompanying drawing 5, and using interval folding, to the XIN of input, YIN, ZIN enters line translation, makes Input angle angle value be adjusted to (0, π/it is 4) interval in, and export quadrant information.

3 be first angle prediction module in accompanying drawing 5, and three parts are split into the radian after adjustment：A high position, middle position is low Position：

To θ^HMathematics conversion is carried out, r is once produced after recodification₁,r₂,…,r₁₀.The determination method of each value is as follows：Due to Radian after adjustment is just, so r always₁=1, r_i=1-2b_i-1.Determine r₁,r₂,…,r₁₀After value, with reference to α₁,α₂,…,α₉'s Value, the angular error produced every time after rotation and θ^M+θ^LAddition forms θ^M′+θ^L′。

α above-mentioned₁,α₂,…,α₉Determined by each postrotational error amount, if the angle actually turned over is more than target Angle, then as previous α_i=-1, if the angle actually turned over is less than angle on target, as previous α_i=1.

2 be that the first rotary module includes 10 rotations in accompanying drawing 5.The angle of first order rotation is θ₁=tan^-1(2^-1+2^-5+ 2^-6), transition matrix is：The angle of second level rotation is θ₂= tan^-1(2^-2+2^-8+2^-10), transition matrix isThe third level is revolved The angle turned is θ₃=tan^-1(2^-3+2^-11+2^-13), transition matrix isLevel Four is followed successively by θ to ten grades of anglecs of rotation below₄= tan^-12^-4, θ₅=tan^-12^-5, θ₆=tan^-12^-6, θ₇=tan^-12^-7, θ₈=tan^-12^-8, θ₉=tan^-12^-9, θ₁₀=tan^-12^-10 =2^-10, transition matrix is traditional matrix.r₁,r₂,…,r₁₀After first angle prediction module is obtained, start by the first order to Tenth grade calculates successively, and the delay of X/Y paths strictly keeps the same.It is optional according to prior art in the calculating of X/Y additions The 4-2 compression adder speed-up computations of CSA adders synthesis are selected, the X/Y values calculated input the second selecting module.

5 be second angle prediction module, the cumulative incoming second angle of remaining angle of first angle prediction module in accompanying drawing 5 Prediction module.The the 10th to the 14th of remaining angle converted with the first prediction module high-value identical equation, is produced r₁₀,r₁₁,r₁₂,r₁₃,r₁₄,r₁₅Value, while producing low N/2 remaining angle, θ^L", according to θ^L" value produce r₁₅,…,r₂₉ Value.r₁₀=1-2b '₉,r₁₁=2b '₁₀-1,r₁₂=2b '₁₁-1,r₁₃=2b '₁₂-1,r₁₄=2b '₁₃-1,r₁₅=2b '₁₄-1。 r₁₅,…,r₂₉With θ^L" true form correspondence, r₁₅=b "₁₅..., r₂₉=b "₂₉。

4 be that the second rotational calculation module includes 6 tradition rotation calculating and a multiplication is calculated in accompanying drawing 5.Second rotation Module receives the X/Y of the first rotary module output, receives the r of second angle prediction module output₉,r₁₀,…,r₁₅,r₁₅,…, r₂₉Value.r₉,r₁₀,…,r₁₅Direction for determining the 10th grade to the 15th grade tradition rotation.r₁₅,…,r₂₉Value determine 15 Multiplier.

6 be post-processing module in accompanying drawing 5, post-processing module to the X/Y of the second rotational calculation module according to exported in 1 as Limit information carries out the recovery of quadrant.Finally export X/Y value.

Claims (4)

1. a kind of coordinate transform parallel calculating method and device, it is characterised in that including step：

(1) input angle angle value, to the XIN of input, YIN, ZIN is pre-processed, and preserves quadrant information, and the system effective accuracy is N；

(2) by the radian value after processing according to binary form indicating value, three sections are divided into, rotation prediction first is carried out to first paragraph, really Determine m rotational direction value of first paragraph, it is then pre- according to remaining radian two sections after the radian residual error of first paragraph is added to Survey remaining direction of rotation；

(3) quadrant information preserved according to pre-processing module, carries out quadrant recovery, exports trigonometric function value.

2. a kind of coordinate transform parallel calculating method according to claim 1 and device, it is characterised in that the step (2) first paragraph rotation prediction and computational methods in are concretely comprised the following steps：

(2.1) the rotation index of first paragraph is 1 Dao m times, and described m is less than [(N-log₂3) maximum integer/3], according to arc It is 2 to spend position power^-1To 2^-mValue, directly obtain 1 to m rotational direction value；

(2.2) radian 2 is rotated for ith target^-i, described i belongs to 1 between m, to the radian 2^-iCarry out at quantization Reason,σ_j∈ { -1,0,1 }, α_i∈ { -1,1 }, σ_jValue meetOn the premise of, take

(2.3) ith spin matrix is transformed to by traditional spin matrix

<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) residual error by first m times rotation is added to low level radian, forms new radian value.

3. a kind of coordinate transform parallel calculating method according to claim 1 and device, it is characterised in that the step (2) second in, three sections of rotation predictions and computational methods are concretely comprised the following steps：

(3.1) computational methods of second segment are that the new radian value formed described in 2.4 directly predicts m according to claim 2 To N/2 rotational direction value, to m to N/2 rotations using traditional spin matrix；

(3.2) rotate the remaining radian error produced to second segment to be added to the 3rd section, produce new radian value, it is then directly pre- It is N/2 to n times rotational direction value to measure remaining radian index.

4. a kind of coordinate transform parallel calculating method according to claim 1 and device, it is characterised in that the step (2) the 3rd section of computational methods in are：

According to claim 3, the N/2 obtained in 3.2 to n times rotational direction value, by N/2 to n times spin matrix boil down to：

r_i∈ { 0,1 }, if remaining radian is negative, s is -1, is otherwise 1.