JP2013250357A - 音響解析装置およびプログラム - Google Patents

Abstract

【課題】音響信号の相異なる音色に対応する調波成分を高精度に解析する。
【解決手段】全極型伝達関数１/|Ａ_f ^j|で表現されて相異なる音色の調波成分に対応するＪ個のスペクトル包絡ＶA_f ^jの各々と、ガウス関数列で表現されて相異なる基本周波数μ_n ^kに対応するＫ個の調波構造Ｇ_n,f ^kの各々との組合せに対応する(Ｊ×Ｋ)個の調波要素ＥA_n ^j,kを要素毎の音量Ｕ_n ^mで混合する音響モデルのスペクトログラムＸ_n,fが、音響信号ＳyのスペクトログラムＹ_n,fに近似するように、全極型伝達関数１/|Ａ_f ^j|の係数α_p ^jと、各調波要素ＥA_n ^j,kの音量Ｕ_n ^mと、各調波構造Ｇ_n,f ^kの基本周波数μ_n ^kとを、反復的な更新で推定する。
【選択図】図２

Description

本発明は、音響信号を解析する技術に関する。

音響信号を要素成分毎（例えば楽器毎）に分離する技術が従来から提案されている。例えば非特許文献１には、非負値行列因子分解（NMF：Non-negative Matrix Factorization）を利用した音源分離が開示されている。非負値行列因子分解を利用した音源分離では、音響信号の各成分の振幅スペクトルに対応する基底ベクトルを配列した基底行列と、各基底ベクトルの加重値の時間変化を示す係数行列とに音響信号が分解される。非特許文献２には、複数のガウス分布を周波数軸上に等間隔に配列した音響モデルを定義し、音響信号の振幅スペクトルを時刻毎に複数の音響モデルに分配する技術（ハーモニッククラスタリング）が開示されている。

P. Smaragdis, et. al., "Non-negative Matrix Factorization for Polyphonic Music Transcription", Proc. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, 2003, p. 170-180 H. Kameoka, et. al., "Extraction of Multiple Fundamental Frequencies from Polyphonic Music Using Harmonic Clustering", In Proceedings of 18th International Congress on Acoustics, 2004, p. I-59-62

非特許文献１の技術では、音色が共通で音高が相違する複数の音響（例えば１種類の楽器が発音した各音高の音響）が相異なる基底ベクトルに分離されるため、基底行列内の複数の基底ベクトルを音色毎（楽器毎）に正確に分類することが困難であるという問題がある。また、非特許文献２の技術では、音響信号の振幅スペクトルが時刻毎に独立に複数の音響モデルに分配されるから、時間的な変動が小さい音響特性（典型的には楽器毎の音色）を推定できず、非特許文献１と同様に、音響信号を音色毎に正確に分離することは困難である。以上の事情を考慮して、本発明は、音響信号の相異なる音色に対応する調波成分を高精度に解析することを目的とする。

以上の課題を解決するために本発明が採用する手段を説明する。なお、本発明の理解を容易にするために、以下の説明では、本発明の要素と後述の実施形態の要素との対応を括弧書で付記するが、本発明の範囲を実施形態の例示に限定する趣旨ではない。

本発明の音響解析装置は、第１全極型伝達関数（例えば全極型伝達関数１/|Ａ_f ^j|）で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡（例えばＪ個のスペクトル包絡ＶA_f ^j）の各々と、ガウス関数列で表現されて相異なる基本周波数（例えば基本周波数μ_n ^k）に対応する複数の調波構造（例えばＫ個の調波構造Ｇ_n,f ^k）の各々との組合せに対応する複数の調波要素（例えば調波要素ＥA_n ^j,k）を、要素毎の音量（例えば音量Ｕ_n ^m）で混合した音響モデルのスペクトログラム（例えばスペクトログラムＸ_n,f）が、対象音響信号（例えば音響信号Ｓy）のスペクトログラム（例えばスペクトログラムＹ_n,f）に近似するように、第１全極型伝達関数の係数（例えば係数α_p ^j）と各調波要素の音量と各調波構造の基本周波数とを反復的な更新で推定する変数解析手段を具備する。以上の構成によれば、調波成分に関連する各変数を高精度に解析することが可能である。なお、本発明の好適な態様において、調波成分に対応する各スペクトル包絡（調波成分の音色）は時不変とされる。以上の構成によれば、例えばガウス関数列を適用した時変のモデルで調波成分の各スペクトル包絡を表現した場合と比較して、調波成分のスペクトル包絡を高精度に推定できるという利点がある。

本発明の好適な態様において、音響モデルは、第２全極型伝達関数（例えば全極型伝達関数１/|Ｂ_f ^l|）でスペクトル包絡（例えばスペクトル包絡ＶB_f ^l）が表現されて相異なる音色に対応する複数の非調波要素（例えばＬ個の非調波要素ＥB^l）と複数の調波要素とを要素毎の音量で混合し、変数解析手段は、音響モデルのスペクトログラムと対象音響信号のスペクトログラムとが相互に近似するように、第１全極型伝達関数および第２全極型伝達関数の各係数と、各調波要素および各非調波要素の音量と、各調波構造の基本周波数とを、反復的な更新で推定する。以上の態様では、調波成分および非調波成分の双方について各変数を高精度に解析できるという利点がある。なお、本発明の好適な態様において、非調波要素に対応する各スペクトル包絡（非調波成分の音色）は時不変とされる。以上の構成によれば、例えばガウス関数列を適用した時変のモデルで非調波成分の各スペクトル包絡を表現した場合と比較して、非調波成分のスペクトル包絡を高精度に推定できるという利点がある。

