double_ema_filter

Exponential moving average (EMA) フィルタを 2 つ直列につないだフィルタを使って AD エンベロープを作ります。

以下に完全なテストコードと、プロットに使ったコードを掲載しています。

以下は ClangSynth で使った C++ による二重 EMA フィルタによる AD エンベロープの実装です。

二重 EMA フィルタ

伝達関数

以下は EMA フィルタの伝達関数です。 \(k\) は範囲 \([0, 1]\) のフィルタ係数です。

以下は EMA フィルタを 2 つ直列につないだときの伝達関数です。

\[ \begin{aligned} H(z) &= \frac{k}{1 + (k - 1)z^{-1}} \times \frac{k}{1 + (k - 1)z^{-1}} \\ &= \frac{k^2}{1 + 2(k - 1) z^{-1} + (k - 1)^2 z^{-2}} \end{aligned} \]

この文章では、上の伝達関数を持つフィルタのことを二重 EMA フィルタと呼ぶことにします。二重 EMA フィルタは、 EMA フィルタと同様にステップ応答がオーバーシュートしないのでエンベロープへの応用に適しています。

class DoubleEMA:
    def __init__(self, k=1, value=0):
        self.k = k
        self.v0 = value
        self.v1 = value

    def reset(self, value=0):
        self.v0 = value
        self.v1 = value

    def process(self, x0):
        self.v0 += self.k * (x0 - self.v0)
        self.v1 += self.k * (self.v0 - self.v1)
        return self.v1

インパルス応答

Wolfram Alpha を使って二重 EMA フィルタの伝達関数を逆 Z 変換したところ、以下の式が得られました。逆 Z 変換によって得られた式はインパルス応答を表しています。

ゲート信号への応答

ここでの目的はエンベロープへの応用なのでゲート信号への応答の形に式を変形します。ここではゲート信号を、時間が 0 サンプルから n サンプルの区間で振幅 1 、それ以外では振幅 0 となる信号と定義しています。以下はゲート信号の図です。

試行錯誤したところ Python3 で実装した以下の関数が、伝達関数の形で実装した二重 EMA フィルタの出力と一致しました。アイデアとしては、インパルスの現れる時間をずらした複数のインパルス応答の総和を取っています。

def doubleEmaClosed(n, k, init=0, gain=1):
    value = 0
    for i in range(0, n + 1):
        value += k * k * (i + 1) * (1 - k)**i
    return init + gain * value

n はゲートが開く時間の長さです。今回の応用ではノートオンからノートオフまでの間隔を表しています。また n が無限大のときにステップ応答となります。

数式にします。 \(V\) は初期値、 \(G\) はゲインです。総和を取り除いた形は Wolfram Alpha から取得しました。

\[ \begin{aligned} x(k, V, G)[n] &= V + G \sum_{i=0}^{n} k^2 (i + 1) (1 - k)^i \\ &= V + G \Big( k^2 (n + 1) (1 - k)^n - k n (1 - k)^n - (1 - k)^n + 1 \Big) \\ &= V + G \Big( 1 + (1 - k)^n (k^2 (n + 1) - k n - 1) \Big) \\ &= V + G \Big( 1 - (1 - k)^{n + 1} (k n + k + 1) \Big) \end{aligned} \]

AD エンベロープ

2 つの二重 EMA フィルタを乗算して、ターゲットとなる AD エンベロープ \(E\) の式を作ります。 \(k_A\) はアタック、 \(k_D\) はディケイのフィルタ係数です。

\[ \begin{aligned} E[n] &= x(k_A, 0, 1)[n] \times x(k_D, 1, -1)[n] \\ &= \Big( 1 - (1 - k_A)^{n + 1} (k_A n + k_A + 1) \Big) \Bigg( 1 - \Big( 1 - (1 - k_D)^{n + 1} (k_D n + k_D + 1) \Big) \Bigg) \\ &= \Big( 1 - (1 - k_A)^{n + 1} (k_A n + k_A + 1) \Big) (1 - k_D)^{n + 1} (k_D n + k_D + 1) \end{aligned} \]

def doubleEmaEnvelopeD0(n, k_A, k_D):
    A = (1 - k_A)**(n + 1) * (k_A * n + k_A + 1)
    D = (1 - k_D)**(n + 1) * (k_D * n + k_D + 1)
    return (1 - A) * D

def samplesToKp(timeInSamples):
    y = 1 - np.cos(2 * np.pi / timeInSamples)
    return -y + np.sqrt(y * (y + 2))

n = np.arange(10000)
k_A = samplesToKp(1000)
k_D = samplesToKp(20000)
envelope = doubleEmaEnvelopeD0(n, k_A, k_D)

以下は \(k_A\) と \(k_D\) を適当に変えたときのエンベロープの出力とコードです。縦軸は振幅、横軸はサンプル数であらわされた時間です。

ピークの探索

エンベロープとして使いたいので \(k_A\) と \(k_D\) の値によらず、ピークの値が 1 となるように正規化します。ここでは数値計算によって \(E\) の最大値を求め、逆数にして乗算することで正規化を行います。 \(E\) を微分した式を使うことも試みたのですが、得られた式が煩雑な形になったので諦めました。

以下は Python3 によるピークを正規化したエンベロープの計算を行うコードです。変数が 1 つだけの任意の関数の最小値を求めるときは scipy.optimize.minimize_scalar が使えます。求まる値はあくまでも最小値なので doubleEmaEnvelopeD0 の正負の符号を逆転した doubleEmaEnvelopeD0Negative を使っています。あとは得られた peakValue の逆数を AD エンベロープの出力に乗算すればピークの値が 1 に正規化された出力が得られます。

