現実世界はアナログかデジタルか
新緑の息吹を感じる季節(図 1)、ワタシはふと思いました。「さて、現実世界はアナログなのかデジタルなのか」。アナログ・デバイセズという社名はその名のとおり、アナログである物理量を取り扱うという意味があります。
「電圧量や電流量はアナログ値」とはいいますが、新緑の息吹を感じるこのアナログな季節に、「電子の電荷量𝑒は
だよな。0.00016pA の電流量は電子
個が、1 秒間に導体中を移動することになるわけだ!」と思いました。1000 個、つまり離散値として数えることができるわけですから、「電流量はデジタル量といえるな!」に至りました。
また PN 接合を電流𝐼 [A] (電子)が流れるときに生じるショット雑音[1]の 1Hz 電流密度
も、電子 1 個が PN 接合を横切るときに発生する電流性ノイズの合算ですので、「電流量はデジタル量」という考えを補強するものといえるでしょう。
ところでアナログ・デバイセズでも電流検出回路のことを「クーロン・カウンタ」(カウンタ!)とか呼んだりしますので、これも驚きをもって受け入れています(笑)。
でも量子力学まで考慮すると、どうなるんでしょうね。それでも上記の考え方でふるまいの連続性がいったん途切れるのでしょう。詳しい方で素人にも分かりやすくご説明いただける方、ぜひご一緒にお茶のみでもいたしましょう!
CORDIC を分かりやすく解説してみたい
CORDIC(COordinate Rotation DIgital Computer)はデジタル演算処理で三角関数など(調べてみると、これがまた、かなり多彩な演算を実現できるようで、びっくりですが)を、簡易な論理回路構造で実現できるアルゴリズムです。関数電卓内のデジタル回路の論理構造(回路でなくソフトウェア演算のようですが…。詳しくは分からず…)にも活用されていますし、マイコンなどにもペリフェラルとして CORDIC 演算回路が内蔵されています。
もう少し言っておくと、CORDIC での演算は足し算、引き算、ビット・シフトを繰り返すことで実現できるものです。デジタル回路では実現しやすいアルゴリズムなのです。
私も CORDIC については、名前と基本的な考え方は知ってはいましたが、具体的なアルゴリズムやインプリメント(実装)方法は分かっていませんでした。これまでも適切な記事を探していましたが、数点ほど読んでみても、いまひとつ理解できませんでした。しかし参考文献[2]を見つけて「なるほど!」となりました。その CORDIC の理解が浅かったころ、とあるサイトに GNU の VHDL ソースがありました。それを使って実際の実験もおこなったこともありました(次の技術ノートでその実験結果も示します)。
今回と次回の技術ノートでは、ネットに多数解説記事は出ているものの、どうも理解不能というか腹落ち感の低い、この CORDIC を、我ながらの手法で解説してみようと思って書き始めております。いつもながら、書き始めた段階では、まだトンネルの出口には到達していない、手探り状態というところです。
CORDIC の計算手順(まずは答えの得かたのご紹介)
これまでいろいろな記事を読んでも良く分からなかった(というよりは読み始めてすぐに挫折していた)CORDIC は、実はとても簡単な計算手順なのでした。
ところで今回の技術ノートでは、文章と数式の説明をビジュアル的にも理解いただけるように、説明の部分付近に、その説明に該当する図を挿入するレイアウトにしてみました。いくばくかでもご理解が深まれば幸いです。
初期値の設定
目的の角度を𝜃とし、その SIN = 正弦、COS = 余弦を生成することを考えます。CORDIC 計算開始の初期値(𝑖 = 0)として
と設定します。詳しくは以降に示しますが、これは図 2 のイメージになります。𝑧(0)に目的の角度𝜃を設定するわけなのですね…(これがポイント)。図 2 では𝑧(0) = 60° と仮定しています。
以降で詳しく説明していきますが、十分な繰り返し演算が行われると、答えとして
が得られるというアルゴリズムです。答えとしては𝑧(𝑁)は不要なものです。
ここで𝐴𝑁はスケーリング係数と呼ばれる定数なのですが、これは以降で示します(繰り返し回数が十分おおきければ 1.64676⋯に収束するものです)。この𝑥(0)に1/𝐴𝑁を設定することなく、
𝑥(0) = 1
を設定して、最後に得られた答えを𝐴𝑁で割っても「数学的」には問題ありません。しかし実デジタル回路での計算回路規模を減らす目的からすれば、このやり方だと「本末転倒的」になってしまいます。
ステップ演算を繰り返す漸化式
つづいて𝑖がゼロからスタートする漸化式(漸化式とは、ひとつ前のステップの演算結果を用いてそのステップでの計算結果を得て、それを繰り返すもの)
を考えます。𝑥(𝑖), 𝑦(𝑖), 𝑧(𝑖)は単なる数値(デジタル値)です。 𝐷(𝑖), 𝑠(𝑖)は以降で示します。図 3 に、1 回目の計算を直交座標としてグラフィカルに表現したようすを、そして図 4 に 2 回目の計算のようすを、それぞれ図 2 の設定[𝑧(0) = 60° ]を例にして示します。
𝐷(𝑖)は CORDIC の回転方向(極性)で
𝑧(𝑖)がマイナスなら𝐷(𝑖) = −1、𝑧(𝑖)がプラスなら𝐷(𝑖) = +1とするだけです。図 3、図 4 においては𝐷(0) = +1,𝐷(1) = +1という設定になっています[𝑧(0)が 45°を超えている条件なので]。
≫ 𝑖は𝑖ビットの右シフト演算です。もっと数学的に書くと
とも表現できます。式(3)の𝑧(𝑖)の式中や、図 3 や図 4 に
がありますが、これは 1 ステップの演算における角度変化量です。これはあらためて以降で説明します。
𝑠(𝑖)のイメージを図 5 に示します。ここでは𝑠(0), 𝑠(1), 𝑠(2)を求めるようすをビジュアルに示しています。
デジタル回路の実装においては、𝑠(𝑖)は計算済みテーブルとして事前に用意しておけばいいだけです。テーブル𝑠(𝑖)の単位(たとえば弧度法なのかラジアンなのか、ビット表記や小数点位置なども含めて)は、𝑧(0)に設定する目的の角度𝜃の単位と同じ にしておく必要があります[𝑠(𝑖)は𝑧(𝑖)と足し引きされるため]。これは重要なポイントです。
繰り返し回数を十分に大きくすれば答えが得られる
ここで𝑖を十分大きく、つまり演算繰り返しを十分な回数𝑁にすると(𝑖 + 1 = 𝑁と表します)、
に収束します。この式を理解しようとすると少し難儀するかもしれませんので、ひとまずスルーして OK です(今回の後半と次回に詳しく説明します)。
実際の CORDIC の使い方として sin/cos を生成する目的であれば、最初の説明のように、それぞれの初期値を𝑥(0) = 1/𝐴𝑁, 𝑦(0) = 0とし、また𝑧(0) = 𝜃として目的の角度𝜃𝜃を設定します。そうすると上記の式(7)は、
となり、𝑁回の繰り返し演算でこれらが得られるというお話です。このようすを図 6 に示します。
繰り返し回数𝑁は、テーブル値𝑠(𝑖) = 0になるところ(デジタル演算として 1 ビット以下になるところ)を目安とします。𝑠(𝑖) は事前計算済みテーブルなので、𝑠(𝑖) = 0になるところは事前に既知、すなわち繰り返し回数𝑁は事前に既知となります。これが CODRIC の計算、うごきなのです。
なお𝑧(0)として対応可能な範囲は−𝜋/2 ≤ 𝑧(0) ≤ 𝜋/2になるので、他の象限部分で計算をさせたい場合は、初期化時におまじないが必要です。詳細は参考文献[2]の p. 3 左コラム中央あたりをご覧ください。
スケーリング係数𝑨𝑵を考える
さきに出てきた定数𝐴𝑁は図 7 に示すように、CORDIC の漸化式を逐次的にステップごとに計算していくうえで生じる、図 7 の緑色破線の長さが増加する係数𝐴(𝑖)の繰り返し積で、
𝑁が十分大きければ 1.64676⋯に収束します。ここで∏という記号は、要素ごとをすべて掛け算する「総乗」という意味です[3]。総和∑の掛け算版といえます。この𝐴𝑁はスケーリング係数と呼ばれます。これは次回詳しく説明します。
実際にエクセルで計算してみる
演算論理ブロックをお見せするまえに(これは次回にご説明します)、この計算がホントに正しく実現できるかエクセルを使ってやってみましょう。
図 8 はエクセルの表計算上で CORDIC を計算してみたようすです。この演算では目的の角度𝜃を 57°として、𝑧(0) = 57° とします(図中の 2 行目を黄色でハイライトしています)。またここまでの説明のように𝑥(0) = 1/𝐴𝑁, 𝑦(0) = 0とします。
なお𝐴𝑁 = 1.64676 ⋯なので、𝑥(0) = 1/1.64676 = 0.607253となります。
「計算済みテーブル𝑠(𝑖)」は、ここではエクセルで
をそのまま計算させています。
𝑁 = 16で 16 回の漸化式演算(ステップごとの演算)を行った結果が表中の一番下、𝑖 = 16のところです(黄色でハイライトしています)。表中の結果𝑥(16)、𝑦(16)が、cos(57°)、 sin(57°) に近くなります。それぞれのセルの欄外下にエクセルで計算した cos(57°)、sin(57°)の結果を示してありますが、かなり近くなっていることが分かります。
これで CORDIC の計算が、説明した方法で正しくできることを確認できました。
CORDIC の理論解説をしてみよう
ネットには各種の理論解説があるが
ネットで CORDIC の理論解説の記事を検索してみると、たくさん出てきます。多くが三角関数の回転行列(これは以降でも使っていきたいと思っていますが)や直角三角形をつなぎ合わせた説明になっています。読む気がない私が悪いのでしょうが(笑)、「簡単にわかる」とか「誰でもわかる」とか書いてある記事のわりには、どれもすっと腑に落ちてきません。
そこでこの WEB ラボでは、独自の切り口(?)と表現方法で、この理論解説をしてみたいと思います。
糸口「その 1」は𝒛(𝒊)を残留角度誤差量と考えること
ものごとを理解していくには、何かを「糸口」として見つけることがポイントです。この CORDIC の理解も同じでしょう。そういう私もこの技術ノートを書き始めながら、まず、そこを探っていきました。
一番重要と思われる糸口が、𝑧(𝑖)は、計算していくうえでの「角度の残留角度誤差量」だと考えることです。𝑧(𝑖)は
であり、引き算される𝑠(𝑖)が 1 演算ステップごとで小さくなり、 𝑧(𝑖)から引き算されていきます。図 8 の最初の 3 ステップの演算を例にして、𝑧(𝑖), 𝑠(𝑖)の変化のようすを図示したものを図 9 に示します。目的とする角度は 57°つまり、𝑧(𝑖) = 57° です。図 8 においても、演算ステップを追うごとで残留角度誤差量である 𝑧(𝑖)が低減していくことが分かります。
𝑧(𝑖)自体の極性がプラスかマイナスかで、𝑠(𝑖)が𝑧(𝑖)から引かれるか足されるか(極性)が𝐷(𝑖)として決まります。この「引かれるか/足されるか」という表現も「残留角度誤差量を低下させていく操作」とみることもできます。
また演算繰り返しを十分な回数𝑁にすると
に収束すると説明しました。ここからもさらに「𝑧(𝑖)は計算していくうえでの残留角度誤差量」ということが理解できます。
糸口「その 2」変化角度𝒔(𝒊)は直角三角形として決まり、その底辺が 1 演算ごとで 2 倍づつになっていく
つづいて減算量𝑠(𝑖)です。𝑠(𝑖)は 1 演算ステップごとで小さくなっていきます。𝑠(𝑖)自体の最初の 3 ステップでの値が、どのようにグラフィカルに表されるかを図 10 に示します。𝑠(𝑖)は
ですから、
● 1 回目の演算ステップ、𝑖 + 1 = 1(𝑖 = 0)では、カッコ内 は 1/1 となり、𝑠(0)は図中の直角二等辺三角形(赤紫色/マゼンタ色)内に示している角の角度(45°)になる
● 2 回目の演算ステップ、𝑖 + 1 = 2(𝑖 = 1)では、カッコ内 は 1/2 となり、𝑠(1)は図中の底辺が 2 倍となる直角三角形内に示している角の角度(26.565°)になる
● 3 回目の演算ステップ、𝑖 + 1 = 3(𝑖 = 2)では、カッコ内 は 1/4 となり、𝑠(2)は図中の底辺が 4 倍となる直角三角形内に示している角の角度(14.036°)になる
このように𝑠(𝑖)は、1 演算ごとで底辺が 2 倍づつになっていく直角三角形の鋭角の角度なのです。
なぜこのように𝑠(𝑖)を「2 倍づつになっていく直角三角形」の鋭角の角度とするのでしょうか。それこそ 45°、(45/2)°、 (45/4)°…というようにしていけばよいのでは?と頭をよぎるかと思います。
このようにする理由は、CORDIC の 1 演算ステップごとで、𝑖ビットのシフトにより演算を実現する(実際は1/2𝑖にする)ために必要な角度設定だからです。
詳しくは次回に示していきますが、このことは「ホント CORDIC、よくできてるな!」と思わせるポイントでもありました。次回の説明をご期待ください!
糸口「その 3」なんとこの計算は逐次比較(SAR)型 ADC の動作と同じではないか!
糸口「その 1」のように、残留角度誤差量𝑧(𝑖)と減算量となるテーブル値𝑠(𝑖)は、
であり、減算量𝑠(𝑖)は 1 演算ステップごとで小さく、そしてその減算極性は𝑧(𝑖)自体の極性がプラスかマイナスかで決まります。
これはまるで図 11 のような天秤と分銅による重量計測(重さが 1/2 の分銅を載せたり下したりする)と同じということに気がつくかと思います。さらにこれは図 12 のような逐次比較(Successive Approximation Register; SAR)型 ADCの動作と同じという気づきに至ることになると思います。
このように考えていけば、「CORDIC もそんなにややこしいものでもないな」と、気持ちの敷居を下げられるのではないでしょうか。
糸口「その 4」は別に見出しを立てて
この糸口「その 4」の説明は少し長くなりますので、いったん別に見出しを立てて別の節として取り扱ってみましょう。
糸口「その 4」高校のときに習った三角関数の加法定理も糸口
高校のときに「三角関数の加法定理 [5]」というのを学習した記憶があるのではないでしょうか。
私が高校生のときに習った数学の先生は、これを「サイン・コス・コス・サイン、コスコス・サイン・サインと覚えなさい」とおっしゃっていました。
これがなんと!このまま CORDIC の繰り返し演算の 1 回分の操作、また𝑁回の繰り返し演算(つまり𝑖 + 1 = 𝑁)で得られる計算結果に相当するのです!
なおこの「三角関数の加法定理」については(一般的な高校数学であることから)本技術ノートではとくに説明や詳細な解説は行いません。証明なども含めて、ネット上に多数の解説がありますので、そちらをご覧いただければと思います。[5]にも解説が載っています。
𝑵回繰り返し演算で得られる計算結果と加法定理との関係性を考える
ここではまず CORDIC での𝑁回の繰り返し演算(つまり𝑖 + 1 = 𝑁)で得られる計算結果𝑥(𝑁),𝑦(𝑁)について、加法定理との関係を説明し、次回の技術ノートで CORDIC の演算ステップ 1 回分について加法定理との関係を説明していきます。
加法定理を行列計算に変形する
さきの式(10)の sin と cos(上下の式)を逆転してみます。
これを行列的に表記しなおしてみると(図 13 がその説明)
これを「回転行列」というようですね。なおここでの行列計算とは
に相当し、図 13 もこの式と同じ変換になります。
加法定理を行列に変形したものをベクトル演算と考える
ここでベクトル
として、これまでの式の𝛽に関連する部分を、「長さ(ノルム) が 1 で、角度𝛽のベクトル𝑽(𝛽)」の X 成分(cos 𝛽)と Y 成分 (sin 𝛽)として、図 14 のように考えます。
また
として回転行列を、ひとつの記号𝑹(𝛼)で表します。ここでは回転行列の回転量パラメータを回転角度𝛼とします。
回転行列でベクトルの回転ができる
この行列的表記により、図 15 のように
と表現すれば、ベクトル𝑽(𝛽)(青の破線矢)をベクトル𝑽(𝛽 + 𝛼) (赤の破線矢)に変換する、ベクトルの角度を回転できる演算になっていることが分かります。
要素ごとに分解して表記しても、回転行列で要素ごとの角度変化ができている
同じことを繰り返しやってみるわけですが…。またこれをもとに戻して、X 成分(cos 𝛽)と Y 成分(sin 𝛽)として要素ごとに分解してみます(数学的にはこれも「ベクトル」と表現されますが…)。
これは長さ(ノルム)が 1 で角度𝛽のベクトルの、X 成分である cos 𝛽と、Y 成分である sin 𝛽を、回転行列𝑹(𝛼)により、変化角度 𝛼で、「角度𝛽 + 𝛼の cos と sin」までその角度を変化できるということです。
この説明をすっ飛ばして回転行列だけとか、複素数とかで CORDIC を解説している記事も多い(読者に理解してもらいたい意思があるのか??)
英語の解説ページが多いのですが、上記にした説明をすっ飛ばして、いきなり回転行列
の式のみを示していて、「CORDIC はこうですよー」とか「CORDIC はこうなりますよー」とか書いてあったり、複素数極座標表記
で CORDIC の回転が表現できますとか、式の関係性の導出なしに書いてあったりする解説があります。
結局、これでは何がなんだか、わけがわかりません(涙)。「読者に理解してもらいたい意思があるのか?」とか思います…。それとも本質が分かってないのかな…、とかも思います…。
ともあれこのように cos と sin(ベクトル)の角度を、高校で習った加法定理での「回転行列」を用いることで変化させることができるのです。これで CORDIC の演算を説明できるのです。
回転行列が CORDIC 演算そのもの
この「回転行列」が CORDIC の演算そのものなのです。 CORDIC で得られる計算結果の一般式は、式(7)のように
𝐴𝑁は図 7 に示したように、CORDIC の漸化式を逐次的にステップごとに計算していくうえで生じる、ベクトルの長さが増加する増加係数の繰り返し積で、𝑁が十分大きければ 1.64676⋯に収束します。
CORDIC 漸化式(一般式)を行列に変形する
この𝑦(𝑁)の式において、項の順序を変更し、
つづいて各項の表記順序を逆にして
そして行列的に
を使って表してみます。そうすると
これは式(16)と全く同じかたちになり、
とすれば、図 16 のようにベクトルの成分𝑥(0), 𝑦(0)が
として、𝑹[𝑧(0)]で図 15 と全く同じように、𝑥(𝑁),𝑦(𝑁)に回転変換されることが分かります。「ベクトル𝑽」の表記で、
とすれば、
とも表記できます。ここでも𝑹[𝑧(0)]で、𝑽(0) が𝑽(𝑁)に回転変換されることが分かります。
CORDIC 漸化式行列(一般式)に初期値を設定する
さらに式(1, 2)で示したように、CORDIC で sin/cos を生成する目的であれば初期値と答えは、
なので、式(17)はもっと視覚的に分かりやすい形で
と表されます。
初期値を角度ゼロのベクトルと考える
さらに図 17 のように、初期値を角度 0°のベクトルだとして
として𝑥(0)と𝑦(0)を設定すれば、式(17)や式(20)は
となります。
角度ゼロから回転行列で角度𝜽のベクトルの成分(CORDIC の答え)が得られる
この CORDIC の式は、0°から角度𝜃への変換を、回転行列
をもちいて
として、回転角度量を𝜃にして計算する演算とまったく同じということが分かります。
まとめ
書き始めた段階では、まだトンネルの出口には到達していない、手探り状態というスタートでしたが、探求していくことで満足できるところまで理解を深めることができました。
CORDIC とは初期値𝑧(0)を目的の角度、𝑧(𝑖)を残留角度誤差量として、加法定理を基礎とした計算アルゴリズムで、ベクトルの sin/cos 値を求めていくものだと分かりました。またそれは SAR ADC の動作とかなり近いしくみだとご説明しました。
次回はこの理解から、CORDIC の 1 ステップ演算がどのように構成されるかなどをより深く考えていき、実際に CORDIC アルゴリズムを FPGA に DDS(Direct Digital Synthesizer)として実装し、DACIC を使って正弦波を出力してみた例をご紹介したいと思います。
しかし「よく出来てるなぁー」としきりに感心する、CORDIC アルゴリズムでありました!