creating_a_clap_plugin_through_simulation.md (10764B)
1 --- 2 title: シミュレーションによるCLAPプラグインの作成 3 tags: 4 - nih-plug 5 private: false 6 updated_at: '2026-07-01T20:48:24+09:00' 7 id: 925a92ee93e20b6e91ba 8 organization_url_name: null 9 slide: false 10 ignorePublish: false 11 --- 12 # シミュレーションによるCLAPプラグインの作成 13 14 普段は大学で物理学徒をしている私が、シミュレーションを用いてCLAPプラグインを作ろうとする冒険譚です。 15 16 # 構想 17 18 少し前に私はAerophoneという楽器を頂く機会がありました。 19 私はDAWをメインに曲を作っているのですが、どうにもAerophoneで演奏した音源をDAWに取り込む方式だと、Interfaceなどの関係で音の質があまりよろしくないということに気が付きました。 20 そこで、楽器の発音と共鳴の機構をシミュレートしたプラグインを作成して、音源よりもリアルな音を出す電子的な楽器を作ってしまおうというのが今回の趣旨です。 21 時代の後退ですし、録音をMidiで流したほうがずっとうまくいくのは理解していますが、SoundFont(sf2など)のようなものにはない表現力を持ってた方が面白いかなと思ったのが始まりです。 22 23 # 設計 24 設計は大きく分けて次のふたつの部分をつなぎ合わせることで実現しています。 25 26 - 発信部(サックスでいうリードの部分) 27 - 共鳴部(サックスでいうボディの部分) 28 29 個人的にはチェット・ベイカーのような優しいトランペットの音を目指したいのですが、サックスだろうがトランペットだろうがの発振は圧力駆動型バルブというモデルで説明できます。 30 フルートやリコーダーは空気の渦を作って発振するので全然違うのですが、流体シミュレーションはまだ理解しきれてないので、今回はちょっと勘弁してもらって。 31 32 # 実装理論(TL;DR) 33 ここから長丁場です。実際、結構な時間を費やしています。 34 35 ## 発信部 36 まず、圧力駆動型バルブの運動方程式を示します。 37 変位はxで、kは弾性、P(t)は息による空気圧、rは振動の減衰のための定数です。 38 $$ m \frac{d^2x}{dt^2} = -k x + f(P(t),x) + r \frac{dx}{dt} $$ 39 rの項の妥当性について一応述べておくと、例えば吹くのを止めて弾性によってのみリードが振動している場合、空気抵抗や構造への負担から振動が減衰していくことは経験的に理解できるでしょう。 40 この方程式を解くうえで厄介なところはP(t)が目まぐるしく変わってしまうことに加え、fがx依存性を持っているところです。 41 もしリードが息に押されてマウスピースに近づくと、その近さの2乗に応じて空気の通る速度は早くなります。(ロケットエンジンが吹き出し口からスカート型に広がっている理由なのですが。) 42 また、リードがマウスピースに触れた瞬間に空気の流れは0になり、十分に弱い圧力だけがリードにかかることになります。するとリードの弾性が優位になってリードはマウスピースから離れていきます。 43 今回はこの空気のモデルを扱うために、f(P(t),x)についてもう少し考えていきます。(k,m,rもリードの噛む強さによって変化する(間接的にtに依存する)が、方程式を解くときは定数として扱います) 44 45 息を吐く強さを$P(t)$とすると、$\rho$を空気の密度、$g(t)$を開口幅、$v_f(t)$を空気の速度とすると、 46 47 $$P(t) = \rho L \frac{dv_f(t)}{dt} + B(t) v_f(t)^2$$ 48 49 ただし、$B(t) = \frac{\rho}{4 g(t)^2}$ 50 51 と表すことができます。これを双一次変換すると 52 53 $$\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}$$ 54 55 から 56 57 $$\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)$$ 58 59 が求まり、 60 61 $$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$$ 62 63 これを整理すると 64 65 $$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$$ 66 67 これを$v_f(t)$についてとくと 68 69 $$v_f[n] = \frac{-A + \sqrt{A^2 + 4 B[n] C[n-1]}}{2 B[n]}$$ 70 71 where, 72 $$A = \left( \frac{2\rho L}{T} \right)$$ 73 $$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]$$ 74 75 となります。 76 77 ところでリードに空気からかかる力を 78 79 $$f[n] = \pm \frac{1}{2} \rho v_f[n]^2 g[n]$$ 80 ($\pm$が正のときリード(サクソフォン等)のモデル、負のときリップリード(トランペット等)) 81 82 と考えます。 83 84 最初に示した 85 86 $$m \frac{d^2 x(t)}{dt^2} + r \frac{dx(t)}{dt} + k x(t) = f(t)$$ 87 88 を$T$をサンプリング周期として一次双変換すると 89 90 $$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}$$ 91 92 where, 93 $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$ 94 95 となります(一次双変換の過程や実際に運動方程式の解の近似になっていることはリポジトリのReadmeのTL;DR Derivation of the simulation formulaセクションで示しています)。 96 97 これで、$x[n]$がnステップ目の位置を表す関数として記述することができました。 98 99 ## 共鳴部 100 101 筒の剛性や、空気の密度の変化から気柱による音波の反射をシミュレートすることも考えましたが、 102 音波の距離によって遅れる反射は認めるものとして、反射によるエネルギー減衰のモデルを考えることで倍音の高音成分がなくなるシミュレートを実装します。 103 エネルギーは位相のt微分の2乗に比例する量であり、減衰はxのt微分の2乗に対して行います。 104 エネルギーが減衰すると、高音成分から現象していくことは理解できるかと思います。 105 式としては、発信部から来た変位をx_n、遅延して送られてくる変位をx_resonance、前回の合計の変位をx_prevとすると、減衰定数aを使って 106 107 $$ x = a (x_xrev - (x_n + x_resonance))^2 $$ 108 109 と書くことができます。ここで、x_resonanceは筒を往復するのにかかった時間だけ前の変位に、開管なら1を、閉管なら-1をかけたものです。 110 111 x_resonanceが管を往復するのにかかった時間だけ前の変位に反射の係数をかけたもので十分な理由を示します。 112 113 音波の伝搬について、 114 115 $$\frac{\partial^2 p}{\partial t^2} - c^2 \frac{\partial^2 p}{\partial x^2} = 0$$ 116 117 この一般解はダランベールの解を用いて 118 119 $$p(x, t) = f(t - x/c) + g(t + x/c)$$ 120 121 と書くことができ、波の振幅の移動として書くことができます。 122 123 また、反射が起きるx=Lの点について書くと 124 125 1. 開放端のとき 126 127 筒の端が大気圧になっているので$p=0$という条件を入れて 128 129 $$f(t - L/c) + g(t + L/c) = 0$$ 130 131 これを $g$ について解くと 132 133 $$g(t + L/c) = -f(t - L/c)$$ 134 135 2. 閉口端のとき 136 137 筒の端で空気の速度が0という条件を入れて 138 139 $$\frac{\partial p}{\partial x} = -\frac{1}{c} f'(t - x/c) + \frac{1}{c} g'(t + x/c) = 0$$ 140 141 $x=L$を代入して 142 143 $$g'(t + L/c) = f'(t - L/c)$$ 144 145 $$g(t + L/c) = f(t - L/c)$$ 146 147 を得ることができ、往復するのにかかった時間だけ前の変位に反射の係数をかければ良いということがわかります。 148 149 補足ですが、楽器の閉管は片方が開放端でもう片方が閉口端、開管は両方が開放端なものを指します。 150 なので、閉管では1往復すると位相が逆->そのままという二回の反射で位相が逆に、開管では1往復すると位相が逆->逆という二回の反射で位相が戻ります。なので上のような係数になっています。 151 152 元も子もない話をしてしまうと、ローパスフィルターをかけて出力すればそれで終わりなのですが、それだとあまりにも面白くないと思うので、このような回りくどい実装になっています。 153 154 実装において重要な往復にかかる時間の導出を行っておきます。 155 156 入力されたmidiノートの周波数を$f$、音速を$c$とすると、波長は 157 158 $$\lambda = \frac{c}{f}$$ 159 160 管の長さLは$\lambda$を用いて 161 162 $$\lambda = 2L$$ 163 164 とかける(閉管の場合は4Lだが、同様の変形を行うだけなので、以下では開管の場合のみを扱う) 165 166 また、Lを往復するのにかかる時間は 167 168 $$\frac{2 L}{c} = \frac{2}{c} \frac{c}{2f} = \frac{1}{f}$$ 169 170 とわかるので、Lやcを実際に使う必要はなく、Lを往復するのにかかる時間が求められます。 171 172 # 実装 173 174 実装はソースコードにあるものが全てです。うまく動かないと理論を見直したり、プラグイン部分の不具合を排除するために理論式をipynbでプロットしながら検証をしたりしました。 175 176 ## コーディングエージェントについて 177 コーディングAgentはうまく動かない場面が多かったです。pythonのプログラムやLaTeX式をRustの関数にする部分で利用しましたが、結局のところ、私が書き換える羽目になりました。 178 よく知られている実装以外をするときは向かないですね。その代わりに向いている部分もありました。 179 パラメーターフィティングです。運動方程式のmやrやkは少し変わるとうまく発振しないので、パラメーターフィッティングが難航していましたが、AIの背景にある理論自体がパラメーターフィッティングであるので、その部分を大いにうまく利用できたのかなと思っています。 180 181 # 現状について 182 183 設計の全ての部分の実装はできました。 184 しかし、目指している音を鳴らすことはできませんでした。 185 反響の共鳴に関わらずリードの固有振動数らしき周波数で鳴り続ける挙動をしています。 186 パラメーターが多すぎるのか、モデルが正確ではないのか。わからないところです。 187 より謎なのは、ipynbでうまく動作していても実際にプラグインに記述するとうまく動かないところです。 188 何はともあれ、私のシミュレーターを作る力が不足していたということでしょう。