正弦波制約微分方程式（SCDE）に基づくピッチ推定法における重み付き最小二乗法（WLS）の導入

はじめに

過去に正弦波制約微分方程式（SCDE）に基づくピッチ推定法の再現実装に関する記事を書いた。

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

本記事はSCDE法を少しだけ拡張し、ピッチ推定精度の向上が期待できる方法を考察する。

はじめに
SCDEの復習 [1], [2]
- 境界値項をゼロにするための窓関数の導入
- 最小二乗法の確認
重み付き最小二乗法（Weighted Least Squares; WLS）への拡張
複数の独立した正弦波を推定するSCDEへの組み込み
- 準備と復習
- 重み付き最小二乗法の導入
実装
実験1
実験2
- 実験結果
- 考察
実験3
- 実験結果
- 考察
実験4
- 実験結果
- 考察
おわりに
参考文献

SCDEの復習 [1], [2]

単一の正弦波 $x(t)$ が満たす、2階微分方程式を考えることから始める。

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

ここで $a$ は求めたい未知の角周波数である。この式は全ての時刻 $t$ で成り立っている。しかし計算機は無限の時刻を扱えない。そこで「全ての時刻で成り立つ」という条件を、積分を使って、一つの式に集約する。

いま観測区間を $[-T/2, T/2]$ とすれば、2階微分方程式は下記の積分表現と等価である。

$\begin{align*} \int_{-T/2}^{T/2} \left[\frac{d^{2}}{dt^{2}} x(t) + a^{2} x(t) \right] w(t) \: dt = 0 \end{align*}$

ここで $w(t)$ は荷重を表す任意の連続関数である。フーリエ変換を見越して、 $w(t) = \exp(-j\omega t)$ としておく。 $j$ は虚数単位である。また $\omega$ は積分評価用の角周波数である。これは 自由に設定できる パラメータであり、求める対象ではない（既知の値）。上記の積分表現は次式となる。

$\begin{align*} \int_{-T/2}^{T/2} \left[\frac{d^{2}}{dt^{2}} x(t) \right] \exp(-j\omega t) \: dt + a^{2} \int_{-T/2}^{T/2} x(t) \exp(-j\omega t) \: dt = 0 \end{align*}$

左辺第2項には $x(t)$ のフーリエ変換が現れる。これを $g_{\omega}$ と置こう：

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

左辺第1項は部分積分を2回適用することで

$\begin{align*} \int_{-T/2}^{T/2} \left[\frac{d^{2}}{dt^{2}} x(t) \right] \exp(-j\omega t) \: dt = -\omega^{2} g_{\omega} + \delta \omega \end{align*}$

を得る。ここで $\delta \omega$ は積分区間の両端 $\pm T/2$ における $x(t)$ や導関数 $x^{\prime}(t)$ で書ける項を表す（境界値）。

積分表現は最終的に、 $a^{2}$ を未知数とする以下の代数方程式となる。

$\begin{align*} (-\omega^{2} g_{\omega} + \delta \omega) + a^{2} g_{\omega} = 0 \end{align*}$

境界値項をゼロにするための窓関数の導入

ちなみに、荷重関数 $w(t)$ を、特殊な窓関数 $p(t)$ を乗じた $w(t) = p(t)\exp(-j\omega t)$ へと変更するならば、境界値項をゼロにできる [1, 2]。この窓関数は、区間の両端 $\pm T/2$ において以下の条件を満たす原点対称の関数である。

関数の値がゼロ　 $p(T/2) = p(-T/2) = 0$
一階微分の値がゼロ　 $p^{\prime}(T/2) = p^{\prime}(-T/2) = 0$

Hann窓は両端で値がゼロになるが、一階微分の値はゼロにならない。文献 [1, 2] には、2階微分もゼロとなる窓関数が紹介されている。窓関数を導入した積分表現は次式となる。

$\begin{align*} \int_{-T/2}^{T/2} \left[\frac{d^{2}}{dt^{2}} x(t) \right] p(t) \exp(-j\omega t) \: dt + a^{2} \int_{-T/2}^{T/2} x(t) p(t) \exp(-j\omega t) \: dt = 0 \end{align*}$

