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

はじめに

前回までの記事で、正弦波制約微分方程式（SCDE）に基づくパラメータ推定法について、その拡張を含めて紹介してきた。これまでのモデルは単一の音源を前提としていたが、多くの実用的な場面では、複数の音源が混在した信号を扱う必要がある。ピアノの和音や機械の動作音などが典型例である。

本記事では、複数の減衰正弦波が混合した信号から各成分のパラメタを推定する問題に取り組む。信号に含まれる成分のうち、最もエネルギーの大きいものから順番にパラメータを推定し、その成分を信号から差し引く処理を反復するのである。本記事ではこれを逐次除去法と呼ぶことにする。推定性能を検証するため、いくつかのシミュレーション実験を行ったので、結果を報告する。

はじめに
逐次除去法
実装
実験
- 実験1
- 実験2
- 実験3
おわりに

逐次除去法

複数の減衰正弦波が混合した信号から、全てのパラメタを一度に推定しようとするのは、誰もがまず最初に考えることである。しかし、多変数かつ複雑な非線形最適化問題を正確に解く必要があり、安定した解を得るのは難しい（ことが容易に予想される）。

そこで、この問題をより扱いやすい形に分解することを考える。具体的には以下の処理の流れを踏む。

1．入力された混合信号の中から、最もエネルギーの大きい単一の減衰正弦波成分を推定する。推定には、前回の記事で実装した、単一成分のパラメタ推定器（関数）を利用する。この関数は信号の周波数 $f_o$ と減衰係数 $\epsilon$ を出力する。

2．1.で得られた $f_o$ と $\epsilon$ を既知として、その成分の振幅 $A$ と位相 $\phi$ を、最小二乗法を用いて計算する。

3．1. および 2. で得られた全てのパラメタ $(f_o, \epsilon, A, \phi)$ を用いて、推定した単一成分の波形を再合成する。そして、その波形を元の混合信号から引き去り、「残差信号」を生成する。

4．この残差信号を新たな入力として、1. に戻る。これらのステップを、あらかじめ指定された音源の数 $N$ に達するまで繰り返す。

この手法は、初期のステップでの推定誤差が、後続のステップで用いる残差信号に影響を及ぼしうる。そのため、各ステップでの単一成分の推定精度の高さが求められる。

実装

以下に置いた。Enjoy!

scde_demo_wls_multi_dumped.py · GitHub

先述した逐次除去法に基づくアルゴリズムをPythonで実装した。ポイントとなる関数たちを挙げておく。

estimate_damped_sinusoid：単一の減衰正弦波の周波数 $f_o$ と減衰係数 $\epsilon$ を出力する。内部的には観測量の計算、周波数ビンの選択、WLSの解計算などが行われる。
estimate_amp_phase_damped： $f_o$ と $\epsilon$ が既知の状態で、その成分の振幅 $A$ と位相 $\phi$ を推定する。内部的には減衰項 $\exp(-\epsilon t)$ を考慮した線形の最小二乗問題を解いている。
synthesize_damped_sinusoid：推定された全パラメタ $(f_o, \epsilon, A, \phi)$ から、除去対象となる単一成分の波形を再合成する。
estimate_multiple_damped_sinusoids：上記関数を順次呼び出して逐次除去法のループ処理を実行する。

単一の減衰正弦波のパラメタを格納するためのNamedTupleである DampedSinusoidParameter を用意した。

class DampedSinusoidParameter(NamedTuple):
    """Represents the parameters of a single damped sinusoid."""

    fo: float  # fundamental frequency
    epsilon: float  # dumping coefficient
    amp: float  # amplitude
    phase: float  # initial phase

実験

逐次除去法に基づく複数減衰正弦波のパラメタ推定法の性能を検証するため、いくつかの条件下でシミュレーション実験を行った。

パラメタが既知の複数の減衰正弦波を合成した信号を生成し、それに対して実装したアルゴリズムを適用した。サンプリング周波数は 44100 Hz で固定する。信号長は各実験で別途指定する。そして推定された各パラメタと、設定した真の値との誤差を比較することで、アルゴリズムの精度を評価する。なお生成された信号には0.01でスケーリングしたホワイトノイズを加えている。

実験1

3つの音源が含まれる条件下で、アルゴリズムが安定して動作し、各成分を正しく分離・推定できるかを確認する。

