(57)【要約】【課題】複数の確率的仮定を同時に考慮したコスト関数を用いて潜在変数を最適化する技術を提供する。【解決手段】潜在変数~w*を最適化する最適化部を含む潜在変数最適化装置であって、最適化部は、vj(1≦j≦J)を行列Djとベクトルbjを用いてvj=Dj~w*+bjと表される潜在変数~w*の補助変数とし、補助変数vjのコスト項は、log-concaveである補助変数vjの確率分布を用いて表されるものであり、最適化部は、補助変数vjのコスト項の和を含むコスト関数の最小化問題を解くことにより潜在変数~w*を最適化する。【選択図】図２
2020126138.tif 000002
【特許請求の範囲】
【請求項１】
潜在変数~w*を最適化する最適化部を含む潜在変数最適化装置であって、
vj(1≦j≦J)を行列Djとベクトルbjを用いてvj=Dj~w*+bjと表される潜在変数~w*の補助変数とし、
補助変数vjのコスト項は、log-concaveである補助変数vjの確率分布を用いて表されるものであり、
前記最適化部は、
補助変数vjのコスト項の和を含むコスト関数の最小化問題を解くことにより潜在変数~w*を最適化する
潜在変数最適化装置。
【請求項２】
請求項１に記載の潜在変数最適化装置であって、
前記コスト関数は
【数３９】
2020126138.tif 000041
（ただし、L0は潜在変数~w*のコスト項、Lj(1≦j≦J)は補助変数vjのコスト項）で定義される
潜在変数最適化装置。
【請求項３】
請求項２に記載の潜在変数最適化装置であって、
uj(1≦j≦J)を補助変数vjの双対変数、^v=(v1T, …, vJT)T, ^D=(D1T, …, DJT)T, ^b=(b1T, …, bJT)T, ^u=(u1T, …, uJT)T、γを所定の定数とし、
前記最適化部は、
次式により、潜在変数~w*を更新する潜在変数更新部と、
【数４０】
2020126138.tif 000042
次式により、補助変数vj(1≦j≦J)を更新する補助変数更新部と、
【数４１】
2020126138.tif 000043
次式により、双対変数uj(1≦j≦J)を更新する双対変数更新部と
【数４２】
2020126138.tif 000044
を含む潜在変数最適化装置。
【請求項４】
ビームフォーマのフィルタ係数~w*=(w1*T, …, wF*T)（ただし、wf*(1≦f≦F)は周波数ビンfのフィルタ係数）を最適化する最適化部を含むフィルタ係数最適化装置であって、
ef,t(1≦f≦F, 1≦t≦T)を時間フレームtにおける周波数ビンfの推定非目的音を表すフィルタ係数~w*の補助変数、yf,t(1≦f≦F, 1≦t≦T)を時間フレームtにおける周波数ビンfの推定目的音を表すフィルタ係数~w*の補助変数、ηf(1≦f≦F-2)をηf=wf*-2wf+1*+wf+2*で定義されるフィルタ係数~w*の補助変数とし、
補助変数ef,t, yf,t, ηfのコスト項は、それぞれlog-concaveである補助変数ef,t, yf,t, ηfの確率分布を用いて表されるものであり、
前記最適化部は、
補助変数ef,t, yf,t, ηfのコスト項の和を含むコスト関数の最小化問題を解くことによりフィルタ係数~w*を最適化する
フィルタ係数最適化装置。
【請求項５】
請求項４に記載のフィルタ係数最適化装置であって、
前記コスト関数は
【数４３】
2020126138.tif 000045
（ただし、Rf(1≦f≦F)は周波数ビンfの非目的音に関する空間相関行列、β(>0)は所定の定数、λは所定の定数）で定義される
フィルタ係数最適化装置。
【請求項６】
請求項５に記載のフィルタ係数最適化装置であって、
ue,f,t(1≦f≦F, 1≦t≦T)を補助変数ef,tの双対変数、uy,f,t(1≦f≦F, 1≦t≦T)を補助変数yf,tの双対変数、uη,f(1≦f≦F-2)を補助変数ηfの双対変数、^u=[ue,1,1, …, ue,F,T, uy,1,1, …, uy,F,T, uη,1, …, uη,F-2]、γを所定の定数とし、
前記最適化部は、
次式により、フィルタ係数~w*を更新するフィルタ係数更新部と、
【数４４】
2020126138.tif 000046
（ただし、^D, ^bは次式により与えられ、
【数４５】
2020126138.tif 000047
zf,t(1≦f≦F, 1≦t≦T)は時間フレームtにおける周波数ビンfの観測音、hf(1≦f≦F)は周波数ビンfにおけるビーム方向のアレイマニフォールドベクトルである。）
次式により、補助変数ef,t(1≦f≦F, 1≦t≦T)、補助変数yf,t(1≦f≦F, 1≦t≦T)、補助変数ηf(1≦f≦F-2)を更新する補助変数更新部と、
【数４６】
2020126138.tif 000048
次式により、双対変数ue,f,t(1≦f≦F, 1≦t≦T)、双対変数uy,f,t(1≦f≦F, 1≦t≦T)、双対変数uη,f(1≦f≦F-2)を更新する双対変数更新部と
【数４７】
2020126138.tif 000049
を含むフィルタ係数最適化装置。
【請求項７】
ビームフォーマのフィルタ係数~w*=(w1*T, …, wF*T)（ただし、wf*(1≦f≦F)は周波数ビンfのフィルタ係数）を最適化する最適化部を含むフィルタ係数最適化装置であって、
ηf(1≦f≦F-2)をηf=wf*-2wf+1*+wf+2*で定義されるフィルタ係数~w*の補助変数とし、
補助変数ηfのコスト項は、log-concaveである補助変数ηfの確率分布を用いて表されるものであり、
前記最適化部は、
補助変数ηfのコスト項の和を含むコスト関数の最小化問題を解くことによりフィルタ係数~w*を最適化する
フィルタ係数最適化装置。
【請求項８】
潜在変数最適化装置が、潜在変数~w*を最適化する最適化ステップを実行する潜在変数最適化方法であって、
vj(1≦j≦J)を行列Djとベクトルbjを用いてvj=Dj~w*+bjと表される潜在変数~w*の補助変数とし、
補助変数vjのコスト項は、log-concaveである補助変数vjの確率分布を用いて表されるものであり、
前記最適化ステップは、
補助変数vjのコスト項の和を含むコスト関数の最小化問題を解くことにより潜在変数~w*を最適化する
潜在変数最適化方法。
【請求項９】
請求項８に記載の潜在変数最適化方法であって、
前記コスト関数は
【数４８】
2020126138.tif 000050
（ただし、L0は潜在変数~w*のコスト項、Lj(1≦j≦J)は補助変数vjのコスト項）で定義される
潜在変数最適化方法。
【請求項１０】
フィルタ係数最適化装置が、ビームフォーマのフィルタ係数~w*=(w1*T, …, wF*T)（ただし、wf*(1≦f≦F)は周波数ビンfのフィルタ係数）を最適化する最適化ステップを実行するフィルタ係数最適化方法であって、
ef,t(1≦f≦F, 1≦t≦T)を時間フレームtにおける周波数ビンfの推定非目的音を表すフィルタ係数~w*の補助変数、yf,t(1≦f≦F, 1≦t≦T)を時間フレームtにおける周波数ビンfの推定目的音を表すフィルタ係数~w*の補助変数、ηf(1≦f≦F-2)をηf=wf*-2wf+1*+wf+2*で定義されるフィルタ係数~w*の補助変数とし、
補助変数ef,t, yf,t, ηfのコスト項は、それぞれlog-concaveである補助変数ef,t, yf,t, ηfの確率分布を用いて表されるものであり、
前記最適化ステップは、
補助変数ef,t, yf,t, ηfのコスト項の和を含むコスト関数の最小化問題を解くことによりフィルタ係数~w*を最適化する
フィルタ係数最適化方法。
【請求項１１】
請求項１０に記載のフィルタ係数最適化方法であって、
前記コスト関数は
【数４９】
2020126138.tif 000051
（ただし、Rf(1≦f≦F)は周波数ビンfの非目的音に関する空間相関行列、β(>0)は所定の定数、λは所定の定数）で定義される
フィルタ係数最適化方法。
【請求項１２】
フィルタ係数最適化装置が、ビームフォーマのフィルタ係数~w*=(w1*T, …, wF*T)（ただし、wf*(1≦f≦F)は周波数ビンfのフィルタ係数）を最適化する最適化ステップを実行するフィルタ係数最適化方法であって、
ηf(1≦f≦F-2)をηf=wf*-2wf+1*+wf+2*で定義されるフィルタ係数~w*の補助変数とし、
補助変数ηfのコスト項は、log-concaveである補助変数ηfの確率分布を用いて表されるものであり、
前記最適化ステップは、
補助変数ηfのコスト項の和を含むコスト関数の最小化問題を解くことによりフィルタ係数~w*を最適化する
フィルタ係数最適化方法。
【請求項１３】
請求項１ないし３のいずれか１項に記載の潜在変数最適化装置、請求項４ないし７のいずれか１項に記載のフィルタ係数最適化装置のいずれかとしてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、目的音強調におけるフィルタ係数など最適化の対象となるモデルの潜在変数を最適化する技術に関する。
【背景技術】
【０００２】
特定の方角から到来する音のみを強調し他方向の雑音を抑圧する信号処理手法として、マイクロホンアレイを用いたビームフォーミングがよく知られている。この手法は電話会議システム、自動車内のコミュニケーションシステム、スマートスピーカー等で実用化されている。ビームフォーミングに関する従来手法の多くは、何らかの制約のもとでコスト関数最小化問題を解くことで最適なフィルタを導出していた。例えば、非特許文献１に記載のMVDRビームフォーマは、出力信号のパワーをコスト関数としてこれを目的音源方角に対する無歪み特性制約のもと最小化することで得られる。また、最尤（ML）ビームフォーマは、出力信号に含まれる雑音のパワーをコスト関数としてこれを最小化することで導出される。その他、ビームフォーマの性能を向上させるために、付加的な制約やコスト項をコスト関数に追加する試みもこれまで行われてきた。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】J. Capon, “High-resolution frequency-wavenumber spectrum analysis”, Proceedings of the IEEE, vol.57, no.8, pp.1408-1418, Aug. 1969.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ビームフォーマを現実の状況に適用する際、複数の特性を同時に持たせることができれば、応用上有用と考えられる。例えば、音声に対して高い強調性能を保ちつつ、低遅延特性も両立したビームフォーマが必要とされるような場面があるだろう。ビームフォーマの特性への要請は、理論的にはビームフォーマのフィルタ係数から定義される補助変数に対する確率的仮定という形式でモデル化できる。例えば、もし強調したい音が人間の音声であるという事前知識があったとすれば、推定される信号は時間周波数領域でラプラス分布等のスパース性の高い分布に従うと仮定するのが妥当であろう。また、フィルタ係数は周波数方向について連続かつなめらかに変化するのが自然であるということは経験的に知られているが、従来手法ではフィルタ係数の周波数方向に関する特性を考慮していなかったため、空間相関行列がランク落ちするような周波数ビンで解が不安定になり、この特性が満たされないという状況がみられた。もし、なめらかさを考慮した設計を行うことができれば、遅延の少ないフィルタが得られるという効果も期待される。以上のような仮定をフィルタ係数の推定に同時に組み込むことができれば、目的音強調だけでなく多様な特性を持ったビームフォーマが構成できると期待される。
【０００５】
しかしながら、これまではコスト関数の最適化に関する数理的手法の検討が十分ではなく、特に複数の確率的仮定を同時に考慮したコスト関数の最適化に関する検討は行われてこなかった。
【０００６】
そこで本発明では、複数の確率的仮定を同時に考慮したコスト関数を用いて潜在変数を最適化する技術を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明の一態様は、潜在変数~w*を最適化する最適化部を含む潜在変数最適化装置であって、vj(1≦j≦J)を行列Djとベクトルbjを用いてvj=Dj~w*+bjと表される潜在変数~w*の補助変数とし、補助変数vjのコスト項は、log-concaveである補助変数vjの確率分布を用いて表されるものであり、前記最適化部は、補助変数vjのコスト項の和を含むコスト関数の最小化問題を解くことにより潜在変数~w*を最適化する。
【発明の効果】
【０００８】
本発明によれば、複数の確率的仮定を同時に考慮したコスト関数を用いて潜在変数を最適化することが可能となる。
【図面の簡単な説明】
【０００９】
【図１】潜在変数最適化アルゴリズムを示す図である。
【図２】潜在変数最適化装置１００（フィルタ係数最適化装置１００）の構成を示すブロック図である。
【図３】潜在変数最適化装置１００（フィルタ係数最適化装置１００）の動作を示すフローチャートである。
【図４】最適化部１２０の構成を示すブロック図である。
【図５】最適化部１２０の動作を示すフローチャートである。
【発明を実施するための形態】
【００１０】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００１１】
各実施形態の説明に先立って、この明細書における表記方法について説明する。
【００１２】
_（アンダースコア）は下付き添字を表す。例えば、xy_zはyzがxに対する上付き添字であり、xy_zはyzがxに対する下付き添字であることを表す。
【００１３】
また、ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、明細書の記載表記の制約上、^xや~xと記載しているものである。
【００１４】
＜技術的背景＞
本発明の実施形態では、フィルタ係数自身やフィルタ係数から定まる補助変数に対する確率的仮定に基づいて設計したコスト関数を用いて、フィルタ係数を最適化（学習）する。ここで、補助変数はフィルタ係数のアフィン変換として表現されるものに限定されるが、例えば推定された目的音、出力信号に含まれる雑音、隣接周波数ビン間のフィルタ係数の差分などはこの範疇に含まれる。フィルタ係数やその補助変数に対して仮定される確率分布がすべてlog-concave（対数凹）であれば、フィルタ係数や補助変数を互いに独立とみなした際の同時分布もまたlog-concaveとなるため、負の対数尤度は凸関数となり、コスト関数最適化問題は線形関係式で制約された補助変数に関する凸関数の最適化問題に帰着される。この最適化問題は、例えば、交互方向乗数法(Alternating Direction Method of Multipliers, ADMM)を用いて解くことができ、最適なフィルタ係数を効率的に計算できる。
【００１５】
以下、上記説明した本発明の実施形態の原理について詳しく説明していく。まず、最初に、確率に基づく最適化の観点からビームフォーミングの問題を定式化し、従来のビームフォーミングの最適化問題がこの定式化の枠で記述できることを説明する。
【００１６】
《ビームフォーミングの問題の定式化》
ここでは、記号・ノーテーションを定義し、問題を定式化する。まず、ビームフォーミングの問題を数理的に記述するための種々の定義を行う。
【００１７】
空間中に単一の目的音源と複数個の妨害音源があり、M個の無指向性マイクからなるマイクロホンアレイでこれらの混合音を収録し、観測されたMチャネル信号をビームフォーミングフィルタに通すことで、特定の方角から到来する目的音を強調するという状況を考える。この状況を記述するモデルを導入するため、まず変数を定義する。以下の議論は、基本的に短時間フーリエ変換(short-time Fourier transform, STFT)ドメインで行う。
【００１８】
af∈CM(f=1, …, F)を周波数ビンfにおける目的音源からマイクロホンアレイへの伝達関数、aikf∈CMを周波数ビンfにおけるk番目の妨害音源からマイクロホンアレイへの伝達関数、sf,t∈Cを周波数ビンf、時間フレームt(t=1, …, T)における目的音信号、nikf,t∈CMを周波数ビンf、時間フレームtにおけるk番目の妨害音信号とする。これらの記号を用いて、マイクロホンアレイが観測する信号zf,tは、瞬時混合仮定のもと
【００１９】
【数１】
2020126138.tif 000003
【００２０】
と表される。ここで、nbf,tは、特定の妨害音源由来と仮定していない雑音信号（例えば、マイクの性能自体に起因する雑音）を表す。
【００２１】
我々が求めたいのは、観測信号zf,tから目的音信号sf,tの精度よい推定値yf,tを与えるような線形フィルタである。以下では、この線形フィルタのフィルタ係数をwf∈CMとする。推定値yf,tの時間フレームを表す添字tを省略して、推定値yfと表すことにすると、zf, yf, wfの関係は
【００２２】
【数２】
2020126138.tif 000004
【００２３】
で与えられる。ここで、Hは複素共役転置を表す。
【００２４】
ここで、フィルタ係数と観測音に依存する変数を導入する。フィルタを用いることで、観測音から目的音を抽出できるならば、目的音を観測音から差し引くことで妨害音源などに起因する非目的音も推定できるはずである。そこで、観測音に含まれる非目的音の推定値ef∈CMを
【００２５】
【数３】
2020126138.tif 000005
【００２６】
で定義する。ここで、hf∈CMは目的音源方角のアレイマニフォールドベクトルである。本来、式(1)のようなモデルのもとではアレイマニフォールドベクトルhfではなく伝達関数afを適用するのが望ましいが、実用上伝達関数を常に正確に知るのは難しい。そこで、式(3)の定義ではアレイマニフォールドベクトルhfを用いる。なお、ビームフォーマが適切に目的音を抽出できていれば、非目的音の推定値efは主に妨害音と背景雑音によって構成されると期待される。
【００２７】
ここで更に、目的音と非目的音の推定値はともにフィルタ係数のアフィン変換として表現できることに着目し、フィルタ係数を用いた変換式として表現する。
【００２８】
なお、以後、表記を簡潔にするため、周波数ビンfに関する添え字を持つ任意の変数xfに関して、すべての周波数ビンに関する情報をまとめたものを~x=(x1T, …, xFT)Tで表すこととする。
【００２９】
また、時間フレームtについて行列FtとGtを
【００３０】
【数４】
2020126138.tif 000006
【００３１】
で定義する。
【００３２】
すると、時間フレームtにおけるビームフォーマの出力である目的音の推定値~yt（以下、推定目的音という）と雑音の推定値~et（以下、推定雑音、推定非目的音という）はそれぞれ以下のフィルタ係数~w*のアフィン変換の形式で表現される。
【００３３】
【数５】
2020126138.tif 000007
【００３４】
（ただし、*は複素共役）
【００３５】
《従来のビームフォーミングの最適化問題》
まず、従来のビームフォーミングの最適化問題を、確率モデルの観点から定義されたコスト関数の最小化問題として記述する。
【００３６】
~y, ~e, ~w*が確率変数であると解釈し、これらの確率分布Py(~y), Pe(~e), Pw(~w*)は既に分かっているものと仮定する。このうち、確率分布Py(~y), Pe(~e)については、音の統計的な性質が反映されているものと期待される。一方、確率分布Pw(~w*)は、目的音源方角に対する周波数応答の仮定を表現するためにしばしば用いられるものである。これらの仮定のもと、観測音の時系列{~zt}t=1Tに対する確率変数~w*の尤度関数は
【００３７】
【数６】
2020126138.tif 000008
【００３８】
と表される。ただし、~yt, ~etは式(6)、式(7)で表される~w*のアフィン変換によって定まる。この尤度を~w*に関して最大化することで、確率モデルの意味で最適なフィルタが導出できる。尤度最大化は負の対数尤度最小化と等価であるので、解くべきは
【００３９】
【数７】
2020126138.tif 000009
【００４０】
という形式の問題になる。
【００４１】
様々な従来のビームフォーミングの最適化問題が式(9)に基づく定式化として解釈できる。以下では、具体例として非特許文献１のMVDRビームフォーマの最適化問題について説明する。
【００４２】
［フィルタ設計フェーズ：フィルタ係数~w*のコスト関数］
周波数ビンfで観測音zfの空間相関行列の推定値Rf=Ez_f[zfzfH](f=1, …, F)が既知であり、また観測音に含まれる推定非目的音~etが正規分布N(0, Rf)に従う（つまり、ef,t〜N(0, Rf)）と仮定する。
【００４３】
このとき、尤度関数の非目的音に対する項ΠtPe(~et; ~w*)は
【００４４】
【数８】
2020126138.tif 000010
【００４５】
と表される。他の項（ΠtPy(~yt; ~w*)やPw(~w*)）にはここでは特に確率分布の仮定をおかない。
【００４６】
［フィルタ設計フェーズ：フィルタ係数~w*の制約条件］
フィルタ係数~w*=(w1*T, …, wF*T)の各wf*には目的音源方角に対する無歪み制約wfHhf=1を課す。
【００４７】
［フィルタ設計フェーズ：最適化問題］
これらの仮定のもと、単純な平方完成によって式(9)に基づく最適化問題は周波数ビンfについて
【００４８】
【数９】
2020126138.tif 000011
【００４９】
（ただし、γf=(hfHRf-1hf)-1）という形に帰着できる。
【００５０】
［フィルタ設計フェーズ：コスト関数最適化］
式(12)の問題の解は、よく知られたMVDRビームフォーマ(すなわち、γfRf-1hf)である。
【００５１】
次に、以上の手続きで得られたフィルタ係数を使用して実際にビームフォーマを動作させる際の計算（フィルタ使用フェーズ）について説明する。
【００５２】
［フィルタ使用フェーズ：フィルタ再設計］
ビームフォーミングの処理では、観測音をフレーム毎に区切り、各フレームに対して離散フーリエ変換を行う必要があるが、リアルタイムでビームフォーミングを行う状況ではフレーム長が長いと遅延が大きくなってしまう。そこで、低遅延のフィルタを再設計する。まず、フィルタ設計フェーズで設計したフィルタ係数~w*に対して逆フーリエ変換を行い、フィルタを時間領域の表現に戻すことで、各マイクロホンm(m=1, …, M)のインパルス応答wm[i]を得る。そして、入力として与えられた指定フレーム長Ntapをもとに、各インパルス応答wm[i]のうち最初のNtap/2成分からなるベクトルw'm1[i]と最後のNtap/2成分からなるベクトルw'm2[i]のみを取り出し（つまり、残りの要素はすべて無視し）、長さがNtapに短縮された新たなインパルス応答
【００５３】
【数１０】
2020126138.tif 000012
【００５４】
を導入する。このインパルス応答w"m[i]を再度離散フーリエ変換することで、要素数がNtapに削減された（再設計された）フィルタ係数~w'*が計算される。
【００５５】
［フィルタ使用フェーズ：離散フーリエ変換（DFT）］
次に、ビームフォーミング処理対象となる観測音を時間方向にNtapサンプルずつ切り出して、切り出した各区間（フレーム）に対して離散フーリエ変換を施し、STFT領域での観測音~zを出力する。
【００５６】
［フィルタ使用フェーズ：畳み込み］
ここでは、STFT領域での観測音~zとフィルタ係数~w'*を入力とし、式(2)のたたみ込みを行い、STFT領域での推定目的音~yを出力する。
【００５７】
［フィルタ使用フェーズ：逆離散フーリエ変換（逆DFT）］
最後に、STFT領域での推定目的音~yに逆離散フーリエ変換を施し、ビームフォーミング処理を施した時間領域波形、つまり時系列の推定目的音を得る。
【００５８】
式(9)は、変数~w*に関する最適化問題という定式化になっており、上述のMVDRビームフォーマのような比較的単純な例に対しては容易に解くことができるが、コスト関数が複雑な式になる場合、同様に最適化を行うのは一般に困難である。このことから、従来の手法は2つの問題点を抱えていることがわかる。一つ目として、これまでの手法では目的音~yや非目的音~eに対して仮定される確率分布が正規分布のように単純なものに限定されてしまいがちであったが、正規分布は必ずしも音源分布の記述として適切ではない点がある。二つ目として、フィルタ係数~w*に対する制約として用いることができるコスト関数も制限されていて、特に多様な確率的仮定を同時に考慮するのは困難である点がある。これまでにもフィルタ係数~w*に対して付加的な制約の導入を検討したものもあったが、一般に多様な確率的仮定を同時に考慮することで構成されるような複雑なコスト関数の最小化は極めて難しい問題であった。特に、低遅延・安定性・高い雑音抑圧性能を同時に達成するようなビームフォーマを構成するのは困難であった。
【００５９】
上記問題を解決するために、式(6)の推定目的音~ytや式(7)の推定雑音~etのような補助変数に対しても確率的仮定をおき、コスト関数を各補助関数に関するコスト項の和として表現することを考える。以下では、この考えに基づくコスト関数の設計方法について説明する。
【００６０】
《複数の特性を考慮したコスト関数に基づく最適化問題》
ビームフォーミングの最適化問題のための、新たなコスト関数は、さまざまな凸関数の項の和として表現する。また、各項の引数は、フィルタ係数~w*や新たに導入される（推定目的音~ytや推定雑音~etのような）フィルタ係数~w*のアフィン変換として定義できる補助変数とする。。これらの要請を満たすコスト関数は、いずれも最適化問題が解きやすい。言い換えれば、これらの要請を満たす範囲内で自由にコスト関数を設計できる。以下、詳しく説明する。
【００６１】
［フィルタ設計フェーズ：フィルタ係数~w*の補助変数］
Jを任意の自然数とし、補助変数vj(j=1, …, J)を導入する。補助変数vjとしてはフィルタ係数~w*と線形の関係、つまり、線形関係式vj=Dj~w*+bjを満たすものを用いる。各補助変数vjとそれがみたす関係式は式(6)、式(7)の一般化になっており、この意味で線形関係式を満たすとの制約は従来の手法を包含するものである。
【００６２】
簡単のため、^v=(v1T, …, vJT)T, ^D=(D1T, …, DJT)T, ^b=(b1T, …, bJT)Tのような表記を以下では採用する。
【００６３】
［フィルタ設計フェーズ：フィルタ係数~w*のコスト項・補助変数vjのコスト項］
コスト関数Lは、フィルタ係数~w*のコスト項L0と補助変数vjのコスト項Lj(j=1, …, J)を用いて
【００６４】
【数１１】
2020126138.tif 000013
【００６５】
（ただし、Lj(j=0, …, J)は凸関数）という形式で表す
【００６６】
凸関数の和は凸関数であるから、コスト関数Lも凸関数となる。補助変数vjのコスト項が凸関数であるという制約は一見突飛かもしれないが、これは実は式(6)の推定目的音~ytや式(7)の推定雑音~etなどの補助変数の確率分布にlog-concaveなものを採用していることに他ならない。ここで、確率分布がlog-concaveであるとはその確率密度関数の対数の-1倍(negative log)が凸関数であるという意味である。正規分布やラプラス分布など音源の確率モデルの記述に一般的に用いられる確率分布の多くもこの性質を満たしている。式(14)のコスト項Ljは補助変数vjの確率密度関数の対数の-1倍と解釈できるため、コスト項の凸性はlog-concaveな確率分布のみを考える限り自動的に満たされる性質である。
【００６７】
［フィルタ設計フェーズ：最適化問題］
以上の検討から、我々が解くべき問題は
【００６８】
【数１２】
2020126138.tif 000014
【００６９】
という典型的な線形制約付き凸最適化問題に帰着できる。式(15)の問題は、フィルタ係数に関する項と補助変数に関する項に分割し、それぞれについて交互に最適化を行うことで解くことができる。この問題を解く具体的なアルゴリズムとして様々なものが知られており、一例として交互方向乗数法（ADMM）を採用したアルゴリズムがある（当該アルゴリズムについて後述する）。
【００７０】
続いて、式(15)の問題定式化において、補助変数に課す確率的仮定として多様なものが使用可能であることを示すため、高い雑音抑圧性能を保ちつつ低遅延かつ音声の強調に適したフィルタを設計する例について説明する。
【００７１】
《複数の特性を考慮したビームフォーマの具体的設計例》
ここでは問題設定として実践的な状況を仮定し、式(15)の枠組みでコスト関数を具体的に設計する例を示す。具体的には、複数の妨害音が鳴っている環境下で、既知の位置から発せられた音声をストリーミング配信する、という状況を仮定する。なお、妨害音源は周波数ビン毎に複素正規分布に従う雑音を発しているものとする。本状況では、目的音源が音声であるという情報を反映した上で、遅延が少なくかつ高い強調性能を保つようなビームフォーマが要望されるだろう。
【００７２】
［フィルタ設計フェーズ：フィルタ係数~w*のコスト項］
上記状況では特にフィルタ係数~w*におくべき制約は存在しないので、フィルタ係数~w*に対するコスト項L0は考えないことにする。つまり、L0(~w*)=0とする。
【００７３】
［フィルタ設計フェーズ：フィルタ係数~w*の補助変数、補助変数のコスト項］
続いて、音源に関して既知の情報やビームフォーマに要望する特性の一つ一つについて検討し、補助変数とそのコスト項の設計を行っていく。
【００７４】
まず、目的音の分布を考える。上記仮定した状況では目的音は音声であるという情報が既知である。音声はスパース性という性質を持つことが知られているため、推定目的音がスパースな確率分布に従うという仮定を考慮したコスト項を設計することで、この既知の情報が活かせると考えられる。そこで、補助変数として推定目的音~ytを採用する。補助変数~ytの定義は式(6)と同一である。そして、補助変数~ytが、次式のラプラス分布に従うという仮定をおく。
【００７５】
【数１３】
2020126138.tif 000015
【００７６】
ここでβ(>0)は分布の形状を定める定数パラメータである。ラプラス分布はスパースな変数の分布の表現にしばしば用いられるものであり、上記仮定した状況では適切であると考えられる。式(16)の仮定のもと、補助変数~yのコスト項Lyはラプラス分布の対数の-1倍である
【００７７】
【数１４】
2020126138.tif 000016
【００７８】
という形になる。ラプラス分布はlog-concaveであるので、コスト項Lyは凸関数になり、式(15)の枠組みのもとで扱える。
【００７９】
次に、非目的音に関しても何らかの確率分布を仮定し、同様に補助変数・コスト項の導入を行う。観測音に含まれる非目的音の推定量として、式(7)で定義される推定非目的音~etを補助変数として導入する。上記仮定した状況では、主に妨害音からなる非目的音は正規分布に従うと仮定する。つまり、補助変数~etは
【００８０】
【数１５】
2020126138.tif 000017
【００８１】
という確率分布に従って出力されるとみなす。ここで、Rfは非目的音に関する空間相関行列であり、観測データから見積もることができる。式(18)の仮定を、補助変数~eに対するコスト項の形式に書き換えると、
【００８２】
【数１６】
2020126138.tif 000018
【００８３】
という形になる。正規分布はlog-concaveであるので、このコスト項Leも凸関数である。
【００８４】
ここで更なる補助変数とコスト項を導入することで、ビームフォーマに低遅延性を持たせることを目指す。そのためフィルタ係数~w*にどのようなコスト項を課せば低遅延なフィルタが実現可能か検討する。従来の広帯域ビームフォーマは周波数ビン毎に独立してフィルタ係数を導出しており、隣接周波数ビン間の関連性は考慮していなかった。しかし、周波数ビン方向に不連続やなめらかでないという周波数特性は、時間領域において裾の長いインパルス応答の原因となる。また、位相遅れを引き起こす群遅延も抑制されることが望ましい。このような特性を持たないフィルタ係数を解として得るために、フィルタ係数の周波数ビン方向に関する差分を新たな補助変数として導入し、これらの補助変数（のノルム）を小さくするようなコスト項を課すのが有効であると考えられる。具体的には、新たに
【００８５】
【数１７】
2020126138.tif 000019
【００８６】
というF-2個の補助変数ηfを定義する。式(20)は、ηfがフィルタの振幅・位相特性の周波数方向に関する2階微分の情報を含むように意図されている。式(20)を用いて、補助変数ηfに対するコスト項Lη_f(ηf)を
【００８７】
【数１８】
2020126138.tif 000020
【００８８】
で定義する。
【００８９】
［フィルタ設計フェーズ：コスト関数最適化］
以上のように、式(18)、式(16)、式(20)のような補助変数に関する仮定をおいたことで、コスト関数Lは、各コスト項の和
【００９０】
【数１９】
2020126138.tif 000021
【００９１】
となる。このコスト関数Lに出現する2FT+F-2個の補助変数はみなフィルタ係数~w*のアフィン変換として表されるので、式(22)の最小化問題は式(15)の具体例である。
【００９２】
ここまではビームフォーミングを対象として最適化問題の議論をしてきたが、これまで説明した数理的な枠組みはより汎用的な適用範囲を持つものであり、音響処理に限られるものではない。本枠組みの汎用性を端的に示すため、上記枠組みを画像処理に適用した例について説明する。
【００９３】
《画像処理における最適化問題の一例》
例えば、同じ形状の物体が大量に写っている画像のように周期的な絵柄の画像（以下、元の画像という）にノイズが重畳された画像が入力として与えられ、当該画像からノイズを除去した画像を得たいという状況を考える。元の画像の各画素の値を表す行列をS=[Sx,y]1≦x≦X,1≦y≦Y、各画素に加わるノイズを表す行列をNと表す。ノイズの値は画素毎に独立に平均0・分散1の正規分布に従って生成されるものとする。我々が観測できるのはノイズを含んだ画像Y=S+Nである。このとき、画像Yから元の画像Sを精度よく推定する問題を考えるため、行列Sを式(15)における~w*とみなし、行列Sや行列Sのアフィン変換によって定まる補助変数に関するコスト項を構成していく。
【００９４】
［フィルタ設計フェーズ：行列Sのコスト項］
まず、推定結果として得られた画像は元の画像と概ね一致していてほしいので、行列Sに関するコスト項として、入力画像Yの各画素との二乗誤差を課すことにする。行列Sに関するコスト項を具体的に書くと、次式になる。
【００９５】
【数２０】
2020126138.tif 000022
【００９６】
式(24)のコスト項Lsは凸関数である。
【００９７】
［フィルタ設計フェーズ：補助変数・補助変数のコスト項］
次に、適切にノイズを除去するための補助変数とそのコスト項を設計する。画像は通常滑らかで、隣接画素間の値の変動は小さいということを我々は経験的に知っている。画素毎に独立なノイズはこの性質に反した不自然な振る舞いを示すため、この不自然さを嫌うようなコスト項を設計することでノイズが除去できると考えられる。そこで、隣接画素間の差分として定義される量D1, D2を次式で与えられる補助変数として導入する。
【００９８】
【数２１】
2020126138.tif 000023
【００９９】
そして、滑らかでより自然性の高い画像ならば補助変数D1, D2の絶対値は傾向として小さくなるはずである。そこで、補助変数D1, D2に対して次のような凸なコスト項を課す。
【０１００】
【数２２】
2020126138.tif 000024
【０１０１】
これらのコスト項LD1, LD2は、ノイズ除去の意味合いを持つコスト項である。
【０１０２】
ここでは更に、元の画像が周期的な構造を有するという事前知識を持っている状況を仮定し、このような事前知識を活かせる補助変数やコスト項を設計する。周期的な画像では、画像を2次元フーリエ変換して得られた空間周波数スペクトルがスパースな構造を持つと期待される。2次元フーリエ変換はアフィン変換として記述できるので、空間周波数スペクトルを補助変数として採用することにし、これらの補助変数をスパースに誘導するようなコスト項を設計すれば我々の目的が達成されると考えられる。具体的には、画像の2次元フーリエ変換R=[Rk,j]を補助変数として導入する。これらは離散フーリエ変換行列Wk,jによって
【０１０３】
【数２３】
2020126138.tif 000025
【０１０４】
という式で定義でき、確かに行列Sのアフィン変換になっている。コスト項としては、
【０１０５】
【数２４】
2020126138.tif 000026
【０１０６】
という形の凸関数を仮定する。
【０１０７】
［フィルタ設計フェーズ：コスト関数最適化］
以上のコスト項の設計により、最適化すべきコスト関数Lは
【０１０８】
【数２５】
2020126138.tif 000027
【０１０９】
という形になる。式(31)の変数のうち行列Sが推定対象の変数であり、その他の変数は行列Sの補助変数である。
【０１１０】
以上の議論から、画像処理の場合においても式(15)の枠組みでコスト関数が設計できることがわかる。
【０１１１】
《ADMMに基づく最適化アルゴリズム》
図１は、式(15)で表される線形制約付き凸最適化問題を実際に解くための反復アルゴリズムを示す図である。当該アルゴリズムは式(15)の問題を効率的に解くアルゴリズムの一つとして知られるADMMに基づくものである。ADMMは元の問題の双対問題に対して最適化を行うようなアルゴリズムであり、補助変数vjと同じ次元を持つ双対変数ujを用いる。
【０１１２】
以下では、ビームフォーマの事例の一つの例として構成したコスト関数(22)に対して図１のアルゴリズムを適用した例について説明する。ここでは、図１の式より具体的な反復更新式を導出する。
【０１１３】
まず、変数~w*の更新則に関してはコスト項L0が式(22)では消えていることから、図１のステップ３の式は
【０１１４】
【数２６】
2020126138.tif 000028
【０１１５】
という形式に帰着される。この式に現れる行列^DH^D=ΣjDjHDjは
【０１１６】
【数２７】
2020126138.tif 000029
【０１１７】
という形式のブロック帯行列となるため、行列^DH^Dのコレスキー分解を行うことにより更新の際必要になる(^DH^D)-1による乗算の計算を効率化できる。
【０１１８】
続いて、補助変数の更新則を求める。この更新則は各コスト項の近接作用素として記述される。ここで、関数fの近接作用素proxfはproxf(x)=argminyf(y)+||x-y||22/2という形で定義される。この形とコスト項を見比べると、補助変数yf,tと補助変数ηfの更新則はl2ノルムの近接作用素
【０１１９】
【数２８】
2020126138.tif 000030
【０１２０】
で表されることがわかる。一方、補助変数ef,tに関するコスト項は単純な二次形式であるので、ef,tの更新式は定義から解析的に容易に導ける。結局、補助変数の更新則は
【０１２１】
【数２９】
2020126138.tif 000031
【０１２２】
という形になる。ここで、Iは単位行列を表す。
【０１２３】
《効果》
本発明の実施形態の原理は、ビームフォーマのフィルタ係数導出をコスト関数最適化問題として解釈した上で、フィルタ係数やその補助変数に対して個別のコスト項に基づく制約を課すことで、所望する複数の特性を兼ね備えたビームフォーマを設計するものである。
【０１２４】
従来法では、事前知識や所望する特性などの様々な要因を考慮した複雑なコスト関数を使用した設計を行うことができなかった。一方、本発明の実施形態の原理に従えば、補助変数という形で新たな変数を複数導入し、これらに関して個別にコスト項を設計するという枠組みでコスト関数を構成する。各コスト項は確率的仮定であるという意味を持っており、特にlog-concaveな確率的仮定を課す場合、線形制約付き凸最適化問題に帰着され、様々な数理的手法で比較的容易に最適化問題を解くことができる。これにより、複数の仮定を同時に考慮したフィルタ設計が可能となる。
【０１２５】
＜第１実施形態＞
以下、図２〜図３を参照して潜在変数最適化装置１００を説明する。図２は、潜在変数最適化装置１００の構成を示すブロック図である。図３は、潜在変数最適化装置１００の動作を示すフローチャートである。図２に示すように潜在変数最適化装置１００は、セットアップデータ計算部１１０と、最適化部１２０と、記録部１９０を含む。記録部１９０は、潜在変数最適化装置１００の処理に必要な情報を適宜記録する構成部である。記録部１９０は、例えば、最適化対象となる潜在変数を記録する。
【０１２６】
潜在変数最適化装置１００は、最適化用データを用いて、最適化の対象となるモデルの潜在変数~w*を最適化する。ここで、モデルとは、入力データを入力とし、出力データを出力とする関数（例えば、観測音を入力データとし、目的音を出力データとするビームフォーマのフィルタ）のことであり、最適化用データとは、潜在変数の最適化に用いる入力データ、または、潜在変数の最適化に用いる入力データと出力データの組のことをいう。
【０１２７】
図３に従い潜在変数最適化装置１００の動作について説明する。
【０１２８】
Ｓ１１０において、セットアップデータ計算部１１０は、最適化用データを用いて、潜在変数~w*を最適化する際に用いるセットアップデータを計算する。例えば、潜在変数~w*を最適化するために用いるコスト関数L
【０１２９】
【数３０】
2020126138.tif 000032
【０１３０】
（ただし、vj(1≦j≦J)は行列Djとベクトルbjを用いてvj=Dj~w*+bjと表される潜在変数~w*の補助変数、L0は潜在変数~w*のコスト項、Lj(1≦j≦J)は補助変数vjのコスト項）で用いるDj(1≦j≦J)、bj(1≦j≦J)、コスト項Li(0≦i≦J)に含まれるパラメータがセットアップデータの一例である。なお、コスト項Li(0≦i≦J)は凸関数とするのが好ましい。
【０１３１】
例えば、補助変数vj(1≦j≦J)のコスト項は、log-concaveである補助変数vjの確率分布を用いて表されるものとすると、コスト項Li(1≦i≦J)は凸関数となる。また、例えば、潜在変数~w*のコスト項L0=0としてもよく、コスト関数Lは補助変数vj(1≦j≦J)のコスト項の和を含むものであればよい。
【０１３２】
Ｓ１２０において、最適化部１２０は、コスト関数Lの最小化問題を解くことにより潜在変数~w*を最適化する。以下、図４〜図５を参照して最適化部１２０について説明する。図４は、最適化部１２０の構成を示すブロック図である。図５は、最適化部１２０の動作を示すフローチャートである。図４に示すように最適化部１２０は、初期化部１２１、潜在変数更新部１２２と、補助変数更新部１２３と、双対変数更新部１２４と、カウンタ更新部１２５と、終了条件判定部１２６を含む。
【０１３３】
図５に従い最適化部１２０の動作について説明する。
【０１３４】
Ｓ１２１において、初期化部１２１は、カウンタnを初期化する。具体的には、n=1とする。また、初期化部１２１は、補助変数^v=(v1T, …, vJT)T、双対変数^u=(u1T, …, uJT)Tを初期化する。さらに、初期化部１２１は、γにも初期値となる定数を設定する。
【０１３５】
Ｓ１２２において、潜在変数更新部１２２は、現時点で得られている補助変数^v、双対変数^uの値を用いて、次式により、潜在変数~w*を更新する。
【０１３６】
【数３１】
2020126138.tif 000033
【０１３７】
ここで、^D=(D1T, …, DJT)T、^b=(b1T, …, bJT)Tである。
【０１３８】
Ｓ１２３において、補助変数更新部１２３は、現時点で得られている潜在変数~w*、双対変数ujの値を用いて、次式により、補助変数vj(1≦j≦J)を更新する。
【０１３９】
【数３２】
2020126138.tif 000034
【０１４０】
Ｓ１２４において、双対変数更新部１２４は、現時点で得られている潜在変数~w*、補助変数vj、双対変数ujの値を用いて、次式により、双対変数uj(1≦j≦J)を更新する。
【０１４１】
【数３３】
2020126138.tif 000035
【０１４２】
Ｓ１２５において、カウンタ更新部１２５は、カウンタnを1だけインクリメントする。具体的には、n←n+1とする。
【０１４３】
Ｓ１２６において、終了条件判定部１２６は、カウンタnが所定の更新回数Niteration（Niterationは1以上の整数であり、例えば10万）に達した場合（つまり、n>Niterationとなり、終了条件が満たされた場合）は、そのときの潜在変数の値~w*を出力して、処理を終了する。それ以外の場合、Ｓ１２２の処理に戻る。つまり、最適化部１２０は、Ｓ１２２〜Ｓ１２６の処理を繰り返す。
【０１４４】
なお、式(*)で定義されるコスト関数Lを用いる場合、Jが2以上であっても最適化することができる。
【０１４５】
また、式(*)のように、コスト関数Lは、補助変数vjのコスト項をlog-concaveである補助変数vjの確率分布を用いて表されるものとして、補助変数vjのコスト項の和を含むものであればよく、例えば、コスト関数Lは、補助変数vjのコスト項をlog-concaveである補助変数vjの確率分布を用いて表されるものとして、補助変数vjのコスト項とlog-concaveな確率分布に基づいて定まるコスト項との和として表されるものであってもよい。
【０１４６】
本実施形態の発明によれば、潜在変数及び潜在変数から定まる補助変数に対する確率的仮定に基づくコスト関数を用いて潜在変数を最適化することが可能となる。
【０１４７】
[適用例]
ここでは、潜在変数最適化装置１００を音源強調に用いるビームフォーマのフィルタ係数の最適化に適用した例について説明する。そこで、以下では潜在変数最適化装置１００のことをフィルタ係数最適化装置１００と呼ぶことにする。フィルタ係数最適化装置１００の最適化対象はビームフォーマのフィルタ係数となる。フィルタ係数最適化装置１００の構成は図２に示す通りである。
【０１４８】
以下、図３に従いフィルタ係数最適化装置１００の動作について説明する。
【０１４９】
Ｓ１１０において、セットアップデータ計算部１１０は、最適化用データを用いて、フィルタ係数~w*=(w1*T, …, wF*T)（ただし、wf*(1≦f≦F)は周波数ビンfのフィルタ係数）を最適化する際に用いるセットアップデータを計算する。例えば、フィルタ係数~w*を最適化するために用いるコスト関数L
【０１５０】
【数３４】
2020126138.tif 000036
【０１５１】
（ただし、ef,t(1≦f≦F, 1≦t≦T)は時間フレームtにおける周波数ビンfの推定非目的音を表すフィルタ係数~w*の補助変数、yf,t(1≦f≦F, 1≦t≦T)は時間フレームtにおける周波数ビンfの推定目的音を表すフィルタ係数~w*の補助変数、ηf(1≦f≦F-2)はηf=wf*-2wf+1*+wf+2*で定義されるフィルタ係数~w*の補助変数、Rf(1≦f≦F)は周波数ビンfの非目的音に関する空間相関行列、β(>0)は所定の定数、λは所定の定数）で用いるコスト項Le,f,t(ef,t)=ef,tHRf-1ef,t(1≦f≦F, 1≦t≦T)、コスト項Ly,f,t(yf,t)=β|yf,t| (1≦f≦F, 1≦t≦T)、コスト項Lη,f(ηf)=λ||ηf||2(1≦f≦F-2)に含まれるパラメータがセットアップデータの一例である。なお、コスト項Le,f,t(ef,t) (1≦f≦F, 1≦t≦T), Ly,f,t(yf,t) (1≦f≦F, 1≦t≦T)、Lη,f(ηf) (1≦f≦F-2)はいずれも凸関数である。
【０１５２】
ただし、コスト項Le,f,t(ef,t), Ly,f,t(yf,t), Lη,f(ηf)は上記コスト項に限るものではなく、例えば、補助変数ef,t, yf,t, ηfのコスト項は、それぞれlog-concaveである補助変数ef,t, yf,t, ηfの確率分布を用いて表されるものであればよい。
【０１５３】
なお、上記コスト関数Lの定義式では、フィルタ係数~w*のコスト項L0=0となっている。
【０１５４】
Ｓ１２０において、最適化部１２０は、コスト関数Lの最小化問題を解くことによりフィルタ係数~w*を最適化する。以下、図４〜図５を参照して最適化部１２０について説明する。図４は、最適化部１２０の構成を示すブロック図である。図５は、最適化部１２０の動作を示すフローチャートである。ここでは、最適化部１２０に含まれる潜在変数更新部１２２のことをフィルタ係数更新部１２２と呼ぶことにする。
【０１５５】
以下、図５に従い最適化部１２０の動作について説明する。
【０１５６】
Ｓ１２１において、初期化部１２１は、カウンタnを初期化する。具体的には、n=1とする。また、初期化部１２１は、補助変数^v=[e1,1, …, eF,T, y1,1, …, yF,T, η1, …, ηF-2]、双対変数^u=[ue,1,1, …, ue,F,T, uy,1,1, …, uy,F,T, uη,1, …, uη,F-2]（ただし、ue,f,t(1≦f≦F, 1≦t≦T)は補助変数ef,tの双対変数、uy,f,t(1≦f≦F, 1≦t≦T)は補助変数yf,tの双対変数、uη,f(1≦f≦F-2)は補助変数ηfの双対変数とする）を初期化する。さらに、初期化部１２１は、γにも初期値となる定数を設定する。
【０１５７】
Ｓ１２２において、フィルタ係数更新部１２２は、現時点で得られている補助変数^v、双対変数^uの値を用いて、次式により、フィルタ係数~w*を更新する。
【０１５８】
【数３５】
2020126138.tif 000037
【０１５９】
ここで、^D, ^bは次式により与えられる。
【０１６０】
【数３６】
2020126138.tif 000038
【０１６１】
Ｓ１２３において、補助変数更新部１２３は、現時点で得られている潜在変数~w*、双対変数ue,f,t, uy,f,t, uη,fの値を用いて、次式により、補助変数ef,t(1≦f≦F, 1≦t≦T)、補助変数yf,t(1≦f≦F, 1≦t≦T)、補助変数ηf(1≦f≦F-2)を更新する。
【０１６２】
【数３７】
2020126138.tif 000039
【０１６３】
（ただし、zf,t(1≦f≦F, 1≦t≦T)は時間フレームtにおける周波数ビンfの観測音、hf(1≦f≦F)は周波数ビンfにおけるビーム方向のアレイマニフォールドベクトル）
【０１６４】
Ｓ１２４において、双対変数更新部１２４は、現時点で得られている潜在変数~w*、補助変数ef,t, yf,t, ηfの値を用いて、次式により、双対変数ue,f,t, uy,f,t, uη,fを更新する。
【０１６５】
【数３８】
2020126138.tif 000040
【０１６６】
Ｓ１２５において、カウンタ更新部１２５は、カウンタnを1だけインクリメントする。具体的には、n←n+1とする。
【０１６７】
Ｓ１２６において、終了条件判定部１２６は、カウンタnが所定の更新回数Niteration（Niterationは1以上の整数であり、例えば10万）に達した場合（つまり、n>Niterationとなり、終了条件が満たされた場合）は、そのときのフィルタ係数の値~w*を出力して、処理を終了する。それ以外の場合、Ｓ１２２の処理に戻る。つまり、最適化部１２０は、Ｓ１２２〜Ｓ１２６の処理を繰り返す。
【０１６８】
なお、式(**)のように、コスト関数Lは、補助変数ef,t, yf,t, ηfのコスト項をそれぞれlog-concaveである補助変数ef,t, yf,t, ηfの確率分布を用いて表されるものとして、補助変数ef,t, yf,t, ηfのコスト項の和を含むものであればよく、例えば、コスト関数Lは、補助変数ef,t, yf,t, ηfのコスト項をそれぞれlog-concaveである補助変数ef,t, yf,t, ηfの確率分布を用いて表されるものとして、補助変数ef,t, yf,t, ηfのコスト項とlog-concaveな確率分布に基づいて定まるコスト項との和として表されるものであってもよい。
【０１６９】
＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ−ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【０１７０】
ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。
【０１７１】
ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。
【０１７２】
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【０１７３】
既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【０１７４】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１７５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１７６】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記憶装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１７７】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【０１７８】
上述の本発明の実施形態の記載は、例証と記載の目的で提示されたものである。網羅的であるという意思はなく、開示された厳密な形式に発明を限定する意思もない。変形やバリエーションは上述の教示から可能である。実施形態は、本発明の原理の最も良い例証を提供するために、そして、この分野の当業者が、熟考された実際の使用に適するように本発明を色々な実施形態で、また、色々な変形を付加して利用できるようにするために、選ばれて表現されたものである。すべてのそのような変形やバリエーションは、公正に合法的に公平に与えられる幅にしたがって解釈された添付の請求項によって定められた本発明のスコープ内である。
【図１】
2020126138.tif 000052
【図２】
2020126138.tif 000053
【図３】
2020126138.tif 000054
【図４】
2020126138.tif 000055
【図５】
2020126138.tif 000056