先と同様に、左辺第1項の計算に部分積分を2回適用する。最終的に得られる代数方程式は

$\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*} a^{2} = \omega^{2} - \frac{k_{\omega} - 2 j\omega h_{\omega}}{g_{\omega}} \end{align*}$

最小二乗法の確認

SCDEで連立方程式を解く状況を考える。 $K$ 個の異なる角周波数 $\omega_k$ $(k=1, \ldots, K)$ を選び、それぞれについて微分方程式の荷重積分を計算すると、 $K$ 本の方程式が得られる。積分表現から得られる代数方程式に誤差項 $\epsilon_{k}$ を導入すると、

$\begin{align*} g_{\omega_{k}} \omega_{k}^{2} - a^{2} g_{\omega_{k}} + \delta \omega_{k} = \epsilon_{k} \end{align*}$

となる。ここで、適当な窓関数を導入して $\delta \omega_{k}$ が消去された（あるいは実質的に無視できる）と仮定する。未知数 $a^{2}$ に関する連立方程式は次式で与えられる。

$\begin{align*} \begin{bmatrix} g_{\omega_{1}} \omega_{1}^{2} \\ g_{\omega_{2}} \omega_{2}^{2} \\ \vdots \\ g_{\omega_{K}} \omega_{K}^{2} \end{bmatrix} - a^{2} \begin{bmatrix} g_{\omega_{1}} \\ g_{\omega_{2}} \\ \vdots \\ g_{\omega_{K}} \end{bmatrix} &= \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{K} \end{bmatrix} \end{align*}$

ベクトル表記を導入すれば

$\begin{align*} \mathbf{b}- a^{2}\mathbf{A} &= \boldsymbol{\epsilon} \end{align*}$

となる。ただし $\mathbf{A} = [g_{\omega_{1}}, g_{\omega_{2}} , \ldots, g_{\omega_{K}} ]^{\top}$ 、 $\mathbf{b} = [ g_{\omega_{1}} \omega_{1}^{2}, g_{\omega_{2}} \omega_{2}^{2}, \ldots, g_{\omega_{K}} \omega_{K}^{2}]^{\top}$ 、 $\boldsymbol{\epsilon} = [\epsilon_{1}, \epsilon_{2}, \ldots, \epsilon_{K} ]^{\top}$ と置いた。

残差ベクトル $\boldsymbol{\epsilon}$ のL2ノルム、あるいはそれと等価な残差の二乗和 $E = \sum_{k=1}^{K} |\epsilon_{k}|^{2}$ を最小にするのが、通常の最小二乗法である。その解は複素形の正規方程式 $\mathbf{A}^{*} \mathbf{A} x = \mathbf{A}^{*} \mathbf{b}$ を解くことで得られる（ $\mathbf{A}^{*}$ は随伴行列）。

$\begin{align*} a^{2} = \frac{\sum_{k=1}^{K} |g_{\omega_k}|^{2} \omega_{k}^{2} } {\sum_{k=1}^{K} |g_{\omega_k}|^{2}} \end{align*}$

上式はそれぞれの周波数ビンにおけるパワースペクトルで重みづけた形で角周波数を推定していると解釈できる。

重み付き最小二乗法（Weighted Least Squares; WLS）への拡張

先の最小二乗法の定式化においては、各方程式の誤差 $\epsilon_{k}$ は等分散、つまり同じ重要度を持つと仮定していた。しかしながら、信号の真の角周波数 $a$ から遠い $\omega_{k}$ で計算した式の信頼度は低いはずである（信念）。安藤氏の論文にも示唆されているように、スペクトルのピークに近い、SN比の高い情報を持つ $\omega_{k}$ （周波数ビン）を重視することで、推定の頑健性を高められる。そこで、各方程式に、その信頼度を表す重み $w_k$ を導入する。

つまり最小化の目的関数が

$\begin{align*} E ( \epsilon_{1}, \epsilon_{2}, \ldots, \epsilon_{K}; w_{1}, w_{2}, \ldots, w_{K}) = \sum_{k=1}^{K} w_k |\epsilon_{k}|^{2} \end{align*}$

へと変わる。ここで重みを対角成分に持つ重み行列 $\mathbf{W}$ を定義する。

