あけましておめでとうございます。今年もよろしくお願いします。オーディオ・フレームワークの開発が一段落してきたので、今年はアプリケーションとアルゴリズムで遊ぶ年にしたいと思っています。
さて最終的な目標であるSDRを作るには、複素発振器が必要になります。これはcos(),sin()関数を同時に計算するものです。
そこで複素発振器用アルゴリズムとして32bit固定小数点三角関数をBlackfinに実装しました。このプログラムは0.75kBの内部データSRAMで動作し、角分解能、振幅分解能がいずれも32bitです。誤差は最大±3LSBに収まっていると考えられ、精度はおおよそ29bitです。開発にはGNUツールチェーンおよびGDB内蔵のシミュレータを使いました。
大まかなアプローチ
大まかな方法として、テーブル参照と三角関数の加法公式を使ったテーブル分割を使うことにしました。これは以前16bit版を作った時に使用したものです。
三角関数をテーブル参照で作ることは誰でも考え付く方法ですが、必要な精度を満たすテーブルを用意できるかどうかは別の問題です。たとえば、角分解能が16bitのテーブルは64kエントリが必要になります。SDRAMになら余裕で搭載できますが、多くの場合これはあまりいい方法とは言えません。
そこで、三角関数の加法公式を使ってテーブルを２分します。64kエントリのテーブルであれば、角度の上位8bitを担当する256エントリと角度の下位8bitを担当する256エントリに分割します。そうして上位エントリで2πの範囲をカバーし、下位エントリで2π/256をカバーします。こうすると、上位エントリのcos/sinテーブルと、下位エントリのcos/sinテーブルを使って、全周16bit分解能64k点による三角関数を生成できます。この方法は乗算時に雑音が混入することから、完全テーブル方式に比べるとややS/Nが劣ります。しかしながら、テーブルサイズはわずか1/128に収まりますので、強力なメリットといえます。たとえばDSPであれば内蔵SRAMに搭載できます。
今回はこの方法を拡張することにしました。目標は32bit精度の角分解能ですので、まともにテーブルにすれば4Gエントリになります。冗談としては面白いのですが手に負えません。2分割してもそれぞれ64kエントリなのでやはりblackfin向きではなく、多段テーブル方式などから考えることにしました。
微小角cos/sinの計算に多項式近似を用いる
32bit角を4分割して256エントリのテーブルを4つもつ方法も考えられます。この場合必要なサイズは8kBです。大きすぎますし芸がありません。
当初角度の上位16bitのみテーブル化することを考えました。実際には符号ビットがありますので17bitテーブル化と同じです。で、残りのLSB15bitについては次の方法を考えました。
- cos( x ) は1とみなす。
- sin( x ) は、Axで代用する(Aは係数)。
よく知られているようにが非常に小さいときには sin( x ) ≒ x とみなすことができます。また、cos( x ) ≒ 1とみなすこともできます。
これをもとに精度評価したところ、以下のようなことが分かりました。
- sin( x )は非常に広い範囲で Ax とみなすことができる。具体的には角度の下位19bitを Axで代用しても、sin( x ) からの誤差は0.5LSB以下。
- cos( x )を1で代用すると、sin/cosの演算誤差は許容できないほど大きくなる。
- cos( x )を近似するときには2次間数にすると誤差が小さくなる
これらのことから、結果的に以下のような構成をとることにしました。
- 符号ビットと続く1bitで4象限を表す。
- 上位テーブルで1象限をカバー。角分解能は5bit 32エントリ。
- 下位テーブルを6bit 64エントリとする。
- 残り19bitのエントリを関数近似で代用する。sin() は線形近似、cos()は2次関数による近似。
大まかな評価結果
全周を128k点に分割してそれぞれのcos/sinを計算し、Libre Office Calcで解析しました。その結果、Libre Officeによる計算結果と比べると次のようになります。
- cos()の誤差の範囲は –5.2 .. 3.0、平均は –1.2、標準偏差は 1.4
- sin()の誤差の範囲は –5.2 .. 5.2、平均は –1.1、標準偏差は 1.8
全点検査したわけではありませんが、ざっくり言って雑音振幅は3bitに収まっていると考えてよさそうです。これが仮に5bitになったとしても、現在考えている用途ではまったく問題になりません。
オフセットが1LSBほどありますが、気にしないことにします。この辺りは丸め誤差や正負のオフセット(-1はあるが+1はない)などのほかに、手作業ゆえの誤差が入っている可能性も否めませんが、商用は考えていませんのでこれで妥協します。
ソースコードはそのうちに公開します。