有限次数調波拘束微分方程式（FOHCDE）に基づく多重ピッチ推定法のPython実装と再現実験

はじめに

従来のピッチ推定手法、特にDFTをベースとする手法では、時間分解能と周波数分解能のトレードオフが本質的な課題であった。短時間の信号から近接した複数の周波数を分離・推定することは依然として困難なタスクである。この課題に対し、山田氏らによる技術報告 [1] および修士論文 [2] では、微分方程式に基づくアプローチを多重音源かつ調波構造を持つ信号へと拡張した「有限次数調波拘束微分方程式（M-FOHCDE）」に基づく手法が提案されている。本手法は、モデルベースのアプローチにより、短時間信号に対しても高精度な推定を可能にすると述べられている。

本記事では、上記論文で提案されているM-FOHCDE法をPythonで実装し、その性能を検証するために論文中で報告されているいくつかの再現実験を行った。

はじめに
M-FOHCDE法の概要
- 多重音源モデルの定式化
- 交互最適化による目的関数の最小化
  - 交互最適化の手順
実装
再現実験
おわりに
参考文献

M-FOHCDE法の概要

多重音源モデルの定式化

信号に、それぞれが独自の基本周波数 $\omega_1, \omega_2$ ... と調波構造を持つ、 $N$ 個の音源が含まれていると仮定する。信号は次式で表される。

$\begin{align*} x(t) = \sum_{n=1}^{N} \sum_{m=1}^{M_n} A_{n, m} \cos (m\omega_n t + \phi_{n, m}) \end{align*}$

$A_{n, m}$ と $\phi_{n, m}$ はそれぞれ、音源 $n$ に関する第 $m$ 調波成分の振幅と初期位相を表す。 $M_n$ は音源 $n$ に含まれる調波成分の数である。

各音源 $n$ に、単体のFOHCDE作用素 $L_n(p)$ が対応している。

$\begin{align*} L_n(p) = \prod_{m=1}^{M_n} (p^2 + m^2\omega_n^{2}) \end{align*}$

ここで $p = d/dt$ である。信号全体は、これらの作用素をすべて乗算した、単一の高階微分方程式を満たすと考える。

$\begin{align*} L_{\mathrm{total}}(p) [x(t)] = (L₁(p) * L₂(p) * \ldots * L_N(p)) x(t) = 0 \end{align*}$

この方程式ができるだけ成り立つように、すなわち残差が最小になるように、目的関数 $J(\omega_{1}, \omega_{2}, \ldots, \omega_{N})$ を定義し、これを最小化する $\omega_n$ の組を求めることが目標である。

目的関数は次式で定義される。

$\begin{align*} J(\omega_{1}, \omega_{2}, \ldots, \omega_{N}) &= \int_{\Gamma} \left( L_{\mathrm{total}}(p) [x(t)] \right)^2 dt \end{align*}$

交互最適化による目的関数の最小化

複数の音源を扱う問題は、最終的に多変数 $(\omega_{1}, \omega_{2}, \ldots, \omega_{N})$ からなる目的関数 $J$ を最小化する問題として定式化される。

$\begin{align*} \mathrm{Minimize} : J(\omega_{1}, \omega_{2}, \ldots, \omega_{N}) \end{align*}$

しかし、目的関数は非常に複雑な非線形関数であり、全ての変数を同時に動かしながら最小値を見出すのは困難である。そこで本論文では交互最適化によるアプローチが採用されている。

交互最適化の手順

1．全ての周波数 $\omega_{1}, \omega_{2}, \ldots, \omega_{N}$ を、何らかの初期値に固定する

2．その状態で、 $\omega_{1}$ だけを動かして $J$ が最小になる点を探す

$\begin{align*} \mathrm{Minimize} \; \mathrm{w.r.t.}\; \omega_{1} : J(\omega_{1}, \mathrm{const}, \mathrm{const}, \ldots) \end{align*}$

3．次に、更新された $\omega_{1}$ と、 $\omega_{3}$ 以降を固定し、 $\omega_{2}$ だけを動かして $J$ が最小になる点を探す

$\begin{align*} \mathrm{Minimize} \; \mathrm{w.r.t.}\; \omega_{2} : J(\mathrm{const}, \omega_{2}, \mathrm{const}, \ldots) \end{align*}$

