前回、BLITを使って鋸波をうまく生成できることが確認できましたので、今度はlogue SDKで実装してみます。
今回のコードおよびバイナリは以下のリポジトリに置いてあります。
実装にあたって、追加で考慮が必要な点が２つあります。
（１）\(Sinc_M\)関数で分母が0の場合
（２）Leaky Integrator
まず、（１）\(Sinc_M\)関数で分母が0の場合です。
\(Sinc_M\)関数は以下のような形をしています。
$$y(n) = \frac{\sin(M\pi \phi)}{P\sin(\pi\phi)} $$
この関数は、\(\phi=0, 1, 2, \ldots\)のとき分母が0になります。ですので上の式を普通に計算すると、ゼロによる除算でエラーになる可能性があります。
分母が0の時の\(Sinc_M\)関数の値は、ロピタルの定理を使って
$$ \lim_{\phi\rightarrow 0} \frac{\sin(M\pi \phi)}{P\sin(\pi\phi)}
= \lim_{\phi\rightarrow 0} \frac{\sin'(M\pi \phi)}{P\sin'(\pi\phi)} \\
= \lim_{\phi\rightarrow 0} \frac{M\cos(M\pi \phi)}{P\cos(\pi\phi)} \\
= \frac{M}{P} $$
のように求められます。
実際には、この値を分母が0の時だけでなく、0に近いごく小さい値（今回の実装では1.0-E5としています。サンプリング間隔が2.08E-5なので、その半分程度です）で用いました。おそらく計算誤差のせいだと思いますが、試作中に音を鳴らしていたところ、この閾値が小さすぎると、数十秒に一回程度、\(Sinc_M\)関数の値がピークになるべきところで、本来よりずっと小さい値になってしまうことがあったためです。
BLITでは、直前のサンプルとの差分をずっと加算しつづけることで信号を生成しますので、その差分の数値に誤差があると、それはずっと信号に残ってしまいます。つまり、信号の中心線（本来は0）が、だんだんプラス側やマイナス側に寄って行ってしまう可能性があります。
これは誤差以外にも、発音中にピッチが変化した場合に起こります。本来ピッチが変化した場合は\(y(n)\)や\(\phi\)をリセットすべきですが、そうするとピッチベンダやポルタメントのときにプチプチノイズが出てしまいます。
（２）Leaky Integratorは、入力の誤差がずっと蓄積されないよう除去するために「もし差分が0の状態が続いたら」\(y(n)\)の値がゼロに近づいていくようにする処理です。
イメージとしては、ちょっと穴が開いている入れ物に水を少しずつ入れていくのだけど、穴からも水が漏れている、みたいな感じです。
通常の積分器は
$$ y(n+1) = y(n) + delta $$
という計算をしますが、Leaky Integratorは
$$ y(n+1) = (1 – k) y(n) + delta \\
(k \ll 1)$$
という計算をします。常に\(ky(n)\)だけ漏れていくので、\(delta=0\)なら\(y(n)\)はだんだんゼロに近づいていきますね。
\(k\)が大きすぎると、鋸波の波形がたわんできますので、0.1%以下くらいにしておくのが良さそうです。
考慮すべき点は以上ですが、そのほか
・\(\phi\)の初期値は0.5（\(Sinc_M\)の振幅が最小のところ）にすると、積分値yの初期値も0としてよい。
・\(\phi\)の値域は0.0~1.0ではなく、0.0～2.0とする。（\(\sin\)のパラメータが0～\(2\pi\)のため）
といったあたりが注意点になります。
下の図は、実装してみたオシレータでC1（32.7Hz）を鳴らしてみたときの周波数分布です。
高いほうの倍音成分まできちんと出ています。もしこれをウエーブテーブルでやると、１波長分で1468サンプルが必要になります。
出音はこんな感じです。（Prologueの出力を録音してMP3で圧縮してあります）
比較のために、マルチエンジンに入っているSawで鳴らしたものはこちらです。これはおそらくサイズの小さいウエーブテーブル（logue SDKで使えるものと同じもの）を使っていると思います。
音量が若干違いますが、それは別にしても、倍音成分が少ないため、かなりダークな音色です。
ちなみに、以前作成したウエーブテーブルオシレータで鳴らすとこうなります。
テーブルサイズは256サンプル/波長ですが、倍音成分はそれなりにあります。しかし、冒頭の低音部で4KHzあたりに僅かにノイズが出てしまっています。これはテーブルサイズが小さすぎて、倍音を128次（=256/2）で打ち切っているため、残留成分（128×32.7Hz＝4185.6Hz）が出てしまっているのだと思います。
実装そのものは短いので、上のリポジトリのコードを見ていただければと思います。
三角関数の計算は普通のsinf()関数を使っていますが、logue SDKのターゲットCPUのSTM32F401は単精度FPUを搭載していますので、１サンプル２回程度の計算なら処理能力に特に問題は無いようです。
FPUが無いマイコンでBLITを使用する場合は、Sinの高速化について工夫が必要になると思います。この計算は精度が重要ですので、logue SDKのosc_sinf()（長さ256のテーブル）ではもちろん実用になりませんが、fastersinfullf()はまあまあ悪くありませんでした。
なお、上記のリポジトリの実装では、おまけとして、
・Leaky Integratorのパラメータ（上の式での\(k\)）
・信号の周波数の上限（\(f_s / 2\)よりも小さくできる）
・倍音の次数（基本周波数のみ～100次の倍音まで設定可能）
・\(Sinc_M\)関数を出力（試験やデバッグ用）
の４つのパラメータを用意しました。倍音の次数をリアルタイムで変更するのは、フーリエ級数を使ったら重すぎて実現できなかったと思います。
BLITについては以上です。
背後の理論を知らないと、ソースコードを見ても、訳が分からない計算をしているようにしか見えないと思います。
オブジェクトコードのサイズが５キロバイト程度と小さいのに、大きなウエーブテーブルを使った場合と同等以上の音質が得られます。
logue SDKのような、メモリ空間の制約が厳しいけれど計算能力はそこそこ高い環境にはピッタリだと思います。
一方で、実際に音を出してみた感想としては、高域についてはもう少し工夫が必要になりそうです。高域の音については、\(Sinc_M\)の分子部分が周波数の小さな変動に対して大きく反応することがあるからです。
\(Sinc_M\)関数の定義から以下のようなことが分かります。
・生成する周波数が高いほど、\(P\)が小さくなるので\(Sinc_M\)のピークが高くなる
・\(Sinc_M\)の分子は、常にサンプリング周波数の1/2以上であり、比較的高い周波数である
また、積分と言っても離散的に計算して\(\delta t\)を掛ける長方形近似ですから、積分の値にはもともと誤差があります。
そして音程が高くなると、サンプルの数が少なくなります。例えばC7（ノート番号96）は4,186Hzで、サンプリングレートが48KHzなら11.47サンプルしかありません。
するとどうなるかというと、結果だけ書くと、周波数が高くなるほど、積分値の誤差が増えます。
\(M\)の値は離散的に変化しますが、その変化の前後において波形が滑らかには変化しませんし、\(Sinc_M\)の分子がナイキスト周波数に近づくと、うなりのような振幅も出てしまいます。そのため、ピッチベンドやポルタメントなど、連続的に周波数を変化させたときに、振幅の変化に伴うプチノイズが出てきます。このノイズは周波数が高くなるほど、（\(P\)が小さくなるので）大きくなります。
例えば下図は、高域の音程でLFOをかけて周波数をわずかに、かつ低速で変動させた場合の波形です。時間軸を圧縮しているので、波形そのものは見えなくなっていますが、ピーク値が周期的に不連続なジャンプをしていることが見て取れます。
対策としては、オーバーサンプリングで波形を求めて（例えばサンプリング周波数の2倍の周波数で計算）からローパスフィルタを通す方法が考えられますが、計算量が増えてしまうのが悩ましいところです。あるいは、近似的には長方形近似の誤差を補正するように、ナイキスト周波数に近いところの倍音は捨てる（\(M\)の値を切り下げる）か、アドホックですが補正用の係数（0.5~1.0）を掛けてピーク値を調整するということも考えられるかもしれません。
しかし原理的には、高域の音程においては波形そのもののサンプリングデータを持っているウエーブテーブル方式の方が有利になると思われます。
コメント