正弦波制約微分方程式（Sinusoidal Constraint Differential Equation; SCDE）によるピッチ推定法をPythonで実装した

はじめに

従来のピッチ推定手法として、FFT（高速フーリエ変換）や自己相関関数などが広く使われている。しかし、これらの手法には共通の課題があった。それは、「短い時間（例えば数十ミリ秒）の信号から、複数の音が同時に鳴っている場合のピッチを正確に推定するのが難しい」という点である。FFTは時間分解能と周波数分解能がトレードオフの関係にあり、短い信号では周波数分解能が低下し、密接した周波数成分を分離できないことがしばしばある。

本記事では、この課題に挑む1つのアプローチ「正弦波制約微分方程式（Sinusoidal Constraint Differential Equation; SCDE）」[1-3] について紹介する。またPython実装も併せて紹介する。

はじめに
正弦波と微分方程式の関係：SCDEの基礎
複数の独立した正弦波を推定するSCDE
実装
実験
- 実験結果
おわりに
参考文献

正弦波と微分方程式の関係：SCDEの基礎

SCDEの核となるのは、単純な正弦波の性質である。我々が知っている正弦波は、以下のような形で表される。

$\begin{align*} x(t) = A \sin(\omega t + \phi) \end{align*}$

ここで、 $A$ は振幅、 $t$ は時間、 $\omega$ は角周波数、 $\phi$ は初期位相である。この正弦波を2回時間で微分する。

$\begin{align*} \frac{d^{2}}{dt^{2}} x(t) = -A\omega^{2} \sin(\omega t + \phi) = -\omega^{2} x(t) \end{align*}$

右辺を左辺に移項して、以下の2階微分方程式を得る。

$\begin{align*} \frac{d^{2}}{dt^{2}} x(t) + \omega^{2} x(t) = 0 \end{align*}$

これが「正弦波制約微分方程式（SCDE）」の最も基本的な形となる。注目すべきは、この方程式が振幅 $A$ や初期位相 $\phi$ に依存しない点である。つまり、信号の形そのものから、その周波数 $\omega$ だけを抽出できる可能性を示唆している。

SCDEでは、この微分方程式が「できるだけ成り立つ」ように、その残差の2乗を最小化するような $\omega$ を求める。つまり、次式で与えられる目的関数を最小化する $\omega$ を求めることで、周波数の推定を行う。

$\begin{align*} J(\omega ) = \int_{\Gamma} \left[ \frac{d^{2}}{dt^{2}} x(t) + \omega ^{2} x(t) \right]^{2} dt \end{align*}$

ここで、 $\Gamma$ は定積分の範囲を表す。この目的関数を $\omega$ で微分して0と置くと、 $\Omega = \omega^{2}$ に関するシンプルな代数方程式が得られるので、その根を求めることで周波数を推定できる。

以下、具体的な代数方程式の形を示す。まず、 $J(\omega)$ を整理して、 $\Omega$ を含む式にする。目的関数が $\Omega$ の関数であることを明示するために $J(\Omega)$ と書いている。

$\begin{align*} J(\Omega) &= \int_{\Gamma} \left[ \left(\frac{d^{2}}{dt^{2}} x(t)\right)^{2} + 2 \Omega \left(\frac{d^{2}}{dt^{2}} x(t)\right) x(t) + \Omega^{2} x(t)^{2} \right] dt \\ &= \int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t)\right)^{2} dt + 2\Omega \int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t) \right) x(t) dt + \Omega^{2} \int_{\Gamma} x(t)^{2} dt \end{align*}$

続いて $\Omega$ で微分する。

$\begin{align*} \frac{J(\Omega)}{ d \Omega} = 2\int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t)\right) x(t) dt + 2\Omega \int_{\Gamma} x(t)^{2} dt \end{align*}$

ここで次式で定義される共分散を導入する。

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

最適な $\Omega$ を求めるための代数方程式は次式で与えられる。

$\begin{align*} \frac{J(\Omega)}{ d \Omega} = 2 S_{0, 0} \Omega + 2S_{1, 0} = 0 \end{align*}$

この方程式は $\Omega$ の1次方程式なので簡単に解けて、最終的に推定される周波数は

