人間の認知はアナログかデジタルか

前回同様、新緑の息吹を感じる季節、ワタシはふと思いました。 「さて、人間の認知はアナログ動作なのかデジタル動作なのか」。 「おお！脳内伝達物質は『分子』だな！」ここもデジタルなのでした…（笑）。

そういえば数値演算回路の本を買っていたぞ

そんな話題はよしとして（汗）、前回と今回は CORDIC （COordinate Rotation DIgital Computer）について解説しています。デジタルの加減算回路と乗算回路は規模が大きく違います（とくに 1 クロックで処理する場合）。タイトルのように CORDIC は足す・引くで計算できてしまいますので、コンパクトな論理回路を実現できます。

さて、今回の出だしの写真は何にしようかと思い巡らしていたところ、「そういえば」と本棚の奥底（？）にあった、図 1の本 [1]をご紹介にしようと思いました。

調べてみるとすでに絶版のようです。結局ほぼ（「ほとんど」までもいかず…）読んでいないまま、約 15 年のときが経過したようです。

ぱらぱらっとめくってみると Verilog で（Verilog は読めず…） RTL 記述があり、今更ながら「読んでカリカリなロジックを作ってみたい！」とか思うものでした。とかいいつつ、読める日はくるのだろうか（汗）。その前に何を作るか？（汗）

さて、今回の技術ノートは

前回は CORDIC で得られる計算結果（最終的に収束した結果） について、加法定理と回転行列との関係、というかアルゴリズムの成り立ちをご説明しました。

今回は CORDIC の漸化式での 1 演算ステップでは、どのような考え方になるかを考えていきましょう。

また今回は実際に DAC を用いて、CORDIC で作った NCO （Numeric Controlled Oscillator）で DDS（Direct Digital Synthesizer）信号を発生させたようすもご紹介します。

CORDIC 漸化式のおさらい

あらためて前回示した、 CORDIC の漸化式を示します。 CORDIC はこの式の計算を複数ステップ実行し、答えを得ます。 𝑖はゼロからスタートするインデックスです。

𝒙(𝟎), 𝒚(𝟎), 𝒛(𝟎) に初期値を入れて、演算ブロックに𝑵ステップ通して、𝒙(𝑵), 𝒚(𝑵)に答えを得る

ここで𝑥(𝑖), 𝑦(𝑖), 𝑧(𝑖)は単なる数値（デジタル値）です。sin/cos を得る CORDIC 演算では、計算開始の初期値（𝑖 = 0）として

と設定します。𝜃は得たい答えとなるsin/cosの目的の角度です。 𝐷(𝑖)は CORDIC の回転方向（極性）で

𝑧(𝑖)がマイナスなら𝐷(𝑖) = −1、𝑧(𝑖)がプラスなら𝐷(𝑖) = +1とするだけのものです。

角度𝑧(𝑖)は、1 ステップずつ漸化式で演算していくうえでの残留角度誤差量に相当します。繰り返し回数𝑁は、事前に設定しておくテーブル値𝑠(𝑖)が、𝑠(𝑖) = 0になるところ（デジタル演算として 1 ビット以下になるところ）を目安とします。

この演算が𝑁回繰り返されると（𝑖 + 1 = 𝑁で）

となり、𝑥(𝑁), 𝑦(𝑁)に目的の角度𝜃における sin/cos の値が得られることになります。

式(1)の漸化式（1 ステップの演算）をブロック図にしたものを 図 2 に示します。この図は、𝑖 + 1回目のステップの結果が得られるブロックになります。

前回の図 17 に示したように（今回も図 3 として再掲しました）、𝑥(0)と𝑦(0)を、角度 0°のベクトルの X 軸 = cos 成分、Y 軸 = sin 成分

だと理解しておくと、以降の説明の腹落ち感があると思います。

ここで𝐴_𝑁は前回も示したスケーリング係数と呼ばれる定数で、繰り返し回数𝑁が十分おおきければ 1.64676⋯に収束するものです。これは改めて以降で説明します。

1 ステップ演算を回転行列にするための準備

式(1)を変形していき、回転行列との関係を考えていきましょう。回転行列は前回の技術ノートのように、加法定理からそのなりたちを得ることができます。詳細はそちらをご覧ください。

まず CORDIC 漸化式、式(1)の第 2 式の項順を変更します。

式中の𝑖ビットの右シフト演算は、図 4 のように、1 ビット LSB 方向／右にシフトするごとに 1/2 になりますから、

と表すことができます。なおこの図は符号なし整数（unsigned int）で表記しています。

これを先の式(1)、式(6)に代入してみると

