exppoly_envelope

WobblingMetalBoard を作っているときに次のエンベロープを思いつきました。

\(t\) は時間、 \(\alpha,\beta\) は任意の定数です。簡単に性質を調べて実装します。

env(t) := t^α * exp(-β * t);
env_d1: diff(env(t), t);
env_d2: diff(env_d1, t);

\[ \begin{aligned} \dot{E}(t) &= (- \beta t^\alpha + \alpha t^{\alpha - 1}) e^{-\beta t} \\ \ddot{E}(t) &= ( \beta^2 t^\alpha - 2 \alpha \beta t^{\alpha - 1} + \left( \alpha - 1 \right) \alpha t^{\alpha - 2} ) e^{-\beta t} \end{aligned} \]

solve(env_d1 = 0, t);

実際に計算してみると \(t = \dfrac{\alpha}{\beta}\) のとき \(E(t)\) が最大となります。つまりアタックの長さは \(\dfrac{\alpha}{\beta}\) です。

アタックの長さから振幅を \([0, 1]\) に正規化した \(\hat{E}(t)\) を求めることができます。

import numpy

def expPoly(α, β, time):
    return time**α * numpy.exp(-β * time)

def expPolyNormalized(α, β, time):
    normalize = expPoly(α, β, α / β)
    return time**α * numpy.exp(-β * time) / normalize

time = numpy.linspace(0, 16, 1024)
envelope = expPolyNormalized(1, 1, time)

\(\beta = 1\) に固定して \(\alpha\) を変えたときの \(\hat{E}(t)\) のプロットです。

\(\alpha = 1\) に固定して \(\beta\) を変えたときの \(\hat{E}(t)\) のプロットです。

\(E(t)\) は \(t\) が \(+\infty\) に到達すると \(0\) になります。適当に計算を打ち切るために \(\hat{E}(t) = x\) となる時間 \(t\) を求めます。 \(x = \dfrac{t^\alpha e^{-\beta t}}{E \left( \dfrac{\alpha}{\beta} \right)}\) を \(t\) について解いてみます。この問題は Maxima では解けなかったので Wolfram Alpha を使いました。

結果です。 \(W\) は Lambert W-function あるいは product log function と呼ばれている関数です。

\[ t = -\frac{\alpha}{\beta} W \left( -\frac{\beta}{\alpha} \left( x \eta \right)^{1 / \alpha} \right), \quad \eta = E \left( \frac{\alpha}{\beta} \right) \]

import numpy
import matplotlib.pyplot as pyplot
import scipy.special.lambertw as lambertw

def getTime(α, β, x, normalize, k=0):
    return -α * lambertw(-β * (x * normalize)**(1 / α) / α, k) / β

def expPoly(α, β, time):
    return time**α * numpy.exp(-β * time)

def expPolyNormalized(α, β, time):
    normalize = expPoly(α, β, α / β)
    return (
        time**α * numpy.exp(-β * time) / normalize,
        normalize,
    )

samplerate = 48000
duration = 8

time = numpy.linspace(0, duration, duration * samplerate)
alpha = 1
beta = 1

curve, normalize = expPolyNormalized(alpha, beta, time)

xx = numpy.linspace(1, 0, 10)
tt0 = getTime(alpha, beta, xx, normalize, 0)
value_tt0, _ = expPolyNormalized(alpha, beta, tt0)

# k = 0 以外の場合は省略。

pyplot.plot(time, curve)
pyplot.scatter(tt0, value_tt0, color="red", label="k=0")
pyplot.show()

結果です。 \(k\) は Lambert W-function の引数で、適当な整数です。図を見ると、振幅から時間を求めるときはアタックの区間で \(k=0\) 、ディケイの区間では \(k=-1\) を使えば良さそうです。

\(E(t)\) に含まれる \(t^\alpha\) は浮動小数点数で計算するとそのうち inf になります。そこで \(t\) がどこまで進むと \(t^\alpha \approx +\infty\) となるのかを調べます。

\(+\infty\) になる時点を \(t_{\mathrm{max}}\) 、浮動小数点数で表現できる最大値を \(F\) とすると、次の式が立ちます。

\(\alpha\) が 1 より小さいとき、 \(t_{\mathrm{max}}\) は \(F\) より大きいので、 \(t\) が nan や inf でなければ取り得る全ての値の範囲で計算に問題が無いことが分かります。 \(\alpha\) が 1 以上のときの \(t_{\mathrm{max}}\) を計算してみます。

