(57)【特許請求の範囲】
【請求項１】
第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素と、第２全極型伝達関数でスペクトル包絡が表現されて相異なる音色に対応する複数の非調波要素とを、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数および前記第２全極型伝達関数の各係数と、前記各調波要素および前記各非調波要素の音量と、前記各調波構造の基本周波数とを、反復的な更新で推定する変数解析手段
を具備する音響解析装置。
【請求項２】
前記調波成分に対応する各スペクトル包絡と前記非調波要素に対応する各スペクトル包絡とは時不変である
請求項１の音響解析装置。
【請求項３】
前記第１全極型伝達関数で表現される調波成分のスペクトル包絡と、当該調波成分の基本周波数の時間変化と、前記第２全極型伝達関数で表現される非調波要素のスペクトル包絡と、当該非調波要素の音量の時間変化とを含む解析結果画像を表示装置に表示させる表示制御手段
を具備する請求項１または請求項２の音響解析装置。
【請求項４】
第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素を、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数の係数と前記各調波要素の音量と前記各調波構造の基本周波数とを反復的な更新で推定する変数解析手段を具備し、
前記変数解析手段は、複数の基本周波数の各々の初期化後に前記音響モデルの各変数の更新処理を反復し、更新処理の反復過程で閾値を下回る音量となった調波構造に対応する各変数の更新を以後の更新処理での更新対象から除外する音響解析装置。
【請求項５】
第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素と、第２全極型伝達関数でスペクトル包絡が表現されて相異なる音色に対応する複数の非調波要素とを、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数および前記第２全極型伝達関数の各係数と、前記各調波要素および前記各非調波要素の音量と、前記各調波構造の基本周波数とを、反復的な更新で推定する解析処理
をコンピュータに実行させるプログラム。
【請求項６】
第１全極型伝達関数で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡の各々と、ガウス関数列で表現されて相異なる基本周波数に対応する複数の調波構造の各々との組合せに対応する複数の調波要素を、要素毎の音量で混合した音響モデルのスペクトログラムが、対象音響信号のスペクトログラムに近似するように、前記第１全極型伝達関数の係数と前記各調波要素の音量と前記各調波構造の基本周波数とを反復的な更新で推定する解析処理
をコンピュータに実行させるプログラムであって、
前記解析処理においては、複数の基本周波数の各々の初期化後に前記音響モデルの各変数の更新処理を反復し、更新処理の反復過程で閾値を下回る音量となった調波構造に対応する各変数の更新を以後の更新処理での更新対象から除外するプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音響信号を解析する技術に関する。
【背景技術】
【０００２】
音響信号を要素成分毎（例えば楽器毎）に分離する技術が従来から提案されている。例えば非特許文献１には、非負値行列因子分解（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
【発明の概要】
【発明が解決しようとする課題】
【０００４】
非特許文献１の技術では、音色が共通で音高が相違する複数の音響（例えば１種類の楽器が発音した各音高の音響）が相異なる基底ベクトルに分離されるため、基底行列内の複数の基底ベクトルを音色毎（楽器毎）に正確に分類することが困難であるという問題がある。また、非特許文献２の技術では、音響信号の振幅スペクトルが時刻毎に独立に複数の音響モデルに分配されるから、時間的な変動が小さい音響特性（典型的には楽器毎の音色）を推定できず、非特許文献１と同様に、音響信号を音色毎に正確に分離することは困難である。以上の事情を考慮して、本発明は、音響信号の相異なる音色に対応する調波成分を高精度に解析することを目的とする。
【課題を解決するための手段】
【０００５】
以上の課題を解決するために本発明が採用する手段を説明する。なお、本発明の理解を容易にするために、以下の説明では、本発明の要素と後述の実施形態の要素との対応を括弧書で付記するが、本発明の範囲を実施形態の例示に限定する趣旨ではない。
【０００６】
本発明の音響解析装置は、第１全極型伝達関数（例えば全極型伝達関数１/|Ａfj|）で表現されて相異なる音色の調波成分に対応する複数のスペクトル包絡（例えばＪ個のスペクトル包絡ＶAfj）の各々と、ガウス関数列で表現されて相異なる基本周波数（例えば基本周波数μnk）に対応する複数の調波構造（例えばＫ個の調波構造Ｇn,fk）の各々との組合せに対応する複数の調波要素（例えば調波要素ＥAnj,k）を、要素毎の音量（例えば音量Ｕnm）で混合した音響モデルのスペクトログラム（例えばスペクトログラムＸn,f）が、対象音響信号（例えば音響信号Ｓy）のスペクトログラム（例えばスペクトログラムＹn,f）に近似するように、第１全極型伝達関数の係数（例えば係数αpj）と各調波要素の音量と各調波構造の基本周波数とを反復的な更新で推定する変数解析手段を具備する。以上の構成によれば、調波成分に関連する各変数を高精度に解析することが可能である。なお、本発明の好適な態様において、調波成分に対応する各スペクトル包絡（調波成分の音色）は時不変とされる。以上の構成によれば、例えばガウス関数列を適用した時変のモデルで調波成分の各スペクトル包絡を表現した場合と比較して、調波成分のスペクトル包絡を高精度に推定できるという利点がある。
【０００７】
本発明の好適な態様において、音響モデルは、第２全極型伝達関数（例えば全極型伝達関数１/|Ｂfl|）でスペクトル包絡（例えばスペクトル包絡ＶBfl）が表現されて相異なる音色に対応する複数の非調波要素（例えばＬ個の非調波要素ＥBl）と複数の調波要素とを要素毎の音量で混合し、変数解析手段は、音響モデルのスペクトログラムと対象音響信号のスペクトログラムとが相互に近似するように、第１全極型伝達関数および第２全極型伝達関数の各係数と、各調波要素および各非調波要素の音量と、各調波構造の基本周波数とを、反復的な更新で推定する。以上の態様では、調波成分および非調波成分の双方について各変数を高精度に解析できるという利点がある。なお、本発明の好適な態様において、非調波要素に対応する各スペクトル包絡（非調波成分の音色）は時不変とされる。以上の構成によれば、例えばガウス関数列を適用した時変のモデルで非調波成分の各スペクトル包絡を表現した場合と比較して、非調波成分のスペクトル包絡を高精度に推定できるという利点がある。
【０００８】
本発明の好適な態様において、変数解析手段は、音響モデルのスペクトログラムと対象音響信号のスペクトログラムとの間のＩダイバージェンスが最小となるように音響モデルの各変数を推定する。
【０００９】
本発明の好適な態様において、変数解析手段は、複数の基本周波数の各々の初期化後に音響モデルの各変数の更新処理を反復し、更新処理の反復過程で閾値を下回る音量となった調波構造に対応する各変数の更新を以後の更新処理での更新対象から除外する。以上の態様では、閾値を下回る音量となった調波構造に対応する各変数の更新が以後の更新処理での更新対象から除外されるから、全部の調波構造について更新処理を最後まで継続する構成と比較して演算量が削減されるという利点がある。
【００１０】
本発明の好適な態様に係る音響解析装置は、第１全極型伝達関数で表現される調波成分のスペクトル包絡と、当該調波成分の基本周波数の時間変化と、第２全極型伝達関数で表現される非調波要素のスペクトル包絡と、当該非調波要素の音量の時間変化とを含む解析結果画像を表示装置に表示させる表示制御手段を具備する。以上の態様では、各調波成分の基本周波数（音高）の時間変化と各非調波成分の音量の時間変化とを利用者が視覚的に容易に把握できるという利点がある。
【００１１】
本発明の好適な態様に係る音響解析装置は、変数解析手段が解析した複数の音量のうち特定の要素成分に対応する音量を変更することで当該要素成分を抑圧するフィルタ（例えばフィルタ（例えばフィルタＦ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のモデルとして想定する。図２に示すように、(Ｊ×Ｋ)個の調波要素ＥAnj,kの各々を要素毎の音量Ｈnj,kに応じて調整するとともにＬ個の非調波要素ＥBlの各々を要素毎の音量Ｉnlに応じて調整し、調整後の各調波要素ＥAnj,kと調整後の各非調波要素ＥBlと（(ＪＫ＋Ｌ)個）を加算する音響モデルでスペクトログラムＸn,fは表現される。
【００２０】
(Ｊ×Ｋ)個の調波要素ＥAnj,kは、相異なる音色（例えば楽器毎）の調波成分に対応するＪ個のスペクトル包絡ＶAfjの各々と、相異なる基本周波数（音高）μnkに対応するＫ個の調波構造Ｇn,fkの各々との(Ｊ×Ｋ)通りの組合せに対応する。１個のスペクトル包絡ＶAfjは、例えば弦楽器や管楽器等の調波性の１種類の楽器が発音する調波音のスペクトルの包絡線に相当する。なお、本実施形態では、各調波成分のスペクトル包絡ＶAfjが時間的に変動しない（すなわち各調波成分の音色が時不変である）と仮定する。他方、調波構造Ｇn,fkは、基本周波数μnkに対応する基音成分と基本周波数μnkの整数倍の周波数に対応する複数の倍音成分とを配列した系列であり、基本周波数μnkに応じて時刻ｎ毎に刻々と変動する。音量Ｈnj,kは、Ｊ個のうち第ｊ番目のスペクトル包絡ＶAfjとＫ個のうち第ｋ番目の調波構造Ｇn,fkとの組合せに対応する調波要素ＥAnj,kの音量（加重値）に相当し、時刻ｎ毎に刻々と変動する。
【００２１】
他方、Ｌ個の非調波要素ＥBlは、相異なる音色の非調波成分に対応するＬ個のスペクトル包絡ＶBflに対応する。１個のスペクトル包絡ＶBflは、例えば打楽器等の非調波性の１種類の楽器が発音する非調波音のスペクトルの包絡線に相当する。調波成分のスペクトル包絡ＶAfjと同様に、本実施形態では、各非調波成分のスペクトル包絡ＶBflが時間的に変動しない（すなわち各非調波成分の音色が時不変である）と仮定する。音量Ｉnlは、Ｌ個のうち第ｌ番目のスペクトル包絡ＶBflに対応する非調波要素ＥBlの音量（加重値）に相当し、時刻ｎ毎に刻々と変動する。
【００２２】
以上の説明から理解されるように、図２の音響モデルで生成されるスペクトログラムＸn,fは以下の数式(1)で定義される。なお、数式(1)の記号「:＝」は定義を意味する。数式(1)の右辺の第１項が調波成分に対応し、第２項が非調波成分に対応する。
【数１】
0006044119.tif 000002
【００２３】
数式(1)の関数１/|Ａfj|は、第ｊ番目の調波成分のスペクトル包絡ＶAfjをＰ個の係数αpj（ｐ＝１〜Ｐ）に応じて表現する数式(2)の全極型伝達関数である。なお、記号ｉは虚数単位を意味する。また、記号ｆ'は、周波数（周波数ビン）ｆに対応する正規化角周波数を意味する。
【数２】
0006044119.tif 000003
【００２４】
同様に、数式(1)の関数１/|Ｂfl|は、第ｌ番目の非調波成分のスペクトル包絡ＶBflをＱ個の係数βql（ｑ＝１〜Ｑ）に応じて表現する数式(3)の全極型伝達関数である。係数αpjの個数Ｐや係数βqlの個数Ｑは例えば１０個程度に設定される。
【数３】
0006044119.tif 000004
【００２５】
数式(1)の調波構造Ｇn,fkは、基本周波数μnkの基音成分と基本周波数μnkの整数倍の周波数（ｈ×μnk）の各倍音成分とに対応するガウス分布（ガウス関数）を基本周波数μnkに応じた間隔で周波数軸上に配列したガウス関数列を意味する以下の数式(4)で表現される。
【数４】
0006044119.tif 000005
数式(4)の記号ｈは倍音成分の次数（整数）を意味し、記号σ2はガウス分布の分散を意味する。分散σ2は、例えば単一の所定値に設定される。数式(4)の調波構造Ｇn,fkによれば、基本周波数μnkに応じてガウス関数列が時刻ｎ毎に周波数軸上で伸縮されるから、ビブラート等の微細な音高の変動も適切に表現できる。
【００２６】
ところで、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（以下「非特許文献３」という）には、調波成分および非調波成分の双方をガウス関数列でモデル化する構成が開示されている。ガウス関数列（各ガウス分布の間隔）は音高に応じて刻々と変動する。すなわち、非特許文献３の構成では、調波成分および非調波成分の双方のスペクトル包絡が時間的に変動する（音色が時変である）ことが前提となる。他方、本実施形態では、全極型伝達関数１/|Ａfj|を適用した時不変のモデルで各調波成分のスペクトル包絡ＶAfjが表現され、全極型伝達関数１/|Ｂfl|を適用した時不変のモデルで各非調波成分のスペクトル包絡ＶBflが表現される。全極型伝達関数は共鳴過程のモデルとして好適であり、かつ、音色（スペクトル包絡）が時不変であるという過程は現実の音響の傾向に充分に整合するから、本実施形態によれば、非特許文献３の構成と比較して、各調波成分のスペクトル包絡ＶAfjや各非調波成分のスペクトル包絡ＶBflを高精度に推定できるという格別の効果が実現される。
【００２７】
説明の便宜のため、(Ｊ×Ｋ)個の調波要素ＥAnj,kとＬ個の非調波要素ＥBlとに対して図２の上方から下方に向けて通し番号（０,１,２,……,ＪＫ＋Ｌ−１）を付与し、任意の１個の要素を変数ｍ（ｍ＝０〜ＪＫ＋Ｌ−１）で表現したうえで、以下の数式(5)のように変数Ｗn,fmおよび変数Ｕnmを定義する。なお、数式(5)の記号modは剰余を意味し、記号〈 〉は床関数を意味する。
【数５】
0006044119.tif 000006
【００２８】
数式(5)の関係を利用すると、前掲の数式(1)は以下の数式(6)のように変形される。
【数６】
0006044119.tif 000007
数式(6)から理解されるように、音響モデルのスペクトログラムＸn,fは、各要素成分（各調波要素ＥAnj,k，各非調波要素ＥBl）に対応するＭ個（(ＪＫ＋Ｌ)個）のスペクトルパターンＷn,fmと各要素成分に対応するＭ個の時変な音量Ｕnmとで表現される。
【００２９】
図１の変数解析部２４は、数式(6)で表現される音響モデルのスペクトログラムＸn,fと周波数分析部２２が算定した音響信号ＳyのスペクトログラムＹn,fとが相互に近似するように音響モデルの各変数を推定する。具体的には、変数解析部２４は、各調波構造Ｇn,fkの基本周波数μnkと、各調波成分のスペクトル包絡ＶAfjを表現する全極型伝達関数１/|Ａfj|の各係数αpjと、各非調波成分のスペクトル包絡ＶBflを表現する全極型伝達関数１/|Ｂfl|の各係数βqlと、各調波要素ＥAnj,kおよび各非調波要素ＥBlの音量Ｕnm（Ｈnj,k，Ｉnl）とを推定する。各変数（μnk，αpj，βql，Ｕnm）は反復的な更新で推定される。
【００３０】
変数解析部２４による各変数の推定は、以下の数式(7)で表現されるように、スペクトログラムＸn,fとスペクトログラムＹn,fとの乖離の度合を表現する評価関数（距離規準）Ｑを各変数｛μnk，αpj，βql，Ｕnm｝に関して（w.r.t.：with respect to）最小化する最適化問題として定式化される。
【数７】
0006044119.tif 000008
【００３１】
本実施形態では、以下の数式(8)で表現されるように、スペクトログラムＸn,fとスペクトログラムＹn,fとのＩダイバージェンスを評価関数Ｑとして採用する。
【数８】
0006044119.tif 000009
【００３２】
＜Ｉダイバージェンスを規準とした全極型伝達関数の係数の推定＞
図２の音響モデルを評価する評価関数Ｑに数式(8)のＩダイバージェンスを適用する場合、全極型伝達関数（１/|Ａfj|，１/|Ｂfl|）の各係数（αpj，βql）を推定するための更新式の導出が問題となる。そこで、変数解析部２４による具体的な処理の説明に先立ち、数式(9)で表現されるように、時間軸上の１個の時刻（したがって時刻ｎは省略される）での振幅スペクトルＹfを全極型伝達関数γ/|Ａf|で近似する場合を仮定して、全極型伝達関数γ/|Ａf|の係数αpを推定するという小課題を便宜的に検討する。
【数９】
0006044119.tif 000010
数式(9)の記号「〜」は近似を意味する。また、数式(9)の記号γは、小課題の検討のために便宜的に導入した音量を意味する。振幅スペクトルＹfと全極型伝達関数γ/|Ａf|との乖離の度合をＩダイバージェンスで規定する評価関数Ｑは、以下の数式(10)で表現される。ただし、数式(10)では、係数αpの推定に関係しない要素を省略した。
【数１０】
0006044119.tif 000011
【００３３】
数式(10)の評価関数Ｑを最小化する係数αpの更新式を検討する。仮に評価関数Ｑが係数αpの２次形式であれば、評価関数Ｑの係数αpによる偏微分がゼロになるときの係数αpの数値が更新値となり、この条件から係数αpの更新式を解析的に導出することが可能である。しかし、数式(10)で表現される評価関数Ｑは係数αpの２次形式ではないから、更新式の解析的な導出は困難である。以上の事情を考慮して、係数αpの２次形式で表現される適切な補助関数を設定する補助関数法を利用して係数αpの更新式を導出する。
【００３４】
補助関数法は、補助変数ξに対する補助関数Ｑ+(θ,ξ)の最小値が本来の最小化の目的となる関数Ｑ(θ)に合致するように補助関数Ｑ+(θ,ξ)を設計し（Ｑ(θ)＝min Ｑ+(θ,ξ)）、補助関数Ｑ+(θ,ξ)について補助変数ξに関する最小化と本来の変数θに関する最小化とを反復することで間接的に本来の関数Ｑ(θ)を単調減少させる手法である。補助関数Ｑ+(θ,ξ)を最小にする変数θおよび変数ξの双方が解析的に解けるように補助関数Ｑ+(θ,ξ)を設計すれば、変数の推定は簡単化される。
【００３５】
数式(10)の括弧内の第１項の対数関数log|Ａf|の非線形性を解消するために以下の数式(11)を想定する。
【数１１】
0006044119.tif 000012
数式(11)の右辺は、変数|Ａf|2が変数ρfとなる地点での接線に相当するから、変数ρfを補助変数とする補助関数として利用できる。数式(11)の等号が成立するのは、補助変数ρfが変数|Ａf|2に合致する場合（ρf←|Ａf|2）である。
【００３６】
次に、数式(10)の括弧内の第２項の逆数を解消するために、以下の数式(12)で表現されるように点τfを中心とする２次のテイラー近似を検討する。
【数１２】
0006044119.tif 000013
数式(12)の右辺は目的関数１/|Ａf|を下回る可能性があるため、補助関数の要件を厳密には充足しないが、変数τfを変数|Ａf|に合致させれば凸関数に対するニュートン法と同形になるから、変数τfを補助変数と見做した効率的かつ安定的な最適化が可能である。
【００３７】
数式(11)および数式(12)を利用することで、数式(10)の評価関数Ｑに対する数式(13)の補助関数Ｑ+が導出される。なお、数式(13)の変数Ｃは、係数αpを含まない要素を意味する。
【数１３】
0006044119.tif 000014
【００３８】
数式(13)は、変数|Ａf|に対して線形であるが、係数αpに関する２次形式には依然として到達していない。そこで、複素数の補助関数ωfを変数|Ａf|に適用した以下の数式(14)を想定する。
【数１４】
0006044119.tif 000015
数式(14)の記号Ｒe［ ］は実部を意味し、記号＊は複素共役を意味する。
【００３９】
数式(14)と前掲の数式(9)とを数式(13)に適用することで、係数αpの２次形式で表現される数式(15)の補助関数Ｑ++が導出される。
【数１５】
0006044119.tif 000016
【００４０】
数式(15)を利用した係数αpの更新を検討する。前述の３種類の補助変数（ρf，τf，ωf）を数式(16)のように更新し、数式(15)を係数αpで偏微分してゼロとすることで以下の数式(17)が導出される。
【数１６】
0006044119.tif 000017
【数１７】
0006044119.tif 000018
【００４１】
変数ｐのＰ個分を連立することで、振幅スペクトルＹfと全極型伝達関数γ/|Ａf|とのＩダイバージェンス（数式(10)の評価関数Ｑ）が最小化されるように全極型伝達関数γ/|Ａf|の係数αpを更新する更新式(18)が導出される。
【数１８】
0006044119.tif 000019
数式(18)は対称テプリッツ（Toeplitz）型の方程式であり、レビンソン-ダービン（Levinson-Durbin）アルゴリズムを利用することで高速に演算することが可能である。
【００４２】
以上の検討を踏まえて、図１の変数解析部２４が音響モデルの各変数（μnk，αpj，βql，Ｕnm）を推定するための更新式を検討する。
【００４３】
＜音量Ｕnm＞
評価関数Ｑを定義する数式(8)のうち括弧内の第１項の対数関数log（１/Ｘn,f）（＝−logＸn,f）に着目する。音響モデルのスペクトログラムＸn,fを表現する数式(6)を考慮すると、対数関数−logＸn,fは、対数関数が総和（Σ）を内包する形式であると理解できる。以上の形式を解消する（対数関数内から総和を除去する）ためにイェンゼン（Jensen）の不等式を適用すると、以下の数式(19)が導出される。
【数１９】
0006044119.tif 000020
数式(19)の変数λn,fmは、任意の変数ｎ,ｆ,ｍについて正数であり（∀ｎ,ｆ,ｍ：λn,fm＞０）、任意の変数ｎおよびｆについて総和が１となる変数（∀ｎ,ｆ：Σλn,fm＝１）である。数式(19)で等号が成立する条件は、ラグランジュ（Lagrange）の未定乗数法を利用して導出される以下の数式(20)で表現される。
【数２０】
0006044119.tif 000021
【００４４】
数式(19)を利用することで、数式(8)の評価関数Ｑに対する数式(21)の補助関数Ｑ+（対数関数が総和を内包しない形式）が導出される。記号Ｃは、音響モデルの変数（μnk，αpj，βql，Ｕnm）を含まない要素を意味する。
【数２１】
0006044119.tif 000022
【００４５】
数式(21)を音量Ｕnmで偏微分することで以下の数式(22)が導出される。
【数２２】
0006044119.tif 000023
数式(22)をゼロとすることで、数式(8)の評価関数Ｑ（スペクトログラムＸn,fとスペクトログラムＹn,fとのＩダイバージェンス）が最小化されるように音量Ｕnmを更新する以下の更新式(23)が導出される。
【数２３】
0006044119.tif 000024
【００４６】
＜全極型伝達関数の係数αpjおよび係数βql＞
前掲の数式(21)を変形すると、各調波成分のスペクトル包絡ＶAfjを表現する全極型伝達関数１/|Ａfj|の係数αpjに関連する要素は以下の数式(24)で表現される。
【数２４】
0006044119.tif 000025
【００４７】
数式(24)が、前述の小課題の検討で想定した数式(10)の右辺と類似する形式であることを考慮すると、数式(10)に対応する更新式(18)を流用することで係数αpjの更新式が導出されると理解できる。すなわち、数式(10)の変数Ｙfを数式(24)の変数Σk,nＹn,fλn,fjK+kに対応させ、数式(10)の変数γを数式(24)の変数Σk,nＧn,fkＨmj,kに対応させて数式(18)を変形することで、数式(8)の評価関数Ｑが最小化されるように係数αpjを更新する以下の更新式(25)が導出される。
【数２５】
0006044119.tif 000026
【００４８】
同様に、数式(10)の変数Ｙfを変数ΣnＹn,fλn,fjK+lに対応させ、数式(10)の変数γを変数ΣnＩnlに対応させて数式(18)を変形することで、数式(8)の評価関数Ｑが最小化されるように係数βqlを更新する以下の更新式(26)が導出される。
【数２６】
0006044119.tif 000027
【００４９】
＜基本周波数μnk＞
各調波構造Ｇn,fkの基本周波数μnkの更新式を導出するために、前掲の数式(21)の第１項のみに着目する。すなわち、数式(21)の第２項Σm,n,fＷn,fmＵnmは、基本周波数μnkに対する依存が無視できるほど微小であると仮定して省略する。数式(21)の第１項のうち基本周波数μnkに関連する要素は以下の数式(27)で表現される。
【数２７】
0006044119.tif 000028
【００５０】
数式(27)にイェンゼンの不等式を適用することで、以下の数式(28)が導出される。
【数２８】
0006044119.tif 000029
【００５１】
数式(28)の変数φn,fh,kは、任意の変数ｈ,ｋ,ｎ,ｆについて正数であり（∀ｈ,ｋ,ｎ,ｆ：φn,fh,k＞０）、任意の変数ｎおよびｆについて総和が１となる変数（∀ｎ,ｆ：Σφn,fh,k＝１）である。数式(28)を利用することで、数式(8)の評価関数Ｑに対する数式(29)の補助関数Ｑ+が導出される。
【数２９】
0006044119.tif 000030
【００５２】
数式(29)を基本周波数μnkで偏微分してゼロとすることで、数式(8)の評価関数Ｑが最小化されるように基本周波数μnkを更新する以下の更新式(30)が導出される。
【数３０】
0006044119.tif 000031
【００５３】
本実施形態の変数解析部２４は、音量Ｕnmを更新する更新式(23)の演算と、係数αpjを更新する更新式(25)の演算と、係数βqlを更新する更新式(26)の演算と、基本周波数μnkを更新する更新式(30)の演算とを反復的に実行することで音響モデルの各変数（μnk，αpj，βql，Ｕnm）を推定する。具体的には、変数解析部２４は図３の解析処理を実行する。解析処理は、例えば入力装置１６に対する利用者からの指示を契機として実行される。図３の解析処理を開始すると、変数解析部２４は、音響モデルの各変数（μnk，αpj，βql，Ｕnm）を初期化する（ＳA）。各変数を初期化する具体的な方法は任意であるが、例えば以下に例示する方法が好適である。
【００５４】
変数解析部２４は、対数軸上で等間隔に配列するＫ個の周波数の各々を各調波構造Ｇn,fkの基本周波数μnkの初期値に設定する（ＳA1）。なお、基本周波数μnkの初期値が適切でない場合（音響信号Ｓyの実際の基本周波数との誤差が大きい場合）、音響信号Ｓyの実際の基本周波数の整数倍または整数分の一の周波数が基本周波数μnkと誤推定される可能性が高いという傾向がある。以上の傾向を考慮して、本実施形態では、調波構造Ｇn,fkの総数Ｋを、音響信号Ｓyの調波成分に想定される最大同時発音数と比較して充分に大きい数値に予備的に設定し、基本周波数μnkの初期値の妥当性が低いと各変数の更新の反復の過程で評価できる調波構造Ｇn,fkを更新対象から順次に除外する方法（後述のステップＳB6）を採用する。
【００５５】
変数解析部２４は、音響信号ＳyのスペクトログラムＹn,fのうちＪ個のフレームの振幅スペクトルを例えばランダムに選択し、各振幅スペクトルの包絡線を近似する全極型伝達関数の係数を音響モデルの係数αpjの初期値に設定する（ＳA2）。同様に、変数解析部２４は、音響信号ＳyのスペクトログラムＹn,fのうちＬ個のフレームの振幅スペクトルを例えばランダムに選択し、各振幅スペクトルの包絡線を近似する全極型伝達関数の係数を音響モデルの係数βqlの初期値に設定する（ＳA3）。また、変数解析部２４は、音量Ｕnmを非負の乱数値に初期化する（ＳA4）。なお、ステップＳA1からステップＳA4の順序は任意に変更される。
【００５６】
以上の手順で音響モデルの各変数を初期化すると、変数解析部２４は、音響信号ＳyのスペクトログラムＹn,fと各変数の現段階での数値とを適用した演算で各変数（μnk，αpj，βql，Ｕnm）を更新する更新処理ＳBを実行する。更新処理ＳBを開始すると、変数解析部２４は、数式(20)の演算で変数λn,fmを算定する（ＳB1）。そして、変数解析部２４は、更新式(23)の演算で音量Ｕnmを更新し（ＳB2）、更新式(30)の演算で基本周波数μnkを更新し（ＳB3）、更新式(25)の演算で係数αpjを更新し（ＳB4）、更新式(26)の演算で係数βqlを更新する（ＳB5）。なお、ステップＳB2からステップＳB5の順序は任意に変更される。
【００５７】
ステップＳA1で基本周波数μnkの初期値に選定されたＫ個の周波数のうち音響信号Ｓyに実際に包含される基本周波数から乖離した周波数に対応する音量Ｕnmは、ステップＳB2での更新毎に順次に減少するという傾向がある。以上の傾向を考慮して、変数解析部２４は、ステップＳB2での更新後の音量Ｕnmが所定の閾値を下回る調波構造Ｇn,fk（すなわち、基本周波数μnkの初期値の妥当性が低いと評価できる調波構造Ｇn,fk）に関連する変数（基本周波数μnkおよび音量Ｕnm）を、以後の更新処理ＳBでの更新対象から除外する（ＳB6）。すなわち、更新処理の反復過程で音量Ｕnmが閾値を下回った調波構造Ｇn,fkは音響モデルから除去される。したがって、Ｋ個の調波構造Ｇn,fkの全部について更新処理Ｓ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には、変数解析部２４が推定した各変数（μnk，αpj，βql，Ｕnm）で定義される音響モデルのスペクトログラムＸn,fが表示される。以上のようにスペクトログラムＹn,fとスペクトログラムＸn,fとが対比的に表示されるから、利用者は、変数解析部２４による解析の精度を視覚的に確認することが可能である。
【００６２】
領域ＤA1および領域ＤA2は、音響信号Ｓyの調波成分に関する解析結果を利用者に提示する画像領域である。領域ＤA1には、変数解析部２４が推定した係数αpjに応じた全極型伝達関数１/|Ａfj|で表現される各調波成分のスペクトル包絡ＶAfjが表示される。領域ＤA2には、変数解析部２４が調波構造Ｇn,fk毎に推定した各基本周波数μnkの時間的な変動（音高の時間軌跡）が表示される。すなわち、領域ＤA2は、縦軸が音高（基本周波数μnk）を示すピアノロール形式の画像である。利用者は、領域ＤA2を視認することで、各調波成分の音高の時間軌跡（例えば楽器毎の旋律）を直観的に把握することが可能である。なお、領域ＤA2内の各調波成分の音高の時間軌跡の表示態様（濃度や色彩等）を、各調波成分について推定された音量Ｕnmに応じて制御する（すなわち、各調波成分の音量Ｕnmを濃度や色彩で表現する）ことも可能である。
【００６３】
他方、領域ＤB1および領域ＤB2は、音響信号Ｓyの非調波成分に関する解析結果を利用者に提示する画像領域である。領域ＤB1には、変数解析部２４が推定した係数βqlに応じた全極型伝達関数１/|Ｂfl|で表現される各非調波成分のスペクトル包絡ＶBflが表示される。領域ＤB2には、変数解析部２４が各非調波成分について推定した音量Ｕnm（すなわち図２の音量Ｉnl）の時間的な変動が非調波成分毎（非調波要素ＥBl毎）に表示される。利用者は、領域ＤB2を視認することで、各非調波成分の発音の時点（例えば各打楽器の発音点）や、領域ＤA2内の各調波成分の基本周波数μnkとの時間的な関係を直観的に把握することが可能である。
【００６４】
図１の信号処理部２８は、変数解析部２４の解析結果（μnk，αpj，βql，Ｕnm）を適用した信号処理（フィルタ処理）を音響信号Ｓyに対して実行することで音響信号Ｓzを生成する。本実施形態の信号処理部２８は、音響信号Ｓyのうち入力装置１６に対する利用者からの指示に応じた要素成分を抑圧した音響信号Ｓzを生成する。
【００６５】
具体的には、信号処理部２８は、周波数分析部２２が算定した音響信号ＳyのスペクトログラムＹn,fについて以下の数式(31)の演算を実行することで音響信号ＳzのスペクトログラムＺn,fを算定する。数式(31)の演算は、変数解析部２４の解析結果に応じたフィルタＦn,fを音響信号ＳyのスペクトログラムＹn,fに作用させる処理を意味する。
【数３１】
0006044119.tif 000032
信号処理部２８は、数式(31)で算定されたスペクトログラムＺn,fを時間領域の音響信号Ｓzに変換する。例えば、信号処理部２８は、スペクトログラムＺn,fと音響信号Ｓyの位相スペクトログラムとを適用した短時間逆フーリエ変換で音響信号Ｓzを生成する。なお、公知の位相復元法で音響信号Ｓzを生成することも可能である。信号処理部２８が生成した音響信号Ｓzが放音装置１８に供給されて音波として再生される。
【００６６】
数式(31)のフィルタＦn,fは、以下の数式(32)で表現される。
【数３２】
0006044119.tif 000033
数式(32)のフィルタＦn,fの分母は、音響モデルのスペクトログラムＸn,f（数式(6)）に相当する。他方、数式(32)の分子の変数ｕnmは、音響モデルにおけるＭ個（(ＪＫ＋Ｌ)個）の要素成分（調波要素ＥAnj,kおよび非調波要素ＥBl）の音量（以下「調整音量」という）に対応する。Ｍ個の調整音量ｕnmのうち利用者からの指示に応じた要素成分に対応する各調整音量ｕnmは所定値εに設定され、残余の各調整音量ｕnmは変数解析部２４が推定した音量Ｕnmに設定される。所定値εは例えばゼロ（またはゼロに近い正数）に設定される。以上の説明から理解されるように、数式(32)のフィルタＦn,fの分子は、音響モデルのスペクトログラムＸn,fのうち利用者からの指示に応じた特定の要素成分の音量Ｕnmを所定値εに変更したスペクトログラムに相当する。したがって、フィルタＦn,fを音響信号Ｓyに作用させる数式(31)の演算により、音響信号Ｓyから特定の要素成分を抑圧（除去）した音響信号Ｓzが生成される。
【００６７】
利用者は、音響信号Ｓyのうち所望の要素成分を入力装置１６の操作で指定することが可能である。例えばＪ個の調波成分のうち特定の調波成分を利用者が選択した場合、信号処理部２８は、利用者が選択した調波成分のスペクトル包絡ＶAfjとＫ個の調波構造Ｇn,fkの各々との組合せに対応するＫ個の調整音量ｕnmを所定値εに設定し、残余（(Ｍ−Ｋ)個）の各調整音量ｕnmを音量Ｕnmに設定する。したがって、音響信号Ｓyのうち利用者が選択した調波成分（例えば特定の楽器の演奏音）を抑圧した音響信号Ｓzが生成される。
【００６８】
Ｋ個の調波構造Ｇn,fkのうち特定の調波構造Ｇn,fkを利用者が選択した場合、信号処理部２８は、利用者が選択した調波構造Ｇn,fkとＪ個のスペクトル包絡ＶAfjの各々との組合せに対応するＪ個の調整音量ｕnmを所定値εに設定し、残余（(Ｍ−Ｊ)個）の各調整音量ｕnmを音量Ｕnmに設定する。したがって、音響信号Ｓyのうち利用者が選択した調波構造Ｇn,fkに対応する基本周波数μnkの調波成分（すなわち特定の音高）を抑圧した音響信号Ｓzが生成される。
【００６９】
また、Ｌ個の非調波成分のうち特定の非調波成分を利用者が選択した場合、信号処理部２８は、利用者が選択した非調波成分（非調波要素ＥBl）に対応する調整音量ｕnmを所定値εに設定し、残余の各調整音量ｕnmを音量Ｕnmに設定する。したがって、音響信号Ｓyのうち利用者が選択した非調波成分（例えば特定の打楽器の演奏音）を抑圧した音響信号Ｓzが生成される。
【００７０】
図５は、以上に説明した音響解析装置１００による処理結果である。図５では、相異なる２種類の調波性の楽器の演奏音を含む音響信号Ｓy（Ｊ＝２，Ｌ＝０）を楽器毎に分離（一方を抑圧）した場合のＳＮ（Signal/Noise）比が、本実施形態の音響解析装置１００を利用した場合と、非負値行列因子分解（NMF）での分離結果をｋ-means法で楽器毎に分類した場合（以下「対比例」という）とについて対比的に図示されている。ＳＮ比が高いほど分離精度が高いことを意味する。評価用の音楽は、ＲＷＣ（Real World Computing） Music Databeseから選択されたクラシックおよびジャズの音楽である。本実施形態によれば、対比例と比較して音響信号Ｓyの各要素成分を高精度に分離できることが図５から理解される。
【００７１】
＜変形例＞
以上に例示した形態には様々な変形が加えられる。例えば、前述の形態では、Ｊ個の調波成分とＬ個の非調波成分とを含む音響モデルを例示したが、Ｌ個の非調波成分を省略することも可能である。
【００７２】
また、前述の形態では、変数解析部２４の解析結果を表示装置１４による表示と信号処理部２８による信号処理とに適用したが、変数解析部２４の解析結果の利用方法は任意である。例えば、音響信号Ｓyのうち特定の楽器に対応する調波成分の基本周波数μnkの解析結果からその楽器の楽譜を作成する構成（自動採譜）や、音響信号Ｓyの特定の要素成分を解析結果に応じて抽出して選択的に音響効果（例えば残響効果）を付与する構成も採用され得る。
【符号の説明】
【００７３】
１００……音響解析装置、１０……演算処理装置、１２……記憶装置、１４……表示装置、１６……入力装置、１８……放音装置、２２……周波数分析部、２４……変数解析部、２６……表示制御部、２８……信号処理部、５０……解析結果画像。
【図１】
0006044119.tif 000034
【図２】
0006044119.tif 000035
【図３】
0006044119.tif 000036
【図４】
0006044119.tif 000037
【図５】
0006044119.tif 000038