$\begin{align*} \mathbf{W} = \mathrm{diag}(w_1, w_2, \ldots, w_K) \end{align*}$

すると上記の二乗和は $E = \boldsymbol{\epsilon}^{*}\mathbf{W} \boldsymbol{\epsilon} = (\mathbf{A} x - \mathbf{b})^{*} \mathbf{W} (\mathbf{A} x - \mathbf{b})$ と書くことができる。ただし $x = a^{2}$ である。これを最小化する解 $x$ は、重み付き正規方程式 $\mathbf{A}^{*} \mathbf{W} \mathbf{A} x = \mathbf{A}^{*} \mathbf{W} \mathbf{b}$ を解くことで得られる。その解は次式で与えられる。

$\begin{align*} a^{2} = \frac{\sum_{k=1}^{K} w_k |g_{\omega_k}|^{2} \omega_{k}^{2} } {\sum_{k=1}^{K} w_k |g_{\omega_k}|^{2}} \end{align*}$

重みの設定にはいくつかの選択肢が考えられる。一番シンプルなものは、振幅スペクトル（パワースペクトル）に基づく重み付けであろう。つまり $w_k = |g_{\omega_k}|$ や $w_k = |g_{\omega_k}|^{2}$ とするのである。一般には、べき乗つき振幅スペクトル $w_k = |g_{\omega_k}|^{p}$ $(p \geq 1)$ でも良い。その他にも設定は考えられるのだが、それはまたの機会としよう（出し惜しみ）。もし $w_k = |g_{\omega_k}|^{2}$ と設定すれば、実質的にパワーの2乗 $|g_{\omega_k}|^{4}$ で重み付けることになる。

ちなみに安藤氏の論文では $k$ F ( $k$ weight frequency) 法とも呼ばれる [1] 。

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

山田氏によって提案された、複数の独立した正弦波を推定するSCDE [3] に対して、重み付き最小二乗法を組み込んでみることを考える。

これには動機がある。山田氏のSCDEでは、時間領域の高階微分を数値的に計算する。そのため、信号に含まれる高周波ノイズを大きく増幅させ、計算を不安定にするおそれがあった。この問題を回避するため、既に上で説明した積分表現およびフーリエ変換を経由し、SCDEで必要となる様々な微分の積和を、フーリエスペクトルと角周波数のべき乗の積和で計算する方針を取った。

準備と復習

まず、信号は $N$ 個の正弦波からなると仮定する。

$\begin{align*} x(t) &= A_1 \cos (a_1 t + \phi_1) + A_2 \cos (a_2 t + \phi_2) + \cdots + A_N \cos (a_N t + \phi_N) \\ &= \sum_{i = 1}^{N} A_i \cos(a_i t + \phi_i) \end{align*}$

未知となるのは $N$ 個の角周波数 $a_1, a_2, \ldots, a_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*}$

$\alpha_n$ は $a_{n}^{2}$ の基本対称式である。この方程式を、 $K$ 個の角周波数で荷重積分することで、 $K$ 本の連立方程式が得られる。

$\begin{align*} (-\omega_{k}^{2})^{N} g_{\omega_k} + \alpha_{1} (-\omega_{k}^{2})^{N-1} g_{\omega_k} + \cdots + \alpha_{N} g_{\omega_{k}^{2}} = \epsilon_k \end{align*}$

これも未知数 $\mathbf{x} = [\alpha_{1}, \ldots, \alpha_{N}]^{\top}$ に関する線形方程式 $\mathbf{A} \mathbf{x} - \mathbf{b} = \boldsymbol{\epsilon}$ の形で書き直すことができる。ここで $\mathbf{A} = (A_{k, n})_{1 \leq k \leq K,\; 1 \leq n \leq N}$ の成分は

$\begin{align*} A_{k, n} = (-\omega_{k}^{2})^{(N-n)} g_{\omega_{k}} \end{align*}$

で与えられる。また $\mathbf{b} = (b_k)_{1 \leq k \leq K}$ の成分は

$\begin{align*} b_{k} = -(-\omega_{k}^{2})^{N} g_{\omega_{k}} \end{align*}$

で与えられる。

