SCDEに基づく単一減衰正弦波のパラメタ推定とPython実装

はじめに

これまでの記事では、正弦波制約微分方程式（SCDE）に基づく周波数推定法と、その拡張を紹介した。特に前回記事では、重み付き最小二乗法（WLS）を導入することで、ノイズ耐性が向上することを確認した。

これらのモデルは、信号が時間的に変化しない定常的な正弦波であることを暗黙の前提としていた。しかし、ピアノやギターといった楽器の音をはじめとする現実世界の音は、発音された瞬間から時間と共に減衰する。波形のエンベロープ（音の響きの時間変化）は音色を特徴づける重要な要素であるが、これらのモデルでは捉えきれていない。

安藤氏による論文 [1] では、SCDEのフレームワークを減衰振動へと拡張する方法が示唆されていた。そこで、本記事では単一の減衰正弦波の周波数と減衰係数を同時に推定するアルゴリズムを考察する。アルゴリズムはPythonで実装し、簡単な検証実験の結果を報告する。

はじめに
減衰振動SCDE
- 減衰振動の微分方程式
- 複素角周波数への拡張
重み付き最小二乗法（WLS）の適用
実装
実験
- 実験結果
- おまけ
考察
おわりに
参考文献

減衰振動SCDE

減衰振動の微分方程式

減衰振動を表す信号は $x(t) = Ae^{-\epsilon t}\cos(at+\phi)$ と書ける。 $A$ は振幅、 $a$ は角周波数であり、 $\epsilon$ は減衰係数である。減衰係数の値が大きいほど、音は速く消えていく。この信号が満たす2階の微分方程式は、次式で与えられる。

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

この方程式には1階の微分項が含まれており、これまでのSCDEと大きく異なる点である。

減衰振動の特性方程式とその解を示しておく。 $\lambda = d/dt$ 、 $\omega_0^{2} = a^{2} + \epsilon^{2}$ とおいて、

$\begin{align*} \lambda ^2 + 2 \epsilon \lambda + \omega_0^{2} = 0 \end{align*}$

が特性方程式である。この方程式の根は $\lambda = -\epsilon \pm j \sqrt{ \omega_0^{2} - \epsilon^{2}}$ により与えられる。2つの基本解は $\exp (\lambda t)$ の形であるから、

$\begin{align*} x_1 (t) &= \exp (-\epsilon t) \exp (j \omega_d t) \\ x_2 (t) &= \exp (-\epsilon t) \exp (-j \omega_d t) \\ \end{align*}$

となる。ここで $\exp(-\epsilon t)$ が減衰を表し、 $\exp(\pm j \omega_d t)$ が振動を表す。 $\omega_d = \sqrt{ \omega_0^{2} - \epsilon^{2}}$ が実際に観測される角周波数であり、推定対象となる。つまり $a^{2} = \omega_{d}^{2}$ である。

複素角周波数への拡張

前回記事にて、境界値項をゼロにする窓関数を導入することで、SCDEは以下の方程式に帰着することを説明した。

$\begin{align*} k_{\omega} - 2j \omega h_{\omega} - \omega^{2}g_{\omega} + a^{2} g_{\omega} = 0 \end{align*}$

となる。ここで

$\begin{align*} g_{\omega} &= \int_{-T/2}^{T/2} x(t) p(t) \exp(-j\omega t)\: dt \\ h_{\omega} &= \int_{-T/2}^{T/2} x(t) p^{\prime}(t) \exp(-j\omega t)\: dt \\ k_{\omega} &= \int_{-T/2}^{T/2} x(t) p^{\prime\prime}(t) \exp(-j\omega t)\: dt \end{align*}$

である。減衰振動の微分方程式に対しても、同様の荷重積分と部分積分を適用できる。最終的に得られる代数方程式は

$\begin{align*} k_{\omega} - (2j \omega + 2\epsilon) h_{\omega} - \omega^{2} g_{\omega} + ( \omega_{0}^{2} + j 2\omega \epsilon) g_{\omega} = 0 \end{align*}$

となる。これを「複素角周波数の2乗」 $\tilde{a}^{2}$ を導入して

$\begin{align*} k_{\omega} - 2j \omega h_{\omega} - \omega^{2} g_{\omega} + \tilde{a}^{2} g_{\omega} = 0 \end{align*}$

とみなすならば、

$\begin{align*} \tilde{a}^{2} g_{\omega} &= 2 \epsilon h_{\omega}+ ( \omega_{0}^{2} + j 2\omega \epsilon) g_{\omega} \\ \tilde{a}^{2} &= 2 \epsilon h_{\omega} / g_{\omega} + \omega_{0}^{2} + j 2\omega \epsilon \\ & = 2 \epsilon h_{\omega} / g_{\omega} + (a^{2} + \epsilon^{2}) + j 2\omega \epsilon \end{align*}$

