有限次数調波拘束微分方程式（Finite-Order Harmonic Constraint Differential Equation; FOHCDE）に基づく基本周波数推定とPython実装

はじめに

本記事は前回記事の続編である。前回は微分方程式の性質を利用して音の周波数を推定する手法「SCDE（Sinusoidal Constraint Differential Equation）」を紹介した。SCDEは、FFTベースの手法では困難であったわずか数ミリ秒という極めて短い信号からも、複数の正弦波の周波数を高精度に推定できる。

本記事では、より現実的な環境で高精度なピッチ推定を可能にする「有限次数調波拘束微分方程式（Finite-Order Harmonic Constraint Differential Equation; FOHCDE）」[1, 2]を紹介する。

はじめに
基本周波数推定における調波構造
SCDEからFOHCDEへ：調波制約の導入
- 係数の具体例
FOHCDEに基づく基本周波数の推定
- 代数方程式の導出
実装
実験
- 実験結果
おわりに
参考文献

基本周波数推定における調波構造

我々が普段耳にする人間の声や多くの楽器の音は、純粋な単一の正弦波ではない。例えば、ギターのA4（440Hz）の音は、440Hzの音だけでなく、その倍である880Hz（2倍音）、1320Hz（3倍音）、1760Hz（4倍音）といった、基本周波数の整数倍の周波数成分を同時に含む。この基音と倍音からなる構成を「調波構造」と呼ぶ。

もし音が3つの倍音（基音、2倍音、3倍音）からなるとしたら、従来のSCDEでは3つの独立した周波数パラメータを推定する必要があった。しかし、調波構造を仮定すれば、推定すべきは「基本周波数」のみになる。これにより、ノイズの中から抽出するべき情報が絞られ、推定の信頼性（頑健性）の向上が期待できる。

以上が調波構造が周波数推定において重要な理由である。

SCDEからFOHCDEへ：調波制約の導入

FOHCDEは、SCDEの原理をこの調波構造に合わせて拡張する。SCDEでは、各正弦波が満たす微分演算子 $\left( \frac{d^{2}}{dt^{2}} + \omega_{k}^{2} \right)$ を連ねることで高次微分方程式を構築した。FOHCDEでは、これを「単一の基本周波数 $\omega_1$ と、その $M$ 次倍音まで」という形で表現する。

したがって、FOHCDEの核となる微分方程式は以下のようになる。

$\begin{align*} \left(\frac{d^{2}}{dt^{2}} + (1\omega_1)^{2}\right)\left(\frac{d^{2}}{dt^{2}} + (2\omega_1)^{2}\right) \cdots \left(\frac{d^{2}}{dt^{2}} + (M\omega_1)^{2}\right) x(t) = 0 \end{align*}$

この高次微分方程式を展開すると、

$\begin{align*} D^{2M} x(t) + \alpha_{1} \Omega_{1} D^{2(M - 1)} x(t) + \cdots + \alpha_{M-1} \Omega_{1}^{M-1} D^{2} x(t) + \alpha_{M} \Omega_{1}^{M} x(t) = 0 \end{align*}$

となる。ここで $D^{k} = d^{k} / d t^{k}$ である。また $\Omega_{1} = \omega_{1}^{2}$ と置いた。係数 $\alpha_{m}$ は次式に示す $x$ の多項式において、各項の係数を左辺と右辺で比較して求めることができる。

$\begin{align*} \sum_{m=0}^{M} \alpha_{m} x^{m} = \prod_{m=0}^{M} (1 + m^{2} x) \end{align*}$

係数 $\alpha_{m}$ の具体例

$M=2$ としてみよう。

$\begin{align*} \prod_{m=0}^{M} (1 + m^{2} x) &= (1 + 0^{2} x) (1 + 1^{2} x) ( 1 + 2^{2} x) \\ &= 1 + 5x + 4x^{2} \end{align*}$

となるから、 $\alpha_{0} = 1$ , $\alpha_{1} = 5$ , $\alpha_{0} = 4$ である。実際の微分方程式を見てみると、

$\begin{align*} \left(\frac{d^{2}}{dt^{2}} + \omega_{1}^{2}\right)\left(\frac{d^{2}}{dt^{2}} + 4\omega_2^{2}\right) x(t) = 0 \end{align*}$

から左辺を整理して

$\begin{align*} \left( \frac{d^{4}}{dt^{4}} + 5 \omega_{1}^{2} \frac{d^{2}}{dt^{2}} + 4\omega_{1}^{4} \right)x(t) & = 0 \\ \Leftrightarrow D^{4} x(t) + 5 \Omega_{1} D^{2} x(t) + 4\Omega_{1}^{2} x(t) & = 0 \end{align*}$