式(3)に示したように、𝐷(𝑖)は CORDIC の回転方向（極性）で、 𝑧(𝑖)がマイナスなら𝐷(𝑖) = −1、𝑧(𝑖)がプラスなら𝐷(𝑖) = +1とするだけの操作です。𝑧(𝑖)が減算される極性となります。以降の式中では、𝐷(𝑖)を単に「±」と表すことにします。

回転行列のかたちに変形していく

𝟐^𝒊と𝐭𝐚𝐧^−𝟏の関係を用いて式変形していく

ここで前回の技術ノートに示したテーブル値

を取りあげます。これは図 5 のイメージです。以降でも改めて示します。この図では𝑖 = 2を例としてみました。アークタンジェントなので、𝑠(𝑖)は角度量になります。これをタンジェントの式に変換してみると（これも図 5 に示してあります）、

になります。これからさきの式(9)は

と変換されます。

行列のかたちの式に変形していく

上記の式(12)を行列のかたちに変形してみると

ここで、1/ cos[𝑠(𝑖)]でくくってみます。なお、

なので

回転角度±𝒔(𝒊)の関数のかたちに変形していく

図 6 上のように

ですから、式(15)は

と書けます。さらに同じく図 6 下のように

ですから、この式の左辺を活用すると

となります。いまだに𝑥(𝑖)と𝑦(𝑖)がいまひとつイメージがつかないかと思いますが、図 7 のように、𝑥(𝑖)と𝑦(𝑖)はある角度𝜑と長さ𝐿をもつ「ベクトル𝑽」のX 成分（cos 成分）と Y 成分（sin 成分）

だと考えると、以降の説明の理解が進みます。ここで𝜑(𝑖)は x(𝑖), 𝑦(𝑖)が形成するベクトルの角度です。

なお𝑥(0)と𝑦(0)も、先の式(5)や図 3 において、長さが 1 で、角 度が 0°のベクトルの X 成分（cos 成分）と Y 成分（sin 成分）だと考えておくとよいと示しました（前回も図 17 で同様に示しました）。

回転行列のかたちになった

上記の式(19)は前回の技術ノートで示した回転行列

と全く同じで、𝛼 = ±𝑠(𝑖)と考え、

となります。そうすると先の式(19)は

と表されます。式(20)を用いて図 7 のベクトルの X 成分（cos 成 分）と Y 成分（sin 成分）の表記にしてみると、

と表されます。「ベクトル𝑽」の表記にすれば、

と表されます。

ステップごとの「演算」を例示してみる

それでは得られた式から、ステップごとの演算のようすをみてみます。

一番最初の演算ステップ、𝑖 + 1 = 1つまり𝑖 = 0の条件を考えてみましょう。ここでも𝑧(0) = +60°、つまり角度𝜃 = +60°を例に します。図 3 で示したように［式(5)再掲］

として𝑥(0)と𝑦(0)を「角度ゼロのベクトルの各成分」だとすれば、ここまでの議論から［式(23)を計算しています］、

右辺を回転行列の考え方を使って整理してみると、

として回転行列𝑹[±𝑠(0)]を使って、ベクトルの cos/sin 成分が角 度±𝑠(0) = +45°で𝜑(0)まで変化（回転）することになります［図 8］。ここで𝑠(0)は事前に用意するテーブル値

これも図 8に示しています。𝑧(0)が+60°なので「±」で表していた𝐷(𝑖)は正となり

になります。上記の式の𝑥(1)については、分母と分子にある cos 同士を消してしまう計算にも変形できますが、このままとしておきます。この 1/cos の項は、𝑥(0), 𝑦(0)そして𝑥(1), 𝑦(1)をベク トルの X 成分（cos 成分）、Y 成分（sin 成分）だとすれば、 CORDIC 演算 1 ステップで𝑖 = 0から𝑖 = 1にベクトルが変換されていくなかで、そのベクトルの長さ（ノルム）が

だけ長くなることに相当します。ベクトルの X成分［cos 成分である𝑥(0), 𝑥(1)］と Y 成分［sin 成分である𝑦(0), 𝑦(1)］も同じ比率で長さが長くなります。以降で𝑥(2)を求めるときに、式(28) で cos を残していた理由が、より理解できると思います。

つづいて𝑖 + 1 = 2つまり𝑖 = 1の条件を考えてみましょう。

ここで𝑠(1)は事前に用意するテーブル値

このようにステップごとに計算されていきます。これを図 9 に示しておきます。繰り返しますが、この図では𝑧(0) = +60°、つまり得たい答えとなる sin/cos の目的の角度を、𝜃 = +60°としています。

