commit 196d4123245789dc52347faaba8ac5ba3bedd1dd
parent 95f9f9fa23b5d607f7e753cfc3d4f452f6a942ce
Author: minerva-jupiter <ryouturn@gmail.com>
Date: Wed, 17 Jun 2026 22:27:57 +0900
docs: add mathematical derivation for reed and resonance simulation
Update documentation with detailed derivation of simulation formulas for air pressure, reed mechanics, and resonance reflection in CLAP plugin development.
Diffstat:
1 file changed, 110 insertions(+), 1 deletion(-)
diff --git a/public/creating_a_clap_plugin_through_simulation.md b/public/creating_a_clap_plugin_through_simulation.md
@@ -1,7 +1,7 @@
---
title: シミュレーションによるCLAPプラグインの作成
tags:
- - ''
+ - nih-plug
private: true
updated_at: ''
id: null
@@ -41,3 +41,112 @@ rの項の妥当性について一応述べておくと、例えば吹くのを
もしリードが息に押されてマウスピースに近づくと、その近さの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ステップ目の位置を表す関数として記述することができた。
+
+## 共鳴部の実装
+
+筒の剛性や、空気の密度関数から気柱による音波の反射をシミュレートすることも考えたが発振よりも膨大になってしまうため、断念した。
+音波の距離によって遅れる反射は認めるものとして、反射によるエネルギー減衰のモデルを考えることで倍音の高音成分がなくなるシミュレートを実装した。
+
+管に与えられる圧力を$p_total$、管の断面積をS、空気の粒子速度を$v$として、これをv|v|に比例する抵抗を受ける有効音響質量M_eとすると、$r_{loss}$を用いて
+
+$$S \cdot p_{total}(t) = M_e \frac{dv(t)}{dt} + r_{loss} v(t) |v(t)|$$
+
+- $M_e = \rho_0 S \Delta L$
+- $r_{loss} = \frac{1}{2} \rho_0 S$
+- $\rho_0$
+
+のように書ける。
+
+また、$p_{total}(t)$は、入射波を$p^+(t)$、反射波を$p^-(t)$とすると
+
+$$p_{total}(t) = p^+(t) + p^-(t)$$
+
+となり、$v(t)$を
+
+$$v(t) = \frac{1}{Z_0} \left( p^+(t) - p^-(t) \right)$$
+
+と書くことができる。これを用いて$p^-(t)$は
+
+$$p^-(t) = p_{total}(t) - p^+(t)$$
+
+となり、元の方程式に代入すると
+
+$$p_{total}(t) = 2p^+(t) - Z_0 v(t)$$
+
+が得られる。これを$T$をサンプリング周期として双一次変換する。
+
+$$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]|$$
+
+これを$v[n]$について解くと、
+
+$$v[n] = \frac{-B_{bc} + \sqrt{B_{bc}^2 + 4 A_{bc} C_{bc}[n]}}{2 A_{bc}}$$
+
+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}$$
+
+これを
+
+$$p^-[n] = p^+[n] - Z_0 v[n]$$
+
+に代入することで、波の反射を書くことができる。
+v[n]=0として反射波の位相を逆位相から同位相に切り替えると固定端反射を実現できそうなものだが、よく考えればv[n]=0なら減衰が起こらないので、このモデルを考えること自体が無意味になる。
+ここらへんの代わりの実装が私には思いつかなかった(元の発想はエネルギーを0より大きく1より小さい定数をかけて減衰させるようなモデルを考えれば高周波の音からの減衰を再現できるのではないかというもの)ので、実装は開口端についてのみになっている。
+元も子もない話をしてしまうと、ローパスフィルターを書けてしまえばそれで終わりなのだけど、それだとあまりにも面白くない。