となる。ここで、 $2 \epsilon h_{\omega} / g_{\omega}$ が $\omega_{0}^{2}$ に比べて十分小さく、無視できるという近似を入れる。すると、

$\begin{align*} \tilde{a}^{2} \approx (a^{2} + \epsilon^{2}) + j 2\omega \epsilon \end{align*}$

となる。形式的には複素角周波数に置き換わるのみで、これまでの定常SCDEと全く同じ代数方程式を解く問題に帰着する。代数方程式を解いて $\tilde{a}^{2}$ が求まれば、減衰係数と角周波数はそれぞれ

$\begin{align*} \epsilon &= \frac{\Im\{\tilde{a}^{2} \}} {2 \omega} \\ a &= \sqrt{ \Re\{\tilde{a}^{2} \} - \epsilon^{2} } \end{align*}$

として推定できる。

重み付き最小二乗法（WLS）の適用

減衰振動を記述する微分方程式は、未知の複素周波数 $\tilde{a}^{2}$ を解くための、以下の代数方程式に帰着することを示した。

$\begin{align*} k_{\omega} - 2j \omega h_{\omega} - \omega^{2} g_{\omega} + \tilde{a}^{2} g_{\omega} = 0 \end{align*}$

この方程式は、評価用の周波数 $\omega$ を選ぶごとに一本得られる。それゆえ、 $\tilde{a}^{2}$ という唯一の未知数を解くためには、原理的には、どこか一つの評価用の周波数 $\omega$ を選んで、以下の方程式を解けば十分である。

$\begin{align*} \tilde{a}^{2} = \omega^{2} - \frac{k_{\omega} - 2 j\omega h_{\omega}}{g_{\omega}} \end{align*}$

しかしながら、現実の信号を扱う上では、ただ1つの方程式だけでは不十分である。現実の信号には、必ずノイズが含まれる。もし、選択したひとつの評価用周波数が、たまたまノイズの影響を強く受けていたりしたら、観測量 $g_{\omega}$ , $h_{\omega}$ , $k_{\omega}$ に含まれる誤差も大きくなる。結果的に $\tilde{a}^{2}$ の推定値も誤差を大きく含むものになってしまう。

そこで、前回記事にもあるように、ノイズに対する推定の頑健性を高めるため、複数の周波数ビンからの情報を統合する方針を取る。具体的には信頼度の高い周波数ビンを重視する重み付き最小二乗法（WLS）を適用する。

具体的には、以下の手順で $\tilde{a}^{2}$ を推定する：

1．スペクトルのピーク周辺から、信頼度の高い周波数ビン $\omega_k$ $(k=1, \ldots, K)$ を選択する

2．各ビンについて、方程式 $g_{\omega} \tilde{a}^{2} \approx \omega^{2} g_{\omega} + 2j \omega h_{\omega} - k_{\omega}$ を立てる。これらは $K$ 本の連立方程式 $\mathbf{A} x \approx \mathbf{b}$ と見なせる（ただし $x = \tilde{a}^{2}$ ）

3．各方程式 $k$ の信頼度を表す重み $w_k$ を導入し、重み付きの正規方程式 $\mathbf{A}^{*} \mathbf{W} \mathbf{A} x = \mathbf{A}^{*} \mathbf{W} \mathbf{b}$ を解く。ここで $\mathbf{A}^{*}$ は随伴行列（エルミート転置）であり、 $\mathbf{W} = \mathrm{diag}(w_1, w_2, \ldots, w_K)$ は重み行列である。

4．この方程式の解として、所望の複素周波数 $\tilde{a}^{2}$ が推定される。

実装

以下に置いた。Enjoy!

scde_demo_wls_dumped.py · GitHub

実装のポイントとなる関数たちに説明を付す。

calculate_observables：観測量 $g_{\omega}$ , $h_{\omega}$ , $k_{\omega}$ を計算
solve_complex_a_sq_wls：WLSによる複素周波数の推定
separate_params_from_complex_a_sq：推定された複素周波数から角周波数 $a$ と減衰係数 $\epsilon$ を抽出
estimate_damped_sinusoid：上記関数たちを呼び出して推定を遂行

実験

標本化周波数 44100 Hz、周波数 440 Hz 、減衰係数 10.0 の減衰信号を対象とする。信号長は 10 ms, 40 ms, 100 ms の3パターンを試す。また、推定はホワイトノイズを加えないオリジナルの信号と、ホワイトノイズを加えた信号に対して行う。ホワイトノイズには 0.05 というスケール定数を掛けている。