ステップごとの「変化率」を例示してみる

1 演算ステップごとの角度変化

この図 8 と図 9 から、𝑥(𝑖), 𝑦(𝑖)によるベクトルの角度𝜑(𝑖)が ［𝑧(0) = +60°を初期値、つまり目的とする角度を、𝜃 = +60°としたとき］

から

に、さらに

に変化していくようすが分かります。なお𝜑(3)も目的とする角度を𝜃 = +60°として計算してみると

ここで±の極性は残留角度誤差𝑧(𝑖)がゼロに収束していく方向として決まります。そのため𝜑(3)での𝑠(2)は引き算になっています。

1 演算ステップごとの大きさの変化率

sin/cos 成分それぞれの大きさ（ベクトルの長さ。ノルム）も

から

に変化していくことが分かります。

CORDIC 漸化式 1 演算ステップを一般化する

式(24)や図 7～図 9から、CORDICの漸化式、𝑖番目から𝑖 + 1番目 を計算する 1 ステップ演算は、

となります。ここで𝐿(𝑖)から𝐿(𝑖 + 1)へ大きさを変化させる係数

ただし

はあるものの、1 演算ステップにおいても、ここまで説明してきた回転行列の考え方そのままです。

つまり CORDIC の繰り返し演算の 1 回分は、デジタル数値計算で（足し算と引き算とビット・シフトで計算量を低減しつつ）ベクトルの sin, cos 成分を回転行列により角度±𝑠(𝑖) で回転させるもの（大きさを変化させる係数は除外して）になるのです。

回転角度は底辺が 2 倍ずつになる直角三角形の鋭角の角度

また式(22)は

となりますから、先にも示しましたが、この式は図 10 のように（前回の図 10 再掲）、𝑖がひとつ大きくなると底辺の長さが 2 倍になっていく、赤紫色／マゼンタ色で示した直角三角形の鋭角の角度𝑠(𝑖)が回転角となる、回転行列になります。

「2 倍」というのは図 4 で示したとおり 1/2 倍、 2倍の演算はデジタル処理においてはシフトで実現でき、CORDIC アルゴリズムとして、コンパクトな回路が実現できる点が根底にあるからです。

なおこの回転方向は、先にも示したように、𝐷(𝑖)でその±極性が決まります。𝐷(𝑖) は 1 ステップずつ漸化式で演算していくうえで、残留角度誤差量𝑧(𝑖)に𝑠(𝑖)を足し引きする「演算極性」になります。

スケーリング係数の元となる 1 ステップ演算での大きさの変化𝑨(𝒊)

式(23)では、回転行列の外に

という係数𝐴(𝑖)があります。これは CORDIC の 1 ステップ演算において、回転行列によるベクトル角度の回転以外の要素ともいえる、ベクトルの大きさ、そしてその成分 sin/cos の大きさが 増加する率です。

CORDIC で𝑁回演算された結果、1 回ごとの増加率𝐴(𝑖)をそれぞれ乗じたものが、先にも示した、また以下でも詳しく示していく「スケーリング係数」𝐴_𝑁になります。

このスケーリング係数𝐴_𝑁を求めるうえで、1 ステップ演算における±𝑠(𝑖)の項にある、±の極性は、図 6 下のとおり

になりますから、計算結果に±極性は無関係になります。

スケーリング係数𝑨_𝑵をより深く考察してみる

ここまでスケーリング係数𝐴_𝑁というものを説明しました。これは CORDIC で逐次的に 1 演算ステップごとに漸化式を計算していくうえで生じる係数（sin/cos の大きさの変化率）

の繰り返し積（総乗）であり

と表されます。ここで𝑁が十分大きければ 1.64676⋯に収束します。

この計算について考えてみましょう。[2]によると、逆三角関数の公式では

とあります。これは図形を描いてみると意味が分かります。これを図 11 に示します。そうすると

となり

となることも分かります。

実際に CORDIC を動かしてみた！

数式やブロック図によるイメージでの説明だけではつまらないので、実際に FPGA で CORDIC を動かし、DAC から信号を出力させてみたいと思います。使用する DAC は 16 ビットの AD5543BRMZ で、これで DDS として動作させます。AD5543 をご紹介しておきましょう。

AD5543 DDS、µSOIC-8 パッケージの 16 ビット DAC

https://www.analog.com/jp/ad5543

【概要】

AD5543/AD5553 は、高精度 16/14 ビット、ロー・パワー、 電流出力、小型の D/Aコンバータ（DAC）です。これらのデバイスは、±10V 乗算リファレンスを使い 5V 単電源で動作するように設計されています。

