Python で CORDIC
概要
CORDIC (COordinate Rotation DIgital Computer) というアルゴリズムをPython で動かしてみた．
動機
FPGAで三角関数の計算をどうしたものかといろいろ調べて，フルデジタル無線機の信号処理という本でCORDICアルゴリズムの存在を知ったため，整数計算が可能か実験してみた． 説明
CORDICは関数電卓等で古くから使われてきたそうで，その一番の特徴はアルゴリズム内に乗算を必要としない点である．
それはそうと，何に使えるのかというと，以下の近似計算ができる．
三角関数($ \cos{z},$ \sin{z})
極座標変換($ \sqrt{x^2+y^2},$ \tan^{-1}{\frac{x}{y}})
しかもこの2項目の2つの計算についてそれそぞれの計算は同時に計算可能である．
つまり，FPGAなど並列演算が可能なデバイスを用いれば完全に同じタイミングで数値を算出でき，複素計算をするとき等に役立つ．
アルゴリズムの解説は，「CORDIC」で検索すれば丁寧な解説が出てくるのでここでは割愛して，整数で計算するときの値の変換方法に重点を置きたい．
ただ，私の理解としては，角度を二分探索的に収束させていく手法だと理解している．
定数について
定数について，ひとつなぜかどこにも説明がないもの($ K)があったため，自分で導出した．
CORDICで N bits の近似値を得たい場合，N個の$ \tan^{-1}{\frac{1}{2^n}},n=0\cdots (N-1)と定数 $ Kが必要である．この$ Kは以下の式で表される．
$ K=k_{N-1}, k_{n+1}=\sqrt{{k_{n}}^2+\left(k_{n}\frac{1}{2^n}\right)^2},k_{0}=1, n=0\cdots N-1
いきなりわけわからんと思うが，直角三角形の一辺が$ k_0=1でそれに直角な辺が$ k_0\frac{1}{2^{n=0}}=k_0=1なら，
斜辺は
$ \sqrt{{k_0}^2+k_0^2}=\sqrt{1^2+1^2}=\sqrt{1+1}=\sqrt{2}\approx1.41421356\cdots =k_{1}
である．
要はピタゴラスの定理である．この三角形の斜辺に対して$ \frac{1}{2^n}倍の辺をもつ直角三角形をかく．
今度の斜辺は，
$ \sqrt{{k_1}^2+\left(k_1\frac{1}{2}\right)^2}=\sqrt{2+\left(\sqrt{2}\times\frac{1}{2}\right)^2}=\sqrt{2+2\times\frac{1}{4}}=\sqrt{2.5}\approx1.581138830\cdots =k_{2}
となる．
これを繰り返すと，$ K=k_{N-1}\approx1.6467602, in N=12が求まる．$ Nが大きければ，だいたい似た数字になる．ここで説明した直角三角形の直角を構成する長辺側の斜辺との角度が$ \tan^{-1}{\frac{1}{2^n}},n=0\cdots (N-1)というわけである．
table:k_n
0 1.4142135623731
1 1.58113883008419
2 1.62980060130066
3 1.64248406575224
4 1.64568891575725
5 1.64649227871248
6 1.64669325427364
7 1.6467435065969
8 1.64675607020488
9 1.64675921113982
10 1.64675999637562
11 1.6467601926847
12 1.64676024176197
13 1.64676025403129
14 1.64676025709862
15 1.64676025786545
16 1.64676025805716
17 1.64676025810509
全く関係ないが，夏目漱石の「こゝろ」という小説にKという登場人物が出てくる．『精神的に向上心のない者は馬鹿だ』はとても印象的である．しかしながら，時には「精神的な向上心」から解き放つことにより得られる「向上」もあること，あまり根詰めても仕方がなく，本能に任せることも大切であろう．あまりオチを語りたくないが，実際Kはこの言葉によって自分の首を絞めることになる．．．
閑話休題，CORDICによってNbit精度に対して上で説明したN+1個の定数をメモリに格納して，アルゴリズムをN回実行することによって求められるのである．
ここで，FPGAで計算することを考える．
角度の表現
象限毎に上位2bitを割り当てる
具体的には例えば12bitなら
$ 0 \mathrm{rad} = \mathrm{0x000}
$ \pi \mathrm{rad} = \mathrm{0x400}
$ 2\pi \mathrm{rad} = \mathrm{0x800}
$ 3\pi \mathrm{rad} = \mathrm{0xC00}
定数Kの整数値換算