シミュレーションによるCLAPプラグインの作成
普段は大学で物理学徒をしている私が、シミュレーションを用いてCLAPプラグインを作ろうとして失敗した話です。
Repository
構想
少し前に私はAerophoneという楽器を頂く機会がありました。
私はDAWをメインに曲を作っているのですが、どうにもAerophoneで演奏した音源をDAWに取り込む方式だと、Interfaceなどの関係で音の質があまりよろしくないということに気が付きました。
そこで、楽器の発音と共鳴の機構をシミュレートしたプラグインを作成して、音源よりもリアルな音を出す電子的な楽器を作ってしまおうというのが今回の趣旨です。
時代の後退ですし、録音をMidiで流したほうがずっとうまくいくのは理解していますが、SoundFont(sf2など)のようなものにはない表現力を持ってた方が面白いかなと思ったのが始まりです。
設計
設計は大きく分けて次のふたつの部分をつなぎ合わせることで実現しています。
- 発信部(サックスでいうリードの部分)
- 共鳴部(サックスでいうボディの部分)
個人的にはチェット・ベイカーのような優しいトランペットの音を目指したいのですが、サックスだろうがトランペットだろうがの発振は圧力駆動型バルブというモデルで説明できます。
フルートやリコーダーは空気の渦を作って発振するので全然違うのですが、流体シミュレーションはまだ理解しきれてないので、今回はちょっと勘弁してもらって。
実装理論(TL;DR)
ここから長丁場です。実際、結構な時間を費やしています。
発信部
まず、圧力駆動型バルブの運動方程式を示します。
変位はxで、kは弾性、P(t)は息による空気圧、rは振動の減衰のための定数です。
$$ m \frac{d^2x}{dt^2} = -k x + f(P(t),x) + r \frac{dx}{dt} $$
rの項の妥当性について一応述べておくと、例えば吹くのを止めて弾性によってのみリードが振動している場合、空気抵抗や構造への負担から振動が減衰していくことは経験的に理解できるでしょう。
この方程式を解くうえで厄介なところはP(t)が目まぐるしく変わってしまうことに加え、fがx依存性を持っているところです。
もしリードが息に押されてマウスピースに近づくと、その近さの2乗に応じて空気の通る速度は早くなります。(ロケットエンジンが吹き出し口からスカート型に広がっている理由なのですが。)
また、リードがマウスピースに触れた瞬間に空気の流れは0になり、十分に弱い圧力だけがリードにかかることになります。するとリードの弾性が優位になってリードはマウスピースから離れていきます。
今回はこの空気のモデルを扱うために、f(P(t),x)についてもう少し考えていきます。(k,m,rもリードの噛む強さによって変化する(間接的にtに依存する)が、方程式を解くときは定数として扱います)
息を吐く強さを$P(t)$とすると、$\rho$を空気の密度、$g(t)$を開口幅、$v_f(t)$を空気の速度とすると、
$$P(t) = \rho L \frac{dv_f(t)}{dt} + B(t) v_f(t)^2$$
ただし、$B(t) = \frac{\rho}{4 g(t)^2}$
と表すことができます。これを双一次変換すると
$$\frac{dv_f(t)}{dt} \approx \frac{2}{T} (v_f[n] - v_f[n-1]) - \left.\frac{dv_f(t)}{dt}\right|_{n-1}$$
から
$$\frac{dv_f(t)}{dt} \approx \frac{2}{T} (v_f[n] - v_f[n-1]) - \frac{1}{\rho L} \left( P[n-1] - B[n-1] v_f[n-1]^2 \right)$$
が求まり、
$$P[n] = $\rho L \left[ \frac{2}{T} (v_f[n] - v_f[n-1]) - \frac{1}{\rho L} \left( P[n-1] - B[n-1] v_f[n-1]^2 \right) \right] + B[n] v_f[n]^2$$
これを整理すると
$$B[n] v_f[n]^2 + \left( \frac{2\rho L}{T} \right) v_f[n] - \left[ P[n] + P[n-1] + \frac{2\rho L}{T} v_f[n-1] - B[n-1] v_f[n-1]^2 \right] = 0$$
これを$v_f(t)$についてとくと
$$v_f[n] = \frac{-A + \sqrt{A^2 + 4 B[n] C[n-1]}}{2 B[n]}$$
where,
$$A = \left( \frac{2\rho L}{T} \right)$$
$$C[n-1] = \left[ P[n] + P[n-1] + \frac{2\rho L}{T} v_f[n-1] - B[n-1] v_f[n-1]^2 \right]$$
となります。
ところでリードに空気からかかる力を
$$f[n] = \pm \frac{1}{2} \rho v_f[n]^2 g[n]$$
($\pm$が正のときリード(サクソフォン等)のモデル、負のときリップリード(トランペット等))
と考えます。
最初に示した
$$m \frac{d^2 x(t)}{dt^2} + r \frac{dx(t)}{dt} + k x(t) = f(t)$$
を$T$をサンプリング周期として一次双変換すると
$$x[n] = \frac{b_0 f[n] + b_1 f[n-1] + b_2 f[n-2] - a_1 x[n-1] - a_2 x[n-2]}{a_0}$$
where,
$b_0 = T^2, \quad b_1 = 2T^2, \quad b_2 = T^2$$a_0 = 4m + 2rT + kT^2$$a_1 = -8m + 2kT^2$$a_2 = 4m - 2rT + kT^2$
となります(一次双変換の過程や実際に運動方程式の解の近似になっていることはリポジトリのReadmeのTL;DR Derivation of the simulation formulaセクションで示しています)。
これで、$x[n]$がnステップ目の位置を表す関数として記述することができました。
共鳴部
筒の剛性や、空気の密度の変化から気柱による音波の反射をシミュレートすることも考えましたが、
音波の距離によって遅れる反射は認めるものとして、反射によるエネルギー減衰のモデルを考えることで倍音の高音成分がなくなるシミュレートを実装します。
エネルギーは位相のt微分の2乗に比例する量であり、減衰はxのt微分の2乗に対して行います。
エネルギーが減衰すると、高音成分から現象していくことは理解できるかと思います。
式としては、発信部から来た変位をx_n、遅延して送られてくる変位をx_resonance、前回の合計の変位をx_prevとすると、減衰定数aを使って
$$ x = a (x_xrev - (x_n + x_resonance))^2 $$
と書くことができます。ここで、x_resonanceは筒を往復するのにかかった時間だけ前の変位に、開管なら1を、閉管なら-1をかけたものです。
x_resonanceが管を往復するのにかかった時間だけ前の変位に反射の係数をかけたもので十分な理由を示します。
音波の伝搬について、
$$\frac{\partial^2 p}{\partial t^2} - c^2 \frac{\partial^2 p}{\partial x^2} = 0$$
この一般解はダランベールの解を用いて
$$p(x, t) = f(t - x/c) + g(t + x/c)$$
と書くことができ、波の振幅の移動として書くことができます。
また、反射が起きるx=Lの点について書くと
- 開放端のとき
筒の端が大気圧になっているので$p=0$という条件を入れて
$$f(t - L/c) + g(t + L/c) = 0$$
これを $g$ について解くと
$$g(t + L/c) = -f(t - L/c)$$
- 閉口端のとき
筒の端で空気の速度が0という条件を入れて
$$\frac{\partial p}{\partial x} = -\frac{1}{c} f'(t - x/c) + \frac{1}{c} g'(t + x/c) = 0$$
$x=L$を代入して
$$g'(t + L/c) = f'(t - L/c)$$
$$g(t + L/c) = f(t - L/c)$$
を得ることができ、往復するのにかかった時間だけ前の変位に反射の係数をかければ良いということがわかります。
補足ですが、楽器の閉管は片方が開放端でもう片方が閉口端、開管は両方が開放端なものを指します。
なので、閉管では1往復すると位相が逆->そのままという二回の反射で位相が逆に、開管では1往復すると位相が逆->逆という二回の反射で位相が戻ります。なので上のような係数になっています。
元も子もない話をしてしまうと、ローパスフィルターをかけて出力すればそれで終わりなのですが、それだとあまりにも面白くないと思うので、このような回りくどい実装になっています。
実装において重要な往復にかかる時間の導出を行っておきます。
入力されたmidiノートの周波数を$f$、音速を$c$とすると、波長は
$$\lambda = \frac{c}{f}$$
管の長さLは$\lambda$を用いて
$$\lambda = 2L$$
とかける(閉管の場合は4Lだが、同様の変形を行うだけなので、以下では開管の場合のみを扱う)
また、Lを往復するのにかかる時間は
$$\frac{2 L}{c} = \frac{2}{c} \frac{c}{2f} = \frac{1}{f}$$
とわかるので、Lやcを実際に使う必要はなく、Lを往復するのにかかる時間が求められます。
実装
実装はソースコードにあるものが全てです。うまく動かないと理論を見直したり、プラグイン部分の不具合を排除するために理論式をipynbでプロットしながら検証をしたりしました。
コーディングエージェントについて
コーディングAgentはうまく動かない場面が多かったです。pythonのプログラムやLaTeX式をRustの関数にする部分で利用しましたが、結局のところ、私が書き換える羽目になりました。
よく知られている実装以外をするときは向かないですね。その代わりに向いている部分もありました。
パラメーターフィティングです。運動方程式のmやrやkは少し変わるとうまく発振しないので、パラメーターフィッティングが難航していましたが、AIの背景にある理論自体がパラメーターフィッティングであるので、その部分を大いにうまく利用できたのかなと思っています。
現状について
設計の全ての部分の実装はできました。
しかし、目指している音を鳴らすことはできませんでした。
反響の共鳴に関わらずリードの固有振動数らしき周波数で鳴り続ける挙動をしています。
パラメーターが多すぎるのか、モデルが正確ではないのか。わからないところです。
より謎なのは、ipynbでうまく動作していても実際にプラグインに記述するとうまく動かないところです。
何はともあれ、私のシミュレーターを作る力が不足していたということでしょう。