# `samplesToKp` と `doubleEmaEnvelopeD0` については一つ前のコードスニペットを参照。

import scipy.optimize as optimize

def doubleEmaEnvelopeD0Negative(n, k_A, k_D):
    A = (1 - k_A)**(n + 1) * (k_A * n + k_A + 1)
    D = (1 - k_D)**(n + 1) * (k_D * n + k_D + 1)
    return (A - 1) * D

k_A = samplesToKp(1000)
k_D = samplesToKp(20000)

result = optimize.minimize_scalar(
    doubleEmaEnvelopeD0Negative,
    args=(k_A, k_D),
)
peakTime = result.x
peakValue = -result.fun

n = np.arange(10000)
envelope = doubleEmaEnvelopeD0(n, k_A, k_D)

以下はピークを正規化したエンベロープの \(k_A\) と \(k_D\) を適当に変えたときの出力です。縦軸は振幅、横軸はサンプル数であらわされた時間です。

正規化の誤差

エンベロープの出力は \([0, 1]\) の範囲内であることが望ましいです。また、伝達関数が表す差分方程式による計算と、 doubleEmaEnvelopeD0 が表す閉じた形の解 (closed-form solution) による計算では浮動小数点数の挙動によって計算結果がずれることが考えられます。そこで気になってくるのが正規化したピークが 1 を超えないかどうかです。

以下は適当に \(k_A\) と \(k_D\) を変えたときの正規化したピーク値のプロットです。縦軸は正規化したピーク値、横軸は \(k_D\) の値です。縦軸のスケールは 1e-8 、オフセットは +9.9999999e-1 です。つまり左上の 1.00 を超えていなければ正規化したピークは 1 以内に収まっているということです。

傾向としては \(k_A\) と \(k_D\) が小さいときは、正規化したピークも小さくなるようです。また、調べた範囲では正規化したピークが 1 を超えるケースは見つかりませんでした。

その他

再帰のない EMA フィルタの計算式

Wolfram Alpha で EMA フィルタの伝達関数を逆 Z 変換したところ、以下の式が得られました。

二重 EMA フィルタの微分

\(x(k, V, G)[n]\) の微分です。とりあえず微分したものの、用途が見当たりませんでした。

\[ \begin{aligned} \frac{dx[n]}{dn} &= G (\log(1-k) (1-k)^n (k^2 (n+1)-k n-1)+(1-k)^n (k^2-k)) \\ &= G (k - 1) (1-k)^n ((k n + k + 1) \log(1-k) + k) \\ &= -G (1 - k)^{n + 1} ((k n + k + 1) \log(1-k) + k) \\ \end{aligned} \]

AD エンベロープの式の微分

AD エンベロープのピークの値を取得するために使えるかと思ったのですが、式の整理が困難だったので、このアプローチは諦めました。

\[ \begin{aligned} & \frac{d}{dn} \Bigg( \Big( 1 - (1 - k_A)^{n + 1} (k_A n + k_A + 1) \Big) (1 - k_D)^{n + 1} (k_D n + k_D + 1) \Bigg) \\ & \quad = k_D (1 - (1 - k_A)^{n + 1} (k_A n + k_A + 1)) (1 - k_D)^{n + 1} \\ & \qquad + (1 - (1 - k_A)^{n + 1} (k_A n + k_A + 1)) (k_D n + k_D + 1) (1 - k_D)^{n + 1} \log(1 - k_D) \\ & \qquad + (k_D n + k_D + 1) (1 - k_D)^{n + 1} (-k_A (1 - k_A)^{n + 1} - (k_A n + k_A + 1) (1 - k_A)^{n + 1} \log(1 - k_A)) \\ & \quad = + D_0 (1 - A_0 A_1) k_D \\ & \qquad + D_0 D_1 (1 - A_0 A_1) D_2 \\ & \qquad - D_0 D_1 A_0 k_A \\ & \qquad - D_0 D_1 A_0 A_1 A_2 \\ & \quad = D_0 \Bigg( (k_D + D_1 D_2) (1 - A_0 A_1) - D_1 A_0 (k_A + A_1 A_2) \Bigg) \end{aligned} \]

\[ \begin{array}{llll} & D_0 = (1 - k_D)^{n + 1} && A_0 = (1 - k_A)^{n + 1} \\ & D_1 = k_D n + k_D + 1 && A_1 = k_A n + k_A + 1 \\ & D_2 = \log(1 - k_D) && A_2 = \log(1 - k_A) \\ \end{array} \]

上の式の \(n\) を 0 について解けば AD エンベロープのピークの時点が得られるはずです。 Wolfram Alpha と Maxima で試したのですが、式の形がとても煩雑になるので整理は困難です。

def doubleEmaEnvelopeD1(n, k_A, k_D):
    D_0 = (1 - k_D)**(n + 1)
    A_0 = (1 - k_A)**(n + 1)
    D_1 = k_D * n + k_D + 1
    A_1 = k_A * n + k_A + 1
    D_2 = np.log(1 - k_D)
    A_2 = np.log(1 - k_A)
    return D_0 * ((k_D + D_1 * D_2) * (1 - A_0 * A_1) - D_1 * A_0 * (k_A + A_1 * A_2))

二重 EMA フィルタによる AD エンベロープ