を得る。多項式の展開から求めた $\alpha_{0}$ , $\alpha_{1}$ , $\alpha_{2}$ が、微分方程式の展開から得られる係数とうまく整合していることが分かる。

FOHCDEに基づく基本周波数の推定

件の微分方程式が「できるだけ成り立つ」ような $\omega_1$ を求めるため、その左辺を二乗して積分した目的関数をSCDEと同様に定義する。

$\begin{align*} J(\Omega_1) = \int_{\Gamma} \left[ \sum_{m=0}^{M}\alpha_m \Omega_{1}^{m} D^{2(M-m)} x(t) \right]^{2} dt \end{align*}$

この目的関数を最小化するような $\Omega_1$ を考えるため、 $J(\Omega_1)$ を $\Omega_1$ で微分して結果を0と置く。

$\begin{align*} \frac{d J (\Omega_1)}{d\Omega_1} = 0 \end{align*}$

上記は $\Omega_1$ に関する代数方程式であり、その根を求めることで直接 $\Omega_1 = \omega_{1}^{2}$ を推定し、ひいては基本周波数 $\omega_{1}$ を得ることができる。

代数方程式の導出

$J(\Omega_1)$ を整理するうえで、係数 $\alpha_m$ や $\Omega_{1}^{m}$ は積分の外に出してよいので

$\begin{align*} J(\Omega_1) &= \int_{\Gamma} \left[ \sum_{i=0}^{M}\alpha_i \Omega_{1}^{i} D^{2(M-i)} x(t) \right] \left[ \sum_{j=0}^{M}\alpha_j \Omega_{1}^{j} D^{2(M-j)} x(t) \right] dt \\ &=\sum_{i=0}^{M} \sum_{j=0}^{M} \alpha_i \alpha_j \Omega_{1}^{i+j} \int_{\Gamma} \left(D^{2(M-i)} x(t)\right) \left(D^{2(M-j)} x(t) \right) dt \\ &= \sum_{i=0}^{M} \sum_{j=0}^{M} \alpha_i \alpha_j S_{M-i, M-j} \Omega_{1}^{i+j} \end{align*}$

となる。ただし $S_{i, j}$ を次式で定義した。

$\begin{align*} S_{i, j} = \int_{\Gamma} \left[ \frac{d^{2i}}{dt^{2i}} x(t) \frac{d^{2j}}{dt^{2j}} x(t) \right] dt \end{align*}$

$J(\Omega_1)$ は $\Omega_{1}$ に関する $2M$ 次の多項式であり、各項の係数を $\beta_m$ と置く：

$\begin{align*} J(\Omega_1) = \sum_{m=0}^{2M} \beta_m \Omega_{1}^{m} = \sum_{i=0}^{M} \sum_{j=0}^{M} \alpha_i \alpha_j S_{M-i, M-j} \Omega_{1}^{i+j} \end{align*}$

ゆえに $\Omega_{1}$ に関する $J(\Omega_1)$ の極値条件として、

$\begin{align*} \frac{d J (\Omega_1)}{d\Omega_1} = \sum_{m=1}^{2M-1} m \beta_m \Omega_{1}^{m-1} = 0 \end{align*}$

を得る。これは $\Omega_{1}$ に関する $2M - 1$ 次の代数方程式である。

実装

以下に置いた。Enjoy!

fohcde_demo.py · GitHub

実装した関数の概要を書いておく。

apply_bandpass_filter: バタワース特性のバンドパスフィルタを適用する関数である。
get_kth_even_derivative : 偶数次数の導関数 $D^{2i} x(t)$ を計算するための関数である。数値微分は中心差分によって実装した。
calc_covmat_mn: 高次導関数どうしの共分散行列の成分を計算するための関数である。
calc_covmat:calc_covmat_mnを呼び出して共分散行列全体を計算するための関数である。
calc_alpha_coeff: 微分方程式に現れる係数 $\alpha_m$ を計算するための関数である。

$\begin{align*} \sum_{m=0}^{M} \alpha_{m} x^{m} = \prod_{m=0}^{M} (1 + m^{2} x) \end{align*}$

という $x$ の多項式を0とおいた方程式は、右辺の形から $x = - 1/m^{2}\; (m = 1, 2, \ldots, M)$ を解に持つ。ゆえにその解を元に、numpy.polynomial.polynomial の polyfromroots 関数を使って多項式の係数が計算できる。このとき、最高次の項の係数は1に正規化されているので、係数の値を元に戻す処理を入れている。つまり次式における右辺の $\prod_{m=1}^{M} m^{2}$ を乗じる処理に相当する。

