統計的信号処理を用いた合奏音からの楽器パート音の分離

全文

(1)Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. 統計的信号処理を用いた合奏音からの楽器パート音の分離 †. 深山幸穂田中大介. ††. †. 高本健生佐々木崇人. 1. はじめに本報では統計的最小二乗規範に基づき，合奏音のモノラル信号を指定された楽器パートごとの演奏音に分離する信号処理法を検討する．よく知られているように，MIDI 標準 [2] が規定するノート番号 i の音高の名目上の周波数は，国際標準ピッチ(1939)に従う. †. i − 69. fi = 440 × 2 12 [Hz] となり， " , f 69 , f 70 , f 71 " は大譜表における中央の C 音の直ぐ上の " , A, A#, B, " 音に対応している．このとき，楽器パートの音域は，例えば，標準的なピアノの全鍵盤を用いれば 21 ≤ i ≤ 108 となる．本研究では，バイオリン，トランペット等，第 p 番目の楽器パートの演奏音を分. 本報では統計的最小二乗規範に基づき，合奏音のモノラル信号を指定された楽器パートごとの演奏音に分離する信号処理法を検討する．前提として，対象となる楽器パートの音域内の音高の異なる各音について，自己相関行列と相互相関行列を既知とした．この信号処理法では合奏音を前述の各音に音量係数を乗じた和と仮定し，まず，入力した合奏音の自己相関行列から各音の音量係数を同定し，次に，合奏音と各音との相互相関行列を算出している．そして，これらの自己相関行列と相互相関行列から，楽器パートごとの音を推定するフィルタの係数を求める．バイオリン音，トランペット音，フルート音による実施例を付す．. 12 音平均律[4]により. 離するため，その楽器の音域に属するノート番号の個数だけ当該音を推定するフィルタを用意し，これらを通過した信号を加算する．上述したフィルタの設計にあたり，合奏の各楽器パートで用いるノート番号の異なる. Separation of each Instrumental part in music sounds with a statistical signal processing. 音を「要素となる音」と称することとし，これらにユニークな音番号. l = 1, 2," , L. を. 付して取扱う．すなわち，合奏音を要素となる音に音量係数を乗じた和としてモデル化する．このとき，要素となる音それぞれは互いに広義の定常な確率過程によるランダムな波形であり，統計量は既知と仮定している．さらに，音量係数は対応する音の波形に比して緩慢に時間変化し，短い時間では定数と考える．以下，次節より統計的信号処理の考え方[3]に基づいて具体的に議論する．. Yukio Fukayama†, Daisuke Tanaka†, Kensei Takamoto† and Takato Sasaki† An optimal filter that separates a tone of particular pitch highness played by particular instrument applying the statistical least square criterion is discussed. The filter is based on a priori autocorrelations of and cross correlations among tones. It firstly identifies current amplitude factors of tones; then, it evaluates cross correlation between current input signal and each tone. Finally it obtains the estimate of seeking tone from the input signal which is mixture of sounds played by instruments and corrupted by white noise. Some numerical examples with violin, trumpet and flute tones are also introduced.. 2. 入力信号のモデル化と音量係数の同定本節では，第. l = l (i , p ) 数. σl. p. 番目の楽器パートで演奏され，ノート番号. 番目の要素となる音. { Sl ( k )} ， (l = 1, 2," , L ). の求め方を検討する．ここに，時点. の下に，時刻. t = t0 + k +t. k. i. の音高である，第. について，その音量係. はサンプリング間隔. +t ，開始時刻 t0. における値を示す．. † 新居浜工業高等専門学校電子制御工学科 Electronic Control Department, Niihama National College of Technology. †† 新居浜工業高等専門学校専攻科 Advanced Study Courses, Niihama National College of Technology. 1. ⓒ 2011 Information Processing Society of Japan.