本発明の好適な態様において、変数解析手段は、音響モデルのスペクトログラムと対象音響信号のスペクトログラムとの間のＩダイバージェンスが最小となるように音響モデルの各変数を推定する。

本発明の好適な態様において、変数解析手段は、複数の基本周波数の各々の初期化後に音響モデルの各変数の更新処理を反復し、更新処理の反復過程で閾値を下回る音量となった調波構造に対応する各変数の更新を以後の更新処理での更新対象から除外する。以上の態様では、閾値を下回る音量となった調波構造に対応する各変数の更新が以後の更新処理での更新対象から除外されるから、全部の調波構造について更新処理を最後まで継続する構成と比較して演算量が削減されるという利点がある。

本発明の好適な態様に係る音響解析装置は、第１全極型伝達関数で表現される調波成分のスペクトル包絡と、当該調波成分の基本周波数の時間変化と、第２全極型伝達関数で表現される非調波要素のスペクトル包絡と、当該非調波要素の音量の時間変化とを含む解析結果画像を表示装置に表示させる表示制御手段を具備する。以上の態様では、各調波成分の基本周波数（音高）の時間変化と各非調波成分の音量の時間変化とを利用者が視覚的に容易に把握できるという利点がある。

本発明の好適な態様に係る音響解析装置は、変数解析手段が解析した複数の音量のうち特定の要素成分に対応する音量を変更することで当該要素成分を抑圧するフィルタ（例えばフィルタ（例えばフィルタＦ_n,f）を設定するとともにフィルタを対象音響信号に作用させる信号処理手段を具備する。本発明の音響解析装置によれば、対象音響信号の各調波成分が高精度に解析されるから、変数解析手段による解析結果に応じたフィルタを対象音響信号に作用させることで、対象音響信号の要素成分を高精度に抑圧することが可能である。

以上の各態様に係る音響解析装置は、音響信号の解析に専用されるＤＳＰ（Digital Signal Processor）などのハードウェア（電子回路）によって実現されるほか、ＣＰＵ（Central Processing Unit）等の汎用の演算処理装置とプログラムとの協働によっても実現される。本発明のプログラムは、第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素を、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、第１全極型伝達関数の係数と各調波要素の音量と各調波構造の基本周波数とを反復的な更新で推定する解析処理をコンピュータに実行させる。以上のプログラムによれば、本発明の音響解析装置と同様の作用および効果が奏される。本発明のプログラムは、コンピュータが読取可能な記録媒体に格納された形態で提供されてコンピュータにインストールされるほか、通信網を介した配信の形態で提供されてコンピュータにインストールされる。

本発明のひとつの実施形態に係る音響解析装置のブロック図である。 音響モデルの説明図である。 変数解析部が実行する解析処理のフローチャートである。 解析結果画像の模式図である。 実施形態の効果の説明図である。

図１は、本発明の好適な実施形態に係る音響解析装置１００のブロック図である。本実施形態の音響解析装置１００は、音色が相違する複数の音響成分（調波成分および非調波成分）が混合された音響信号Ｓyを解析する信号処理装置であり、図１に示すように、演算処理装置１０と記憶装置１２と表示装置１４と入力装置１６と放音装置１８とを具備するコンピュータシステムで実現される。

演算処理装置１０は、記憶装置１２に格納されたプログラムＰGMを実行することで、音響信号Ｓyを解析するための複数の機能（周波数分析部２２，変数解析部２４，表示制御部２６，信号処理部２８）を実現する。なお、演算処理装置１０の各機能を複数の装置に分散した構成や、専用の電子回路（DSP）が一部の機能を実現する構成も採用され得る。

記憶装置１２は、演算処理装置１０が実行するプログラムＰGMや演算処理装置１０が使用する各種のデータを記憶する。半導体記録媒体や磁気記録媒体等の公知の記録媒体や複数種の記録媒体の組合せが記憶装置１２として任意に採用され得る。本実施形態の記憶装置１２は音響信号Ｓyを記憶する。なお、可搬型または内蔵型の記録媒体を再生する外部再生装置（図示略）から音響解析装置１００が音響信号Ｓyを取得することも可能である。

表示装置１４（例えば液晶表示パネル）は、演算処理装置１０による解析結果を表示する。入力装置１６は、利用者からの指示を受付ける機器であり、例えば複数の操作子を含んで構成される。放音装置１８（例えばスピーカやヘッドホン）は、演算処理装置１０から指示された音波を再生する。

周波数分析部２２は、音響信号ＳyのスペクトログラムＹ_n,fを算定する。スペクトログラムＹ_n,fは、時間軸上のフレーム毎に算定された振幅スペクトルの時系列である。記号ｎは、時間軸上に離散的に設定された任意の時点（フレームの番号）を意味し、記号ｆは、周波数軸上に離散的に設定された任意の周波数（周波数ビン）を意味する。スペクトログラムＹ_n,fの算定には、短時間フーリエ変換等の公知の周波数解析が任意に採用される。

本実施形態では、図２の音響モデルで生成されるスペクトログラムＸ_n,fを音響信号ＳyのスペクトログラムＹ_n,fのモデルとして想定する。図２に示すように、(Ｊ×Ｋ)個の調波要素ＥA_n ^j,kの各々を要素毎の音量Ｈ_n ^j,kに応じて調整するとともにＬ個の非調波要素ＥB^lの各々を要素毎の音量Ｉ_n ^lに応じて調整し、調整後の各調波要素ＥA_n ^j,kと調整後の各非調波要素ＥB^lと（(ＪＫ＋Ｌ)個）を加算する音響モデルでスペクトログラムＸ_n,fは表現される。

