MPMアルゴリズムのNSDFの求めかた

結論だけ知りたい人向け

NSDF の式

$\displaystyle n'_t(\tau) = \frac{2r'_t(\tau)}{m'_t(\tau)} \quad (0 \leq \tau < W)$

の $m'_t(\tau)$ について、前から求めると

$\displaystyle m'_t(\tau) = \begin{cases} 2r'_t(0) & (\tau=0) \\ m'_t(\tau-1) - x_{t+\tau-1}^2 - x_{t+W-\tau}^2 & (0 < \tau < W) \end{cases}$

後ろから求めると

$\displaystyle m'_t(\tau) = \begin{cases} 0 & (\tau=W) \\ m'_t(\tau+1) + x_{t+\tau}^2 + x_{t+W-\tau-1}^2 & (0 \leq \tau < W) \end{cases}$

となるので、好きな方でやってください。

実装（C#）

https://github.com/azyobuzin/PitchDetector/blob/81b5b4e3ff5cf20dce2b1f01fd36294360323722/PitchDetector/PitchAccord.cs#L35-L46

はじめに

ピッチ検出アルゴリズムについて検索すると、まぁ大体このページ（音声の波形からピッチを検出するアルゴリズム - まめめも）に行き着くと思うのですが、 NSDF のインクリメンタルな求め方について特に説明もなく、クワイン用の well compressed なプログラムしか書いていないので、基の論文であるところの『A Smarter Way to Find Pitch』に書いてある式を変形しながら、本当にインクリメンタルに求まるのか検証していきます。

ここでは Normalized Square Difference Function の略として NSDF と言っていますが、 McLeod 氏の 2008 年の論文では specially-normalized autocorrelation function で略して SNAC と呼んでいますね。まぁいいや。

前から求めていく

まず、 NSDF は自己相関から得られる値で、最大値は 1 です。 $\tau=0$ のとき、元の信号と完全に一致するわけですから、 $n'_t(0)=1$ になるはずで、 $n'_t(\tau) = \frac{2r'_t(\tau)}{m'_t(\tau)}$ だから $m'_t(0) = 2r'_t(0)$ になります。

次に $m'_t(\tau)$ のもともとの定義を見てみましょう。

$\displaystyle m'_t(\tau) = \sum_{j=t}^{t+W-\tau-1} (x_j^2 + x_{j+\tau}^2)$

話を簡単にするために $t=0,W=4$ としてみます。そしてシグマを展開すると次のようになります。

$\displaystyle \begin{eqnarray} m'_0(0) &=& \sum_{j=0}^3 (x_j^2 + x_{j+\tau}^2) \\ &=& (x_0^2 + x_0^2) + (x_1^2 + x_1^2) + (x_2^2 + x_2^2) + (x_3^2 + x_3^2) \\ &=& 2x_0^2 + 2x_1^2 + 2x_2^2 + 2x_3^2 \\ m'_0(1) &=& \sum_{j=0}^2 (x_j^2 + x_{j+\tau}^2) \\ &=& (x_0^2 + x_1^2) + (x_1^2 + x_2^2) + (x_2^2 + x_3^2) \\ &=& x_0^2 + 2x_1^2 + 2x_2^2 + x_3^2 \\ &=& m'_0(0) - x_0^2 - x_3^2 \\ m'_0(2) &=& \sum_{j=0}^1 (x_j^2 + x_{j+\tau}^2) \\ &=& (x_0^2 + x_2^2) + (x_1^2 + x_3^2) \\ &=& x_0^2 + x_1^2 + x_2^2 + x_3^2 \\ &=& m'_0(1) - x_1^2 - x_2^2 \\ m'_0(3) &=& \sum_{j=0}^0 (x_j^2 + x_{j+\tau}^2) \\ &=& x_0^2 + x_3^2 \\ &=& m'_0(2) - x_2^2 - x_1^2 \end{eqnarray}$

で、これを眺めて雑に一般化してしまえば

$\displaystyle m'_0(\tau) = m'_0(\tau-1) - x_{\tau-1}^2 - x_{W-\tau}^2$

になるなぁというのがわかると思います。というわけで

$\displaystyle m'_t(\tau) = \begin{cases} 2r'_t(0) & (\tau=0) \\ m'_t(\tau-1) - x_{t+\tau-1}^2 - x_{t+W-\tau}^2 & (0 < \tau < W) \end{cases}$

が得られました。 NSDF の値自体は $n'_t(\tau) = \frac{2r'_t(\tau)}{m'_t(\tau)}$ に求まった $m'_t(\tau)$ を突っ込むだけなので、特に言うことはないと思います。

後ろから求めていく

普通実装するにあたって、すでに信号は読み込まれているはずなので、逆から舐めていっても問題ないはずです。そこで先ほど展開した $m'_0(\tau)$ をもう一度眺めてみると

$\displaystyle \begin{eqnarray} m'_0(3) &=& x_0^2 + x_3^2 \\ m'_0(2) &=& x_0^2 + x_1^2 + x_2^2 + x_3^2 \\ &=& m'_0(3) + x_2^2 + x_1^2 \\ m'_0(1) &=& x_0^2 + 2x_1^2 + 2x_2^2 + x_3^2 \\ &=& m'_0(2) + x_1^2 + x_2^2 \\ m'_0(0) &=& 2x_0^2 + 2x_1^2 + 2x_2^2 + 2x_3^2 \\ &=& m'_0(1) + x_0^2 + x_3^2 \end{eqnarray}$