第7回

Sin 波を作る――テーブル参照のテクニック

堀江 誠一

音声処理に欠かせない関数として，Sin と Cos がある．メモリ容量や演算速度に余裕のある環境であれば，巨大なテーブルや複雑な演算式で計算できる．しかし，組み込み機器では限られたメモリ容量と CPU パワーで現実的な精度の演算結果を得る必要がある．今回は DSP で Sin と Cos を実現するための方法について解説する．
（編集部）

Sin 波と Cos 波は，発振器のみならず信号処理を行う上仕組みはすぐに思い付くアイデアですが，あっという間にであちこちに登場する重要な信号です．ディジタル信号処壁にぶつかってしまいます．それは，精度を高めようとす

理に限らず，実際の電子回路においても発振器の設計は重るとテーブルが大きくなり過ぎるということです．円周を要なポイントとして扱われます．

16 ビット精度で分割すると，必要なテーブルの大きさはSin（正弦）関数，Cos（余弦）関数といった三角関数は広64K エントリになります．これでは話になりません．

く使われるので，短い時間で精度良く計算することが求めSDRAM に格納することが前提ならば，不可能な数字ではられます．三角関数の計算方法にはいろいろありますが，ありませんが，容量が小さな高速内部メモリに置くことは今回はテーブルを使うものについて紹介します．また，三できません．SDRAM に置いた場合もキャッシュは使える

角関数を作りながら，関連する話として丸め誤差についてものの，かなりの速度低下を甘んじて受け入れなければなも考えてみます．

りません．

最後に，出来上がった発振器を使って少し遊んでみま
しょう

注1

しかし，三角関数の公式の一つを使えば，巨大なテーブルを分割して小さくし，内部メモリに配置できます．その．

方法を紹介しましょう．

1. テーブル参照による三角関数の計算● 三角関数の和の公式を使ってテーブルを縮小三角関数については，高校の数学の時間に教わった以下

テーブル参照による三角関数の計算について考えてみましょう．
テーブルに三角関数を入れておいてそれを引く，というの公式を覚えているかと思います．
cos( a + b) = cos( a)cos(b) − sin( a)sin(b)…………………（1）sin( a + b) = sin( a)cos(b) + cos(b)sin( a)…………………（2）

これらの公式を使うと，64K エントリの巨大なテーブル
円周の256分の1をさらに256等分したテーブルsi_l[256]とco_l[256]を作る
円周を256等分したテーブル
si_h[256]とco_h[256]を作るを 1K エントリに圧縮できます．その原理を図 1 を見ながら説明しましょう．ポイントは式（1）の a+b を使って，三角関数の引き数を分割するというアイデアです．
先に仮定したように円周を 16 ビット精度で表すとしま拡大す．このとき，図のように，円周を 256 分割する荒っぽいテーブルと，円周の 256 分の 1 をさらに 256 分割する細かいテーブルがあると考えてください．それぞれのエントリのインデックスを A，B とすると，それぞれのインデック
スは円周を表す 16 ビット・インデックスの上位，下位にsi_h[A]*co_l[B]+co_h[A]*si_l[B]を計算すれば，sin((A*256+B)*2*π/65536)を得ることができる．
必要なテーブルのサイズは256*4

図 1 テーブルによる三角関数の計算

注 1 ：サンプル・プログラムは，CQ 出版社の Web サイトからダウンロードできる．http://www.cqpub.co.jp/interface/download/

166 KEYWORD ―― Blackfin，Sin 波，発信機，シングル・トーン，DSB 変調，AM 変調，FM 変調，
フランジャ

Oct. 2007

リスト 1 Scilab によるテーブル生成

16ビット符号なし整数（0〜65535）で0〜2πを表現する

index=[0:255]*2*%pi/256;

X：