(Ｊ×Ｋ)個の調波要素ＥA_n ^j,kは、相異なる音色（例えば楽器毎）の調波成分に対応するＪ個のスペクトル包絡ＶA_f ^jの各々と、相異なる基本周波数（音高）μ_n ^kに対応するＫ個の調波構造Ｇ_n,f ^kの各々との(Ｊ×Ｋ)通りの組合せに対応する。１個のスペクトル包絡ＶA_f ^jは、例えば弦楽器や管楽器等の調波性の１種類の楽器が発音する調波音のスペクトルの包絡線に相当する。なお、本実施形態では、各調波成分のスペクトル包絡ＶA_f ^jが時間的に変動しない（すなわち各調波成分の音色が時不変である）と仮定する。他方、調波構造Ｇ_n,f ^kは、基本周波数μ_n ^kに対応する基音成分と基本周波数μ_n ^kの整数倍の周波数に対応する複数の倍音成分とを配列した系列であり、基本周波数μ_n ^kに応じて時刻ｎ毎に刻々と変動する。音量Ｈ_n ^j,kは、Ｊ個のうち第ｊ番目のスペクトル包絡ＶA_f ^jとＫ個のうち第ｋ番目の調波構造Ｇ_n,f ^kとの組合せに対応する調波要素ＥA_n ^j,kの音量（加重値）に相当し、時刻ｎ毎に刻々と変動する。

他方、Ｌ個の非調波要素ＥB^lは、相異なる音色の非調波成分に対応するＬ個のスペクトル包絡ＶB_f ^lに対応する。１個のスペクトル包絡ＶB_f ^lは、例えば打楽器等の非調波性の１種類の楽器が発音する非調波音のスペクトルの包絡線に相当する。調波成分のスペクトル包絡ＶA_f ^jと同様に、本実施形態では、各非調波成分のスペクトル包絡ＶB_f ^lが時間的に変動しない（すなわち各非調波成分の音色が時不変である）と仮定する。音量Ｉ_n ^lは、Ｌ個のうち第ｌ番目のスペクトル包絡ＶB_f ^lに対応する非調波要素ＥB^lの音量（加重値）に相当し、時刻ｎ毎に刻々と変動する。

以上の説明から理解されるように、図２の音響モデルで生成されるスペクトログラムＸ_n,fは以下の数式(1)で定義される。なお、数式(1)の記号「:＝」は定義を意味する。数式(1)の右辺の第１項が調波成分に対応し、第２項が非調波成分に対応する。

数式(1)の関数１/|Ａ_f ^j|は、第ｊ番目の調波成分のスペクトル包絡ＶA_f ^jをＰ個の係数α_p ^j（ｐ＝１〜Ｐ）に応じて表現する数式(2)の全極型伝達関数である。なお、記号ｉは虚数単位を意味する。また、記号ｆ'は、周波数（周波数ビン）ｆに対応する正規化角周波数を意味する。

同様に、数式(1)の関数１/|Ｂ_f ^l|は、第ｌ番目の非調波成分のスペクトル包絡ＶB_f ^lをＱ個の係数β_q ^l（ｑ＝１〜Ｑ）に応じて表現する数式(3)の全極型伝達関数である。係数α_p ^jの個数Ｐや係数β_q ^lの個数Ｑは例えば１０個程度に設定される。

数式(1)の調波構造Ｇ_n,f ^kは、基本周波数μ_n ^kの基音成分と基本周波数μ_n ^kの整数倍の周波数（ｈ×μ_n ^k）の各倍音成分とに対応するガウス分布（ガウス関数）を基本周波数μ_n ^kに応じた間隔で周波数軸上に配列したガウス関数列を意味する以下の数式(4)で表現される。

数式(4)の記号ｈは倍音成分の次数（整数）を意味し、記号σ²はガウス分布の分散を意味する。分散σ²は、例えば単一の所定値に設定される。数式(4)の調波構造Ｇ_n,f ^kによれば、基本周波数μ_n ^kに応じてガウス関数列が時刻ｎ毎に周波数軸上で伸縮されるから、ビブラート等の微細な音高の変動も適切に表現できる。

ところで、H. Kameoka, et. al., "Speech Spectrum Modeling for Joint Estimation of Spectral Envelope and Fundamental Frequency", IEEE Trans. on Audio, Speech and Language Processing, Vol. 18, No.6, p. 1507-1516, 2010（以下「非特許文献３」という）には、調波成分および非調波成分の双方をガウス関数列でモデル化する構成が開示されている。ガウス関数列（各ガウス分布の間隔）は音高に応じて刻々と変動する。すなわち、非特許文献３の構成では、調波成分および非調波成分の双方のスペクトル包絡が時間的に変動する（音色が時変である）ことが前提となる。他方、本実施形態では、全極型伝達関数１/|Ａ_f ^j|を適用した時不変のモデルで各調波成分のスペクトル包絡ＶA_f ^jが表現され、全極型伝達関数１/|Ｂ_f ^l|を適用した時不変のモデルで各非調波成分のスペクトル包絡ＶB_f ^lが表現される。全極型伝達関数は共鳴過程のモデルとして好適であり、かつ、音色（スペクトル包絡）が時不変であるという過程は現実の音響の傾向に充分に整合するから、本実施形態によれば、非特許文献３の構成と比較して、各調波成分のスペクトル包絡ＶA_f ^jや各非調波成分のスペクトル包絡ＶB_f ^lを高精度に推定できるという格別の効果が実現される。

説明の便宜のため、(Ｊ×Ｋ)個の調波要素ＥA_n ^j,kとＬ個の非調波要素ＥB^lとに対して図２の上方から下方に向けて通し番号（０,１,２,……,ＪＫ＋Ｌ−１）を付与し、任意の１個の要素を変数ｍ（ｍ＝０〜ＪＫ＋Ｌ−１）で表現したうえで、以下の数式(5)のように変数Ｗ_n,f ^mおよび変数Ｕ_n ^mを定義する。なお、数式(5)の記号modは剰余を意味し、記号〈 〉は床関数を意味する。

