fract16型の余弦波と正弦波を生成するC++のクラスです。
まずはクラス宣言から。
// すべてのDDS型の先祖になる抽象型
class COscillator{
protected:
unsigned int state[2];
// state[0]:位相
// state[1]:周波数
// 位相は32bitで2πをあらわす
// 周波数は1サンプルごとの位相進角をあらわす。
public:
COscillator( unsigned int freq );
virtual void setFreq( unsigned int freq ); // 位相の強制変更
virtual void setPhase( unsigned int phase ); // 周波数の強制変更
virtual void run ( fract16 sin[], fract16 cos[], unsigned int count )=0;
};
// 三角関数の加法公式によるDDS
class CDDS: public COscillator{
public:
CDDS( unsigned int freq );
virtual void run ( fract16 sin[], fract16 cos[], unsigned int count );
};
COscillatorは抽象型なので実装するのはCDDSの方です。
CDDS::runメソッドはアセンブリ言語によるルーチンを呼び出すだけです。ルーチンのアルゴリズムは2191空挺団の関数表を工夫するで説明した2分割表と三角関数の和の公式を使うものです。今回、せっかくBlackfinを使うのでデュアル・ロードとデュアル演算を活用してみました。
また、32ビット・ロード１回でcosとsinの表の値を取得できるよう表の形を工夫しています。ループの中にnop;が二回続けて現れますが、これはストール・サイクルを明示的に埋めるものです。このようにしておけばあとで最適化を行うときに便利です。なお、3ヶ月間無料で使えるお試し版ならば、パイプラインのどこでどのようなストールが発生しているのか簡単に知ることができます。
// parameter list
// 符号なし整数を位相として受け取り、振幅を整数で返す
#define COS 1 // 余弦波出力テーブル
#define SIN 2 // 正弦波出力テーブル
#define COUNT 3 // 生成するサンプル数
#define STATE 4 // 位相の内部状態と周波数を含むunsigned int * 2の構造体
// void ddsbody( fract16 cos[], fract16 sin[], int count, unsigned int satate[] );
// cossinの要素はそれぞれ32ビット。各ビットは下位にcos、上位にsin
// の値を符号付16ビット固定小数点の値として収めている。
.section L1_data_a;
.byte2 cos_sina[512]="cossina.txt";
.var extmask=0x1808; // pos:24, length:8
.section L1_data_b;
.byte2 cos_sinb[512]="cossinb.txt";
// M0 : cos_sina 上位8ビットによる配列ベース
// M1 : cos_sinb 下位8ビットによる配列ベース
// I2 : 上位8ビットによる配列オフセット
// I3 : 下位8ビットによる配列オフセット
// I0 : cos
// I1 : sin
// R0 : phase
// R1 : freq
// R7 : 0x00ff
.section program;
.global _ddsbody;
_ddsbody:
link 0;
[--sp] = (r7:4, p5:3); // レジスタ退避
m0.H =cos_sina;
m0.L =cos_sina;
m1.H =cos_sinb;
m1.L =cos_sinb;
p0=[fp+20]; // 内部状態のアドレス
p5=r2; // count
r7=0x01fe; // r7 : 0x00ff
i0=r0; // p5 : COS出力配列
i1=r1; // i1 : SIN出力配列
p4.L=extmask;
p4.H=extmask;
r7.L=0x1008; // pos:16, length 8;
loop DDS lc0=p5;
r0=[p0]; // phase
r1=[p0+4]; // freq
r2.H = 0;
r6 = [p4];
loop_begin DDS;
// 16bit位相量を8bitのインデックスに分解
r2 = extract(r0, r6.L)(z); // 上位8bit
r3 = extract(r0, r7.L)(z); // 下位8bit
r2 = r2 << 2; // オフセットに変換
r3 = r3 << 2; // オフセットに変換
i2 = r2; // 上位8bitによるオフセット
i3 = r3; // 下位8bitによるオフセット
r0 = r0 + r1(ns); // 位相のアップデート
nop;
nop;
i2 += m0; // 上位8bitによる要素アドレス
i3 += m1; // 下位8bitによる要素アドレス
r4 = [i2] || r5 = [i3]; // cossin(h), cossin(L)
// cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
// sin(a+b)=sin(a)cos(b)+cos(a)sin(b)
a1 =r4.L*r5.L, a0 =r4.H*r5.L; // a1:cos, a0:sin
r6.H=(a1-=r4.H*r5.H), r6.L=(a0+=r4.L*r5.H); //.H:cos, .L:sin
w[i0++] = r6.H;
w[i1++] = r6.L || r6 = [p4];
loop_end DDS;
// 更新したPHASEを戻す
[p0] = r0; // phase
(r7:4, p5:3) = [sp++]; // レジスタ復帰
unlink;
rts;
_ddsbody.end:
ループボディーのサイズは16行です。デュアル・アクセスが効果的に行えるようにメモリー配置を工夫すれば、16サイクルで実行できます。600MHzで動作させると、このルーチンはほぼ40MHzの信号を生成できます。