アジョブジ星通信

進捗が出た頃に更新されるブログ。

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}

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

はじめに

ピッチ検出アルゴリズムについて検索すると、まぁ大体このページ(音声の波形からピッチを検出するアルゴリズム - まめめも)に行き着くと思うのですが、 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}

という感じなので、頭悪く m'_0(4)=0 と置いてしまえば

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

になります。というわけで

\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}

になります。

私は 0 \leq \tau < W のループできれいに書ける逆順のほうが好きです。こちらからは以上です。