基本周波数系列の推定結果を凸最適化によりポストフィルタする話

はじめに

CA社の記事があった．

AI Lab、音声・音響信号処理分野のトップカンファレンス「ICASSP 2026」にて2本の論文採択 | 株式会社サイバーエージェント

同記事では下記論文が採択された旨が書いてあった．

[2602.01727] Voting-based Pitch Estimation with Temporal and Frequential Alignment and Correlation Aware Selection

基本周波数系列を推定する手法は数多く提案されている．それら手法の推定結果の median を取るヒューリスティックな後処理（多数決，voting）によって，推定誤りを補完し，頑健な推定が可能となることが知られている．本論文では，voting により推定精度が改善する根拠および voting の妥当性を， median や mode に着目して理論的に示しつつ，さらなる精度改善の手法を提案している．

わたしは複数の推定結果の voting という後処理、いわばポストフィルタに興味を惹かれた．特に median による統合処理には凸最適化が使えそうと分かり，ちょっと手元で推定を試してみたのであった．本記事はその記録である．

はじめに
定式化
- 補題
- なぜ L1 ノルムが適切なのか？
実験
- 実験その1
  - 実験条件
  - 実験結果
- 実験その2
  - 実験条件
  - 実験結果
おわりに

定式化

まず多数の推定器の結果をアンサンブルするために，基本周波数の推定に関する voting では median を，また有声/無声判定の推定では mode を取っている．その根拠・妥当性は論文に書かれているので詳細は割愛する（サボり）．一応，簡易メモを残しておくと，

median（基本周波数）：推定の誤差分散を具体的に導出し，相関の低い推定器を取り込むことで誤差分散が $1/n$ オーダー（ $n$ は手法数）で減少することを示している．相関の高い推定器を取り込む場合（＝どの推定器も同じ間違いをする）は，誤差分散は単独の手法（ $n = 1$ ）と変わらない．
mode（有声/無声判定）：二値分類問題として扱われるが，mode による統合はコンドルセの陪審定理（Condorcet's Jury Theorem）によって正当化される．多数決が正解する確率を具体的に示し，推定器の個数 $n$ が増えるにつれて，その確率が 1 に収束することを述べている．

論文ではさらに推定結果を向上させるために「時間シフト」と「周波数シフト」を適用して誤差分布を中心化し，補正済みのデータに対して median を取ることで，誤差分散の縮小を狙っている．

本記事は，研究の発展を見すえて，これらシフトを適用する前段階の medianを取る処理（基本周波数が対象）を，凸最適化に置き換える試みである．

各手法の推定結果に対する median の計算は，L1ノルム（差の絶対値の和）の最小化に等しい（後述の「補題」参照）．凸最適化問題として定式化すると，次式となる．

$\begin{align*} \hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left[ \sum_{t=1}^{T} \sum_{i=1}^{n} | y_{i,t} - \theta_t | + \lambda \sum_{t=2}^{T} | \theta_t - \theta_{t-1} | \right] \end{align*}$

$y_{it}$ は手法 $i$ による時刻 $t$ の観測値（推定値）， $\theta_{t}$ は推定対象となる時刻 $t$ の基本周波数， $n$ は手法数， $T$ はフレーム数である． $\lambda$ は重みパラメタである．系列全体をベクトルを使って $\boldsymbol{\theta} = (\theta_{1}, \theta_{2}, \ldots, \theta_{T} )^{\top}$ と書いている．

第1項は観測誤差をL1ノルム（絶対誤差の和）で測っており，推定値を観測データの median に近づける効果がある．データの忠実性（fidelity）を評価しているとも見なせる．また第2項は正則化項であり，フレーム間におけるピッチの変化量を測りつつ，総和を取ることで，フレーム全体の総変動（total variation）を評価している．この項の導入により，系列の不要な振動を抑えて平坦にしつつ，急激な変化も許容する効果を持つ．

ちなみにベクトル記法に移れば，次式となる．

