JP6535112B2 - マスク推定装置、マスク推定方法及びマスク推定プログラム - Google Patents

マスク推定装置、マスク推定方法及びマスク推定プログラム Download PDF

Info

Publication number
JP6535112B2
JP6535112B2 JP2017567967A JP2017567967A JP6535112B2 JP 6535112 B2 JP6535112 B2 JP 6535112B2 JP 2017567967 A JP2017567967 A JP 2017567967A JP 2017567967 A JP2017567967 A JP 2017567967A JP 6535112 B2 JP6535112 B2 JP 6535112B2
Authority
JP
Japan
Prior art keywords
mask
distribution
time frequency
shape
feature quantity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2017567967A
Other languages
English (en)
Other versions
JPWO2017141542A1 (ja
Inventor
信貴 伊藤
信貴 伊藤
荒木　章子
章子 荒木
中谷　智広
智広 中谷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nippon Telegraph and Telephone Corp
Original Assignee
Nippon Telegraph and Telephone Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Telegraph and Telephone Corp filed Critical Nippon Telegraph and Telephone Corp
Publication of JPWO2017141542A1 publication Critical patent/JPWO2017141542A1/ja
Application granted granted Critical
Publication of JP6535112B2 publication Critical patent/JP6535112B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0272—Voice signal separating
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2134—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0272—Voice signal separating
    • G10L21/028—Voice signal separating using properties of sound source
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0272—Voice signal separating
    • G10L21/0308—Voice signal separating characterised by the type of parameter measurement, e.g. correlation techniques, zero crossing techniques or predictive techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computational Linguistics (AREA)
  • Quality & Reliability (AREA)
  • Signal Processing (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Circuit For Audible Band Transducer (AREA)
  • Complex Calculations (AREA)

Description

本発明は、マスク推定装置、マスク推定方法及びマスク推定プログラムに関する。
複数のマイクロホンにより取得された観測信号から目的音を推定する際に用いられる、各時間周波数点における各音響信号の寄与度を表すマスクを推定する技術がある。マスクは、観測信号の雑音除去や音源分離等に用いられる。
例えば、各音響信号の到来方向が互いに異なることに着目し、各時間周波数点において各音響信号音がどの方向から到来するかの情報を表す特徴量ベクトルをクラスタリングすることにより、マスクを推定する技術がある。
M. Souden, S. Araki, K. Kinoshita, T. Nakatani, and H. Sawada, "A Multichannel MMSE-Based Framework for Speech Source Separation and Noise Reduction," IEEE TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 21, NO. 9, SEPTEMBER 2013, pp. 1913 - 1928.
しかしながら、上記技術では、特徴量ベクトルをモデル化する混合分布の位置パラメータ及び広がりパラメータを更新することにより、混合分布の要素分布の位置及び広がりを調整する。そのため、要素分布は、ある軸に関して回転対称な分布形状を表すに過ぎない。特徴量ベクトルの分布形状は、マイクロホン位置や音響空間の伝達特性等の様々な影響を受け、必ずしも回転対称になるとは限らず、例えば楕円状の分布形状となる場合がある。そのため、上記技術では、混合分布を特徴量ベクトルの分布に十分近づけることができず、マスク推定の精度が必ずしも高くないという問題がある。
本願が開示する実施形態の一例は、上記に鑑みてなされたものであって、例えば、マスク推定の精度をより向上させるマスク推定装置、マスク推定方法及びマスク推定プログラムを提供することを目的とする。
本願の実施形態の一例において、例えば、マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する。マスク推定装置は、特徴量ベクトルと、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を示すマスクを更新する。マスク推定装置は、更新されたマスクから、混合重みを更新する。マスク推定装置は、特徴量ベクトルと、マスクとから、形状パラメータを更新する。
本願が開示する実施形態の一例によれば、例えば、特徴量ベクトルの分布が回転対称でない場合でも高精度なマスク推定が可能になる。
図１は、実施形態１（実施形態２）に係るマスク推定装置の一例を示す図である。 図２は、実施形態１（実施形態２）に係るマスク推定処理の一例を示すフローチャートである。 図３は、実施形態３に係る目的音推定システムの一例を示す図である。 図４は、実施形態３に係る目的音推定処理の一例を示すフローチャートである。 図５は、プログラムが実行されることにより、実施形態に係るマスク推定装置及び目的音推定装置を含む目的音推定システムが実現されるコンピュータの一例を示す図である。 図６は、背景技術に係るマスク推定装置の一例を示す図である。
以下、本願が開示する技術に関する実施形態の一例について、図面を参照して説明する。なお、以下の実施形態により、本願の開示技術が限定されるものではない。また、以下の実施形態は、矛盾しない範囲で適宜組合せてもよい。
なお、以下の記載において、“Ａ”に対して“＾Ａ”と表記する場合は、「Ａの直上に＾が付された記号」と同等であるとする。また、“Ａ”がベクトルである場合には「ベクトルＡ」、“Ａ”が行列である場合には「行列Ａ」、“Ａ”が集合である場合には「集合Ａ」と表記する。
［背景技術について］
本願の実施形態の開示に先立ち、本願が開示する技術の基礎となる背景技術について説明する。背景技術は、文献１「M. Souden, S. Araki, K. Kinoshita, T. Nakatani, and H. Sawada, “A Multichannel MMSE-Based Framework for Speech Source Separation and Noise Reduction,” IEEE TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 21, NO. 9, SEPTEMBER 2013, pp. 1913 - 1928.」に開示されている。
図６は、背景技術に係るマスク推定装置の一例を示す図である。背景技術に係るマスク推定装置１０Ａは、記憶部２０Ａが接続される。マスク推定装置１０Ａは、特徴抽出部１１Ａ、マスク更新部１２Ａ、混合重み更新部１３Ａ、位置パラメータ更新部１４Ａ、広がりパラメータ更新部１５Ａを有する。
記憶部２０Ａは、後述するように特徴抽出部１１Ａにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの、各要素分布の混合重み、各要素分布の位置を表す位置パラメータと、各要素分布の広がりを表す広がりパラメータを記憶する。
特徴抽出部１１Ａは、Ｎ個（Ｎ＞１）の音響信号が混在する状況において、異なる位置で収録されたＭ個（Ｍ＞１）の観測信号を受け取り、時間周波数点毎にすべての観測信号の時間周波数成分をまとめたＭ次元縦ベクトルを構成し、Ｍ次元縦ベクトルを用いて所定のノルムを持つＭ次元の特徴量ベクトルを抽出する。
マスク更新部１２Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルを特徴抽出部１１Ａから受け取り、記憶部２０Ａに保存されている混合重みと、位置パラメータと、広がりパラメータを記憶部２０Ａから受け取り、各時間周波数点に各要素分布がどの程度寄与しているかの割合を表すマスクを更新する。
混合重み更新部１３Ａは、マスク更新部１２Ａにより更新されたマスクを受け取り、混合重みを更新する。位置パラメータ更新部１４Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルと、マスク更新部１２Ａにより更新されたマスクを受け取り、位置パラメータを更新する。広がりパラメータ更新部１５Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルと、マスク更新部１２Ａにより更新されたマスクを受け取り、広がりパラメータを更新する。
上記文献１では、音響信号の数Ｎを既知と仮定すると共に、Ｎ−１個が目的音であり、１個が背景雑音であると仮定する。Ｎ＝２の場合は、例えば、背景雑音の存在下で一人の音声をＭ個のマイクロホンで収録する場合に対応する。また、Ｎ＞２の場合は、例えば、背景雑音の存在下でＮ−１人による会話をＭ個のマイクロホンで収録する場合に対応する。
いま、Ｍ個の観測信号のそれぞれに短時間フーリエ変換等の時間周波数分析を適用して得られる各観測信号の時間周波数成分を、時間周波数点毎にまとめたＭ次元縦ベクトル（以下、観測信号ベクトルと呼ぶ）をｙ(t,f)で表す。ｔ（ｔ＝１，・・・Ｔ。ｔは正整数）は時刻、ｆ（ｆ＝１，・・・Ｆ。ｆは正整数）は周波数ビン番号である。
上記文献１では、目的音はスパース性を有し、各時間周波数点において高々一つの目的音のみが存在すると仮定され、また、背景雑音はすべての時間周波数点において存在すると仮定されている。このとき、観測信号ベクトルｙ(t,f)は、各時間周波数点において次のＮ個の状態のうちのいずれかを取る。
状態Ｓn：目的音のうちｎ番目の目的音のみが存在する状態（ｎ＝１，・・・，Ｎ−１。ｎは正整数）
状態ＳN：目的音が一つも存在しない状態
このため、観測信号ベクトルｙ(t,f)は、下記（１−１）式又は（１−２）式のいずれかでモデル化できるとされている。
Figure 0006535112
上記（１−１）式は、当該時間周波数点においてｎ番目の目的音のみが存在する場合を表し、上記（１−２）式は、当該時間周波数点において目的音が一つも存在しない場合を表す。上記（１−１）式及び（１−２）式におけるベクトルｓn(t,f)はｎ番目の目的音に対応する成分を表し、ベクトルｖ(t,f)は背景雑音に対応する成分を表す。
マスク推定装置１０Ａの目的は、観測信号ベクトルｙ(t,f)が与えられたときに、各時間周波数点(t,f)に各状態Ｓn（ｎ＝１，・・・，Ｎ）がどの程度寄与しているかの割合を表すマスクγn(t,f)を推定することである。マスクγn(t,f)は、γ1(t,f)＋・・・＋γN(t,f)＝１，０≦γn(t,f)≦１を満たす。マスクγn(t,f)は、例えば、当該時間周波数点での状態がＳnの場合にγn(t,f)＝１、そうでない場合にγn(t,f)＝０と定義することができる。あるいは、より一般的には、マスクγn(t,f)を当該時間周波数点における状態Ｓnの事後確率と定義することもできる。以下の説明では、マスクγn(t,f)は、当該時間周波数点における状態Ｓnの事後確率と定義する。
マスクγn(t,f)が推定できると、各目的音ｓn(t,f)の推定等、さまざまな応用に用いることができる。例えば、マスクγn(t,f)を用いて当該目的音が鳴っている時間周波数成分を集めてくることにより、当該目的音の性質を表す統計量を推定することができる。この統計量を用いて設計したフィルタにより各目的音ｓn(t,f)を推定することができる。
上記文献１では、各音響信号の到来方向は互いに異なることに着目し、到来方向をクラスタリングすることにより、マスクγn(t,f)を推定する。特徴抽出部１１Ａは、各時間周波数点で音がどの方向から到来するかの情報を表す特徴量ベクトルを抽出する。観測信号ベクトルｙ(t,f)の大きさ（ノルム）は音響信号自体に依存するが、観測信号ベクトルｙ(t,f)の方向はその音源位置のみで決まると仮定する。この仮定に基づき、音源位置を表す特徴量ベクトルとして、一定のノルムを持つＭ次元の特徴量ベクトルｚ(t,f)を観測信号ベクトルｙ(t,f)から抽出する。特徴量ベクトルｚ(t,f)は、例えば下記（２）式により抽出できる。
Figure 0006535112
ここで、||・||は２ノルムを表す。当該時間周波数点がどの状態Ｓnを取るかによって、音の到来方向が異なるため、特徴量ベクトルｚ(t,f)は、状態Ｓn毎に異なる確率分布を持つ。
そこで、特徴量ベクトルｚ(t,f)に対してクラスタリング技術を適用することにより、各時間周波数点に各状態がどの程度寄与しているかの割合を表すマスクを推定することができる。これが、上記文献１における基本的な考え方である。
上記文献１では、周波数毎に、特徴量ベクトルz(t,f)の確率分布を、下記（３）式の混合分布でモデル化する。
Figure 0006535112
上記（３）式の混合分布が、特徴量ベクトルｚ(t,f)の分布を近似するようにモデルパラメータ（集合）Θを推定することにより、クラスタリングを行う。上記（３）式の混合分布の要素分布は、下記（４）式で表される。
Figure 0006535112
上記（４）式の要素分布は、状態が既知であるという条件下での特徴量ベクトルｚ(t,f)の条件付き確率分布をモデル化し、α1(f)＋・・・＋αK(f)＝１，０≦αk(f)≦１を満たすαk(f)は、ｋ番目の要素分布の混合重みである。また、||ａk(f)||＝１を満たす位置パラメータ（ベクトル）ａk(f)は、ｋ番目の要素分布の位置を表す。また、広がりパラメータσ2 k(f)は、ｋ番目の要素分布の広がりを表す。なお、Ｋは混合数であり、モデルパラメータ（集合）Θは上記（３）式の混合分布のモデルパラメータの集合｛αk(f),ａk(f),σ2 k(f)｝を表し、・Hは行列のエルミート転置を表す。上記（３）式における要素分布ｐW(ｚ(t,f);ａk(f),σ2 k(f))を山に例えると、上記（３）式の混合分布のモデルパラメータである混合重みαk(f)は山の高さ、位置パラメータａk(f)は山の位置、広がりパラメータσ2 k(f)は山の裾野の広さに該当する。上記文献１では、混合数Ｋは、音響信号の数Ｎ（既知と仮定）と等しいとする。
モデルパラメータΘが求まれば、ベイズの定理より、下記（５）式により、特徴量ベクトルｚ(t,f)が観測された条件下で当該時間周波数点がｋ番目の要素分布に対応する事後確率（すなわちマスク）が求まる。
Figure 0006535112
ここで、モデルパラメータΘをいかに推定するか問題となる。上記文献１では、モデルパラメータΘを用いて、上記（５）式によりマスクγk(t,f)を更新するステップと、マスクγk(t,f)を用いて、モデルパラメータΘを更新するステップの２つのステップを交互に反復することにより、モデルパラメータΘ及びマスクγk(t,f)を推定する。モデルパラメータΘを更新するステップでは、マスクγk(t,f)を用いて、下記（６−１）式〜（６−３）式により、モデルパラメータΘを更新する。
Figure 0006535112
ここで、行列Ｒk(f)は、下記（７）式により計算される。
Figure 0006535112
上記の反復する処理は、期待値最大化法による対数尤度の最適化として、理論的に導出することができる。なお、上記（５）式はマスク更新部１２Ａがマスクγk(t,f)を更新する処理に相当し、上記（６−１）式は混合重み更新部１３Ａが混合重みαk(f)を更新する処理に相当し、上記（６−２）式は位置パラメータ更新部１４Ａが位置パラメータａk(f)を更新する処理に相当し、上記（６−３）式は広がりパラメータ更新部１５Ａが広がりパラメータσ2 k(f)を更新する処理に相当する。
ここで、上記文献１の方法では、上記（３）式の混合分布に基づき上記（５）式によりマスクを推定するため、上記（３）式の混合分布が特徴量ベクトルｚ(t,f)の分布をどの程度正確に近似できるかが、マスクの推定精度に大きく影響する。文献１の方法では、上記（３）式の混合分布を特徴量ベクトルｚ(t,f)の分布に近づけるために、上記（３）式の混合分布の位置パラメータａk(f)と広がりパラメータσ2 k(f)を更新することにより、上記（４）式の要素分布の位置と広がりを調整する。そのため、上記（４）式の要素分布は、ある軸に関して回転対称な分布形状しか表せないという制限がある。
特徴量ベクトルｚ(t,f)の分布形状は、マイクロホン配置や部屋の音響伝達特性等のさまざまな影響を受け、必ずしも回転対称になるとは限らない。そのため、上記文献１の方法では、必ずしも上記（３）式の混合分布が特徴量ベクトルｚ(t,f)の分布を十分近似することができず、上記（５）式によるマスク推定の精度が高くないという問題がある。
そこで、開示の実施形態の第１の態様のマスク推定装置は、Ｎ個（Ｎ＞１）の音響信号が混在する状況で収録したＭ個の観測信号（Ｍ＞１）に基づくＭ次元の特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの、各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータをもとにマスクを更新する。なお、形状とは、図形の属性のうち、位置及び広がりを示す属性を除いた属性とする。形状には、例えば、楕円の長軸及び短軸に相当する情報等がある。
第１の態様のマスク推定装置によれば、形状パラメータの更新により要素分布の形状を調整することで、上記文献１の方法と比べ、混合分布を用いて特徴量ベクトルの分布をより正確に近似することができ、マスクをより正確に推定することができる。
また、開示の実施形態の第２の態様のマスク推定装置は、第１の態様のマスク推定装置において、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの各要素分布を複素ビンガム分布とし、複素ビンガム分布のパラメータ行列を形状パラメータとする。このとき、特徴量ベクトルｚ(t,f)の確率分布ｐ(ｚ(t,f);Θ)は、複素ビンガム分布を要素分布とする、下記（８）式に示す混合分布（以下、混合複素ビンガム分布と呼ぶ）によりモデル化される。
Figure 0006535112
なお、上記（８）式におけるｐB(ｚ;Ｂ)は、下記（９）式により定義される複素ビンガム分布である。なお、exp・は、指数関数を表す。
Figure 0006535112
上記（９）式における行列Ｂは、複素ビンガム分布のパラメータ行列であり、ｃ(Ｂ)は規格化定数である。パラメータΘは、上記（８）式の混合複素ビンガム分布のモデルパラメータの集合｛αk(f),Ｂk(f)｝を表す。
なお、上記（８）式の混合分布のモデルパラメータである混合重みαk(f)は要素分布ｋの高さを表し、行列Ｂk(f)は要素分布ｋの位置、広がりに加えて分布形状（分布形状を表す楕円の各軸方向の分布広がり等）も制御可能な形状パラメータである。行列Ｂk(f)の第一固有ベクトルは要素分布ｋの位置を表し、行列Ｂk(f)の第一固有値と第二固有値の差の絶対値は要素分布ｋの広がりの小ささを表し、行列Ｂk(f)の第一固有値と第ｍ固有値（３≦ｍ≦Ｍ）の差の絶対値は、要素分布ｋの分布形状を表す楕円の第（ｍ−１）軸方向の分布広がりの小ささを表す。
いま、上記（９）式の複素ビンガム分布は、パラメータ行列Ｂがσ2とａの２つのパラメータのみを用いてＢ＝(１/σ2)・ａａHと表される特殊な場合（分布形状を表す楕円の第一軸方向の分布広がり、第二軸方向の分布広がり、・・・・・・、第（Ｍ−１）軸方向の分布広がりがすべて等しい場合）を考えると、下記（１０）式のようになる。
Figure 0006535112
上記（１０）式は、上記（４）式と同じ形をしており、ベクトルａに平行な軸に関して回転対称となる。これに対し、上記（９）式の複素ビンガム分布は、パラメータ行列Ｂに特段の制約は設けられておらず、パラメータ行列Ｂを更新することにより、上記（９）式の複素ビンガム分布の分布形状を調整して回転対称でない分布形状も表現することができる。
これにより、上記（９）式の複素ビンガム分布を用いれば、特徴量ベクトルｚ(t,f)の分布形状が回転対称でない場合でも、上記（８）式の混合分布で特徴量ベクトルｚ(t,f)の分布を十分近似することができる。よって、上記文献１に記載の方法よりも高精度にマスクを推定することが可能になる。
上記（８）式によるモデル化の下、期待値最大化法に基づき、モデルパラメータΘを用いて、マスクγk(t,f)を更新するステップと、マスクγk(t,f)を用いて、モデルパラメータΘを更新するステップの２つのステップを交互に反復することにより、マスクγk(t,f)とモデルパラメータΘを推定するアルゴリズム（上述）を導出することができる。
また、開示の実施形態の第３の態様のマスク推定装置は、第１の態様のマスク推定装置において、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの各要素分布をComplex Angular Central Gaussian (cACG) distribution（以下、複素角度ガウス分布と呼ぶ）とし、複素角度ガウス分布のパラメータ行列を形状パラメータとする。
このとき、特徴量ベクトルｚ(t,f)の確率分布ｐ(ｚ(t,f);Θ)は、複素角度ガウス分布を要素分布とする、下記（１１）式に示す混合分布（以下、混合複素角度ガウス分布と呼ぶ）によりモデル化される。
Figure 0006535112
ここで、上記（１１）式におけるｐA(ｚ;Σ)は、パラメータ行列を行列Σとする、下記（１２）式により定義される複素角度ガウス分布である。
Figure 0006535112
行列Σk(f)は、要素分布ｋの位置、広がりに加えて分布形状（分布形状を表す楕円の各軸方向の分布広がり等）も制御可能な形状パラメータである。行列Σk(f)の第一固有ベクトルは要素分布ｋの位置を表し、行列Σk(f)の第一固有値を第二固有値で割った比は要素分布ｋの広がりの小ささを表し、行列Σk(f)の第一固有値を第ｍ固有値（３≦ｍ≦Ｍ）で割った比は、要素分布ｋの分布形状を表す楕円の第（ｍ−１）軸方向の分布広がりの小ささを表す。モデルパラメータΘは、上記（１１）式における混合複素角度ガウス分布のモデルパラメータの集合｛αk(f);Σk(f)｝を表す。
上記の第２の態様のマスク推定装置における上記（９）式の複素ビンガム分布と同様に、上記（１２）式の複素角度ガウス分布においてもパラメータ行列Σに特段の制約は設けられておらず、回転対称でない分布形状を表現することができる。そのため、上記（１２）式の複素角度ガウス分布を用いれば、特徴量ベクトルｚ(t,f)の分布形状が回転対称でない場合でも分布形状をより正確にモデル化でき、上記文献１の方法よりもより高精度にクラスタリングを行うことが可能になる。
［実施形態１］
実施形態１では、要素分布として複素ビンガム分布を用いる。実施形態１では、上記文献１と同様に、音響信号の数Ｎを既知と仮定し、Ｎ−１個の目的音と、１個の背景雑音が混ざったものをＭ個のマイクロホンで収録するものとする。
（実施形態１に係るマスク推定装置）
図１は、実施形態１に係るマスク推定装置の一例を示す図である。実施形態１に係るマスク推定装置１０は、記憶部２０が接続されている。マスク推定装置１０は、特徴抽出部１１、マスク更新部１２、混合重み更新部１３、パラメータ更新部１４を有する。τを時間を表す番号（サンプル番号）として、マイクロホンｍで収録された時間領域の観測信号をｙm(τ)と書くことにする。ｙm(τ)は、ｎ番目の目的音源信号ｓmn(τ)(ｎ＝１〜Ｎ−１。ｎは正整数)と、背景雑音ｖm(τ)の和として、下記（１３）式によりモデル化できる。
Figure 0006535112
特徴抽出部１１は、複数のマイクロホンで収録された観測信号を受け取り、各観測信号ｙm(τ)毎に時間周波数分析を適用して各観測信号の時間周波数成分ｙm(t,f)を求める（ｍはマイクロホンの番号を表し、１〜Ｍの整数とする）。時間周波数分析としては、短時間フーリエ変換や短時間離散コサイン変換などさまざまな方法を用いることができる。
次に、特徴抽出部１１は、時間周波数点毎にすべての観測信号の時間周波数成分をまとめた、下記（１４）式により表されるＭ次元縦ベクトル（観測信号べクトルと呼ぶ）ｙ(t,f)を構成する。
Figure 0006535112
次に、特徴抽出部１１は、観測信号ベクトルｙ(t,f)を用いて所定のノルムを持つＭ次元の特徴量ベクトルｚ(t,f)を抽出する。特徴量ベクトルｚ(t,f)としては、文献２「H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind sourceseparation via frequency bin-wise clustering and permutation alignment,” IEEE Trans-actions on Audio, Speech and Language Processing(ASLP), vol. 19, no. 3, pp. 516 - 527, Mar. 2011.」、文献３「D.H. Tran Vu and R. Haeb-Umbach, “Blind speech separation em-ploying directional statistics in an expectation maximization framework,” Proceedings of the International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 241 - 244, Mar. 2010.」、文献４「S. Araki, H. Sawada, R. Mukai, and S. Makino, “Underdeter-mined blind sparse source separation for arbitrarily arranged multiple sensors,” Signal Processing, vol. 87, no. 8, pp. 1833 - 1847, Aug. 2007.」に記載の特徴量ベクトル等、さまざまなものを用いることができる。例えば、上記文献２に倣うと、特徴量ベクトルｚ(t,f)は、下記（１５）式により抽出することができる。
Figure 0006535112
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)を受け取り、各時間周波数点(t,f)にｋ番目の要素分布がどの程度寄与しているかの割合を表すマスクγk(t,f)を計算し出力する。
以下、マスク更新部１２における処理の流れを説明する。まず、マスク更新部１２は、記憶部２０にあらかじめ設定し保持しておいたモデルパラメータΘ:=｛αk(f),Ｂk(f)｝の初期値を記憶部２０から読み出す。モデルパラメータΘの初期値の設定は、乱数を用いる等、さまざまな方法により行うことができる。
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)と、モデルパラメータΘの現在の推定値を受け取って、ｋ番目の要素分布に対応するマスクを、例えば、事後確率γk(t,f)として、下記（１６）式により更新する。
Figure 0006535112
ここで、上記（１６）式における行列Ｂk(f)の固有値をｅk1(f),・・・, ｅkM(f)とすると、規格化定数ｃ（Ｂk(f)）は、下記（１７）式のようになる。また、混合数Ｋは、Ｋ＝Ｎとする。
Figure 0006535112
次に、混合重み更新部１３は、事後確率γk(t,f)を受け取って、下記（１８）式により、混合重みを更新する。
Figure 0006535112
次に、パラメータ更新部１４は、特徴量ベクトルｚ(t,f)と、事後確率γk(t,f)を受け取って、正定値エルミート行列Ｒk(f)を、下記（１９）式により計算する。
Figure 0006535112
次に、パラメータ更新部１４は、下記（２０）式のように、行列Ｒk(f)の固有値分解を行う。
Figure 0006535112
上記（２０）式の固有値分解によれば、行列Ｒk(f)の固有ベクトルにより構成されるユニタリ行列Ｕk(f)と、行列Ｒk(f)の固有値により構成される対角行列Ｄk(f)が求まる。ここで、行列Ｄk(f)のｍ番目の対角成分（固有値）ｄkm(f)は、ｄk1(f)≦・・・≦ｄkM(f)のように、昇順に並べてあるものとする。
次に、パラメータ更新部１４は、下記（２１）式により、パラメータ行列Ｂk(f)を更新する。
Figure 0006535112
ここで、上記（２１）式における行列Ｅk(f)は、ｍ番目の対角成分がｅkm(f)である対角行列である。ｅkm(f)は、下記（２２）式で与えられる。
Figure 0006535112
マスク更新部１２、混合重み更新部１３、パラメータ更新部１４での処理は、終了条件が満たされるまで反復して行われる。終了条件としては、例えば、「所定の回数（例えば２０回）反復が行われた」、「更新前後での（後述する）対数尤度関数の増加量が所定の閾値以下であった」等、さまざまな条件を用いることができる。
なお、記憶部２０は、混合重み更新部１３により更新された混合重み、パラメータ更新部１４により更新された形状パラメータを保存し、次回、マスク更新部１２、パラメータ更新部１４での処理の際に、保存した混合重み、形状パラメータを提供する。
（マスク推定装置１０における処理の理論的背景）
以下、マスク推定装置１０における処理の理論的背景について説明する。マスク推定装置１０における処理は、下記（２３）式の対数尤度関数Ｌ(Θ)をモデルパラメータΘに関して最大化する最適化問題を解くことにより導かれるものである。
Figure 0006535112
上記（２３）式の対数尤度関数Ｌ(Θ)は、期待値最大化法に基づいて、以下に述べるＥステップとＭステップを交互に反復することで最適化することができる。
||ｚ||＝１より、任意の実数ｈに対してｐB(ｚ;Ｂ+ｈＩ)＝ｐB(ｚ;Ｂ)（ただし、Ｉは単位行列）が成り立つことから、以下では、一般性を失わず、行列Ｂk(f)の最大固有値ｅkM(f)を０とする。これは、複素ビンガム分布は、各固有値の値そのものではなく、各固有値の間の差分により分布形状が定まるものであることを意味する。このことから、複素ビンガム分布は、最大固有値を指定値に固定する条件を追加したとしても、固有値間の値の差が保たれれば不変であるため、最大固有値ｅkM(f)＝０と固定することができる。
Ｅステップでは、モデルパラメータΘの現在の推定値を用いて、ｋ番目の要素分布の事後確率γk(t,f)を、下記（２４−２）式により更新する。ただし、規格化定数ｃ(・)は、上記（１７）式により定義される。
Figure 0006535112
Ｍステップでは、Ｅステップで更新した事後確率γk(t,f)を用いて、下記（２５−１）式及び（２５−２）式のように定義されるＱ関数を最大化することにより、モデルパラメータΘを更新する。なお、行列Ｒk(f)は、上記（１９）式により定義され、trは行列のトレースを表す。
Figure 0006535112
すなわち、ラグランジュの未定乗数法により、制約条件α1(f)＋・・・＋αK(f)＝１の下で、関数Ｑ(Θ)を最大化するαk(f)が、上記（１８）式のように得られる。以下では、要素分布である複素ビンガム分布の形状パラメータＢk(f)の更新式を導出する。行列Ｒk(f)の固有値分解を下記（２６−１）式のように定義し、行列Ｂk(f)の固有値分解を下記（２６−２）式のように定義する。
Figure 0006535112
ここで、上記（２６−１）式における行列Ｕk(f)は行列Ｒk(f)の固有ベクトルからなるユニタリ行列であり、上記（２６−２）式における行列Ｖk(f)は行列Ｂk(f)の固有ベクトルからなるユニタリ行列である。また、上記（２６−１）式における行列Ｄk(f)は行列Ｒk(f)の固有値からなる対角行列であり、上記（２６−２）式における行列Ｅk(f)は行列Ｂk(f)の固有値からなる対角行列である。
そして、行列Ｒk(f)の固有値をｄk1(f)≦・・・≦ｄkM(f)、行列Ｂ(f)の固有値をｅk1(f)≦・・・≦ｅkM(f)＝０とするとき、行列Ｄk(f)及び行列Ｅk(f)は、それぞれ下記（２７−１）式及び（２７−２）式のように表される。なお、行列Ｂ(f)の最大固有値ｅkM(f)を０とおいた。
なお、ある音源に着目して空間相関行列Ｒk(f)と複素ビンガム分布のパラメータ行列Ｂ(f)を固有値分解してマイクロホンの個数分の固有値を求める際、ｍ番目に大きい第ｍ固有値（ｍ＝１〜Ｍ）は、次の意味を持つ。すなわち、第１固有ベクトルは、着目音源の分布のピークの位置を表す。また、第ｍ固有ベクトル（ｍ＝２〜Ｍ）は、着目音源の分布を表す楕円の第（ｍ−１）軸の方向（行列Ｒk(f)の第ｍ固有値が軸の長さに相当）を表す。
Figure 0006535112
ここで、文献５「K.V. Mardia and P.E. Jupp, “Directional Statistics,” John Wiley & Sons, WestSussex, 2000.」に開示されている結果より、行列Ｂk(f)が関数Ｑ(Θ)を最大化するとき、Ｖk(f)＝Ｕk(f)となる。したがって、上記（２５−２）式において、下記（２８）式のようになることを用いて、ｅkm(f)（ｍ＝１,・・・,Ｍ−１）に関する関数Ｑ(Θ)の偏微分を０と置くと、下記（２９）式を得る。
Figure 0006535112
Figure 0006535112
上記（２９）式をｅkm(f)について解くために、行列Ｒ(f)の固有値ｄk1(f),・・・,ｄkM(f)はすべて異なると仮定する。この仮定は、現実のデータに対してもほとんど常に成り立つ。このとき、上記（１７）式により、規格化定数ｃ（Ｂk(f)）が与えられる。
いま、各音響信号が特定の方向から到来するため複素ビンガム分布は十分に集中していると仮定し、ｅkm(f)＝κk(f)ｅkm0(f)とおいてκk(f)→∞の極限を取ると、下記（３０−１）式〜（３０−３）式のように、規格化定数ｃ（Ｂk(f)）が近似できる。
Figure 0006535112
すなわち、複素ビンガム分布では、「音源の位置が大きく変化せずに音を発生するなら、音源を表現する特徴量ベクトルの分布形状が、比較的急峻なピークを持つこととなる（分布が十分に集中することになる）」ことに着目すると、上記（３０−２）式でのκk(f)は非常に大きな値を有するとの仮定を導入することができ、この結果、規格化定数ｃ（Ｂk(f)）を（３０−３）式の形式で表現できるようになる。したがって、上記（２９）式をｅkm(f)について解くことにより、上記（２２）式が得られる。
（実施形態１に係るマスク推定処理）
図２は、実施形態１に係るマスク推定処理の一例を示すフローチャートである。図２に示すように、先ず、特徴抽出部１１は、複数のマイクロホンで収録された観測信号から、Ｍ次元の特徴量ベクトルｚ(t,f)を抽出する（ステップＳ１１）。
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)、混合重み、及び形状パラメータからマスクγk(t,f)を計算し更新する（ステップＳ１２）。次に、混合重み更新部１３は、事後確率γk(t,f)から混合重みを更新する（ステップＳ１３）。次に、パラメータ更新部１４は、特徴量ベクトルｚ(t,f)と、事後確率γk(t,f)から、パラメータ行列を更新する（ステップＳ１４）。
次に、マスク更新部１２は、終了条件充足か否かを判定する（ステップＳ１５）。マスク更新部１２は、終了条件充足の場合（ステップＳ１５：Ｙｅｓ）、マスク推定処理を終了し、終了条件充足でない場合（ステップＳ１５：Ｎｏ）、ステップＳ１２へ処理を移す。
［実施形態２］
実施形態２では、要素分布として複素角度ガウス分布を用いる。実施形態１では、上記文献１と同様に、音響信号の数Ｎを既知と仮定し、Ｎ−１個の目的音と、１個の背景雑音が混ざったものをＭ個のマイクロホンで収録するものとする。
図１を参照して、実施形態２に係るマスク推定装置１０−２の一例を説明する。マスク推定装置１０−２は、記憶部２０−２が接続されている。マスク推定装置１０−２は、特徴抽出部１１、マスク更新部１２−２、混合重み更新部１３−２、パラメータ更新部１４−２を有する。
特徴抽出部１１は、実施形態１と同様である。マスク更新部１２−２は、特徴量ベクトルｚ(t,f)を受け取り、各時間周波数点(t,f)に、ｋ番目の要素分布がどの程度寄与しているかの割合を表すマスクγk(t,f)を計算し出力する。
以下、マスク更新部１２−２における処理の流れを説明する。まず、マスク更新部１２−２は、記憶部２０−２にあらかじめ設定し保持しておいたモデルパラメータΘ:=｛αk(f),Σk(f)｝の初期値を記憶部２０−２から読み出す。行列Σk(f)は、複素角度ガウス分布のパラメータ行列である。
次に、マスク更新部１２−２は、特徴量ベクトルｚ(t,f)とモデルパラメータΘの現在の推定値を受け取って、ｋ番目の要素分布に対応するマスクを、例えば、事後確率γk(t,f)として、下記（３１）式により更新する。なお、混合数ＫはＫ＝Ｎとする。
Figure 0006535112
次に、混合重み更新部１３−２は、事後確率γk(t,f)を受け取って、下記（３２）式により、混合重みを更新する。
Figure 0006535112
次に、パラメータ更新部１４−２は、特徴量ベクトルｚ(t,f)と、パラメータ行列Σk(t,f)と、事後確率γk(t,f)を受け取って、パラメータ行列Σk(t,f)を、下記（３３）式により更新する。
Figure 0006535112
マスク更新部１２−２、混合重み更新部１３−２、パラメータ更新部１４−２での処理は、実施形態１と同様に、終了条件が満たされるまで反復して行われる。
なお、記憶部２０−２は、混合重み更新部１３−２により更新された混合重み、パラメータ更新部１４−２により更新された形状パラメータを保持し、次回、マスク更新部１２−２、パラメータ更新部１４−２での処理の際に、保存したマスク、形状パラメータを提供する。
（マスク推定装置１０−２における処理の理論的背景）
以下、マスク推定装置１０−２における処理の理論的背景について説明する。マスク推定装置１０−２における処理は、下記（３４）式の対数尤度関数Ｌ(Θ)をモデルパラメータΘに関して最大化する最適化問題を解くことにより導かれるものである。
Figure 0006535112
上記（３４）式の対数尤度関数Ｌ(Θ)は、期待値最大化法に基づいて、以下に述べるＥステップとＭステップを交互に反復することで最適化することができる。
Ｅステップでは、モデルパラメータΘの現在の推定値を用いて、ｋ番目の要素分布の事後確率γk(t,f)を、下記（３５−２）式により更新する。
Figure 0006535112
Ｍステップでは、Ｅステップで更新した事後確率γk(t,f)を用いて、下記（３６−１）式及び（３６−２）式のように定義されるＱ関数を最大化することにより、モデルパラメータΘを更新する。なお、ln・は、対数関数を表す。
Figure 0006535112
すなわち、ラグランジュの未定乗数法により、制約条件α1(f)＋・・・＋αK(f)＝１の下で、関数Ｑ(Θ)を最大化するαk(f)が、上記（３２）式のように得られる。また、関数Ｑ(Θ)のΣk(f)に関する偏微分は、下記（３７）式のようになる。
Figure 0006535112
そして、上記（３７）式で表される関数Ｑ(Θ)のΣk(f)に関する偏微分を０とおくと、上記（３３）式が得られる。
なお、実施形態２に係るマスク推定処理は、実施形態１に係るマスク推定処理と同様に、図２に示すマスク推定処理の一例を示すフローチャートの処理順序に沿って実行される。
［実施形態１及び実施形態２の変形例］
実施形態１及び実施形態２では、Ｎ個の音響信号のうちＮ−１個が目的音であり、１個が背景雑音であるとする。しかし、これに限らず、実施形態１及び実施形態２の変形例では、Ｎ個の音響信号がすべて目的音である場合にも、同一の処理により各音響信号（目的音）に対応するマスクを推定することができる。
また、実施形態１及び実施形態２では、パラメータΘは確定的であるとみなし、最尤法によりパラメータΘを推定する。しかし、これに限らず、実施形態１及び実施形態２の変形例、例えば、パラメータΘを確率変数とみなして事前分布を与え、事後確率最大化によりパラメータΘを推定してもよい。
［実施形態３］
実施形態３は、実施形態１のマスク推定装置１０又は実施形態２のマスク推定装置１０−２のいずれかを用いて目的音を推定する目的音推定システムを開示する。図３は、実施形態３に係る目的音推定システムの一例を示す図である。
実施形態３に係る目的音推定システム１００は、実施形態１のマスク推定装置１０又は実施形態２のマスク推定装置１０−２のいずれか１つと、目的音推定装置３０とを有する。目的音推定装置３０は、観測信号ベクトルｙ(t,f)と、マスク推定装置１０又はマスク推定装置１０−２からのマスクγk(t,f)を入力とし、ｎ番目の目的音の推定値＾ｓn(t,f)（ｎ＝１〜Ｎ−１）を出力する。
目的音推定装置３０は、行列計算部３１、ウィーナーフィルタ構成部３２、目的音推定部３３を有する。
実施形態３では、周波数ｆごとにマスクを求めるため、同一番号ｎの音響信号に対応するクラスタが、異なる周波数ｆでは異なるクラスタ番号を持つということが起こる場合がある。この問題を、パーミュテーション問題という。このパーミュテーション問題のため、実施形態１又は実施形態２において求めたマスクγk(t,f)を用いて各目的音を推定するために、同一番号ｎの音響信号に対応するクラスタが、周波数ｆによらず同じクラスタ番号を持つように対応を取る。これをパーミュテーション解決という。パーミュテーション解決は、上記文献２「H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment,” IEEE Transactions on Audio, Speech and Language Processing(ASLP), vol. 19, no. 3, pp. 516 - 527, Mar. 2011.」に記載の方法等、さまざまな方法により行うことができる。
いま、ｎ番目の音響信号に対応するパーミュテーション解決後のマスクを、改めてγn(t,f)と書くことにする。また、Ｎ個のマスクγn(t,f)のうち、いずれが背景雑音に対応するかが分かっていると仮定する。例えば、マスクを用いてマスキングにより作成した分離音を聞き比べることにより、人為的に背景雑音に対応するクラスタを同定することができる。
以下では、一般性を失わずに、背景雑音に対応するクラスタがＮ番目のクラスタであると仮定する。行列計算部３１は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、下記（３８）式により雑音下目的音共分散行列Φn+v(f)を計算する。
Figure 0006535112
次に、行列計算部３１は、下記（３９）式により雑音共分散行列Φv(f)を計算する。
Figure 0006535112
次に、行列計算部３１は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、下記（４０）式により目的音共分散行列Φn(f)を求める。
Figure 0006535112
次に、行列計算部３１は、下記（４１）式により、観測共分散行列Φy(f)を求める。
Figure 0006535112
次に、ウィーナーフィルタ構成部３２は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、マルチチャネル・ウィーナー・フィルタＷn(f)を、下記（４２）式により求める。
Figure 0006535112
そして、目的音推定部３３は、下記（４３）式のように、マルチチャネル・ウィーナー・フィルタＷn(f)を観測信号ベクトルｙ(t,f)に適用することで、背景雑音や、ｎ番目以外の目的音の成分を抑圧し、目的音ｎの成分の推定値＾ｓn(t,f)を得ることができる。
Figure 0006535112
（実施形態３に係る目的音推定処理）
図４は、実施形態３に係る目的音推定処理の一例を示すフローチャートである。先ず、目的音推定装置３０の行列計算部３１は、マスク推定装置１０（１０−２）からマスク情報を取得する（ステップＳ２１）。
次に、行列計算部３１は、雑音下目的音共分散行列Φn+v(f)を計算する（ステップＳ２２）。次に、行列計算部３１は、雑音共分散行列Φv(f)を計算する（ステップＳ２３）。次に、行列計算部３１は、目的音共分散行列Φn(f)を計算する（ステップＳ２４）。次に、行列計算部３１は、観測共分散行列Φy(f)を計算する（ステップＳ２５）。
次に、ウィーナーフィルタ構成部３２は、マルチチャネル・ウィーナー・フィルタＷn(f)を構成する（ステップＳ２６）。次に、目的音推定部３３は、ステップＳ２６で構成したマルチチャネル・ウィーナー・フィルタＷn(f) を観測信号ベクトルｙ(t,f)に適用することで、目的音ｎの成分の推定値＾ｓn(t,f)を得て、出力する（ステップＳ２７）。
（性能実証結果）
実施形態３の性能実証として、実施形態１のマスク推定装置１０を含む目的音推定システム１００を用い、バスの中、カフェなどの背景雑音が存在する環境において、１人の話者がタブレットに向かって文章を読み上げている状況（実施形態１においてＮ＝２の場合）で、タブレットに装着されたＭ＝６個のマイクロホンで収録した信号に対して音声認識を行った。音声認識性能（単語誤り率）は、次のとおりであった。
すなわち、マスク推定を行わず音声認識をした場合の単語誤り率は、14.29(%)であった。また、混合複素ワトソン分布でマスク推定をした後、マルチチャネル・ウィーナー・フィルタを適用した音声認識の場合の単語誤り率は、9.51(%)であった。実施形態１のマスク推定装置１０を含む目的音推定システム１００において、混合複素ビンガム分布でマスク推定した後、マルチチャネル・ウィーナー・フィルタを適用した場合の単語誤り率は、8.53(%)であった。以上から、実施形態３は、従来技術と比較して、音声認識性能が向上することが分かる。
上記の実施形態のマスク推定装置１０（１０−２）及び目的音推定装置３０を含む目的音推定システム１００において行われる各処理は、全部又は任意の一部が、ＣＰＵ（Central Processing Unit）等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、マスク推定装置１０（１０−２）及び目的音推定装置３０を含む目的音推定システム１００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。
また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。また、実施形態においてフローチャートを参照して説明した各処理は、最終的な実行結果が同一である限りにおいて、適宜実行順序の入れ替えや並列実行を行うことができる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。
（プログラムについて）
図５は、プログラムが実行されることにより、実施形態に係るマスク推定装置及び目的音推定装置を含む目的音推定システムが実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。
メモリ１０１０は、ＲＯＭ１０１１及びＲＡＭ１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。
ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわちマスク推定装置１０（１０−２）及び目的音推定装置３０の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、マスク推定装置１０（１０−２）及び目的音推定装置３０における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。
また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
１０、１０−２、１０Ａ マスク推定装置
１１、１１Ａ 特徴抽出部
１２、１２−２、１２Ａ マスク更新部
１３、１３−２、１３Ａ 混合重み更新部
１４、１４−２ パラメータ更新部
１４Ａ 位置パラメータ更新部
１５Ａ 広がりパラメータ更新部
２０、２０−２、２０Ａ 記憶部
３０ 目的音推定装置
３１ 行列計算部
３２ ウィーナーフィルタ構成部
３３ 目的音推定部
１０００ コンピュータ
１０１０ メモリ
１０２０ ＣＰＵ

