タケダノヲト-vol.18 ソフトウエアLFO

ソフトウエアでLFO

LFOは低い周波数の発振器。タケダノヲト的なコンテキストでは、シンセの音色を変調するためのCVを発生させるために使う。たとえば、ビブラートや、シーケンサーのテンポの設定とかに使うのだけど、マイコンを使って、コレを生成してみようというのが狙い。

まずは、グラフ用紙を取り出してみる。これに2点を定め直線で結ぶ。線引きで、ぴゃっと結ばずに、マス目を塗りつぶして線を描くようにする。で、あるXの時のYの値をCVとして出せばいいんじゃないか？一定の値を超えたら変化の符号を反転させてやれば、発振状態になるだろう？というところからスタート。
縦軸としてのCVの出力は、256レベルの変化を出すのなら256数える。256になったら、今度は0までカウントダウン。コレを繰り返せばいい。カウントアップ/ダウンするタイミングをタイマーの管理し、このタイミングを変えることで周波数を変えることができそうだ。
ただ、このアルゴリズムでは、貴重なリソースとしてのタイマーをLFOの出力だけの為に占有してしまうことになる。たとえば、すでに、タイマーはPWMで出力で使っていて、新たにタイミングを生成するタイマーが無いこともありうる。
いわば、発振周波数を変えるために、グラフ用紙の横軸を伸び縮みさせている状態ともいえる。たとえば、1枚のグラフ用紙の上に、複数のLFOや、EGの出力を載せようとすると無理が生じてくる。

ブレゼンハムのアルゴリズム

システム全体のクロックは共通とする。誰かの都合で横軸を伸ばしたり縮めたりしない。つまり、グラフ用紙は1枚。この1枚のグラフ用紙から、複数のLFOや、EGの出力を取り出すことができたら素敵じゃないか？というわけで、1枚のグラフ用紙の変わりに、コンピュータのディスプレイを当てはめてみると、コンピュータグラフィックのアルゴリズムの中でも基本中の基本であるブレゼンハムのアルゴリズムを使えば、縦横が伸び縮みしないディスプレイの上に、自由な線分を描画することができる。
これを応用してソフトウエアでLFOを書いてみた。

直線の描画の考え方