未知数が $N$ 個なので、方程式は少なくとも $N$ 本必要である。つまり $K\geq N$ である。通常のSCDEでは、 $K=N$ となるように周波数ビンを選び、方程式を解いているのであった。

重み付き最小二乗法の導入

$K > N$ となるよう、周波数ビンを多く選ぶ。つまり方程式の数が未知数より多い状況を考えて、重み付き最小二乗解を求める。解くべき重み付き正規方程式は、

$\begin{align*} \mathbf{A}^{*} \mathbf{W} \mathbf{A} x = \mathbf{A}^{*} \mathbf{W} \mathbf{b} \end{align*}$

である。最適な係数 $\alpha_1, \ldots, \alpha_N$ が求まれば、後の流れは全く同じである。すなわち、以下の $N$ 次代数方程式の根 $\Omega_n$ を求め、 $a_n = \sqrt{\Omega_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*}$

以降、この手法を WLS-SCDE と呼ぶことにする。

実装

以下に置いた。Enjoy!

scde_demo_wls.py · GitHub

実装のポイントをまとめておく。

1．スペクトル計算 (calculate_spectrum)

この関数は、入力された時間領域の信号フレームを受け取り、後続の処理で必要となる複素スペクトルと周波数軸を計算する。フレーム境界の不連続性から生じるスペクトル漏れを抑制するため、窓の両端が滑らかにゼロに収束するブラックマンハリス窓を適用した。安藤氏の論文で示唆されている「境界項の効果」を抑制する働きも期待できる。

2．周波数ビン選択 (select_bins_robust)

計算されたスペクトルの中から、 $N$ 個の未知の周波数を解くために、推定に有利な情報を持つ $K$ 個の周波数ビンを選択する。単純にパワースペクトルピークの上位 $K$ 個を選ぶだけでは、近接した周波数が単一のピークにくっついてしまうケースに対応できない。そこで、scipy.signal.find_peaks を用いて大まかなピーク位置を特定し、各ピークの周辺領域からバランス良く周波数ビンを収集するようにした。

安藤氏の論文 [2] では、スペクトルピークの周辺から、連続した周波数ビンを選ぶのが良い戦略である旨が述べられている。しかし、具体的にどのようにビンを選択するかについては述べられているわけではない。実は具体的な周波数ビン選択については、いくつかアルゴリズムを試行錯誤しており、今回はその中で一番結果が安定したものを実装として紹介している。後日、アルゴリズム間の性能を比較する記事を書く予定。

実験1

WLS-SCDEが、山田氏の時間領域SCDEと比較して、どの程度の性能向上を達成できるかを検証する。

周波数 440 Hz, 460 Hz, 480 Hzという、周波数が近接する3つの正弦波の和で書ける信号の周波数を推定する実験を行った。標本化周波数は 44100 kHz である。信号にホワイトノイズを加えたうえで推定させる実験では、ノイズレベルは 10 dB、ノイズの周波数帯域は3つの周波数を覆う 390 Hz〜580 Hz とした。

実験結果その1（信号長が 0.2 sec の場合）

信号長が0.2 secの場合の結果を示す。まず、ノイズが一切観測されない場合について、時間領域SCDEと性能差はあまり無いように見える。

--- Test without noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 439.928  459.918  479.906 ] Hz
Estimation errors: [ 0.072  0.082  0.094 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.052  460.274  479.831 ] Hz
Estimation errors: [ -0.052  -0.274  0.169 ] Hz

続いて、信号に帯域制限つきホワイトノイズを加えた場合の結果を示す。時間領域SCDEは推定精度を落としている。

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 442.329  474.256  515.471 ] Hz
Estimation errors: [ -2.329  -14.256  -35.471 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.050  460.270  479.832 ] Hz
Estimation errors: [ -0.050  -0.270  0.168 ] Hz

実験結果その2（信号長が 1.0 sec の場合）

信号長が1.0 secの場合の結果を示す。まずノイズがない場合において、WLS-SCDEの推定精度の向上が確認できた。

--- Test without noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 439.928  459.918  479.906 ] Hz
Estimation errors: [ 0.072  0.082  0.094 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.003  460.094  479.959 ] Hz
Estimation errors: [ -0.003  -0.094  0.041 ] Hz

