MPMによる基本周波数の測定 - TadaoYamaokaの開発日記

AndroidとiOSでリリースしている「Vocal Pitch Monitor」では、自己相関関数と決定木のアルゴリズムで基本周波数(ピッチ)の測定を行っている。

Vocal Pitch Monitor - Apps on Google Play

‎「VocalPitchMonitor」をApp Storeで

自己相関関数は、ノイズに強く比較的安定して高速にピッチを測定できるが、

スペクトル漏れの問題
サブハーモニックエラー(高調波をピッチと誤る)

が起きる。

前者は窓関数で軽減できるが、高調波がサンプリング間隔の定数倍で、基本周波数が定数倍から外れる場合、高調波が強調されて、サブハーモニックエラーが起きやすくなる。

また、自己相関関数を固定ウィンドウ幅で計算する場合、テーパリング(先細り)が発生するため、基本周波数のピークが減衰して高調波よりも低い値になる場合がある。

アルゴリズム選定の経緯

アプリの元となったWindows用のソフトを開発したのは、2009年で、現在のようにインターネットで簡単に論文を検索できなかった。
はじめナイーブにFFTからピークを抽出していたが、精度が出なかったため、図書館で「コンピュータ音楽: 歴史・テクノロジー・アート」という本を見つけて、紹介されていたアルゴリズムをいくつか実装して、自己相関関数が高速で安定していたため、自己相関関数を採用した。

ただ、上記の通り、サブハーモニックエラーが起きるため、決定木を機械学習して精度を上げている。
このときに機械学習を勉強したのが今の将棋AIの開発にもつながっている。

LLMで調査

もっと精度のよいアルゴリズムがないか気になったので、LLMで調べてみた。

YIN

多数の手法が回答されたが、スマホで動かす際に使えるノンパラメトリックの手法では、YINというアルゴリズムが見つかった。
YINは、自己相関関数をベースにして、オクターブ誤りを大幅に減らした手法ということである。

インターネット上で、論文が見つかったので、読んでみた。
http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf

2002年の論文で、比較的古いが、上記の書籍は2001年出版なので、まだ紹介されていなかった。

自己相関関数の代わりに、累積平均正規化差分関数(Cumulative mean normalized difference function)を使うことで、自己相関よりも振幅変換に強く、ピークの検索範囲上限設定がなくなるというメリットがある。
ただし、固定ウィンドウ幅の自己相関のようにFFTが使えず、Nの2乗の計算が必要になるため、スマホで動かすには向いていなさそうである。
FFTを使う場合は誤差が生じるため、既存手法の方が良いと論文にも書かれていた。

McLeod Pitch Method (MPM)

YINを検索したら、似たような手法で、McLeod Pitch Method (MPM)というものが見つかったので、論文を読んでみた。
https://www.cs.otago.ac.nz/graphics/Geoff/tartini/papers/A_Smarter_Way_to_Find_Pitch.pdf

こちらは2005年の論文である。
オクターブエラーが起こりにくく、窓関数がなくてもよいため、小さなウィンドウ幅にできるため、時間分解能を高くできると提案されている。
固定ウィンドウの自己相関関数で計算できるので、スマホでも高速に計算できそうである。

アルゴリズムの詳細は割愛するが、固定ウィンドウ幅の自己相関関数をエネルギーで正規化してテーパリングを打ち消した手法である。
正規化した関数は、Normalised Square Difference Function (NSDF)と呼ばれる。

Pythonで実装

MPMアルゴリズムをPythonで実装して可視化してみた。

高調波とノイズを含む信号の、FFTと自己相関とNSDを表示する。

import numpy as np
import matplotlib.pyplot as plt

# サンプリング周波数とサンプル数の設定
fs = 8000            # サンプリング周波数 [Hz]
N = 4096             # サンプル数
t = np.arange(N) / fs  # 時間軸

# 基本周波数の設定
f0 = 59              # 基本周波数 [Hz]

# 波形の生成：基本周波数 + 2倍、3倍の高調波 + ホワイトノイズ
y = (0.2 * np.sin(2 * np.pi * f0 * t) +
     0.4 * np.sin(2 * np.pi * 2 * f0 * t) +
     4.0 * np.sin(2 * np.pi * 3 * f0 * t) +
     0.2 * np.random.randn(N))

# FFTの計算
Y = np.fft.fft(y)
freqs = np.fft.fftfreq(len(y), d=1/fs)

# プロット用に正の周波数部分のみを抽出
half_n = len(Y) // 2
positive_freqs = freqs[:half_n]
Y_positive = Y[:half_n]

# 自己相関の計算（全範囲）
autocorr_full = np.correlate(y, y, mode='full')
lags_full = np.arange(-len(y) + 1, len(y))

# 正のラグ部分のみを抽出（0ラグ以上）
mask = lags_full >= 0
autocorr = autocorr_full[mask]
lags = lags_full[mask]

#MPM
m1 = (y * y)[::-1].cumsum()
m1 = np.roll(m1, 1)
m1[0] = 0
m1 = autocorr[0] - m1

m2 = (y * y).cumsum()
m2 = np.roll(m2, 1)
m2[0] = 0
m2 = autocorr[0] - m2

nsd = 2 * autocorr / (m1 + m2)

# プロットの設定
fig, axes = plt.subplots(4, 1, figsize=(10, 12))

# 1. 波形のプロット
axes[0].plot(t, y)
axes[0].set_title("Waveform")
axes[0].set_xlabel("Time [s]")
axes[0].set_ylabel("Amplitude")
axes[0].xaxis.set_major_locator(plt.MaxNLocator(20))  # 目盛りを細かく設定
axes[0].minorticks_on()  # 小目盛りを表示

# 2. FFTのプロット（振幅スペクトル）
axes[1].plot(positive_freqs, np.abs(Y_positive))
axes[1].set_title("FFT Spectrum")
axes[1].set_xlabel("Frequency [Hz]")
axes[1].set_ylabel("Amplitude")
axes[1].xaxis.set_major_locator(plt.MaxNLocator(20))
axes[1].minorticks_on()

# 3. 自己相関のプロット（正のラグ部分）
axes[2].plot(lags, autocorr)
axes[2].set_title("Autocorrelation (Positive Lags)")
axes[2].set_xlabel("Lag")
axes[2].set_ylabel("Correlation")
axes[2].xaxis.set_major_locator(plt.MaxNLocator(20))
axes[2].minorticks_on()

# 4. NSDのプロット
axes[3].plot(lags, nsd)
axes[3].set_title("Normalized Square Difference")
axes[3].set_xlabel("Lag")
axes[3].set_ylabel("NSD")
axes[3].xaxis.set_major_locator(plt.MaxNLocator(20))
axes[3].minorticks_on()

plt.tight_layout()
plt.show()