CORDIC を使用した平方根の計算

この例では、MATLAB® で CORDIC カーネル アルゴリズムを使用して平方根を計算する方法を説明します。CORDIC ベースのアルゴリズムは、モーター制御、ナビゲーション、信号処理、ワイヤレス通信などのさまざまな組み込みアプリケーションにとって重要です。

CORDIC は、COordinate Rotation DIgital Computer の略語です。ギブンス回転に基づく CORDIC アルゴリズムは、反復的なシフト加算演算だけが必要であるため、ハードウェア効率が最も高いアルゴリズムのうちの 1 つです ([1] および [2])。CORDIC アルゴリズムは、明示的な乗数を必要とせず、正弦関数、余弦関数、逆正弦関数、逆余弦関数、逆正接関数、ベクトル振幅関数、除算関数、平方根関数、双曲線関数、対数関数など、さまざまな関数の計算に適しています。

固定小数点 CORDIC アルゴリズムには次の演算が必要です。

• 反復あたり 1 回のテーブル ルックアップ

• 反復あたり 2 回のシフト

• 反復あたり 3 回の加算

平方根など、CORDIC ベースの双曲線アルゴリズムでは、結果の収束を得るために特定の反復 ($$i=4,13,40,121,\dots ,k,3k+1,\dots$$) が繰り返されることに注意してください ([2])。

双曲線計算モードを使用した CORDIC カーネル アルゴリズム

CORDIC 計算モード アルゴリズムを使用して、双曲線三角、平方根、対数、指数などの双曲線関数を計算できます。

双曲線ベクトル モードでの CORDIC 方程式

双曲線ベクトル モードは平方根を計算するために使用されます。

ベクトル モードでの CORDIC 方程式は以下のとおりです。

$$\begin{array}{lcl}{x}_{i+1}& =& x_i+y_i\cdot d_i\cdot 2^{-i}\\ y_{i+1}& =& y_i+x_i\cdot d_i\cdot 2^{-i}\\ z_{i+1}& =& z_i-d_i\cdot \text{atanh}(2^{-i})\end{array}$$

ここで、

$$i=\{1,2,3,4,4,5,\dots \}$$ は、$$3k+1$$ 反復ごとにステップが繰り返されるシーケンスです。

$y_i<0$ の場合は $d_i=+1$、それ以外の場合は $-1$ です。

このモードにより、$N$ が $+\infty$ に近づくときに以下の結果が得られます。

$$\begin{array}{lcl}{x}_N& \approx & A_N\sqrt{x_0^2-y_0^2}\\ y_N& \approx & 0\\ z_N& \approx & z_0+\text{atanh}(y_0/x_0)\end{array}$$

ここで、

$A_N={\prod }_{i=\{1,2,3,4,4,5,\dots \}}\sqrt{1-2^{-2i}}$.

積内のシーケンスは、$$3k+1$$ ステップごとに繰り返されることに注意してください ([2])。

一般的には、$N$ には大きな定数値が選ばれます。そのため、$A_N$ は事前に計算できる場合があります。

また、平方根には $x_N$ の結果のみが使用されることに注意してください。

MATLAB への CORDIC 双曲線ベクトル モード アルゴリズムの実装

CORDIC 双曲線ベクトル モード カーネル アルゴリズムを MATLAB コードに実装する例は、スカラー x、y および z の場合、次のようになります。固定小数点および浮動小数点のデータ型の両方にこの同じコードを使用できます。

CORDIC 双曲線ベクトル モード カーネル

k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:N
    xtmp = bitsra(x, idx); % multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % multiply by 2^(-idx)
    if y < 0
        x(:) = x + ytmp;
        y(:) = y + xtmp;
        z(:) = z - atanhLookupTable(idx);
    else
        x(:) = x - ytmp;
        y(:) = y - xtmp;
        z(:) = z + atanhLookupTable(idx);
    end
    if idx==k
        xtmp = bitsra(x, idx); % multiply by 2^(-idx)
        ytmp = bitsra(y, idx); % multiply by 2^(-idx)
        if y < 0
            x(:) = x + ytmp;
            y(:) = y + xtmp;
            z(:) = z - atanhLookupTable(idx);
        else
            x(:) = x - ytmp;
            y(:) = y - xtmp;
            z(:) = z + atanhLookupTable(idx);
        end
        k = 3*k + 1;
    end