数式(5)の関係を利用すると、前掲の数式(1)は以下の数式(6)のように変形される。

数式(6)から理解されるように、音響モデルのスペクトログラムＸ_n,fは、各要素成分（各調波要素ＥA_n ^j,k，各非調波要素ＥB^l）に対応するＭ個（(ＪＫ＋Ｌ)個）のスペクトルパターンＷ_n,f ^mと各要素成分に対応するＭ個の時変な音量Ｕ_n ^mとで表現される。

図１の変数解析部２４は、数式(6)で表現される音響モデルのスペクトログラムＸ_n,fと周波数分析部２２が算定した音響信号ＳyのスペクトログラムＹ_n,fとが相互に近似するように音響モデルの各変数を推定する。具体的には、変数解析部２４は、各調波構造Ｇ_n,f ^kの基本周波数μ_n ^kと、各調波成分のスペクトル包絡ＶA_f ^jを表現する全極型伝達関数１/|Ａ_f ^j|の各係数α_p ^jと、各非調波成分のスペクトル包絡ＶB_f ^lを表現する全極型伝達関数１/|Ｂ_f ^l|の各係数β_q ^lと、各調波要素ＥA_n ^j,kおよび各非調波要素ＥB^lの音量Ｕ_n ^m（Ｈ_n ^j,k，Ｉ_n ^l）とを推定する。各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）は反復的な更新で推定される。

変数解析部２４による各変数の推定は、以下の数式(7)で表現されるように、スペクトログラムＸ_n,fとスペクトログラムＹ_n,fとの乖離の度合を表現する評価関数（距離規準）Ｑを各変数｛μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m｝に関して（w.r.t.：with respect to）最小化する最適化問題として定式化される。

本実施形態では、以下の数式(8)で表現されるように、スペクトログラムＸ_n,fとスペクトログラムＹ_n,fとのＩダイバージェンスを評価関数Ｑとして採用する。

＜Ｉダイバージェンスを規準とした全極型伝達関数の係数の推定＞
図２の音響モデルを評価する評価関数Ｑに数式(8)のＩダイバージェンスを適用する場合、全極型伝達関数（１/|Ａ_f ^j|，１/|Ｂ_f ^l|）の各係数（α_p ^j，β_q ^l）を推定するための更新式の導出が問題となる。そこで、変数解析部２４による具体的な処理の説明に先立ち、数式(9)で表現されるように、時間軸上の１個の時刻（したがって時刻ｎは省略される）での振幅スペクトルＹ_fを全極型伝達関数γ/|Ａ_f|で近似する場合を仮定して、全極型伝達関数γ/|Ａ_f|の係数α_pを推定するという小課題を便宜的に検討する。

数式(9)の記号「〜」は近似を意味する。また、数式(9)の記号γは、小課題の検討のために便宜的に導入した音量を意味する。振幅スペクトルＹ_fと全極型伝達関数γ/|Ａ_f|との乖離の度合をＩダイバージェンスで規定する評価関数Ｑは、以下の数式(10)で表現される。ただし、数式(10)では、係数α_pの推定に関係しない要素を省略した。

数式(10)の評価関数Ｑを最小化する係数α_pの更新式を検討する。仮に評価関数Ｑが係数α_pの２次形式であれば、評価関数Ｑの係数α_pによる偏微分がゼロになるときの係数α_pの数値が更新値となり、この条件から係数α_pの更新式を解析的に導出することが可能である。しかし、数式(10)で表現される評価関数Ｑは係数α_pの２次形式ではないから、更新式の解析的な導出は困難である。以上の事情を考慮して、係数α_pの２次形式で表現される適切な補助関数を設定する補助関数法を利用して係数α_pの更新式を導出する。

補助関数法は、補助変数ξに対する補助関数Ｑ⁺(θ,ξ)の最小値が本来の最小化の目的となる関数Ｑ(θ)に合致するように補助関数Ｑ⁺(θ,ξ)を設計し（Ｑ(θ)＝min Ｑ⁺(θ,ξ)）、補助関数Ｑ⁺(θ,ξ)について補助変数ξに関する最小化と本来の変数θに関する最小化とを反復することで間接的に本来の関数Ｑ(θ)を単調減少させる手法である。補助関数Ｑ⁺(θ,ξ)を最小にする変数θおよび変数ξの双方が解析的に解けるように補助関数Ｑ⁺(θ,ξ)を設計すれば、変数の推定は簡単化される。

数式(10)の括弧内の第１項の対数関数log|Ａ_f|の非線形性を解消するために以下の数式(11)を想定する。

数式(11)の右辺は、変数|Ａ_f|²が変数ρ_fとなる地点での接線に相当するから、変数ρ_fを補助変数とする補助関数として利用できる。数式(11)の等号が成立するのは、補助変数ρ_fが変数|Ａ_f|²に合致する場合（ρ_f←|Ａ_f|²）である。

次に、数式(10)の括弧内の第２項の逆数を解消するために、以下の数式(12)で表現されるように点τ_fを中心とする２次のテイラー近似を検討する。

数式(12)の右辺は目的関数１/|Ａ_f|を下回る可能性があるため、補助関数の要件を厳密には充足しないが、変数τ_fを変数|Ａ_f|に合致させれば凸関数に対するニュートン法と同形になるから、変数τ_fを補助変数と見做した効率的かつ安定的な最適化が可能である。

