(57)【特許請求の範囲】
【請求項１】
雑音は1つ以上の雑音テンプレートの混合で表現されるものとし、周波数領域に変換された雑音の音信号x(L)j,0を用いて、1つ以上の雑音テンプレートを含み、確率モデルであらわされる雑音モデルを学習する雑音モデル学習部と、
前記音信号x(L)j,0と、周波数領域に変換された既定の動作音kと雑音とが混ざった音信号x(L)i,kとを用いて、前記動作音kにおける周波数ωの重要度γω,kを計算する周波数重要度計算部と、
h=1,2,…,Hkとし、前記重要度γω,kが高いと判断された周波数fk,hと、前記周波数fk,hの総数Hkとを求める正弦波数及び周波数計算部と、
前記音信号x(L)j,0と前記音信号x(L)i,kとを用いて、前記周波数fk,hにおける正弦波強度Ak,hを計算する正弦波強度計算部とを含む、
学習装置。
【請求項２】
請求項１の学習装置で得た雑音モデル及び周波数fk,h、総数Hk、正弦波強度Ak,hを用いて、識別対象の音響信号に含まれる動作音を識別する識別装置であって、
周波数領域に変換された識別対象の音響信号Xω,τと前記雑音モデルに含まれる雑音テンプレートとを用いて、前記雑音テンプレートの混合比を求め、前記雑音テンプレートと前記混合比とから雑音の推定値Nω,τを計算する雑音計算部と、
前記周波数fk,hと前記総数Hkと前記正弦波強度Ak,hとを用いて、正弦波重畳モデルにより、動作音kの近似値を求め、前記音響信号Xω,τと前記近似値と前記推定値Nω,τとから、ある時刻において前記音響信号Xω,τに前記動作音kが含まれる確率であるフレームワイズ動作音確率を計算するフレームワイズ動作音確率計算部とを含む、
識別装置。
【請求項３】
請求項２の識別装置であって、
識別対象の音響信号に含まれると想定される動作音の種類の総数をKとし、k=1,2,…,K、k'=1,2,…,Kとし、ある時刻の動作音k'から次の時刻の動作音kへ移り変わる遷移確率と前記フレームワイズ動作音確率とを用いて、一つ前の時刻の動作音を考慮した、ある時刻において前記音響信号Xω,τに前記動作音kが含まれる確率である動作音確率を計算する動作音確率計算部を含む、
識別装置。
【請求項４】
請求項２または請求項３の識別装置であって、
前記フレームワイズ動作音確率に基づき得られる動作音の推定結果zτと前記近似値と前記推定値Nω,τとを用いて、前記音響信号Xω,τに含まれる動作音を強調するためのフィルタを設計するフィルタ設計部と、
前記音響信号Xω,τと前記フィルタとを用いて、前記音響信号Xω,τに含まれる動作音を強調するフィルタリング部とを含む、
識別装置。
【請求項５】
請求項２から請求項４の何れかの識別装置であって、
前記推定結果zτと、推定結果が得られた時刻とを用いて、動作音のログを生成する、
識別装置。
【請求項６】
雑音は1つ以上の雑音テンプレートの混合で表現されるものとし、周波数領域に変換された雑音の音信号x(L)j,0を用いて、1つ以上の雑音テンプレートを含み、確率モデルであらわされる雑音モデルを学習する雑音モデル学習ステップと、
前記音信号x(L)j,0と、周波数領域に変換された既定の動作音kと雑音とが混ざった音信号x(L)i,kとを用いて、前記動作音kにおける周波数ωの重要度γω,kを計算する周波数重要度計算ステップと、
h=1,2,…,Hkとし、前記重要度γω,kが高いと判断された周波数fk,hと、前記周波数fk,hの総数Hkとを求める正弦波数及び周波数計算ステップと、
前記音信号x(L)j,0と前記音信号x(L)i,kとを用いて、前記周波数fk,hにおける正弦波強度Ak,hを計算する正弦波強度計算ステップとを含む、
学習方法。
【請求項７】
請求項６の学習方法で得た雑音モデル及び周波数fk,h、総数Hk、正弦波強度Ak,hを用いて、識別対象の音響信号に含まれる動作音を識別する識別方法であって、
周波数領域に変換された識別対象の音響信号Xω,τと前記雑音モデルに含まれる雑音テンプレートとを用いて、前記雑音テンプレートの混合比を求め、前記雑音テンプレートと前記混合比とから雑音の推定値Nω,τを計算する雑音計算ステップと、
前記周波数fk,hと前記総数Hkと前記正弦波強度Ak,hとを用いて、正弦波重畳モデルにより、動作音kの近似値を求め、前記音響信号Xω,τと前記近似値と前記推定値Nω,τとから、ある時刻において前記音響信号Xω,τに前記動作音kが含まれる確率であるフレームワイズ動作音確率を計算するフレームワイズ動作音確率計算ステップとを含む、
識別方法。
【請求項８】
請求項１の学習装置、または、請求項２から請求項５の何れかの識別装置として、コンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
音響信号に含まれる所定の音を識別する技術に関する。
【背景技術】
【０００２】
音響信号が何の音であるか（足音、水が流れる音など）を分析する技術として、非特許文献１が知られている。非特許文献１では、音響イベント列の生成過程を、大規模かつ複雑な統計モデルでモデル化することで、音響特徴量の生成モデルや識別器を構築している。生成モデルから、長時間の音響信号毎の状況の生成確率を分析したり、新たに入力された音響イベント列と生成モデルの距離を利用して、音響信号が示す状況を推定することを可能とする。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】井本佳右ほか,“音響イベント列の確率的生成モデルを利用した音響トピックとユーザ行動との関係分析”,日本音響学会，2013年春期音響学会
【発明の概要】
【発明が解決しようとする課題】
【０００４】
音響信号に含まれる所定の音を識別する技術を、工場内の機器の動作音の識別に利用することを考える。工場などに設置された大型の製造機および造型機などの業務用機器は、故障により稼働がストップするだけで、業務に大きな支障をもたらす。そのため、その動作状況を日常的に監視し、故障の発生を事前に防がなくてはならない。解決案として、業務用機器の管理業者が、定期的に現場へ整備員を派遣し、パーツの摩耗などを確認する方法がある。しかし、多大な人件費や移動費、労力がかかるため、すべての業務用機器や工場でこれを実施するのは難しい。代替案として、その機械を構成するパーツごとに、動作状況ログを残す方法がある。例として、図１のような、ベルトコンベア１、塗料噴射器２、アーム３のパーツで構成される機器の動作ログを残すことを考える。「動作ログ１」のように、「いつ」「どのパーツが」動作したかの“詳細な動作ログ”を残せれば、稼働状況や使用頻度から、各パーツの消耗・摩耗度合いなどを推定できる。これにより、故障時期の予測や、パーツ交換の提案などが可能になり、故障率を大幅に落とせる。ところが現在のほとんどの動作ログは、「動作ログ２」のような大まかなログ（作業開始時刻と作業終了時刻だけなど）を残す機能しか備えておらず、詳細な動作ログを残す機能を備えてるものは少ない。
【０００５】
そこで、動作音から詳細な動作ログを作成する方法を考える。例えば、各パーツの動作音に着目し、その動作音から詳細な動作ログを作成する。機器の動作音を収録し、その動作音がどのパーツから発せられていたかを識別することで、これを実現する。動作音から動作ログを作成するメリットは、“詳細な動作ログ”を容易に作成できることに加え、動作ログと動作音を紐づけることで、人間が遠隔地から各パーツごとの動作音を聞き、異常診断ができる点である（図２）。そのための要素技術として、音響信号から動作音の種類を識別し、所望の動作音だけを個別に強調することが必要である。
【０００６】
ところが、非特許文献１では、家庭内などの比較的静かな環境を想定し、MFCCやLPCといった一般的な音響特徴量を用いて識別をしている。大型の業務用機械が動いている環境はえてして、他の機械（例えば、冷却ファン）が大騒音のノイズを出していることが多い。また、該当機器に取り付けられた、一部のパーツが大騒音のノイズを出していることもある（コンプレッサーなど）。そのため、各パーツの動作音は信号雑音比（ＳＮ比）が悪く、MFCCやLPCなどの一般的な特徴量では動作音の識別は難しい。また同様に、音響信号から所望の動作音だけを個別に強調することも困難である。
【０００７】
本発明は、大騒音のノイズが存在する状況においても、音響信号に含まれる動作音を識別することができる識別装置、識別装置において用いられるパラメータを学習する学習装置、その方法及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の課題を解決するために、本発明の一態様によれば、学習装置は、雑音は1つ以上の雑音テンプレートの混合で表現されるものとし、周波数領域に変換された雑音の音信号x(L)j,0を用いて、1つ以上の雑音テンプレートを含み、確率モデルであらわされる雑音モデルを学習する雑音モデル学習部と、音信号x(L)j,0と、周波数領域に変換された既定の動作音kと雑音とが混ざった音信号x(L)i,kとを用いて、動作音kにおける周波数ωの重要度γω,kを計算する周波数重要度計算部と、h=1,2,…,Hkとし、重要度γω,kが高いと判断された周波数fk,hと、周波数fk,hの総数Hkとを求める正弦波数及び周波数計算部と、音信号x(L)j,0と音信号x(L)i,kとを用いて、周波数fk,hにおける正弦波強度Ak,hを計算する正弦波強度計算部とを含む。
【０００９】
上記の課題を解決するために、本発明の他の態様によれば、学習方法は、雑音は1つ以上の雑音テンプレートの混合で表現されるものとし、周波数領域に変換された雑音の音信号x(L)j,0を用いて、1つ以上の雑音テンプレートを含み、確率モデルであらわされる雑音モデルを学習する雑音モデル学習ステップと、音信号x(L)j,0と、周波数領域に変換された既定の動作音kと雑音とが混ざった音信号x(L)i,kとを用いて、動作音kにおける周波数ωの重要度γω,kを計算する周波数重要度計算ステップと、h=1,2,…,Hkとし、重要度γω,kが高いと判断された周波数fk,hと、周波数fk,hの総数Hkとを求める正弦波数及び周波数計算ステップと、音信号x(L)j,0と音信号x(L)i,kとを用いて、周波数fk,hにおける正弦波強度Ak,hを計算する正弦波強度計算ステップとを含む。
【発明の効果】
【００１０】
本発明によれば、大騒音のノイズが存在する状況においても、音響信号に含まれる動作音を識別することができるという効果を奏する。
【図面の簡単な説明】
【００１１】
【図１】動作ログを説明するための図。
【図２】動作ログと動作音を紐づけ、遠隔地で動作音を人間が聞く構成を説明するための図。
【図３】動作音の時間周波数解析結果の例を示す図。
【図４】機器がとりうる状態遷移確率θk,k'の一例を図解した図。
【図５】式(4)の意味を説明するための図。
【図６】学習の手順を示す図。
【図７】識別,強調の手順を示す図。
【図８】第一実施形態に係る学習装置の機能ブロック図。
【図９】動作音について周波数重要度を計算した結果の例を示す図。
【図１０】第一実施形態に係る識別装置の機能ブロック図。
【発明を実施するための形態】
【００１２】
以下、本発明の実施形態について、説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号[-]「~」等は、本来直前の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直後に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【００１３】
＜第一実施形態のポイント＞
観測信号（音響信号）から機械音の動作音の種類を識別する。より詳しくは、識別したい動作音毎に与えられた学習データを用いて、大騒音下で動作音の識別と強調を行う技術に関する。
【００１４】
本実施形態では、大騒音下で特定の動作音の種類を識別し、強調する機械学習ベースの技術を提供する。機械の動作音に特化してこれを実現するために、機械の動作音の以下の特徴に着目する。
【００１５】
＜機械の動作音の4つの特徴＞
1. 各動作音は、毎回同じ音が鳴る。
2. 機械の動作順には法則性がある。
3. 工場内の騒音は、ほぼ定常（冷却ファンなど）か、突発音が繰り返して鳴る（打鋲音など）。
4. 観測音は、動作音と騒音が混ざった音である。
【００１６】
以下では、この特徴について詳しく説明する。
＜1.各動作音は、毎回同じ音が鳴る＞
各パーツの動作音は、モーター等に起因する「ウィーン」のような音と、アラームなどに起因する「ピー」という音の2種類に大別できると考える。製造機器は、製品ごとの個体差をなくすために、毎回同じ動きをする。つまり、各パーツでは、動作音として常に同じような音が鳴ると考えられる。また、「ウィーン」や「ピー」といった動作音は、特徴的な周波数を持つことが多い。図３は、比較的静かな環境で収録した動作音の時間周波数解析結果である。この動作音は、「ピー」という音であり、3000Hzと4000Hzのあたりに強いピークを確認できる。このことから動作音st,kは、以下の正弦波重畳モデルで近似的に表現できるものとする。
【００１７】
【数１】
0006301891.tif 000002
【００１８】
ここでπは円周率、tは時間領域の時間を表すインデックス、k∈{1,2,…,K}は動作音の種類を表すインデックス（ただしKは音響信号に含まれる動作音の種類の総数）、Hkはk番目の動作音を構成する特徴的な周波数の数、h∈{1,2,…,Hk}はk番目の動作音を構成するHk個の特徴的な周波数のうちのh番目の特徴的な周波数を表すインデックス、fk,hとAk,hとはそれぞれk番目の動作音のh番目の特徴的な周波数(Hz)とその強度とを表す。図３を例に挙げれば式（1）は、Hk=2、fk,1=3000、fk,2=4000であり、動作音が2つの正弦波の足し算で表現できることを示している。
【００１９】
＜2.機械の動作には順番に法則性がある＞
機械による製品の製造は、アルゴリズムで記述できる決められた「手順」に基づき行われる。図１を例に挙げれば、「ベルトコンベアによる製品の移動」→「塗料の噴射」→「製品の回収」→「休止」、という手順に従って、機械は動作している。これは現在動作しているパーツは、1時刻前の動作しているパーツに応じて推測できることを示している。時間フレームをτと置き、zτ∈{1,2,…,K}を時刻τでの動作の種類と置くと、zτはzτ-1に依存した1次のマルコフモデルで記述できる。図１を例に挙げると、「ベルトコンベアによる製品の移動」をzτ=1、「塗料の噴射」をzτ=2、「製品の回収」をzτ=3、「休止」をzτ=4と置いたとき、zτ-1=1ならば、zτは1か2しか取りえない。このことを確率を用いて記述すると、以下のように書ける。
【００２０】
【数２】
0006301891.tif 000003
【００２１】
ここでθk,k≧0は、時刻τ-1の状態がk'(ただし、k'∈{1,2,…,K}である)だったとき、時刻τの状態がkである、動作が移り変わる確率（状態遷移確率）を示す。図４は、図１の機器がとりうる状態遷移確率θk,k'の一例を図解したものである。k'=1だったとき、θ3,1θ4,1=0であることが、zτ-1=1ならば、zτは1か2しか取りえないことを示している。
【００２２】
＜3.工場内の騒音は、ほぼ定常か、突発音が繰り返して鳴る＞
工場内の騒音(背景雑音、単に雑音ともいう)は、ほぼ定常（冷却ファンなど）か数個のスペクトル状態で記述できる突発音（打鋲音など）の2種類に大別できると考える。ここで定常な雑音とは、送風機のような「ゴー」という雑音である。また数個のスペクトル状態で記述できる突発音とは、ドリルやハンマーのような「ドドドドドド」など一定のリズムで大きな打撃音が鳴るものである。また、実際の騒音(背景雑音)は、これらの音、つまりM種類の雑音のテンプレートμω,mが、混合比率wm,τ（ただしΣMmwm,τ=1）で混ぜ合わさったものと考え、時刻τでの雑音の音響信号Nω,τを以下のように記述する。
【００２３】
【数３】
0006301891.tif 000004
【００２４】
ここでωは周波数領域の周波数インデックスであり、ω∈{1,2,…,Ω}である。つまり、本実施形態では雑音を1つ以上の雑音テンプレートの混合で表現する。
【００２５】
＜4.観測音は、動作音と騒音が混ざった音である＞
実際にマイクロホンで観測される音（振幅スペクトル）Xω,τは、1.から3.で説明される音が足しあわされて観測されるものと考える。
【００２６】
【数４】
0006301891.tif 000005
【００２７】
ここでSω,z_τ(ただし、下付き添え字のz_τはzτを意味する。)は式(1)のst,kをフーリエ変換し、絶対値を取った、動作音の振幅スペクトルである。式(4)の意味を、図５を使って説明する。まずSω,z_τは、時刻τではzτ番目のパーツが動作しており、その動作音はSω,z_τであることを示している。つまり、式(1)(2)にしたがって動作音が生成されていることを表す。続いてNω,τは時刻τでは、それぞれの雑音源mがwω,τの混合比で混ぜ合わさっている（式(3)）。そして、上述の手順で生成されたSω,z_τとNω,τとが重畳され、観測信号Xω,τが生成される。
【００２８】
本実施形態は、1.から4.の説明に基づき、式(4)の逆問題を解いてzτ（時刻τで動作しているパーツの種類）を識別するものである。また、zτを解くことにより動作音の振幅スペクトルSω,z_τが求まり、動作音だけの強調が可能になる。この逆問題を解くために必要なパラメータを、学習データから学習する。本実施形態の学習と識別,強調の手順の例を、それぞれ図６と図７に示す。以降、まず図６に沿って学習の手順を説明し、次に図７に沿って識別/強調の手順を説明する。
【００２９】
＜学習手順＞
本実施形態の学習手順を実現する学習装置１００について説明する。
【００３０】
学習装置１００は、学習データ(音響信号x(L)t、ラベルファイルlabel)を入力とし、雑音モデルのパラメータΥ={αm,μ1,m,…,μΩ,m,Σm}Mm=1、各動作音kの正弦波重畳モデルパラメータ(K個のHkと、(ΣKk=1Hk)個の周波数fk,hと、(ΣKk=1Hk)個の正弦波強度Ak,h)を学習し、出力する。
【００３１】
学習装置１００は、周波数領域変換部１０２と、データ切り出し部１０３と、雑音モデル学習部１０４と、K個の特徴周波数計算部１０５−ｋとを含む。図８は、学習装置１００の機能ブロック図を示す。
【００３２】
＜学習データ＞
識別対象の機器の近くに（内部でも可能）マイクホロンを設置し、音響信号（以下、マイクロホンで収音された音響信号を「観測信号」ともいう）x(L)tを収集する。なお、上付き添え字（L）は学習データ自体、または、学習データから得られるデータを表すインデックスである。収集した音響信号x(L)tに、生成したい詳細な動作ログと同じ粒度でラベルファイル（動作ログ）labelを人手で作成する。ラベルファイルlabelは、時刻tに対応する音響信号x(L)tがどの動作音に対応するかを示す情報であり、例えば、各時刻tに対応する各音響信号x(L)tに各動作音を示すラベルを付与してもよいし、対象となる学習データ(音響信号x(L)1，x(L)2，…)全てに対して、各動作音を示すラベルとその動作音の開始時刻及び終了時刻を付与してもよい。
【００３３】
このときのサンプリング周波数は、動作音の特徴的な周波数を記録できるよう設定すればよい。例えば16kHzなどに設定する。
【００３４】
さらに可能であれば、識別対象の機器の近くに複数のマイクロホンを設置したり、雑音発生源の近くにマイクロホンを設置してもよい。複数のマイクロホンを設置した場合、ビームフォーミングやスペクトルサブトラクションなどの技術が使用可能になるため、動作音、雑音の推定精度が向上する。
【００３５】
＜周波数領域変換部１０２＞
周波数領域変換部１０２は、音響信号x(L)tを受け取り、音響信号x(L)tを短時間フーリエ変換（STFT）などの手法で周波数領域に変換する。その後、複素スペクトルの絶対値を取り、振幅スペクトル(以下、この振幅スペクトルを単に「周波数領域に変換された音響信号」または「音響信号」ともいう)X(L)τ=(X(L)1,τ，X(L)2,τ，…,X(L)ω,τ，…,X(L)Ω,τ，)Tを得（Ｓ１０２）、出力する。ただし、Tは転置を表す。ここで、フーリエ変換長は32ms、シフト幅は16msなどに設定できる。
【００３６】
＜データ切り出し部１０３＞
データ切り出し部１０３は、周波数領域に変換された音響信号X(L)τと、ラベルファイルlabelとを受け取り、ラベルファイルlabelを元に、周波数領域に変換された音響信号X(L)τを、雑音と、動作音の種類ごとに分割し（Ｓ１０３）、雑音の音響信号x(L)j,0と動作音kの音響信号x(L)i,kと出力する。ここで、雑音の音響信号x(L)j,0には雑音のみが、動作音kの音響信号x(L)i,kには動作音と雑音が重畳された状態となっている。なお、音響信号x(L)j,0に含まれる雑音と音響信号x(L)i,kとに含まれる雑音は、同じ種類の雑音である。ただし、そのレベルは異なってもよい。
【００３７】
なお、音響信号X(L)τから切り出されたk番目の動作音のIkフレーム分の音響信号を
x(L)i,k=(X(L)1,i,k，X(L)2,i,k，…,X(L)ω,i,k,…,X(L)Ω,i,k)T 
とおき、音響信号X(L)τから切り出された雑音のJフレーム分の音響信号を
x(L)j,0=(X(L)1,j,0，X(L)2,j,0，…,X(L)ω,j,0,…,X(L)Ω,j,0)T 
とおく。ただし、Ikは動作音k毎に異なり、i=1,2,…,Ikである。
【００３８】
なお、ある時刻τの音響信号X(L)τに2つ以上の動作音が含まれる(2つ以上の動作音が重畳された状態の)場合には、推定精度を高めるために、その部分を排除し、以下の学習には用いない構成としてもよい。
【００３９】
＜雑音モデル学習部１０４＞
雑音モデル学習部１０４は、音響信号X(L)τから切り出された雑音のJフレーム分の音響信号x(L)j,0を受け取り、この値を用いて、雑音モデルを学習し（Ｓ１０４）、学習装置１００の出力値として雑音モデルパラメータΥを出力する。ここで、雑音モデルは、M個(Mは1以上の整数)の雑音テンプレートμ1,m,μ2,m,…,μΩ,m(ただし、m=1,2,…,M)を含み、確率モデルであらわされる。言い換えると、学習時には、雑音のみの学習データ(雑音のJフレーム分の音響信号x(L)j,0)を用いて、何らかの確率モデルにより雑音モデルを事前学習する。GMMの場合は、何らかの確率モデルとして「混合ガウス分布（GMM）」を用いる。GMMの場合、EMアルゴリズムなどの手法を用いて雑音モデルを学習する(例えば、参考文献１参照)。
[非特許文献２]小西貞則、“多変量解析入門、付録C EMアルゴリズム”、岩波書店、pp.294-298、2010年
この場合、雑音データ（ここでは、雑音のJフレーム分の音響信号x(L)j,0=(X(L)1,j,0，X(L)2,j,0，…,X(L)Ω,j,0)Tであり、X(L)ω,j,0をN(L)ω,τとも記載する(X(L)ω,j,0=N(L)ω,τ)。ここでｊは、時刻フレームτの何れかに対応する）は以下の確率分布に従って生成されたものとしてモデル化する。
【００４０】
【数５】
0006301891.tif 000006
【００４１】
ここでαm、μω,m、Σmはそれぞれ、m番目の混合分布の混合比、周波数ωにおける平均ベクトル(雑音テンプレート)、共分散行列である。この場合、雑音モデルパラメータは、雑音GMMパラメータであり、Υ={αm,μ1,m,…,μΩ,m,Σm}Mm=1で表される。なお、共分散行列Σmは対角行列に制限してもよいし、m番目の混合分布の分散σ2mを用いて、σ2mIと制限してもよい。ここでIは単位行列を表す。
【００４２】
GMMの学習は、雑音の音響信号N(L)ω,τを構成する雑音をM種類にクラスタリングしていることに相当する。つまり、工場の中で鳴っているM種類の雑音を推定しており、μω,mはm番目の雑音テンプレートとみなすこともできる。そのため、GMMの混合数Mは、工場内で想定される雑音源の数を設定すべきであり、例えば8などに設定できる。
【００４３】
なお、学習データを得る際に、複数のマイクロホンを配置した場合、GMMの学習と平均ベクトル(雑音テンプレート)μω,mの作成手順が異なる。複数チャネルを用いてビームフォーミングやスペクトラルサブトラクションを施した音響特徴量ベクトルをN(add)τとする。この音響特徴量ベクトルN(add)τを、音響信号X(L)ω,τから切り出された雑音のJフレーム分の音響信号X(L)ω,j,0=N(L)ω,τと結合させて、以下のGMMを学習する。
【００４４】
【数６】
0006301891.tif 000007
【００４５】
つまり雑音GMMパラメータは、Υ={αm,μ-m,Σm}Mm=1であり、雑音テンプレートはμ-mの上位Ω次元となる。
【００４６】
なお、識別・強調時に雑音を推定するには、事前学習した雑音モデルを用いて、観測信号から雑音を推定する。推定法には、最小平均二乗誤差(MMSE)推定や最大事後確率(MAP)推定などを用いることができる。どの推定法を用いるかは任意であるが、雑音モデルに用いた確率分布の種類に応じて、推定精度と計算速度のトレードオフを考慮しながら決めることが望ましい。なお、GMMの場合は、このトレードオフを考慮して、MMSE推定を使うことが一般的である。
【００４７】
＜特徴周波数計算部１０５−ｋ＞
Ｋ個の特徴周波数計算部１０５−ｋはそれぞれ以下の処理を行う。特徴周波数計算部１０５−ｋは、雑音の音響信号x(L)j,0と動作音kの音響信号x(L)i,kとを受け取り、動作音kの正弦波重畳モデルパラメータ(ここでは、正弦波数Hk,その周波数fk,h,正弦波強度Ak,h)を学習し（Ｓ１０５）、出力する。例えば、特徴周波数計算部１０５−ｋは、周波数重要度計算部１０５Ａと、正弦波数及び周波数計算部１０５Ｂと正弦波強度計算部１０５Ｃとを含む（図８参照）。
【００４８】
(周波数重要度計算部１０５Ａ)
周波数重要度計算部１０５Ａは、雑音の音響信号x(L)j,0と動作音kの音響信号x(L)i,kとを受け取り、これらの値を用いて、動作音kにおける周波数ωの重要度γω,kを計算し（Ｓ１０５Ａ）、出力する。なお、重要度γω,kとは、周波数ωの成分が前記動作音kを識別する際に重要な成分か否かを示す指標である。例えば、以下のようにして重要度γω,kを計算する。
【００４９】
まず、雑音の音響信号x(L)j,0と動作音kの音響信号x(L)i,kとを用いて、ロジスティック回帰モデルを学習する（例えば参考文献２参照）。
(参考文献２)小西貞則、“多変量解析入門、4章ロジスティック回帰モデル”、岩波書店、pp.79-96、2010年。
【００５０】
説明変数行列Xkと応答変数ベクトルYkをそれぞれ
【００５１】
【数７】
0006301891.tif 000008
と置き、以下のロジスティック回帰モデルを学習する。
【００５２】
【数８】
0006301891.tif 000009
【００５３】
ここでβk=(β1,k,β2,k,…,βΩ,k)Tは回帰ベクトルである。
【００５４】
ロジスティック回帰モデルの回帰ベクトルの特徴として、識別に有効な次元に対応する回帰係数の絶対値が大きくなる。すなわち、|βω,k|が大きな周波数ほど、動作音kの識別に有効な周波数、すなわち重要な周波数とみなすことができる。そこで本実施形態では周波数ωに対する重要度γω,kを以下のように定義する。
【００５５】
【数９】
0006301891.tif 000010
【００５６】
この値が大きな周波数ほど、その動作音の特徴をよく表現する周波数と考える。2つの異なる動作音について周波数重要度を計算した結果の例を図９に示す。いくつかの周波数において、周波数重要度のピークが確認できる。なお、回帰ベクトルの過学習を避け、安定して周波数重要度を計算するために、ロジスティック回帰の学習時に正則化項を付与した方がよい。この正則化項には、回帰ベクトルのL1ノルムを用いるLASSOや回帰ベクトルのL2ノルムを用いるリッヂ回帰を使うことができる。
【００５７】
(正弦波数及び周波数計算部１０５Ｂ)
正弦波数及び周波数計算部１０５Ｂは、重要度γω,kを受け取り、重要度γω,kが高いと判断された周波数fk,hと、周波数fk,hの総数Hkとを求め（Ｓ１０５Ｂ）、出力する。なお、重要度γω,kが高いとは、周波数ωの成分が動作音kを識別する際に重要な成分であることを意味する。例えば、周波数重要度に対する閾値γTHを設定する。重要度γω,kと閾値γTHとの大小関係により、重要度γω,kが高いか否かを判断する。重要度γω,kが式(12)で表される場合には、重要度γω,kが大きければ大きいほど、動作音kに対して周波数ωが重要であると考えられるので、γω,k>γTHとなる周波数が動作音kを識別する際に重要な成分であると考える。そこで、γω,k>γTHとなる周波数をfk,hとする。ただし、fk,hは周波数の値自体であってもよいし、それと等価な値（例えば周波数のインデックス）であってもよい。さらに、γω,k>γTHとなる周波数fk,hの個数を数え、その総数をHkとする。総数Hkと周波数fk,hの例を図９に示す。動作音１を例に挙げれば、Hk=2であり、fk,1は2100Hz付近、fk,2は4200Hz付近となる。なお閾値γTHは総数Hkが大きくなりすぎないよう経験的に決定すべきであり、たとえば0.5に設定できる。
【００５８】
(正弦波強度計算部１０５Ｃ)
正弦波強度計算部１０５Ｃは、総数Hkと、Hk個の周波数fk,hと、雑音の音響信号x(L)j,0と動作音kの音響信号x(L)i,kとを受け取り、周波数fk,hにおける音響信号の正弦波強度Ak,hを計算し（Ｓ１０５Ｃ）、出力する。例えば、正弦波強度Ak,hは以下のように計算する。
【００５９】
【数１０】
0006301891.tif 000011
【００６０】
X(L)f_k,h,i,k(ただし、下付添え字f_k,hはfk,hを意味する。)はフレームiの音響信号x(L)i,k=(X(L)1,i,k，X(L)2,i,k，…,X(L)Ω,i,k)Tのうちの周波数fk,hに対応する周波数成分を意味し、X(L)f_k,h,j,0は、フレームjのx(L)j,0=(X(L)1,j,0，X(L)2,j,0，…,X(L)Ω,j,0)Tのうちの周波数fk,hに対応する周波数成分を意味する。また万が一、Ak,hが0を下回った場合、Ak,h=0とする。ここで、式(13)の第一項は動作音kがあるときの周波数fk,hの周波数成分(振幅スペクトル値)の平均値、第二項は雑音の周波数fk,hの周波数成分(振幅スペクトル値)の平均値である。
【００６１】
学習装置１００は、以上の処理により、学習データから雑音モデルパラメータΥ={αm,μ1,m,…,μΩ,m,Σm}Mm=1と正弦波重畳モデルパラメータ(K個のHkと、(ΣKk=1Hk)個の周波数fk,hと、(ΣKk=1Hk)個の正弦波強度Ak,h)とを求め、出力する。
【００６２】
＜識別及び強調手順＞
図７を用いて、本実施形態の識別及び強調手順を実現する識別装置２００について説明する。
【００６３】
識別装置２００は、雑音モデルパラメータΥ={αm,μ1,m,…,μΩ,m,Σm}Mm=1と正弦波重畳モデルパラメータ(K個のHkと、(ΣKk=1Hk)個の周波数fk,hと、(ΣKk=1Hk)個の正弦波強度Ak,h)とを受け取り、識別及び強調処理を行う前に予め、設定しておく。
【００６４】
識別装置２００は、識別対象の音響信号xtを受け取り、音響信号xtに含まれる動作音を識別し（例えば、時刻tの音響信号xtにどの動作音kが含まれるかを識別する）、動作音の種類zτと、動作音を強調した音響信号x(E)tを出力する。
【００６５】
なお、識別対象の音響信号xtは、学習データと同じ状況で収集される。
【００６６】
図８は、識別装置２００の機能ブロック図の例を示す。識別装置２００は、周波数領域変換部２０２、雑音計算部２０３、フレームワイズ動作音確率計算部２０４、動作音確率計算部２０５、動作音推定部２０６、ウィナーフィルタ設計部２０７、ウィナーフィルタリング部２０８及び時間領域変換部２０９を含む。
【００６７】
＜周波数領域変換部２０２＞
周波数領域変換部２０２は、識別対象の音響信号xtを受け取り、音響信号xtを短時間フーリエ変換（STFT）などの手法で周波数領域に変換する。その後、複素スペクトルの絶対値を取り、振幅スペクトル(以下、この振幅スペクトルを単に「周波数領域に変換された音響信号」または「音響信号」ともいう)Xτ=(X1,τ,X2,τ,…,Xω,τ,…,XΩ,τ)Tを得（Ｓ２０２）、出力する。ここで、変換パラメータは、学習時と同じものを用いる。
【００６８】
＜雑音計算部２０３＞
雑音計算部２０３は、音響信号Xτと雑音モデルパラメータΥとを受け取り、雑音の推定値Nτ=(N1,τ,N2,τ,…,Nω,τ,…,NΩ,τ)Tを式(3)で計算し（Ｓ２０３）、出力する。例えば、音響信号Xτと雑音モデルパラメータΥに含まれる雑音テンプレートμ1,m,…,μΩ,mとを用いて、雑音テンプレートμ1,m,…,μΩ,mの混合比wm,τを求め、雑音テンプレートμ1,m,…,μΩ,mと混合比wm,τとから雑音の推定値Nτを計算する。雑音GMMを用いた場合、まず、混合比wm,τを以下の式で計算する。
【００６９】
【数１１】
0006301891.tif 000012
【００７０】
その後、推定値Nτ=(N1,τ,N2,τ,…,Nω,τ,…,NΩ,τ)Tを式(3)で計算し、出力する。
【００７１】
【数１２】
0006301891.tif 000013
【００７２】
なお、複数のマイクロホンを用いて音響信号を収録している際は、複数チャネルを用いてビームフォーミングやスペクトラルサブトラクションを施した音響特徴量ベクトルX(add)τを、観測データと結合させて、式(14)を計算すればよい。例えば、式(14)〜(16)において、
【００７３】
【数１３】
0006301891.tif 000014
とする。
【００７４】
＜フレームワイズ動作音確率計算部２０４＞
フレームワイズ動作音確率計算部２０４は、音響信号Xτと、雑音の推定値Nτと、正弦波重畳モデルパラメータ(ここでは、K個のHkと、(ΣKk=1Hk)個の周波数fk,hと、(ΣKk=1Hk)個の正弦波強度Ak,h)とを受け取り、まず、正弦波重畳モデルにより、動作音kの近似値を求める。ここでは、正弦波重畳モデルパラメータを用いて、式(1)により時間領域での動作音kの近似値st,kを計算する。
【００７５】
【数１４】
0006301891.tif 000015
【００７６】
計算した近似値st,kを窓関数で切り出し、短時間フーリエ変換（STFT）などの手法で周波数領域に変換する。その後、複素スペクトルの絶対値を取り、振幅スペクトル(以下、この振幅スペクトルを単に「周波数領域に変換された動作音の近似値」または「動作音の近似値」ともいう)Sτ=(S1,τ,S2,τ,…,Sω,τ,…,SΩ,τ)Tを得る。ここで、変換パラメータは、学習時と同じものを用いる。なお、動作音の近似値Sτは時刻τごとに計算するのではなく、正弦波重畳モデルパラメータを取得した段階で、識別及び強調処理を行う前に計算しておく。このため、上述の通り、サンプリング周波数は、動作音の特徴的な周波数を記録できるように設定すればよい。動作音の近似値Sτは、フレームワイズ動作音確率計算部２０４において以降の処理で用いられるとともに、ウィナーフィルタ設計部２０７に出力される。
【００７７】
次に、雑音の推定値Nω,τと、動作音kの近似値Sω,kとを重畳し、以下の変数を計算する。
【００７８】
【数１５】
0006301891.tif 000016
【００７９】
この変数X~ω,τ,kは、時刻τで動作音kが鳴っていた場合の、音響信号Xω,τの推定値である。
【００８０】
そして、音響信号Xω,τと、推定値X~ω,τ,kとを用いて、時刻τ-1の状態と関係なく時刻τだけの情報から、ある時刻τにおいて音響信号Xω,τに動作音kが含まれる確率であるフレームワイズ動作音確率p(zτ|Xω,τ,Nω,τ)を以下のように計算し(Ｓ２０４)、出力する。
【００８１】
【数１６】
0006301891.tif 000017
【００８２】
＜動作音確率計算部２０５＞
動作音確率計算部２０５は、フレームワイズ動作音確率p(zτ|Xω,τ,Nω,τ)と遷移確率θk,k'とを受け取り、これらの値を用いて、時刻τ-1の動作音kを考慮した、ある時刻τにおいて音響信号Xω,τに動作音kが含まれる確率である動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1)を計算し（Ｓ２０５）、出力する。動作の状態を逐次的に推定するために、隠れマルコフモデルの前向き状態推定アルゴリズムに基づき、以下のように動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1)を求める。
【００８３】
【数１７】
0006301891.tif 000018
【００８４】
ここで遷移確率θk,k'は、学習時に最尤推定で求めてもよいし、機械の動作法則に従って人手で定めてもよい。なお、実用上は、動作音の学習データがほとんど集まらないことが多いので、遷移確率θk,k'は人手で与えることが多い。より詳しくは、遷移確率θk,k'を学習するのに十分なラベル量が、雑音GMMパラメータや正弦波重畳モデルパラメータを学習するのに必要なラベル量よりも圧倒的に多いため、費用対効果や人件費を考慮すると、人手で遷移確率θk,k'を与えた方（チューニングした方が）が効率的な場合が多い。
【００８５】
＜動作音推定部２０６＞
動作音推定部２０６は、動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1)を受け取り、時刻τでの動作音の種類zτを、以下の式に従って推定し（Ｓ２０６）、推定結果である種類zτを出力する。
【００８６】
【数１８】
0006301891.tif 000019
【００８７】
式(22)による状態推定は、最大事後確率（MAP）推定と呼ばれ、時刻τまでのデータを得た下で、最も事後確率(動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1))の高い状態kを時刻τでの動作音の種類zτと推定することを示している。推定結果zτは、フレームワイズ動作音確率p(zτ|Xω,τ,Nω,τ)に基づき得られる動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1)を利用するため、推定結果zτは、間接的に、フレームワイズ動作音確率p(zτ|Xω,τ,Nω,τ)に基づき得られる値と言える。
【００８８】
なお、図２のように動作ログと動作音を紐づける必要がなければ、以降の処理を省略し、推定結果である種類zτのみを識別装置２００の出力値として出力すればよい。
【００８９】
＜ウィナーフィルタ設計部２０７＞
ウィナーフィルタ設計部２０７は、推定結果である種類zτと雑音の推定値Nτ(雑音計算部２０３の出力値)と動作音の近似値Sk（フレームワイズ動作音確率計算部２０４においてフレームワイズ動作音確率の計算過程で得られる値）とを受け取り、次式により、音響信号Xω,τに含まれる動作音k=zτだけを強調するためのウィナーフィルタGτ=(G1,τ,G2,τ,…,Gω,τ,…,GΩ,τ)Tを設計し（Ｓ２０７）、出力する。
【００９０】
【数１９】
0006301891.tif 000020
【００９１】
＜ウィナーフィルタリング部２０８＞
ウィナーフィルタリング部２０８は、音響信号XτとウィナーフィルタGτとを受け取り、音響信号XτをウィナーフィルタGτでフィルタリングし（Ｓ２０８）、音響信号Xτに含まれる動作音k=zτを強調し、強調した音響信号X(E)τ=(X(E)1,τ,X(E)2,τ,…,X(E)ω,τ,…,X(E)Ω,τ)Tを出力する。例えば、次式のように、音響信号XτにウィナーフィルタGτを乗算することでフィルタリングする。
X(E)ω,τ=Gω,τXω,τ
【００９２】
＜時間領域変換部２０９＞
時間領域変換部２０９は、音響信号X(E)τを受け取り、周波数領域変換部２０２で行った周波数領域への変換方法に対応する時間領域への変換方法(例えばIFFT)を用いて、音響信号X(E)τを時間領域の音響信号x(E)tに変換し（Ｓ２０９）、識別装置２００の出力値として出力する。
【００９３】
上述の方法によって、識別装置２００は、動作音の種類zτと強調した音響信号x(E)tを出力する。なお、強調した音響信号x(E)tを必要としない場合には、ウィナーフィルタ設計部２０７、ウィナーフィルタリング部２０８、時間領域変換部２０９を設けず、その処理（Ｓ２０７〜Ｓ２０９）を省略してもよい。
【００９４】
＜効果＞
以上の構成により、大騒音のノイズが存在する状況においても、音響信号に含まれる動作音を識別することができる。
【００９５】
本実施形態では、機械音に特化して大騒音下で特定の動作音の種類を識別し強調するために、上述の機械の動作音の4つの特徴に着目した、識別/強調アルゴリズムとその学習法を提供する。なお、識別対象の動作法則に着目した、隠れマルコフモデルの前向き状態推定アルゴリズム自体は、古くからある技術である(非特許文献１等)。しかし、1.動作音を正弦波重畳モデルで近似計算（式(1)）し、2.雑音を、雑音テンプレートの混合で表現することで（式(3)(14)）、3.音響信号からの動作音の種類を、少量の学習データから効率よく、学習・識別/強調するアルゴリズム（システム全体）については従来技術にはなく、これにより、大騒音のノイズが存在する状況においても、音響信号に含まれる動作音を識別することができる。
【００９６】
また、正弦波重畳モデルのパラメータ推定で用いたロジスティック回帰のパラメータ推定は公知の技術であるが、本実施形態では、この公知技術を利用して、回帰係数の絶対値比率（式(12)）を周波数重要度とみなし、正弦波重畳モデルのパラメータ推定に用いている(Ｓ１０５Ａ〜Ｓ１０５Ｃ)。
【００９７】
また、統計的パターン認識に基づく状態推定は、大量の学習データ必要である。例えば、音声認識などは、人間の発話だけで100時間以上の学習データが必要であった。本実施形態では、機械の動作音の動作特性に着目し、学習法を工夫することで、少ない学習データ量でパラメータ推定が可能である。具体的には、各パーツの動作が50回程度（機械の種類にもよるが、15分程度）のデータで学習が可能である。また、動作音を識別するだけでなく、強調する方法も同時に提供することで、遠隔地から動作音に基づく異常診断を行えるようになり、移動時間・移動経費の削減など、業務効率化にも寄与する。
【００９８】
＜変形例＞
また、推定結果zτと、推定結果が得られた時刻τとを用いて、動作音のログ(図１参照)を生成してもよい。このような構成により詳細な動作ログを容易に残すことができる。工場などに設置された大型の製造機・造型機などの、詳細な動作ログを音響信号から作成することで、機械の稼働状況や使用頻度から、各パーツの消耗・摩耗度合いなどを推定でき、故障時期の予測や交換時期の提案などが可能になり、故障率を大幅に落とせる。産業、特に製造業の効率化に寄与する技術である。
【００９９】
本実施形態では、動作音推定部２０６において、動作音確率p(zτ|Xω,τ,Nω,τ,zτ-1)を用いて、時刻τでの動作音の種類zτを推定しているが、フレームワイズ動作音確率p(zτ|Xω,τ,Nω,τ)を用いて、次式により、時刻τでの動作音の種類zτを推定してもよい。
【０１００】
【数２０】
0006301891.tif 000021
【０１０１】
ただし、この場合、機械の動作音の4つの特徴のうちの1つである、「2. 機械の動作順には法則性がある。」を利用していないため、推定精度は落ちる。
【０１０２】
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０１０３】
＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
【０１０４】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【０１０５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。
【０１０６】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１０７】
また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
0006301891.tif 000022
【図２】
0006301891.tif 000023
【図３】
0006301891.tif 000024
【図４】
0006301891.tif 000025
【図５】
0006301891.tif 000026
【図６】
0006301891.tif 000027
【図７】
0006301891.tif 000028
【図８】
0006301891.tif 000029
【図９】
0006301891.tif 000030
【図１０】
0006301891.tif 000031