(57)【要約】【課題】ビームフォーミングのフィルタ群を導出する際に用いるPSD分離行列を解析的に設計する分離行列設計装置を提供する。【解決手段】分離行列設計装置は、Mを2以上の整数とし、M個の収音信号を入力とし、L個のビームフォーミング出力値を得るビームフォーミングにおいて、l=1,2,…,Lとし、i=1,2,…,Lとし、PSD分離行列GL,ωが、l番目のビームフォーマの方向Θiに対する感度を要素とするL×Lの非負行列であるPSD混合行列DL,ωの逆行列となり、かつ、PSD分離行列GL,ωがM行列となるように、PSD分離行列GL,ωを設計する分離行列設計部を含む。【選択図】図７
2018121108.tif 000002
【特許請求の範囲】
【請求項１】
Mを2以上の整数とし、M個の収音信号を入力とし、L個のビームフォーミング出力値を得るビームフォーミングにおいて、l=1,2,…,Lとし、i=1,2,…,Lとし、PSD分離行列GL,ωが、l番目のビームフォーマの方向Θiに対する感度を要素とするL×Lの非負行列であるPSD混合行列DL,ωの逆行列となり、かつ、前記PSD分離行列GL,ωがM行列となるように、前記PSD分離行列GL,ωを設計する分離行列設計部を含む、
分離行列設計装置。
【請求項２】
請求項１の分離行列設計装置であって、
M個のマイクロホンの空間的配置の中心を通るある回転軸に対して1回転させたときに、M個の前記マイクロホンの空間的配置と同一の空間的配置が複数回現れるようにM個の前記マイクロホンが配置されるものとし、前記分離行列設計部は、前記PSD分離行列GL,ωが巡回性を持つM行列となるように前記PSD分離行列GL,ωを設計する、
分離行列設計装置。
【請求項３】
請求項１または請求項２の分離行列設計装置であって、
AHは行列Aの共役転置、ΛL,ωは正の対角行列、FLはDFT基底行列、ILは単位行列、前記PSD分離行列GL,ωは、
【数７７】
2018121108.tif 000079
であり、
A*はAの複素共役であり、前記PSD分離行列GL,ωのパラメータΨL,ωが以下の条件を満たすようにパラメータΨL,ωを設計するパラメータ設計部を含む、
【数７８】
2018121108.tif 000080
分離行列設計装置。
【請求項４】
請求項３の分離行列設計装置であって、
ユーザによって指定された収音領域Aを受け取り、広さを表すためのパラメータμを求める収音領域指定部を含み、
前記パラメータ設計部は、前記パラメータμが大きいほどパラメータψ1,ωの値が大きくなるように設計する、
分離行列設計装置。
【請求項５】
請求項１から請求項４の何れかの分離行列設計装置が設計した前記PSD分離行列GL,ωの逆行列である前記PSD混合行列DL,ωを用いてビームフォーミングのフィルタ係数WL,ωを算出するフィルタ係数算出部を含み、
ΔL,ωは前記PSD混合行列DL,ωの各要素の平方根として求められる行列、HL,ωはアレイマニフォールドベクトル、ΥΔ,M,ωはΔM,ωのDFT出力値により構成される行列であり、FMはDFT基底行列、ΥΔ,M,ωはΔM,ωのDFT出力値により構成される行列、ΔM,ω,HM,ωはΔL,ω,HL,ωに対してそれぞれアップコンバージョンを行った行列であり、前記フィルタ係数算出部は、前記フィルタ係数Wを、M=Lのとき、
【数７９】
2018121108.tif 000081
として算出し、M>Lのとき、
【数８０】
2018121108.tif 000082
として算出する、
フィルタ係数算出装置。
【請求項６】
Mを2以上の整数とし、M個の収音信号を入力とし、L個のビームフォーミング出力値を得るビームフォーミングにおいて、l=1,2,…,Lとし、i=1,2,…,Lとし、PSD分離行列GL,ωが、l番目のビームフォーマの方向Θiに対する感度を要素とするL×Lの非負行列であるPSD混合行列DL,ωの逆行列となり、かつ、前記PSD分離行列GL,ωがM行列となるように、前記PSD分離行列GL,ωを設計する分離行列設計ステップを含む、
分離行列設計装置が実行する分離行列設計方法。
【請求項７】
請求項６の分離行列設計方法で設計した前記PSD分離行列GL,ωの逆行列である前記PSD混合行列DL,ωを用いてビームフォーミングのフィルタ係数WL,ωを算出するフィルタ係数算出ステップを含み、
ΔL,ωは前記PSD混合行列DL,ωの各要素の平方根として求められる行列、HL,ωはアレイマニフォールドベクトル、ΥΔ,M,ωはΔM,ωのDFT出力値により構成される行列であり、FMはDFT基底行列、ΥΔ,M,ωはΔM,ωのDFT出力値により構成される行列、ΔM,ω,HM,ωはΔL,ω,HL,ωに対してそれぞれアップコンバージョンを行った行列であり、前記フィルタ係数算出ステップは、前記フィルタ係数Wを、M=Lのとき、
【数８１】
2018121108.tif 000083
として算出し、M>Lのとき、
【数８２】
2018121108.tif 000084
として算出する、
フィルタ係数算出装置が実行するフィルタ係数算出方法。
【請求項８】
請求項１から請求項４の何れかの分離行列設計装置、または、請求項５のフィルタ係数算出装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数のマイクロホンで構成された音響アレイ装置を用いて、任意の音源信号や位置といった情報を推定、抽出するビームフォーミング技術に関する。
【背景技術】
【０００２】
ビームフォーミングとウィーナーフィルタリングを組み合わせた音源強調方式について研究を進めている。ウィーナーフィルタを推定するための方式として、局所PSD(Power Spectral Density、パワースペクトル密度)推定法がある。局所PSD推定法では、線形性の逆写像により、複数のビームフォーミングの出力信号のPSD群から音源のPSD群を推定し、音源のPSD群からウィーナーフィルタを構成する(非特許文献１参照)。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】K. Niwa, T. Kawase, K. Kobayashi and Y. Hioka, "PSD estimation in beamspace using property of M-matrix", IWAENC2016, 2016.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
従来技術では、ヒューリスティックにビームフォーミングを設計している。そのため、ウィーナーフィルタを構成する過程において、設計したビームフォーミングの特性によって、音源のPSD群の推定精度が劣化してしまうことが分かっている。
【０００５】
本発明では、ビームフォーミングのフィルタ群を導出する際に用いるPSD分離行列を解析的に設計する分離行列設計装置、解析的に設計されたPSD分離行列を用いて解析的にフィルタ係数を算出するフィルタ係数算出装置、その方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記の課題を解決するために、本発明の一態様によれば、分離行列設計装置は、Mを2以上の整数とし、M個の収音信号を入力とし、L個のビームフォーミング出力値を得るビームフォーミングにおいて、l=1,2,…,Lとし、i=1,2,…,Lとし、PSD分離行列GL,ωが、l番目のビームフォーマの方向Θiに対する感度を要素とするL×Lの非負行列であるPSD混合行列DL,ωの逆行列となり、かつ、PSD分離行列GL,ωがM行列となるように、PSD分離行列GL,ωを設計する分離行列設計部を含む。
【０００７】
上記の課題を解決するために、本発明の他の態様によれば、分離行列設計方法は、Mを2以上の整数とし、M個の収音信号を入力とし、L個のビームフォーミング出力値を得るビームフォーミングにおいて、l=1,2,…,Lとし、i=1,2,…,Lとし、PSD分離行列GL,ωが、l番目のビームフォーマの方向Θiに対する感度を要素とするL×Lの非負行列であるPSD混合行列DL,ωの逆行列となり、かつ、PSD分離行列GL,ωがM行列となるように、PSD分離行列GL,ωを設計する分離行列設計ステップを含む。
【発明の効果】
【０００８】
本発明によれば、音源のPSD群の推定精度が高まるようなビームフォーミングのフィルタ群等を解析的に導出することができるという効果を奏する。
【図面の簡単な説明】
【０００９】
【図１】ビームフォーミングとウィーナーポストフィルターを使った音源強調系の機能ブロック図。
【図２】ビームフォーマと到来方向との関係を示す図。
【図３】音源PSDを推定方法を説明するための図。
【図４】従来のビームフォーミングフィルタ係数群、PSD混合行列、PSD分離行列を推定するフローを示す図。
【図５】等間隔に離散化されたM方向から平面波が到来することを仮定したとき、HM,ωを説明するための図。
【図６】ΔM,ωを多重極展開したときの具体的な例を示す図。
【図７】第一実施形態に係るフィルタ係数算出装置の機能ブロック図。
【図８】第一実施形態に係るフィルタ係数算出装置の処理フローの例を示す図。
【図９】第二実施形態に係る音声強調装置の機能ブロック図。
【図１０】第二実施形態に係る音声強調装置の処理フローの例を示す図。
【図１１】第三実施形態に係る音声強調装置の機能ブロック図。
【図１２】第三実施形態に係る音声強調装置の処理フローの例を示す図。
【図１３】ΨL,ωを変えながら作成したフィルタの指向特性を示す図。
【図１４】ΨL,ωを変えながら作成したフィルタの指向特性を示す図。
【図１５】ΨL,ωを変えながら作成したフィルタの指向特性を示す図。
【図１６】ΨL,ωを変えながら作成したフィルタの指向特性を示す図。
【発明を実施するための形態】
【００１０】
以下、本発明の実施形態について、説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【００１１】
＜局所PSD推定に基づく音源強調＞
非特許文献１を従来研究とし、局所PSD(Power Spectral Density: パワースペクトル密度)推定に基づく音源強調についての原理とその課題についてより詳細に説明する。
【００１２】
＜ビームフォーミングとウィーナーポストフィルターを使った音源強調＞
図１は、ビームフォーミングとウィーナーポストフィルターを使った音源強調系の機能ブロック図を示す。ビームフォーミング部９１は、M個の観測信号Xm,ω,τを用いて、ビームフォーミングによりビームフォーミング出力信号Yl,ω,τを求める。ここで、ωとτはそれぞれ周波数、フレーム時間のインデックスを表す。ウィーナーフィルタ推定部９２は、M個のビームフォーミング出力信号Yl,ω,τを用いて、PSD推定に基づきウィーナーフィルタを推定する。ウィーナーフィルタリング部９３は、ビームフォーミング出力信号Yl,ω,τに対してウィーナーフィルタを用いてフィルタリングを行い、雑音成分を抑圧した出力信号Zl,ω,τを得る。
【００１３】
まず、観測信号とビームフォーミング出力信号のモデルについて説明する。次に、局所PSD推定に基づくウィーナーフィルタの推定について説明する。
【００１４】
(ビームフォーミング出力信号のモデル化)
M(≧2)本のマイクロホンを用いて観測することを想定する。音場には、K個の音源がある。k番目の音源信号(時間周波数成分)をSk,ω,τと表す。k番目の音源とm番目のマイクロホンの間の伝達特性をAm,k,ωと表すとき、m番目の観測信号Xm,ω,τは以下でモデル化される。
【００１５】
【数１】
2018121108.tif 000003
【００１６】
ビームフォーミングは、指向性を制御するためによく使われる技術である。l番目のビームフォーマのm番目のマイクロホンに対する係数をWl,m,ωと書くとき、l番目のビームフォーマの出力信号Yl,ω,τは以下で計算される。
【００１７】
【数２】
2018121108.tif 000004
【００１８】
ここで、Dl,k,ωと＊は、l番目のビームフォーマのk番目の音源に対する感度と複素共役演算をそれぞれ表す。
【００１９】
(局所PSD推定に基づくウィーナーフィルタ推定)
音源信号が互いに無相関であることを仮定すると、l番目のビームフォーマの出力信号のPSDは以下でモデル化される。
【００２０】
【数３】
2018121108.tif 000005
【００２１】
ここで、<・>とφS_k,ωは、それぞれ期待値演算子(例えば、時間平均化処理で置き換えることができる)、k番目の音源信号のPSDを表す。なお、下付き添え字S_k,ωはSk,ωを意味する。式(3)の関係が、L個の方向(ある程度の角度幅をもった領域){Θ1,…,ΘL}において音源群をグループ化した際にも成立するのであれば、φY_l,ω(下付き添え字Y_l,ωはYl,ωを意味する）は以下で近似して表される。
【００２２】
【数４】
2018121108.tif 000006
【００２３】
ここで、Dl,Θ_i,ω(下付き添え字l,Θ_i,ωはl,Θi,ωを意味する）はl番目のビームフォーマの方向Θiに対する感度であり、φS_Θ_i,ω(下付き添え字S_Θ_i,ωはSΘ_i,ωを意味し、この中の下付き添え字Θ_iはΘiを意味する）は方向Θiを中心としたある程度の角度幅における音源群の加算値のPSD(音源PSDと以後、呼ぶ)である。φS_Θ_i,ω、φY_l,ωの関係は以下で表される。
【００２４】
【数５】
2018121108.tif 000007
【００２５】
行列DL,ωのことを以後、PSD混合行列と呼ぶことにする。PSD混合行列DL,ωは、対角成分が0より大きく、非対角成分が0以上のL×Lの非負行列である。
【００２６】
音源PSD(ΦS_ω(下付き添え字S_ωはSωを意味する）)を求めるために、式(5)の逆問題を解く必要がある。局所PSD推定法では、ΦS_ωは以下で推定される。
【００２７】
【数６】
2018121108.tif 000008
【００２８】
以後、GL,ωをPSD分離行列と呼ぶことにする。PSD分離行列の構造は以下でモデル化される。
【００２９】
【数７】
2018121108.tif 000009
【００３０】
つまり、GL,ω・DL,ω＝ILであり、PSD分離行列GL,ωは、対角成分が0より大きく、非対角成分が0以下のL×Lの行列である。実用上、効果のある方法として、音源PSDを時間フレームごとに求めることもできる。これは、音源の時間スパース性が高いことを想定しており、以下で実行される。
【００３１】
【数８】
2018121108.tif 000010
【００３２】
l番目の音源群を強調するためのウィーナーポストフィルタは以下で計算される。
【００３３】
【数９】
2018121108.tif 000011
【００３４】
これをビームフォーミング出力信号Yl,ω,τに畳み込むことにより、出力信号Zl,ω,τを得る。
【００３５】
【数１０】
2018121108.tif 000012
【００３６】
(音源PSDを高精度に推定するためにPSD分離行列に求められる条件)
音源PSDを高精度に推定するためにPSD分離行列GL,ωに求められる条件について、まず仮説を立てる。次に、その仮説を満たすためにPSD分離行列GL,ωに求められる条件について明確化する。
【００３７】
(仮説)
図２に示すように、l番目のビームフォーマは、方向Θlから到来する音源を強調すると想定する。式(6)と(7)より、l番目の音源PSDを計算する過程は以下で書き下すことができる。
【００３８】
【数１１】
2018121108.tif 000013
【００３９】
ビームフォーミング出力信号Yl,ω,τには、l番目の音源群を強調した信号と抑圧はされているものの周囲の雑音とが混在する。そのとき、φY_l,ω,τ(下付添え字Y_l,ω,τはYl,ω,τを意味する)を構成する要素のうちl番目の音源PSDが支配的であり、式(11)の第二項目の演算は、周囲の雑音を減算する役割を担う。図３に示すように、φY_l,ω,τから周囲雑音を減算するためには、PSD分離行列の対角項λl,ωが正の実数で(強調する役割)、非対角項γl,i,ωは0以下の実数(減算する役割)になるべきである。
【００４０】
【数１２】
2018121108.tif 000014
【００４１】
だから、高精度に音源PSDを推定するための仮説は、式(12)を満たすPSD分離行列GL,ωを設計することで表される。
【００４２】
PSD分離行列GL,ωの対角項は正、非対角項は負であり、その逆行列(つまり、PSD混合行列DL,ω)が定義通りに非負の行列となるとき、PSD分離行列GL,ωは、M行列と呼ばれる。故に、上述した高精度に音源PSDを推定するための仮説は、GL,ωをM行列になるように設計するという言葉で言い換え可能である。
【００４３】
(PSD分離行列がM行列となるための条件の導出)
PSD分離行列GL,ωがM行列となるための条件を導出する。PSD分離行列GL,ωの逆行列を解析的にモデル化するために、以下のように行列を分解する。
【００４４】
【数１３】
2018121108.tif 000015
【００４５】
ここで、ΛL,ω正の対角行列、BL,ωは非負の非対角行列、ILは単位行列である。
【００４６】
【数１４】
2018121108.tif 000016
【００４７】
PSD分離行列GL,ωの逆行列は以下で計算される。
【００４８】
【数１５】
2018121108.tif 000017
【００４９】
Λ-1L,ωは非負行列であるから、PSD混合行列DL,ωが定義通りに非負行列になるための条件は、(IL-BL,ω)-1が非負行列になることである。
【００５０】
(IL-BL,ω)-1が非負行列になるための条件を導出する。ノイマン級数展開することで以下のように展開される。
【００５１】
【数１６】
2018121108.tif 000018
【００５２】
なお、式(18)が成り立つ条件は、BL,ωのスペクトルノルム(行列の特異値の絶対値の最大値)が1未満になることである。このことは、以下の式で表される。
【００５３】
【数１７】
2018121108.tif 000019
【００５４】
ここで、σi,ωは、BL,ωのi番目の特異値を表す。まとめると、PSD分離行列GL,ωがM行列になるための条件は、式(19)で表される。つまり、BL,ωのスペクトルノルムが1未満になればいい。
【００５５】
(従来法の処理フロー)
非特許文献１では、BL,ωの解析的なモデル化が難しいことから、式(19)の関係が成り立つと仮定して、PSD混合行列DL,ωとPSD分離行列GL,ωの算出方法について提案した。簡略的ではあるが、ビームフォーミングフィルタ係数群WL,ω、PSD混合行列DL,ω、PSD分離行列GL,ωを推定するフローについて図４に載せる。
【００５６】
まず、最小分散無歪応答法(参考文献１)により、ビームフォーミングフィルタを設計する（Ｓ１００）。
(参考文献１)浅野太，「音のアレイ信号処理」,コロナ社,2011, pp.86-90.
なお、L個のビームフォーマの焦点はそれぞれ異なる方向{Θ1,…,ΘL}であるとし、それらの係数は行列WL,ωに格納されるとする。
【００５７】
【数１８】
2018121108.tif 000020
【００５８】
次に、行列WL,ωを用いてPSD混合行列を計算する(Ｓ１０１)。
【００５９】
【数１９】
2018121108.tif 000021
【００６０】
ここで、HL,ωはアレイマニフォールドベクトル(自由音場にアレイを設置したときに得られるであろう伝達特性の模擬値)である。
【００６１】
【数２０】
2018121108.tif 000022
【００６２】
なお、Hm,Θ_l,ω(下付き添え字m,Θ_l,ωはm,Θl,ωを意味する)は、m番目のマイクロホンに対して方向Θlから音が到来したときの自由音場における伝達特性の模擬値を表し、事前に計算可能である。計算されたΔL,ωのパワーを式(5)のように計算することで、PSD混合行列DL,ωを求め、出力する。
【００６３】
【数２１】
2018121108.tif 000023
【００６４】
最後に、PSD分離行列GL,ωを近似計算する(Ｓ１０２)。式(19)の関係が成り立つと仮定し、ΛL,ω=ILとするのであれば、PSD混合行列DL,ωを以下のように近似的にモデル化することができる。
【００６５】
【数２２】
2018121108.tif 000024
【００６６】
PSD混合行列DL,ωは既に計算されているのだから、ΓL,ωを近似的に求めることはできる。その近似したΓL,ωを使えば、PSD分離行列GL,ωもまた近似的に求めることができる。
【００６７】
【数２３】
2018121108.tif 000025
【００６８】
一連の操作により、目的の量(WL,ω、DL,ω、GL,ω)を求めることはできたが、近似操作を伴うために、PSDの推定誤差を増大させてしまうという欠点があった。これが、従来技術の課題の詳細である。
【００６９】
(本実施形態の原理)
まず、音源PSDを高精度に推定するためのPSD混合/分離行列の設計方式に関する原理を述べる。次に、設計したビームパターンに対応したフィルタ係数を解析的に導出するための方法について述べる。
【００７０】
(音源PSDを高精度に推定するためのPSD混合/分離行列の設計)
BL,ωの特異値を解析的にモデル化するために、PSD混合/分離行列の構造を巡回性を持つ行列（以下、巡回行列ともいう）になるように拘束をかける。この物理的意味合いについて、まず説明する。次に、PSD分離行列GL,ωが巡回M行列となるような条件を導出する。
【００７１】
(対称性ビームフォーマを設計した場合のPSD混合/分離行列のモデル化)
BL,ωの特異値を解析的にモデル化するための一つの方法として、PSD混合行列DL,ωが巡回行列となるように制約を課すこととした。
【００７２】
【数２４】
2018121108.tif 000026
【００７３】
式(25)を構成するL個のビームフォーマの特性は基本的に同じで、回転させているだけであることが分かる。式(25)のように、PSD混合行列DL,ωが巡回行列になるためには、M個のマイクロホンの空間的配置の中心を通るある回転軸に対して1回転させたときに、M個のマイクロホンの空間的配置と同一の空間的配置が複数回現れるようにM個のマイクロホンが配置される必要がある。言い換えると、回転対称性を満たすように、M個のマイクロホンが配置される必要がある。例えば、等間隔の円形アレイを用いて音場を観測し、対称性を持たせてビームフォーマを設計することで達成されるだろう。他にも、正多面体の頂点にマイクロホンを配置する方法や、それらを多重化して配置する(例えば、二つの同心円を用いて、大きい円形アレイと、小さい円形アレイとを用いて二重化して配置する)方法が考えられる。このとき、マイクロホンは、バッフルに取り付けてもよいし、バッフルを設けずに、例えば、糸で吊るしたり、細い棒で固定して、空中に配置してもよい。
【００７４】
式(19)の条件が成立することをここでは仮定すると、式(17)、(18)より、DL,ωは以下のようにモデル化される。
【００７５】
【数２５】
2018121108.tif 000027
【００７６】
PSD混合行列DL,ωが巡回行列であるとき、PSD混合行列DL,ωを構成するBL,ωとΛL,ωもまた巡回行列になることが知られている。
【００７７】
【数２６】
2018121108.tif 000028
【００７８】
ここで、αωは任意の正の実数である。
【００７９】
巡回行列の基本的な性質より、解析対象の行列BL,ωの固有ベクトルは、離散フーリエ変換の基底行列(DFT基底行列)で構成されることが知られている。
【００８０】
【数２７】
2018121108.tif 000029
【００８１】
ここで、Hは共役転置を表す。FLはDFT基底行列であり、ωとは独立である。
【００８２】
【数２８】
2018121108.tif 000030
【００８３】
固有値行列ΨL,ωはDFTの出力値で構成されることになる。
【００８４】
【数２９】
2018121108.tif 000031
【００８５】
式(19)の条件を満たすとき、BL,ωとΛL,ωは解析的にモデル化できることが分かった。その時、PSD混合行列DL,ωやPSD分離行列GL,ωもまた解析的にモデル化できる。式(29)より、PSD混合行列DL,ωの構成要素である(IL-BL,ω)-1は、以下でモデル化できる。
【００８６】
【数３０】
2018121108.tif 000032
【００８７】
ここで、
【００８８】
【数３１】
2018121108.tif 000033
【００８９】
である。式(13)と(26)より、PSD混合行列DL,ωとPSD分離行列GL,ωは以下でモデル化される。
【００９０】
【数３２】
2018121108.tif 000034
【００９１】
式(36)より、PSD分離行列GL,ωは巡回性を持つM行列(巡回M行列)になることが分かる。
【００９２】
(PSD分離行列GL,ωが巡回M行列となるための条件の導出)
PSD分離行列GL,ωが巡回M行列となるための条件を導出する。BL,ωは非負非対角行列であるから、DFT出力成分の一つである直流成分ψ1,ωが{ψ1,ω,…,ψL,ω}の中で最大値をであることは明らかである。
【００９３】
【数３３】
2018121108.tif 000035
【００９４】
f0L=1であるので、ψ1,ωは解析的に計算できる。
【００９５】
【数３４】
2018121108.tif 000036
【００９６】
よって、式(19)にあるBL,ωのスペクトルノルムは以下で書ける。
【００９７】
【数３５】
2018121108.tif 000037
【００９８】
故に、PSD分離行列GL,ωを巡回M行列にするという条件は、ψ1,ωを1未満にすることで達成される。
【００９９】
ψ1,ωを1未満にするための条件を導出する。BL,ωの対角項は全て0であることと、FLBL,ωFHL=ΨL,ωであることから、１つの条件は以下で書き表すことができる。
【０１００】
【数３６】
2018121108.tif 000038
【０１０１】
BL,ωのスペクトルノルムは1未満であるべきなので、式(40)の条件は以下のように書き直すことができる。
【０１０２】
【数３７】
2018121108.tif 000039
【０１０３】
その上で、{ψ1,ω,…,ψL,ω}は、実数系列{0,B2,ω,…,BL,ω}を入力した際にのDFT出力群であることから、式(42)の関係は成り立つ。
【０１０４】
【数３８】
2018121108.tif 000040
【０１０５】
まとめると、GL,ωが巡回M行列となるための条件は、式(41)と(42)で表される。式(42)の複素共役の関係より、ΨL,ωを決めるための自由度は
【０１０６】
【数３９】
2018121108.tif 000041
【０１０７】
に限定される。
【０１０８】
数値例を示すことで理解の促進を図りたい。L=6として、PSD混合行列DL,ωやPSD分離行列GL,ωを具体的に計算する。式(41)の条件を満たす中で、
【０１０９】
【数４０】
2018121108.tif 000042
【０１１０】
個の要素を任意に設計する。
【０１１１】
【数４１】
2018121108.tif 000043
【０１１２】
ここで、ψ2,ω,ψ3,ω,ψ4,ωの選択肢として複素数でも構わないのであるが、ここでは実数を使用した。式(42)より、残りの成分のいくつかは自動的に決まる。
【０１１３】
【数４２】
2018121108.tif 000044
【０１１４】
最後に、式(41)によりψ1,ωを計算する。
【０１１５】
【数４３】
2018121108.tif 000045
【０１１６】
ψ1,ωが1未満であるという条件を満たすので、式(35)や(36)を使って、DL,ωやGL,ωを解析的に設計することが可能である。ΛL,ω=ILとして、DL,ω、GL,ωはそれぞれ以下で計算される。
【０１１７】
【数４４】
2018121108.tif 000046
【０１１８】
ここで、diag(・)は対角行列化する演算子である。PSD分離行列GL,ωが巡回M行列であることが確認できる。
【０１１９】
(解析的なビームフォーミングのフィルタ係数の導出方法)
まず、式(35)に従うようなビームパターンを得るためのフィルタ係数の導出方法について、説明する。次に設計したビームフォーミングの物理的な性質について調査を行う。
【０１２０】
(ビームフォーミングのフィルタ係数の導出)
ビームフォーミングのフィルタ係数を解析的に導出するために、MがL以上であること(M≧L)を前提とする。ビームパターン(ビームフォーマの感度)ΔL,ω、フィルタ係数行列WL,ω、アレイマニフォールドベクトル(自遊空間における伝達特性の模擬)HL,ωの関係は以下で定義される。
【０１２１】
【数４５】
2018121108.tif 000047
【０１２２】
ここで、
【０１２３】
【数４６】
2018121108.tif 000048
【０１２４】
である。ΔL,ωは、PSD混合行列DL,ωの各要素の平方根を求めることで求められる。また、PSD混合行列DL,ωが巡回行列であれば、ΔL,ωもまた巡回行列であるので、以下のように表現できる。
【０１２５】
【数４７】
2018121108.tif 000049
【０１２６】
ここで、
【０１２７】
【数４８】
2018121108.tif 000050
【０１２８】
である。L=M、かつHL,ωがフルランクの行列であるとき、フィルタ係数行列WL,ωは以下で計算される。
【０１２９】
【数４９】
2018121108.tif 000051
【０１３０】
しかしながら、L<Mであるとき、少し工夫が必要である。それは、HL,ωの疑似逆行列は正確に求まらないためである。この問題を解決するための基本的なアイデアは、L個のビームパターンの特性を保持したままで、ΔL,ωの次元をLからMへ上げることである。このアップコンバージョン処理を通じて、WL,ω、ΔL,ω、HL,ωはM次元の正方行列WM,ω、ΔM,ω、HM,ωになる。アップコンバージョンされたΔM,ωは以下でモデル化される。
【０１３１】
【数５０】
2018121108.tif 000052
【０１３２】
ここで、
【０１３３】
【数５１】
2018121108.tif 000053
【０１３４】
である。ΥΔ,L,ωやΥΔ,M,ωは、ΔL,ωやΔM,ωのDFT出力値により構成されるのだから、L個のビームパターンの特性を維持したままに、それらを変換することは容易に達成できる。Lが偶数であるとき、ΔM,ωは以下でモデル化される。
【０１３５】
【数５２】
2018121108.tif 000054
【０１３６】
Lが奇数であるとき、以下の変換が成り立つ。
【０１３７】
【数５３】
2018121108.tif 000055
【０１３８】
上述の(対称性ビームフォーマを設計した場合のPSD混合/分離行列のモデル化)で述べたように、ここでは、ΔM,ωを巡回行列とするために、等間隔の円形アレイで観測することを前提とする。図５に示すように、等間隔に離散化されたM方向から平面波が到来することを仮定したとき、HM,ωは以下でモデル化される。
【０１３９】
【数５４】
2018121108.tif 000056
【０１４０】
ここで、τl,mはl番目の方向から平面波が到来したときのm番目のマイクロホンと円形アレイの中心の間の遅延時間差であり、以下で計算可能である。
【０１４１】
【数５５】
2018121108.tif 000057
【０１４２】
ここで、dとcはそれぞれ、アレイ半径と音速を表す。τl,mはインデックスの差(l-m)にのみ依存するので、HM,ωは巡回行列になる。故に以下のように分解することが可能である。
【０１４３】
【数５６】
2018121108.tif 000058
【０１４４】
ここで、
【０１４５】
【数５７】
2018121108.tif 000059
【０１４６】
である。式(57)、(58)、(64)より、L<Mであるときのビームフォーミングフィルタは以下により解析的に計算することができる。
【０１４７】
【数５８】
2018121108.tif 000060
【０１４８】
(設計したビームフォーマの物理特性の調査)
設計したビームフォーマの特性を調査するために、ΔM,ωを以下のように、多重極展開した。
【０１４９】
【数５９】
2018121108.tif 000061
【０１５０】
ここで、Mは偶数であることを想定した。Mが奇数であるときは、以下で展開される。
【０１５１】
【数６０】
2018121108.tif 000062
【０１５２】
ΩM,l,ωに含まれるビームパターンは、2(l-1)番目の多重極に対応する。l=1であるとき、ΩM,1,ωはモノポールの指向性で構成される。
【０１５３】
【数６１】
2018121108.tif 000063
【０１５４】
ここで、
【０１５５】
【数６２】
2018121108.tif 000064
【０１５６】
である。ΩM,2,ωはダイポールで構成される。
【０１５７】
【数６３】
2018121108.tif 000065
【０１５８】
ここで、
【０１５９】
【数６４】
2018121108.tif 000066
【０１６０】
である。ΩM,l,ωの一般形は以下で書き表すことができる。
【０１６１】
【数６５】
2018121108.tif 000067
【０１６２】
ここで、νl,ωは、Mが偶数の場合、奇数の場合のそれぞれにおいて、式(74)、または式(75)で計算することができる。
【０１６３】
【数６６】
2018121108.tif 000068
【０１６４】
図６に、ΔM,ωを多重極展開したときの具体的な例を示す。なお、ここで使用したPSD混合行列DL,ωの数値は、式(49)と同じである。この図より、多重極の重みづけ加算により、設計したフィルタの指向特性が書き表されることが分かる。また、離散ではなく、連続的な指向特性が予測可能であることも分かる。多重極に対応した指向特性を広い周波数バンドに渡って得るために、マイクロホン間の距離を小さくしたほうがいい。
【０１６５】
＜第一実施形態＞
提案法の原理をまとめると、式(7)のようにDL,ωGL,ω=ILという拘束下で、PSD分離行列GL,ωがM行列となるようなPSD混合/分離行列の設計法とそれに対応したビームフォーミングフィルタの設計法であった。なお、従来方式では、DL,ωGL,ω=ILという関係は近似されているので成立しないし、ビームフォーミングフィルタの導出も解析的ではなかった。
【０１６６】
図７は第一実施形態に係るフィルタ係数算出装置２００の機能ブロック図を、図８はその処理フローを示す。
【０１６７】
フィルタ係数算出装置２００は、CPUと、RAMと、以下の処理を実行するためのプログラムを記録したROMを備えたコンピュータで構成され、機能的には次に示すように構成されている。
【０１６８】
フィルタ係数算出装置２００は、パラメータ設計部２０１と、分離行列設計部２０２と、フィルタ係数算出部２０３とを含む。
【０１６９】
＜パラメータ設計部２０１＞
パラメータ設計部２０１は、PSD分離行列GのパラメータΨL,ωが以下の条件を満たすようにパラメータΨL,ωを設計し（Ｓ２０１）、出力する。
【０１７０】
【数６７】
2018121108.tif 000069
【０１７１】
例えば、人手により、
【０１７２】
【数６８】
2018121108.tif 000070
【０１７３】
を、パラメータ設計部２０１が設計する。
【０１７４】
＜分離行列設計部２０２＞
分離行列設計部２０２は、パラメータΨL,ωを受け取り、PSD分離行列GL,ωがPSD混合行列DL,ωの逆行列となり、かつ、PSD分離行列GL,ωがM行列となるように、PSD分離行列GL,ωを設計し（Ｓ２０２）、出力する。例えば、式(35)、(36)により、DL,ω、GL,ωを解析的に設計する。
【０１７５】
【数６９】
2018121108.tif 000071
【０１７６】
＜フィルタ係数算出部２０３＞
フィルタ係数算出部２０３は、PSD混合行列DL,ωを受け取り、ビームフォーミングのフィルタ係数を解析的に
M=Lのとき、
【０１７７】
【数７０】
2018121108.tif 000072
【０１７８】
として算出し、M>Lのとき、
【０１７９】
【数７１】
2018121108.tif 000073
【０１８０】
として算出し（Ｓ２０３）、出力する。なお、式(57)のΔL,ωはPSD混合行列DL,ωの各要素の平方根として求められる行列、HL,ωはアレイマニフォールドベクトル(自由音場にアレイを設置したときに得られるであろう伝達特性の模擬値)である。式(65)のFMはDFT基底行列、ΥΔ,M,ωはΔM,ωのDFT出力値により構成される行列、ΔM,ω,HM,ωはΔL,ω,HL,ωに対してそれぞれアップコンバージョンを行った行列である。
【０１８１】
＜効果＞
以上の構成により、回転対称性を満たすように、M個のマイクロホンが配置される（例えば、等間隔円形アレイを用いて音場を観測する）ことを前提として、高精度に音源PSDを推定するために適した性質を持つPSD混合/分離行列、及びビームフォーミングのフィルタ係数を解析的に導出する方法について提案した。これまで、ヒューリスティックに決めてきた変数(PSD混合/分離行列、ビームフォーミングのフィルタ係数)を解析的に最適化することで、音源PSDの推定誤差を軽減し、クリアな目的音強調が可能になる。
【０１８２】
＜変形例＞
フィルタ係数算出部２０３のフィルタ係数を算出する式は、本実施形態の式に限定されない。本実施形態で求めたPSD混合/分離行列を用いて、他の式等によりフィルタ係数を算出してもよい。また、パラメータ設計部２０１と分離行列設計部２０２とを別装置とし、分離行列設計装置としてもよい。さらに、パラメータ設計部２０１を備えず、分離行列設計部２０２のみからなる離行列設計装置としてもよい。このときパラメータΨは人手により入力してもよいし、別装置から与えられてもよい。
【０１８３】
＜第二実施形態＞
第一実施形態のフィルタ係数算出装置２００で求めたフィルタ係数WL,ω及びPSD分離行列GL,ωを用いた音声強調装置３００について説明する。
【０１８４】
図９は第二実施形態に係る音声強調装置３００の機能ブロック図、図１０はその処理フローの例を示す図である。
【０１８５】
音声強調装置３００は、ビームフォーミング部３０１、ウィーナーフィルタ推定部３０２、ウィーナーフィルタリング部３０３、記憶部３０４を含む。
【０１８６】
音声強調装置３００は、M個のマイクロホンで収音した観測信号Xm,ω,τを入力とし、方向Θlから到来する音を強調した出力信号Zl,ω,τを求め、出力する。なお、M個のマイクロホンは、回転対称性を満たすように配置されるものとする。
【０１８７】
記憶部３０４には、第一実施形態で求めたフィルタ係数WL,ω及びPSD分離行列GL,ωを格納しておく。
【０１８８】
ビームフォーミング部３０１は、観測信号Xm,ω,τとフィルタ係数Wl,ω,τとを受け取り、次式により、ビームフォーミング出力信号Yl,ω,τを求め（Ｓ３０１）、出力する。
【０１８９】
【数７２】
2018121108.tif 000074
【０１９０】
ウィーナーフィルタ推定部３０２は、ビームフォーミング出力信号Yl,ω,τとPSD分離行列GL,ωとを受け取り、次式により、ウィーナーフィルタUl,ω,τを算出し（Ｓ３０２）、出力する。
【０１９１】
【数７３】
2018121108.tif 000075
【０１９２】
ウィーナーフィルタリング部３０３は、ビームフォーミング出力信号Yl,ω,τとウィーナーフィルタUl,ω,τとを受け取り、ビームフォーミング出力信号Yl,ω,τをウィーナーフィルタUl,ω,τでフィルタリングし（Ｓ３０３）、出力信号Zl,ω,τを求め、出力する。
【０１９３】
【数７４】
2018121108.tif 000076
【０１９４】
＜第三実施形態＞
第二実施形態と異なる部分を中心に説明する。
【０１９５】
第一実施形態の技術を応用すると、ユーザの操作に応じて強調したい音源の収音領域の角度幅を連続的に変えることができる。例えば、雑音となる成分が近傍にある場合には、収音領域を狭めることで、出力SN比を向上する。
【０１９６】
図１１は第三実施形態に係る音声強調装置４００の機能ブロック図、図１２はその処理フローの例を示す図である。
【０１９７】
音声強調装置３００は、収音領域指定部４０１、パラメータ設計部４０２、分離行列設計部２０２、フィルタ係数算出部２０３、ビームフォーミング部３０１、ウィーナーフィルタ推定部３０２、ウィーナーフィルタリング部３０３を含む。
【０１９８】
第一実施形態及び第二実施形態と異なる収音領域指定部４０１、パラメータ設計部４０２について説明する。
【０１９９】
収音領域指定部４０１は、図示しない入力手段(キーボート、マウス、タッチパネル等)を介して、ユーザによって指定された収音領域Aを受け取り、広さを表すためのパラメータμを求め(Ｓ４０１)、出力する。例えば、0≦μ＜1とし、μが小さいほど収音領域を狭め、μが1に近づくほど収音領域を広げることに対応する。収音領域Aを0≦μ＜1とし、収音領域Aを広さを表すためのパラメータμとしてそのまま出力してもよい。この値は音源強調の処理中に変えることが可能で、連続的に変えることで収音領域をスムーズに広げたり狭めたりすることができる。
【０２００】
パラメータ設計部４０２は、パラメータμを入力とし、パラメータμに対応するパラメータΨL,ωを設計し（Ｓ４０２）、出力する。収音領域のユーザ要求にマッチして、ΨL,ωを決める方法には任意性がある。大きな方針としては、μが大きくなるほど、ψ1,ωが大きな値になるようにすればよい。上述の＜設計したビームフォーマの物理特性の調査＞で説明したように、ψ1,ωは無指向性で収音した信号に掛かる係数に対応する。詳細は述べないが、ψ1,ωが大きくなるほど、フィルタの特性における無指向性収音信号の重みが大きくなるため、収音領域幅が大きくなる。例えば、ψ1,ω=μとし、他の値は式(41)と(42)を満たすようにランダムに値を割り振ったり同じ値を代入してもよい。
【０２０１】
【数７５】
2018121108.tif 000077
【０２０２】
例えば、１つ以上の閾値とパラメータμとの大小関係に基づき、予め用意しておいたパラメータΨL,ωを選択する構成としてもよい。例えば、2個の閾値T1,T2(T1＜T2)を用意し、
【０２０３】
【数７６】
2018121108.tif 000078
【０２０４】
とする。ただし、a1＜a2＜a3とする。また、所定の計算式に基づき、μを用いて、パラメータΨL,ωを計算により求める構成としてもよい。前述の通り、μが大きくなるほど、ψ1,ωが大きな値になるように、これらの閾値とパラメータΨL,ω、または所定の計算式を用意すればよい。
【０２０５】
図１３−図１６に、ψ1,ωを変えながら作成したフィルタの指向特性を載せる。ψ1,ωを大きくすることで、ビーム幅が広がることが分かる。なお、あまり、感度差がつかないようなビームフォーマを作成しても、正確に音源PSDを推定できないことから、ψ1,ωの最大値は0.8程度までにしたほうがいいことも数値シミュレーションより分かった。
【０２０６】
＜効果＞
このような構成とすることで、第二実施形態と同様の効果を得ることができる。さらに、柔軟に収音領域を変更できる。
【０２０７】
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０２０８】
＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
【０２０９】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【０２１０】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。
【０２１１】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０２１２】
また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
2018121108.tif 000085
【図２】
2018121108.tif 000086
【図３】
2018121108.tif 000087
【図４】
2018121108.tif 000088
【図５】
2018121108.tif 000089
【図６】
2018121108.tif 000090
【図７】
2018121108.tif 000091
【図８】
2018121108.tif 000092
【図９】
2018121108.tif 000093
【図１０】
2018121108.tif 000094
【図１１】
2018121108.tif 000095
【図１２】
2018121108.tif 000096
【図１３】
2018121108.tif 000097
【図１４】
2018121108.tif 000098
【図１５】
2018121108.tif 000099
【図１６】
2018121108.tif 000100