end % idx loop

逆双曲線正接ルックアップ テーブル atanhLookupTable は次のように定義されます。これはコード生成時に事前に計算され、定数になります。

atanhLookupTable = cast(atanh(2.^-(1:N)),'like',z);

CORDIC 双曲線ベクトル モード カーネルを使用した平方根の計算

適切な初期値を選択することにより、CORDIC カーネル双曲線ベクトル モード アルゴリズムで平方根を計算できます。

最初に、次の初期化手順を実行します。

• $x_0$ を $v+0.25$ に設定します。

• $y_0$ を $v-0.25$ に設定します。

$N$ 回の反復の後、これらの初期値は $N$ が $+\infty$ に近づくときの以下の出力を導きます。

$$x_N\approx A_N\sqrt{(v+0.25{)}^2-(v-0.25{)}^2}$$

これは、以下のようにさらに簡素化できる場合があります。

$$x_N\approx A_N\sqrt{v}$$

ここで、$A_N$ は上記で定義された CORDIC ゲインです。

なお、平方根については、$z$ および atanhLookupTable は結果に影響しません。したがって、$z$ および atanhLookupTable は使用されません。

MATLAB への CORDIC 平方根カーネルの実装

CORDIC 平方根カーネル アルゴリズムを MATLAB コードに実装する例は、スカラー x および y の場合、次のようになります。固定小数点および浮動小数点のデータ型の両方にこの同じコードを使用できます。

CORDIC 平方根カーネル

k = 4; % Used for the repeated (3*k + 1) iteration steps
for idx = 1:N
    xtmp = bitsra(x, idx); % multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % multiply by 2^(-idx)
    if y < 0
        x(:) = x + ytmp;
        y(:) = y + xtmp;
    else
        x(:) = x - ytmp;
        y(:) = y - xtmp;
    end
    if idx==k
         xtmp = bitsra(x, idx); % multiply by 2^(-idx)
         ytmp = bitsra(y, idx); % multiply by 2^(-idx)
         if y < 0
             x(:) = x + ytmp;
             y(:) = y + xtmp;
         else
             x(:) = x - ytmp;
             y(:) = y - xtmp;
         end
         k = 3*k + 1;
      end
 end % idx loop

このコードは、z と atanhLookupTable が使用されていない点を除けば、上記の CORDIC 双曲線ベクトル モード カーネルの実装と同じです。これにより、1 つのテーブル ルックアップと 1 つの加算を反復ごとに削減できます。

例

10 回の CORDIC カーネル反復を使用し、関数 cordicsqrt を使用して v_fix の近似平方根を計算します。

step  = 2^-7;
v_fix = fi(0.5:step:(2-step), 1, 20); % fixed-point inputs in range [.5, 2)
niter = 10; % number of CORDIC iterations
x_sqr = cordicsqrt(v_fix, niter);

比較のために CORDIC 出力の実際値 (RWV) を取得し、MATLAB 参照値と CORDIC 平方根値との誤差をプロットします。