最初の音源の設定では、周波数を互いに十分に離した（220.0 Hz, 415.3 Hz, 622.6 Hz）。415 Hzは220 Hzの倍音ではなく、622 Hzも他の倍音ではない。それぞれ異なる減衰率を持ち（5.0, 12.0, 8.0）、極端な強弱の差もない（1.0, 0.8, 0.9）。信号長は 0.1 s である。

DampedSinusoidParameter(fo=220.0, epsilon=5.0, amp=1.0, phase=np.pi / 4),
DampedSinusoidParameter(fo=415.3, epsilon=12.0, amp=0.8, phase=-np.pi / 3),
DampedSinusoidParameter(fo=622.6, epsilon=8.0, amp=0.9, phase=np.pi / 2)

推定結果を以下に示す。期待通り、高精度に推定できたことが確認できる。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	220.00	219.99	5.00	5.00	1.00	1.01
音源 2	415.30	415.30	12.00	12.06	0.80	0.80
音源 3	622.60	622.60	8.00	7.99	0.90	0.89

続いて、周波数がほぼ倍音関係にある状況を設定する。周波数は 220 Hz、そのほぼ2倍音である 441 Hz、そして660 Hz である。高い音ほど速く減衰するようにした（5.0, 8.0, 12.0）。振幅は倍音の方が少し小さくなるようにした（1.0, 0.7, 0.5）。

DampedSinusoidParameter(fo=220.0, epsilon=5.0, amp=1.0, phase=np.pi / 4),
DampedSinusoidParameter(fo=441.0, epsilon=8.0, amp=0.7, phase=-np.pi / 3),
DampedSinusoidParameter(fo=660.0, epsilon=12.0, amp=0.5, phase=np.pi / 2)

推定結果を以下に示す。この状況でも、うまく推定できたことが確認できる。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	220.00	220.00	5.00	5.00	1.00	1.00
音源 2	441.00	440.99	8.00	8.02	0.70	0.70
音源 3	660.00	659.99	12.00	12.05	0.50	0.50

実験2

実験1と比較して、より困難な条件下での性能を評価したい。そこで、(1) 成分間に極端な強弱差がある場合、(2) 信号長が短い場合でアルゴリズムの挙動を検証する。強度差をつけたのは、いわゆるマスキングを想定している。

(1) のマスキング条件では、3音源のうち1つの振幅を、他の音源よりも大幅に大きく設定した信号を生成した（400 Hz の音源が振幅 1.0, 他2つの音源は振幅 0.1 と 0.2）。信号長は 0.1 s である。

DampedSinusoidParameter(fo=250.0, epsilon=10.0, amp=0.1, phase=np.pi / 4),
DampedSinusoidParameter(fo=400.0, epsilon=5.0, amp=1.0, phase=-np.pi / 3),
DampedSinusoidParameter(fo=600.0, epsilon=8.0, amp=0.2, phase=np.pi / 2)

(2) の条件では、4音源が含まれる信号に対し、信号長を 0.1 s と 0.04 s の2つの条件で推定した。

DampedSinusoidParameter(fo=261.63, epsilon=4.0, amp=1.0, phase=np.pi / 4),
DampedSinusoidParameter(fo=329.63, epsilon=6.0, amp=0.7, phase=-np.pi / 3),
DampedSinusoidParameter(fo=392.00, epsilon=8.0, amp=0.6, phase=np.pi / 2),
DampedSinusoidParameter(fo=493.88, epsilon=10.0, amp=0.5, phase=np.pi / 6)

(1) の実験結果を以下に示す。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	250.00	249.94	10.00	10.08	0.10	0.10
音源 2	400.00	400.00	5.00	4.98	1.00	1.00
音源 3	600.00	599.99	8.00	7.79	0.20	0.20

信号長が 0.1s の場合、振幅に 10 倍の差がある $N = 3$ の信号に対しても、まず最も強い成分を正確に除去し、その下に隠れていた弱い成分を高い精度で推定することができている。

(2) の実験結果を以下に示す。まずは信号長 0.04 s である。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	261.63	276.70	4.00	18.03	1.00	0.65
音源 2	329.63	284.38	6.00	7.04	0.70	0.43
音源 3	392.00	284.29	8.00	21.25	0.60	0.11
音源 4	493.88	284.13	10.00	20.79	0.50	0.01