さらにノイズを加えた場合についても、WLS-SCDEは推定精度をキープしている一方で、時間領域SCDEは推定精度が低下した。

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 442.271  473.846  508.942 ] Hz
Estimation errors: [ -2.271  -13.846  -28.942 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.003  460.098  479.962 ] Hz
Estimation errors: [ -0.003  -0.098  0.038 ] Hz

実験結果その3（信号長が 40 msec の場合）

信号長が 40 msec の場合の結果を示す。ノイズの有無に関わらず、WLS-SCDEは推定精度が低下した。時間領域SCDEはノイズが加えられる場合に推定精度が低下した。

Sampling Frequency: 44100 Hz
Duration: 0.04 sec
--- Test without noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 439.928  459.918  479.906 ] Hz
Estimation errors: [ 0.072  0.082  0.094 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 427.653  465.748  492.976 ] Hz
Estimation errors: [ 12.347  -5.748  -12.976 ] Hz

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.829  476.498  519.118 ] Hz
Estimation errors: [ -0.829  -16.498  -39.118 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 427.668  464.290  492.962 ] Hz
Estimation errors: [ 12.332  -4.290  -12.962 ] Hz

考察

時間領域SCDEは、いずれの信号長のときもノイズに対する脆弱性を示した。WLS-SCDEは、信号長が 0.2 sec と 1.0 sec のとき、ノイズに対して高い頑健性を示した。これは周波数領域での計算とWLSによる重み付けが効果的に働いたことを示している。また信号長が 1.0 sec のような十分に長い信号長の下では、WLS-SCDEの推定誤差（平均 0.05 Hz）は、時間領域SCDEの推定誤差（平均 15 Hz）と比較して約 300 倍ほど良い結果となった。これはノイズの有無に関わらない。ただし信号長が 40 msec のとき、WLS-SCDEの推定精度は顕著に低下しており、信号長と推定精度にはトレードオフが存在することが分かった。

時間領域SCDEは、40 msec のような短い信号に対しても推定が可能である一方で、やはりノイズや数値誤差に敏感であると考えられる。1.0 sec のようなある程度の長さの信号において精度向上が見られなかったことから、数値微分の誤差蓄積が性能上限を決めている可能性はある。
WLS-SCDEは、FFTをベースとするため、信号長が長くなるほど周波数分解能が向上し、推定精度が向上する。ただ、40 msec での推定精度低下はむしろ推定対象となる正弦波の周波数が近接しているためである。詳しくは後述する実験2で明らかとなる。

実験2

WLS-SCDEのノイズ耐性を調べるために、さらに実験を行った。実験1と同様、3つの正弦波の和で書ける信号の周波数を推定する。ただしここでは、周波数は 100 Hz, 200 Hz, 300 Hz のように、実験1よりも周波数間の距離を離した。また信号長は時間領域SCDEに有利な 40 msec とした。標本化周波数は 44100 kHz である。信号にホワイトノイズを加えたうえで推定させる実験では、ノイズレベルは 10 dB、ノイズの周波数帯域は3つの周波数を覆う50 Hz〜500 Hz とした。

実験結果

まず、ノイズを加える前の推定では、WLS-SCDEが若干推定精度が低い結果となった。これは周波数分解能の点からやむを得ない。

--- Test without noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000 ] Hz
Estimated frequencies with noise: [ 99.999  199.993  299.977 ] Hz
Estimation errors: [ 0.001  0.007  0.023 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 100.000  200.000  300.000 ] Hz
Estimated frequencies with noise: [ 99.880  199.832  298.992 ] Hz
Estimation errors: [ 0.120  0.168  1.008 ] Hz

一方、ノイズを加えた後の推定では、時間領域SCDEは大きく推定精度を落としている。WLS-SCDEは依然として推定精度をキープしている。

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000 ] Hz
Estimated frequencies with noise: [ 121.236  282.338  362.077 ] Hz
Estimation errors: [ -21.236  -82.338  -62.077 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 100.000  200.000  300.000 ] Hz
Estimated frequencies with noise: [ 99.930  199.854  301.623 ] Hz
Estimation errors: [ 0.070  0.146  -1.623 ] Hz

考察