実験結果

信号長 10 ms の結果を示す。周波数推定は安定しているが、ノイズありの信号では、減衰係数の推定誤差が増えている。

--- Damped Sinusoid Estimation (original) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 440.00 Hz, epsilon = 10.22

--- Damped Sinusoid Estimation (noise) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 440.38 Hz, epsilon = 16.29

続いて、信号長 40 ms の結果を示す。ノイズの有無によらず、周波数と減衰係数がほぼ同じ誤差レベルで推定できたことが分かる。

--- Damped Sinusoid Estimation (original) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 439.99 Hz, epsilon = 10.22

--- Damped Sinusoid Estimation (noise) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 439.98 Hz, epsilon = 10.34

最後に、信号長 100 ms の結果を示す。減衰係数の推定誤差が小さくなったことが分かる。

--- Damped Sinusoid Estimation (original) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 439.99 Hz, epsilon = 9.97

--- Damped Sinusoid Estimation (noise) ---
True values:      f0 = 440.00 Hz, epsilon = 10.00
Estimated values: f0 = 440.01 Hz, epsilon = 10.08

おまけ

減衰振動SCDEの導出にあたり、 $2 \epsilon h_{\omega} / g_{\omega}$ が $\omega_{0}^{2}$ に比べて無視できるとしていた。この項が具体的にどれくらいのオーダーだったのかを確認しておく必要がある。

関数 solve_complex_a_sq_wls において、以下のコードを使って出力させてみた。

print(f"{np.max(np.abs(hwk / gwk)):.3f}")

結果は

信号長が 10 ms のとき 1680.162
信号長が 40 ms のとき 464.477
信号長が 100 ms のとき 111.616

であった。今回の実験では減衰係数 $\epsilon$ の値は 10 であるから、 $2 \epsilon h_{\omega} / g_{\omega}$ は

信号長が 10 ms のとき 33603.24
信号長が 40 ms のとき 9289.54
信号長が 100 ms のとき 2232.32

である。 $\omega_{0}$ はほぼ $2\pi * 440 \approx 2764$ であるので、 $\omega_{0}^{2} \approx (2764)^{2} \approx 7,640,000$ である。ここから、 $\omega_{0}^{2}$ の何%に相当するかを計算すると

信号長が 10 ms のとき 33603.24 / 7640000 = 約 0.44 %
信号長が 40 ms のとき 9289.54 / 7640000 = 約 0.12 %
信号長が 100 ms のとき 2232.32 / 7640000 = 約 0.03 %

となる。信号長に比例してオーダーは小さくなる。少なくとも、1 %を大きく下回り、0.1 % に迫るオーダーを確保できれば、十分な推定精度が得られるようである。

考察

40 ms 以上の信号長において、周波数と減衰係数を同時にかつ高精度で推定できた。ノイズを加えた後も推定精度が劣化しなかったことから、高い頑健性を示している。

信号長を 40 ms から100 ms へと長くすると、推定誤差はさらに減少し、真の値にほぼ一致するレベルであった。信号長の増加はつまり利用できる情報量の増加であり、これを推定精度の向上に直接結びつけることができたといえる。

しかし、信号長が 10 ms という短い信号長では、減衰係数の推定精度が、ノイズによって大きく低下する現象が観測された。これは本手法の限界を示唆している。減衰係数には、信号波形のエンベロープの情報が含まれるが、10 ms という短時間では、エンベロープの変化自体が微小である。そのため、ノイズの影響を受けやすくなったと考えられる。その一方で、周波数の推定は比較的安定していた。これは短い信号の中にも、推定に十分な周期構造の情報が含まれていると解釈できる。

おわりに

本記事では、単一の減衰正弦波に含まれる周波数と減衰係数という2つのパラメータを、SCDEに基づいて推定するアルゴリズム（減衰振動SCDE）を考察した。

Pythonによる簡単な検証実験により、ある程度の信号長（本実験では 40 ms 以上）が確保できれば、ノイズの多い環境でも、周波数と減衰を同時に推定できることが分かった。信号長が極端に短い場合（10 ms）には、減衰係数の推定が、ノイズに対して敏感になることも分かった。

今後は、減衰振動SCDEを複数の減衰正弦波へと拡張することを予定している。

追記：複数の減衰正弦波のパラメタ推定に応用した記事は以下。

tam5917.hatenablog.com

参考文献

[1] S. Ando and T. Nara, "An Exact Direct Method of Sinusoidal Parameter Estimation Derived From Finite Fourier Integral of Differential Equation," in IEEE Transactions on Signal Processing, vol. 57, no. 9, pp. 3317-3329, Sept. 2009.