Claims (11)

  1. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記要素分布は、回転対称ではない分布形状も表現することができることを特徴とするマスク推定装置。
  2. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記要素分布の確率密度関数は、ｚを独立変数、ａを単位ベクトル、・Ｈをエルミート転置として、ａとｚの内積ａＨｚの絶対値｜ａＨｚ｜の関数ではない確率密度関数も表現することができることを特徴とするマスク推定装置。
  3. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記要素分布の確率密度関数は、ｚを独立変数、Ａを行列、・Ｈをエルミート転置として、二次形式ｚＨＡｚの関数であり、
    前記形状パラメータは、前記行列Ａを定めるパラメータである
    ことを特徴とするマスク推定装置。
  4. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記形状パラメータは、各要素分布の分布形状を表す楕円の各軸方向の分布広がりの情報を含む
    ことを特徴とするマスク推定装置。
  5. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記要素分布は、複素ビンガム分布であり、
    前記形状パラメータは、複素ビンガム分布のパラメータ行列である
    ことを特徴とするマスク推定装置。
  6. 複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
    前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
    前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
    を備え、
    前記要素分布は、複素角度ガウス分布であり、
    前記形状パラメータは、複素角度ガウス分布のパラメータ行列である
    ことを特徴とするマスク推定装置。
  7. マスク推定装置が実行するマスク推定方法であって、
    前記マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめることにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部を備え、
    前記特徴量ベクトルを抽出する特徴抽出ステップと、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新ステップと、
    前記マスク更新ステップにより更新されたマスクから、前記混合重みを更新する混合重み更新ステップと、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新ステップと
    を含み、
    前記要素分布は、回転対称ではない分布形状も表現することができることを特徴とするマスク推定方法。
  8. マスク推定装置が実行するマスク推定方法であって、
    前記マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめることにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部を備え、
    前記特徴量ベクトルを抽出する特徴抽出ステップと、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新ステップと、
    前記マスク更新ステップにより更新されたマスクから、前記混合重みを更新する混合重み更新ステップと、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新ステップと
    を含み、
    前記要素分布の確率密度関数は、ｚを独立変数、ａを単位ベクトル、・Ｈをエルミート転置として、ａとｚの内積ａＨｚの絶対値｜ａＨｚ｜の関数ではない確率密度関数も表現することができることを特徴とするマスク推定方法。
  9. マスク推定装置が実行するマスク推定方法であって、
    前記マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめることにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部を備え、
    前記特徴量ベクトルを抽出する特徴抽出ステップと、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新ステップと、
    前記マスク更新ステップにより更新されたマスクから、前記混合重みを更新する混合重み更新ステップと、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新ステップと
    を含み、
    前記要素分布の確率密度関数は、ｚを独立変数、Ａを行列、・Ｈをエルミート転置として、二次形式ｚＨＡｚの関数であり、
    前記形状パラメータは、前記行列Ａを定めるパラメータである
    ことを特徴とするマスク推定方法。
  10. マスク推定装置が実行するマスク推定方法であって、
    前記マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめることにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部を備え、
    前記特徴量ベクトルを抽出する特徴抽出ステップと、
    前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新ステップと、
    前記マスク更新ステップにより更新されたマスクから、前記混合重みを更新する混合重み更新ステップと、
    前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新ステップと
    を含み、
    前記形状パラメータは、各要素分布の分布形状を表す楕円の各軸方向の分布広がりの情報を含む
    ことを特徴とするマスク推定方法。
  11. 請求項１〜６のいずれか１つに記載のマスク推定装置としてコンピュータを機能させるためのマスク推定プログラム。