4．このプロセスを、 $\omega_{3}, \omega_{4}, \ldots, \omega_{N}$ と順番に、全ての周波数について実行する

5．1から4までのサイクルを、各周波数の値がほとんど変化しなくなる（収束する）まで繰り返す

各ステップでは1変数のみの最適化（1次元の探索問題）なので、多変数の同時最適化よりも問題はずっと解きやすい。

次節では、この交互最適化の各ステップ、特に「1変数最適化」をどのように実装したかについて、具体的に説明する。

実装

以下に置いた。Enjoy!

m-fohcde_demo.py · GitHub

前節では、多変数の最適化問題を1変数の最適化問題の繰り返しに帰着させる「交互最適化」について述べたが、これをどのように実装したかについて説明する。

交互最適化のアルゴリズム全体を動かすための estimate_multiple_pitch 関数を用意した。主に以下の処理からなる：

共分散行列の事前計算：最適化のループを開始する前に、計算コストの高い共分散行列を事前に一度だけ計算しておく
反復ループの実行：各反復において、全ての音源 $n = 1, 2, \ldots, N$ について、その時点での最適な周波数 $\omega_n$ を計算し、値を更新する
収束判定：反復の前後で周波数の更新量が十分に小さくなった場合、解が収束したとみなして計算を打ち切る

def estimate_multiple_pitch(...):    
    # 1. 共分散行列の事前計算
    S_matrix = # ... S_ij を計算する処理
    # 初期化
    omegas_sq = # ... 初期値の準備
    # 2. 反復ループ
    for i in range(max_iter):
        prev_omegas_sq = np.copy(omegas_sq)
        for n in range(n_sources):
            # n番目の周波数以外を固定し、1変数最適化を実行
            new_omega_sq = solve_for_one_omega(...)
            # 値を更新
            if new_omega_sq is not None:
                omegas_sq[n] = new_omega_sq
        # 3. 収束判定
        if np.linalg.norm(omegas_sq - prev_omegas_sq) < tolerance:
            break
    return # ... 最終的な周波数をHzで返す

反復ループ内では、1次元の最適化を行う solve_for_one_omega 関数が呼び出される。この関数は、他の $N-1$ 個の周波数が固定された状態で、残る1つの周波数 $ω_{n}$ に関する目的関数 $J(ω_{n})$ を最小化する。しかし $J(ω_{n})$ は $ω_{n}$ の複雑な多項式となるため、解析的に係数を導出・計算するのはコストが高い。そこで、 $J(ω_{n})$ の点をサンプリングし、それらサンプル点から多項式をフィッティングすることで、未知の関数を「復元」する方針をとった。

以下、solve_for_one_omega 関数における処理の概要を示す。

1． $J(ω_{n})$ の値のサンプリング：最適化対象の変数 $Ω_{n} = \omega_{n}^2 > 0$ を複数点サンプリングする。サンプリングの範囲は物理的にあり得る範囲で広めに取る。それぞれのサンプル値における目的関数の値を具体的に計算することで、「点群データ」が手に入る。

2．多項式フィッティングによる関数の復元：得られた点群データを最もよく表現する多項式を、Numpyの polyfit 関数を用いてフィッティングする。ここでいう「最もよく表現する」とは、各サンプル点における実際の値と、多項式による予測値との差の二乗和（残差二乗和）を最小化することを意味する。

3．導関数の根の計算：復元した多項式の係数から導関数 $dJ/d\Omega_{n}$ の係数を求める。そしてNumpyの roots 関数を用いて根（一般に複数個）を求める。

4．最適解の選択：上のステップで得られた根は目的関数を最小化する候補たちであり、妥当なものを選択する必要がある。極端に値の大きすぎる根や複素根、負値の実根を捨てたうえで、正の値の実根たちで絞り込む。その中で目的関数の値を最も小さくするものを最適解として出力する。

ちなみにサンプリングした点群データから未知の関数を復元する、というやり方はかつて実装したことがあった（当時はガウス過程状態空間モデルだったが）。

再現実験

実験その1（文献 [1], 4.2節の再現）

提案手法 (M-FOHCDE) が、DFTをベースとする従来手法と比較して、どの程度の優位性を持つかを評価する。

実験方法

基本周波数が200Hzと440Hzの2つの音源（それぞれ2次倍音までを含む）を合成した信号を生成する。信号長を1msから50msまで変化させ、各信号長において以下の4手法で200Hzの基本周波数を推定する。

