(57)【要約】  信号分析装置（１）は、Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列である同時無相関化行列Ｐまたは／及び、そのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る信号パラメータ分析部（２０）を有する。
2019163487.tif 000002
【特許請求の範囲】
【請求項１】
Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列である同時無相関化行列Ｐまたは／及び、そのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る信号パラメータ分析部
を有することを特徴とする信号分析装置。
【請求項２】
請求項１に記載の信号分析装置であって、
前記信号パラメータ分析部は、前記空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列に近づくように求めた行列である同時無相関化行列Ｐまたは／及びそのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得ることを特徴とする信号分析装置。
【請求項３】
請求項１に記載の信号分析装置であって、
前記信号パラメータ分析部は、不動点法、勾配法、自然勾配法、または、ニュートン法に基づいて、同時無相関化行列Ｐの推定値を更新することを特徴とする信号分析装置。
【請求項４】
請求項１に記載の信号分析装置であって、
前記信号パラメータ分析部が得た前記同時無相関化行列Ｐのエルミート転置ＰＨを、入力されたＩ個の観測信号による観測信号ベクトルに掛けることにより、Ｊ個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る同時無相関化部
をさらに有することを特徴とする信号分析装置。
【請求項５】
請求項４に記載の信号分析装置であって、
前記同時無相関化部が得た前記同時無相関化済観測信号ベクトルの各要素に、
前記同時無相関化行列Ｐに基づく行列ＰＨＲｊＰ（ｊは１以上Ｊ以下の整数）の対角成分に基づくシングルチャネルウィーナーフィルタを適用して同時無相関化済源信号成分推定値によるベクトルを得るウィーナーフィルタ部
をさらに有することを特徴とする信号分析装置。
【請求項６】
請求項５に記載の信号分析装置であって、
前記ウィーナーフィルタ部が得た前記同時無相関化済源信号成分推定値によるベクトルに、前記同時無相関化行列Ｐのエルミート転置ＰＨの逆行列を掛けることにより、源信号成分の推定値を得る同時無相関化逆変換部
をさらに有することを特徴とする信号分析装置。
【請求項７】
２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒ１，Ｒ２に対するｐｉＨＲ２ｐｍ＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする行列である同時無相関化行列Ｐまたはそのエルミート転置ＰＨを、前記Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る信号パラメータ分析部
を有することを特徴とする信号分析装置。
【請求項８】
請求項７に記載の信号分析装置であって、
前記信号パラメータ分析部が得た前記同時無相関化行列Ｐのエルミート転置ＰＨを、入力されたＩ個の観測信号による観測信号ベクトルに掛けることにより、２個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る同時無相関化部
をさらに有することを特徴とする信号分析装置。
【請求項９】
請求項８に記載の信号分析装置であって、
前記同時無相関化部が得た前記同時無相関化済観測信号ベクトルの各要素に、
前記同時無相関化行列Ｐに基づく行列ＰＨＲｊＰ（ｊは１または２）の対角成分に基づくシングルチャネルウィーナーフィルタを適用して同時無相関化済源信号成分推定値によるベクトルを得るウィーナーフィルタ部
をさらに有することを特徴とする信号分析装置。
【請求項１０】
請求項９に記載の信号分析装置であって、
前記ウィーナーフィルタ部が得た前記同時無相関化済源信号成分推定値によるベクトルに、前記同時無相関化行列Ｐのエルミート転置ＰＨの逆行列を掛けることにより、源信号成分の推定値を得る同時無相関化逆変換部
をさらに有することを特徴とする信号分析装置。
【請求項１１】
Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列Ｐを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するための同時無相関化行列として、前記同時無相関化行列のエルミート転置を、前記観測信号ベクトルに掛けることにより、Ｊ個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る同時無相関化部
を有することを特徴とする信号分析装置。
【請求項１２】
請求項１１に記載の信号分析装置であって、
前記同時無相関化部が得た前記同時無相関化済観測信号ベクトルの各要素に、
前記同時無相関化行列Ｐに基づく行列ＰＨＲｊＰ（ｊは１以上Ｊ以下の整数）の対角成分に基づくシングルチャネルウィーナーフィルタを適用して同時無相関化済源信号成分推定値によるベクトルを得るウィーナーフィルタ部
をさらに有することを特徴とする信号分析装置。
【請求項１３】
請求項１２に記載の信号分析装置であって、
前記ウィーナーフィルタ部が得た前記同時無相関化済源信号成分推定値によるベクトルに、前記同時無相関化行列Ｐのエルミート転置ＰＨの逆行列を掛けることにより、源信号成分の推定値を得る同時無相関化逆変換部
をさらに有することを特徴とする信号分析装置。
【請求項１４】
Ｊ個（Ｊは２以上の整数）の混在する源信号から得た同時無相関化された共分散行列に基づいて、源信号の空間的特性をモデル化するパラメータである空間共分散行列を求める空間共分散行列計算部、または／及び前記同時無相関化された共分散行列に基づいて、源信号のサンプル点ごとの分散をモデル化するパラメータである分散パラメータを求める分散パラメータ計算部
を有することを特徴とする信号分析装置。
【請求項１５】
請求項１４に記載の信号分析装置であって、
補助関数法を用いて、前記空間共分散行列または／及び分散パラメータを求めることを特徴とする信号分析装置。
【請求項１６】
信号分析装置が実行する信号分析方法であって、
Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列である同時無相関化行列Ｐまたは／及び、そのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る工程
を含んだことを特徴とする信号分析方法。
【請求項１７】
信号分析装置が実行する信号分析方法であって、
２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒ１，Ｒ２に対するｐｉＨＲ２ｐｍ＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする行列である同時無相関化行列Ｐまたはそのエルミート転置ＰＨを、前記Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る工程
を含んだことを特徴とする信号分析方法。
【請求項１８】
信号分析装置が実行する信号分析方法であって、
Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列Ｐを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するための同時無相関化行列として、前記同時無相関化行列のエルミート転置を、前記観測信号ベクトルに掛けることにより、Ｊ個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る工程
を含んだことを特徴とする信号分析方法。
【請求項１９】
信号分析装置が実行する信号分析方法であって、
Ｊ個（Ｊは２以上の整数）の混在する源信号から得た同時無相関化された共分散行列に基づいて、源信号の空間的特性をモデル化するパラメータである空間共分散行列を求める工程
を含んだことを特徴とする信号分析方法。
【請求項２０】
コンピュータを、請求項１〜１５のいずれか一つに記載の信号分析装置として機能させるための信号分析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号分析装置、信号分析方法及び信号分析プログラムに関する。
【背景技術】
【０００２】
従来、２個の源信号が混在する状況において、それぞれ異なる位置で取得された複数の観測信号から、各源信号に由来する成分である源信号成分を推定する方法が提案されている。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】N.Q.K.  Duong,  E.  Vincent  and  R.  Gribonval,  “Under-Determined  Reverberant  Audio  Source  Separation  Using  a  Full-Rank  Spatial  Covariance  Model”,  IEEE  Transactions  on  Audio,  Speech,  and  Language  Processing,  vol.18,  no.  7,  pp.  1830−1840,  Sep.  2010.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ここで、図２４を用いて、従来の信号分析装置について説明する。図２４は、従来の信号分析装置の構成を示す図である。図２４に示すように、従来の信号分析装置１Ｐは、マルチチャネルウィーナーフィルタ部１０ａを有する。なお、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。【０００５】
マルチチャネルウィーナーフィルタ部１０ａは、観測信号ｙｉ（ｋ）（ｉは観測信号のインデックス、ｋはサンプル点のインデックス）と、源信号のサンプル点ごとの分散をモデル化するパラメータである分散パラメータｖｊ（ｋ）（ｊは源信号のインデックス）と、源信号の空間的特性をモデル化するパラメータである空間共分散行列Ｒｊと、に基づいて、源信号成分の推定値を計算する。【０００６】
具体的には、マルチチャネルウィーナーフィルタ部１０ａは、すべての観測信号を要素とするＩ次元縦ベクトルである観測信号ベクトルｙ（ｋ）を次の（１）式により作成する。なお、以降において、上付きのＴは、転置を表す。【０００７】
【数１】
2019163487.tif 000003
【０００８】
観測信号ベクトルｙ（ｋ）は、２個の源信号に由来する成分（源信号成分）ｘ１（ｋ），ｘ２（ｋ）の和として、次の（２）式によってモデル化される。【０００９】
【数２】
2019163487.tif 000004
【００１０】
マルチチャネルウィーナーフィルタ部１０ａは、さらに、観測信号ベクトルｙ（ｋ）にＩ×Ｉの行列であるマルチチャネルウィーナーフィルタＷｊ（ｋ）を適用することにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を次の（３）式により計算する。【００１１】
【数３】
2019163487.tif 000005
【００１２】
しかしながら、従来の信号分析装置１Ｐでは、大きい計算量を必要とする逆行列演算をサンプル点ごとに行う必要があった。具体的には、従来の信号分析装置１Ｐでは、ｋごとに逆行列（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１を計算する必要があった。このため、従来の信号分析装置１Ｐでは、源信号に由来する成分である源信号成分の推定値を計算するために大きい計算量が必要であった。【００１３】
本発明は、上記に鑑みてなされたものであって、源信号成分の推定値を少ない計算量で計算することができる信号分析装置、信号分析方法及び信号分析プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１４】
上述した課題を解決し、目的を達成するために、本発明の信号分析装置は、Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列である同時無相関化行列Ｐまたは／及び、そのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る信号パラメータ分析部を有することを特徴とする。【００１５】
また、本発明の信号分析装置は、２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒ１，Ｒ２に対するｐｉＨＲ２ｐｍ＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする行列である同時無相関化行列Ｐまたはそのエルミート転置ＰＨを、Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る信号パラメータ分析部を有することを特徴とする。【００１６】
また、本発明の信号分析装置は、Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列Ｐを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するための同時無相関化行列として、同時無相関化行列のエルミート転置を、観測信号ベクトルに掛けることにより、Ｊ個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る同時無相関化部を有することを特徴とする。【００１７】
また、本発明の信号分析装置は、Ｊ個（Ｊは２以上の整数）の混在する源信号から得た同時無相関化された共分散行列に基づいて、源信号の空間的特性をモデル化するパラメータである空間共分散行列を求める空間共分散行列計算部、または／及び同時無相関化された共分散行列に基づいて、源信号のサンプル点ごとの分散をモデル化するパラメータである分散パラメータを求める分散パラメータ計算部を有することを特徴とする。
【発明の効果】
【００１８】
本発明によれば、源信号成分の推定値を少ない計算量で計算することができる。
【図面の簡単な説明】
【００１９】
【図１】図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。
【図２】図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図３】図３は、第２の実施形態に係る信号分析装置の構成の一例を示す図である。
【図４】図４は、第２の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図５】図５は、第３の実施形態に係る信号分析装置の構成の一例を示す図である。
【図６】図６は、第３の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図７】図７は、従来の方法及び第３の実施形態を用いた音源分離実験の実験条件を示す図である。
【図８】図８は、従来の方法及び第３の実施形態を用いた音源分離実験結果（real  time  factor）を示す図である。
【図９】図９は、従来の方法及び第３の実施形態を用いた音源分離実験結果（音源分離性能）を示す図である。
【図１０】図１０は、第４の実施形態に係る信号分析装置の構成の一例を示す図である。
【図１１】図１１は、第４の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図１２】図１２は、第５の実施形態に係る信号分析装置の構成の一例を示す図である。
【図１３】図１３は、第５の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図１４】図１４は、従来の方法及び第５の実施形態を用いた音源分離実験の実験条件を示す図である。
【図１５】図１５は、従来の方法及び第５の実施形態を用いた音源分離実験結果（real  time  factor）を示す図である。
【図１６】図１６は、従来の方法及び第５の実施形態を用いた音源分離実験結果（音源分離性能）を示す図である。
【図１７】図１７は、第６の実施形態に係る信号分析装置の構成の一例を示す図である。
【図１８】図１８は、第６の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図１９】図１９は、周波数領域信号分離部の構成の一例を示す図である。
【図２０】図２０は、周波数領域信号分離処理の処理手順の一例を示すフローチャートである。
【図２１】図２１は、第８の実施形態に係る信号分析装置の構成の一例を示す図である。
【図２２】図２２は、第８の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【図２３】図２３は、プログラムが実行されることにより、信号分析装置が実現されるコンピュータの一例を示す図である。
【図２４】図２４は、従来の信号分析装置の構成を示す図である。
【発明を実施するための形態】
【００２０】
以下に、本願に係る信号分析装置、信号分析方法及び信号分析プログラムの実施形態を図面に基づいて詳細に説明する。また、本発明は、以下に説明する実施形態により限定されるものではない。なお、以下では、ベクトル、行列又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同じであるとする。ベクトル、行列又はスカラーであるＡに対し、“~Ａ”と記載する場合は「“Ａ”の直上に“~”が記された記号」と同じであるとする。また、以下では、信号分離は、信号分析の下位概念とする。
【００２１】
［第１の実施形態］
まず、第１の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。なお、第１の実施形態においては、２個の源信号が混在する状況において、それぞれ異なる位置でマイクロホンなどのセンサにより取得されたＩ個（Ｉは２以上の整数）の観測信号ｙｉ（ｋ）（ｉ＝１，・・・，Ｉは観測信号のインデックス、ｋはサンプル点のインデックス）と、源信号のサンプル点ごとの分散（パワー）をモデル化するパラメータである分散パラメータｖｊ（ｋ）（ｊ＝１，２は源信号のインデックス）と、源信号の空間的特性をモデル化するパラメータである空間共分散行列Ｒｊと、が信号分析装置に入力されるものとする。
【００２２】
なお、本第１の実施形態における「信号」の例としては、音、脳波、無線信号などが挙げられる。本第１の実施形態における「観測信号」とは、複数のセンサ（センサアレイ）により取得された信号を指す。センサアレイの例としては、マイクロホンアレイ、脳波計・脳磁計、アンテナアレイなどが挙げられる。また、本第１の実施形態における「源信号」とは、観測信号の構成要素である信号を指す。観測信号は２個の源信号の和としてモデル化される。
【００２３】
また、本第１の実施形態における「信号分離」の問題設定は、２個の目的信号（たとえば音声）が混在する状況においてこれらに対応する成分を推定する問題設定（たとえば音源分離）に加えて、目的信号（たとえば音声）とノイズ（たとえばテレビから流れる音楽や雑踏における喧噪）が混在する状況においてこれらに対応する成分（または目的信号に対応する成分のみ）を推定する問題設定（ノイズ除去とも呼ばれる）をも含むものとする。前者の場合における源信号は、目的信号（たとえば音声）である。後者の場合における源信号は、目的信号（たとえば音声）とノイズ（たとえばテレビから流れる音楽や雑踏における喧噪）とである。なお、「雑踏における喧噪」は、多数のノイズが混ざったものであるが、この場合まとめて１つの源信号とみなす。
【００２４】
［信号分析装置の構成及び処理］
まず、図１及び図２を用いて、第１の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図１は、第１の実施形態に係る信号分析装置の構成の一例を示す図である。図２は、第１の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。第１の実施形態に係る信号分析装置１は、例えば、ＲＯＭ（Read  Only  Memory）、ＲＡＭ（Random  Access  Memory）、ＣＰＵ（Central  Processing  Unit）等を含むコンピュータ等に所定のプログラムが読み込まれて、ＣＰＵが所定のプログラムを実行することで実現される。
【００２５】
図１に示すように、信号分析装置１は、観測信号ベクトル作成部１０、信号パラメータ分析部２０、同時無相関化部３０、シングルチャネルウィーナーフィルタ部４０及び同時無相関化逆変換部５０を有する。
【００２６】
まず、信号分析装置１の各部の概要について説明する。観測信号ベクトル作成部１０は、入力された観測信号ｙｉ（ｋ）を取得し（ステップＳ１０）、（４）式に示す、すべての観測信号からなるＩ次元縦ベクトルである観測信号ベクトルをサンプル点ごとに作成する（ステップＳ１１）。
【００２７】
【数４】
2019163487.tif 000006
【００２８】
信号パラメータ分析部２０は、入力された空間共分散行列Ｒｊを取得し、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Ｐを計算する（ステップＳ１２）。
【００２９】
同時無相関化部３０は、観測信号ベクトル作成部１０からの観測信号ベクトルｙ（ｋ）と、信号パラメータ分析部２０からの一般化固有ベクトル行列Ｐと、に基づいて、同時無相関化された観測信号ベクトルｙ´（ｋ）を計算する（ステップＳ１３）。具体的には、同時無相関化部３０は、（５）式に示すように、観測信号ベクトルｙ（ｋ）に、一般化固有ベクトル行列Ｐのエルミート転置ＰＨを掛けることにより、同時無相関化された観測信号ベクトルｙ´（ｋ）を計算する。
【００３０】
【数５】
2019163487.tif 000007
【００３１】
ここで、「´（上付きのプライム）」は、同時無相関化を施された変数であることを表す。なお、同時無相関化とは、観測信号ベクトルｙ（ｋ）に含まれる、２個の源信号に対応する成分（源信号成分）を同時に無相関化する線型変換を指す。
【００３２】
シングルチャネルウィーナーフィルタ部４０は、入力された分散パラメータｖｊ（ｋ）を取得し、信号パラメータ分析部２０からの一般化固有値行列Λと、同時無相関化部３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、に基づいて、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を計算する（ステップＳ１４）。具体的には、シングルチャネルウィーナーフィルタ部４０は、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を（６），（７）式により計算する。
【００３３】
【数６】
2019163487.tif 000008
【数７】
2019163487.tif 000009
【００３４】
ただし、ベクトルαに対し、［α］ｌはαの第ｌ要素を表し、行列Ａに対し、［Ａ］ｌｍはＡの（ｌ，ｍ）要素を表す。また、変数の上の∧（ハット）は、推定値であることを表す。
【００３５】
同時無相関化逆変換部５０は、信号パラメータ分析部２０からの一般化固有ベクトル行列Ｐと、シングルチャネルウィーナーフィルタ部４０からの同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）と、に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算し、出力する（ステップＳ１５）。具体的には、同時無相関化逆変換部５０は、（８）式に示すように、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）に対し、同時無相関化部３０で掛けた一般化固有ベクトル行列Ｐのエルミート転置ＰＨの逆行列（ＰＨ）−１を掛けることにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算し、出力する。
【００３６】
【数８】
2019163487.tif 000010
【００３７】
［第１の実施形態の処理の背景］
次に、第１の実施形態の処理の背景について説明する。なお、以下では、観測信号が複素数で表される場合について説明する（観測信号が実数で表される場合の説明は、ほぼ同様だから省略する）。第１の実施形態においては、観測信号ベクトルｙ（ｋ）は、２個の源信号に対応する成分（源信号成分）ｘ１（ｋ），ｘ２（ｋ）の和として、次に示す（９）式でモデル化される。
【００３８】
【数９】
2019163487.tif 000011
【００３９】
源信号成分ｘｊ（ｋ）は、平均が０、共分散行列がＳｊ（ｋ）であるとし、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）は互いに無相関であるとする。２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）は、ともに目的信号（たとえば音声）に対応する源信号成分であっても良いし、一方が目的信号（たとえば音声）に対応する源信号成分で他方がノイズ（たとえばテレビからの雑音や雑踏における喧噪）に対応する源信号成分であっても良い。
【００４０】
観測信号ベクトルｙ（ｋ）から源信号成分ｘｊ（ｋ）を推定する方法としては、マルチチャネルウィーナーフィルタに基づく方法がある。この方法では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を（１０）式により計算する。
【００４１】
【数１０】
2019163487.tif 000012
【００４２】
すなわち、この方法では、観測信号ベクトルｙ（ｋ）に対し、Ｉ×Ｉの行列であるマルチチャネルウィーナーフィルタＳｊ（ｋ）（Ｓ１（ｋ）＋Ｓ２（ｋ））−１を掛けることにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する。
【００４３】
［従来の信号分析方法］
ここで、従来の信号分析方法（例えば、非特許文献１に記載の信号分析方法）では、（１０）式に基づいて源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する。その際、源信号成分ｘｊ（ｋ）の共分散行列Ｓｊ（ｋ）は、次の（１１）式によりモデル化される。なお、（１１）式におけるＥは、期待値演算である。
【００４４】
【数１１】
2019163487.tif 000013
【００４５】
ここで、ｖｊ（ｋ）は、源信号のサンプル点ごと（ｋごと）の分散をモデル化するパラメータである分散パラメータであり、Ｒｊは、源信号の空間的特性をモデル化するパラメータである空間共分散行列である。従来の信号分析方法では、（１０）式に（１１）式を代入して得られる（１２）式を用いて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算していた。
【００４６】
【数１２】
2019163487.tif 000014
【００４７】
上述の従来の信号分析方法の問題点として、膨大な計算量を要するため、大規模なデータセットの処理（サンプル点ｋの数が膨大になる）や、リアルタイム処理（短時間での処理が要求される）、計算能力の低い機器などへの適用が困難だった。
【００４８】
実際、従来技術では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する際に、逆行列計算（具体的には、（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１の計算）をサンプル点毎に（つまり多数回）行う必要があった。逆行列計算は、Ｉ次行列に対しＯ（Ｉ３）という大きな計算量を要する演算である。また、サンプル点の数は、典型的には数百〜数千万程にもなる。そのため、従来の信号分析方法では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する際に、膨大な計算量を必要としていた。
【００４９】
［第１の実施の形態に係る信号分析方法］
これに対し、第１の実施形態では、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）の同時無相関化に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を少ない計算量で計算する方法を提供する。以下、第１の実施形態における考え方について説明する。
【００５０】
各サンプル点ｋに対し信号α（ｋ）の共分散行列Ｅ（α（ｋ）α（ｋ）Ｈ）が対角行列であるとき、α（ｋ）は無相関であるという。これは、簡単に言えば、「ベクトルα（ｋ）の異なる要素同士は相関を持たない」ということである。いま仮に、（１１）式の源信号成分ｘｊ（ｋ）の共分散行列ｖｊ（ｋ）Ｒｊ（ｊは１または２）がいずれも対角行列であるとすると、（１２）式における逆行列演算（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１は単なる対角要素ごとの逆数演算により実現できるため、少ない計算量で計算できる。しかしながら、通常、源信号成分ｘｊ（ｋ）の異なる要素同士は互いに強い相関を持つため、その共分散行列ｖｊ（ｋ）Ｒｊは対角行列ではない。そこで、本第１の実施形態では、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）を同時無相関化することを考える。
【００５１】
２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）の同時無相関化とは、（１３）式及び（１４）式に示すように、ｘ１（ｋ），ｘ２（ｋ）に正則行列Ｖを掛ける線型変換を施して、線型変換後の源信号成分ｘ´１（ｋ），ｘ´２（ｋ）が同時に無相関になるようにする（すなわち、ｘ´１（ｋ），ｘ´２（ｋ）の共分散行列Ｓ´１（ｋ），Ｓ´２（ｋ）が同時に対角行列になるようにする）ことを指すものとする。この同時無相関化により、源信号成分ｘｊ（ｋ）の共分散行列ｖｊ（ｋ）Ｒｊはいずれも対角行列に変換される。したがって、同時無相関化により、（１２）式における逆行列演算（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１は、対角行列に対する逆行列演算に変わり、単なる対角要素ごとの逆数演算により実現できるようになるため、少ない計算量で計算できるようになる。
【００５２】
【数１３】
2019163487.tif 000015
【数１４】
2019163487.tif 000016
【００５３】
ｘ´１（ｋ），ｘ´２（ｋ）の共分散行列Ｓ´１（ｋ），Ｓ´２（ｋ）を計算すると、次の（１５）〜（１９）式のようになる。
【００５４】
【数１５】
2019163487.tif 000017
【数１６】
2019163487.tif 000018
【００５５】
同時無相関化の条件は、これらが同時に対角行列になること、すなわち以下の（２０）式及び（２１）式が同時に成り立つことである。なお、Λ１，Λ２は或る対角行列である。
【００５６】
【数１７】
2019163487.tif 000019
【数１８】
2019163487.tif 000020
【００５７】
すなわち、同時無相関化の目標は、空間共分散行列Ｒ１，Ｒ２に対し、左から正則行列Ｖを、右からＶのエルミート転置ＶＨを、それぞれ掛けたときに、これらが同時に対角行列になるような正則行列Ｖを求めること（同時対角化）である。
【００５８】
このような行列Ｖは、一般化固有値問題に基づいて、以下のようにして求めることができる。（２２）式に示す、空間共分散行列Ｒ１，Ｒ２に対する一般化固有値問題を解く。
【００５９】
【数１９】
2019163487.tif 000021
【００６０】
すると、Ｉ個の一般化固有値λ１，・・・，λＩと、λ１，・・・，λＩにそれぞれ対応するＩ個の一般化固有ベクトルｐ１，・・・，ｐＩと、が得られる。ここで、一般化固有ベクトルｐ１，・・・，ｐＩは、次の（２３）式を満たすように取ることができる。
【００６１】
【数２０】
2019163487.tif 000022
【００６２】
これらの一般化固有ベクトルｐ１，・・・，ｐＩを用いて、行列Ｐを次の（２４）式により定めると、これは、次の（２５）式及び（２６）式を満たす（例えば、参考文献１「G.  Strang,  “Linear  Algebra  and  Its  Applications”,  Harcourt  Brace  Jovanovich,  1988.」参照）。
【００６３】
【数２１】
2019163487.tif 000023
【００６４】
【数２２】
2019163487.tif 000024
【数２３】
2019163487.tif 000025
【００６５】
ただし、Ｉは単位行列、Λは一般化固有値λ１，・・・，λＩからなる（２７）式の対角行列である。
【００６６】
【数２４】
2019163487.tif 000026
【００６７】
したがって、行列Ｖを次の（２８）式により定めると、これは同時無相関化の条件である（２０）式及び（２１）式を満たす（Λ１=Λ、Λ２=Ｉである）。
【００６８】
【数２５】
2019163487.tif 000027
【００６９】
以下、（２７）式における一般化固有値からなる対角行列Λを一般化固有値行列と呼び、（２４）式における一般化固有ベクトルからなる行列Ｐを一般化固有ベクトル行列と呼ぶことがある。また、上述のように、Ｐは同時無相関化を実現する行列でもあるから、その意味でＰを同時無相関化行列とも呼ぶ。
【００７０】
信号パラメータ分析部２０は、空間共分散行列Ｒ１，Ｒ２に対する一般化固有値λ１，・・・，λＩと、λ１，・・・，λＩにそれぞれ対応する一般化固有ベクトルｐ１，・・・，ｐＩであって（２３）式を満たすものと、を計算し、一般化固有値λ１，・・・，λＩを対角要素とする対角行列である一般化固有値行列Λと、一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする正則行列である一般化固有ベクトル行列Ｐと、を、（２７）式と（２４）式とにより作成する。
【００７１】
上述のように、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）は、（２４）式で定義される行列Ｐのエルミート転置ＰＨを掛ける線型変換により、同時無相関化される。観測信号ベクトルのモデルを示す（９）式の両辺に対し、この線形変換を施すと、次の（２９）式を得る。
【００７２】
【数２６】
2019163487.tif 000028
【００７３】
ここで、（１８）式、（１９）式、（２５）式、（２６）式及び（２８）式より、ｘ´１（ｋ），ｘ´２（ｋ）の共分散行列Ｓ´１（ｋ），Ｓ´２（ｋ）は、次の（３０）式及び（３１）式のように、ともに対角行列である。
【００７４】
【数２７】
2019163487.tif 000029
【数２８】
2019163487.tif 000030
【００７５】
ただし、Λは（２７）式で定義される対角行列である。（２９）式にマルチチャネルウィーナーフィルタを適用する場合、ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）は、次の（３２）式で与えられる。
【００７６】
【数２９】
2019163487.tif 000031
【００７７】
（３２）式における共分散行列Ｓ´１（ｋ），Ｓ´２（ｋ）は、（３０）式及び（３１）式に示すように対角行列である。したがって、（３２）式は、以下の（３６）式及び（４０）式に示すように非常に簡単になる。このように、同時無相関化により、（１２）式における逆行列演算（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１は、（３３）式と（３７）式のように対角行列に対する逆行列演算（ｖ１（ｋ）Λ＋ｖ２（ｋ）Ｉ）−１に変わり、単なる対角要素ごとの逆数演算により実現できるようになるため、少ない計算量で計算できるようになる。
【００７８】
【数３０】
2019163487.tif 000032
【００７９】
【数３１】
2019163487.tif 000033
【００８０】
シングルチャネルウィーナーフィルタ部４０は、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を（３６）式及び（４０）式により計算する。
【００８１】
（３６）式及び（４０）式の処理は、シングルチャネルウィーナーフィルタに基づく処理と解釈することもできる。そのことを説明するために、まず、シングルチャネルウィーナーフィルタとは何であるかを説明する。マルチチャネルウィーナーフィルタが複数（Ｉ個）の観測信号を用いて推定を行うのに対し、シングルチャネルウィーナーフィルタは１個の観測信号のみを用いて推定を行う。シングルチャネルウィーナーフィルタでは、１個の観測信号（スカラー）ｙ（ｋ）を、２個の源信号成分（スカラー）ｘ１（ｋ），ｘ２（ｋ）の和として次の（４１）式でモデル化する。
【００８２】
【数３２】
2019163487.tif 000034
【００８３】
ここで、観測信号が１個しかないことに対応して、ｙ（ｋ），ｘ１（ｋ），ｘ２（ｋ）はいずれもスカラーである。源信号成分ｘｊ（ｋ）は、平均が０、分散がσ２ｊ（ｋ）であるとし、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）は互いに無相関であるとする。シングルチャネルウィーナーフィルタに基づく方法では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）は、次の（４２）式により計算される。
【００８４】
【数３３】
2019163487.tif 000035
【００８５】
すなわち、この方法では、観測信号ｙ（ｋ）に対し、シングルチャネルウィーナーフィルタσ２ｊ（ｋ）／（σ２１（ｋ）＋σ２２（ｋ））を掛けることにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する。
【００８６】
ここで、（２９）式を要素ごとに見た場合、第ｉ要素より次の（４３）式を得る。
【００８７】
【数３４】
2019163487.tif 000036
【００８８】
上式に基づくシングルチャネルウィーナーフィルタを用いて、［ｙ´（ｋ）］ｉから［ｘ´ｊ（ｋ）］ｉを推定することを考える。一般に、確率ベクトルαの共分散行列がＡであるとき、αの第ｌ要素［α］ｌの分散は、Ａの（ｌ，ｌ）要素［Ａ］ｌｌであることから、［ｘ´ｊ（ｋ）］ｉの分散は［Ｓ´ｊ（ｋ）］ｉｉである。よって、（３０）式及び（３１）式より、［ｘ´１（ｋ）］ｉの分散はｖ１（ｋ）λｉであり、［ｘ´２（ｋ）］ｉの分散はｖ２（ｋ）である。したがって、シングルチャネルウィーナーフィルタに基づく［ｘ´ｊ（ｋ）］ｉの推定値［＾ｘ´ｊ（ｋ）］ｉは、次の（４４）式及び（４５）式で与えられる。
【００８９】
【数３５】
2019163487.tif 000037
【数３６】
2019163487.tif 000038
【００９０】
これらはそれぞれ、（３６）式及び（４０）式の第ｉ要素に他ならない。
【００９１】
さて、上記のようにして得られた同時無相関化された源信号成分ｘ´ｊ（ｋ）＝ＰＨｘｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）は、推定が上手く行っていれば、次の（４６）式のように真値に近い値を取ると期待される。
【００９２】
【数３７】
2019163487.tif 000039
【００９３】
上式の両辺に行列（ＰＨ）−１を掛けると、次の（４７）式を得る。
【００９４】
【数３８】
2019163487.tif 000040
【００９５】
そこで、次の（４８）式により、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する。
【００９６】
【数３９】
2019163487.tif 000041
【００９７】
このように、第１の実施形態によれば、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）の計算において、サンプル点ごとの逆行列演算を行う必要がないため、少ない計算量で信号分離を実現することができる。なお、上式は逆行列（ＰＨ）−１を含むが、これはｋに依存しないため、サンプル点ごとではなく１回計算するだけで良いことに注意する。
【００９８】
従来の信号分析方法における源信号成分ｘｊ（ｋ）の推定値を＾ｘｊｃｏｎｖ（ｋ）で表し、第１の実施形態における源信号成分ｘｊ（ｋ）の推定値を＾ｘｊｐｒｏｐ（ｋ）で表す。＾ｘｊｃｏｎｖ（ｋ）と＾ｘｊｐｒｏｐ（ｋ）とは理論上一致することが示される。ここで、＾ｘｊｃｏｎｖ（ｋ）は、（４９）式によって与えられる。
【００９９】
【数４０】
2019163487.tif 000042
【０１００】
また、（５）式と、（３３）式〜（４０）式と、（４８）式と、により、＾ｘｊｐｒｏｐ（ｋ）は次の（５０）式及び（５１）式に等しい。
【０１０１】
【数４１】
2019163487.tif 000043
【数４２】
2019163487.tif 000044
【０１０２】
以下では、ｊ＝１の場合について説明する。なお、ｊ＝２の場合については、ｊ＝１の場合と同様であるから、説明を省略する。（２５）式及び（２６）式によって、以下の（５２）式及び（５３）式が導き出される。
【０１０３】
【数４３】
2019163487.tif 000045
【数４４】
2019163487.tif 000046
【０１０４】
これらを、（４９）式に代入すると、次の（５４）式〜（５７）式を得る。
【０１０５】
【数４５】
2019163487.tif 000047
【０１０６】
上述のことから、従来の信号分析方法と第１の実施形態とで、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）の計算結果は理論上一致する。ただし、もちろんコンピュータによる計算誤差が生じることはある。
【０１０７】
このように、本実施の形態１に係る信号分析装置１では、信号パラメータ分析部２０は、２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒ１，Ｒ２に対するｐｉＨＲ２ｐｍ＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする行列である同時無相関化行列Ｐまたはそのエルミート転置ＰＨを、Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る。信号分析装置１では、この無相関化するためのパラメータを用いて同時無相関化を行うことによって、フィルタリングでは、計算量の少ない要素ごとの掛け算を行うだけで足り、源信号成分の推定値の計算量を低減することができる。
【０１０８】
同時無相関化部３０は、信号パラメータ分析部２０が得た同時無相関化行列Ｐのエルミート転置ＰＨを、入力されたＩ個の観測信号による観測信号ベクトルに掛けることにより、２個の源信号に対応する成分が無相関化された同時無相関化済観測信号ベクトルを得る。
【０１０９】
そして、シングルチャネルウィーナーフィルタ部４０は、同時無相関化部３０が得た同時無相関化済観測信号ベクトルの各要素に、同時無相関化行列Ｐに基づく行列ＰＨＲｊＰ（ｊは１以上２以下の整数）の対角成分に基づくシングルチャネルウィーナーフィルタを適用して同時無相関化済源信号成分推定値によるベクトルを得る。
【０１１０】
さらに、同時無相関化逆変換部５０は、シングルチャネルウィーナーフィルタ部４０が得た同時無相関化済源信号成分推定値によるベクトルに、同時無相関化行列Ｐのエルミート転置ＰＨの逆行列（ＰＨ）−１を掛けることにより、源信号成分の推定値を得る。
【０１１１】
［第１の実施形態の効果］
このような構成を有する第１の実施形態によれば、上記のように、従来の信号分析方法と比べて、信号分離性能を低下させることなく、より少ない計算量で信号分離を実現することができる。すなわち、第１の実施形態では、観測信号ベクトルに対し、同時無相関化を行ってから、フィルタリングを行い、最後に同時無相関化の逆変換を行うことにより、源信号成分の推定値を計算する。ここで、同時無相関化を実現する行列Ｐは、空間共分散行列に対する一般化固有値問題を解くことにより、求めることができる。
【０１１２】
この結果、第１の実施形態では、従来技術と比べて少ない計算量で源信号成分の推定値を計算できる。これは主に、以下の理由による。
【０１１３】
第一に、第１の実施の形態では、同時無相関化により、フィルタリングは要素ごとの掛け算という計算量の小さい処理になり、サンプル点ごとの逆行列演算という計算量の大きい処理が不要になる。第二に、第１の実施の形態では、サンプル点によらない空間共分散行列に対する一般化固有値問題を解いて得られる行列Ｐを用いて同時無相関化を行うため、一般化固有値問題はサンプル点ごとではなく１回だけ解けば良い。これは、源信号成分の共分散行列が、サンプル点によらない空間共分散行列のスカラー倍であるとモデル化したことにより可能になったことである。第三に、第１の実施の形態では、同時無相関化の逆変換を行う際には、同時無相関化の行列ＰＨの逆行列演算を行う必要があるが、同じ理由により、これはサンプル点ごとではなく１回だけ行えば良い。
【０１１４】
さらに、第１の実施の形態では、観測信号ベクトルに対して、同時無相関化を行ってから、フィルタリングを行い、最後に同時無相関化の逆変換を行う、第１の実施形態における処理と、観測信号ベクトルに対して直接フィルタリングを行う従来技術における処理とは、等価である。このため、第１の実施形態によれば、従来技術と比べて信号分離性能を劣化させることなく、従来技術より少ない計算量で源信号成分の推定値を計算することができる。
【０１１５】
［第１の実施形態の変形例１］
次に、第１の実施形態の変形例１について説明する。第１の実施形態では、一般化固有ベクトルｐ１，・・・，ｐＩは、（２３）式の制約条件を満たすものとした。本第１の実施形態の変形例１では、一般化固有ベクトルｐ１，・・・，ｐＩが（２３）式以外の制約条件を満たすものとした場合の信号分析装置について説明する。本第１の実施形態の変形例１に係る信号分析装置における信号パラメータ分析部では、空間共分散行列Ｒ１，Ｒ２に対する（５８）式に示す一般化固有値問題を解く。
【０１１６】
【数４６】
2019163487.tif 000048
【０１１７】
この結果、Ｉ個の一般化固有値λ１，・・・，λＩとλ１，・・・，λＩにそれぞれ対応するＩ個の一般化固有ベクトルｐ１，・・・，ｐＩと、が得られる。ここで、上に述べた実施の形態１の例では、一般化固有ベクトルｐ１，・・・，ｐＩは、次の（５９）式を満たすように取る。
【０１１８】
【数４７】
2019163487.tif 000049
【０１１９】
これに対して、第１の実施形態の変形例１では、一般化固有ベクトルｐ１，・・・，ｐＩを、次の（６０）式を満たすように取る（各σｌは所定の正数）。
【０１２０】
【数４８】
2019163487.tif 000050
【０１２１】
これらの一般化固有ベクトルｐ１，・・・，ｐＩを用いて、行列Ｐを次の（６１）式により定めると、これは、（６２）式及び（６３）式を満たす。
【０１２２】
【数４９】
2019163487.tif 000051
【数５０】
2019163487.tif 000052
【数５１】
2019163487.tif 000053
【０１２３】
ただし、Λ１，Λ２は、次の（６４）式及び（６５）式に示す対角行列である。
【０１２４】
【数５２】
2019163487.tif 000054
【数５３】
2019163487.tif 000055
【０１２５】
したがって、行列Ｖを次の（６６）式により定めると、これは、同時無相関化の条件（２０）式及び（２１）式を満たす。ただし、（６４）式及び（６５）式のΛ１及びΛ２が、（２０）式及び（２１）式におけるΛ１及びΛ２に相当する。
【０１２６】
【数５４】
2019163487.tif 000056
【０１２７】
以下、（２７）式における一般化固有値からなる対角行列Λを一般化固有値行列、（６１）式における一般化固有ベクトルからなる行列Ｐを一般化固有ベクトル行列と呼ぶことがある。このようにして求めた行列Ｐを、同時無相関化行列として用いることができる。本第１の実施形態の変形例１に係る信号分析装置におけるシングルチャネルウィーナーフィルタ部では、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を、次の（６７）式及び（６８）式により計算する。
【０１２８】
【数５５】
2019163487.tif 000057
【数５６】
2019163487.tif 000058
【０１２９】
他の各部における処理は、上記の第１の実施形態と同様である。
【０１３０】
このように、第１の実施形態の変形例１では、信号パラメータ分析部２０は、２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒ１，Ｒ２に対するｐｉＨＲ２ｐｍ＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，・・・，ｐＩを列ベクトルとする行列である同時無相関化行列Ｐまたはそのエルミート転置ＰＨを、Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る。
【０１３１】
［第１の実施形態の変形例２］
次に、第１の実施形態の変形例２について説明する。第１の実施形態の変形例２では、信号パラメータ分析部における処理のもう一つの変形例について説明する。第１の実施形態では、信号パラメータ分析部２０は、固有値分析に基づいて、ＰＨＲ１ＰとＰＨＲ２Ｐとが対角行列になる行列Ｐとそのときの対角行列Λ＝ＰＨＲ１Ｐを求めていた。これに対し、本第１の実施形態の変形例２では、信号パラメータ分析部２０は、最適化に基づいて、Λ１＝ＰＨＲ１Ｐと、Λ２＝ＰＨＲ２Ｐとができるだけ、ある対角行列Σ１，Σ２に近づくような行列ＰとそのときのΛ１＝ＰＨＲ１Ｐ及びΛ２＝ＰＨＲ２Ｐを求める。
【０１３２】
すなわち、第１の実施形態の変形例２では、信号パラメータ分析部２０は、次のような行列ＰとそのときのΛ１＝ＰＨＲ１Ｐ及びΛ２＝ＰＨＲ２Ｐを求める（（６９）式及び（７０）式参照）。
【０１３３】
【数５７】
2019163487.tif 000059
【数５８】
2019163487.tif 000060
【０１３４】
具体的には、次の（７１）式のコスト関数を行列Ｐと対角行列Σ１とΣ２に関して最小化することにより、行列Ｐを求めることができる。
【０１３５】
【数５９】
2019163487.tif 000061
【０１３６】
ここで、||Ａ||Ｆは、行列Ａのフロベニウスノルムを表す。例えば、上式のコスト関数の最小化は、参考文献２「A.  Yeredor,  “Non-Orthogonal  Joint  Diagonalization  in  the  Least-Squares  Sense  with  Application  in  Blind  Source  Separation”,  IEEE  Transactions  on  Signal  Processing,  vol.  50,  no.  7,  pp.  1545−1553,  Jul.  2002.」に記載のアルゴリズムに従って行うことができる。このようにして求めた行列Ｐを、同時無相関化行列として用いることができる。
【０１３７】
また、Λ１＝ＰＨＲ１Ｐ及びΛ２＝ＰＨＲ２Ｐに基づき、第１の実施形態の変形例１と同様にして、シングルチャネルウィーナーフィルタ部４０における同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を計算することができる。他の各部における処理は、上に述べた第１の実施形態と同様である。
【０１３８】
このように、第１の実施形態の変形例２では、信号パラメータ分析部２０は、Ｊ個（この変形例ではＪは２）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列に近づくように求めた行列である同時無相関化行列Ｐまたは／及びそのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る。
【０１３９】
［第１の実施形態の変形例３］
次に、第１の実施形態の変形例３について説明する。本第１の実施形態の変形例３では、第１の実施形態の変形例２の方法に基づき、源信号が２個の場合に限らず、源信号が一般にＪ個の場合に拡張する場合について説明する。
【０１４０】
なお、第１の実施形態の変形例３においては、Ｊ個の源信号が混在する状況において、それぞれ異なる位置で取得されたＩ個の観測信号ｙｉ（ｋ）と、源信号のサンプル点ごとの分散（パワー）をモデル化するパラメータである分散パラメータｖｊ（ｋ）（ｊ＝１，・・・，Ｊは源信号のインデックス）と、源信号の空間的特性をモデル化するパラメータである空間共分散行列Ｒｊと、が信号分析装置１に入力されるものとする。
【０１４１】
第１の実施の形態の変形例３では、信号パラメータ分析部２０は、最適化に基づいて、Λ１＝ＰＨＲ１Ｐ，・・・，ΛＪ＝ＰＨＲＪＰが、できるだけ、ある対角行列Σ１，・・・，ΣＪに近づくような行列Ｐと、そのときのΛ１＝ＰＨＲ１Ｐ，・・・，ΛＪ＝ＰＨＲＪＰを求める。
【０１４２】
すなわち、第１の実施形態の変形例３では、信号パラメータ分析部２０は、次の（７２）式〜（７３）式に示すような行列Ｐと、そのときのΛ１＝ＰＨＲ１Ｐ，・・・，ΛＪ＝ＰＨＲＪＰを求める。
【０１４３】
【数６０】
2019163487.tif 000062
【０１４４】
具体的には、次の（７４）式のコスト関数を行列Ｐと各対角行列Σｊ（ｊは１以上Ｊ以下の整数）に関して最小化することにより、行列Ｐを求めることができる。
【０１４５】
【数６１】
2019163487.tif 000063
【０１４６】
例えば、上式のコスト関数の最小化は、参考文献２に記載のアルゴリズムに従って行うことができる。本第１の実施形態の変形例３では、このようにして求めたＰを、同時無相関化行列として用いることができる。また、Λ１＝ＰＨＲ１Ｐ，・・・，ΛＪ＝ＰＨＲＪＰに基づき、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を計算することができる。
【０１４７】
次に、源信号がＪ個の場合の信号分析装置１の構成と処理について説明する。信号分析装置１は、観測信号ベクトル作成部１０、信号パラメータ分析部２０、同時無相関化部３０、シングルチャネルウィーナーフィルタ部４０及び同時無相関化逆変換部５０を有する。観測信号ベクトル作成部１０は、入力された観測信号ｙｉ（ｋ）を取得し、すべての観測信号からなるＩ次元縦ベクトルである観測信号ベクトルｙ（ｋ）をサンプル点ごとに作成する。観測信号ベクトルｙ（ｋ）は、次の（７５）式で示される。
【０１４８】
【数６２】
2019163487.tif 000064
【０１４９】
信号パラメータ分析部２０は、入力された空間共分散行列Ｒｊを取得する。
【０１５０】
【数６３】
2019163487.tif 000065
【０１５１】
そして、上述の（７６）〜（７８）式で定義される各行列Λｊ（ｊは１以上Ｊ以下の整数）が対角行列になるような同時無相関化行列Ｐと、そのときのΛ１，・・・，ΛＪを計算する。
【０１５２】
同時無相関化部３０は、観測信号ベクトル作成部１０からの観測信号ベクトルｙ（ｋ）と、信号パラメータ分析部２０からの同時無相関化行列Ｐと、に基づいて、同時無相関化された観測信号ベクトルｙ´（ｋ）を計算する。具体的には、（７９）式に示すように、同時無相関化部３０は、観測信号ベクトルｙ（ｋ）に同時無相関化行列Ｐのエルミート転置ＰＨを掛けることにより、同時無相関化された観測信号ベクトルｙ´（ｋ）を計算する。
【０１５３】
【数６４】
2019163487.tif 000066
【０１５４】
なお、同時無相関化とは、観測信号ベクトルｙ（ｋ）に含まれる、Ｊ個の源信号に対応する成分（源信号成分）を同時に無相関化する線型変換を指す。
【０１５５】
シングルチャネルウィーナーフィルタ部４０は、入力された分散パラメータｖｊ（ｋ）を取得し、信号パラメータ分析部２０からの行列Λｊと、同時無相関化部３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、に基づいて、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を計算する。具体的には、シングルチャネルウィーナーフィルタ部４０は、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を次の（８０）式により計算する。
【０１５６】
【数６５】
2019163487.tif 000067
【０１５７】
同時無相関化逆変換部５０は、信号パラメータ分析部２０からの同時無相関化行列Ｐと、シングルチャネルウィーナーフィルタ部４０からの同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）と、に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算し、出力する。具体的には、同時無相関化逆変換部５０は、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）に対し、（８１）式のように、同時無相関化部３０で掛けた行列ＰＨの逆行列（ＰＨ）−１を掛けることにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算し、出力する。
【０１５８】
【数６６】
2019163487.tif 000068
【０１５９】
このように、第１の実施形態の変形例３では、信号パラメータ分析部２０は、Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列に近づくように求めた行列である同時無相関化行列Ｐまたは／及びそのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る。
【０１６０】
［第１の実施形態及び第１の実施形態の変形例１〜３の効果］
このように、本第１の実施形態及び本第１の実施形態の変形例１〜３によれば、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を少ない計算量で計算することができる。特に、第１の実施形態の固有値分析に基づく構成によれば、従来技術と比べて信号分離性能を劣化させることなく、従来技術より少ない計算量で信号分離を行うことができる。
【０１６１】
なお、第１の実施形態の変形例２，３では、対角行列Σｊに近づくように行列Λｊを求める。ただし、行列Λｊは実際には対角行列にならないことがある。この場合には、対角行列に近い行列Λｊ＝行列ＰＨＲｊＰ」の「対角成分」をλとして使ってシングルチャネルウィーナーフィルタの処理を行う。したがって、第１の実施形態の変形例２，３のシングルチャネルウィーナーフィルタの処理では、行列Λｊの非対角成分は用いず、対角成分のみを用いるので、非対角成分は特定しない構成としても良い。
【０１６２】
［第２の実施形態］
次に、第２の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
【０１６３】
第１の実施形態で説明したように、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが既知の場合には、これらに基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算できる。一方、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の場合には、これらを推定する必要がある。
【０１６４】
そこで、本第２の実施形態においては、２個の源信号が混在する状況において、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊと、源信号成分ｘｊ（ｋ）と、を反復法により同時に推定する方法について説明する。第２の実施形態においては、それぞれ異なる位置で取得されたＩ個の観測信号ｙｉ（ｋ）が信号分析装置に入力されるものとする。
【０１６５】
［信号分析装置の構成及び処理］
図３及び図４を用いて、第２の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図３は、第２の実施形態に係る信号分析装置の構成の一例を示す図である。図４は、第２の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０１６６】
図３に示すように、信号分析装置２０１は、初期化部（図示しない）、観測信号ベクトル作成部２１０、信号パラメータ分析部２２０、同時無相関化部２３０、シングルチャネルウィーナーフィルタ部２４０、事後共分散行列更新部２５０、分散パラメータ更新部２６０、空間共分散行列更新部２７０、収束判定部（図示しない）及び同時無相関化逆変換部２８０を有する。
【０１６７】
まず、信号分析装置２０１の各部の概要について説明する。図示しない初期化部は、分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）と、空間共分散行列Ｒｊの初期推定値＾Ｒｊ（０）と、を設定する（ステップＳ２０）。これらは、例えば乱数に基づいて設定すればよい。
【０１６８】
観測信号ベクトル作成部２１０は、入力された観測信号ｙｉ（ｋ）を取得し（ステップＳ２１）、すべての観測信号からなるＩ次元縦ベクトルである観測信号ベクトルｙ（ｋ）をサンプル点ごとに作成する（ステップＳ２２）。
【０１６９】
信号パラメータ分析部２２０は、空間共分散行列更新部２７０からの同時無相関化された空間共分散行列Ｒ´ｊの推定値＾Ｒ´ｊ（ただし、例外として、信号パラメータ分析部２２０における最初の処理の際には、入力された空間共分散行列Ｒｊの初期推定値＾Ｒｊ（０））に基づいて、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Ｐを更新する（ステップＳ２３）。
【０１７０】
同時無相関化部２３０は、観測信号ベクトル作成部２１０からの観測信号ベクトルｙ（ｋ）と、信号パラメータ分析部２２０からの一般化固有ベクトル行列Ｐと、に基づいて、同時無相関化された観測信号ベクトルｙ´（ｋ）を更新する（ステップＳ２４）。
【０１７１】
シングルチャネルウィーナーフィルタ部２４０は、信号パラメータ分析部２２０からの一般化固有値行列Λと、同時無相関化部２３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、分散パラメータ更新部２６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）（ただし、例外として、シングルチャネルウィーナーフィルタ部２４０における最初の処理の際には、入力された分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）と、に基づいて、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を更新する（ステップＳ２５）。
【０１７２】
事後共分散行列更新部２５０は、信号パラメータ分析部２２０からの一般化固有値行列Λと、分散パラメータ更新部２６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）（ただし、例外として、事後共分散行列更新部２５０における最初の処理の際には、入力された分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ））と、に基づいて、同時無相関化された事後共分散行列Σ´（ｋ）を更新する（ステップＳ２６）。
【０１７３】
分散パラメータ更新部２６０は、信号パラメータ分析部２２０からの一般化固有値行列Λと、シングルチャネルウィーナーフィルタ部２４０からの同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）と、事後共分散行列更新部２５０からの同時無相関化された事後共分散行列Σ´（ｋ）と、に基づいて、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を更新する（ステップＳ２７）。
【０１７４】
空間共分散行列更新部２７０は、シングルチャネルウィーナーフィルタ部２４０からの同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）と、事後共分散行列更新部２５０からの同時無相関化された事後共分散行列Σ´（ｋ）と、分散パラメータ更新部２６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、に基づいて、同時無相関化された空間共分散行列Ｒ´ｊの推定値＾Ｒ´ｊを更新する（ステップＳ２８）。
【０１７５】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ２９）。収束判定部が収束していないと判定した場合（ステップＳ２９：Ｎｏ）、ステップＳ２３の信号パラメータ分析部２２０での処理に戻って、処理が継続される。
【０１７６】
一方、収束判定部が収束したと判定した場合（ステップＳ２９：Ｙｅｓ）、同時無相関化逆変換部２８０での以下の処理が行われる。
【０１７７】
具体的には、同時無相関化逆変換部２８０は、信号パラメータ分析部２２０からの一般化固有ベクトル行列Ｐと、シングルチャネルウィーナーフィルタ部２４０からの同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）と、に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算し、出力する（ステップＳ３０）。
【０１７８】
次に、信号分析装置２０１の各部の詳細について説明する。なお、初期化部については上述の通りである。
【０１７９】
観測信号ベクトル作成部２１０は、入力された観測信号ｙｉ（ｋ）を取得し、すべての観測信号からなるＩ次元縦ベクトルである観測信号ベクトルｙ（ｋ）を次の（８２）式により作成する。
【０１８０】
【数６７】
2019163487.tif 000069
【０１８１】
信号パラメータ分析部２２０は、信号パラメータ分析部２２０における最初の処理の際に、初期化部からの空間共分散行列Ｒｊの初期推定値＾Ｒｊ（０）を受け取る。そして、信号パラメータ分析部２２０は、＾Ｒ１（０），＾Ｒ２（０）に対する一般化固有値問題を解くことにより、一般化固有値λ１，・・・，λＩと、λ１，・・・，λＩにそれぞれ対応する一般化固有ベクトルｐ１，・・・，ｐＩであって、（８３）式を満たすものを計算する。
【０１８２】
【数６８】
2019163487.tif 000070
【０１８３】
これによって、信号パラメータ分析部２２０は、一般化固有値行列Λ及び一般化固有ベクトル行列Ｐを次の（８４）式及び（８５）式により更新する。
【０１８４】
【数６９】
2019163487.tif 000071
【数７０】
2019163487.tif 000072
【０１８５】
信号パラメータ分析部２２０は、信号パラメータ分析部２２０における２回目以降の処理の際に、空間共分散行列更新部２７０からの同時無相関化された空間共分散行列Ｒ´ｊの推定値＾Ｒ´ｊを受け取って、＾Ｒ´１，＾Ｒ´２に対する一般化固有値問題を解くことにより、一般化固有値λ１，・・・，λＩと、λ１，・・・，λＩにそれぞれ対応する一般化固有ベクトルｑ１，・・・，ｑＩであって、（８６）式を満たすものを計算する。
【０１８６】
【数７１】
2019163487.tif 000073
【０１８７】
そして、信号パラメータ分析部２２０は、一般化固有値行列Λ及び行列Ｑを（８７）式及び（８８）式により更新する。
【０１８８】
【数７２】
2019163487.tif 000074
【数７３】
2019163487.tif 000075
【０１８９】
そして、信号パラメータ分析部２２０は、一般化固有ベクトル行列Ｐを次の（８９）式により更新する。
【０１９０】
【数７４】
2019163487.tif 000076
【０１９１】
同時無相関化部２３０は、（９０）式に示すように、観測信号ベクトルｙ（ｋ）に一般化固有ベクトル行列Ｐのエルミート転置ＰＨを掛けることにより、同時無相関化された観測信号ベクトルｙ´（ｋ）を更新する。
【０１９２】
【数７５】
2019163487.tif 000077
【０１９３】
シングルチャネルウィーナーフィルタ部２４０は、同時無相関化された源信号成分ｘ´ｊ（ｋ）の推定値＾ｘ´ｊ（ｋ）を（９１）式及び（９２）式により更新する。
【０１９４】
【数７６】
2019163487.tif 000078
【数７７】
2019163487.tif 000079
【０１９５】
ただし、例外として、シングルチャネルウィーナーフィルタ部２４０における最初の処理の際に、上の２式において、分散パラメータ更新部２６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）の代わりに、初期化部からの分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）を用いる。
【０１９６】
事後共分散行列更新部２５０は、同時無相関化された事後共分散行列Σ´（ｋ）を、次の（９３）式により更新する。
【０１９７】
【数７８】
2019163487.tif 000080
【０１９８】
ただし、例外として、事後共分散行列更新部２５０における最初の処理の際には、上式において、分散パラメータ更新部２６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）の代わりに、初期化部からの分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）を用いる。
【０１９９】
分散パラメータ更新部２６０は、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を次の（９４）式により更新する。
【０２００】
【数７９】
2019163487.tif 000081
【０２０１】
空間共分散行列更新部２７０は、同時無相関化された空間共分散行列Ｒ´ｊの推定値＾Ｒ´ｊを次の（９５）式により更新する。（９５）式において、Ｋはサンプル点の数である。
【０２０２】
【数８０】
2019163487.tif 000082
【０２０３】
同時無相関化逆変換部２８０は、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を次の（９６）式により計算し、出力する。
【０２０４】
【数８１】
2019163487.tif 000083
【０２０５】
［第２の実施形態に係る信号分析装置の処理の背景］
［従来の信号分析方法］
次に、信号分析装置２０１の処理の背景について説明する。まず、従来技術について説明する。２個の源信号が混在する状況において、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊと、源信号成分ｘｊ（ｋ）と、を反復法により同時に推定する方法が、非特許文献１に開示されている。
【０２０６】
ここで、従来技術における観測信号ベクトルｙ（ｋ）のモデル化について説明する。（９）式における源信号成分ｘｊ（ｋ）は、次の（９７）式のように平均０、共分散行列ｖｊ（ｋ）Ｒｊの複素ガウス分布に従うものとする。
【０２０７】
【数８２】
2019163487.tif 000084
【０２０８】
ここで、Ｎ（ｚ；μ，Σ）は、平均がμ、共分散行列がΣの複素ガウス分布（（９８）式）を表す。
【０２０９】
【数８３】
2019163487.tif 000085
【０２１０】
また、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）は互いに無相関であると仮定する。観測信号ベクトルｙ（ｋ）が与えられた条件下での源信号成分ｘｊ（ｋ）の事後確率密度は、次の（９９）式で与えられる（例えば、非特許文献１を参照）。
【０２１１】
【数８４】
2019163487.tif 000086
【０２１２】
ここで、平均＾ｘｊ（ｋ）及び共分散行列Σｊ（ｋ）は、次の（１００）式及び（１０１）式により与えられる。
【０２１３】
【数８５】
2019163487.tif 000087
【数８６】
2019163487.tif 000088
【０２１４】
（１００）式における平均＾ｘｊ（ｋ）は、（１０）式におけるマルチチャネルウィーナーフィルタに基づく源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）に他ならない（両者に同じ記号を用いているのはそのためである）。したがって、（１０）式におけるマルチチャネルウィーナーフィルタに基づく源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）は、観測信号ベクトルｙ（ｋ）が与えられた条件下での源信号成分ｘｊ（ｋ）の事後確率密度を最大化する点と解釈できることが分かった。（１０１）式はｊに依らないから、（１０２）式と書く。
【０２１５】
【数８７】
2019163487.tif 000089
【０２１６】
従来技術では、下記第１の処理〜第４の処理により、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊと、源信号成分ｘｊ（ｋ）と、を同時に推定する。
【０２１７】
第１の処理では、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、空間共分散行列Ｒｊの推定値＾Ｒｊと、を初期化する。これらは、例えば、乱数に基づいて初期化すればよい。
【０２１８】
第２の処理では、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、空間共分散行列Ｒｊの推定値＾Ｒｊと、に基づいて、源信号成分ｘｊ（ｋ）を推定する。第１の実施形態では、源信号成分ｘｊ（ｋ）の１つの推定値＾ｘｊ（ｋ）を計算したのに対し、ここでは源信号成分ｘｊ（ｋ）の事後確率密度（（９９）式参照）を更新する。具体的には、（９９）式の平均＾ｘｊ（ｋ）と共分散行列Σ（ｋ）（（１０２）式の通りΣ１（ｋ）及びΣ２（ｋ）に等しい。（１０１）式も参照。）を次の（１０３）式及び（１０４）式により更新する。
【０２１９】
【数８８】
2019163487.tif 000090
【数８９】
2019163487.tif 000091
【０２２０】
第３の処理では、源信号成分ｘｊ（ｋ）の事後確率密度の平均＾ｘｊ（ｋ）と共分散行列Σ（ｋ）に基づいて、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、空間共分散行列Ｒｊの推定値＾Ｒｊと、を次の（１０５）式及び（１０６）式により更新する。ここで、ｔｒはトレースである。
【０２２１】
【数９０】
2019163487.tif 000092
【数９１】
2019163487.tif 000093
【０２２２】
第４の処理では、第２及び第３の処理を、収束するまで交互に反復し、収束した場合に、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を出力する。
【０２２３】
従来技術の問題点として、膨大な計算量を要するため、サンプル点ｋの数が膨大になる大規模なデータセットの処理、短時間での処理が要求されるリアルタイム処理、或いは、計算能力の低い機器などへの適用が困難だった。
【０２２４】
実際、従来技術では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）と事後共分散行列Σ（ｋ）を計算する際に、逆行列計算（具体的には（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１の計算）と行列乗算（具体的には（ｖ１（ｋ）Ｒ１）（ｖ１（ｋ）Ｒ１＋ｖ２（ｋ）Ｒ２）−１（ｖ２（ｋ）Ｒ２）の計算）とを、サンプル点毎に（つまり多数回）行う必要があった。逆行列計算と行列乗算とはともに、Ｉ次行列に対し、Ｏ（Ｉ３）という大きな計算量を要する演算である。また、サンプル点の数は、典型的には数百〜数千万程にもなる。
【０２２５】
このため、従来の信号分析方法では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）と事後共分散行列Σ（ｋ）を計算する際に、膨大な計算量を必要としていた。さらに、従来では、信号分析装置全体として、上記の第４の処理における反復を行う必要がある（反復回数は典型的には数十回にもなる）ため、計算量はさらに膨大となる。
【０２２６】
［第２の実施の形態に係る信号分析方法］
これに対し、第２の実施形態では、２個の源信号成分ｘ１（ｋ），ｘ２（ｋ）の同時無相関化に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）と事後共分散行列Σ（ｋ）と空間共分散行列Ｒと分散パラメータｖとを少ない計算量で計算する方法を提供する。この結果、信号分析装置２０１全体としても、少ない計算量で信号分析を実現することができる。以下、第２の実施形態における考え方について説明する。
【０２２７】
従来技術（非特許文献１参照）の更新則（１０３）式〜（１０６）式において、何回目の反復であるかを明示すると、以下の（１０７）式〜（１１０）式のようになる。なお、ｎ＝１，・・・，Ｎは何回目の反復であるかを表すインデックス、Ｎは反復回数である。
【０２２８】
【数９２】
2019163487.tif 000094
【数９３】
2019163487.tif 000095
【数９４】
2019163487.tif 000096
【数９５】
2019163487.tif 000097
【０２２９】
空間共分散行列Ｒ１，Ｒ２の推定値＾Ｒ１（ｎ−１），＾Ｒ２（ｎ−１）の一般化固有値問題に対して、第１の実施形態と同様にして定義される一般化固有値行列及び一般化固有ベクトル行列を、Λ（ｎ−１）及びＰ（ｎ−１）で表す。これらは、（１１１）式及び（１１２）式を満たす。
【０２３０】
【数９６】
2019163487.tif 000098
【数９７】
2019163487.tif 000099
【０２３１】
すなわち、Λ（ｎ−１）及びＰ（ｎ−１）は（１１３）式及び（１１４）式を満たす。Ｐ（ｎ−１）は同時無相関化（同時対角化）を実現する行列であるから、その意味でＰ（ｎ−１）を同時無相関化行列とも呼ぶ。
【０２３２】
【数９８】
2019163487.tif 000100
【数９９】
2019163487.tif 000101
【０２３３】
（１０７）式からサンプル点ごとの逆行列演算を消すために、（１０７）式に（１１３）式と（１１４）式とを代入して整理すると、第１の実施形態と同様に、次の（１１５）式、（１１６）式を得る。
【０２３４】
【数１００】
2019163487.tif 000102
【０２３５】
同様に、（１０８）式からサンプル点ごとの逆行列演算を消すために、（１０８）式に、（１１３）式と（１１４）式とを代入すると、次の（１１８）式〜（１２１）式を得る。
【０２３６】
【数１０１】
2019163487.tif 000103
【０２３７】
これで、サンプル点ごとの逆行列演算を消すことができた。しかしながら、（１２１）式は、まだ、サンプル点ごとの行列の乗算を含んでいる。そこで、Σ（ｎ）（ｋ）の代わりに、（１２２）式により定義される、Σ´（ｎ）（ｋ）を更新することを考える。
【０２３８】
【数１０２】
2019163487.tif 000104
【０２３９】
（１２１）式より、Σ´（ｎ）（ｋ）は、（１２３）式で表せる。
【０２４０】
【数１０３】
2019163487.tif 000105
【０２４１】
（１２３）式には、サンプル点ごとの行列の乗算が含まれない。Σ´（ｎ）（ｋ）は同時無相関化の線型変換を施された事後共分散行列とみなせるから、以下、同時無相関化された事後共分散行列と呼ぶ。
【０２４２】
そして、（１１５）式と（１１６）式とにおいても、＾ｘｊ（ｎ）（ｋ）の代わりに、（１２４）式で定義される＾ｘ´ｊ（ｎ）（ｋ）を更新することを考える。
【０２４３】
【数１０４】
2019163487.tif 000106
【０２４４】
（１１５）式及び（１１６）式より、（１２５）式及び（１２６）式となる。
【０２４５】
【数１０５】
2019163487.tif 000107
【数１０６】
2019163487.tif 000108
＾ｘ´ｊ（ｎ）（ｋ）は同時無相関化された源信号成分である。
【０２４６】
次に、（１０９）式に着目すると、これは、＾ｘｊ（ｎ）（ｋ）とΣ（ｎ）（ｋ）に依存している。上述のように、第２の実施形態では、＾ｘｊ（ｎ）（ｋ）とΣ（ｎ）（ｋ）の代わりに＾ｘ´ｊ（ｎ）（ｋ）とΣ´（ｎ）（ｋ）とを更新する。そこで、（１０９）式も、＾ｘ´ｊ（ｎ）（ｋ）とΣ´（ｎ）（ｋ）とを用いて計算できる形に書き換える必要がある。（１２４）式及び（１２２）式より、＾ｘｊ（ｎ）（ｋ）とΣ（ｎ）（ｋ）は（１２７）式及び（１２８）式と表せる。
【０２４７】
【数１０７】
2019163487.tif 000109
【数１０８】
2019163487.tif 000110
【０２４８】
これらを（１０９）式に代入すると（１２９）式〜（１３３）式を経て（１３４）式が得られる。
【０２４９】
【数１０９】
2019163487.tif 000111
【０２５０】
（１１０）式についても同様に、＾ｘ´ｊ（ｎ）（ｋ）とΣ´（ｎ）（ｋ）とを用いて計算できる形に書き換えるために、（１２７）式と（１２８）式とを代入すると、次の（１３５）式及び（１３６）式のようになる。
【０２５１】
【数１１０】
2019163487.tif 000112
【０２５２】
ここでも、＾Ｒｊ（ｎ）の代わりに、（１３７）式で定義される＾Ｒ´ｊ（ｎ）を更新することにする。
【０２５３】
【数１１１】
2019163487.tif 000113
【０２５４】
これにより、逆行列演算と行列の乗算とをいくつか減らせるため、計算量を削減することができる。また、＾Ｒ´ｊ（ｎ）は、（１３６）式より、（１３８）式のように表せる。
【０２５５】
【数１１２】
2019163487.tif 000114
【０２５６】
＾Ｒ´ｊ（ｎ）は同時無相関化の線型変換を施された空間共分散行列とみなせるから、以下、同時無相関化された空間共分散行列と呼ぶ。
【０２５７】
最後に、一般化固有値行列Λ及び一般化固有ベクトル行列Ｐの更新について説明する。最初の反復においては、空間共分散行列Ｒ１，Ｒ２の初期推定値＾Ｒ１（０），＾Ｒ２（０）の一般化固有値問題を解くことにより、Λ（０）及びＰ（０）を求めればよい。
【０２５８】
２回目以降の反復においては、同時無相関化された空間共分散行列の推定値＾Ｒ´１（ｎ），＾Ｒ´２（ｎ）が得られている。これを同時無相関化の逆変換により＾Ｒ１（ｎ），＾Ｒ２（ｎ）に戻してから、一般化固有値問題を適用することも可能だが、やや無駄な計算が生じる。
【０２５９】
そこで、以下のようにすれば、より効率的に計算できる。まず、＾Ｒ´１（ｎ），＾Ｒ´２（ｎ）の一般化固有値問題を解いて、一般化固有値を対角成分とする対角行列Λ（ｎ）と、一般化固有ベクトルを列ベクトルとする行列Ｑ（ｎ）と、を求める。これらは、（Ｑ（ｎ））Ｈ（＾Ｒ´１（ｎ））Ｑ（ｎ）＝Λ（ｎ）及び（Ｑ（ｎ））Ｈ（＾Ｒ´２（ｎ））Ｑ（ｎ）＝Ｉを満たすものとする。
【０２６０】
次に、Ｐ（ｎ）を（１３９）式により更新する。
【０２６１】
【数１１３】
2019163487.tif 000115
【０２６２】
このとき、Λ（ｎ）及びＰ（ｎ）は、＾Ｒ１（ｎ），＾Ｒ２（ｎ）に対する一般化固有値行列及び一般化固有ベクトル行列になる。これは、次の（１４０）式〜（１４７）式のようにして確認できる。
【０２６３】
【数１１４】
2019163487.tif 000116
【数１１５】
2019163487.tif 000117
【０２６４】
［第２の実施形態の効果］
このように、第２の実施形態によれば、（同時無相関化された）源信号成分の推定値及び（同時無相関化された）事後共分散行列と空間共分散行列と分散パラメータとを従来技術と比べて少ない計算量で計算できるため、従来技術と比べて少ない計算量で信号分離を実現できる。
【０２６５】
これは以下の理由による。第一に、第１の実施形態と同様に、同時無相関化に基づくことで、サンプル点ごとの逆行列演算が不要になる。第二に、事後共分散行列Σ（ｎ）（ｋ）そのものではなく、同時無相関化された事後共分散行列Σ´（ｎ）（ｋ）を更新するようにすることで、サンプル点ごとの行列の乗算が不要になる。
【０２６６】
また、上記のように、従来の信号分析方法の処理と第２の実施形態の処理とは理論上等価であるため、従来の信号分析方法と第２の実施形態とで信号分析結果は理論上一致する（ただし、もちろんコンピュータによる計算誤差が生じることはある）。したがって、第２の実施形態によれば、従来の信号分析方法と比べて、信号分析性能を低下させることなく、より少ない計算量で信号分析を実現することができる。
【０２６７】
また、第１の実施形態とは異なり、第２の実施形態によれば、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の場合でも、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊと、源信号成分ｘｊ（ｋ）と、を同時に推定することができる。
【０２６８】
［第３の実施形態］
次に、第３の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
【０２６９】
第３の実施形態においては、２個の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたＩ個（Ｉは２以上の整数）の観測信号ｙｉ（ｔ）が信号分析装置に入力されるものとする。観測信号ｙｉ（ｔ）は、周波数分解され、周波数ごとに第２の実施形態に記載の方法で信号分離がなされ、最後に時間領域に変換される。
【０２７０】
［信号分析装置の構成及び処理］
図５及び図６を用いて、第３の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図５は、第３の実施形態に係る信号分析装置の構成の一例を示す図である。図６は、第３の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０２７１】
図５に示すように、信号分析装置３０１は、周波数領域変換部３１０、周波数領域信号分離部３２０及び時間領域変換部３３０を有する。
【０２７２】
まず、信号分析装置３０１の各部の概要について説明する。周波数領域変換部３１０は、入力された観測信号ｙｉ（ｔ）を受け取り（ステップＳ４０）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙｉ（ｋ，ｆ）を計算する（ステップＳ４１）。ここで、ｋは時間フレームのインデックスであり、ｆは周波数ビンのインデックスである。以下、ｙｉ（ｋ，ｆ）をｙｉｆ（ｋ）で表す。
【０２７３】
周波数領域信号分離部３２０は、ｆごとに第２の実施形態の方法を適用して、周波数領域での信号分離を行う（ステップＳ４２）。すなわち、ｆ＝１に対しては、ｙｉ１（ｋ）を入力として、第２の実施形態の方法により源信号成分の推定値＾ｘ１１（ｋ），＾ｘ２１（ｋ）を得る。
【０２７４】
そして、ｆ＝２に対しては、ｙｉ２（ｋ）を入力として、第２の実施形態の方法により、源信号成分の推定値＾ｘ１２（ｋ），＾ｘ２２（ｋ）を得る。各ｆについてこれらを繰り返す。
【０２７５】
そして、ｆ＝Ｆに対しては、ｙｉＦ（ｋ）を入力として、第２の実施形態の方法により、源信号成分の推定値＾ｘ１Ｆ（ｋ），＾ｘ２Ｆ（ｋ）を得る。以下、＾ｘｊｆ（ｋ）を、＾ｘｊ（ｋ，ｆ）で表す。
【０２７６】
時間領域変換部３３０は、時間周波数領域の源信号成分の推定値＾ｘｊ（ｋ，ｆ）に基づいて、逆短時間フーリエ変換などにより、時間領域の源信号成分の推定値＾ｘｊ（ｔ）を計算し、出力する（ステップＳ４３）。
【０２７７】
このように、第３の実施形態では、時間周波数領域において、周波数ビンｆごとに第２の実施形態に係る方法を適用する。これに対し、非特許文献１には、時間周波数領域において、周波数ビンｆごとに第２の実施形態に記載した従来技術を適用する方法が開示されている。
【０２７８】
この従来技術の問題点として、膨大な計算量を要するため、大規模なデータセットの処理や、リアルタイム処理（短時間での処理が要求される）、計算能力の低い機器などへの適用が困難だった。実際、従来技術では、逆行列計算と行列乗算とを、周波数ビンごと、フレームごと、反復ごとに（つまり多数回）行う必要があった。このため、従来の信号分析方法では、膨大な計算量を必要としていた。
【０２７９】
例えば、後述の実験では、８秒の観測信号に対して、周波数ビンの数は２５６、フレームの数は２４９、反復数は１０であり、処理に７５０秒を要した。これが８００００秒（およそ２２時間）のデータセットになれば、処理に７５０００００秒（およそ８７日）掛かる計算になる。
【０２８０】
［第３の実施形態の効果］
本発明の効果を確認するために、従来の方法及び第３の実施形態を用いた音源分離実験について説明する。
【０２８１】
図７は、従来の方法及び第３の実施形態を用いた音源分離実験の実験条件を示す図である。図８は、従来の方法及び第３の実施形態を用いた音源分離実験結果（real  time  factor）を示す図である。図９は、従来の方法及び第３の実施形態を用いた音源分離実験結果（音源分離性能）を示す図である。
【０２８２】
８秒間の英語音声に、図７に示す実験室で測定された残響インパルス応答を畳み込んで残響音声を生成し、複数人の残響音声を加算することで観測信号を生成した。このとき、観測信号に対して、各方法を用いて音源分離を行った場合の性能を図８及び図９に示す。図８は、real  time  factor（観測信号長を１としたときの、処理に要した時間）である。また、図９は、音源分離性能である信号対歪み比（値が大きいほど音源分離性能が高い）である。
【０２８３】
図８及び図９に示すように、本第３の実施形態に係る方法によれば、従来の方法と比べて音源分離性能を低下させることなく、２５０倍以上の計算時間削減が実現できたことが分かる。なお、残響時間が１３０ミリ秒及び４４０ミリ秒の場合には、本第３の実施形態に係る方法による信号対歪み比は、従来の方法と比べてわずかに低下しているが、これはコンピュータによる計算の誤差によるものである。このように、本第３の実施形態に係る方法は従来の方法と理論的には等価だが、実際には計算誤差に起因して性能差が生じることもあるものの、それは図８のようにごく僅かな差に過ぎない。
【０２８４】
このように、第３の実施形態によれば、従来の方法と比べて音源分離性能を低下させることなく、計算時間の大きな削減が実現できる。
【０２８５】
［第４の実施形態］
次に、第４の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。
【０２８６】
第４の実施形態は、第２の実施形態において、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊを、補助関数法を用いて推定する構成に変更したものである。
【０２８７】
［信号分析装置の構成及び処理］
図１０及び図１１を用いて、第４の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図１０は、第４の実施形態に係る信号分析装置の構成の一例を示す図である。図１１は、第４の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０２８８】
図１０に示すように、信号分析装置４０１は、初期化部（図示しない）、観測信号ベクトル作成部４１０、信号パラメータ分析部４２０、同時無相関化部４３０、分散パラメータ更新部４４０、空間共分散行列更新部４５０、収束判定部（図示しない）を有する。
【０２８９】
まず、信号分析装置４０１の各部の概要について説明する。図示しない初期化部は、分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）と、空間共分散行列Ｒｊの初期推定値＾Ｒｊ（０）と、を設定する（ステップＳ５１）。これらは、例えば乱数に基づいて設定すればよい。
【０２９０】
観測信号ベクトル作成部４１０は、入力された観測信号ｙｉ（ｋ）を取得し（ステップＳ５２）、すべての観測信号からなるＩ次元縦ベクトルである観測信号ベクトルｙ（ｋ）をサンプル点ごとに作成する（ステップＳ５３）。
【０２９１】
信号パラメータ分析部４２０は、空間共分散行列更新部４５０からの空間共分散行列Ｒｊの推定値＾Ｒｊに基づいて、一般化固有値問題を解くことにより、一般化固有値行列Λ及び一般化固有ベクトル行列Ｐを更新する(ステップＳ５４)。具体的には、信号パラメータ分析部４２０は、空間共分散行列Ｒ１，Ｒ２の推定値＾Ｒ１，＾Ｒ２の一般化固有値問題に対して、第１の実施形態と同様にして定義される一般化固有値行列Λ及び一般化固有ベクトル行列Ｐを求める。これらは、（１４８）式及び（１４９）式を満たすものとする。
【０２９２】
【数１１６】
2019163487.tif 000118
【数１１７】
2019163487.tif 000119
【０２９３】
ただし、例外として、信号パラメータ分析部４２０における最初の処理の際には、信号パラメータ分析部４２０における上述の処理において、空間共分散行列更新部４５０からの空間共分散行列Ｒｊの推定値＾Ｒｊの代わりに、入力された空間共分散行列Ｒｊの初期推定値^Ｒｊ（０）を用いる。
【０２９４】
同時無相関化部４３０は、観測信号ベクトル作成部４１０からの観測信号ベクトルｙ（ｋ）と、信号パラメータ分析部４２０からの一般化固有ベクトル行列Ｐと、に基づいて、同時無相関化された観測信号ベクトルｙ´（ｋ）を（９０）式により更新する（ステップＳ５５）。
【０２９５】
分散パラメータ更新部４４０は、同時無相関化部４３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、信号パラメータ分析部４２０からの一般化固有値行列Λと、に基づいて、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を次の（１５０）式及び（１５１）式により更新する（ステップＳ５６）。
【０２９６】
【数１１８】
2019163487.tif 000120
【数１１９】
2019163487.tif 000121
【０２９７】
ここで、右辺の＾ｖ１（ｋ）及び＾ｖ２（ｋ）は、分散パラメータ更新部４４０における前回の処理で得られた分散パラメータの推定値（ただし、例外として、分散パラメータ更新部４４０における最初の処理の際には、初期化部からの分散パラメータの初期推定値）を表す。
【０２９８】
空間共分散行列更新部４５０は、同時無相関化部４３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、信号パラメータ分析部４２０からの一般化固有ベクトル行列Ｐ及び一般化固有値行列Λと、分散パラメータ更新部４４０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、に基づいて、空間共分散行列Ｒｊの推定値＾Ｒｊを更新する（ステップＳ５７）。具体的には、空間共分散行列更新部４５０は、まず次の（１５２）式により行列Ｃｊを更新する。
【０２９９】
【数１２０】
2019163487.tif 000122
【０３００】
次に、空間共分散行列更新部４５０は、次の（１５３）式によりベクトルｚ（ｋ）を更新する。
【０３０１】
【数１２１】
2019163487.tif 000123
【０３０２】
次に、空間共分散行列更新部４５０は、次の（１５４）式により行列Ｄｊを更新する。
【０３０３】
【数１２２】
2019163487.tif 000124
【０３０４】
次に、空間共分散行列更新部４５０は、次の（１５５）式により空間共分散行列Ｒｊの推定値＾Ｒｊを更新する。
【０３０５】
【数１２３】
2019163487.tif 000125
【０３０６】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ５８）。収束判定部が収束していないと判定した場合（ステップＳ５８：Ｎｏ）、信号パラメータ分析部４２０での処理（ステップＳ５４）に戻って、処理が継続される。
【０３０７】
一方、収束判定部が収束したと判定した場合（ステップＳ５８：Ｙｅｓ）、分散パラメータ更新部４４０が分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を出力し、空間共分散行列更新部４５０が空間共分散行列Ｒｊの推定値＾Ｒｊを出力する（ステップＳ５９）。
【０３０８】
[第４の実施形態の背景]
第４の実施形態の背景について説明する。補助関数法に基づけば、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊは、次のように推定できる。分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値はあらかじめ初期化され、以下に述べる反復処理により更新される。分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）は、次の（１５６）式により更新される。
【０３０９】
【数１２４】
2019163487.tif 000126
【０３１０】
空間共分散行列Ｒｊの更新においては、まず（１５７）式及び（１５８）式により行列Ｃｊ及び行列Ｄｊが更新され、次に、（１５９）式により空間共分散行列Ｒｊの推定値＾Ｒｊが更新される。
【０３１１】
【数１２５】
2019163487.tif 000127
【数１２６】
2019163487.tif 000128
【数１２７】
2019163487.tif 000129
【０３１２】
ここで、半正定値エルミート行列Ａの平方根Ａ１／２は、（１６０）式をＡの固有値分解とすると、（１６１）式により定義される。ここで、ＵはＡの固有ベクトルからなるユニタリ行列、ΛはＡの固有値からなる対角行列、Σ１／２はΣの各対角要素[Σ]ｉｉをその平方根（[Σ]ｉｉ）１／２で置き換えることにより得られる対角行列である。
【０３１３】
【数１２８】
2019163487.tif 000130
【数１２９】
2019163487.tif 000131
【０３１４】
上記の分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値＾Ｒｊの更新は、収束するまで反復される。
【０３１５】
上記の分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値＾Ｒｊの更新則は、参考文献３「Kazuyoshi  Yoshii，Ryota  Tomioka，Daichi  Mochihashi，and  Masataka  Goto，“Infnite  Positive  Semidefinite  Tensor  Factorization  for  Source  Separation  of  Mixture  Signals”，Proceedings  of  the  30th  International  Conference  on  Machine  Learning，2013．」、参考文献４「Hiroshi  Sawada，Hirokazu  Kameoka，Shoko  Araki，and  Naonori  Ueda，“Multichannel  Extensions  of  Non−Negative  Matrix  Factorization  With  Complex−Valued  Data”，IEEE  TRANSACTIONS  ON  AUDIO，  SPEECH，  AND  LANGUAGE  PROCESSING，VOL．21，No.5，MAY  2013．」に開示されている方法と同様の方法により導出することができる。
【０３１６】
上記の方法の問題点として、（１５６）式〜（１５８）式の各処理において、サンプル点ごとの逆行列演算を行う必要があるため、膨大な計算量を必要とする。
【０３１７】
これに対し、本第４の実施形態では、同時無相関化（同時対角化）に基づいて、サンプル点ごとの逆行列演算を回避することにより、計算量を削減することができる。実際、（１４８）式及び（１４９）式を＾Ｒ１、＾Ｒ２について解いたものを（１５６）式〜（１５８）式に代入することにより、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値の更新則が以下のように得られる。
【０３１８】
分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）の更新則は、次の（１６２）式及び（１６３）式のようになる。
【０３１９】
【数１３０】
2019163487.tif 000132
【０３２０】
ここで、ｙ´（ｋ）は、（１６４）式の同時無相関化済観測信号ベクトルであり、Λ１は一般化固有値行列Λであり、Λ２は単位行列Ｉである。
【０３２１】
【数１３１】
2019163487.tif 000133
【０３２２】
空間共分散行列Ｒｊの推定値＾Ｒｊの更新則は次の（１６５）式及び（１６６）式のようになる。まず、行列Ｃｊに関しては、次式のようになる。
【０３２３】
【数１３２】
2019163487.tif 000134
【０３２４】
行列Ｄｊに関しては次の（１６７）式及び（１６８）式のようになる。
【０３２５】
【数１３３】
2019163487.tif 000135
【０３２６】
ただし、ベクトルｚ（ｋ）は、次の（１６９）式〜（１７１）式により定義される。
【０３２７】
【数１３４】
2019163487.tif 000136
【０３２８】
［第４の実施形態の効果］
このように、同時無相関化（同時対角化）により、サンプル点ごとの逆行列演算を不要にし、計算量を削減することができる。具体的には、同時無相関化（同時対角化）に基づき、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）の更新式を（１６２）式のように対角行列Λｊ＝ＰＨＲｊＰを用いて表すことができ、その結果、（１６３）式のように逆行列演算を不要にすることができる。空間共分散行列Ｒｊの推定値＾Ｒｊの更新則についても同様である。
【０３２９】
なお、（１５５）式の処理では、（１５９）式の処理と同様、逆行列演算を行う必要があるが、これはサンプル点ごとではなく各反復にて１回だけ行えば良いため、その計算量は比較的小さい。同様に、（１５５）式の処理では、（１５９）式の処理と同様、行列の平方根の計算のために固有値分解を行う必要があるが、これはサンプル点ごとではなく各反復にて２回だけ行えば良いため、その計算量は比較的小さい。また、（１６６）式、（１６８）式、及び（１５５）式の処理では、逆行列演算と同様に計算量の大きい演算である行列の乗算を含んでいるが、これもサンプル点ごとに行う必要はないため、その計算量は比較的小さい。
【０３３０】
［第４の実施形態の変形例］
なお、本実施形態においては、源信号が２個の場合について説明したが、第１の実施形態の変形例３と同様の考え方に基づけば、源信号が２個の場合に限らず、源信号が一般にＪ個の場合に拡張することができる。
【０３３１】
上記のように、第１〜第４の実施形態では、信号パラメータ分析部が、Ｊ個（Ｊは２以上の整数）の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列Ｒｊ（ｊは１以上Ｊ以下の整数）について、ＰＨＲｊＰがすべて対角行列になるような行列である同時無相関化行列Ｐまたは／及び、そのエルミート転置ＰＨを、それぞれ異なるＩ個（Ｉは２以上の整数）の位置で取得される観測信号による観測信号ベクトルについてＪ個の源信号に対応する成分を無相関化するためのパラメータとして得る。信号分析装置１，２０１，３０１，４０１では、この無相関化するためのパラメータを用いて同時無相関化を行うことによって、フィルタリングでは、計算量の少ない要素ごとの掛け算を行うだけで足り、源信号成分の推定値の計算量を低減することができる。
【０３３２】
［第５の実施形態］
次に、第５の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。第５の実施形態は、第１の実施形態の変形例３のバリエーションである。第１の実施形態の変形例３では、分散パラメータｖｊ（ｋ）と空間共分散行列Ｒｊが信号分析装置に入力されるとの前提条件がある。これに対し、第５の実施形態では、分散パラメータｖｊ（ｋ）と空間共分散行列Ｒｊとが信号分析装置に入力される必要はなく、分散パラメータｖｊ（ｋ）と、空間共分散行列Ｒｊをモデル化するパラメータである同時無相関化行列Ｐ及び同時無相関化された空間共分散行列Λｊと、源信号成分ｘｊ（ｋ）との同時推定を実現する。
【０３３３】
本第５の実施形態では、一般にＪ個の源信号が混在し、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の状況で、分散パラメータｖｊ（ｋ）と、空間共分散行列Ｒｊをモデル化するパラメータである同時無相関化行列Ｐ及び同時無相関化された空間共分散行列Λｊと、源信号成分ｘｊ（ｋ）と、を同時に推定することにより、信号分離を実現する方法について説明する。本第５の実施形態では、それぞれ異なる位置で取得されたＩ個の観測信号ｙｉ（ｋ）が信号分析装置に入力されるものとする。
【０３３４】
［信号分析装置の構成及び処理］
図１２及び図１３を用いて、第５の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図１２は、第５の実施形態に係る信号分析装置の構成の一例を示す図である。図１３は、第５の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０３３５】
図１２に示すように、信号分析装置５０１は、初期化部（図示しない）、観測信号ベクトル作成部５１０、同時無相関化部５３０、シングルチャネルウィーナーフィルタ部５４０、事後共分散行列更新部５５０、分散パラメータ更新部５６０、空間共分散行列更新部５７０、信号パラメータ分析部５２０、収束判定部（図示しない）及び同時無相関化逆変換部５８０を有する。
【０３３６】
まず、信号分析装置５０１の各部の概要について説明する。図示しない初期化部は、分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ）と、同時無相関化行列Ｐの初期推定値＾Ｐ（０）と、同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０）と、を設定する（ステップＳ６１）。これらは、例えば乱数に基づいて設定すればよい。
【０３３７】
観測信号ベクトル作成部５１０は、入力されたそれぞれ異なる位置で取得されたＩ個の観測信号ｙｉ（ｋ）を取得し（ステップＳ６２）、Ｉ次元縦ベクトルである観測信号ベクトルｙ（ｋ）をサンプル点ごとに（８２）式により作成する（ステップＳ６３）。
【０３３８】
同時無相関化部５３０は、観測信号ベクトル作成部５１０からの観測信号ベクトルｙ（ｋ）と、信号パラメータ分析部５２０からの同時無相関化行列Ｐの推定値＾Ｐ（ただし例外として、同時無相関化部５３０における最初の処理の際には、初期化部からの同時無相関化行列Ｐの初期推定値＾Ｐ（０））と、に基づいて、同時無相関化された観測信号ベクトルｙ´（ｋ）を次の（１７２）式により更新する（ステップＳ６４）。
【０３３９】
【数１３５】
2019163487.tif 000137
【０３４０】
ただし、同時無相関化部５３０における最初の処理の際には、（１７２）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。
【０３４１】
シングルチャネルウィーナーフィルタ部５４０は、同時無相関化部５３０からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、分散パラメータ更新部５６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）（ただし、例外として、シングルチャネルウィーナーフィルタ部５４０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ））と、空間共分散行列更新部５７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊ（ただし、例外として、シングルチャネルウィーナーフィルタ部５４０における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０））と、に基づいて、同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）を次の（１７３）式により更新する（ステップＳ６５）。
【０３４２】
【数１３６】
2019163487.tif 000138
【０３４３】
ただし、シングルチャネルウィーナーフィルタ部５４０における最初の処理の際には、（１７３）式の右辺において、＾Λｊ（ｊ＝１，・・・，Ｊ）を＾Λｊ（０）（ｊ＝１，・・・，Ｊ）に、＾ｖｊ（ｋ）（ｊ＝１，・・・，Ｊ；ｋ＝１，・・・，Ｋ）を＾ｖｊ（０）（ｋ）（ｊ＝１，・・・，Ｊ；ｋ＝１，・・・，Ｋ）に、それぞれ置き換える。
【０３４４】
事後共分散行列更新部５５０は、分散パラメータ更新部５６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）（ただし、例外として、事後共分散行列更新部５５０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ））と、空間共分散行列更新部５７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊ（ただし、例外として、事後共分散行列更新部５５０における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０））と、に基づいて、同時無相関化された事後共分散行列Σｊ´（ｋ）を次の（１７４）式により更新する（ステップＳ６６）。
【０３４５】
【数１３７】
2019163487.tif 000139
【０３４６】
ただし、事後共分散行列更新部５５０における最初の処理の際には、（１７４）式の右辺において、＾Λｊを＾Λｊ（０）に、＾ｖｊ（ｋ）を＾ｖｊ（０）（ｋ）に、それぞれ置き換える。
【０３４７】
分散パラメータ更新部５６０は、シングルチャネルウィーナーフィルタ部５４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、事後共分散行列更新部５５０からの同時無相関化された事後共分散行列Σｊ´（ｋ）と、空間共分散行列更新部５７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊ（ただし、例外として、分散パラメータ更新部５６０における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０））と、に基づいて、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を次の（１７５）式により更新する（ステップＳ６７）。
【０３４８】
【数１３８】
2019163487.tif 000140
【０３４９】
ただし、分散パラメータ更新部５６０における最初の処理の際には、（１７５）式の右辺において、＾Λｊを＾Λｊ（０）に置き換える。
【０３５０】
空間共分散行列更新部５７０は、シングルチャネルウィーナーフィルタ部５４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、事後共分散行列更新部５５０からの同時無相関化された事後共分散行列Σｊ´（ｋ）と、分散パラメータ更新部５６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、に基づいて、同時無相関化された空間共分散行列Λｊの推定値＾Λｊを次の（１７６）式により更新する（ステップＳ６８）。
【０３５１】
【数１３９】
2019163487.tif 000141
【０３５２】
信号パラメータ分析部５２０は、シングルチャネルウィーナーフィルタ部５４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、事後共分散行列更新部５５０からの同時無相関化された事後共分散行列Σｊ´（ｋ）と、分散パラメータ更新部５６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、空間共分散行列更新部５７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊと、（信号パラメータ分析部５２０における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、更に、初期化部からの同時無相関化行列Ｐの初期推定値＾Ｐ（０）と、）を受け取って、次の（１７７）式により同時無相関化行列Ｐの推定値＾Ｐを更新する（ステップＳ６９）。ただし、ｅｉは単位行列の第ｉ列を表し、行列Ａに対して［Ａ］ｉは、Ａの第ｉ列を表す。
【０３５３】
【数１４０】
2019163487.tif 000142
【０３５４】
ただし、信号パラメータ分析部５２０における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、（１７７）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。なお、信号パラメータ分析部５２０は、（１７７）式を複数回連続して適用することにより同時無相関化行列Ｐの推定値＾Ｐを更新してもよい。
【０３５５】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ７０）。収束判定部が収束していないと判定した場合（ステップＳ７０：Ｎｏ）、同時無相関化部５３０によるステップＳ６４の処理に戻って、処理が継続される。
【０３５６】
一方、収束判定部が収束したと判定した場合（ステップＳ７０：Ｙｅｓ）、同時無相関化部５３０における処理（ステップＳ７１）と、シングルチャネルウィーナーフィルタ部５４０における処理（ステップＳ７２）と、同時無相関化逆変換部５８０における処理（ステップＳ７３）と、が行われる。同時無相関化部５３０及びシングルチャネルウィーナーフィルタ部５４０における処理については上述の通りであるから、同時無相関化逆変換部５８０での処理について説明する。同時無相関化逆変換部５８０は、信号パラメータ分析部５２０からの同時無相関化行列Ｐの推定値＾Ｐと、シングルチャネルウィーナーフィルタ部５４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、に基づいて、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を次の（１７８）式により計算し、出力する（ステップＳ７３）。
【０３５７】
【数１４１】
2019163487.tif 000143
【０３５８】
［信号分析装置の処理の背景］
次に、本第５の実施形態に係る信号分析装置５０１の処理の背景について説明する。本第５の実施形態においては、観測信号ベクトルｙ（ｋ）は、Ｊ個の源信号成分ｘ１（ｋ），・・・，ｘＪ（ｋ）の和として、次の（１７９）式でモデル化される。
【０３５９】
【数１４２】
2019163487.tif 000144
【０３６０】
源信号成分ｘｊ（ｋ）は、（９７）式のように平均０、共分散行列Ｓｊ（ｋ）＝ｖｊ（ｋ）Ｒｊの複素ガウス分布に従うとする。ここで、ｖｊ（ｋ）はｊ番目の源信号のサンプル点ごとの分散をモデル化する分散パラメータであり、Ｒｊはｊ番目の源信号の空間的特性をモデル化するパラメータである空間共分散行列である。分散パラメータｖｊ（ｋ）は正であり、空間共分散行列Ｒｊは正定値エルミート行列であると仮定する。Ｊ個の源信号成分ｘ１（ｋ），・・・，ｘＪ（ｋ）は互いに無相関であるとする。
【０３６１】
観測信号ベクトルｙ（ｋ）から源信号成分ｘｊ（ｋ）を推定する方法としては、マルチチャネルウィーナーフィルタに基づく方法がある。この方法では、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を次の（１８０）式により計算する。
【０３６２】
【数１４３】
2019163487.tif 000145
【０３６３】
本第５の実施形態のように、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の状況で（１８０）式により源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する（すなわち信号分離を実現する）ためには、これらのパラメータを推定する必要がある。
【０３６４】
［従来の信号分析方法］
分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の状況で、分散パラメータｖｊ（ｋ）と、空間共分散行列Ｒｊと、源信号成分ｘｊ（ｋ）と、を同時に推定することにより、信号分離を実現する方法が、非特許文献１に開示されている。
【０３６５】
従来技術では、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊを最尤推定する。この最尤推定は、下記の（１８１）式の最適化問題として定式化される。
【０３６６】
【数１４４】
2019163487.tif 000146
【０３６７】
ここで、Ψはパラメータの全体を表し、分散パラメータｖｊ（ｋ）（ｊ＝１，・・・，Ｊ；ｋ＝１，・・・，Ｋ）及び空間共分散行列Ｒｊ（ｊ＝１，・・・，Ｊ）からなる。Ｌ１（Ψ）は、次の（１８２）式の対数尤度関数を表す。
【０３６８】
【数１４５】
2019163487.tif 000147
【０３６９】
下の（１８３）式は、行列Ａが正定値エルミート行列であることを表す。
【０３７０】
【数１４６】
2019163487.tif 000148
【０３７１】
従来技術では、ＥＭアルゴリズムに基づいて最尤推定を行う。このＥＭアルゴリズムでは、ＥステップとＭステップとが交互に反復される。
【０３７２】
Ｅステップでは、パラメータΨの現在の推定値＾Ψに基づく、源信号成分ｘｊ（ｋ）の事後確率ｐ（ｘｊ（ｋ）｜ｙ（ｋ），＾Ψ）を計算する。この事後確率は、次の（１８４）式のように複素ガウス分布であることが示される。
【０３７３】
【数１４７】
2019163487.tif 000149
【０３７４】
ただし、平均＾ｘｊ（ｋ）及び共分散行列Σｊ（ｋ）は、パラメータの現在の推定値＾Ψを用いて次の（１８５）式及び（１８６）式のように表される（＾ｖｊ（ｋ）は分散パラメータｖｊ（ｋ）の現在の推定値であり、＾Ｒｊは空間共分散行列Ｒｊの現在の推定値である）。
【０３７５】
【数１４８】
2019163487.tif 000150
【数１４９】
2019163487.tif 000151
【０３７６】
したがって、Ｅステップでは、平均＾ｘｊ（ｋ）及び共分散行列Σｊ（ｋ）を、（１８５）式及び（１８６）式により更新する。なお、（１８５）式は、（１８０）式のマルチチャネルウィーナーフィルタによる源信号成分ｘｊ（ｋ）の推定値に他ならない。以下、＾ｘｊ（ｋ）を源信号成分ｘｊ（ｋ）の推定値と呼び、Σｊ（ｋ）を事後共分散行列と呼ぶ。
【０３７７】
Ｍステップでは、Ｅステップで更新された源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）及び事後共分散行列Σｊ（ｋ）に基づいて、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値＾Ｒｊを更新する。これらは、下記の更新則（（１８７）式及び（１８８）式）により更新される。
【０３７８】
【数１５０】
2019163487.tif 000152
【数１５１】
2019163487.tif 000153
【０３７９】
しかしながら、従来技術では、（１８５）式の源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）の更新と、（１８６）式の事後共分散行列Σｊ（ｋ）の更新と、において、逆行列計算と行列乗算とをサンプル点ごとに行う必要があった。具体的には、（１８９）式に示す逆行列計算（反復ごとサンプル点ごとに１回）と、（１９０）式に示す行列乗算（反復ごとサンプル点ごとに２Ｊ回）とを行う必要があった。
【０３８０】
【数１５２】
2019163487.tif 000154
【数１５３】
2019163487.tif 000155
【０３８１】
そのため、従来技術では膨大な計算量を要していたため、リアルタイムでの信号分離や、センサネットワーク等における分散処理など、サンプル点の数Ｋが大きい（観測信号長が長い）場合や計算資源が限られている状況への適用が困難であった。
【０３８２】
［第５の実施形態に係る信号分析方法］
これに対し、本第５の実施形態では、Ｊ個の源信号成分ｘ１（ｋ），・・・，ｘＪ（ｋ）の同時無相関化に基づいて、Ｅステップにおける源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）及び事後共分散行列Σｊ（ｋ）の更新を少ない計算量で実現することにより、少ない計算量での信号分離を可能にする。以下、本第５の実施形態における考え方について説明する。
【０３８３】
いま仮に、（１８５）式及び（１８６）式における源信号成分ｘｊ（ｋ）の共分散行列ｖｊ（ｋ）Ｒｊがいずれも対角行列であるとすると、（１８５）式及び（１８６）式における逆行列演算及び行列乗算は単なる対角要素ごとの逆数演算及び乗算に帰着するため、少ない計算量で計算できる。しかしながら、通常、ベクトルｘｊ（ｋ）の異なる要素（つまり異なる位置で観測されたｊ番目の源信号）は互いに相関を持つため、その共分散行列ｖｊ（ｋ）Ｒｊの非対角要素は０ではない、すなわち共分散行列ｖｊ（ｋ）Ｒｊは対角行列ではない。
【０３８４】
そこで本第５の実施形態では、Ｊ個の源信号成分ｘ１（ｋ），・・・，ｘＪ（ｋ）を、正則行列Ｐのエルミート転置ＰＨによる次の（１９１）式の変換により同時無相関化することを考える。
【０３８５】
【数１５４】
2019163487.tif 000156
【０３８６】
変換後のｘｊ´（ｋ）の共分散行列は、次の（１９２）式及び（１９３）式で与えられる。
【０３８７】
【数１５５】
2019163487.tif 000157
【０３８８】
変換後のｘｊ´（ｋ）（ｊ＝１，・・・，Ｊ）がすべて無相関になる条件、すなわち、（１９３）式で与えられるこれらの共分散行列がすべて対角行列になる条件は、次の（１９４）式〜（１９６）式である。
【０３８９】
【数１５６】
2019163487.tif 000158
【０３９０】
既に述べたように、信号数がＪ＝２の場合には、一般化固有値問題を解くことにより、（１９４）式〜（１９６）式を満たす正則行列Ｐ及び対角行列Λ１，・・・，ΛＪを得ることができる。これに対し、本第５の実施形態では、Ｊ＝２の場合に限らず一般の信号数Ｊの場合に適用可能なアプローチとして、未知のＪ個の空間共分散行列Ｒ１，・・・，ＲＪが、ある正則行列Ｐにより同時対角化されるという制約条件、すなわち、ある正則行列Ｐ及びあるＪ個の対角行列Λ１，・・・，ΛＪに対し、（１９４）式〜（１９６）式が成立するという制約条件の下、正則行列Ｐ及び対角行列Λ１，・・・，ΛＪを最尤推定する。上述の制約条件は、ある正則行列Ｐ及びあるＪ個の対角行列Λ１，・・・，ΛＪに対し、次の各式（（１９７）式〜（１９９）式）が成立することと同値である。
【０３９１】
【数１５７】
2019163487.tif 000159
【０３９２】
本第５の実施形態では、空間共分散行列Ｒ１，・・・，ＲＪは、パラメータである正則行列Ｐとパラメータである対角行列Λ１，・・・，ΛＪとを用いて、（１９７）式〜（１９９）式のように表される（モデル化される）と仮定する。この場合、空間共分散行列Ｒ１，・・・，ＲＪを推定することは、行列Ｐ（同時無相関化を実現する行列であるから、以下、同時無相関化行列と呼ぶ）及び対角行列Λ１，・・・，ΛＪ（同時無相関化行列Ｐを用いて（１９４）式〜（１９６）式のように空間共分散行列Ｒ１，・・・，ＲＪを変換して得られる行列であるから、以下、同時無相関化された空間共分散行列と呼ぶ）を推定することに帰着する。
【０３９３】
そこで、本第５の実施形態では、同時無相関化行列Ｐと、同時無相関化された空間共分散行列Λｊと、分散パラメータｖｊ（ｋ）と、を最尤法により同時推定する。これにより、一般の信号数Ｊに対しても、源信号成分の同時無相関化（源信号成分の共分散行列の同時対角化）が実現し、少ない計算量での信号分離が可能になる。
【０３９４】
上述の最尤推定は、次の最適化問題（（２００）式）として定式化される。
【０３９５】
【数１５８】
2019163487.tif 000160
【０３９６】
ここで、Θはすべてのパラメータの全体を表し、分散パラメータｖｊ（ｋ）（ｊ＝１，・・・，Ｊ；ｋ＝１，・・・，Ｋ）と、同時無相関化行列Ｐと、同時無相関化された空間共分散行列Λｊ（ｊ＝１，・・・，Ｊ）と、からなる。Ｌ２（Θ）は、次の（２０１）式の対数尤度関数であり、（１８２）式に（１９７）式〜（１９９）式を代入することにより得られる。なお、ＧＬ（Ｉ，Ｃ）はすべての正則なＩ次複素行列からなる集合である。
【０３９７】
【数１５９】
2019163487.tif 000161
【０３９８】
本第５の実施形態では、ＥＭアルゴリズムにより、上述の最尤推定を実現する。このＥＭアルゴリズムはＥステップとＭステップとからなる。
【０３９９】
Ｅステップでは、パラメータΘの現在の推定値＾Θに基づいて、源信号成分の事後確率を更新する。従来技術の場合と同様に、ｘｊ（ｋ）の事後確率ｐ（ｘｊ（ｋ）｜ｙ（ｋ），＾Θ）は、（１８４）式の右辺の複素ガウス分布であり、その平均は、（１８５）式で、その共分散行列は、（１８６）式で与えられることが示せる。ただし、従来技術の場合とは異なり、（１８５）式及び（１８６）式における空間共分散行列Ｒｊの推定値＾Ｒｊは、同時無相関化行列Ｐの推定値＾Ｐ及び同時無相関化された空間共分散行列Λｊの推定値＾Λｊを用いて、次の（２０２）式〜（２０４）式のように表される。
【０４００】
【数１６０】
2019163487.tif 000162
【０４０１】
（１８５）式及び（１８６）式に（２０２）式〜（２０４）式を代入することにより、次の（２０５）式〜（２０８）式を得る。
【０４０２】
【数１６１】
2019163487.tif 000163
【数１６２】
2019163487.tif 000164
【０４０３】
したがって、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を同時無相関化行列Ｐの推定値＾Ｐにより基底変換して得られる＾ｘｊ´（ｋ）（以下、同時無相関化された源信号成分の推定値と呼ぶ）は、（１７３）式により更新でき、事後共分散行列Σｊ（ｋ）を同時無相関化行列Ｐの推定値＾Ｐにより基底変換して得られるΣｊ´（ｋ）（以下、同時無相関化された事後共分散行列と呼ぶ）は、（１７４）式により更新できる。その際、同時対角化により、逆行列演算及び行列乗算は、対角行列に対する逆行列演算及び行列乗算になっているため、その計算量は、Ｏ（Ｉ３）ではなくＯ（Ｉ）であり少ない計算量で計算できる。
【０４０４】
Ｍステップでは、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、同時無相関化された空間共分散行列Λｊの推定値＾Λｊと、同時無相関化行列Ｐの推定値＾Ｐと、を下記のＱ関数の最大化に基づいて更新する（（２０９）式及び（２１０）式参照）。ただし、（２０９）式及び（２１０）式ではΘに依らない定数項は無視してあり、したがって等号はその両辺がΘに依らない定数の差を除いて等しいことを表している。
【０４０５】
【数１６３】
2019163487.tif 000165
【０４０６】
（２１０）式のＱ関数を分散パラメータｖｊ（ｋ）に関して偏微分して０とおくと、次の（２１１）式及び（２１２）式を得る。
【０４０７】
【数１６４】
2019163487.tif 000166
【０４０８】
これで、（１７５）式の分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）の更新式が導かれた。
【０４０９】
また、（２１０）式のＱ関数を同時無相関化された空間共分散行列Λｊに関して偏微分したものが零行列に等しいと置くと、次の（２１３）式及び（２１４）式を得る。
【０４１０】
【数１６５】
2019163487.tif 000167
【０４１１】
これで、（１７６）式の同時無相関化された空間共分散行列Λｊの推定値＾Λｊの更新則が導かれた。なお、（２１３）式において、絶対値二乗｜・｜２はベクトルの要素ごとに計算するものとし、ベクトルαに対してｄｉａｇ（α）は、αの要素を対角要素とする対角行列を表す。
【０４１２】
（２０９）式のＱ関数を同時無相関化行列Ｐの複素共役Ｐ＊に関して偏微分すると、次の（２１５）式を得る。
【０４１３】
【数１６６】
2019163487.tif 000168
【０４１４】
これが零行列に等しいと置いて両辺のｖｅｃを取ると、次の（２１６）式を得る。
【０４１５】
【数１６７】
2019163487.tif 000169
【０４１６】
ここで、行列Ａに対し、ｖｅｃ（Ａ）は、行列Ａの列ベクトルを縦に並べて得られる列ベクトルを表す。これに関し、次の（２１７）式の公式が成り立つ。
【０４１７】
【数１６８】
2019163487.tif 000170
【０４１８】
また、「×」を丸で囲んだ記号は、クロネッカ積を表す。ブロック対角行列の性質より、（２１６）式は次の（２１８）式〜（２２０）式と等価である。
【０４１９】
【数１６９】
2019163487.tif 000171
【０４２０】
そこで、不動点法に基づいて、（１７７）式により、同時無相関化行列Ｐの推定値＾Ｐを更新すれば良い。このように不動点法に基づくことで、勾配法や自然勾配法とは異なり学習率の調整が不要であるため、安定な最適化が可能である。
【０４２１】
以下、事後確率の詳細な導出を述べる。隠れ変数ｘ１（ｋ），・・・，ｘＪ−１（ｋ）の対数事後確率は、次の（２２１）式〜（２２３）式のように書き下せる。
【０４２２】
【数１７０】
2019163487.tif 000172
【０４２３】
ここで、等号「＝」の上に「ｃ」を付した記号は、両辺がｘ１（ｋ），・・・，ｘＪ−１（ｋ）に依らない定数の差を除いて等しいことを表す。行列Ａは、Ｉ次の単位行列ＩをＪ−１個横に並べた次の（２２４）式の行列を表す。また、空間共分散行列Ｒｊ（ｊ＝１，・・・，Ｊ）は（１９７）式〜（１９９）式で与えられる。なお、本明細書では、Ｉを観測信号の個数（スカラー）と単位行列（行列）との２通りの意味で用いているが、Ｉがスカラーであるか行列であるかは文脈より明らかであるから、これにより混乱が生じるおそれはない。
【０４２４】
【数１７１】
2019163487.tif 000173
【０４２５】
~ｘ（ｋ）は、ｘ１（ｋ），・・・，ｘＪ−１（ｋ）を縦に並べた次の（２２５）式のベクトルを表す。
【０４２６】
【数１７２】
2019163487.tif 000174
【０４２７】
（２２３）式より、隠れ変数ｘ１（ｋ），・・・，ｘＪ−１（ｋ）の事後確率が次の（２２６）式の複素ガウス分布であることが導かれる。
【０４２８】
【数１７３】
2019163487.tif 000175
【０４２９】
ここで、共分散行列~Σ（ｋ）及び平均＾~ｘ（ｋ）は次の（２２７）式〜（２３０）式及び（２３１）式〜（２３３）式で与えられる。但し、（２２７）式から（２２８）式への式変形の際に逆行列補題を用いた。
【０４３０】
【数１７４】
2019163487.tif 000176
【０４３１】
【数１７５】
2019163487.tif 000177
【０４３２】
後に分かるように、Ｅステップでは、行列~Σ全体を計算する必要は実はなく、行列~ΣをＩ次の正方行列を要素とするＪ−１次の正方ブロック行列とみなしたときの対角ブロックのみを計算すれば良い。ｊ番目の対角ブロックをΣｊ（ｋ）と書くことにすると、これは次の（２３４）式のように計算できる。明らかに、Ｊ＝２の場合には、（２３４）式は（１０１）式と等価である。
【０４３３】
【数１７６】
2019163487.tif 000178
【０４３４】
また、＾ｘｊ（ｋ）を、＾~ｘ（ｋ）をＩ次元縦ベクトルを要素とするサイズ（Ｊ−１）×１のブロック行列とみなしたときの（ｊ，１）ブロックと定義すると、これは次の（２３５）式のように表される。
【０４３５】
【数１７７】
2019163487.tif 000179
【０４３６】
＾ｘｊ（ｋ）及びΣｊ（ｋ）は次の（２３６）式及び（２３７）式のようにも表されることに注意する。
【０４３７】
【数１７８】
2019163487.tif 000180
【数１７９】
2019163487.tif 000181
【０４３８】
ただし、上式においては、ｊ＝１，・・・，Ｊ−１である。また、Ｅ（・｜ｙ（ｋ））は、隠れ変数の事後確率ｐ（ｘ１（ｋ），・・・，ｘＪ−１（ｋ）｜ｙ（ｋ），Θ）に関する期待値（（２３８）式）である。
【０４３９】
【数１８０】
2019163487.tif 000182
【０４４０】
これにならい、ｊ＝Ｊに対しても、＾ｘＪ（ｋ）及びΣＪ（ｋ）を次の（２３９）式〜（２４３）式及び（２４４）式〜（２５１）式のように定義する。
【０４４１】
【数１８１】
2019163487.tif 000183
【０４４２】
【数１８２】
2019163487.tif 000184
【０４４３】
よって、以下の（２５２）式及び（２５３）式は、ｊ＝１，・・・，Ｊに対して成立する。
【０４４４】
【数１８３】
2019163487.tif 000185
【数１８４】
2019163487.tif 000186
【０４４５】
以下、Ｑ関数の詳細な導出について説明する。観測信号ベクトルｙ（ｋ）と隠れ変数ｘ１（ｋ），・・・，ｘＪ−１（ｋ）との同時分布の対数は、次の（２５４）式〜（２５６）式のように書き下される。ただし、（２５４）式〜（２５６）式では定数項は無視してあり、したがって等号はその両辺が定数の差を除いて等しいことを表している。
【０４４６】
【数１８５】
2019163487.tif 000187
【０４４７】
ここで、次の各式（（２５７）式〜（２６２）式）が成立することに注意する。
【０４４８】
【数１８６】
2019163487.tif 000188
【０４４９】
したがって、Ｑ関数は次式（（２６３）式〜（２６６）式）で与えられる。（２６４）式は（２０９）式に他ならない。
【０４５０】
【数１８７】
2019163487.tif 000189
【０４５１】
ここで、（２６６）式における絶対値の二乗｜・｜２は、要素ごとに取るものとする。また、等号「＝」の上に「ｃ」を付した記号は、両辺がΘに依らない定数の差を除いて等しいことを表す（先と異なる意味で記号が用いられていることに注意されたい）。
【０４５２】
［第５の実施形態の変形例１］
次に、第５の実施形態の変形例１について説明する。本第５の実施形態の変形例１では、信号パラメータ分析部における処理の変形例について説明する。第５の実施形態では、信号パラメータ分析部は、不動点法に基づいて同時無相関化行列Ｐの推定値＾Ｐを更新していた。これに対し、第５の実施形態の変形例１に掛かる信号パラメータ分析部は、不動点法以外の任意の最適化法（例えば、勾配法、自然勾配法、ニュートン法など）に基づいて、同時無相関化行列Ｐの推定値＾Ｐを更新する。
【０４５３】
例えば、勾配法の場合、本第５の実施形態の変形例１に係る信号パラメータ分析部は、同時無相関化行列Ｐの推定値＾Ｐを次の（２６７）式により更新する（これは、（２１５）式から直ちに導かれる）。
【０４５４】
【数１８８】
2019163487.tif 000190
【０４５５】
ただし、信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、（２６７）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。μは学習率と呼ばれる正の数であり、例えばμ＝０．０５などと定めればよい。
【０４５６】
また、自然勾配法の場合、本第５の実施形態の変形例１に係る信号パラメータ分析部は、同時無相関化行列Ｐの推定値を次の（２６８）式により更新する（これも（２１５）式から直ちに導かれる）。
【０４５７】
【数１８９】
2019163487.tif 000191
【０４５８】
ただし、信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、（２６８）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。μは、学習率と呼ばれる正の数であり、例えば、μ＝０．０５などと定めればよい。自然勾配法によれば、空間の計量が適切に考慮されるため、一般に勾配法よりも安定な更新が可能であることが知られている。
【０４５９】
［第５の実施形態の変形例２］
次に、第５の実施形態の変形例２について説明する。本第５の実施形態の変形例２では、信号パラメータ分析部における処理の変形例について説明する。第５の実施形態や第５の実施形態の変形例１では、信号パラメータ分析部は、Ｑ関数に基づいて同時無相関化行列Ｐの推定値＾Ｐを更新していた。これに対し、本第５の実施形態の変形例２に係る信号パラメータ分析部は、直接、対数尤度関数Ｌ２（Θ）に基づいて、同時無相関化行列Ｐの推定値＾Ｐを更新する。
【０４６０】
本第５の実施形態の変形例２に係る信号パラメータ分析部は、同時無相関化部からの同時無相関化された観測信号ベクトルｙ´（ｋ）と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊと、分散パラメータ更新部からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、（信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、更に、初期化部からの同時無相関化行列Ｐの初期推定値＾Ｐ（０）と、）に基づいて、同時無相関化行列Ｐの推定値＾Ｐを次式（（２６９）式）により更新する（Ｅは単位行列を表すものとする）。
【０４６１】
【数１９０】
2019163487.tif 000192
【０４６２】
ただし、信号パラメータ分析部における＾Ｐの最初の更新の際には、（２６９）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。なお、信号パラメータ分析部は、（２６９）式を複数回連続して適用することにより、同時無相関化行列Ｐの推定値＾Ｐを更新してもよい。
【０４６３】
（２６９）式の対数尤度関数Ｌ２（Θ）に基づく同時無相関化行列Ｐの推定値＾Ｐの更新則の導出について説明する。（２０１）式の対数尤度関数を同時無相関化行列Ｐの複素共役Ｐ＊に関して偏微分すると次の（２７０）式を得る。
【０４６４】
【数１９１】
2019163487.tif 000193
【０４６５】
（２７０）式が零行列に等しいと置いて、両辺のｖｅｃを取ると、次式（（２７１）式）を得る。
【０４６６】
【数１９２】
2019163487.tif 000194
【０４６７】
（２７１）式の右辺の行列は、次の（２７２）式及び（２７３）式のように変形できる。
【０４６８】
【数１９３】
2019163487.tif 000195
【０４６９】
よって次式（（２７４）式及び（２７５）式）を得る。
【０４７０】
【数１９４】
2019163487.tif 000196
【０４７１】
（２７５）式に（２７６）式を代入することにより、（２６９）式の不動点法による同時無相関化行列Ｐの推定値＾Ｐの更新則が得られる。
【０４７２】
【数１９５】
2019163487.tif 000197
【０４７３】
このように、不動点法に基づくことで、勾配法や自然勾配法とは異なり学習率の調整が不要であるため、安定な最適化が可能である。
【０４７４】
また、第５の実施形態の変形例１と同様に、本第５の実施形態の変形例２においても、不動点法に限らず任意の最適化法（勾配法、自然勾配法、ニュートン法など）により同時無相関化行列Ｐの推定値＾Ｐを更新してよい。以下、そのような例について説明する。
【０４７５】
勾配法の場合、本第５の実施形態の変形例２に係る信号パラメータ分析部は、観測信号ベクトル作成部からの観測信号ベクトルｙ（ｋ）と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊと、分散パラメータ更新部からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、（信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、更に、初期化部からの同時無相関化行列Ｐの初期推定値＾Ｐ（０）と、）に基づいて、同時無相関化行列Ｐの推定値＾Ｐを次の（２７７）式により更新する（これは、（２７０）式より直ちに導かれる）。
【０４７６】
【数１９６】
2019163487.tif 000198
【０４７７】
ただし、信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、（２７７）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。なお、信号パラメータ分析部は、（２７７）式を複数回連続して適用することにより同時無相関化行列Ｐの推定値＾Ｐを更新してもよい。また、信号パラメータ分析部は、（２７７）式に（２７６）式を代入することにより得られる、同時無相関化された観測信号ベクトルｙ´（ｋ）を用いて表した式を、（２７７）式の代わりに用いても良い。μは、学習率と呼ばれる所定の正の数であり、例えば、μ＝０．０５などと定めればよい。この勾配法の場合には、（２７７）式の１回の適用につき逆行列計算が１回のみで行列乗算は不要（対角行列に対する逆行列計算はカウントしていない）であるため、不動点法の場合と比べてより少ない計算量で同時無相関化行列Ｐの推定値＾Ｐを更新できる。
【０４７８】
また、自然勾配法の場合、本第５の実施形態の変形例２に係る信号パラメータ分析部は、観測信号ベクトル作成部からの観測信号ベクトルｙ（ｋ）と、空間共分散行列更新部からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊと、分散パラメータ更新部からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、（信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、更に、初期化部からの同時無相関化行列Ｐの初期推定値＾Ｐ（０）と、）に基づいて、同時無相関化行列Ｐの推定値＾Ｐを次の（２７８）式により更新する（これも（２７０）式より直ちに導かれる）。
【０４７９】
【数１９７】
2019163487.tif 000199
【０４８０】
ただし、信号パラメータ分析部における同時無相関化行列Ｐの推定値＾Ｐの最初の更新の際には、（２７８）式の右辺において、＾Ｐを＾Ｐ（０）で置き換える。なお、信号パラメータ分析部は、（２７８）式を複数回連続して適用することにより同時無相関化行列Ｐの推定値＾Ｐを更新してもよい。また、信号パラメータ分析部は、（２７８）式に（２７６）式を代入することにより得られる、同時無相関化された観測信号ベクトルｙ´（ｋ）を用いて表した式を、（２７８）式の代わりに用いても良い。μは、学習率と呼ばれる所定の正の数であり、例えば、μ＝０．０５などと定めればよい。また、自然勾配法によれば、空間の計量が適切に考慮されるため、一般に勾配法よりも安定な更新が可能であることが知られている。
【０４８１】
［第５の実施形態の変形例３］
次に、第５の実施形態の変形例３について説明する。第５の実施形態の変形例３では、第３の実施形態にならい、周波数ごとに第５の実施形態を適用する。
【０４８２】
第５の実施形態の変形例３においては、Ｊ個の音源信号が混在する状況において、それぞれ異なる位置でマイクロホンにより取得されたＩ個（Ｉは２以上の整数）の観測信号ｙｉ（ｔ）が信号分析装置に入力されるものとする。観測信号ｙｉ（ｔ）は、周波数分解され、周波数ごとに第５の実施形態に記載の方法で信号分離がなされ、最後に時間領域に変換される。
【０４８３】
図５及び図６を用いて、第５の実施形態の変形例３に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。
【０４８４】
周波数領域変換部３１０（ステップＳ４０、ステップＳ４１）及び時間領域変換部３３０（ステップＳ４３）での処理は、第３の実施形態と同じである。以下、相違点である周波数領域信号分離部３２０（ステップＳ４２）での処理について説明する。
【０４８５】
周波数領域信号分離部３２０は、周波数領域変換部３１０から（２７９）式に示す時間周波数領域の観測信号を受け取り、ｆごとに第５の実施形態の方法により周波数領域での信号分離を行うことにより、（２８０）式に示す時間周波数領域の源信号成分の推定値を計算し、時間領域変換部３３０に渡す（ステップＳ４２）。
【０４８６】
【数１９８】
2019163487.tif 000200
【数１９９】
2019163487.tif 000201
【０４８７】
すなわち、ｆ＝１に対しては、（２８１）式に示す観測信号を入力として、第５の実施形態の方法により（２８２）式の源信号成分の推定値を得る。
【０４８８】
【数２００】
2019163487.tif 000202
【数２０１】
2019163487.tif 000203
【０４８９】
そして、ｆ＝２に対しては、（２８３）式に示す観測信号を入力として、第５の実施形形態の方法により（２８４）式の源信号成分の推定値を得る。
【０４９０】
【数２０２】
2019163487.tif 000204
【数２０３】
2019163487.tif 000205
【０４９１】
以下、同じ処理を各ｆに対して繰り返し、最後に、ｆ＝Ｆに対しては、（２８５）式に示す観測信号を入力として、第５の実施形態の方法により（２８６）式の源信号成分の推定値を得る。
【０４９２】
【数２０４】
2019163487.tif 000206
【数２０５】
2019163487.tif 000207
【０４９３】
こうして得られた時間周波数領域の源信号成分の推定値（（２８７）式参照）を時間領域変換部３３０に渡す。
【０４９４】
【数２０６】
2019163487.tif 000208
【０４９５】
［第５の実施形態の効果］
本発明の効果を確認するために、従来の方法及び第５の実施形態を用いた音源分離実験について説明する。
【０４９６】
図１４は、従来の方法及び第５の実施形態を用いた音源分離実験の実験条件を示す図である。図１５は、従来の方法及び第５の実施形態を用いた音源分離実験結果（real  time  factor）を示す図である。図１６は、従来の方法及び第５の実施形態を用いた音源分離実験結果（音源分離性能）を示す図である。
【０４９７】
８秒の英語音声に、図１４に示す実験室で測定された残響インパルス応答を畳み込んで残響音声を生成し、複数人の残響音声を加算することで観測信号を生成した。このとき、観測信号に対して、各方法を用いて信号分離を行った場合の性能を図１５及び図１６に示す。図１５は、real  time  factor（観測信号長を１としたときの、処理に要した時間）である。また、図１６は、音源分離性能である信号対歪み比（値が大きいほど音源分離性能が高い）である。
【０４９８】
図１５及び図１６に示すように、本第５の実施形態に係る方法によれば、音源数が３個以上である場合にも、４２０倍以上の計算時間削減が実現でき、従来の方法を上回る音源分離性能を得ることができた。
【０４９９】
このように、本第５の実施形態によれば、従来の方法と比べて計算時間を大きく削減できるのに加え、従来の方法よりも高い信号分離性能を実現できる。本第５の実施形態は、信号数が２個でない場合にも適用でき、信号数が３個以上の場合にも適用できる。
【０５００】
［第６の実施形態］
次に、第６の実施形態に係る信号分析装置の構成、処理の流れ及び効果について説明する。第６の実施形態では、第４の実施形態の変形例を実施するための具体的な手法について説明する。
【０５０１】
本第６の実施形態では、第５の実施形態と同様に、一般にＪ個の源信号が混在し、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の状況で、分散パラメータｖｊ（ｋ）と、空間共分散行列Ｒｊをモデル化するパラメータである同時無相関化された空間共分散行列Λｊ及び同時無相関化行列Ｐと、源信号成分ｘｊ（ｋ）と、を同時に推定することにより、信号分離を実現する方法について説明する。
【０５０２】
なお、第５の実施形態は、ＥＭアルゴリズムに基づいていたのに対し、本第６の実施形態は補助関数法に基づいている。本第６の実施形態では、それぞれ異なる位置で取得されたＩ個の観測信号ｙｉ（ｋ）が信号分析装置に入力されるものとする。
【０５０３】
［信号分析装置の構成及び処理］
図１７及び図１８を用いて、第６の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図１７は、第６の実施形態に係る信号分析装置の構成の一例を示す図である。図１８は、第６の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０５０４】
図１７に示すように、信号分析装置６０１は、初期化部（図示しない）、観測信号ベクトル作成部６１０、同時無相関化部６３０、シングルチャネルウィーナーフィルタ部６４０、観測共分散行列更新部６５０、分散パラメータ更新部６６０、空間共分散行列更新部６７０、信号パラメータ分析部６２０、収束判定部（図示しない）及び同時無相関化逆変換部６８０を有する。
【０５０５】
まず、信号分析装置の各部の概要について説明する。初期化部（図示しない）の処理（ステップＳ８１）、観測信号ベクトル作成部６１０の処理（ステップＳ８２，Ｓ８３）、同時無相関化部６３０の処理（ステップＳ８４）及びシングルチャネルウィーナーフィルタ部６４０の処理（ステップＳ８５）については、第５の実施形態と同様であるから、説明を省略する。
【０５０６】
観測共分散行列更新部６５０は、分散パラメータ更新部６６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）（ただし例外として、観測共分散行列更新部６５０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ）の初期推定値＾ｖｊ（０）（ｋ））と、空間共分散行列更新部６７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊ（ただし例外として、観測共分散行列更新部６５０における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０））と、を受け取って、同時無相関化された観測共分散行列Ｕ´（ｋ）を次式（２８８）式により更新する（ステップＳ８６）。
【０５０７】
【数２０７】
2019163487.tif 000209
【０５０８】
ただし、観測共分散行列更新部６５０における最初の処理の際には、上式において、＾ｖｊ（ｋ）を＾ｖｊ（０）（ｋ）に、＾Λｊを＾Λｊ（０）に、それぞれ置き換える。
【０５０９】
分散パラメータ更新部６６０は、シングルチャネルウィーナーフィルタ部６４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、観測共分散行列更新部６５０からの同時無相関化された観測共分散行列Ｕ´（ｋ）と、空間共分散行列更新部６７０からの同時無相関化された空間共分散行列Λｊの推定値＾Λｊ（ただし例外として、分散パラメータ更新部６６０における最初の処理の際には、初期化部からの同時無相関化された空間共分散行列Λｊの初期推定値＾Λｊ（０））と、を受け取って、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）を次式（（２８９）式）により更新する（ステップＳ８７）。
【０５１０】
【数２０８】
2019163487.tif 000210
【０５１１】
ただし、分散パラメータ更新部６６０における最初の処理の際には、上式において、＾Λｊを＾Λｊ（０）で置き換える。
【０５１２】
空間共分散行列更新部６７０は、シングルチャネルウィーナーフィルタ部６４０からの同時無相関化された源信号成分ｘｊ´（ｋ）の推定値＾ｘｊ´（ｋ）と、観測共分散行列更新部６５０からの同時無相関化された観測共分散行列Ｕ´（ｋ）と、分散パラメータ更新部６６０からの分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、を受け取って、同時無相関化された空間共分散行列Λｊの推定値＾Λｊを次の（２９０）式により更新する（ステップＳ８８）。
【０５１３】
【数２０９】
2019163487.tif 000211
【０５１４】
信号パラメータ分析部６２０は、第５の実施形態の変形例２と同様の処理を行って、同時無相関化行列Ｐの推定値＾Ｐを更新する（ステップＳ８９）。
【０５１５】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ９０）。収束判定部が収束していないと判定した場合（ステップＳ９０：Ｎｏ）、同時無相関化部６３０での処理（ステップＳ８４）に戻って、処理が継続される。
【０５１６】
一方、収束判定部が収束したと判定した場合（ステップＳ９０：Ｙｅｓ）、同時無相関化部６３０における処理（ステップＳ９１）、シングルチャネルウィーナーフィルタ部６４０における処理（ステップＳ９２）、及び、同時無相関化逆変換部６８０における処理（ステップＳ９３）が、第５の実施形態と同様にして行われることにより、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）が計算され、出力される。
【０５１７】
［信号分析装置の処理の背景］
次に、第６の実施形態に係る信号分析装置の背景について説明する。第５の実施形態で述べたように、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊが未知の状況で、（１８０）式により、源信号成分ｘｊ（ｋ）の推定値＾ｘｊ（ｋ）を計算する（すなわち信号分離を実現する）ためには、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊを推定する必要がある。
【０５１８】
補助関数法に基づいて、（１８１）式の最適化問題を解くことにより、分散パラメータｖｊ（ｋ）及び空間共分散行列Ｒｊを最尤推定することができる。
【０５１９】
この補助関数法では、補助変数Φ及び次式を満たす補助関数Ｌ１−（Ψ，Φ）に基づいて、補助関数Ｌ１−（Ψ，Φ）を補助変数Φについて最大化することによりΦを更新するステップと、補助関数Ｌ１−（Ψ，Φ）が減少しないようにパラメータΨを更新するステップと、が交互に反復される（（２９１）式参照）。
【０５２０】
【数２１０】
2019163487.tif 000212
【０５２１】
本第６の実施形態では、補助関数として（２９２）式のＬ１−（Ψ，Φ）を用いる。
【０５２２】
【数２１１】
2019163487.tif 000213
【０５２３】
ただし、右辺における「定数」は、パラメータΨ及び補助変数Φによらない定数である。（２９２）式が（２９１）式の補助関数の条件を満たすことは、参考文献３及び参考文献４と同様に示すことができる。なお、補助変数Φは、Ｕ（ｋ）及びＧｊ（ｋ）からなり、Ｕ（ｋ）は正定値エルミート行列であり、Ｇｊ（ｋ）は（２９３）式を満たすものとする。
【０５２４】
【数２１２】
2019163487.tif 000214
【０５２５】
補助関数Ｌ１−（Ψ，Φ）を補助変数Φについて最大化することによりΦを更新するステップでは、分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び空間共分散行列Ｒｊの推定値＾Ｒｊに基づいて、Ｕ（ｋ）及びＧｊ（ｋ）を（２９４）式及び（２９５）式により更新する。Ｕ（ｋ）は観測信号ベクトルｙ（ｋ）の共分散行列に他ならないから、観測共分散行列と呼ぶ。
【０５２６】
【数２１３】
2019163487.tif 000215
【数２１４】
2019163487.tif 000216
【０５２７】
補助関数Ｌ１−（Ψ，Φ）が減少しないようにパラメータΨを更新するステップでは、上で更新した補助変数に基づいて、次の（２９６）式及び（２９７）式により分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）と、空間共分散行列Ｒｊの推定値＾Ｒｊと、を更新する。
【０５２８】
【数２１５】
2019163487.tif 000217
【数２１６】
2019163487.tif 000218
【０５２９】
ここで、Ａ＃Ｂは、（２９８）式で定義されるＡとＢとの幾何平均である。
【０５３０】
【数２１７】
2019163487.tif 000219
【０５３１】
また、行列Ｃｊ及び行列Ｄｊは、それぞれ（２９９）式及び（３００）式により計算される。
【０５３２】
【数２１８】
2019163487.tif 000220
【数２１９】
2019163487.tif 000221
【０５３３】
（２９４）式〜（３００）式の更新式の導出は、参考文献３及び参考文献４と同様にできる。更に、補助変数Ｇｊ（ｋ）を代入消去することにより、補助変数Ｇｊ（ｋ）が現れない更新則を導くことができる。この場合の更新則（１回の反復）を以下の（３０１）式〜（３０５）式に示す。明らかに、Ｊ＝２の場合には、（３０１）式〜（３０５）式による更新は、（１５６）式〜（１５９）式による更新と等価である。
【０５３４】
【数２２０】
2019163487.tif 000222
【数２２１】
2019163487.tif 000223
【数２２２】
2019163487.tif 000224
【数２２３】
2019163487.tif 000225
【数２２４】
2019163487.tif 000226
【０５３５】
しかしながら、上記の方法では、（３０２）式〜（３０４）式に含まれる、逆行列計算Ｕ（ｋ）−１をサンプル点ごとに行う必要があった。そのため、上記の方法は膨大な計算量を要する。
【０５３６】
これに対し、本第６の実施形態では、第５の実施形態と同様に、Ｊ個の源信号成分ｘ１（ｋ），・・・，ｘＪ（ｋ）の同時無相関化に基づいて、逆行列計算Ｕ（ｋ）−１を少ない計算量で計算できる対角行列に対する逆行列計算に帰着させることにより、少ない計算量での信号分離を実現する。具体的には、本第６の実施形態では、第５の実施形態と同様に、Ｊ個の空間共分散行列が（１９７）式〜（１９９）式の形に表されるという制約条件の下、同時無相関化行列Ｐ及び同時無相関化された空間共分散行列Λ１，・・・，ΛＪを最尤推定する。この最尤推定は、（２００）式の最適化問題として定式化される。本第６の実施形態では、補助関数法に基づいて、上述の最尤推定を実現する。この方法では、補助関数法の１回の反復を適用することにより分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び同時無相関化された空間共分散行列Λｊの推定値＾Λｊを更新するステップと、同時無相関化行列Ｐの推定値＾Ｐを更新するステップと、が交互に反復される。
【０５３７】
補助関数法に基づいて分散パラメータｖｊ（ｋ）の推定値＾ｖｊ（ｋ）及び同時無相関化された空間共分散行列Λｊの推定値＾Λｊを更新するステップは、（３０６）式の補助関数Ｌ２−（Θ，Φ）に基づいている。
【０５３８】
【数２２５】
2019163487.tif 000227
【０５３９】
ただし、右辺の「定数」はパラメータΘ及び補助変数Φによらない定数である。
【０５４０】
先と同様に補助変数Ｇｊ（ｋ）を代入消去すると、更新則は下記の（３０７）式〜（３１１）式のように与えられる。ただし、（３１０）式の右辺における斜体のＤは、行列の非対角成分を０で置き換える演算子である。
【０５４１】
【数２２６】
2019163487.tif 000228
【数２２７】
2019163487.tif 000229
【数２２８】
2019163487.tif 000230
【数２２９】
2019163487.tif 000231
【数２３０】
2019163487.tif 000232
【０５４２】
これは、更に次の（３１２）式〜（３１５）式に書き直せる。
【０５４３】
【数２３１】
2019163487.tif 000233
【数２３２】
2019163487.tif 000234
【数２３３】
2019163487.tif 000235
【数２３４】
2019163487.tif 000236
【０５４４】
同時無相関化行列Ｐの推定値＾Ｐの更新則の導出については、第５の実施形態の変形例２と同様であるから説明を省略する。以上で、第６の実施形態における更新則が導かれた。
【０５４５】
［第６の実施形態の変形例］
第６の実施形態では対数尤度に基づいて同時無相関化行列Ｐの推定値＾Ｐを更新するが、補助関数に基づいて同時無相関化行列Ｐの推定値＾Ｐを更新してもよい。詳細は、第５の実施形態及び第５の実施形態の変形例１と同様だから省略する。
【０５４６】
［第６の実施形態の効果］
第６の実施形態によれば、同時無相関化に基づき、計算量を大幅に削減できる。また、ＥＭアルゴリズムは反復回数が少ない場合には必ずしも高い信号分離性能を実現できないが、第６の実施形態によれば、補助関数法に基づき、反復回数が少ない場合にも高い信号分離性能を実現できる。本第６の実施形態は、信号数が３個以上でも適用できる。
【０５４７】
［第７の実施形態］
本第７の実施形態では、観測信号に含まれる源信号の総数ＬがＬ≧３の場合でも少ない計算量で信号分離を可能にする方法を示す。第５の実施形態や第６の実施形態でもそのような方法を示したが、本実施形態ではもう一つの実現方法を示す。
【０５４８】
第５の実施形態や第６の実施形態は、「すべての時間周波数点において観測信号は、Ｌ個すべての源信号を含む」というモデル化（（１７９）式参照）の下、Ｌ個すべての源信号の同時無相関化（Ｌ個すべての源信号の空間共分散行列の同時対角化）に基づいて、Ｌ≧３の場合でも少ない計算量で信号分離を実現していた。この同時無相関化（同時対角化）は、「Ｌ個すべての源信号の空間共分散行列が１つの同時無相関化行列によって同時対角化できる」という制約条件（（１９７）式〜（１９９）式参照）に基づいていた。しかしながら、実際にはこの制約条件は現実に即していない場合もあり、そのような場合には、第５の実施形態や第６の実施形態ではモデル化精度が低下し、結果として信号分離性能が低下するおそれがあった。なお、第５の実施形態や第６の実施形態では、Ｌ個すべての源信号の同時無相関化（Ｌ個すべての源信号の空間共分散行列の同時対角化）に基づいているため、同時無相関化される源信号の最大個数（同時対角化される源信号の空間共分散行列の最大個数）をＪとすれば、Ｊ＝Ｌである。したがって、第５の実施形態や第６の実施形態では、これらをともにＪで表していた。
【０５４９】
一方、本第７の実施形態では、源信号の時間周波数領域におけるスパース性に着目する。源信号は、時間周波数領域におけるスパース性、すなわち、「少数の時間周波数成分にエネルギーが集中する」という性質を持つことが多い（例えば、参考文献５「O.  Yilmaz  and  S.  Rickard,  “Blind  separation  of  speech  mixtures  via  time-frequency  masking”,  IEEE  Transactions  on  Signal  Processing,  vol.  52,  no.  7,  pp.  1830−1847,  Jul.  2004.」参照）。この場合、各時間周波数点において観測信号は、Ｌ個すべての源信号を含むのではなく、Ｌ個よりも少ない源信号しか含まないとみなせることが多い。
【０５５０】
そこで、本第７の実施形態では、「各時間周波数点において観測信号は高々２個の源信号しか含まない」というモデル化の下、各時間周波数点において観測信号に含まれる高々２個の源信号の同時無相関化（各時間周波数点において観測信号に含まれる高々２個の源信号の空間共分散行列の同時対角化）に基づいて、Ｌ≧３の場合でも少ない計算量で信号分離を実現する。一般化固有値問題に基づいて２個以下の源信号は常に一つの同時無相関化行列によって同時無相関化できる（２個以下の源信号の空間共分散行列は常に一つの同時無相関化行列によって同時対角化できる）から、本第７の実施形態では、第５の実施形態や第６の実施形態とは異なり、源信号の空間共分散行列に特段の制約条件を課す必要がない。
【０５５１】
したがって、本第７の実施形態によれば、「Ｌ個すべての源信号の空間共分散行列が１つの同時無相関化行列によって同時対角化できる」という制約条件が現実に即していない場合でも高精度なモデル化が可能であり、結果として、そのような場合でも高い信号分離性能を実現可能である。なお、本第７の実施形態では、同時無相関化される源信号の最大個数Ｊ（同時対角化される源信号の空間共分散行列の最大個数Ｊ）はＪ＝２である。
【０５５２】
本第７の実施形態では、Ｌ個（Ｌは２以上の自然数）の源信号が混在する状況において、それぞれ異なる位置で取得されたＩ個（Ｉは２以上の自然数）の観測信号ｙｉ（ｔ）が信号分析装置に入力されるものとする。
【０５５３】
［信号分析装置の構成及び処理］
図５及び図６を用いて、本第７の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。図５は、本第７の実施形態に係る信号分析装置の構成の一例を示す図である。図６は、本第７の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。
【０５５４】
図５に示すように、信号分析装置３０１は、周波数領域変換部３１０、周波数領域信号分離部３２０、時間領域変換部３３０からなる。
【０５５５】
まず、信号分析装置３０１の各部の概要について説明する。周波数領域変換部３１０は、入力された観測信号ｙｉ（ｔ）を受け取り（ステップＳ４０）、短時間フーリエ変換などにより時間周波数領域の観測信号ｙｉ（ｋ，ｆ）を計算する（ステップＳ４１）。ここで、ｋは時間フレームのインデックスであり、ｆは周波数ビンのインデックスである。ｙｉ（ｋ，ｆ）をｙｉｆ（ｋ）とも表す。
【０５５６】
周波数領域信号分離部３２０は、ｆごとに信号分離を行って（ステップＳ４２）、時間周波数領域の源信号成分ｘｊ（ｋ，ｆ）の推定値＾ｘｊ（ｋ，ｆ）を得る。＾ｘｊ（ｋ，ｆ）を＾ｘｊｆ（ｋ）とも表す。
【０５５７】
時間領域変換部３３０は、時間周波数領域の源信号成分ｘｊ（ｋ，ｆ）の推定値＾ｘｊ（ｋ，ｆ）に基づいて、逆短時間フーリエ変換などにより、時間領域の源信号成分ｘｊ（ｔ）の推定値＾ｘｊ（ｔ）を計算し（ステップＳ４３）、出力する。
【０５５８】
次に、周波数領域信号分離部３２０についてより詳細に説明する。図１９及び図２０を用いて、周波数領域信号分離部３２０の構成及び周波数領域信号分離処理の処理手順について説明する。図１９は、周波数領域信号分離部３２０の構成の一例を示す図である。図２０は、周波数領域信号分離処理の処理手順の一例を示すフローチャートである。
【０５５９】
図１９に示すように、周波数領域信号分離部３２０は、初期化部（図示しない）、観測信号ベクトル作成部７１０、信号パラメータ分析部７２０、同時無相関化部７３０、信号存在事後確率更新部７４０、シングルチャネルウィーナーフィルタ部７５０、事後共分散行列更新部７６０、混合重み更新部７７０、分散パラメータ更新部７８０、空間共分散行列更新部７９０、収束判定部（図示しない）及び源信号成分推定部８００を有する。これら周波数領域信号分離部３２０の各部での処理は、周波数ビン毎に行う。
【０５６０】
周波数領域信号分離部３２０は、ｆ←０により周波数ビンのインデックスｆを初期化する（ステップＳ１０１）。周波数領域信号分離部３２０は、ｆ←ｆ＋１により周波数ビンのインデックスｆを１増やす（ステップＳ１０２）。
【０５６１】
図示しない初期化部は、混合重みαＢ（ｆ）の初期推定値＾αＢ（０）（ｆ）と、分散パラメータｖｊ（ｋ，ｆ）の初期推定値＾ｖｊ（０）（ｋ，ｆ）と、空間共分散行列Ｒｊ（ｆ）の初期推定値＾Ｒｊ（０）（ｆ）と、を設定する（ステップＳ１０３）。これらは、例えば乱数に基づいて設定すればよい。
【０５６２】
観測信号ベクトル作成部７１０は、周波数領域変換部からの時間周波数領域の観測信号ｙｉ（ｋ，ｆ）を受け取り（ステップＳ１０４）、観測信号ベクトルｙ（ｋ，ｆ）を（３１６）式により作成する（ステップＳ１０５）。
【０５６３】
【数２３５】
2019163487.tif 000237
【０５６４】
ここで、本第７の実施形態における観測信号ベクトルｙ（ｋ，ｆ）のモデル化について説明する。本実施形態では、各時間周波数点（ｋ，ｆ）において観測信号ベクトルｙ（ｋ，ｆ）は高々２個の源信号成分しか含まないと仮定する。各時間周波数点（ｋ，ｆ）において観測信号ベクトルｙ（ｋ，ｆ）に含まれる源信号成分は、そのような源信号成分のインデックスの全体である集合Ｂ（ｋ，ｆ）により表すことができる（ただし、集合Ｂ（ｋ，ｆ）の元の個数＃Ｂ（ｋ，ｆ）は＃Ｂ（ｋ，ｆ）≦２を満たすものとする）。このとき、観測信号ベクトルｙ（ｋ，ｆ）は、Ｂ（ｋ，ｆ）のすべての元ｊに対する源信号成分ｘｊ（ｋ，ｆ）の和として、（３１７）式によりモデル化される。
【０５６５】
【数２３６】
2019163487.tif 000238
【０５６６】
例えば、ある時間周波数点（ｋ，ｆ）でＢ（ｋ，ｆ）＝｛１，３｝である場合、（ｋ，ｆ）では、（３１７）式は、（３１８）式のようになる。
【０５６７】
【数２３７】
2019163487.tif 000239
【０５６８】
Ｂ（ｋ，ｆ）が取りうる値（集合）の全体を、集合｛１，・・・，Ｌ｝の空でない部分集合からなる集合系（すなわち、冪集合２｛１，・・・，Ｌ｝の部分集合であって空集合を元として含まないもの）Ｈで表す。例えば、Ｂ（ｋ，ｆ）が＃Ｂ≦２なる任意の空でない集合Ｂ⊂｛１，・・・，Ｌ｝になり得る場合、Ｈは、（３１９）式及び（３２０）式で与えられる。
【０５６９】
【数２３８】
2019163487.tif 000240
【０５７０】
（３１９）式及び（３２０）式の集合系ＨをＨｍａｘで表す。本第７の実施形態では、集合系Ｈは、あらかじめ与えられており、Ｈ⊂Ｈｍａｘであるとする。また、＃Ｂ＝２なる集合Ｂ∈Ｈが少なくとも一つはあるものとする。
【０５７１】
すなわち、集合系Ｈの元（集合）Ｂが存在して、＃Ｂ＝２であるとする。
【０５７２】
信号パラメータ分析部７２０は、空間共分散行列更新部７９０からの空間共分散行列Ｒｊ（ｆ）の推定値＾Ｒｊ（ｆ）（ただし、例外として、信号パラメータ分析部７２０での最初の処理の際には、初期化部からの空間共分散行列Ｒｊ（ｆ）の初期推定値＾Ｒｊ（０）（ｆ））を受け取って、Ｌ個の源信号成分のうちの任意の２個ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）（ｊ，ｌは１以上Ｌ以下の相異なる任意の自然数）に対して、同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）と対角行列Λｊ，｛ｊ，ｌ｝（ｆ）を計算する（ステップＳ１０６）。ここで、同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）は２個の源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）を同時無相関化する行列（すなわち、２個の源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）の空間共分散行列を同時対角化する行列）であり、添え字である集合｛ｊ，ｌ｝は、同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）によって同時無相関化される源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）のインデックスを表している。また、対角行列Λｊ，｛ｊ，ｌ｝（ｆ）は、源信号成分ｘｊ（ｋ，ｆ）の空間共分散行列を同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）によって対角化した際の対角行列である。
【０５７３】
具体的には、信号パラメータ分析部７２０は以下の処理を行う。まず、２個の源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）（ｊ，ｌは、１≦ｊ＜ｌ≦Ｌなる任意の２つの自然数）の空間共分散行列の推定値に対する一般化固有値問題＾Ｒｌ（ｆ）ｐ＝λ＾Ｒｊ（ｆ）ｐを解くことにより、（ただし、例外として、信号パラメータ分析部７２０での最初の処理の際には、これらの源信号成分の空間共分散行列の初期推定値に対する一般化固有値問題＾Ｒｌ（０）（ｆ）ｐ＝λ＾Ｒｊ（０）（ｆ）ｐを解くことにより、）一般化固有値λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）と、一般化固有値λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）にそれぞれ対応する一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）と、を求める。ただし、一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）は、（３２１）式を満たすように取るものとする。
【０５７４】
【数２３９】
2019163487.tif 000241
【０５７５】
なお、信号パラメータ分析部７２０での最初の処理の際には、（３２１）式において、＾Ｒｊ（ｆ）を＾Ｒｊ（０）（ｆ）で置き換える。信号パラメータ分析部７２０は、次に、一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）を列ベクトルとする行列を同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）とし、Ｉ次単位行列を対角行列Λｊ，｛ｊ，ｌ｝（ｆ）とし、一般化固有値λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）を対角要素とする対角行列を対角行列Λｌ，｛ｊ，ｌ｝（ｆ）とする。
【０５７６】
同時無相関化部７３０は、観測信号ベクトル作成部７１０からの観測信号ベクトルｙ（ｋ，ｆ）と、信号パラメータ分析部７２０からの同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）と、を受け取って、源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）（ｊ，ｌは、１≦ｊ＜ｌ≦Ｌなる任意の２つの自然数）について同時無相関化された観測信号ベクトルｙ´｛ｊ，ｌ｝（ｋ，ｆ）を次の（３２２）式により計算する（ステップＳ１０７）。
【０５７７】
【数２４０】
2019163487.tif 000242
【０５７８】
ここで、添え字である集合｛ｊ，ｌ｝は、同時無相関化される源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）のインデックスを表す。
【０５７９】
信号存在事後確率更新部７４０は、信号パラメータ分析部７２０からの同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）及び対角行列Λｊ，｛ｊ，ｌ｝（ｆ）と、同時無相関化部７３０からの同時無相関化された観測信号ベクトルｙ´｛ｊ，ｌ｝（ｋ，ｆ）と、混合重み更新部７７０からの混合重みαＢ（ｆ）の推定値＾αＢ（ｆ）（ただし、例外として、信号存在事後確率更新部７４０における最初の処理の際には、初期化部からの混合重み＾αＢ（ｆ）の初期推定値＾αＢ（０）（ｆ））と、分散パラメータ更新部７８０からの分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）（ただし、例外として、信号存在事後確率更新部７４０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ，ｆ）の初期推定値＾ｖｊ（０）（ｋ，ｆ））と、を受け取って、集合系Ｈの各元（集合）Ｂに対して、信号存在事後確率γＢ（ｋ，ｆ）を（３２３）式により更新する（ステップＳ１０８）。
【０５８０】
【数２４１】
2019163487.tif 000243
【０５８１】
なお、集合系Ｈの定義によっては、（３２３）式において、ただ一つの元からなる集合（一元集合）｛ｊ｝に対する同時無相関化行列Ｐ｛ｊ｝（ｆ）、対角行列Λｊ，｛ｊ｝（ｆ）及び同時無相関化された観測信号ベクトルｙ´｛ｊ｝（ｋ，ｆ）が現れることがあるが、これらは二元集合に対するＰ｛ｊ，ｌ｝（ｆ）、Λｊ，｛ｊ，ｌ｝（ｆ）及びｙ´｛ｊ，ｌ｝（ｋ，ｆ）を用いて、（３２４）式、（３２５）式及び（３２６）式により定義される。
【０５８２】
【数２４２】
2019163487.tif 000244
【数２４３】
2019163487.tif 000245
【数２４４】
2019163487.tif 000246
【０５８３】
ここで、ρ：｛１，．．．，Ｌ｝→｛１，．．．，Ｌ｝は、「すべてのｊ∈｛１，．．．，Ｌ｝に対してρ（ｊ）≠ｊ」なる所定の関数である。このように定義する理由は、次の通りである。
【０５８４】
すなわち、１個の源信号成分ｘｊ（ｋ，ｆ）の同時無相関化とは、単なるｘｊ（ｋ，ｆ）の無相関化に他ならず、それは２個の源信号成分ｘｊ（ｋ，ｆ），ｘρ（ｊ）（ｋ，ｆ）の同時無相関化により当然実現されるからである。
【０５８５】
シングルチャネルウィーナーフィルタ部７５０は、同時無相関化部７３０からの同時無相関化された観測信号ベクトルｙ´｛ｊ，ｌ｝（ｋ，ｆ）と、信号パラメータ分析部７２０からの対角行列Λｊ，｛ｊ，ｌ｝（ｆ）と、分散パラメータ更新部７８０からの分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）（ただし、例外として、シングルチャネルウィーナーフィルタ部７５０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ，ｆ）の初期推定値＾ｖｊ（０）（ｋ，ｆ））と、を受け取って、同時無相関化された事後平均ｍ´ｊ，｛ｊ，ｌ｝（ｋ，ｆ）（ｊ，ｌは１以上Ｌ以下の相異なる任意の自然数）を（３２７）式により更新する（ステップＳ１０９）。
【０５８６】
【数２４５】
2019163487.tif 000247
【０５８７】
事後共分散行列更新部７６０は、信号パラメータ分析部７２０からの対角行列Λｊ，｛ｊ，ｌ｝（ｆ）と、分散パラメータ更新部７８０からの分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）（ただし、例外として、事後共分散行列更新部７６０における最初の処理の際には、初期化部からの分散パラメータｖｊ（ｋ，ｆ）の初期推定値＾ｖｊ（０）（ｋ，ｆ））と、を受け取って、同時無相関化された事後共分散行列Σ´｛ｊ，ｌ｝（ｋ，ｆ）（ｊ，ｌは１≦ｊ＜ｌ≦Ｌなる任意の２つの自然数）を次の（３２８）式により更新する（ステップＳ１１０）。
【０５８８】
【数２４６】
2019163487.tif 000248
【０５８９】
混合重み更新部７７０は、信号存在事後確率更新部７４０からの信号存在事後確率γＢ（ｋ，ｆ）を受け取って、（３２９）式により混合重みαＢ（ｆ）の推定値＾αＢ（ｆ）（Ｂは集合系Ｈの任意の元（集合））を更新する（ステップＳ１１１）。
【０５９０】
【数２４７】
2019163487.tif 000249
【０５９１】
分散パラメータ更新部７８０は、信号パラメータ分析部７２０からの対角行列Λｊ，｛ｊ，ｌ｝（ｆ）と、信号存在事後確率更新部７４０からの信号存在事後確率γＢ（ｋ，ｆ）と、シングルチャネルウィーナーフィルタ部７５０からの同時無相関化された事後平均ｍ´ｊ，｛ｊ，ｌ｝（ｋ，ｆ）と、事後共分散行列更新部７６０からの同時無相関化された事後共分散行列Σ´｛ｊ，ｌ｝（ｋ，ｆ）とを受け取って、分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）を（３３０）式により更新する（ステップＳ１１２）。
【０５９２】
【数２４８】
2019163487.tif 000250
【０５９３】
なお、集合系Ｈの定義によっては、（３３０）式において、一元集合｛ｊ｝に対する同時無相関化された事後平均ｍ´ｊ，｛ｊ｝（ｋ，ｆ）及び同時無相関化された事後共分散行列Σ´｛ｊ｝（ｋ，ｆ）が現れることがあるが、これらは（３３１）式及び（３３２）式により定義される。
【０５９４】
【数２４９】
2019163487.tif 000251
【数２５０】
2019163487.tif 000252
【０５９５】
空間共分散行列更新部７９０は、信号パラメータ分析部７２０からの同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）と、信号存在事後確率更新部７４０からの信号存在事後確率γＢ（ｋ，ｆ）と、シングルチャネルウィーナーフィルタ部７５０からの同時無相関化された事後平均ｍ´ｊ，｛ｊ，ｌ｝（ｋ，ｆ）と、事後共分散行列更新部７６０からの同時無相関化された事後共分散行列Σ´｛ｊ，ｌ｝（ｋ，ｆ）と、分散パラメータ更新部７８０からの分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）と、を受け取って、空間共分散行列Ｒｊ（ｆ）の推定値＾Ｒｊ（ｆ）を（３３３）式により更新する（ステップＳ１１３）。
【０５９６】
【数２５１】
2019163487.tif 000253
【０５９７】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ１１４）。収束判定部が収束していないと判定した場合（ステップＳ１１４：Ｎｏ）、信号パラメータ分析部７２０での処理に戻って、処理が継続される。
【０５９８】
一方、収束判定部が収束したと判定した場合（ステップＳ１１４：Ｙｅｓ）、源信号成分推定部８００での処理が行われる。源信号成分推定部８００は、信号パラメータ分析部７２０からの同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）と、信号存在事後確率更新部７４０からの信号存在事後確率γＢ（ｋ，ｆ）と、シングルチャネルウィーナーフィルタ部７５０からの同時無相関化された事後平均ｍ´ｊ，｛ｊ，ｌ｝（ｋ，ｆ）と、を受け取って、源信号成分ｘｊ（ｋ，ｆ）の推定値＾ｘｊ（ｋ，ｆ）を（３３４）式により推定し（ステップＳ１１５）、出力する。
【０５９９】
【数２５２】
2019163487.tif 000254
【０６００】
周波数領域信号分離部３２０は、ｆ＝Ｆであるかどうか判定する（ステップＳ１１６）。ｆ＝Ｆでないと判定された場合（ステップＳ１１６：Ｎｏ）、周波数領域信号分離部３２０でのｆ←ｆ＋１の処理に戻って、処理が継続される。
【０６０１】
一方、ｆ＝Ｆであると判定された場合（ステップＳ１１６：Ｙｅｓ）、周波数領域信号分離部３２０での処理は終了する。
【０６０２】
［信号分析装置の処理の背景］
次に、本第７の実施形態に係る信号分析装置の処理の背景について説明する。本第７の実施形態では、観測信号ベクトルｙ（ｋ，ｆ）は、（３１７）式でモデル化される。
【０６０３】
本第７の実施形態では、源信号成分ｘｊ（ｋ，ｆ）が平均０、共分散行列Ｓｊ（ｋ，ｆ）＝ｖｊ（ｋ，ｆ）Ｒｊ（ｆ）の複素ガウス分布に従うと仮定する。このとき、Ｂ（ｋ，ｆ）が値Ｂ∈Ｈを取るという条件下での観測信号ベクトルｙ（ｋ，ｆ）の確率分布は、次の（３３５）式で与えられる。
【０６０４】
【数２５３】
2019163487.tif 000255
【０６０５】
ここで、Φ１はすべてのモデルパラメータをまとめて表したものである。また、集合Ｂ（ｋ，ｆ）の確率分布は次の（３３６）式で与えられる。
【０６０６】
【数２５４】
2019163487.tif 000256
【０６０７】
ただし、αＢ（ｆ）は、次の（３３７）式の制約条件を満たす。
【０６０８】
【数２５５】
2019163487.tif 000257
【０６０９】
（３３８）式及び（３３６）式より、観測信号ベクトルｙ（ｋ，ｆ）の確率分布は、次の（３３８）式で与えられる。
【０６１０】
【数２５６】
2019163487.tif 000258
【０６１１】
以下、αＢ（ｆ）を混合重みと呼ぶ。パラメータΦ１は、混合重みαＢ（ｆ）、分散パラメータｖｊ（ｋ，ｆ）及び空間共分散行列Ｒｊ（ｆ）からなる。パラメータΦ１は、次の（３３９）式の尤度関数の最大化により推定できる。
【０６１２】
【数２５７】
2019163487.tif 000259
【０６１３】
参考文献６「R.  Ikeshita,  M.  Togami,  Y.  Kawaguchi,  Y.  Fujita,  and  K.  Nagamatsu,  “Local  Gaussian  model  with  source-set  constraints  in  audio  source  separation”,  IEEE  International  Workshop  on  Machine  Learning  for  Signal  Processing,  Sept.  2017.」には、ＥＭアルゴリズムに基づいて上式の尤度関数を最大化することにより、パラメータΦ１を推定する方法が開示されている。この従来技術では、ＥステップとＭステップを交互に反復することにより、パラメータΦ１を推定する。
【０６１４】
Ｅステップでは、信号存在事後確率γＢ（ｋ，ｆ）、事後平均ｍｊ，Ｂ（ｋ，ｆ）及び事後共分散行列ΣＢ（ｋ，ｆ）を次の（３４０）式〜（３４４）式により計算する。
【０６１５】
【数２５８】
2019163487.tif 000260
【数２５９】
2019163487.tif 000261
【数２６０】
2019163487.tif 000262
【数２６１】
2019163487.tif 000263
【数２６２】
2019163487.tif 000264
【０６１６】
ここで、信号存在事後確率γＢ（ｋ，ｆ）は、観測信号ベクトルｙ（ｋ，ｆ）が与えられたという条件下での集合Ｂ（ｋ，ｆ）の条件付き確率分布である事後分布として、次の（３４５）式で定義される。
【０６１７】
【数２６３】
2019163487.tif 000265
【０６１８】
ここで、＾Φ１はパラメータΦ１の現在の推定値を表す。また、事後平均ｍｊ，Ｂ（ｋ，ｆ）及び事後共分散行列ΣＢ（ｋ，ｆ）とは、観測信号ベクトルｙ（ｋ，ｆ）と集合Ｂ（ｋ，ｆ）が与えられたという条件下での源信号成分ｘｊ（ｋ，ｆ）の条件付き確率分布である事後分布（複素ガウス分布であることが示される）の平均及び共分散行列として定義される。この事後分布は次の（３４６）式のように表される。
【０６１９】
【数２６４】
2019163487.tif 000266
【０６２０】
Ｍステップでは、混合重みαＢ（ｆ）の推定値＾αＢ（ｆ）、分散パラメータｖｊ（ｋ，ｆ）の推定値＾ｖｊ（ｋ，ｆ）及び空間共分散行列Ｒｊ（ｆ）の推定値＾Ｒｊ（ｆ）を、次の（３４７）式〜（３４９）式の更新則により更新する。
【０６２１】
【数２６５】
2019163487.tif 000267
【数２６６】
2019163487.tif 000268
【数２６７】
2019163487.tif 000269
【０６２２】
しかしながら、この従来技術では、（３４０）式の信号存在事後確率γＢ（ｋ，ｆ）の更新、（３４２）式の事後平均ｍｊ，Ｂ（ｋ，ｆ）の更新、及び（３４４）式の事後共分散行列Σ｛ｊ，ｌ｝（ｋ，ｆ）の更新において、行列式計算、逆行列計算及び行列乗算を時間周波数点ごとに行う必要があった。具体的には、（３４０）式における（３５０）式の形の行列式計算と、（３４０）式及び（３４２）式及び（３４４）式における（３５１）式の形の逆行列計算と、（３４４）式における（３５２）式の形の行列乗算と、を時間周波数点（ｋ，ｆ）ごとに行う必要があった。
【０６２３】
【数２６８】
2019163487.tif 000270
【数２６９】
2019163487.tif 000271
【数２７０】
2019163487.tif 000272
【０６２４】
行列式計算、逆行列計算及び行列乗算はいずれもＯ（Ｉ３）のオーダーの計算量を必要とする行列演算であり、かつ時間周波数点の数は通常、非常に大きいため、この従来技術では膨大な計算量を要していた。
【０６２５】
いま仮に、（３５０）式〜（３５２）式における共分散行列＾ｖｊ（ｋ，ｆ）＾Ｒｊ（ｆ）及び＾ｖｌ（ｋ，ｆ）＾Ｒｌ（ｆ）がともに対角行列であるとすると、（３５０）式〜（３５２）式における行列式計算、逆行列計算及び行列乗算は単なる対角要素に対するスカラー演算に帰着するため、少ない計算量で計算できる。しかしながら、通常、源信号成分（ベクトル）の異なる要素は互いに相関を持つため、共分散行列＾ｖｊ（ｋ，ｆ）＾Ｒｊ（ｆ）及び＾ｖｌ（ｋ，ｆ）＾Ｒｌ（ｆ）は対角行列ではない。
【０６２６】
そこで、本第７の実施形態では、一般化固有値問題に基づいて２つの空間共分散行列＾Ｒｊ（ｆ），＾Ｒｌ（ｆ）を同時対角化することを考える。λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）を空間共分散行列＾Ｒｌ（ｆ），＾Ｒｊ（ｆ）に対する一般化固有値問題（（３５３）式）の一般化固有値とし、ｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）を一般化固有値λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）に対応する一般化固有ベクトルとする。
【０６２７】
【数２７１】
2019163487.tif 000273
【０６２８】
ただし、一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）は、次の（３５４）式を満たすように取る。
【０６２９】
【数２７２】
2019163487.tif 000274
【０６３０】
ｐ｛ｊ，ｌ｝（ｆ）を一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）を列ベクトルとする行列、Λ｛ｊ，ｌ｝（ｆ）を一般化固有値λ１，｛ｊ，ｌ｝（ｆ），・・・，λＩ，｛ｊ，ｌ｝（ｆ）を対角要素とする対角行列とすると、次の（３５５）式が成り立つ。
【０６３１】
【数２７３】
2019163487.tif 000275
【０６３２】
（３５５）式は、行列Ｐ｛ｊ，ｌ｝（ｆ）が２つの空間共分散行列＾Ｒｊ（ｆ），＾Ｒｌ（ｆ）を同時対角化する（２つの源信号成分ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）を同時無相関化する）ことを示している。そこで、行列Ｐ｛ｊ，ｌ｝（ｆ）を同時無相関化行列と呼ぶ。
【０６３３】
（３５５）式を空間共分散行列＾Ｒｊ（ｆ），＾Ｒｌ（ｆ）について解いて、（３４０）式、（３４２）式及び（３４４）式に代入し、（３５６）式〜（３５８）式と定義する（ｊ，ｌは１以上Ｌ以下の互いに異なる任意の整数）と、（３６２）、（３６３）式及び（３６４）式（後述）を得る。
【０６３４】
【数２７４】
2019163487.tif 000276
【数２７５】
2019163487.tif 000277
【数２７６】
2019163487.tif 000278
【０６３５】
（３５６）式及び（３５７）式を事後平均ｍｊ，｛ｊ，ｌ｝（ｋ，ｆ）及び事後共分散行列Σ｛ｊ，ｌ｝（ｋ，ｆ）について解き、（３４８）式及び（３４９）式に代入すると（３６６）式及び（３６８）式を得る。
【０６３６】
これで、本第７の実施形態に係る信号分析装置における更新則が導出された。
【０６３７】
このように、本第７の実施形態に係る信号分析装置では、信号パラメータ分析部７２０は、２個の混在する源信号それぞれの空間的特性をモデル化する空間共分散行列＾Ｒｌ（ｆ），＾Ｒｊ（ｆ）に対するｐｉ，｛ｊ，ｌ｝（ｆ）Ｈ＾Ｒｊ（ｆ）ｐｍ，｛ｊ，ｌ｝（ｆ）＝０（ｉ，ｍは１以上Ｉ以下の互いに異なる任意の整数）を満たす、それぞれ異なるＩ個（Ｉは２以上の整数）の一般化固有ベクトルｐ１，｛ｊ，ｌ｝（ｆ），・・・，ｐＩ，｛ｊ，ｌ｝（ｆ）を列ベクトルとする行列である同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）またはそのエルミート転置Ｐ｛ｊ，ｌ｝（ｆ）Ｈを、Ｉ個の位置で取得される観測信号による観測信号ベクトルについて２個の源信号に対応する成分を無相関化するためのパラメータとして得る。
【０６３８】
［第７の実施形態の実施例］
本実施例では、集合系Ｈの具体例について説明する。本実施例では、集合系Ｈを（３５９）式とする。
【０６３９】
【数２７７】
2019163487.tif 000279
【０６４０】
これは、各時間周波数点において観測信号ベクトルは、１番目の源信号成分をつねに含み、２番目〜Ｊ番目の源信号成分のうちのちょうど１つを含むことを意味する。例えば、１番目の源信号成分は背景雑音のようなスパースでない源信号成分をモデル化し、２番目〜Ｊ番目の源信号成分は音声のようなスパースな源信号成分をモデル化する。
【０６４１】
信号パラメータ分析部７２０、同時無相関化部７３０、信号存在事後確率更新部７４０、シングルチャネルウィーナーフィルタ部７５０、事後共分散行列更新部７６０、混合重み更新部７７０、分散パラメータ更新部７８０、空間共分散行列更新部７９０及び源信号成分推定部８００における処理の更新式を（３６０）式〜（３７０）式に示す。本第７の実施形態では、Ｈの元である集合｛１，２｝，｛１，３｝，・・・，｛１，Ｊ｝は、それぞれインデックス２，３，・・・，Ｊと１対１に対応しているから、簡単のため以下の記法を用いる。すなわち、ｙ´｛１，ｊ｝（ｋ，ｆ）をｙ´ｊ（ｋ，ｆ）と書き、Ｐ｛１，ｊ｝（ｆ）をＰｊ（ｆ）と書き、γ｛１，ｊ｝（ｋ，ｆ）をγｊ（ｋ，ｆ）と書き、＾α｛１，ｊ｝（ｆ）を＾αｊ（ｆ）と書き、Λｌ，｛１，ｊ｝（ｆ）をΛｌｊ（ｆ）と書き、ｍ´ｌ，｛１，ｊ｝（ｋ，ｆ）をｍ´ｌｊ（ｋ，ｆ）と書き、Σ´｛１，ｊ｝（ｋ，ｆ）をΣ´ｊ（ｋ，ｆ）と書く。ただし、ｊ＝２，・・・，Ｊとし、ｌ＝１，ｊとする（下記の更新式でも同様）。なお、本実施例では、信号パラメータ分析部７２０は、Ｌ個の源信号成分のうちの任意の２個ｘｊ（ｋ，ｆ），ｘｌ（ｋ，ｆ）（ｊ，ｌは１以上Ｌ以下の相異なる任意の自然数）に対して、同時無相関化行列Ｐ｛ｊ，ｌ｝（ｆ）と対角行列Λｊ，｛ｊ，ｌ｝（ｆ）を計算する必要はなく、（３６０）式にあるように、ｘ１（ｋ，ｆ），ｘｊ（ｋ，ｆ）（ｊは２以上Ｌ以下の任意の自然数）に対して、同時無相関化行列Ｐ｛１，ｊ｝（ｆ）（すなわちＰｊ（ｆ））と対角行列Λｌ，｛１，ｊ｝（ｆ）（すなわちΛｌｊ（ｆ））を計算すれば十分である。
【０６４２】
【数２７８】
2019163487.tif 000280
【０６４３】
【数２７９】
2019163487.tif 000281
【０６４４】
【数２８０】
2019163487.tif 000282
【０６４５】
【数２８１】
2019163487.tif 000283
【０６４６】
【数２８２】
2019163487.tif 000284
【０６４７】
【数２８３】
2019163487.tif 000285
【０６４８】
【数２８４】
2019163487.tif 000286
【０６４９】
【数２８５】
2019163487.tif 000287
【０６５０】
【数２８６】
2019163487.tif 000288
【０６５１】
【数２８７】
2019163487.tif 000289
【０６５２】
【数２８８】
2019163487.tif 000290
【０６５３】
［第８の実施形態］
本第８の実施形態では、第３の実施形態に係る信号分析装置の他の構成例について説明する。図２１は、第８の実施形態に係る信号分析装置の構成の一例を示す図である。図２２は、第８の実施形態に係る信号分析処理の処理手順の一例を示すフローチャートである。図２１に示すように、第８の実施の形態に係る信号分析装置８０１は、周波数領域変換部３１０、周波数領域信号分離部３２０−１及び時間領域変換部３３０を有する。
【０６５４】
信号分析装置８０１では、周波数領域信号分離部３２０−１において、スペクトル基底ｔｎ（ｆ）と、アクティベーション基底ｗｎ（ｋ）と、分配関数ｚｎｊと、同時無相関化された空間共分散行列Λｊ（ｆ）と、同時無相関化行列Ｐ（ｆ）との推定値とを計算する。
【０６５５】
スペクトル基底ｔｎ（ｆ）と、アクティベーション基底ｗｎ（ｋ）と、分配関数ｚｎｊとは、分散パラメータｖｊ（ｋ，ｆ）をモデル化するパラメータである。分散パラメータｖｊ（ｋ，ｆ）は、これらのパラメータを用いて（３７１）式によりモデル化される。
【０６５６】
【数２８９】
2019163487.tif 000291
【０６５７】
ここで、ｎは、スペクトル基底およびアクティベーション基底のインデックスである。Ｎは、スペクトル基底およびアクティベーション基底の個数である。
【０６５８】
周波数領域信号分離部３２０−１は、初期化部（図示しない）、観測信号ベクトル作成部８１０、同時無相関化８３０、シングルチャネルウィーナーフィルタ部８４０、観測共分散行列更新部８５０、分散パラメータ更新部８６０、空間共分散行列更新部８７０、信号パラメータ分析部８２０、収束判定部（図示しない）及び同時無相関化逆変換部８８０を有する。
【０６５９】
［信号分析装置の構成及び処理］
図２１及び図２２を用いて、第８の実施形態に係る信号分析装置の構成及び信号分析処理の処理手順について説明する。
【０６６０】
図示しない初期化部は、分散パラメータｖｊ（ｋ，ｆ）と、スペクトル基底ｔｎ（ｆ）と、アクティベーション基底ｗｎ（ｋ）と、分配関数ｚｎｊと、同時無相関化された空間共分散行列Λｊ（ｆ）と、同時無相関化行列Ｐ（ｆ）と、の初期推定値を計算する初期化を行う（ステップＳ１２１）。初期推定値は、上付きの（０）を付して、^ｖｊ（０）（ｋ，ｆ）などと表す。
【０６６１】
周波数領域変換部３１０は、観測信号を取得し（ステップＳ１２２）、観測信号を周波数領域に変換する（ステップＳ１２３）ことにより、時間周波数領域の観測信号ｙｉ（ｋ，ｆ）を計算する。
【０６６２】
観測信号ベクトル作成部８１０は、観測信号ベクトルｙ（ｋ，ｆ）を第７の実施形態と同様にして作成する（ステップＳ１２４）。
【０６６３】
同時無相関化部８３０は、同時無相関化された観測信号ベクトルｙ´（ｋ，ｆ）を（３７２）式により計算し、更新する（ステップＳ１２５）。
【０６６４】
【数２９０】
2019163487.tif 000292
【０６６５】
ここで、＾Ｐ（ｆ）は信号パラメータ分析部８２０からの同時無相関化行列の推定値である。
【０６６６】
観測共分散行列更新部８５０は、同時無相関化された観測共分散行列Ｕ´（ｋ，ｆ）を（３７３）式により更新する（ステップＳ１２６）。
【０６６７】
【数２９１】
2019163487.tif 000293
【０６６８】
ここで、＾ｖｊ（ｋ，ｆ）は、分散パラメータ更新部８６０からの分散パラメータの推定値である。＾Λｊ（ｆ）は、空間共分散行列更新部８７０からの同時無相関化された空間共分散行列の推定値である。
【０６６９】
分散パラメータ更新部８６０は、次の手順により分散パラメータの推定値＾ｖｊ（ｋ，ｆ）を更新する（ステップＳ１２７）。
【０６７０】
まず、分散パラメータ更新部８６０は、ｎ番目の基底に対応する同時無相関化された空間共分散行列Σｎ（ｆ）を（３７４）式により更新する。（３７４）式において、＾ｚｎｊは、分配関数の推定値である。
【０６７１】
【数２９２】
2019163487.tif 000294
【０６７２】
次に、分散パラメータ更新部はスペクトル基底の推定値＾ｔｎ（ｆ）を（３７５）式により更新する。（３７５）式において、＾ｗｎ（ｋ）はアクティベーション基底の推定値である。
【０６７３】
【数２９３】
2019163487.tif 000295
【０６７４】
次に、分散パラメータ更新部８６０は、アクティベーション基底の推定値＾ｗｎ（ｋ）を（３７６）式により更新する。
【０６７５】
【数２９４】
2019163487.tif 000296
【０６７６】
次に、分散パラメータ更新部８６０は、分配関数の推定値＾ｚｎｊを（３７７）式により更新する。
【０６７７】
【数２９５】
2019163487.tif 000297
【０６７８】
最後に、分散パラメータ更新部８６０は、分散パラメータの推定値＾ｖｊ（ｋ，ｆ）を（３７８）式により更新する。
【０６７９】
【数２９６】
2019163487.tif 000298
【０６８０】
空間共分散行列更新部８７０は、（３７９）式により、同時無相関化された空間共分散行列の推定値＾Λｊ（ｆ）を更新する（ステップＳ１２８）。
【０６８１】
【数２９７】
2019163487.tif 000299
【０６８２】
ここで、＃は、行列の幾何平均であり、Ａｊ（ｆ）およびＢｊ（ｆ）は、（３８０）式および（３８１）式で計算される行列である。
【０６８３】
【数２９８】
2019163487.tif 000300
【０６８４】
【数２９９】
2019163487.tif 000301
【０６８５】
信号パラメータ分析部８２０は、同時無相関化行列の推定値＾Ｐ（ｆ）を（３８２）式により更新する（ステップＳ１２９）。
【０６８６】
【数３００】
2019163487.tif 000302
【０６８７】
図示しない収束判定部は、収束したかどうかの判定を行う（ステップＳ１３０）。収束判定部が収束していないと判定した場合（ステップＳ１３０：Ｎｏ）、同時無相関化部８３０での処理（ステップＳ１２５）に戻って、処理が継続される。
【０６８８】
一方、収束判定部が収束したと判定した場合（ステップＳ１３０：Ｙｅｓ）、シングルチャネルウィーナーフィルタ部８４０は、同時無相関化された源信号成分の推定値＾ｘ´ｊ（ｋ，ｆ）を（３８３）式により計算する（ステップＳ１３１）。
【０６８９】
【数３０１】
2019163487.tif 000303
【０６９０】
シングルチャネルウィーナーフィルタ部８４０は、同時無相関化された源信号成分の推定値＾ｘ´ｊ（ｋ，ｆ）を更新する（ステップＳ１３２）。
【０６９１】
同時無相関化逆変換部８８０は、源信号成分の推定値＾ｘｊ（ｋ，ｆ）を（３８４）式により計算する（ステップＳ１３３）。
【０６９２】
【数３０２】
2019163487.tif 000304
【０６９３】
時間領域変換部３３０は、第７の実施形態と同様の処理を行って、時間領域の源信号成分の推定値＾ｘｊ（ｔ）を計算し、出力する（ステップＳ１３４）。
【０６９４】
［第８の実施形態の効果］
第３の実施形態では、音の方向性のみを基に信号処理を行っていた。これに対して、本第８の実施形態では、音の構造（パーツ）も考慮して信号分析を行うため、音源分離の精度を向上させることができる。
【０６９５】
そして、第８の実施形態では、パーミテーション問題が原理的に起こらないという効果を奏する。本第８の実施形態では、全ての周波数において、どのようなパーツがあるかを考慮した処理を行うため、周波数ごとに独立に処理しないことから、パーミテーション問題が原理的に起こらない。
【０６９６】
［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、ＣＰＵ及び当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【０６９７】
また、本実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。すなわち、上記学習方法及び音声認識方法において説明した処理は、記載の順にしたがって時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。
【０６９８】
［プログラム］
図２３は、プログラムが実行されることにより、信号分析装置１，２０１，３０１，４０１，５０１，６０１が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。
【０６９９】
メモリ１０１０は、ＲＯＭ１０１１及びＲＡＭ１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。
【０７００】
ハードディスクドライブ１０９０は、例えば、ＯＳ（Operating  System）１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、信号分析装置１，２０１，３０１，４０１，５０１，６０１の各処理を規定するプログラムは、コンピュータ１０００により実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、信号分析装置１，２０１，３０１，４０１，５０１，６０１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid  State  Drive）により代替されてもよい。
【０７０１】
また、上述した実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０７０２】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３及びプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３及びプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【０７０３】
以上、本発明者によってなされた発明を適用した実施形態について説明したが、本実施形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。すなわち、本実施形態に基づいて当業者等によりなされる他の実施形態、実施例及び運用技術等はすべて本発明の範疇に含まれる。
【符号の説明】
【０７０４】
１，２０１，３０１，４０１，５０１，６０１，８０１  信号分析装置
１０，２１０，４１０，５１０，６１０，７１０，８１０  観測信号ベクトル作成部
２０，２２０，４２０，５２０，６２０，７２０，８２０  信号パラメータ分析部
３０，２３０，４３０，５３０，６３０，７３０，８３０  同時無相関化部
４０，２４０，５４０，６４０，７５０，８４０  シングルチャネルウィーナーフィルタ部
５０，２８０，５８０，６８０，８８０  同時無相関化逆変換部
２５０，５５０，７６０  事後共分散行列更新部
２６０，４４０，５６０，７８０，８６０  分散パラメータ更新部
２７０，４５０，５７０，７９０，８７０  空間共分散行列更新部
３１０  周波数領域変換部
３２０  周波数領域信号分離部
３３０  時間領域変換部
６５０，８５０  観測共分散行列更新部
７４０  信号存在事後確率更新部
７７０  混合重み更新部
８００  源信号成分推定部
【図１】
2019163487.tif 000305
【図２】
2019163487.tif 000306
【図３】
2019163487.tif 000307
【図４】
2019163487.tif 000308
【図５】
2019163487.tif 000309
【図６】
2019163487.tif 000310
【図７】
2019163487.tif 000311
【図８】
2019163487.tif 000312
【図９】
2019163487.tif 000313
【図１０】
2019163487.tif 000314
【図１１】
2019163487.tif 000315
【図１２】
2019163487.tif 000316
【図１３】
2019163487.tif 000317
【図１４】
2019163487.tif 000318
【図１５】
2019163487.tif 000319
【図１６】
2019163487.tif 000320
【図１７】
2019163487.tif 000321
【図１８】
2019163487.tif 000322
【図１９】
2019163487.tif 000323
【図２０】
2019163487.tif 000324
【図２１】
2019163487.tif 000325
【図２２】
2019163487.tif 000326
【図２３】
2019163487.tif 000327
【図２４】
2019163487.tif 000328
【国際調査報告】
2019163487.tif 000329
2019163487.tif 000330
2019163487.tif 000331
2019163487.tif 000332