(2) Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. まず，入力信号. {Y ( k )}. は確率過程であり，次式にモデル化する:. ここに，音量係数ベクトル. Y ( k ) = ∑ σ l Sl ( k ) + W ( k ) ，. σ = (σ 1 " σ L ) ．次に，入力信号を 2M 近傍においてベクトルで表記する：  S1 ( k )    T Y ( k ) = (σ ⊗ I 2 M +1 )  #  + W ( k )  S (k )   L . (1). l =1. ここに，要素となる音と音量係数との積和を合奏音とし，白色雑音て入力信号となった．また，音量係数. σl. は時点. k. {W ( k )}. が混入し. と共に緩慢に変化し，後述する. 近傍では定数とみなせる.. このとき. { Sl ( k )}. と. {W ( k )}. は互いに広義の定常(w.s.s.)な確率過程である．すな. わち，統計量として，平均は零，自己相関と相互相関. (l , j = 1, 2," , L ). m の関数と仮定する： E{ Sl ( k )} = E{W ( k )} = 0 ，. は時点. k. に. = (σ 1 I 2 M +1. よらず，時点の差. rlj ( m), ( m ≤ 2 M ) ， E{ Sl ( k + m)S j ( k )} =   0, ( m ≥ 2M + 1)  w2 , ( m = 0) ， E{W ( k + m)W ( k )} =   0, ( m ≠ 0) E{ Sl ( k + m)W ( k )} = 0, (for all m) ，ここに，その外側において. k. の. 2M. { Sl ( k )}. (2). また，信号のエネルギ. (3). 番目の要素となる音のベクトル音ベクトル. e. ， (l. (7). (8). Y ( k ) ∈ R 2 M +1 ,. 第. l. = 1, 2," , L ) ，および，白色雑. を定義した：. (5). Sl ( k ) = ( Sl ( k + M ) " Sl ( k − M ) ). T. (9). ， T. W ( k ) = (W ( k + M ) " W ( k − M ) ). の自己相関や相互相関が無視できるように，時点. T. ，. (10). ．. (11). RYY ∈ R ( 2 M +1)×( 2 M +1) は {Y ( k )} から直接に評価でき, さらに，緩慢に時間変化する σ を用いて次式にモデル化できる： RYY = E{Y ( k )Y T ( k )}. このとき，自己相関行列. は (1) 式より次式に評価できる：.  E{ S1T ( k )S1 ( k )} " E{ S1T ( k )S L ( k )}    2 = σT  # % #  σ + E{(W ( k )) }  E{ S LT ( k )S1 ( k )} " E{ S LT ( k )S L ( k )}     r11 (0) " r1L (0)    =σ  # % #  σ + w2 ，  r ( 0) " r ( 0)  LL  L1 . Sl ( k ) ∈ R. Y ( k ) = (Y ( k + M ) " Y ( k − M ) ). 2. T. W (k ) ∈ R. 2 M +1. 2 M +1. (4). e = E{(Y ( k )) } 2.  S1 ( k )    " σ L I 2 M +1 )  #  + W ( k ) ，  S (k )   L . ここに，⊗ はクロネッカ積であり，次のとおり，信号ベクトル. 近傍を設定した． 2. を定義した：. T. L. 2M. σ ∈ RL.  R11 " R1L    = (σ ⊗ I 2 M +1 )  # % #  (σ ⊗ I 2 M +1 ) + w2 I 2 M +1 ， (12) R " R  LL   L1 ( 2 M +1)×( 2 M +1) ここに， Rlj ∈ R は，要素となる Sl ( k ) および S j ( k ) ， (l , j = 1, 2," , L ) について，自己相関行列 (l = j ) または相互相関行列 (l ≠ j ) であり，時点 k によらないため，予め評価しておける： T. (6). 2. ⓒ 2011 Information Processing Society of Japan.

(3) Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. bm (≥ 0) , ( m = 1, 2," , 2 M + 1). " rlj ( 2 M )   rlj (0)   # % # ． Rlj = E{ Sl ( k ) S Tj ( k )} =  (13)  r ( −2M ) " r (0)  lj  lj  このとき， RYY は正定対称行列であるから次のとおり規格化直交系をなす固有ベクトル. v m ∈ R 2 M +1. J (σ ) =. vmT. いま，固有ベクトル. および. F (σ ) = =. をそれぞれ左右から. (12) 式に乗じると次式を得る：. さらに，(6) 式を用いると，. w2. f (σ ) ∈ R L. (21) とヘッセ行列. F ( σ ) ∈ R L× L. は以下とな. (22). m. (Vm + VmT )[σσ T (Vm + VmT ) + (σ TVm σ − γ m ) I L ] ． σ (0). の下，次式の収束計算. (i = 1, 2,"). (23) により. (16). σ ( i ) = σ ( i −1) − F −1 (σ ( i −1) ) f (σ ( i −1) ) ， (i ) ここに， σ = σ は次式の収束判定条件 ε (i ) ( i −1) σ −σ <ε .. を消去できる：. をもって与える：. (25). 3. 楽器パート音分離のためのフィルタ (18). 本節では，入力信号ノート番号. (19). ここに，ロバスト性を得るため，次式のあてはめ誤差. {Y ( k )}. に含まれている，第. の音高である，第. まず，合奏音中の第. を最小二乗. l = l (i , p ). p. 番目の音. 番目の楽器パートで演奏され，. {U l ( k )} ， (l = 1, 2," , L ). {Uˆ l ( k )} を得るフィルタを導出する． l 番目の音は次式にモデル化できる：. U l ( k ) = σ l Sl ( k ) ， (20). J (σ ). i. について，その推定値. 解として同定できる：. σ = arg min{ J (σ )} ，. を決. (24). を定義した：. T T  r11 (0) " r1L (0)   v m R11v m " v m R1L vm      Vm =  # % # − # % # ， T T  r ( 0) " r ( 0)   v R v " v R v  LL m LL m   L1   m L1 m 2 γ m = e − λm ．従って，通常は L 2 M + 1 であるから，(17) 式を満たす σ (≥ 0 L ). σ. 定できる：. (17). γm. ∑ 2b. よって，適当な初期条件. γ m = σ TVm σ ，ここに， Vm および. (σ TVm σ − γ m )2 ．. ∂ 2 J (σ ) ∂f (σ ) = ∂σ ∂σ T ∂σ T. 2 M +1 m =1.  vmT R11vm " vmT R1L vm    2 % # λm = σ T  # σ + w ，  vmT RL1vm " vmT RLL vm   . λm に応じて与えればよい:. ∂J (σ ) 2 M +1 = ∑ 2bm (Vm + VmT )σ (σ TVm σ − γ m ) ， ∂σ m =1. f (σ ) =. (15). v m ，( m = 1, 2," , 2 M + 1). を固有値. る：. ここに，対応する固有値. RYY vm = λm vm ,. m. このとき，勾配ベクトル. (14). λm は以下のように置ける： (λ1 ≥ λ2 ≥ " ≥ λ2 M +1 ≥ 0) ．. ∑b m =1. を選ぶことができる：. vmT vn = δ mn ， ( m, n = 1, 2," , 2M + 1) ，. 2 M +1. ここに，音量係数. において重み係数. (26). σ l ，および，要素となる音 { Sl ( k )} ， (l = 1, 2," , L ). に係わる. 事項は前節と同様である．. 3. ⓒ 2011 Information Processing Society of Japan.

(4) Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. いま，推定値. {Uˆ l ( k )}. を得るため，インパルス応答. { hl ( m )} ，( − M ≤ m ≤ M ). なる音どうしの相互相関行列. Uˆ l ( k ) =. ∑. m =− M. hl ( m)Y ( k − m) = hlTY ( k ) ，. (27). ここに，インパルス応答ベクトルを定義した：. hl = ( hl ( − M ) hl ( − M + 1) " hl ( M ) ) このとき，統計的最小二乗規範により，第. l. T. ．. (33). 4. 統計量の評価. 番目の音の推定に最適なインパルス応. hl ， (l = 1, 2," , L ) は次式の評価関数を最小にする： C( hl ) = E{(U l ( k ) − Uˆ l ( k ))2 }. 前節までに述べた. = σ l2 rll (0) − 2σ l hlT rYl + hlT RYY hl. まず，第. −1 −1 −1 = ( hl − σ l RYY rYl )T RYY ( hl − σ l RYY rYl ) + σ l2 [rll (0) − rYTl RYY rYl ]. (30). 相互相関ベクトル. RYY. rYl ∈ R. は正定対称で逆行列を有する．また，入力信号のと次式の. 2 M +1. {Uˆ l ( k )}. の推定には， Rlj および. RYY. が必要で. l. 番目と第. j. 番目の要素となる音. { Sl ( k )}. と. { S j ( k )}. Rlj ，(l , j = 1, 2," , L ). は確率過程で. は次の構成である：. " rlj ( 2M )   rlj (0)   Rlj =  # % # ， (34)  r ( −2 M ) " r ( 0 )  lj  lj  これは時点 k に依存しないため，事前にそれぞれの音を単独に演奏して見本過程 { sl (κ )} と { s j (κ )} に用いれば，上式中の rlj ( m) ， ( m = 0, ±1," , ±2 M ) を N 個のデータの時間平均として以下のとおり評価できる： 1 N rlj ( m) = E{ Sl ( k + m)S j ( k )} = ∑ sl (κ + m)s j (κ ) ， (35) N κ =1 ここに，時点 κ は事前の演奏の開始時刻 τ 0 にて，時刻 τ = τ 0 + κ +t の値を示す．一方，確率過程 {Y ( k )} について，自己相関行列 RYY は時点 k により緩慢に変化し，見本過程として実測された入力信号 { y( k )} を成分とするベクトル y ( k ) ∈ R 2 M +1 の時間平均として評価できる．さらに，忘却係数 α , (0 < α < 1) の. 以上の展開において，前節に述べた互いに広義の定常な性質を適用した．従って，(29) 式の等号成立より最適なインパルス応答ベクトルが得られる：. ここに，自己相関行列. の同定や. あり，自己相関行列または相互相関行列 (29). −1 hl = σ l RYY rYl ，. σ. あり，本節では，二次の定常確率過程におけるエルゴード性を仮定して，これらの期待値を見本過程の時間平均から評価する．. = E{(σ l Sl ( k ) − hlTY ( k ))(σ l Sl ( k ) − hlTY ( k ))T }. ，. ．. (28). 答ベクトル. −1 ≥ σ l2 [rll (0) − rYTl RYY rYl ]. を用いて次式に求められる：.  E{ S1 ( k ) SlT ( k )}    L # RYl = (σ T ⊗ I 2 M +1 )   = ∑ σ j R jl  E{ S L ( k ) SlT ( k )}  j =1  . の非因果型フィルタを考える: M. Rlj ， (l , j = 1, 2," , L ). を用いた：. rYl = ( rYl ( M ) " rYl ( − M ))T ． (31) このとき， rYl は次式に示すとおり，入力信号と第 l 番目の要素となる音との相互相関行列 RYl , (l = 1, 2," , L ) の中央の列として求められる:. " rYl ( 2M )   rYl (0)   RYl = E{Y ( k ) S ( k )} =  # % # (32) ，  r ( −2M ) " r (0)  Yl  Yl  ここに，(8) 式より， RYl は音量係数ベクトル σ ，および，予め評価が可能な要素と T l. べき乗を重みとした時間平均で近似すれば，以下のとおり漸化式に帰着させて計算量を. 4. ⓒ 2011 Information Processing Society of Japan.

(5) Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. 低減できる:. リング周波数.  rYY (0) " rYY ( 2M )    T # % # RYY ( k ) =   = E{Y ( k )Y ( k )}  r ( 2 M ) " r ( 0)  YY  YY  1 = N. fs = 44.1 [kHz]. で発生した際の自己相関を，横軸として時点の差にと. って示している．これらは同一音高であるが音色の差異をもたらす自己相関は明瞭に異なる．なお，実際の楽器を用いても，自己相関には再現性があり，音色の差異に依存することを確認している．. −1.  ∞ j  ∞ j  T ≈ y ( k ) y ( k )  ∑ α   ∑ α y( k − j ) y ( k − j )  ∑ k =1  j =0   j =0  N. T. ∞. = (1 − α )∑ α j y( k − j ) y T ( k − j ) j =0. ∞. = α (1 − α )∑ α j y( k − 1 − j ) y T ( k − 1 − j ) + (1 − α ) y( k ) y T ( k ) j =0. = α RYY ( k − 1) + (1 − α ) y( k ) y T ( k ) ，ここに，. y(k ). (36). は次式の構成であり，入力信号のエネルギ. e. 2. は. RYY ( k ). の対角 Fig.1 Autocorrelation of violin tone (440Hz). 要素として評価できる：. y( k ) = ( y( k + M ) " y( k − M ) ). T. (37). ，. e = E{(Y ( k )) } = rYY (0) ， 2. 2. さらに,. −1 YY. R. (38). の算出は逆行列の補題[6]により漸化式に帰着できる：. R ( k ) = (α RYY ( k − 1) + (1 − α ) y( k ) y T ( k ))−1 −1 YY. =. 1. α. −1 RYY ( k − 1) −. −1 −1 (1 − α ) RYY ( k − 1) y( k ) y T ( k ) RYY ( k − 1) ． T −1 α [α + (1 − α ) y ( k ) RYY ( k − 1) y( k )]. (39). 5. 実施例図 1 から図 3 には，順に，バイオリン (Vn) 音，トランペット (Tp) 音，および，フルート (Fl) 音により，大譜表の中央の周波数. f 69 = 440 [Hz] ，を. C. 音の直ぐ上の. A. Fig.2 Autocorrelation of trumpet tone (440Hz). 音，すなわち，名目. MIDI システム (YAMAHA 社 XGworks) によりサンプ. 5. ⓒ 2011 Information Processing Society of Japan.

(6) Vol.2011-MUS-89 No.20 2011/2/12. 情報処理学会研究報告 IPSJ SIG Technical Report. 間をかけて，図 4 の左右とも Vn 音の音量を増加させ，同時に左側のケースでは Tp 音を，右側のケースでは Fl 音の音量を減少した場合である．両ケースを通じ，2 節に述べた手法が妥当に機能したことを例証できる．これらの同定結果を，さらに 3 節に述べた推定手法に適用することにより，Vn 音，Tp 音，Fl 音の混合について，互いに音高が異なる場合には極めて明瞭に，同一音高の場合でも，例えば，楽器パートごとの採譜に有効な程度に，それぞれの楽器音に分離できることを確認している．. 6. 結論合奏音のモノラル信号を指定された楽器パートごとの演奏音に分離する際に適用できる下記の特徴をもつ信号処理法を提案した： 1) 楽器パートごとの各音高を要素となる音に選び，事前に評価したこれらの自己相関，相互相関に基づき，統計的最小二乗規範を適用 2) 合奏音を各要素となる音に音量係数を乗じた和と仮定し，これに白色雑音を加えて入力信号としてモデル化 3) まず，各要素となる音の音量係数を同定した後，入力信号と各要素となる音と相互相関を評価して，楽器パートごとの演奏音を推定 4) バイオリン音，トランペット音，フルート音を混合し，同一音高でも分離が可能な信号処理性能を確認. Fig.3 Autocorrelation of flute tone (440Hz). 参考文献 [1] Y. Fukayama et al.: Identification of Instruments and Keys for Music on the Time-frequency plane; Proceedings of 41st international symposium on stochastic systems theory and its applications, pp.19-24, ISCIE, Kobe JP (2009) [2] MIDI manufacturers association; http://www.midi.org/aboutmidi/tutorials.php [3] R. A. Roberts and C. T. Mullis: Digital Signal Processing; Addison Wesley, Reading MA (1987) [4] M. Kennedy: Oxford dictionary of music; Oxford, New York NY (1985) [5] A. Klapuri and M. Davy, editors: Signal processing methods for music transcription, Springer, New York NY (2006) [6] G. Golub and C. Van Loan: Matrix computations; Johns Hopkins, Baltimore MA (1983). Fig.4 Identification of amplitude factor (440Hz) 図 4 には 2 節に述べた手法の適用例として，同一音高の 2 音の混合より音量係数を同定した結果を示す．これらは 4 節に従い， f s. = 44.1 [kHz] ， f 69 = 440 [Hz]. の音高で MIDI システム (YAMAHA 社 XGworks) を用いて Vn 音，Tp 音，Fl 音の統計量を事前に評価した後，別の MIDI システム (ROLAND 社 VSC) にて，すべて 2 秒. 6. ⓒ 2011 Information Processing Society of Japan.

(7)