M-FOHCDE (提案手法)
SCDE (調波構造を考慮しない微分方程式ベースの手法)
Spectral Reassignment (高度なDFTベースの手法)
DFT + Parabolic Interpolation (基本的なDFTベースの手法)

ちなみにDFTベースの2手法は文献 [2] に推定アルゴリズムが書かれており、実装のうえで参考になる。

実験結果

図1: 各手法における信号長と推定誤差の関係；横軸は入力信号の長さ（ミリ秒）、縦軸は200Hzの基本周波数に対する推定誤差（Hz）を示す

図1に各手法の推定結果を示す。論文で報告された傾向と、よく一致している結果が得られた。提案手法であるM-FOHCDE（およびSCDE）による推定は、信号長にほとんど依存せず、2msという短い信号においても誤差を低く保つことができた。一方、DFTをベースとする2つの従来法は、信号長が短い領域（約12ms以下）では周波数分解能が低くなったことにより、推定誤差が非常に大きくなった。

実験その2（文献 [1], 4.3節の再現）

ここでは提案手法の性能が、入力信号に含まれる周波数の組み合わせによってどのように変化するかを評価し、その適用限界を明らかにすることが目的である。

実験方法

一方の音源の基本周波数を440Hzに固定し、もう一方の基本周波数を1Hzから1000Hzまで変化させる。各音源は2次倍音までを含む。信号長は10msに固定し、各周波数の組み合わせについて提案手法M-FOHCDEを適用する。

実験結果

図2に実験結果を示す。ほとんどの周波数の組み合わせにおいて推定値と真の値が一致しており、高い推定精度を持つことを示した。論文中では、特定の周波数（220, 440, 880Hz; 一方の周波数がもう一方の倍音と重複する条件）において誤差の増大が報告されていた。しかし、今回の再現実験では、そのような誤差は観測されなかった。

今回の実装において、真値に近い初期値を与えて最適化を開始したことは、推定誤差が観測されなかった要因の一つと考えられる。また、 solve_for_one_omega 関数が安定して動作したことにより、アルゴリズムが局所解に陥ることなく、大域的最適解を見つけ出すことができたためとも考えられる。ノイズのない理想的な信号を用いたことも、今回の結果に寄与した可能性がある。

実験その3（文献 [1], 4.4節の再現）

ここでは、周波数が時々刻々と変化する信号に対して、提案手法がどの程度追従できるかを評価する。

実験方法

一方の音源の基本周波数を正弦波状に変動させ（ビブラートを模擬）、もう一方を固定周波数とした信号を生成する。この信号に対し、5msの分析窓を1msずつ移動させながら、各フレームでM-FOHCDEを適用し、2つの基本周波数を推定する。

実験結果

図3に実験結果を示す。本図が示すように、推定された周波数（青い点群）は、真の周波数の時間変動（赤い破線）によく追従している。M-FOHCDEは高い時間分解能を持ち、音楽信号におけるビブラートのような微細な周波数変化の分析にも有効であることが確認できた。ただし2つの周波数が近接・交差する領域では、推定値にばらつきが見られる。これは手法の特性を反映したものと考えられ、適用限界を示している。

おわりに

本記事では、山田氏らにより提案された、有限次数調波拘束微分方程式（M-FOHCDE）に基づく多重ピッチ推定法をPythonで実装し、その性能を再現実験によって検証した。

実験の結果、提案手法は特に短時間信号の分析において、DFTベースの従来手法に対する優位性を持つことが確認された。また、周波数が近接・重複する条件下では推定精度が低下する挙動も見られ、手法の特性と適用限界についても把握できた。

ちなみに今回の実装は、論文で主に扱われている「2音源・2次倍音」という特定のケースを前提としていない。任意の音源数 $N$ と、各音源の任意の倍音次数 $M_n$ の組み合わせに対応可能な設計となっている。

参考文献

[1] 山田健太, 升山義紀, 山岡洸瑛, 植野夏樹, 小野順貴, "微分方程式に基づく有限次数調波信号の多重ピッチ推定," 電子情報通信学会技術報告, vol. 123, no. 401, pp. 315-320, 2024.

[2] 山田健太, "微分方程式に基づく有限次数調波信号の短時間基本周波数推定", 東京都立大学修士論文, 2024. Link