前回確定したアルゴリズムを使って実用的な関数に仕立て、クラスで包んでオブジェクト化します。
新たに作った関数は以下のとおりです。
// void ddsbody(
// int sin[], //
// int cos[], // cosが不要のときは0を渡す
// int step, // 配列内に値を格納する時の間隔
// int count, // 計算する要素数
// int * phase // phaseは内部位相(L,H)と周波数(L,H)の4ワードからなる。
// それぞれは32bit数で、位相の上位16びっとを表引きに使う
// 順番は下位からPL,PH,FL,FH。
// );
sinとcosを同時に計算し、与えられた配列に格納していきます。cosがNULLのときには計算だけ行って数値は返しません。値はcount組計算し、stepワード毎にメモリーに配置します。phase引数は4ワードを渡します。最初の2ワードは32bitの内部位相です。次の2ワードは32ビットの周波数(サンプルごとの位相増分)です。いずれもリトルエンディアンで格納します。
ルーチン本体を以下にしめします。ループ内部に20ステップと比較的重い関数ですが、行っている計算一つ一つは非常に単純です。
// parameter list
#define SIN 1
#define COS 2
#define STEP 3
#define COUNT 4
#define PHASE 5
// 符号なし整数を位相として受け取り、振幅を整数で返す
#include <asm_sprt.h>
.section/dm data1;
.var sina[256]="sina.txt";
.var sinb[256]="sinb.txt";
.var cosa[256]="cosa.txt";
.var cosb[256]="cosb.txt";
// I0 : sina
// I1 : cosa
// I2(R) : sinb
// I3(R) : cosb
// M1 : indexH
// M3 : indexL
// I6 : sin
// I7(R),SI : cos
// M7 : step
// AX1 : phaseH
// AX0 : phaseL
// AY1 : freqH
// AY0 : freqL
// AF : 0x00ff
.section/pm program;
.global _ddsbody;
_ddsbody:
entry; // スタックフレーム確立
pushs( i2 );
pushs( i3 );
pushs( i7 );
ax0=0x00ff;
af=ax0+0; // af : 0x00ff
get_arg(i7, PHASE); // 引数を取得
m7=1;
ax0=dm(i7+=m7); // phaseL
ax1=dm(i7+=m7); // phaseH
ay0=dm(i7+=m7); // freqL
ay1=dm(i7+=m7); // freqH
get_arg(i6, SIN);
get_arg(i7, COS);
get_arg(m7, STEP);
get_arg(cntr, COUNT);
i0=sina;
i1=cosa;
i2=sinb;
i3=cosb;
dis m_mode; // 固定小数点モード
si=i7; // cos
ar= pass si; // cos 引数は0か?
if ne jump doloop; // 0でなければそのまま
i7=i6; // 0のときは、sinと同じアドレスにする
doloop:
do calc until ce;
// 16bit位相量を8bitのインデックスに分解
ar=ax1 and af; // phaseH and 0x00ff
m3=ar; // 下位8bit;
sr = lshift ax1 by -8(lo);
m1=sr0; // 上位8bit;
// 位相のアップデート
ar=ax0+ay0; // PHASE_L += FREQ_L
ax0=ar;
ar=ax1+ay1+c; // PHASE_H += FREQ_H
ax1=ar;
// 正弦、余弦表をアクセス
mx0=dm(i0+m1); // sin(h)
my0=dm(i1+m1); // cos(h)
mx1=dm(i2+m3); // sin(l)
my1=dm(i3+m3); // cos(l)
// cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
mr=my0*my1(ss);
mr=mr+mx1*mx0(rnd);
sat mr; // 計算値を飽和させる
dm(i7+=m7) = mr1; // 余弦を格納
// cos 引数が0のとき、余弦はsin引数に代入される。
// 次のステップでsin引数は上書きされるので、結局
// cos引数が0なら値は捨てられる
// sin(a+b)=sin(a)cos(b)+cos(a)sin(b)
mr=mx0*my1(ss);
mr=mr+mx1*my0(rnd);
sat mr; // 計算値を飽和させる
calc: dm(i6+=m7) = mr1; // 正弦を格納
// 更新したPHASEを戻す
m7=1;
get_arg( i7, PHASE );
dm(i7+=m7)=ax0; // PHASE_Lを戻す
dm(i7+=m7)=ax1; // PHASE_Hを戻す
ena m_mode; // 整数モードに戻す
pops( i7 );
and_pops( i3 );
and_pops( i2 );
exit; // スタックフレームを廃棄
ラッパークラスはDDSをカプセル化します。クラスに包み込むことで内部変数をユーザーから遮蔽し、複数インスタンスの取り扱いを容易にします。なお、位相や周波数を変更できるように設計しました。
namespace dds{
// すべてのDDS型の先祖になる抽象型
class CDDS{
protected:
int phase[4]; // phaseは内部位相(L,H)と周波数(L,H)の4ワードからなる
// それぞれは32bit数で、位相の上位16ビットを表引きに使う
// 順番は下位からPL,PH,FL,FH
// 位相は符号なしで、32bitで2πをあらわす
// 周波数は1サンプルごとの位相進角をあらわす。
public:
CDDS( int freqH, int freqL );
virtual void setFreq( int freqH, int freqL ); // 位相の強制変更
virtual void setPhase( int phaseH, int phaseL ); // 周波数の強制変更
virtual void run ( int sin[], int step, int count )=0;
virtual void run ( int sin[], int cos[], int step, int count )=0;
};
// 三角関数の加法公式によるDDS
class CSinCosDDS: public CDDS{
public:
CSinCosDDS( int freqH, int freqL );
virtual void run ( int sin[], int step, int count ); // 正弦波のみの生成
virtual void run ( int sin[], int cos[], int step, int count ); // 正弦波と余弦波の生成
};
};
#include "dds.h"
extern "C" void ddsbody(
int sin[], // 1
int cos[], // 2 cosが不要のときは0を渡す
int step, // 3
int count, // 4
int * phase // 5 phaseは内部位相(L,H)と周波数(L,H)の4ワードからなる
// それぞれは32bit数で、位相の上位16びっとを表引きに使う
// 順番は下位からPL,PH,FL,FH
);
namespace dds{
CDDS::CDDS( int freqH, int freqL )
{
setPhase( 0, 0 );
setFreq( freqH, freqL );
}
void CDDS::setFreq( int freqH, int freqL )
{
phase[2] = freqL;
phase[3] = freqH;
}
void CDDS::setPhase( int phaseH, int phaseL )
{
phase[0] = phaseL;
phase[1] = phaseH;
}
CSinCosDDS::CSinCosDDS( int freqH, int freqL ) : CDDS( freqH, freqL )
{
}
void CSinCosDDS::run( int sin[], int step, int count )
{
ddsbody( sin, 0, step, count, phase );
}
void CSinCosDDS::run( int sin[], int cos[], int step, int count )
{
ddsbody( sin, cos, step, count, phase );
}
};
比較的小さなテーブルで精度の高い正弦波を生成するプログラムに挑戦してみました。関数表を2分することで大幅にメモリーを節約できます。3分すればさらに小さくなりますが、すでにレジスタを使い切っていますから、この先はオーバーヘッドも増大します。
計算に必要な時間はsin/cos組で1サンプルあたり20サイクルです。これを大きいと見るか小さいと見るかですが、オーディオ帯域ではサンプル周波数48KHzのとき3000サイクルを1サンプルの処理時間に割り当てることができます。このことを考えれば1サンプルあたり20サイクルは十分小さな値だと言えます。
プログラムのダウンロードはこちら: ダウンロード 12KB
最後にsin/cosを計算したときのスクリーンショットを紹介します。
図1 スクリーンショット
次は⇒音を聞いてみる