JP2017567967A 2016-02-16 2016-12-20 マスク推定装置、マスク推定方法及びマスク推定プログラム Active JP6535112B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2016027424 2016-02-16
JP2016027424 2016-02-16
PCT/JP2016/087996 WO2017141542A1 (ja) 2016-02-16 2016-12-20 マスク推定装置、マスク推定方法及びマスク推定プログラム

Publications (2)

Publication Number Publication Date
JPWO2017141542A1 JPWO2017141542A1 (ja) 2018-07-12
JP6535112B2 true JP6535112B2 (ja) 2019-06-26

Family

ID=59625834

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2017567967A Active JP6535112B2 (ja) 2016-02-16 2016-12-20 マスク推定装置、マスク推定方法及びマスク推定プログラム

Country Status (4)

Country Link
US (1) US10878832B2 (ja)
JP (1) JP6535112B2 (ja)
CN (1) CN108701468B (ja)
WO (1) WO2017141542A1 (ja)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11423924B2 (en) * 2018-02-23 2022-08-23 Nippon Telegraph And Telephone Corporation Signal analysis device for modeling spatial characteristics of source signals, signal analysis method, and recording medium
JP6915579B2 (ja) * 2018-04-06 2021-08-04 日本電信電話株式会社 信号分析装置、信号分析方法および信号分析プログラム
JP6992709B2 (ja) * 2018-08-31 2022-01-13 日本電信電話株式会社 マスク推定装置、マスク推定方法及びマスク推定プログラム
CN109859769B (zh) * 2019-01-30 2021-09-17 西安讯飞超脑信息科技有限公司 一种掩码估计方法及装置
CN110674528B (zh) * 2019-09-20 2024-04-09 深圳前海微众银行股份有限公司 联邦学习隐私数据处理方法、设备、系统及存储介质
CN113539290B (zh) * 2020-04-22 2024-04-12 华为技术有限公司 语音降噪方法和装置
CN112564885B (zh) * 2020-11-26 2022-07-12 南京农业大学 基于掩码变量最大概率密度函数分布的侧信道测试分析方法
JP7487795B2 (ja) 2020-12-14 2024-05-21 日本電信電話株式会社 音源信号生成装置、音源信号生成方法、プログラム
US11755888B1 (en) * 2023-01-09 2023-09-12 Fudan University Method and system for accelerating score-based generative models with preconditioned diffusion sampling

