三角関数の高速化


オドメトリなどを使うことを考えると,\sin,\cosなどの三角関数を実装する必要がある.SHなどの浮動小数点演算ユニットを持たない組み込みコントローラを用い,1~10msec周期でオドメトリを実行するには高速な三角関数の実装が必須です.そこで,本稿ではテーブル引きとテイラー展開を用いた三角関数の実装をまとめる.

公式

準備として以下に本稿を読むに必要な数学知識をまとめる.

三角関数の公式

以下のような公式が知られている.

(1)    \begin{align*} \sin (x+\frac{\pi}{2})&=\cos x\\ \sin (-x)&=-\sin x\\ \sin (\pi-x)&=\sin x\\ \sin (a+b)&=\sin a\cos b+\cos a\sin b\\ \end{align*}

正弦関数のテイラー展開

ある定数値a\in\mathbb{R}\displaystylea\displaystyleと十分近い値x\displaystyleに対し,以下のような近似式が成立する.

(2)    \begin{align*} \sin (x)&=\sin(a)+(x-a)\cos(a)-\frac{1}{2!}(x-a)^2\sin(a)-\frac{1}{3!}(x-a)^3\cos(a)+\cdots \end{align*}

x-a\displaystyleは大きくなるほど,誤差が大きくなる点に注意しよう.

定義域の制限

テーブルを使って三角関数を求めることを考えよう.より細かい分解能でテーブル用意すれば,精度は向上する.しかし,使用するROMの容量が増大する問題がある.そこで,三角関数の定義域を制限し,不要なテーブルを用意しないことを考える.

\sin (x+\frac{\pi}{2})=\cos x\displaystyleより,高速な\sin\displaystyle関数を用意すれば,加法演算1回で高速な\cos\displaystyle関数を用意することができる.

\sin (-x)=-\sin x\displaystyleより,x\geq 0\displaystyleの場合のみを考えればよい.さらに,\sin (\pi-x)=\sin x\displaystyleより,\frac{\pi}{2}\geq x\geq 0\displaystyleの場合のみテーブルを用意すればよい.以上より,\frac{\pi}{2}\geq x\geq 0\displaystyleの範囲で\sin x\displaystyle関数のテーブルを用意すれば十分なことが分かる.たとえば,\sin (0),\sin (1\cdot 2\pi/360),\sin (2\cdot 2\pi/360),...,\sin (90\cdot 2\pi/360)\displaystyleの値をあらかじめテーブルとして用意する.そして,\sin (10.5\cdot 2\pi/360)\simeq\sin (10\cdot 2\pi/360)\displaystyleのように,もっとも引数に近い値のテーブル値を近似値として出力すれば良い.

本稿で紹介した手法で,作成するテーブルの規模は小さくすることが出来た.しかし,計算精度とテーブル規模に関するトレードオフ問題は回避できていないことが分かる.

テイラー展開を使ったテーブル値の補正

前節では,\sin (10.5\cdot 2\pi/360)\simeq\sin (10\cdot 2\pi/360)\displaystyleのように,テーブルの量子化幅に満たない部分を無視する近似を行っていた.そこで,量子化幅は0に近いことに注目する.

(3)    \begin{align*} \sin (x)&=\sin(a)+(x-a)\cos(a)               -\frac{1}{2!}(x-a)^2\sin(a)-\frac{1}{3!}(x-a)^3\cos(a)+\cdots\\ [crayon-59061598549cf617951091/] \end{align*}

より,\sin a,\cos a\displaystyleをテーブルで求め,この近似式に代入すれば,テーブルを引くだけよりかは高精度で演算を行うことができる.ただし,\cos a\displaystyleをテーブル引きだけで求めるためには,任意のテーブルの引数a\displaystyleに対し,a+\frac{\pi}{2}\displaystyleがテーブルの引数として用意されている必要がある.これは,a=0,\pi/2\displaystyleがテーブルの引数として用意されており,分割幅が一定であることと等価である.

多角形近似

テイラー展開を使った補正近似の他に,多角形近似による補正が考えられる.幾何的に考えると,理論上はテイラー展開の方が優位なはず.そのうち検討.

固定小数点方式で計算

三角関数の計算では,極端に大きい数や小さい数を扱わない.そのため,固定小数点方式(プログラミングの小ネタ#固定小数点方式で解説)で演算を行うことが考えられる.採用した場合,double演算は浮動小数点方式から固定小数点方式への変換部分のみで必要となる.

実装

本稿では,SH7147のような,int型の演算は十分早いが,double型の演算がきわめて遅いプロセッサを想定したコーディングを行う.たとえば,int型の演算が10回増えようが,double型の演算が1回でも減れば計算量の削減とみなす.

以下は128分割テーブルの場合の実装である.

テーブルデータの作成

以下のようにしてテーブルを作る.

テーブルを引くだけの実装

前節のプログラムで作られたテーブルを用い,以下のように三角関数を実装することができる.

254/M_PIを定数として別に保持しておけば,mysin関数はdouble演算を2回のみ行うことがわかる.このプログラムの実行結果を以下に示す.真値とテーブル引きの間では,最大0.005程度の誤差がおこることが分かる.

本質的には三角関数のテーブル引きによる高速化とあまり変わらない.

テイラー展開を用いた補間をする場合

以下に1次のテイラー展開でテーブル値を補間した実装例を示す.

254/M_PIを定数として別に保持しておけば,mysin関数はdouble演算を6回行うことがわかる.このプログラムの実行結果を以下に示す.真値とテーブル引きの間では,最大0.00001程度の誤差がおこることが分かる.このように,劇的な精度の向上が確認できた.なお,テイラー展開を2次で打ち切ると,double演算が9回になるものの,誤差が0になった.(printfで表示できる限界まで正確になる)

sh7147に実装してみたところ,テーブル引きのみだと実行時間を100usec以下に抑えることができた.テイラー展開による補正を加えると100usec以上の計算時間がかかった.

美しいコーディング

ロボット工作研究室wikiに美しいコード例があります.

付録

マクローリン展開(テイラー展開の系)

十分小さなx\displaystyleに対し,以下のような近似式が成立する.

(4)    \begin{align*} \sin (x)&=x-\frac{1}{3!}x^3+\frac{1}{5!}x^5+\cdots\\ \cos (x)&=1-\frac{1}{2!}x^2+\frac{1}{4!}x^4+\cdots \end{align*}

参考文献

外部リンク