与えられた外部リファレンス VREF によって、フルスケール出力電流が決定されます。内部帰還抵抗（RFB）により、外付けオペアンプと組み合わせて、電圧変換用の R-2R と温度トラッキング機能を簡単に実現できます。（後略）

Arduino で使える FPGA 基板が欲しかった

CORDIC による cos/sin 値生成と、AD5543 へのデータ伝送には FPGA の使用が必要です。またもともと、Arduino ともつないでみたいという希望もありました。

まずは「Arduino シールドとなる FPGA 基板はないものか」と探してみました。ネットでいろいろサーチしてみると、Arduino 自体の代わりとなる FPGA 基板（Arduino ピン互換といえるもの）は販売されていますが、「Arduino シールドとなる」つまり、 Arduino と接続して Arduino を USB ブリッジなどに活用して FPGA 基板を制御する製品は見つけることができませんでした。

パネル de ボードの機能モジュールとして作ってみた

そこで P板.com の「パネル deボード」[3]の機能モジュールとして、「Arduino シールド FPGA 基板 IEEC-DI008A」を設計することにしました。「パネル de ボード」は P 板.com とアナログ・デバイセズで共同で立ち上げたサービスです。

図 12 は Arduino と、この FPGA 基板 IEEC-DI008A、そして別に開発した AD5543BRMZ を搭載した DAC 基板を 3 段重ねで接続したようすです。一番上に見える AD5543BRMZ の DAC 基板は 「ユーザ基板」となっており、ここは他の回路に乗せ換え、つまり FPGA で動作させたい回路を搭載できるようになっています。

FPGA 基板は「パネル de ボード」の規格に適合しています。また複数あるピンヘッダのレイアウトは、全て2.54mmピッチの穴位置になっているので、パネル de ボードを用いなくとも、市販のユニバーサル・ボードをこの FPGA 基板の上に「ユーザ基板」として構成することができます。

CORDIC の RTL を実装してみる

CORDIC は Web に落ちていた GNU の VHDL ソースを活用してみました（残念ながらご紹介するユニットはだいぶ前に作ったもので、この GNU の RTL ソースは現在では web に掲載されていません）。

このソースでのいちばんの心配は「DAC を駆動するための数値（たとえば 16 ビット DAC なら sin(-π/2) = －1 で 0x0000、sin(+π/2) = +1 で 0xffff）が出力されるソースになっているか」というところでした。

ソースを解析していくと、テーブル𝑠(𝑖)として数値を用意しているテーブル・ファイルでは、数値情報 32 ビットのうち、上位 2 ビットが小数点以上に相当するビットであること、残りの 30 ビットが小数点以下となっていること、ラジアン（π=3.1415) で設定されていることが分かりました。これを

sin(-π/2) = －1 で 0x0000

sin(+π/2) = +1 で 0xffff

になるようにテーブルの数値をスケーリングして修正しました。

今回使用した GNU の VHDL ソースは 32 ビット演算のものですが、一定ステップで変化する角度値カウンタに 24 ビット幅を用意し、これから 24 ビットの sin/cos 値を CORDIC アルゴリズムで生成します。その上位 16 ビットを AD5543BRMZ に転送し、アナログ正弦波信号を生成します。AD5543BRMZ は 2 個使用し、それぞれで sin と cos の 90°ずれた波形を出力するようにしました。

ちゃんと動いたぞ！

この時間波形、周波数を 1kHz にしたものを図 13 に示します。 図 14 は周波数スペクトラムのようすです。AD5543BRMZ は 16 ビット DAC ということで、高調波ひずみが-89dBc（キャリア比） に抑えられ、変なスプリアスもなく、非常に良好な性能であることが分かります。

まとめ

技術記事とは「分かっている人が知らない人に説明する」というのが一般的です。しかしたまに「分かっている人が、分かっている人にしか分からない説明をする」なんてこともありますし、「中途半端にしか分かっていない人が、知らない人に対して煙をまくように説明する」なんてこともあるでしょう。

前回、今回で説明した CORDIC を解説するネット上の記事も、上記の後半の 2 例のようなものもあるのでは？とか思っていたところです。そういうこの WEB ラボも「後半の 2 例」にならないよう気をつけねばなりません…。

そうそう、今しがたも企業の財務会計に関する簡単な疑問を解決しようと、とある会計系サイトの記事を見ていましたが、表現自体が無茶苦茶で、これでは分かっている人にも分からないぞ、と感じました（当然、私のような分かっていない人には全く理解できないということで…）。

分かりやすい解説、分かりやすい日本語表現、それぞれ難しいものです。