$\begin{align*} \hat{\omega} = \sqrt{\Omega} = \sqrt{-\frac{ S_{1, 0} } { S_{0, 0} }} \end{align*}$

となる。

複数の独立した正弦波を推定するSCDE

基本のSCDEは1つの正弦波に有効だが、実際の音には複数の周波数成分が含まれる。SCDEはこれを拡張し、複数の独立した正弦波の和にも適用できる。一般に、 $N$ 個の正弦波からなる場合、

$\begin{align*} x(t) &= A_1 \cos (\omega_1 t + \phi_1) + A_2 \cos (\omega_2 t + \phi_2) + \cdots + A_N \cos (\omega_N t + \phi_N) \\ &= \sum_{i = 1}^{N} A_i \cos(\omega_i t + \phi_i) \end{align*}$

と書ける。この信号は、それぞれの正弦波が満たす微分演算子を連ねた、より高次の微分方程式を満たす：

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

実際、上記の微分演算子により、 $A_1 \cos (\omega_1 t + \phi_1)$ から $A_N \cos (\omega_N t + \phi_N)$ までの全ての項が方程式を満たしている。左辺を展開して以下のような形に整理できる。

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

ここで $D^{k} = d^{k} / d t^{k}$ という記号を導入した。それぞれの $\alpha_n$ は $\omega_1^{2}, \ldots, \omega_N^{2}$ の基本対称式として書ける。逆に、 $\omega_n$ は $\alpha_0$ から $\alpha_N$ までの係数の情報を使って求めることができる（代数的に求められるとは限らないという意味を込めて、「情報」と表現した）。

1つの正弦波のときと同様に、上式ができるだけ成り立つような $\alpha_n$ を求めるべく、左辺を二乗したうえで定義される目的関数

$\begin{align*} J (\alpha_0, \alpha_1, \ldots, \alpha_N) = \int_{\Gamma} \left[ \sum_{n=0}^{N} \alpha_{n} D^{2(N-n)} x(t) \right]^{2} dt \end{align*}$

の最小化を考える。 $\alpha_1, \ldots, \alpha_N$ の極値条件からなる連立方程式を解くことで、各 $\alpha_n$ を求めることができる。 $J (\alpha_0, \ldots, \alpha_N)$ を $\alpha_1, \ldots, \alpha_N$ のそれぞれについて偏微分して0と置いた結果は

$\begin{align*} S_{N,N-1} &+ \alpha_1 S_{N-1,N-1} + \dots + \alpha_N S_{0, N-1} = 0 \\ S_{N, N-2} &+ \alpha_1 S_{N-1,N-2} + \dots + \alpha_N S_{0, N-2} = 0 \\ &\vdots \nonumber \\ S_{N, 0} &+ \alpha_1 S_{N-1, 0} \hspace{5.1mm} + \dots + \alpha_N S_{0,0}\hspace{5.3mm} = 0 \end{align*}$

となる。ただし、展開式より $\alpha_0 = 1$ と常に定数であるから、 $\alpha_0$ の極値条件は省略している。 $S_{i, j}$ の定義は既に述べた通りであるが、再掲する。

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

方程式を少し整理して

$\begin{align*} \textbf{B} \textbf{a} = - \textbf{b} \end{align*}$

とする。ここで $\textbf{a} = \left( \alpha_1, \alpha_2, \ldots, \alpha_N\right)^{\top}$ 、 $\textbf{b} = \left( S_{N,N-1}, S_{N,N-2} \ldots, S_{N, 0} \right)^{\top}$ 、

$\begin{align*} \textbf{B} = \begin{pmatrix} S_{N-1,N-1} & \cdots & S_{0,N-1} \\ \vdots & \ddots & \vdots \\ S_{N-1, 0} & \cdots & S_{0, 0} \end{pmatrix} \end{align*}$

とした。 $\textbf{B}$ は共分散行列に相当する。最適な $\alpha_1, \ldots, \alpha_N$ は次式で与えられる。

$\begin{align*} \textbf{a} = - \textbf{B}^{-1} \textbf{b} \end{align*}$

目的関数の係数が求まれば、これらは基本対称式であるから、 $\omega_1^{2}, \ldots, \omega_N^{2}$ は以下の $N$ 次方程式の解として推定される。