import numpy

def printTimeMax(dtype):
    fmax = numpy.finfo(dtype).max
    alpha = numpy.array([2**i for i in range(1, 16)], dtype=dtype)
    timeMax = numpy.power(fmax, dtype(1) / alpha, dtype=dtype)
    print(dtype)
    for a, t in zip(alpha, timeMax):
        print(f"{a:9.1f}, {t}")

printTimeMax(numpy.float32)
printTimeMax(numpy.float64)

整形した出力です。左の列が \(\alpha\) 、右の列が \(t_{\mathrm{max}}\) です。単精度のときは \(\alpha\) が 64.0 のときに 4.0 秒しかレンダリングできないことが分かります。

正規化係数 \(E \left( \dfrac{\alpha}{\beta} \right)\) を求めることができる \(\beta\) の範囲です。

\[ \begin{aligned} t_{\mathrm{max}} &\geq \dfrac{\alpha}{\beta} \\ \beta &\geq \dfrac{\alpha}{F^{1 / \alpha}} \end{aligned} \]

#include <cmath>
#include <cfloat>

class ExpPolyEnvelope {
public:
  // attack の単位は秒。
  // curve は任意の値。 β に相当。
  // attack と curve が大きいと計算結果が inf になるときがあるので注意。
  void reset(double sampleRate, double attack, double curve)
  {
    alpha = attack * curve;

    auto betaMin = alpha / getTerminationTime();
    if (curve < betaMin) curve = betaMin;

    peak = pow(alpha / curve, alpha) * exp(-alpha);
    gamma = exp(-curve / sampleRate);
    tick = 1.0 / sampleRate;

    time = 0.0;
    value = 1.0;
  }

  bool isReleasing() { return time >= attack; }
  double getTerminationTime() { return pow(DBL_MAX, 1.0 / alpha); }

  double process()
  {
    auto output = pow(time, alpha) * value / peak;
    if (!std::isfinite(output)) return 0.0; // 念のため。
    time += tick;
    value *= gamma;
    return output;
  }

protected:
  double value = 0;
  double peak = 1;
  double gamma = 0;
  double attack = 0;
  double tick = 0;
  double alpha = 0;
  double time = 0;
};

ExpPoly エンベロープを GenericDrum の衝突の計算に応用できそうな気がしたので、以下の区間積分が解けるかどうか試しました。

\[ \begin{equation} \int_\tau^\infty t^\alpha e^{-\beta t} dt, \quad \tau \geq 0. \label{a} \end{equation} \]

ここでの応用はエネルギーを時間方向に薄く延ばすことです。つまり、任意の正の実数 \(x\) が与えられたときに、式 \(\ref{a}\) と \(x\) が等しくなるような \(\alpha, \beta, \tau\) の組を知りたいです。

SymPy で積分します。 Maxima の integrate では解けず、 Wolfram Alpha では計算時間切れになりました。

import sympy

t = sympy.Symbol("t", real=True, positive=True)
α = sympy.Symbol("α", real=True, positive=True)
β = sympy.Symbol("β", real=True, positive=True)
τ = sympy.Symbol("τ", real=True, positive=True)

result = sympy.integrate(t**α * sympy.E ** (-β * t), (t, τ, sympy.oo))
print(sympy.latex(result))

\[ τ τ^{α} \left(\frac{β^{- α - 1} τ^{- α - 1} \left(α + 1\right) \Gamma\left(- α - 1\right) \gamma\left(α + 1, β τ\right)}{\Gamma\left(- α\right)} + β^{- α - 1} τ^{- α - 1} \Gamma\left(α + 1\right)\right) \]

\[ \begin{aligned} \int_\tau^\infty t^\alpha e^{-\beta t} dt &= τ^{c} \left( \frac{ s c \Gamma(-c) \gamma(c, β τ) }{ \Gamma(- α) } + s \Gamma(c) \right), \\ c &= \alpha + 1,\\ s &= β^{-c} τ^{-c}. \end{aligned} \]

上の SymPy のコードに以下をつけ足して解こうとしたのですが、 NotImplementedError が出て行き詰りました。

x = sympy.Symbol("x", real=True, positive=True)
sympy.solve(sympy.Eq(result, x), τ)