続いて信号長 0.1 s の結果を示す。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	261.63	261.63	4.00	3.93	1.00	1.02
音源 2	329.63	329.63	6.00	6.05	0.70	0.67
音源 3	392.00	391.99	8.00	7.93	0.60	0.61
音源 4	493.88	493.89	10.00	9.93	0.50	0.50

$N = 4$ の場合、信号長 0.1 s では推定に成功したが、信号長を 0.04 s まで短くすると、特に後から推定される成分の誤差が増大し、推定が破綻する結果となった。

各ステップでの推定に利用できる信号長が不足することによって、初期のステップでの推定精度が低下し、その除去誤差が残差信号にも残るため、後続ステップの推定精度も連鎖的に低下したと考えられる。

実験3

最後に、複数の周波数が互いに非常に近接している場合の性能評価を行う。440 Hzを中心に、15 Hz から20 Hz の間隔で周波数を密集させた（420 Hz, 440 Hz, 465 Hz, 480 Hz）。減衰係数はバラバラであり（10.0, 5.0, 12.0, 8.0）、振幅はほぼ同じレベルである（1.0, 0.9, 1.0, 0.8）。信号長は 0.1 s （短時間）と 1.0 s （長時間）の2つの条件とした。

DampedSinusoidParameter(fo=420.0, epsilon=10.0, amp=1.0, phase=np.pi / 4),
DampedSinusoidParameter(fo=440.0, epsilon=5.0, amp=0.9, phase=-np.pi / 3),
DampedSinusoidParameter(fo=465.0, epsilon=12.0, amp=1.0, phase=np.pi / 2),
DampedSinusoidParameter(fo=480.0, epsilon=8.0, amp=0.8, phase=np.pi / 6)

(3) の実験結果を示す。まずは信号長 0.1 s である。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	420.00	423.60	10.00	-27.56	1.00	0.05
音源 2	440.00	436.02	5.00	-45.99	0.90	0.01
音源 3	465.00	442.09	12.00	-28.29	1.00	0.06
音源 4	480.00	482.49	8.00	-13.05	0.80	0.21

続いて信号長 1.0 s である。

音源	真のFo (Hz)	推定Fo (Hz)	真の減衰係数	推定減衰係数	真の振幅	推定振幅
音源 1	420.00	419.99	10.00	10.03	1.00	0.95
音源 2	440.00	440.00	5.00	5.01	0.90	0.91
音源 3	465.00	464.98	12.00	11.99	1.00	0.95
音源 4	480.00	480.00	8.00	7.97	0.80	0.97

周波数が互いに近接している場合、信号長が 0.1 s であってもアルゴリズムは4つの近接した成分を分離できず、減衰係数の推定が破綻した。特に減衰係数が負の値を取っており、物理的にはありえない値となっている。信号長が 1.0 s まで確保されると、良好な推定結果が得られた。

手法の周波数分解能は、根本的にFFTの分析窓の長さ（信号長）に依存する。FFTが分離できないほど近接した周波数は、この手法でも分離することができない。裏を返せば、十分な信号長さえ与えられれば、高い周波数分解能を達成できることが示されたわけである。

おわりに

本記事では、SCDEの枠組みを拡張し、複数の減衰正弦波が混合した信号から、各成分のパラメタを推定するアルゴリズムを考察・実装した。

本記事で採用した逐次除去法では、信号中から最もエネルギーの大きい成分を一つずつ推定し、それを信号から差し引いていく処理を反復する。各ステップにおける単一成分のパラメタ推定には、前回記事で実装した推定器を利用した。

シミュレーション実験を通じて、実装したアルゴリズムの性能を検証した。その結果、

音源の数が3や4といった、より複雑な状況でも安定して動作した。
周波数が倍音関係にある場合や、成分間に大きな振幅差（マスキング）がある場合でも、各成分を良好に分離・推定できた。

一方で、分析に利用できる信号長に性能が強く依存することも明らかになった。信号長が十分に確保されていれば、非常に近接した周波数成分であっても高い精度で分離可能であったが、信号長が短い場合には、誤差の伝播が原因で推定が不安定になるケースも見られた。

今回実装した逐次除去法に基づく減衰振動SCDEは、複雑な音響信号から、その構成要素である減衰正弦波のパラメタを高精度に抽出できるするための有用なツールとなりうる。ただし性能限界を理解したうえで適切な信号長のもとで使用するならば、であるが。

今後は実際の楽器音の分析に応用したり、また逐次除去法以外の推定アルゴリズムを考えていきたい。