ノイズが加わった場合、時間領域SCDEは大きく推定精度を落としたが、まさにノイズに対する脆弱性であり、その原因は高階の数値微分にあると考えられる。今回の実験のように、周波数間の距離が十分に開いている場合（> 50 Hz）には、WLS-SCDEは 40 msec という短い信号長であっても強いノイズ耐性を示した。実験1のように各正弦波の周波数が近接している場合（< 50 Hz）には、推定精度が低下する。これは今回の実装で採用したスペクトルのピーク抽出のアルゴリズムが、近接した周波数の検出を苦手としている、と解釈できる。スペクトルピーク位置の推定アルゴリズムにはまだ検討の余地がある。

実験3

ここでは推定する正弦波の個数を増やして実験してみる。具体的には5つの正弦波の周波数を推定する。周波数はそれぞれ 100 Hz, 200 Hz, 300 Hz, 400 Hz, 500 Hz である。信号長は 40 msec と 1.0 sec である。帯域制限ノイズの周波数範囲は 50 Hz 〜 550 Hz である。ノイズレベルはこれまでと同様に 10 dB とした。

実験結果

信号長が 40 msec の実験結果を以下に示す。時間領域SCDEは、数値微分の計算誤差の影響を大きく受け、推定精度が大幅に低下した。WLS-SCDEはそれほど精度は低下していない。

Sampling Frequency: 44100 Hz
Duration: 0.04 sec
--- Test without noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 117.954  264.010  393.173  499.638  2218.970 ] Hz
Estimation errors: [ -17.954  -64.010  -93.173  -99.638  -1718.970 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 104.301  215.761  312.210  403.963  500.891 ] Hz
Estimation errors: [ -4.301  -15.761  -12.210  -3.963  -0.891 ] Hz

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 121.114  295.739  405.727  502.748  3732.291 ] Hz
Estimation errors: [ -21.114  -95.739  -105.727  -102.748  -3232.291 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 105.157  216.438  315.815  401.758  500.894 ] Hz
Estimation errors: [ -5.157  -16.438  -15.815  -1.758  -0.894 ] Hz

続いて信号長が 1.0 secのときの結果を以下に示す。時間領域SCDEでは 12438.928 Hz や 4648.440 Hzという推定結果が現れており、推定が破綻している。WLS-SCDEは十分な信号長を確保したおかげで、高精度かつ安定な推定を実現していた。

Sampling Frequency: 44100 Hz
Duration: 1.0 sec
--- Test without noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 118.075  264.323  393.269  499.652  12438.928 ] Hz
Estimation errors: [ -18.075  -64.323  -93.269  -99.652  -11938.928 ] Hz

--- Test without noise (WLS) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 100.018  200.044  300.048  399.966  499.975 ] Hz
Estimation errors: [ -0.018  -0.044  -0.048  0.034  0.025 ] Hz

--- Test with band-limited noise (YAMADA) ---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 127.063  290.810  412.217  501.550  4648.440 ] Hz
Estimation errors: [ -27.063  -90.810  -112.217  -101.550  -4148.440 ] Hz

--- Test with band-limited noise (WLS)---
True frequencies: [ 100.000  200.000  300.000  400.000  500.000 ] Hz
Estimated frequencies with noise: [ 100.018  200.045  300.047  399.966  499.975 ] Hz
Estimation errors: [ -0.018  -0.045  -0.047  0.034  0.025 ] Hz

考察

$N$ 個の正弦波を推定するため、時間領域SCDEは $2N$ 階の微分方程式を解く必要がある。実験1と2では、 $N = 3$ だったが、これは 6 階の微分が必要であった。実験3では、 $N = 5$ であり、最高で 10 階の微分が必要である。今回の時間領域SCDEの実装では数値微分の計算に中心差分法を採用したが、微分計算のたびに誤差が蓄積し、指数関数的に増幅する。結果的に連立方程式の係数行列が不安定になり（いわゆるill-conditionedな行列）、周波数の推定精度が極端に悪くなったと考えられる。WLS-SCDEは時間領域で数値微分を計算せず、周波数領域での代数計算、すなわち周波数のべき乗計算と積和になっており、数値誤差は発生しない。

