commit ca2e1595ff132f454c0c94045291b513fba0fd96
parent deeecc96f3e5d20d51c5c3c1f7e954e231a02e65
Author: minerva-jupiter <ryouturn@gmail.com>
Date: Fri, 19 Jun 2026 03:14:30 +0000
docs: update simulation theory and implementation details in article
- Rename section header "実装(TL;DR)" to "実装理論(TL;DR)".
- Rename "発信部" and "共鳴部" headers.
- Refine the explanation of the wave reflection model for the resonance section, including d'Alembert's solution and derivation of return time based on frequency.
- Simplify the implementation section to point towards the source code.
Diffstat:
1 file changed, 57 insertions(+), 35 deletions(-)
diff --git a/public/creating_a_clap_plugin_through_simulation.md b/public/creating_a_clap_plugin_through_simulation.md
@@ -29,10 +29,10 @@ ignorePublish: false
個人的にはチェット・ベイカーのような優しいトランペットの音を目指したいのですが、サックスだろうがトランペットだろうがの発振は圧力駆動型バルブというモデルで説明できます。
フルートやリコーダーは空気の渦を作って発振するので全然違うのですが、流体シミュレーションはまだ理解しきれてないので、今回はちょっと勘弁してもらって。
-# 実装(TL;DR)
+# 実装理論(TL;DR)
ここから長丁場です。実際、結構な時間を費やしています。
-## 発信部の実装
+## 発信部
まず、圧力駆動型バルブの運動方程式を示します。
変位はxで、kは弾性、P(t)は息による空気圧、rは振動の減衰のための定数です。
$$ m \frac{d^2x}{dt^2} = -k x + f(P(t),x) + r \frac{dx}{dt} $$
@@ -92,61 +92,83 @@ $$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セクションに掲載)。
+となります(一次双変換の過程や実際に運動方程式の解の近似になっていることはリポジトリのReadmeのTL;DR Derivation of the simulation formulaセクションで示しています)。
これで、$x[n]$がnステップ目の位置を表す関数として記述することができました。
-## 共鳴部の実装
+## 共鳴部
-筒の剛性や、空気の密度関数から気柱による音波の反射をシミュレートすることも考えたが発振よりも膨大になってしまうため、断念しました。
+筒の剛性や、空気の密度の変化から気柱による音波の反射をシミュレートすることも考えましたが、
音波の距離によって遅れる反射は認めるものとして、反射によるエネルギー減衰のモデルを考えることで倍音の高音成分がなくなるシミュレートを実装します。
+エネルギーは位相のt微分の2乗に比例する量であり、減衰はxのt微分の2乗に対して行います。
+エネルギーが減衰すると、高音成分から現象していくことは理解できるかと思います。
+式としては、発信部から来た変位をx_n、遅延して送られてくる変位をx_resonance、前回の合計の変位をx_prevとすると、減衰定数aを使って
-管に与えられる圧力を$p_{total}$、管の断面積をS、空気の粒子速度を$v$として、これをv|v|に比例する抵抗を受ける有効音響質量M_eとすると、$r_{loss}$を用いて
+$$ x = a (x_xrev - (x_n + x_resonance))^2 $$
-$$S \cdot p_{total}(t) = M_e \frac{dv(t)}{dt} + r_{loss} v(t) |v(t)|$$
+と書くことができます。ここで、x_resonanceは筒を往復するのにかかった時間だけ前の変位に、開管なら1を、閉管なら-1をかけたものです。
-- $M_e = \rho_0 S \Delta L$
-- $r_{loss} = \frac{1}{2} \rho_0 S$
-- $\rho_0$
+x_resonanceが管を往復するのにかかった時間だけ前の変位に反射の係数をかけたもので十分な理由を示します。
-のように書けます。
+音波の伝搬について、
-また、$p_{total}(t)$は、入射波を$p^+(t)$、反射波を$p^-(t)$とすると
+$$\frac{\partial^2 p}{\partial t^2} - c^2 \frac{\partial^2 p}{\partial x^2} = 0$$
-$$p_{total}(t) = p^+(t) + p^-(t)$$
+この一般解はダランベールの解を用いて
-となり、$v(t)$を
+$$p(x, t) = f(t - x/c) + g(t + x/c)$$
-$$v(t) = \frac{1}{Z_0} \left( p^+(t) - p^-(t) \right)$$
+と書くことができ、波の振幅の移動として書くことができます。
-と書くことができます。これを用いて$p^-(t)$は
+また、反射が起きるx=Lの点について書くと
-$$p^-(t) = p_{total}(t) - p^+(t)$$
+1. 開放端のとき
-となり、元の方程式に代入すると
+筒の端が大気圧になっているので$p=0$という条件を入れて
-$$p_{total}(t) = 2p^+(t) - Z_0 v(t)$$
+$$f(t - L/c) + g(t + L/c) = 0$$
-が得られます。これを$T$をサンプリング周期として双一次変換します。
+これを $g$ について解くと
-$$S \left( 2p^+[n] - Z_0 v[n] \right) = M_e \left( \frac{v[n] - v[n-1]}{T} \right) + r_{loss} v[n] |v[n]|$$
+$$g(t + L/c) = -f(t - L/c)$$
-これを$v[n]$について解くと、
+2. 閉口端のとき
-$$v[n] = \frac{-B_{bc} + \sqrt{B_{bc}^2 + 4 A_{bc} C_{bc}[n]}}{2 A_{bc}}$$
+筒の端で空気の速度が0という条件を入れて
-where,
-$$\begin{aligned}
-A_{bc} &= r_{loss} \\
-B_{bc} &= \frac{M_e}{T} + S Z_0 \\
-C_{bc}[n] &= 2S p^+[n] + \frac{M_e}{T} v[n-1]
-\end{aligned}$$
+$$\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を往復するのにかかる時間が求められます。
-$$p^-[n] = p^+[n] - Z_0 v[n]$$
+# 実装
-に代入することで、波の反射を書くことができます。
-v[n]=0として反射波の位相を逆位相から同位相に切り替えると固定端反射を実現できそうなものだが、よく考えればv[n]=0なら減衰が起こらないので、このモデルを考えること自体が無意味になことはすぐにわかります。
-ここらへんの代わりの実装が私には思いつかなかった(元の発想はエネルギーを0より大きく1より小さい定数をかけて減衰させるようなモデルを考えれば高周波の音からの減衰を再現できるのではないかというもの)ので、実装は開口端についてのみを行っています。
-元も子もない話をしてしまうと、ローパスフィルターを書けてしまえばそれで終わりなのだけど、それだとあまりにも面白くないと思っています。
+実装はソースコードにあるものが全てです。うまく動かないと理論を見直したり、プラグイン部分の不具合を排除するために理論式をipynbでプロットしながら検証をしたりしました。