Family Cites Families (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6816632B1 (en) * 2000-02-17 2004-11-09 Wake Forest University Health Sciences Geometric motion analysis
KR100647286B1 (ko) * 2004-08-14 2006-11-23 삼성전자주식회사 교차채널 간섭을 제거하기 위한 후처리장치 및 방법과이를 이용한 다채널 음원 분리장치 및 방법
TWI257482B (en) * 2004-12-15 2006-07-01 Spirox Corp Method and apparatus for measuring jitter of signal
JP2006337851A (ja) * 2005-06-03 2006-12-14 Sony Corp 音声信号分離装置及び方法
JP4675177B2 (ja) * 2005-07-26 2011-04-20 株式会社神戸製鋼所 音源分離装置，音源分離プログラム及び音源分離方法
JP2007156300A (ja) * 2005-12-08 2007-06-21 Kobe Steel Ltd 音源分離装置、音源分離プログラム及び音源分離方法
US8433199B2 (en) * 2008-03-18 2013-04-30 Princeton University System and method for nonlinear self-filtering via dynamical stochastic resonance
JP5480496B2 (ja) * 2008-03-25 2014-04-23 株式会社ニューフレアテクノロジー 荷電粒子ビーム描画方法及び荷電粒子ビーム描画装置
JP5277887B2 (ja) * 2008-11-14 2013-08-28 ヤマハ株式会社 信号処理装置およびプログラム
KR20110072630A (ko) * 2009-12-23 2011-06-29 삼성전자주식회사 빔의 위치 측정 장치 및 방법
JP5728888B2 (ja) * 2010-10-29 2015-06-03 ソニー株式会社 信号処理装置および方法、並びにプログラム
JP5621637B2 (ja) * 2011-02-04 2014-11-12 ヤマハ株式会社 音響処理装置
US9291725B2 (en) * 2012-05-16 2016-03-22 Kabushiki Kaisha Toshiba Random coincidence reduction in positron emission tomography using tangential time-of-flight mask
CN103594093A (zh) * 2012-08-15 2014-02-19 王景芳 基于信噪比软掩蔽语音增强方法
JP6253226B2 (ja) * 2012-10-29 2017-12-27 三菱電機株式会社 音源分離装置
US10475440B2 (en) * 2013-02-14 2019-11-12 Sony Corporation Voice segment detection for extraction of sound source
JP6059072B2 (ja) * 2013-04-24 2017-01-11 日本電信電話株式会社 モデル推定装置、音源分離装置、モデル推定方法、音源分離方法及びプログラム
US9497528B2 (en) * 2013-11-07 2016-11-15 Continental Automotive Systems, Inc. Cotalker nulling based on multi super directional beamformer
CN105096961B (zh) * 2014-05-06 2019-02-01 华为技术有限公司 语音分离方法和装置
US20160216384A1 (en) * 2015-01-26 2016-07-28 Brimrose Technology Corporation Detection of nuclear radiation via mercurous halides

Also Published As

Publication number Publication date
WO2017141542A1 (ja) 2017-08-24
CN108701468A (zh) 2018-10-23
US10878832B2 (en) 2020-12-29
CN108701468B (zh) 2023-06-02
US20190267019A1 (en) 2019-08-29
JPWO2017141542A1 (ja) 2018-07-12

Similar Documents

Publication Publication Date Title
JP6535112B2 (ja) マスク推定装置、マスク推定方法及びマスク推定プログラム
US10643633B2 (en) Spatial correlation matrix estimation device, spatial correlation matrix estimation method, and spatial correlation matrix estimation program
US10650841B2 (en) Sound source separation apparatus and method
US11456003B2 (en) Estimation device, learning device, estimation method, learning method, and recording medium
CN111415676A (zh) 一种基于分离矩阵初始化频点选择的盲源分离方法及系统
JP6517760B2 (ja) マスク推定用パラメータ推定装置、マスク推定用パラメータ推定方法およびマスク推定用パラメータ推定プログラム
Drude et al. Unsupervised training of neural mask-based beamforming
JP7176627B2 (ja) 信号抽出システム、信号抽出学習方法および信号抽出学習プログラム
JP2018141922A (ja) ステアリングベクトル推定装置、ステアリングベクトル推定方法およびステアリングベクトル推定プログラム
JP6538624B2 (ja) 信号処理装置、信号処理方法および信号処理プログラム
EP3550565A1 (en) Audio source separation with source direction determination based on iterative weighting
JP5726790B2 (ja) 音源分離装置、音源分離方法、およびプログラム
JP6910609B2 (ja) 信号解析装置、方法、及びプログラム
JP6370751B2 (ja) ガウス混合モデルパラメータ計算装置、情報推定装置、音強調装置、これらの方法及びプログラム
JP6636973B2 (ja) マスク推定装置、マスク推定方法およびマスク推定プログラム
JP2016156944A (ja) モデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラム
JP6734237B2 (ja) 目的音源推定装置、目的音源推定方法及び目的音源推定プログラム
JP6114053B2 (ja) 音源分離装置、音源分離方法、およびプログラム
KR101647059B1 (ko) 독립 벡터 분석 및 모델 기반 특징 향상을 이용한 강인한 음성 인식 방법
Loweimi et al. On the usefulness of statistical normalisation of bottleneck features for speech recognition
JP6915579B2 (ja) 信号分析装置、信号分析方法および信号分析プログラム
Rafique et al. Speech source separation using the IVA algorithm with multivariate mixed super gaussian student's t source prior in real room environment
CN109074811A (zh) 音频源分离
Vincent Advances in audio source seperation and multisource audio content retrieval
JP2023039288A (ja) 音源分離モデル学習装置、音源分離装置、音源分離モデル学習方法、音源分離方法及びプログラム

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20180221

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20181218

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20190213

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20190528

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20190530

R150 Certificate of patent or registration of utility model

Ref document number: 6535112

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150