数式(11)および数式(12)を利用することで、数式(10)の評価関数Ｑに対する数式(13)の補助関数Ｑ⁺が導出される。なお、数式(13)の変数Ｃは、係数α_pを含まない要素を意味する。

数式(13)は、変数|Ａ_f|に対して線形であるが、係数α_pに関する２次形式には依然として到達していない。そこで、複素数の補助関数ω_fを変数|Ａ_f|に適用した以下の数式(14)を想定する。

数式(14)の記号Ｒe［ ］は実部を意味し、記号＊は複素共役を意味する。

数式(14)と前掲の数式(9)とを数式(13)に適用することで、係数α_pの２次形式で表現される数式(15)の補助関数Ｑ⁺⁺が導出される。

数式(15)を利用した係数α_pの更新を検討する。前述の３種類の補助変数（ρ_f，τ_f，ω_f）を数式(16)のように更新し、数式(15)を係数α_pで偏微分してゼロとすることで以下の数式(17)が導出される。

変数ｐのＰ個分を連立することで、振幅スペクトルＹ_fと全極型伝達関数γ/|Ａ_f|とのＩダイバージェンス（数式(10)の評価関数Ｑ）が最小化されるように全極型伝達関数γ/|Ａ_f|の係数α_pを更新する更新式(18)が導出される。

数式(18)は対称テプリッツ（Toeplitz）型の方程式であり、レビンソン-ダービン（Levinson-Durbin）アルゴリズムを利用することで高速に演算することが可能である。

以上の検討を踏まえて、図１の変数解析部２４が音響モデルの各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を推定するための更新式を検討する。

＜音量Ｕ_n ^m＞
評価関数Ｑを定義する数式(8)のうち括弧内の第１項の対数関数log（１/Ｘ_n,f）（＝−logＸ_n,f）に着目する。音響モデルのスペクトログラムＸ_n,fを表現する数式(6)を考慮すると、対数関数−logＸ_n,fは、対数関数が総和（Σ）を内包する形式であると理解できる。以上の形式を解消する（対数関数内から総和を除去する）ためにイェンゼン（Jensen）の不等式を適用すると、以下の数式(19)が導出される。

数式(19)の変数λ_n,f ^mは、任意の変数ｎ,ｆ,ｍについて正数であり（∀ｎ,ｆ,ｍ：λ_n,f ^m＞０）、任意の変数ｎおよびｆについて総和が１となる変数（∀ｎ,ｆ：Σλ_n,f ^m＝１）である。数式(19)で等号が成立する条件は、ラグランジュ（Lagrange）の未定乗数法を利用して導出される以下の数式(20)で表現される。

数式(19)を利用することで、数式(8)の評価関数Ｑに対する数式(21)の補助関数Ｑ⁺（対数関数が総和を内包しない形式）が導出される。記号Ｃは、音響モデルの変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を含まない要素を意味する。

数式(21)を音量Ｕ_n ^mで偏微分することで以下の数式(22)が導出される。

数式(22)をゼロとすることで、数式(8)の評価関数Ｑ（スペクトログラムＸ_n,fとスペクトログラムＹ_n,fとのＩダイバージェンス）が最小化されるように音量Ｕ_n ^mを更新する以下の更新式(23)が導出される。

＜全極型伝達関数の係数α_p ^jおよび係数β_q ^l＞
前掲の数式(21)を変形すると、各調波成分のスペクトル包絡ＶA_f ^jを表現する全極型伝達関数１/|Ａ_f ^j|の係数α_p ^jに関連する要素は以下の数式(24)で表現される。

数式(24)が、前述の小課題の検討で想定した数式(10)の右辺と類似する形式であることを考慮すると、数式(10)に対応する更新式(18)を流用することで係数α_p ^jの更新式が導出されると理解できる。すなわち、数式(10)の変数Ｙ_fを数式(24)の変数Σ_k,nＹ_n,fλ_n,f ^jK+kに対応させ、数式(10)の変数γを数式(24)の変数Σ_k,nＧ_n,f ^kＨ_m ^j,kに対応させて数式(18)を変形することで、数式(8)の評価関数Ｑが最小化されるように係数α_p ^jを更新する以下の更新式(25)が導出される。

同様に、数式(10)の変数Ｙ_fを変数Σ_nＹ_n,fλ_n,f ^jK+lに対応させ、数式(10)の変数γを変数Σ_nＩ_n ^lに対応させて数式(18)を変形することで、数式(8)の評価関数Ｑが最小化されるように係数β_q ^lを更新する以下の更新式(26)が導出される。

＜基本周波数μ_n ^k＞
各調波構造Ｇ_n,f ^kの基本周波数μ_n ^kの更新式を導出するために、前掲の数式(21)の第１項のみに着目する。すなわち、数式(21)の第２項Σ_m,n,fＷ_n,f ^mＵ_n ^mは、基本周波数μ_n ^kに対する依存が無視できるほど微小であると仮定して省略する。数式(21)の第１項のうち基本周波数μ_n ^kに関連する要素は以下の数式(27)で表現される。

数式(27)にイェンゼンの不等式を適用することで、以下の数式(28)が導出される。

数式(28)の変数φ_n,f ^h,kは、任意の変数ｈ,ｋ,ｎ,ｆについて正数であり（∀ｈ,ｋ,ｎ,ｆ：φ_n,f ^h,k＞０）、任意の変数ｎおよびｆについて総和が１となる変数（∀ｎ,ｆ：Σφ_n,f ^h,k＝１）である。数式(28)を利用することで、数式(8)の評価関数Ｑに対する数式(29)の補助関数Ｑ⁺が導出される。

数式(29)を基本周波数μ_n ^kで偏微分してゼロとすることで、数式(8)の評価関数Ｑが最小化されるように基本周波数μ_n ^kを更新する以下の更新式(30)が導出される。