x_cdc = double(x_sqr); % CORDIC results (scaled by An_hp)
v_ref = double(v_fix); % Reference floating-point input values
x_ref = sqrt(v_ref);   % MATLAB reference floating-point results
figure;
subplot(211);
plot(v_ref, x_cdc, 'r.', v_ref, x_ref, 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (In-Range) and MATLAB Reference Results');
subplot(212);
absErr = abs(x_ref - x_cdc);
plot(v_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

アルゴリズム入力範囲制限の回避

多くの平方根アルゴリズムは、入力値 $v$ を [0.5, 2) の範囲内に正規化します。通常、この前処理は固定語長正規化を使用して行われ、小さな入力値範囲だけでなく大きな入力値範囲をサポートできるようになります。

CORDIC ベースの平方根アルゴリズムの実装は、この範囲外の入力には特に敏感です。関数 cordicsqrt は、以下の数学的関係に基づく正規化手法でアルゴリズムの範囲制限に対処します。

ある $0.5<=u<2$ およびある偶数の整数 $n$ について $v=u\cdot 2^n$。

したがって、次のようになります。

$$\sqrt{v}=\sqrt{u}\cdot 2^{n/2}$$

関数 cordicsqrt 内で、前述のとおり、$u$ および $n$ の値は入力 $v$ の正規化中に検索されます。$n$ は、入力 $v$ の 2 進数表現の最上位ビット (MSB) の先頭のゼロの数です。これらの値は一連のビット単位の論理およびシフトによって検索されます。$n$ は偶数でなければならないので、MSB の先行するゼロの数が奇数の場合、$n$ を偶数にするために 1 つのビット シフトが行われることに注意してください。これらのシフト後の結果の値は $0.5<=u<2$ です。

$u$ は CORDIC ベースの平方根カーネルの入力になります。ここで、$\sqrt{u}$ の近似が計算されます。結果はその後、正しい出力範囲に戻されるよう、$2^{n/2}$ によってスケーリングされます。これは、$n/2$ ビットによる簡単なビット シフトによって達成されます。シフトの方向 (左または右) は $n$ の符号に依存します。

例

CORDIC を使用して、小さな非負の範囲をもつ 10 ビット固定小数点入力データの平方根を計算します。CORDIC ベースのアルゴリズムの結果を同じ入力範囲に対する浮動小数点 MATLAB 参照の結果と比較します。

step     = 2^-8;
u_ref    = 0:step:(0.5-step); % Input array (small range of values)
u_in_arb = fi(u_ref,0,10);    % 10-bit unsigned fixed-point input data values
u_len    = numel(u_ref);
sqrt_ref = sqrt(double(u_in_arb)); % MATLAB sqrt reference results
niter    = 10;
results  = zeros(u_len, 2);
results(:,2) = sqrt_ref(:);

同等の実際値 (RWV) の結果を計算します。CORDIC の RWV と MATLAB 参照の結果をプロットします。

x_out = cordicsqrt(u_in_arb, niter);
results(:,1) = double(x_out);
figure;
subplot(211);
plot(u_ref, results(:,1), 'r.', u_ref, results(:,2), 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (Small Input Range) and MATLAB Reference Results');
axis([0 0.5 0 0.75]);
subplot(212);
absErr = abs(results(:,2) - results(:,1));
plot(u_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

例

CORDIC を使用して、大きい正の範囲をもつ 16 ビット固定小数点入力データの平方根を計算します。CORDIC ベースのアルゴリズムの結果を同じ入力範囲に対する浮動小数点 MATLAB 参照の結果と比較します。

u_ref    = 0:5:2500;       % Input array (larger range of values)
u_in_arb = fi(u_ref,0,16); % 16-bit unsigned fixed-point input data values
u_len    = numel(u_ref);
sqrt_ref = sqrt(double(u_in_arb)); % MATLAB sqrt reference results
niter    = 16;
results  = zeros(u_len, 2);
results(:,2) = sqrt_ref(:);

同等の実際値 (RWV) の結果を計算します。CORDIC の RWV と MATLAB 参照の結果をプロットします。

x_out = cordicsqrt(u_in_arb, niter);
results(:,1) = double(x_out);
figure;
subplot(211);
plot(u_ref, results(:,1), 'r.', u_ref, results(:,2), 'b-');
legend('CORDIC', 'Reference', 'Location', 'SouthEast');
title('CORDIC Square Root (Large Input Range) and MATLAB Reference Results');
axis([0 2500 0 55]);
subplot(212);
absErr = abs(results(:,2) - results(:,1));
plot(u_ref, absErr);
title('Absolute Error (vs. MATLAB SQRT Reference Results)');

参考文献

1. Jack E. Volder, "The CORDIC Trigonometric Computing Technique," IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp. 330-334.

2. J.S. Walther, "A Unified Algorithm for Elementary Functions," Conference Proceedings, Spring Joint Computer Conference, May 1971, pp. 379-385.

参考

CORDIC Square Root HDL Optimized