もともと、時間も空間もアナログ量。言い換えれば、連続系。これをコンピュータのマス目に区切られた画面や、タイマーで区切られる一定時間という離散系へ丸め込む、と考えてみる。最終的に、マス目や、タイマーのタイミングで丸められてしまうのだから整数の計算のみで処理したい。
当然、無理やり当てはめるのだから、誤差は出る。一番最初に紹介した、1LFO/1タイマーメソッドのほうが精度は高い。このメソッドでは、誤差は次のステップに引き継がれ、集積されていくので、ミクロ的には誤差はたんまりあるが、波形数個分にまで視野をひろげ、トータルでみれば、そこそこの精度がでるはず。ここでは、ビブラート用のLFOと用途を絞って考えてみる。
図は、XY平面にX1,Y1、X2,Y2の2点を結ぶ直線を描画している。この直線の傾きは、hh/ww。これは、Xが１進む時のYの増分はhh/wwという事。Xが１増えるごとに、この増分を累積して行き、これが1を超えたとき、Yは１つ進む。1より下は、小数点だから、グラフ用紙には反映されない。ようするに、Yの値は変わらない。Yの増分を累積して、これが１を超えるかどうかだけを見ればOKという形。
１を超えるかどうかの判断は、座標軸とはまったく無関係なので、増分にwwを掛けて約分してしまい、wwを超えるかどうかを判断の基準にしてしまう。１つXが変わるたびに必要な計算は、足し算１つ、ある数字との比較１回、必要な場合には、Yを変化させて、現在のYを出力。ループ1回アタリの演算の数もケチっている。
ポイントは、hhとくらべて、ww(分母)がでかいかどうか。でかい時には、Xが一つ進むごとに、Yが変化するかどうかを判断することができるけど、hhとくらべて、ww(分母)が小さいのなら、hhとwwを入れ替えて、ｙが一つ変化するとき、xが変化するかどうかを監視するパターンに切り替えなければいけない。hhとwwが一緒だったら、Xが一つ進むごとに必ず、Yは一つ変わることになる。
ここでは、hhとくらべて、wwはホボ一緒または、wwのほうがでかい場合(傾きが１以下の時）に限定できるように各条件をそろえることで、実現させてみる。

もともとコンピュータグラフィックスの基本アルゴリズムとして大変有名なものなので、検索すればもっと素敵な解説はたくさん見つけることができる。ここでは傾きの数を小数点の付いた数字にせず、傾きを作る分数のままで考え、分子に当たる数字を累積して、分母を超えるかどうかと言う判断に展開する。（というか、展開しないというか..）
これで計算はすべて、整数のみで行える。

空間の変化を時間の変化へと変換

直線の描画の考え方2

先の図は、X,Y平面だったけど、今度は時間Tと電圧Vの変化のグラフに書き換えてみる。といっても、軸の名前が変わっただけで、ホボ一緒に見えるが、一定時間変わったときに変化する電圧が表現されていると考える。
T1から、T2の間は、傾きの符号は+。T2で傾きの符号だけが変わる。T2からT3までは符号の向きが変わり、右下がりの直線になる。T3で符号が反転すれば、三角波になるはず。
最大値、最小値ぞれぞれの点で増分の累積はリセットされ、Tが１つ変わったとき、Vを１足すか引くかを切り替えて傾きの符号の符号を変える。いうまでもなく、T1からT3間での時間から周波数が読める。
Vは、出力される電圧。たとえば、8bit出力なら256段階、上がって降りて(T1～T3)で１周期だから、すべての段を1回ずつ触って回る(つまり最高速の1周期は)512ステップとなる。1ステップあたり１回CVが変化するチャンスが来ると言う事。
このアルゴリズムでは、決められたインターバルでこのLFO波形を出力する関数を呼び出し、そのタイミングでのCVを取り出す形なので、どういうインターバルでこの関数を呼び出すかで、出力する周波数が変わる。
ここでは、PWMの出力に使うタイマーのオーバーフローの割り込みを基準タイミングとしてみる。

CLK	8MHz	9.5MHz	20MHz
PWM fc	32us(31.25kHz)	26.95uS(37.11kHz)	12.8uS(78.13kHz)
LFO MAX fc	16.39mS(61.01Hz)	13.80mS(72.46Hz)	6.55mS(152.67Hz)

クロックごとの最高速。8bit幅PWM出力時

8bitのPWMを使うのなら、クロックを256回数えると１波形生成が1タイミングとなる。このタイミングごとに現在の傾きの値を計算して、次の段に進むかどうかをチェックするということだ。毎回、Vの値が更新される場合が、仮想的に描画している線分の傾きが1であり、これが最高速となる。
3種類のクロックで8BitのPWMのオーバーフローのタイミングで、波形更新する場合の最高速を一覧表にしてみた。
PWMの波形1サイクルが終わらないと、パルスワイズの更新が起きないので、コレより速いスピードにするには、別途DAコンバーターの増設が必要になる。このあたりの周波数では、凶悪なFM系モジュレーションというより、ビブラート用のLFOぐらいにしか使えないが、そこが狙いならツボともいえる。
この手の各種タイミングの計算や、検討には何時も使っている以下のスクリプトが便利。

TimerOverFlow Count Value

最大/最小周波数の検討

120bpm	時間	周波数
全音符	2S	0.5Hz
2分符	1S	1Hz
4分音符	500mS	2Hz
8分音符	250mS	4Hz
16分音符	125mS	8Hz
32分音符	62.5mS	16Hz

テンポ120のときの音符の長さと周波数

最高速は、ハードウエアの制約で決まってしまう部分がある。9.5MHzで8bit幅のPWMなら、最高速は72Hz程度でるけど、ビブラート専用として考えるなら、テンポ120bpmの32分音符分程度のスピードがあれば十分かもしれない。とりあえず、余裕を見て最高速は20Hz、最低速を0.3Hzとして検討してみる。
先に、最高の周波数を検討する。9.5MHzで走る8BitのPWM波形生成を最小タイミングで最高速で動かせば、72.46Hz。4回に１度LFOを呼び出せば、周波数は1/4、となり18.12Hz。ちょっとたりないかもしれないけど、最高速のデータを出しているときに、仮想的に描画してる線分の傾きを1にすることで、各種計算を楽にする方針なので、ここはぐっと我慢。
一方、最低速のほうは、比較的自由に決められる。たとえば、0.3Hzとするなら、スピードの変化率は、(18.12/0.3=)60.4。適当に整数化して62とする。(61とすべきだけど、ちと、計算違いがあって...ここではアルゴリズムを示すと言うことで数字には細かいことは言わないで頂戴!)。
最高速を１とすれば、最低速度は、最高速度の、1/62。仮想的に描画する直線の傾きが1/62ということになる。これ、分母と分子に、適当な数字をかけて、256に近づけておけば、スピードを変えるパラメーターのDレンジを広げられる。
ここでは、それぞれ4を掛けて、4/248とする。コレが一番のんびりした時。最高速は、1なんだけど、コレを248/248と考え、分子を、4から248まで変化させて、仮想的に描画する線分の傾き徐々にきつくして、最終的には45度で、最高速になると言う寸法。

ここまでのところをソースにまとめてみる。入力はスピードのパラメーター。4から248までが期待されてる。決まったタイミング、(ここでは、107.8uS毎）で呼び出されることが前提で、出力は問い合わせを受けたときのLFOのレベル。

void lfo(uint8_t speed)		//	speed 4 .. 248 -> 18.12Hz .. 0.3Hz
{
	static int current_e = (248-4）/2;	//増分の集積クリア。
	static int current_level = 128 ;	//最大振幅0から255の真ん中からスタート
	static int direction = 1;			//変化の向き

	current_e += speed;				//増分の集積
	if (current_e >= 248){			//分母を超えたかどうか
		current_e -= 248 ;			//はみ出した誤差分を次のステップに送る
		current_level += direction;
	}

	if (current_level > 255 ){		//レベルが最大値に達したときの処理
		direction = -1;				//傾きの反転
		current_e = (248-4）/2;		//増分の集積クリア。
		current_level -= 1;			//新しいレベルの設定
	}
	if (current_level == 0){		//レベルが最大小に達したときの処理
		direction = +1;				//傾きの反転
		current_e = (248-4）/2;		//増分の集積クリア
		current_level += 1;			//新しいレベルの設定
	}
	return current_level;
}

「増分の集積クリア」は0としないで、判断の値の半分からスタートすることで、仮想的に描画される線分の位置を変化の真ん中に合わせる事ができる。

発振周波数を変化させる

この関数に4から、248の数字を入れて、一定タイミングで何度も呼び出し、得た数字を並べると三角波の波形と同じ変化を取り出すことができる。スピードを決める数字は、結果として分子を変えているので、得られる周波数の変化はリニアな変化となる。
たとえば、ADコンバーターから取り出した0から、255までの数字を4から248までの数字にコンバートしていれれば、リニアに周波数の変化するLFOとなる。

	OUT = (IN X ((248-4)) / 256) + 4

ここでは指数的に周波数が変化するようにしたほうがLFOの周波数の設定としては使いやすい、というか、そうしないと使い物にならないと思ったほうがいい。で、もう一工夫してみる。

指数カーブを得たいのなら普通には y=exp(x)で得られる。 y=exp(x)のｙとして、4から248までが得られるｘを求めて、y=exp(x)に、このリニアな変化をいれて、4から248の指数的な変化をする値を取りだし、LFOの周波数をコントロールすればよさげ。
4と、248、それぞれ、指数関数の逆関数としての対数関数に突っ込んでみる。エクセルではlog(4)=0.60206、log(248)=2.3944が得られる。これは、10を底としたときの数字で。指数関数と対数関数の、計算の基準と成る底が違ったらそりゃ結果も違う。エクセルで数字を得るには、指数関数のかわりに、pwoer(n, x)を、対数関数は底を指定して、log(x, n) (それぞれnが底）として計算する。底は10でも2.718でも、もちろん2でも、このカーブの描く形は変わらない。

直線の描画の考え方2

底を2にしたときの4と248の対数は、2と7.9541。先ほどは0から、255までの数字を2から248までにコンバートしたけど、0から、255までの数字を2から7.9541までの数字にコンバートし、

	OUT = (IN X ((7.9541-2)) / 256) + 2

これを2を底とする指数関数（y=exp(x,2))の入力とする。
(4と248のの対数である2や、7.9541は、プログラムの中に定数として埋め込まれる。ライブで計算する必要はなくエクセルなどで1度計算してしまえばOK)
さて、exp(x)がちょっと難関。2と7.9541の間の間にある、256個の小数点つきの数字の指数を計算するのだけど、指数関数の性格として、整数部分と、小数点以下の部分に分割して計算することができる。たとえばこんな感じ。


	log(2.5) = log(2) x Log(0.5)

この法則を使って計算の比較的簡単な整数部分と、面倒げな小数点以下の部分とに分けて計算する。前半、整数部分は普通に掛け算で得られる。しかも、2を底として計算すれば、さらにお得なボーナスがついて計算はシフトで済ませられる。
後半の小数点乗部分を計算してグラフ化してみたのが右の表。ほんの少したわんではいるけど、ほとんど直線ジャン？！そのまま行ってもバレねえよという発想が、小数点以下の数字の直線近似法。
あれ、ちとまてよ、小数点の計算がでてきちゃった。

固定小数点による演算

小さなCPUには高コストな小数点の演算をさけるつもりだったのに、たっぷりでてきちゃった。浮動小数点の計算は、まじに高コストなので、ここでは、固定小数点という演算方法で必要な精度に限定した計算をやってみる。
理屈は簡単。

A X B = C

だったら、

A X B X　D = C'

としてC'を計算して後から、Dで割りゃ良いんじゃん？ということ。Dを2の乗数にすれば、最後の割り算はシフトだけで済ませる。で、Dをどんな値すれば良いかの検討がミソになる。

この計算の結果最終的に得たい数字は4から248なので、結果は8ビットで取り出せればOK。
いまここでターゲットにしている数字は、2から、7.9541までの値。普通に引き算すると5.9541。これを、今求められている精度、8ビットで表現した時、最終的に切り捨てられる小数点以下の数字に、なんビット必要になるかを検討する。
5.954の幅を、8ビット幅で表現するのだから、255で割って1ビットの重みとして0.023を得る。１は、(1/0.023=)43.47。これが、1以下で必要な解像度。とりあえず6ビット(64段階)ないと必要な解像度が得られない。ここではマジックナンバー、Dを64とする。
この考察から、先ほどでた、小数点以下の面倒な部分の計算を直線で近似せず、あらかじめ計算しておいた値を使う場合のテーブルサイズが、64個になると考えることができる。
というわけで、小数点を6ビット目にするのなら、最小単位は、(1/64=)0.0156となり、この系で7.9541を変換すると(7.9541/0.0156=)509.87。同様に最小値は、2は(2/0.0156=)128.20。それぞれざっと正規化して0から、255の値Nを、128から510に変換する。

赤が、きちんと計算したカーブで、黄色が直線で近似のカーブ。

(N X (510-128) / 256) + 128

これで得た数字の6ビット目以上は整数部分。以下は、小数点以下部分となる。
さて、このアルゴリズムをプログラミング言語に実装することを考えると、perlなどの言語では心配しなくて済むが、C言語の場合、データのサイズ(扱える数の最大値）は、自分で設定しなければいけない。たとえば、8ビット同士の足し算なら結果を収めるデータ領域は、9bit、掛け算なら16ビット必要になる。
計算をローコストに抑えるために、可能な限り8bitで済ませたいが、どうしても8bitで足りないときは16bitで押さえるように検討するべき。ここでは、(510-128=) 382に入力の最大値255を掛けると、16ビットを超えてしまう。しかも、この計算のあと、また8Bitの数で割るのだから、トータルで16bitでは収まるはずの数だったりすることを考えると、ここはできるならケチりたい。さらに、計算の中に出てくる割り算は、ホンチャンの割り算は使わずにシフトを使って高速に済ませたいので、式を展開した。このとき、掛け算の桁あふれは怖いが、割り算を一番最後にしないと精度を落とすことになるので要注意。ここでは、定数の掛け算と、定数の割り算が入るので、約分して定数を16bit以内に収めて精度も確保した。


(N X ((510-128)/2) / (256/2)) + 128

すでにマジックナンバーDが掛っている状態で、0から255までのリニアな変化を、2から、7.9541まで変化する値を取り出せた。最後に、マジックナンバーで割るのを忘れないようにしなくっちゃ...
この式は、まま、C言語のソースに突っ込めるのだけど、N以外の部分はずべて定数になるので、正規化されて、プログラムの中での演算は、単純化されている。あとから、見直したとき、可能な限り式は原型のまま残しておいたほうが、あとからの理解/応用に便利だ。
さて、指数の計算。まず、ターゲットの数字を6ビット右にシフトして整数部分で、この指数計算の底である2を左にシフトして、整数部分を得る。
下から6ビット分が小数点以下。本来なら、微妙にたわんだカーブになるところを直線で近似しちゃうわけだけど、単純に、1から2まで変化する値に変換するには１(ここでは、固定小数点化した1ということで64）足せばOK。
こうして得た、整数部分の指数と、小数部分の指数をかければ、求める値の指数を得ることができる。


unsigned char calc_LFO_speed(unsigned char value)
{
	unsigned int temp;		/*16bit */
	unsigned int speed;		/*16bit */
	unsigned char exp;

//	(N X ((510-128)/2) / （256/2）) + 128 プラス１は切り上げ。
	temp = value * ((510-128)/2+1);
	temp >>= 7;
	temp += 128;
	speed = 2 << ( (temp >> 6 )-1) ;				// 整数部分計算	 pow(2, temp/64 ) 
										// 乗数が0の時の処理がない。要注意  
	exp = ((unsigned int) temp & 0x003F) + 64 ;  		// 小数点以下は直線近似
	speed = (speed * (unsigned int)exp );			// 整数部分と小数点以下を合成
	speed >>= 6;							// 小数点以下切捨てて整数化
	return speed;
}

実装と評価

ATTiny13にソフトウエアLFOを実装してみた。

結果論だけど、このアルゴリズムでのLFOの実装は、途中でスピードが変化した場合もそれなりに動作してくれる。スピードのパラメーターをADCから取り出す場合、精度的な問題で取り出した値にワウフラッタ(にも似たデータのブレ）があっても問題なく追従する。LFOの動作の途中に関係なく自由にスピードを換え、変わったところでそのスピードが発振に反映される。
逆に言えば、デジタルのクセに、アナログライクなADコンバータの誤差を含んだ出力が出てくる。デジタルだから絶対正確というわけではない。回路はもうこれ以上ないぐらいシンプル。組み立てる楽しみは殆どない、面倒くさいだけ。

Tiny13LFO-TypeA.zip

AVRstudioのプロジェクトファイルで2種類のソフトが同梱されている。一つは可能な限りこのページで紹介したアルゴリズムに忠実に再現した版。
もう一つは、このアルゴリズムの利点を生かし、タイマー一つで、複数の周波数を出力できるようにした。
LFOが2ch分。三角波が出るほか、片方は、スイッチで矩形波にも出る。
このままで、2.5Vのオフセットを持っているので、オペアンプなどのバッファを介してセンターを0Vへシフトしたほうがいいかもしれない。シンプルなLFOならば、ファンクションジェネレーターの回路を使えば、このLFOにはかなわないまでも、8pinDIPのデュアルオペアンプ一つ他、数点の部品で組み立てることができる。
このままでは、

2ch、違うスピードで取り出せる
スピードの設定が電圧で可変できる

の２点しかメリットがないと言う微妙な仕様。
1ch出力として、さまざまな波形が取り出せるぐらいの機能を持たせないと、アナログチップの強力なリプレイスとしてはあまり意味がないかもしれない。むしろ、別のソフトウエアの中で、必要になるLFO的に変化するデータとして使うほうがさまざまなメリットが生まれそうだ。
本ソフトは、あくまで上記アルゴリズムの紹介の一環として、実用性はどなのよバージョン。

指数関数的に周波数が変わるようにしたのは、一定角度ごとに周波数が倍(1オクターブの変化)が出るようにしたかったから。本ソフトでは、0.3Hzから18Hz強の周波数帯域を持つ訳だけど、0.3Hzの2倍が0.6Hz(これで1オクターブ)さらに倍が1.2Hzと倍々にしていくと6オクターブ弱ある。300度の回転角を持つボリュームなら、50度ごとにオクターブアップしていくはず。言い換えれば、全音符、二分音譜、四分音符と、一定角度毎変化する。

このページは、「阿部さんの掲示板」こと、アナログ震世界のアナログシンセ掲示板に投稿させていただいた記事「試しながら学ぶAVR入門」というスレッド内の、「ローコストLFO実装」を中心に加筆・修正してまとめなおしました。

ソフトウエアでLFOを実現する