本実施形態の変数解析部２４は、音量Ｕ_n ^mを更新する更新式(23)の演算と、係数α_p ^jを更新する更新式(25)の演算と、係数β_q ^lを更新する更新式(26)の演算と、基本周波数μ_n ^kを更新する更新式(30)の演算とを反復的に実行することで音響モデルの各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を推定する。具体的には、変数解析部２４は図３の解析処理を実行する。解析処理は、例えば入力装置１６に対する利用者からの指示を契機として実行される。図３の解析処理を開始すると、変数解析部２４は、音響モデルの各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を初期化する（ＳA）。各変数を初期化する具体的な方法は任意であるが、例えば以下に例示する方法が好適である。

変数解析部２４は、対数軸上で等間隔に配列するＫ個の周波数の各々を各調波構造Ｇ_n,f ^kの基本周波数μ_n ^kの初期値に設定する（ＳA1）。なお、基本周波数μ_n ^kの初期値が適切でない場合（音響信号Ｓyの実際の基本周波数との誤差が大きい場合）、音響信号Ｓyの実際の基本周波数の整数倍または整数分の一の周波数が基本周波数μ_n ^kと誤推定される可能性が高いという傾向がある。以上の傾向を考慮して、本実施形態では、調波構造Ｇ_n,f ^kの総数Ｋを、音響信号Ｓyの調波成分に想定される最大同時発音数と比較して充分に大きい数値に予備的に設定し、基本周波数μ_n ^kの初期値の妥当性が低いと各変数の更新の反復の過程で評価できる調波構造Ｇ_n,f ^kを更新対象から順次に除外する方法（後述のステップＳB6）を採用する。

変数解析部２４は、音響信号ＳyのスペクトログラムＹ_n,fのうちＪ個のフレームの振幅スペクトルを例えばランダムに選択し、各振幅スペクトルの包絡線を近似する全極型伝達関数の係数を音響モデルの係数α_p ^jの初期値に設定する（ＳA2）。同様に、変数解析部２４は、音響信号ＳyのスペクトログラムＹ_n,fのうちＬ個のフレームの振幅スペクトルを例えばランダムに選択し、各振幅スペクトルの包絡線を近似する全極型伝達関数の係数を音響モデルの係数β_q ^lの初期値に設定する（ＳA3）。また、変数解析部２４は、音量Ｕ_n ^mを非負の乱数値に初期化する（ＳA4）。なお、ステップＳA1からステップＳA4の順序は任意に変更される。

以上の手順で音響モデルの各変数を初期化すると、変数解析部２４は、音響信号ＳyのスペクトログラムＹ_n,fと各変数の現段階での数値とを適用した演算で各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を更新する更新処理ＳBを実行する。更新処理ＳBを開始すると、変数解析部２４は、数式(20)の演算で変数λ_n,f ^mを算定する（ＳB1）。そして、変数解析部２４は、更新式(23)の演算で音量Ｕ_n ^mを更新し（ＳB2）、更新式(30)の演算で基本周波数μ_n ^kを更新し（ＳB3）、更新式(25)の演算で係数α_p ^jを更新し（ＳB4）、更新式(26)の演算で係数β_q ^lを更新する（ＳB5）。なお、ステップＳB2からステップＳB5の順序は任意に変更される。

ステップＳA1で基本周波数μ_n ^kの初期値に選定されたＫ個の周波数のうち音響信号Ｓyに実際に包含される基本周波数から乖離した周波数に対応する音量Ｕ_n ^mは、ステップＳB2での更新毎に順次に減少するという傾向がある。以上の傾向を考慮して、変数解析部２４は、ステップＳB2での更新後の音量Ｕ_n ^mが所定の閾値を下回る調波構造Ｇ_n,f ^k（すなわち、基本周波数μ_n ^kの初期値の妥当性が低いと評価できる調波構造Ｇ_n,f ^k）に関連する変数（基本周波数μ_n ^kおよび音量Ｕ_n ^m）を、以後の更新処理ＳBでの更新対象から除外する（ＳB6）。すなわち、更新処理の反復過程で音量Ｕ_n ^mが閾値を下回った調波構造Ｇ_n,f ^kは音響モデルから除去される。したがって、Ｋ個の調波構造Ｇ_n,f ^kの全部について更新処理ＳBを最後まで継続する構成と比較して変数解析部２４の演算量が削減されるという利点がある。

変数解析部２４は、更新処理ＳBの反復を終了する条件（以下「反復停止条件」という）が成立したか否かを判定する（ＳC1）。例えば変数解析部２４は、現段階までの更新処理ＳBの反復回数が所定回数に到達した場合に反復停止条件が成立したと判定し、反復回数が所定回数を下回る場合には反復停止条件が成立していないと判定する。なお、反復停止条件の判定方法は任意である。例えば、音響モデルの各変数の収束の有無を評価（収束判定）することも可能である。すなわち、変数解析部２４は、各変数が収束した場合に反復停止条件が成立したと判定し、各変数が収束していない場合には反復停止条件が成立していないと判定する。各変数の収束判定には公知の技術が任意に採用される。

反復停止条件が成立していない場合（ＳC1：NO）、変数解析部２４は、直前の更新処理ＳBでの更新後の各変数を適用した更新処理ＳBを実行する。すなわち、反復停止条件が成立するまで更新処理ＳBが順次に実行されて各変数が累積的に更新される。他方、反復停止条件が成立した場合（ＳC1：YES）、変数解析部２４は、直前の更新処理ＳBでの更新後の各変数を最終的な解析結果として確定して記憶装置１２に格納する（ＳC2）。変数解析部２４が実行する解析処理の具体的な内容は以上の通りである。