co=cos(index)*32768;
write( '¥hco.txt', int(co)', '(i8)');

上位8ビット（A）

下位8ビット（B）

co_h[A], si_h[A]

co_l[B], si_l[B]si=sin(index)*32768;
write( '¥hsi.txt', int(si)', '(i8)');si_h[A]*co_l[B]+co_h[A]*si_l[B]⇒ sin(X)index=[0:255]*2*%pi/65536;
co=cos(index)*32768;
write( '¥lco.txt', int(co)', '(i8)');co_h[A]*co_l[B]-si_h[A]*si_l[B]⇒ cos(X)si=sin(index)*32768;
write( '¥lsi.txt', int(si)', '(i8)');

図 2 添え字の分解のテーブルを引きます（図 2）．最後に先の三角関数の和の相当することになります．
この分割したインデックスでおのおののテーブルを引き，式（1）を使って値を合成すれば，16 ビット精度で円周を分公式を使って精密な正弦関数の値を計算します．

● 誤差を x86 上のシミュレーション機能で評価する割した三角関数を作り出すことができます．このとき必要関数を作ることができたので，その誤差を評価してみまなテーブルは 256 エントリのものが四つなので，合計で 1024しょう．プログラムをシミュレータで走らせ，余弦関数のエントリとなります．16 ビット DSP なら 2K バイトです．

計算結果を printf()文を使って 65536 点すべての計算結● テーブル参照によるプログラム

果をファイルに出力します．

四つのテーブルはどのようにして作っても構いませんが，
注2

この計算には VisualDSP++ の compiled simulation 機能を使ってみまを使いました．この機能は Blackfin プロセッサのアセンブした．リスト 1 にその様子を示します．生成されるファイリ・プログラムを x86 のアセンブリ・プログラムに変換しルは，整数が 1 行に一つある形式です．これに行末のコンて実行するもので，シミュレーション速度が 100 倍程度にマを追加して，C 言語で読める形式にします．こうして作っなる非常に高速なシミュレーション機能です．評価キット

たテキスト・ファイルをソース・コード中の #include 文のライセンスでは，インストールから 90 日間だけ使えるで読み込み，配列の初期値として利用します．

機能です．計算だけならば EZ-KIT lite 上でやっても構わ

今回は科学技術計算用ソフトウェア Scilab

shortfract 型を使っているので，整数ではなく固定小ないのですが，EZ-KIT Lite は printf 関数の実行が遅い

数点数表記のまま読み込ませればいいのですが，あえて整ため，結果的に compiled simulation が最速になります．

数を使っているのには理由があります．Analog Devices 社

計算結果を Excel に取り込んで cos()の計算結果と比較のコンパイラには制限があり，shortfract 型変数を静的したのが図 3 です．縦軸は誤差を LSB で表しています．一に初期化すると，C++ の初期化用のコードが生成されて見して，正負に振れる大きな誤差があることが分かりましまうからです．これは小規模のテーブルでは問題あります．実はこの誤差は関数計算の過程ではなく，係数テーブせんが，今回のようにそれ相当の規模になるとプログラムルを作る過程で入り込んだものです．

が膨れ上がってしまうのです．

● より誤差の少ない表を作る

Analog Devices 社のツール・グループには改善を提案

このプログラムの係数を作るために使った Scilab のプロしていますが，残念ながら現時点では大きなテーブルの初グラムは，計算結果を整数化するために int()関数を使っ期化は整数を使ってごまかすしかありません．ています．この関数は小数部を単純に切り捨てるため，誤プログラムの流れは至極単純です．与えられた引き数差の絶対値が最大で 1 になります．話はそれだけではすみ（正弦波の位相にあたる）は 16 ビット符号なし整数です．

ません．切り捨ての方向は 0 に向かうため，正負で値の切先に説明したように円周は 65536 分割し，16 ビットで全周り捨て方向が違うのです．図の誤差に途中から大きなオフを表します．このように決めた引き数を上位と下位の 8ビットに分割し，それぞれを添え字として Sin および CosOct. 2007

注 2 ：Scilab に関しては Interface 誌 2006 年 9 月号の特集記事で扱っている．

167