$\begin{align*} \prod_{m=0}^{M} (1 + m^{2} x) = \prod_{m=1}^{M} m^{2} \cdot \prod_{m=1}^{M} \left( \frac{1}{m^{2}} + x\right) \end{align*}$

calc_beta_coeff: 目的関数 $J(\Omega_1)$ を $\Omega_1$ の多項式とみなしたときに現れる、各項の係数 $\beta_m$ を計算するための関数である。

$\begin{align*} \sum_{m=0}^{2M} \beta_m \Omega_{1}^{m} = \sum_{i=0}^{M} \sum_{j=0}^{M} \alpha_i \alpha_j S_{M-i, M-j} \Omega_{1}^{i+j} \end{align*}$

という関係式から、 $\Omega_{1}$ のべき乗が $m = i + j$ となる項を右辺から集めていけば、係数は容易に計算できる。

calc_poly_coeff: 目的関数の極値条件 $\frac{d J (\Omega_1)}{d\Omega_1} =0$ から得られる $\Omega_{1}$ の代数方程式に関して、各項の係数を計算するための関数である。

$\begin{align*} \sum_{m=1}^{2M-1} m \beta_m \Omega_{1}^{m-1} = 0 \end{align*}$

という方程式であるから、計算済の $\beta_m$ を使い、係数は $m \beta_m\; (m = 1, 2, \ldots, M)$ と計算できる。

solve_poly_getf0: $\Omega_{1}$ の代数方程式を解き、得られた $\Omega_{1}$ の解から最終的な基本周波数を推定するための関数である。代数方程式の求解には numpy.polynomial.polynomial の polyroots 関数を使って実装した。基本角周波数は平方根 $\omega_1 = \sqrt{\Omega_1}$ から計算でき、 $2\pi$ で割って Hz 単位の基本周波数にしている。

実験

基本周波数 200 Hz で3次倍音まで含む信号を対象に、基本周波数の推定実験を行った。標本化周波数は 44.1 kHz，信号長は 40 ms と 2 ms である。信号に帯域制限されたホワイトノイズを加えたうえで推定させるパターンも試した。このときのノイズレベルは 10 dB、帯域は 180 Hz 〜 660 Hz である。

本来はSCDEとの比較も行うべきだが、今回はサボった。本記事読者の追実験に期待する。

実験結果

スクリプトを実行した結果は以下の通りである。

--- FOHCDE Single Source Demo: f0=200.0 Hz with M=3 harmonics ---
True fundamental frequency: 200.0 Hz
Signal contains harmonics up to 3x f0 (600.0 Hz)

--- Test without noise (40ms signal) ---
Estimated fundamental frequency: 199.9403 Hz
Estimation error: 0.0597 Hz

--- Test with band-limited noise (40ms signal, SNR 10dB) ---
Estimated fundamental frequency with noise: 200.4510 Hz
Estimation error: 0.4510 Hz

--- Test with very short signal duration (2ms) ---
Estimated fundamental frequency for 2ms signal: 200.3918 Hz
Estimation error:  0.3918 Hz

--- Test with band-limited noise (2ms signal, SNR 10dB) ---
Estimated fundamental frequency with noise: 193.6471 Hz
Estimation error: 6.3529 Hz

40 ms のノイズなし信号については、基本周波数は精度よく推定された。広い帯域でノイズが加わっても、推定は概ねうまくいっているようである。また 2 ms の信号についても、ノイズがない場合は良好な推定値が得られたが、ノイズが加わった場合はさすがに推定精度が低下するようである。ノイズのランダム性を考慮すると、ノイズありの場合は複数回の推定を行った際の統計量（平均と標準偏差）を示すべきだろうが、それもサボった。

おわりに

本記事では調波構造を考慮した SCDE である FOHCDE に基づく基本周波数推定手法を紹介した。検証実験の結果、推定には概ね成功したといえるだろう。信号が短時間であっても高精度に推定できることが手法のウリだが、調波構造の制約を導入することで、精度改善が期待できる。

文献 [2] にはFOHCDEに基づく多重ピッチ推定手法が新たに提案されており、今後の記事にて実装を紹介予定である。ちなみにFOHCDEに関して、文献 [2] に掲載されている一部の式が誤っており（添字）、デバッグに少々時間を費やした。

参考文献

[1] K. Yamada, Y. Masuyama, K. Yamaoka and N. Ono, "Fundamental Frequency Estimation Based on Finite-Order Harmonic Constraint Differential Equation," 2023 Asia Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), pp. 868-872, 2023. https://ieeexplore.ieee.org/document/10317393

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