図１の表示制御部２６は、変数解析部２４の解析結果に応じた画像（以下「解析結果画像」という）を生成して表示装置１４に表示させる。図４に例示されるように、本実施形態の解析結果画像５０は、複数の領域（ＤY，ＤX，ＤA1，ＤA2，ＤB1，ＤB2）を含んで構成される。領域ＤYと領域ＤXと領域ＤA2と領域ＤB2とは時間軸が共通する。

領域ＤYには、周波数分析部２２が算定した音響信号ＳyのスペクトログラムＹ_n,fが表示され、領域ＤXには、変数解析部２４が推定した各変数（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）で定義される音響モデルのスペクトログラムＸ_n,fが表示される。以上のようにスペクトログラムＹ_n,fとスペクトログラムＸ_n,fとが対比的に表示されるから、利用者は、変数解析部２４による解析の精度を視覚的に確認することが可能である。

領域ＤA1および領域ＤA2は、音響信号Ｓyの調波成分に関する解析結果を利用者に提示する画像領域である。領域ＤA1には、変数解析部２４が推定した係数α_p ^jに応じた全極型伝達関数１/|Ａ_f ^j|で表現される各調波成分のスペクトル包絡ＶA_f ^jが表示される。領域ＤA2には、変数解析部２４が調波構造Ｇ_n,f ^k毎に推定した各基本周波数μ_n ^kの時間的な変動（音高の時間軌跡）が表示される。すなわち、領域ＤA2は、縦軸が音高（基本周波数μ_n ^k）を示すピアノロール形式の画像である。利用者は、領域ＤA2を視認することで、各調波成分の音高の時間軌跡（例えば楽器毎の旋律）を直観的に把握することが可能である。なお、領域ＤA2内の各調波成分の音高の時間軌跡の表示態様（濃度や色彩等）を、各調波成分について推定された音量Ｕ_n ^mに応じて制御する（すなわち、各調波成分の音量Ｕ_n ^mを濃度や色彩で表現する）ことも可能である。

他方、領域ＤB1および領域ＤB2は、音響信号Ｓyの非調波成分に関する解析結果を利用者に提示する画像領域である。領域ＤB1には、変数解析部２４が推定した係数β_q ^lに応じた全極型伝達関数１/|Ｂ_f ^l|で表現される各非調波成分のスペクトル包絡ＶB_f ^lが表示される。領域ＤB2には、変数解析部２４が各非調波成分について推定した音量Ｕ_n ^m（すなわち図２の音量Ｉ_n ^l）の時間的な変動が非調波成分毎（非調波要素ＥB^l毎）に表示される。利用者は、領域ＤB2を視認することで、各非調波成分の発音の時点（例えば各打楽器の発音点）や、領域ＤA2内の各調波成分の基本周波数μ_n ^kとの時間的な関係を直観的に把握することが可能である。

図１の信号処理部２８は、変数解析部２４の解析結果（μ_n ^k，α_p ^j，β_q ^l，Ｕ_n ^m）を適用した信号処理（フィルタ処理）を音響信号Ｓyに対して実行することで音響信号Ｓzを生成する。本実施形態の信号処理部２８は、音響信号Ｓyのうち入力装置１６に対する利用者からの指示に応じた要素成分を抑圧した音響信号Ｓzを生成する。

具体的には、信号処理部２８は、周波数分析部２２が算定した音響信号ＳyのスペクトログラムＹ_n,fについて以下の数式(31)の演算を実行することで音響信号ＳzのスペクトログラムＺ_n,fを算定する。数式(31)の演算は、変数解析部２４の解析結果に応じたフィルタＦ_n,fを音響信号ＳyのスペクトログラムＹ_n,fに作用させる処理を意味する。

信号処理部２８は、数式(31)で算定されたスペクトログラムＺ_n,fを時間領域の音響信号Ｓzに変換する。例えば、信号処理部２８は、スペクトログラムＺ_n,fと音響信号Ｓyの位相スペクトログラムとを適用した短時間逆フーリエ変換で音響信号Ｓzを生成する。なお、公知の位相復元法で音響信号Ｓzを生成することも可能である。信号処理部２８が生成した音響信号Ｓzが放音装置１８に供給されて音波として再生される。

数式(31)のフィルタＦ_n,fは、以下の数式(32)で表現される。

数式(32)のフィルタＦ_n,fの分母は、音響モデルのスペクトログラムＸ_n,f（数式(6)）に相当する。他方、数式(32)の分子の変数ｕ_n ^mは、音響モデルにおけるＭ個（(ＪＫ＋Ｌ)個）の要素成分（調波要素ＥA_n ^j,kおよび非調波要素ＥB^l）の音量（以下「調整音量」という）に対応する。Ｍ個の調整音量ｕ_n ^mのうち利用者からの指示に応じた要素成分に対応する各調整音量ｕ_n ^mは所定値εに設定され、残余の各調整音量ｕ_n ^mは変数解析部２４が推定した音量Ｕ_n ^mに設定される。所定値εは例えばゼロ（またはゼロに近い正数）に設定される。以上の説明から理解されるように、数式(32)のフィルタＦ_n,fの分子は、音響モデルのスペクトログラムＸ_n,fのうち利用者からの指示に応じた特定の要素成分の音量Ｕ_n ^mを所定値εに変更したスペクトログラムに相当する。したがって、フィルタＦ_n,fを音響信号Ｓyに作用させる数式(31)の演算により、音響信号Ｓyから特定の要素成分を抑圧（除去）した音響信号Ｓzが生成される。