今回の実験を通して、時間領域SCDEは、中心差分法に基づく数値微分では $N$ が小さい（2か3）という限定された条件でしか機能しないことが分かった。山田氏の修士論文 [3] によれば、リチャードソン加速のような、より高精度な数値微分手法を導入することで、時間領域SCDEの推定精度が改善されることが報告されている。この点は、時間領域SCDEのさらなる改善の可能性を示唆している。しかし、そのような高精度化を行っても、高階微分に伴う根本的な誤差増幅の問題がどの程度解決されるかは、さらなる検証が必要である。 WLS-SCDEは $N$ が増えても推定は破綻せず、汎用性と頑健性を備えていることが示された。これは多数の正弦波が混在する、より現実的な信号の分析において、大きな利点を持つといえる。

実験4

実験1では、440 Hz 、 460 Hz 、480 Hz という正弦波の周波数に対して、ノイズの帯域が 390 Hz 〜 580 Hz という広帯域だった。ここでは山田氏の修士論文に合わせて、ノイズの帯域を 430 Hz 〜 490 Hz という狭帯域にして実験する。ノイズレベルは 10 dB である。信号長は 0.2 sec とした。

実験結果

比較のため、広帯域の結果も合わせて示す。時間領域SCDEは狭帯域ノイズという状況では優れた耐雑音性能を示している。

--- Test with band-limited noise (430.0 Hz - 490.0 Hz) (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 439.929  459.903  479.911 ] Hz
Estimation errors: [ 0.071  0.097  0.089 ] Hz

--- Test with band-limited noise (390.0 Hz - 580.0 Hz) (YAMADA) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 442.571  474.250  521.167 ] Hz
Estimation errors: [ -2.571  -14.250  -41.167 ] Hz

参考まで、WLS-SCDEの結果も示す。

--- Test with band-limited noise (430.0 Hz - 490.0 Hz) (WLS) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.001  460.093  479.959 ] Hz
Estimation errors: [ -0.001  -0.093  0.041 ] Hz

--- Test with band-limited noise (390.0 Hz - 580.0 Hz) (WLS) ---
True frequencies: [ 440.000  460.000  480.000 ] Hz
Estimated frequencies with noise: [ 440.054  460.280  479.832 ] Hz
Estimation errors: [ -0.054  -0.280  0.168 ] Hz

考察

狭帯域ノイズ下では時間領域SCDEは優れた耐雑音性能を示したが、この条件は「信号の周波数帯域が事前にかなり正確に分かっている」という、強い事前知識を仮定している。つまり、高性能なバンドパスフィルタを、前処理として適用できることを前提としている。広帯域ノイズ下では時間領域SCDEの性能は大きく劣化したが、WLS-SCDEは帯域が広がっても性能がそれほど変わらなかった。その理由は大きく2つの要素が考えられる。一つは、WLS-SCDEの実装で用いたブラックマンハリス窓である。この窓は非常に急峻なフィルタ特性を持っており、分析対象の周波数帯域の外側からの影響、いわゆる「スペクトル漏れ」を効果的に抑制したと考えられる。もう一つは今回の実装で採用したビン選択のアルゴリズムである。スペクトルのピーク周辺（高SN比の領域）にうまくフォーカスしたことで、580Hz付近のノイズを分析対象外として効果的に無視したと考えられる。

時間領域SCDEは、(1) $N$ （音源数）が小さく、 (2) ノイズがない、あるいは極めて狭い帯域に限定されている、という理想的な状況下で性能を発揮するアルゴリズムである。WLS-SCDEは、 $N$ が増えても破綻せず、ノイズの帯域が広がっても性能が大きく劣化しない。十分な信号長があれば、精度においても時間領域SCDEを超えており、より広い範囲の、より現実的な条件下で機能する手法といえる。

おわりに

本記事では正弦波制約微分方程式（SCDE）に基づくピッチ推定法に重み付き最小二乗法（WLS）を導入した。実験の結果、WLS-SCDEは近接した正弦波の抽出をそれほど得意としないものの、優れたノイズ耐性を示し、WLSの有効性は確認できた。

参考文献

[1] 安藤繁, 短時間正弦波パラメータ推定の厳密直接解法 : 雑音と干渉正弦波の統計的解析, 信学技報, vol. 107, no. 62, EA2007-11, pp. 25-30, 2007.

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

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