$\begin{align*} \hat{\boldsymbol{\theta}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left( \sum_{i=1}^{n} \| \mathbf{y}_i - \boldsymbol{\theta} \|_1 + \lambda \| \mathbf{D} \boldsymbol{\theta} \|_1 \right) \end{align*}$

$\| \cdot \|_1$ はL1ノルムである．また $\mathbf{D}$ は差分行列である：

$\begin{align*} \mathbf{D} \boldsymbol{\theta} = (\theta_2 - \theta_1, \ldots, \theta_T - \theta_{T-1})^{\top} \end{align*}$

補題

median を取る処理と，L1ノルムの最小化が等しいことを示しておく．

推定したい値を $\theta$ ，観測値を $x_1, \ldots, x_n$ とする．目的関数を「L1ノルム」（絶対誤差の和）と置く．

$\begin{align*} J (\theta) = \sum_{i=1}^{n} |x_i - \theta| \end{align*}$

この $J(\theta)$ を最小にする $\theta$ を求めたい．そこで $\theta$ で微分してみる．

$\begin{align*} \frac{d }{d\theta} J (\theta) = \sum_{i=1}^{n} \frac{d }{d\theta} |x_i - \theta| = \sum_{i=1}^{n} - \textrm{sign}(x_i - \theta) \end{align*}$

ここで絶対値 $|x|$ の微分が符号関数 $\textrm{sign} (x)$ になることを使った（ $x \gt 0$ なら $+1$ ， $x \lt 0$ なら $-1$ ）．上記の導関数を $=0$ と置いてみる．

$\begin{align*} \sum_{i=1}^{n} \textrm{sign}(x_i - \theta) = 0 \end{align*}$

つまり最適な $\theta$ は，「 $x_{i} \gt \theta$ となるデータの個数」（プラスの符号）と「 $x_{i} \lt \theta$ となるデータの個数」（マイナスの符号）が釣り合うような位置に相当する．これは median の定義そのものであり，まさに分布の「中央」に位置しているというわけである．

なぜ L1 ノルムが適切なのか？

基本周波数（ピッチ）の推定において，L2（平均値, mean）ではなくL1（中央値, median）を選ぶ理由は何だろうか．つまるところ、「外れ値への頑健性」にあるといえる．

ピッチ推定では倍ピッチ・半ピッチと呼ばれる推定誤りがしばしば発生する．平均値を計算すると，これら推定誤りの影響を受けやすい．例えば，各推定器の出力が $[100, 102, 200]$ (Hz) だとすると，平均値は 134 (Hz) になり，倍ピッチ誤り（＝外れ値）200 (Hz) の影響を強く受けている．これは L2 の明確な弱点である．それに対して中央値は 102 (Hz) になり，外れ値の影響を受けづらい．

実験

実験その1

実験用にダミーデータを生成し，凸最適化が期待通りに動作しているかを確認する．

実験条件

ダミーデータ（真のピッチ）の仕様：

総フレーム数 200，実区間は $[0, 4 \pi ]$
第50フレームまでは $100 + 10 \sin(t)$
第51フレーム〜第100フレームは $150$ （定数）→ジャンプを表現
第101フレーム〜第150フレームは $150 + 20 (t - t_{100}) / 10$ 　→ $t_{100}$ はほぼ $2\pi$
第151フレーム〜第200フレームは $120$ （定数）

観測データ（つまり各推定器の推定値）は真のピッチにノイズを加える．具体的にはラプラスノイズ（loc=0.0, scale=2.0）を加えつつ，オクターブエラーを表現するため 10 %の確率で +1 オクターブ（2倍）もしくは -1 オクターブ（0.5倍）のミスをするようにノイズ（外れ値）を加えた．以下にデータ生成の実装を抜粋したものを示す．

# 1. 真のピッチ軌跡
t = np.linspace(0, 4 * np.pi, n_frames)
true_pitch = np.zeros(n_frames)
true_pitch[:50] = 100 + 10 * np.sin(t[:50])  # Vibrato
true_pitch[50:100] = 150  # Jump
true_pitch[100:150] = 150 + 20 * (t[100:150] - t[100]) / 10  # Glissando
true_pitch[150:] = 120  # Jump down

# 2. 各手法による推定値
observations = []
for i in range(n_methods):
    # 基本ノイズ (ラプラス分布：裾が重い)
    noise = np.random.laplace(loc=0.0, scale=2.0, size=n_frames)

    # オクターブエラー (確率的に大幅にズレる外れ値を追加)
    # 10%の確率で、+1オクターブ(2倍) or -1オクターブ(0.5倍) のミスをする
    outliers = np.random.choice([1.0, 2.0, 0.5], size=n_frames, p=[0.9, 0.05, 0.05])

    obs = true_pitch * outliers + noise
    observations.append(obs)

観測系列の数 n_methods は 5 とした．

凸最適化の実装には，cvxpy を用いた．実装部分を抜粋する：

import cvxpy as cp

def solve_l1_filtering(observations, lam=5.0):
    n_methods, n_frames = observations.shape

    # 最適化変数 (推定したい真のピッチ)
    x = cp.Variable(n_frames)

    # 1. データ忠実項
    # L1ノルムを使うことで、観測データの中央値（median）に近づこうとする
    # sum_{i} |x - y_i|
    data_fidelity = 0
    for i in range(n_methods):
        data_fidelity += cp.norm(x - observations[i], 1)

    # 2. 正則化項
    # total variation (TV) = 隣接フレームの差分のL1ノルム
    # これにより、急激な変化を許容しつつ、不要な振動を抑える
    regularization = cp.norm(cp.diff(x), 1)

    # 目的関数
    objective = cp.Minimize(data_fidelity + lam * regularization)

    # 制約なし最適化問題として定義
    problem = cp.Problem(objective)

    # ソルバーで解く (ECOSやOSQPなどが自動選択される）
    problem.solve()

    return x.value

lam は正則化の重みである．最適化対象となる目的関数に応じたデータ形式で用意し，あとは Minimize を呼び出せばよいので，とても直感的に使うことができる．

実験結果

図1にピッチ推定の結果を示す．正則化の重みは 5.0 としている．真のピッチ（ground truth, 黒点線）に対して，ノイズありデータ系列はかなり激しく「暴れている」様子が見て取れる．青線が凸最適化を行わず単純に median だけで推定した結果を示し，赤線が今回の凸最適化で推定した結果を示している．median のみでもノイズの影響をそれほど受けずに推定できているが，一部でパルス的な観測値がそのまま残った形になっている．隣接フレームの影響を考慮していないから，といえるだろう．それに対して凸最適化の結果，正則化項の効果により，パルス的な推定値も消え，真のピッチに近い推定値が得られた．また前半部分のビブラート的な波打つ変動にも追従できている．

図2は正則化項の重みを 100.0 にしたときの結果を示している．スパース性が一層強く出ており，よりノート（音符）らしい推定結果，つまり一部の急激なジャンプ（第50フレーム，第150フレーム）を除いて平坦な系列が得られている．

実験その2

実際の音声データを対象にして，推定実験を行う．

実験条件

音声データは The LJ Speech Dataset の "LJ001-0001.wav" とした．

ピッチ推定には librosa から使える pYIN のみを使い，異なるピッチ推定器を複数個使わない．ただしピッチの抽出条件を変えて，擬似的な複数の推定器を用意した．具体的には分析窓の長さを 2048をデフォルトとして，他を1024，4096 としている．窓を短くすれば（1024）時間分解能が上がるが周波数分解能が下がるので，低音部分の抽出が不安定になる．また長くすれば（4096）安定するが細かい変化には追従しづらい．

実験結果

図3にピッチ推定の結果を示す．正則化の重みは 1.0 とした．スペクトログラムに見える縞模様のうち「一番下」のものをトレースできていれば良い．0.3 〜 0.4 秒付近にかけて，pYIN (win=4096) が推定誤りを起こしていることが分かる．また 1.8 〜 1.9 秒付近では pYIN (win=1096) のみが正しいと思われるピッチ系列を推定しているが，他の推定器では追随できておらず，多数決の結果として平坦な系列が得られている．その他のフレームを見ても，概ね各推定器の結果を統合して滑らかにしたような軌跡が得られている．

正則化の重みを 10.0 に変えて推定したときの結果を図4に示す．細かな上下動が抑制され，滑らかなピッチ軌跡が得られている．

最後に，正則化の重みを 50.0 に変えて推定したときの結果を図5に示す．大きな上下動を起こす箇所はわずかであり，ほぼ平坦な軌跡が得られている．これはデメリットではなく，採譜の際にはむしろ有効に働くことが期待される．

おわりに

凸最適化による基本周波数（ピッチ）推定の後処理はうまく動いたといえる．今回は pYIN の結果のみを統合したが，論文で試されていた他の推定器の結果を統合した実験は実施しておきたいところである（アライメントの問題は適切に対処する前提で）．

ピッチ推定手法をさらに発展させるために，以下の方向性が考えられる．

各推定器の絶対誤差に対して，推定器の信頼度を反映する重みを導入する方法である．重みの決め方には色々なものが考えられる．有声/無声判定にも重みは導入できる（重み付き多数決）．
現状の枠組みは，真値の周りに，単一のラプラス分布に従う誤差が発生しているモデルとして解釈できる．ラプラス分布の尤度を最大化することで，median が得られたのであった．これを拡張するため，ラプラス混合分布を導入する．真値の分布以外に，倍ピッチ・半ピッチ誤りに対応する分布を置く．これにより，推定値が真値から生成される過程をうまくモデル化できることが期待される（生成モデル）．
現状の枠組みは，フレームごとに独立して計算し，後処理的にスムージングをしていた．これを時系列モデルや状態空間モデルに組み込む方向性も考えられる（例：粒子フィルタ，隠れマルコフモデル）．論文で議論されていた時間方向・周波数方向のアライメントの問題を，確率モデルによってうまく解決できる可能性がある．
論文のセクション3.2では「Greedy Algorithm による相関を考慮した手法選択」を提案していた．これは事前の静的な選択であるが，「曲のAメロでは Harvest が良いが、Bメロの高音では pYIN が良い」といった時間変化する信頼度には対応できない．そこで Group Lasso による動的な重み付けを導入する方向性が考えられる．
ピッチの軌跡を，「ゆっくり変化する曲のメロディ」と「ビブラートなどの周期的な揺れ」，そして「突発的なエラー（スパースなノイズ）」に分解して最適化する方向性が考えられる．Robust PCA と類似の凸最適化になる．
現状の実装は cvxpy のソルバに依存している（今回は ECOS を採用）．当然，フレーム数や手法数が増えると最適化に多くの計算時間を要する．そこで ADMM が使えるように問題を再定式化し，高速化を図る方向性が考えられる．実際にそれは可能であって，手元の ADMM 実装を用いた数十秒の歌声データに対する簡易的な実験では，最適化の計算時間を削減しつつ（3割〜4割程度の削減; ハイパパラメタ依存），ほぼ同等のピッチ推定が得られることを確認している．