利用者は、音響信号Ｓyのうち所望の要素成分を入力装置１６の操作で指定することが可能である。例えばＪ個の調波成分のうち特定の調波成分を利用者が選択した場合、信号処理部２８は、利用者が選択した調波成分のスペクトル包絡ＶA_f ^jとＫ個の調波構造Ｇ_n,f ^kの各々との組合せに対応するＫ個の調整音量ｕ_n ^mを所定値εに設定し、残余（(Ｍ−Ｋ)個）の各調整音量ｕ_n ^mを音量Ｕ_n ^mに設定する。したがって、音響信号Ｓyのうち利用者が選択した調波成分（例えば特定の楽器の演奏音）を抑圧した音響信号Ｓzが生成される。

Ｋ個の調波構造Ｇ_n,f ^kのうち特定の調波構造Ｇ_n,f ^kを利用者が選択した場合、信号処理部２８は、利用者が選択した調波構造Ｇ_n,f ^kとＪ個のスペクトル包絡ＶA_f ^jの各々との組合せに対応するＪ個の調整音量ｕ_n ^mを所定値εに設定し、残余（(Ｍ−Ｊ)個）の各調整音量ｕ_n ^mを音量Ｕ_n ^mに設定する。したがって、音響信号Ｓyのうち利用者が選択した調波構造Ｇ_n,f ^kに対応する基本周波数μ_n ^kの調波成分（すなわち特定の音高）を抑圧した音響信号Ｓzが生成される。

また、Ｌ個の非調波成分のうち特定の非調波成分を利用者が選択した場合、信号処理部２８は、利用者が選択した非調波成分（非調波要素ＥB^l）に対応する調整音量ｕ_n ^mを所定値εに設定し、残余の各調整音量ｕ_n ^mを音量Ｕ_n ^mに設定する。したがって、音響信号Ｓyのうち利用者が選択した非調波成分（例えば特定の打楽器の演奏音）を抑圧した音響信号Ｓzが生成される。

図５は、以上に説明した音響解析装置１００による処理結果である。図５では、相異なる２種類の調波性の楽器の演奏音を含む音響信号Ｓy（Ｊ＝２，Ｌ＝０）を楽器毎に分離（一方を抑圧）した場合のＳＮ（Signal/Noise）比が、本実施形態の音響解析装置１００を利用した場合と、非負値行列因子分解（NMF）での分離結果をｋ-means法で楽器毎に分類した場合（以下「対比例」という）とについて対比的に図示されている。ＳＮ比が高いほど分離精度が高いことを意味する。評価用の音楽は、ＲＷＣ（Real World Computing） Music Databeseから選択されたクラシックおよびジャズの音楽である。本実施形態によれば、対比例と比較して音響信号Ｓyの各要素成分を高精度に分離できることが図５から理解される。

＜変形例＞
以上に例示した形態には様々な変形が加えられる。例えば、前述の形態では、Ｊ個の調波成分とＬ個の非調波成分とを含む音響モデルを例示したが、Ｌ個の非調波成分を省略することも可能である。

また、前述の形態では、変数解析部２４の解析結果を表示装置１４による表示と信号処理部２８による信号処理とに適用したが、変数解析部２４の解析結果の利用方法は任意である。例えば、音響信号Ｓyのうち特定の楽器に対応する調波成分の基本周波数μ_n ^kの解析結果からその楽器の楽譜を作成する構成（自動採譜）や、音響信号Ｓyの特定の要素成分を解析結果に応じて抽出して選択的に音響効果（例えば残響効果）を付与する構成も採用され得る。

１００……音響解析装置、１０……演算処理装置、１２……記憶装置、１４……表示装置、１６……入力装置、１８……放音装置、２２……周波数分析部、２４……変数解析部、２６……表示制御部、２８……信号処理部、５０……解析結果画像。

Claims (7)

1. 第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素を、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数の係数と前記各調波要素の音量と前記各調波構造の基本周波数とを反復的な更新で推定する変数解析手段
   を具備する音響解析装置。
2. 前記音響モデルは、第２全極型伝達関数でスペクトル包絡が表現されて相異なる音色に対応する複数の非調波要素と前記複数の調波要素とを要素毎の音量で混合し、
   前記変数解析手段は、前記音響モデルのスペクトログラムと前記対象音響信号のスペクトログラムとが相互に近似するように、前記第１全極型伝達関数および前記第２全極型伝達関数の各係数と、前記各調波要素および前記各非調波要素の音量と、前記各調波構造の基本周波数とを、反復的な更新で推定する
   請求項１の音響解析装置。
3. 前記調波成分に対応する各スペクトル包絡と前記非調波要素に対応する各スペクトル包絡とは時不変である
   請求項２の音響解析装置。
4. 前記変数解析手段は、前記音響モデルのスペクトログラムと前記対象音響信号のスペクトログラムとの間のＩダイバージェンスが最小となるように前記音響モデルの各変数を推定する
   請求項１から請求項３の何れかの音響解析装置。
5. 前記変数解析手段は、複数の基本周波数の各々の初期化後に前記音響モデルの各変数の更新処理を反復し、更新処理の反復過程で閾値を下回る音量となった調波構造に対応する各変数の更新を以後の更新処理での更新対象から除外する
   請求項１から請求項４の何れかの音響解析装置。
6. 前記第１全極型伝達関数で表現される調波成分のスペクトル包絡と、当該調波成分の基本周波数の時間変化と、前記第２全極型伝達関数で表現される非調波要素のスペクトル包絡と、当該非調波要素の音量の時間変化とを含む解析結果画像を表示装置に表示させる表示制御手段
   を具備する請求項１から請求項５の何れかの音響解析装置。
7. 第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素を、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数の係数と前記各調波要素の音量と前記各調波構造の基本周波数とを反復的な更新で推定する解析処理
   をコンピュータに実行させるプログラム。

Images