$\begin{align*} \Omega^{N} + (-1)^{1} \alpha_1 \Omega^{N-1} + (-1)^{2} \alpha_2 \Omega^{N-2} + \cdots + (-1)^{N} \alpha_N = 0 \end{align*}$

4次方程式までは代数的に解くことができるが、5次以上の方程式は代数的に解くことはできない。実装ではいずれの方程式も数値的に解いている。

実装

以下に置いた。Enjoy!

A demonstration of Frequency Estimation Based on SCDE. · GitHub

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

apply_bandpass_filter: バタワース特性のバンドパスフィルタを適用する関数である。
get_kth_even_derivative : 導関数 $D^{2i} x(t)$ を計算するための関数である。数値微分は中心差分によって実装した。
calc_covmat_mn: 共分散行列 $\textbf{B}$ の各成分を計算するための関数である。数値微分の度に信号長が短くなるので、高階微分どうしの積和（数値積分）を計算するために信号長を考慮してトリミングする必要がある。
solve_poly: 周波数推定のための $N$ 次方程式の求解を行うための関数である。
estimate_frequencies_scde: calc_cov_mnによる共分散行列の計算結果を用いて、係数 $\alpha_n$ らの連立方程式の求解を行う。続いてそれら係数に基づくsolve_polyによる $N$ 次方程式の求解を実施し、最終的な周波数の推定を行うための関数である。本実験では $N = 3$ としている。
main: 検証実験を実施するための関数である。

実験

周波数 440 Hz, 460 Hz, 480 Hzという、3つの正弦波の和で書ける信号の周波数を推定する実験を行った。標本化周波数は 44.1 kHz である。また信号長は 40 ms と 2 ms の2種類とした。40 ms に関しては、信号にホワイトノイズを加えたうえで推定させるパターンも試した。このときのノイズレベルは 10 dB、ノイズの周波数帯域は3つの周波数を覆う430 Hz〜490 Hz とした。ノイズの帯域制限はバタワース特性をもつバンドパスフィルタにより実現した。

実験結果

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

--- Test without noise ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies: [439.92795114 459.91768358 479.90645314] Hz
Estimation errors: [-0.07204886 -0.08231642 -0.09354686] Hz

--- Test with band-limited noise ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies with noise: [439.86240043 459.95433428 479.8400905 ] Hz
Estimation errors: [-0.13759957 -0.04566572 -0.1599095 ] Hz

--- Test with very short signal duration (2ms) ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies for 2ms signal: [440.91214252 460.90604637 480.98345977] Hz
Estimation errors: [0.91214252 0.90604637 0.98345977] Hz

"Test without noise" は3つの正弦波の和に対して、ノイズを加えずに推定した結果を示す。推定誤差はかなり小さく、ピッチ推定がうまくいっていることが分かる。

"Test with band-limited noise" は帯域制限したホワイトノイズを加えて推定した結果を示す。ノイズレベルが 10 dB であっても推定結果は良好である。

"Test with very short signal duration" は 2 ms という、とても短い信号に対して推定した結果を示す。誤差は 40 ms のときと比較して増大しているが、誤差は 1 Hz 程度に留まった。

おわりに

本記事ではSCDEによるピッチ推定法とその実装を紹介した。簡単な検証実験の結果、ピッチ推定はうまくいったようである。

音声など倍音構造を持つ信号に対してもSCDEは適用可能であるが、各周波数は独立に推定され、周波数間の調波構造は考慮されていない。SCDEに調波構造の制約を課すことで、推定精度のさらなる向上が見込める。またSCDEはノイズに対する頑健性が低いという問題も指摘されている。文献 [2] では、SCDEに調波構造の制約を導入した、有限次数調波拘束微分方程式（Finite-Order Harmonic Constraint Differential Equation; FOHCDE）に基づく手法が提案されている。すでに FOHCDE 自体の実装は済んでおり、記事を後日公開予定である。

追記：FOHCDEの記事を書いた。 tam5917.hatenablog.com

参考文献

[1] K. Yamada, Y. Masuyama, Y. Wakabayashi and N. Ono, "Simultaneous Frequency Estimation for Three or More Sinusoids Based on Sinusoidal Constraint Differential Equation," 2022 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), pp. 975-978, 2022. https://ieeexplore.ieee.org/document/9980228

[2] 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

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