nabeの雑記帳

漸化式*1によりsin信号を生成する方法を紹介します。これにより、マイコンのような機器でも高速にsin信号を生成することができます。FPGA等にも容易に応用できます。少し工夫すると円も描けます。

• 特徴
• 生成アルゴリズム
• サンプル
• 原理
  • 導出
• 参考文献
• 注意

※高精度の演算には使えません。

■特徴

• 任意の周波数、任意のサンプリング周波数でsin波を生成できます。
• 初期値のみ、あらかじめ計算しプログラムに埋め込むか、別の方法で計算する必要があります。
• 2次のIIR型フィルタにより生成するため、時間が経つと計算誤差が累積します。

生成周波数／サンプリング周波数が整数ならば、誤差の累積は非常に小さくて済みますが、それだったらテーブルで良いような気もします（でもメモリの少ないマイコンなら多少は有利かな）。

■生成アルゴリズム

至って簡単です。

初期値計算。

y[0] = sin(phi);
y[1] = sin(phi + 2*PI*f/k);
ar1  = 2*cos(2*PI*f/k);

生成式(n>2)。

y[n] = ar1*y[n-1] - y[n-2];

誤差が累積します。float型でためしたところ、（パラメータにもよりますが）1000回の計算で0.00012程度、10000回の計算で0.0012程度の誤差が出ました。doubleで生成すれば相当低誤差で使えそうです。

もしfがkで割り切れるなら１周した段階で（予め計算してある）初期値に戻してあげる方が良いです。割り切れない場合は、（もしsinを計算可能な環境ならば）適当な回数を演算したところで初期値計算を再度行うと良いでしょう。

永遠に計算し続けると、誤差の累積により発散するか0に収束すると思われます。

■サンプル

Excelファイル。

Cのソース。

#define PI 3.14159265358979323846
float ar1,yz1,yz2;
sin_init(float freq, float sps, float phi) {
	yz1 = sin(phi);
	yz2 = sin(phi+2*PI*freq/sps);
	ar1 = 2*cos(2*PI*freq/sps);
}
float sin_get() {
	float g;
	g   = yz1;
	yz1 = yz2;
	yz2 = ar1*yz1 - g;
	return g;
}
main() {
	sin_init(1000, 10000, 0);
	while(1)
		printf("%f\n", sin_get());
}

■原理

2次IIRフィルタを振動させているだけです。Z平面上での極を単位円上に配置してパラメーターを算出しています。

右の図に示すようにIIRによる振動現象は、単位円上を動く点と捉えることができます*2。きまった角度φずつ単位円の上をスライドする点を観測して生成されるものがsin波です。*3図で全体を100sample/secとして、各ポイントに対応する周波数を書きました。左に回転量が増えるほど周波数が高くなります。*4

IIRフィルタでの回転量を決めるのは、右図に示した「極」です。IIRフィルタで計算すること*5は複素平面上で極を単純にかけ算することに等しいのです。複素数の授業でやったと思いますが、ある点に対して単位円上の点pをかけ算すると、その点は点pの角度だけ回転します。これが、IIRフィルタ（やFIRフィルタ、ARモデルやMAモデル）を貫く基本概念になります。

本当なら1次の極だけで足りるのですが、それですとIIRフィルタの係数が複素数になってややこしいので、一般には実軸対称の位置に（複素共役の位置に）もう1つ点を取ります。すると、2次の実パラメータを持つIIRフィルタが構成できます。

導出

フィルタの極p1は、周波数fとサンプリング周期kから次のように決まります。

φ = 2πf/k 
p1 = cos(φ) + jsin(φ)

もうひとつの極p2は実軸対称の位置に取りますから

p2 = cos(φ) - jsin(φ)

です。極が２つあるときのIIRフィルタ伝達関数は

となり、ここから真面目に計算すると、

IIR(z)の分母 = 1-p1p2(z^-1)+(z^-2)
p1p2 = 2cos(φ)

が導かれます。これをIIRフィルタの定義式に戻すと、

になり、アルゴリズムの式が導かれます。

■参考文献

• ディジタルフィルタとｚ変換（要Java）

Z平面や極がどのような影響を与えるか、どのような意味があるのかを、目で見て、手で触って掴むことができます。FIR/IIRフィルタやARモデルなどを考えるうえで非常にためになるサイトです。

■注意

複数の波形合成を行ってその中から寄与率の低い波形の演算をやめるという手段を取ると、SONYの特許4019824に抵触する恐れがあります（「特許、登録」を選択して「[RN]登録番号」で検索してください）。まあ普通は触れようがないと思いますけど。

この特許、出願時（2002年）の資料を見るとsin波の生成だけで取ろうとした形跡がありますが、大昔から自明なこと（Z平面の単位円上に極を持つIIRフィルタは振動することと同値）なのでそりゃ無理です。