(57)【特許請求の範囲】
【請求項１】
システムからの少なくとも１つの所望の信号を推定する信号処理装置であって、
該信号処理装置は、
第１のノイズ信号、第２のノイズ信号および第３のノイズ信号を含む複数のノイズ信号を生成するランダムノイズ生成器と、
該第１のノイズ信号を受信して、
【数１５０】
0004402879.tif 000164
に従って少なくとも１つの第１の時間変動パラメータを生成する第１のパラメータ生成器であって、
【数１５１】
0004402879.tif 000165
は、時間ｔ＋１におけるソース番号ｉの１つの第１の時間変動パラメータを表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、
【数１５２】
0004402879.tif 000166
は、時間ｔにおけるソース番号ｉの１つの第１の時間変動パラメータを表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、
【数１５３】
0004402879.tif 000167
は、時間ｔ＋１におけるソース番号ｉの該第１のノイズ信号を表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表す、第１のパラメータ生成器と、
該第２のノイズ信号を受信して、
【数１５４】
0004402879.tif 000168
に従って、少なくとも１つの通信チャンネルの通信経路を表す少なくとも１つの第２のパラメータを生成する第２のパラメータ生成器であって、
【数１５５】
0004402879.tif 000169
は、時間ｔ＋１におけるマイクロフォンｊでのソース番号ｉの１つの第２のパラメータを表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、
【数１５６】
0004402879.tif 000170
は、時間ｔにおけるマイクロフォンｊでのソース番号ｉの１つの第２のパラメータを表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、
【数１５７】
0004402879.tif 000171
は、時間ｔ＋１におけるマイクロフォンｊでのソース番号ｉの該第２のノイズ信号を表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、該システムのマイクロフォンｊは、ソース番号ｉとマイクロフォンｊとの間の該少なくとも１つの通信チャンネルを介して少なくとも１つの汚染信号を受信および測定し、該測定された少なくとも１つの汚染信号を該信号処理装置に提供し、該少なくとも１つの汚染信号は、該少なくとも１つの所望の信号と該複数のノイズ信号とを含む、第２のパラメータ生成器と、
該第３のノイズ信号と該少なくとも１つの第１の時間変動パラメータとを受信して、
【数１５８】
0004402879.tif 000172
に従って少なくとも１つの信号状態空間値を生成するサンプル生成器であって、
【数１５９】
0004402879.tif 000173
は、時間ｔ＋１におけるソース番号ｉの１つの信号状態空間値を表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表し、
【数１６０】
0004402879.tif 000174
は、時間ｔにおけるソース番号ｉの信号状態空間値のベクトル形式を表し、ｐｉは、時間ｔにおけるソース番号ｉの線形の組み合わせにおける以前の信号状態空間値の個数であり、ｋは、Ｎ個のランダムに選択された粒子のうちのｋ番目の粒子を表し、
【数１６１】
0004402879.tif 000175
は、時間ｔ＋１におけるソース番号ｉの該第３のノイズ信号を表し、ｋは、Ｎ個のランダムな粒子のうちのｋ番目の粒子を表す、サンプル生成器と、
該少なくとも１つの信号状態空間値と、該少なくとも１つの第１の時間変動パラメータと、該少なくとも１つの第２のパラメータとを受信して、
【数１６２】
0004402879.tif 000176
に従って該少なくとも１つの所望の信号を推定する信号推定手段であって、
Ｉは、ｆｔの関数としての該少なくとも１つの所望の信号のモンテカルロ推定を表し、Ｌは、平滑化遅延を表し、
【数１６３】
0004402879.tif 000177
は、該信号状態空間値を表す該少なくとも１つの所望の信号を表し、
【数１６４】
0004402879.tif 000178
は、該少なくとも１つの第１の時間変動パラメータおよび該少なくとも１つの第２のパラメータを表し、Ｎは、粒子の数を表す、信号推定手段と
を備えた、信号処理装置。
【請求項２】
前記少なくとも１つの所望の信号は、物理的プロセスにより生成される、請求項１に記載の信号処理装置。
【請求項３】
前記少なくとも１つの第１の時間変動パラメータは、前記物理的プロセスをモデリングする、請求項２に記載の信号処理装置。
【請求項４】
前記少なくとも１つの第２のパラメータは、固定され、かつ未知である、請求項１〜３のいずれか１つに記載の信号処理装置。
【請求項５】
前記少なくとも１つの第２のパラメータは、時間変動である、請求項１〜３のいずれか１つに記載の信号処理装置。
【請求項６】
前記少なくとも１つの第１の時間変動パラメータは、前記少なくとも１つの第２の時間変動パラメータの変化のレートとは異なる変化のレートを有する、請求項５に記載の信号処理装置。
【請求項７】
前記少なくとも１つの第１の時間変動パラメータは、平均して、前記少なくとも１つの第２の時間変動パラメータよりも早く変動する、請求項６に記載の信号処理装置。
【請求項８】
前記少なくとも１つの所望の信号は、複数の所望の信号を含み、該複数の所望の信号の各々は、それぞれの状態空間システムとしてモデリングされる、請求項１〜７のいずれか１つに記載の信号処理装置。
【請求項９】
前記複数の所望の信号の全部でなく、少なくとも１つは、少なくとも１つの第３のパラメータを用いてモデリングされ、該第３のパラメータまたは該第３のパラメータの各々は、固定され、かつ未知である、請求項８に記載の信号処理装置。
【請求項１０】
前記複数の所望の信号の全部でなく、少なくとも１つは、少なくとも１つの第４のパラメータを用いてモデリングされ、該第４のパラメータまたは該第４のパラメータの各々は、固定され、かつ公知である、請求項８または９に記載の信号処理装置。
【請求項１１】
前記第２のパラメータまたは第２のパラメータの各々は、公知である、請求項１〜１０のいずれか１つに記載の信号処理装置。
【請求項１２】
前記少なくとも１つの通信チャンネルは、複数の通信チャンネルを備え、前記少なくとも１つの汚染信号は、複数の汚染信号を含む、請求項１に記載の信号処理装置。
【請求項１３】
前記通信チャンネルの数は、所望の信号の数以上である、請求項１２に記載の信号処理装置。
【請求項１４】
前記少なくとも１つの汚染信号は、前記所望の信号の少なくともいくつかの時間遅延バージョンの線形の組み合わせを含む、請求項８に記載の信号処理装置。
【請求項１５】
前記少なくとも１つの汚染信号は、雑音で汚染された前記少なくとも１つの所望の信号を含む、請求項１〜１４のいずれか１つに記載の信号処理装置。
【請求項１６】
前記少なくとも１つのチャンネルは、異なった長さの複数の信号伝播経路を含む、請求項１〜１４のいずれか１つに記載の信号処理装置。
【請求項１７】
前記少なくとも１つの所望の信号は、アナログ信号を含む、請求項１〜１６のいずれか１つに記載の信号処理装置。
【請求項１８】
前記アナログ信号をサンプリングすることにより、時間的にサンプリングされたアナログ信号を生成する手段をさらに含む、請求項１７に記載の信号処理装置。
【請求項１９】
前記少なくとも１つの所望の信号は、少なくとも１つの音声信号を含む、請求項１７または１８に記載の信号処理装置。
【請求項２０】
前記少なくとも１つの音声信号は、音声を含む、請求項１９に記載の信号処理装置。
【請求項２１】
前記汚染信号は、複数の信号受信装置により測定される、請求項１９に記載の信号処理装置。
【請求項２２】
前記少なくとも１つの第１の時間変動パラメータは、ランダムエラーまたは外乱の時間変動パラメータを含む、請求項１９に記載の信号処理装置。
【請求項２３】
前記少なくとも１つの第１の時間変動パラメータは、音響信号の周波数スペクトルを構成するパラメータを含む、請求項１９に記載の信号処理装置。
【請求項２４】
前記少なくとも１つの所望の信号は、時間変動自己回帰としてモデリングされる、請求項１〜２３のいずれか１つに記載の信号処理装置。
【請求項２５】
前記少なくとも１つの所望の信号は、移動平均モデルとしてモデリングされる、請求項１〜２３のいずれか１つに記載の信号処理装置。
【請求項２６】
前記少なくとも１つの所望の信号は、非線形時間変動モデルとしてモデリングされる、請求項１〜２３のいずれか１つに記載の信号処理装置。
【請求項２７】
前記少なくとも１つの通信チャンネルは、有限インパルス応答モデルとしてモデリングされる、請求項１〜２６のいずれか１つに記載の信号処理装置。
【請求項２８】
前記少なくとも１つの通信チャンネルは、無限インパルス応答モデルとしてモデリングされる、請求項１〜２６のいずれか１つに記載の信号処理装置。
【請求項２９】
前記少なくとも１つの通信チャンネルは、非線形時間変動モデルとしてモデリングされる、請求項１〜２６のいずれか１つに記載の信号処理装置。
【請求項３０】
前記少なくとも１つの第１の時間変動パラメータは、確率モデルを用いてモデリングされる、請求項１〜２９のいずれか１つに記載の信号処理装置。
【請求項３１】
前記信号推定手段は、ベイズの推定により前記少なくとも１つの所望の信号をさらに推定する、請求項３０に記載の信号処理装置。
【請求項３２】
前記信号推定手段は、最尤法により前記少なくとも１つの所望の信号をさらに推定する、請求項３０に記載の信号処理装置。
【請求項３３】
前記信号推定手段は、最小２乗法により前記少なくとも１つの所望の信号をさらに推定する、請求項３０に記載の信号処理装置。
【請求項３４】
前記信号推定手段は、モンテカルロ法を前記少なくとも１つの汚染信号に適用することによって前記少なくとも１つの所望の信号をさらに推定する、請求項３１〜３３のいずれか１つに記載の信号処理装置。
【請求項３５】
前記信号推定手段は、カルマンフィルタまたは拡張カルマンフィルタにより前記少なくとも１つの所望の信号をさらに推定する、請求項３１〜３４のいずれか１つに記載の信号処理装置。
【請求項３６】
請求項１〜３５のいずれか１つに記載の信号処理装置を制御するプログラム。
【請求項３７】
請求項３６に記載のプログラムを含むコンピュータ読み取り可能な格納媒体。
【請求項３８】
請求項３６に記載のプログラムによりプログラミングされているコンピュータ。
【発明の詳細な説明】
【０００１】
本発明は、信号を抽出する方法に関する。この方法は、それぞれの通信チャンネルを介して受け取られる１つ以上の汚染信号から、１つ以上の所望の信号を抽出するために用いられ得る。信号は、例えば、ノイズ、マルチパス伝播の場合はその遅延バージョン、所望の信号であり得るまたはあり得ない他の信号、または、これらの組み合わせによって汚染され得る。
【０００２】
通信経路（単数または複数）は、ケーブルを介したり、電磁伝播および音響伝播などの任意の形態であり得る。また、所望の信号は、原則的に任意の形態であり得る。この方法の特定のアプリケーションの１つは、汚染信号（例えば、音響的に伝播するノイズまたは他の音信号など）からの音声などの音信号を抽出することが望ましいシステムである。
【０００３】
ＷＯ  ９９／６６６３８は、状態空間モデリングに基づく信号分離技術を開示する。状態空間モデルは、信号混合プロセスのために想定され、さらなる状態空間モデルは、非混合システムを意図する。信号の通信環境が時間変動であることが提案されるが、これはこの開示技術においてはモデリングされない。
【０００４】
ＵＳ  ５  ８７０  ００１は、セルラー無線システムに用いられる較正技術を開示する。この技術は、カルマンフィルタを用いる従来例である。
【０００５】
ＵＳ  ５  ８４５  ２０８は、セルラー無線システムにおいて受け取った電力を推定する技術を開示する。この技術は、パラメータが固定されておりかつ事前に推定されるまたは「既知」である状態空間自己回帰モデルを利用する。
【０００６】
ＵＳ  ５  ５８１  ５８０は、レイリーフェージング環境においてチャンネルをフェージングする、モデルに基づくチャンネル推定アルゴリズムを開示する。この推定法則は、時間変動通信チャンネル係数に自己回帰モデルを用いる。
【０００７】
ＫｏｔｅｃｈａおよびＤｊｕｒｉｃの「Ｓｅｑｕｅｎｔｉａｌ  Ｍｏｎｔｅ  Ｃａｒｌｏ  ｓａｍｐｌｉｎｇ  ｄｅｔｅｃｔｏｒ  ｆｏｒ  Ｒａｙｌｅｉｇｈ  ｆａｓｔ−ｆａｄｉｎｇ  ｃｈａｎｎｅｌ」、Ｐｒｏｃ．２０００  ＩＥＥＥ  Ｉｎｔ．Ｃｏｎｆ．Ａｃｏｕｓｔｉｃｓ  Ｓｐｅｅｃｈ  ａｎｄ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ、Ｖｏｌ．１、６１〜６４ページは、デジタルまたは離散レベルの信号を処理する技術を開示する。この技術は、動状態空間モデルとしてシステムをモデリングすることを利用する。動状態空間モデルにおいて、チャンネル推定および送信データ検出は、モンテカルロサンプリングフィルタに基づく。チャンネルフェージング係数および送信変数は、隠れ変数として扱われ、チャンネル係数は、自己回帰プロセスとしてモデリングされる。隠れ変数の粒子は、過去の観測に基づいて、インポータンスサンプリング密度から連続して生成される。これらの粒子は、必要とされる条件つきの事後分布によって伝播され、かつ、重さを計られる。これらの粒子およびその重さによって隠れ変数が推定される。この技術は、固定パラメータによるソースのモデリングに限定される。
【０００８】
Ｃｈｉｎ−Ｗｅｉ  ＬｉｎおよびＢｏｒ−Ｓｅｎ  Ｃｈｅｎの「Ｓｔａｔｅ  Ｓｐａｃｅ  Ｍｏｄｅｌ  ａｎｄ  Ｎｏｉｓｅ  Ｆｉｌｔｅｒｉｎｇ  Ｄｅｓｉｇｎ  ｉｎ  Ｔｒａｎｓｍｕｌｔｉｐｌｅｘｅｒ  Ｓｙｓｔｅｍｓ」、Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ、Ｖｏｌ．４３、Ｎｏ．１、１９９５、６５〜７８ページは、通信システムに適用される別の状態空間モデリング技術を開示する。「トランスマルチプレクサ」方式において、所望の信号は、既知かつ固定のパラメータによって非時間変動自己回帰プロセスとしてモデリングされる。
【０００９】
本発明の第１の局面によって、少なくとも１つの測定された汚染信号およびその少なくとも１つの汚染信号を測定する少なくとも１つの通信チャンネルを含むシステムから少なくとも１つの所望の信号を抽出する方法が提供され、この方法は、少なくとも１つの第１の時間変動パラメータによって第１の動状態空間システムとして少なくとも１つの所望の信号をモデリングし、かつ、少なくとも１つの第２のパラメータを有する第２の状態空間システムとして少なくとも１つの通信チャンネルをモデリングすることを含む。
【００１０】
状態空間システムは数学的に既知であり、いくつかの実用的な問題の解決策に適用されてきた。状態空間システムにおいて、システムが推定または抽出を行うことが望ましいシステムの潜在的な状態がある。この状態は、以前の状態値の既知関数および不規則なエラーまたは外乱状態として生成されることが想定される。
【００１１】
利用可能な測定値もまた、現在の状態および別の不規則なエラーまたはノイズ状態の既知関数であることが想定される。
【００１２】
状態空間システムが、１つ以上の測定された汚染信号および通信チャンネルを含むシステムから１つ以上の所望の信号を抽出する問題にうまく適用されることが不意に発見された。この技術は、扱いやすい方法で、かつ、リアルタイムまたはオンラインで、１つ以上の所望の信号を抽出することを可能にする。この技術のさらなる利点は、所望の信号（単数または複数）のサンプルを抽出するために未来のサンプルを必要としないことであるが、いくつかの実施形態において、限られた未来のサンプルを用いることが有利であり得る。後者の場合、所望の信号（単数または複数）のサンプルは幾分か遅延されているが、少なくとも測定された信号のサンプリングレートにおいて利用可能であり、かつ、非常に大量なメモリを必要としない。
【００１３】
少なくとも１つの所望の信号は、物理的プロセスによって生成され得る。少なくとも１つの第１のパラメータは、物理的生成プロセスをモデリングし得る。このようなモデリングは、一般的には実際の物理的プロセスに近いものであるが、信号抽出に十分であることが分かった。
【００１４】
少なくとも１つの第２のパラメータは、固定かつ未知であり得る。少なくとも１つの第２のパラメータは、時間変動であり得る。少なくとも１つの第１の時間変動パラメータは、少なくとも１つの第２の時間変動パラメータの変化率とは異なる変化率を有し得る。少なくとも１つの第１の時間変動パラメータの変化率は、平均して、少なくとも１つの第２の時間変動パラメータの変化率より大きくあり得る。
【００１５】
多くのシステムにおいて、所望の信号（単数または複数）の特性は比較的急速に変化するが、通信チャンネル（単数または複数）の特性はより遅く変化する。チャンネルの特性は突然変化し得るが、このような変化は比較的珍しい。しかし、音声などの信号は、比較的急速に変化する特性を有する。異なる割合の変化がモデリングされるようにこれらの特性をモデリングすることにより、１つ以上の信号の抽出が促進される。
【００１６】
少なくとも１つの所望の信号は、複数の所望の信号を含み得、それぞれの所望の信号は、それぞれの状態空間システムとしてモデリングされる。複数の所望の信号の少なくとも１つ（全てではない）は、少なくとも１つの第３のパラメータによってモデリングされ得、第３のパラメータ（単数または複数）は、固定かつ未知である。
【００１７】
複数の所望の信号の少なくとも１つ（全てではない）は、少なくとも１つの第４のパラメータによってモデリングされ得、第４のパラメータ（単数または複数）は既知である。
【００１８】
第２のパラメータ（単数または複数）は既知であり得る。
【００１９】
少なくとも１つの通信チャンネルは、複数の通信チャンネルを含み得、少なくとも１つの汚染信号は、複数の汚染信号を含み得る。少なくとも１つの所望の信号は、複数の所望の信号を含み得る。通信チャンネルの数は、所望の信号の数以上であり得る。必要ではないが、一般的に、測定された信号の数が、抽出されるべき所望の信号の数以上であることが望ましい。これにより、この方法によって所望の信号が再現される効率が改善され、特に、所望の信号（単数または複数）の再構築または抽出の正確さが改善される。
【００２０】
少なくとも１つの汚染信号は、少なくともいくつかの所望の信号の時間遅延バージョンの一次結合を含み得る。この方法は、従って、マルチパス伝播の場合、信号が互いに汚染し合っている場合、およびこれらの効果の組み合わせの場合に、所望の信号を抽出することが可能である。
【００２１】
少なくとも１つの汚染信号は、ノイズによって汚染されている少なくとも１つの所望の信号を含み得る。従って、この方法は、所望の信号（単数または複数）をノイズから抽出することが出来る。少なくとも１つのチャンネルは、異なる長さの複数の信号伝播経路を含み得る。
【００２２】
少なくとも１つの所望の信号は、アナログ信号を含み得る。アナログ信号は、一時的にサンプリングされたアナログ信号であり得る。
【００２３】
少なくとも１つの所望の信号は、音信号を含み得る。少なくとも１つの音信号は、音声を含み得る。汚染信号は、音領域を空間的にサンプリングすることにより測定される。少なくとも１つの第１のパラメータは、ノイズ生成モデリングパラメータを含み得る。少なくとも１つの第１のパラメータは、ホルマントモデリングパラメータを含み得る。例えば、マイクロフォンなどの音響電気変換器が、例えば、部屋または他の空間に空間的に分布され得、出力信号が、この方法によって処理され得、その結果、音声の他のソースまたは他の情報を持つ音のソースなどのバックグラウンドノイズまたは信号の存在下で、１つのソースから音声が抽出または分離される。
【００２４】
少なくとも１つの所望の信号は、時間変動自己回帰としてモデリングされ得る。この型のモデリングは、多くの種類の所望の信号に適しており、音声を抽出することに特に適している。また、少なくとも１つの所望の信号は、移動平均モデルとしてモデリングされ得る。さらにまた、少なくとも１つの所望の信号は、非線形時間変動モデルとしてモデリングされ得る。
【００２５】
少なくとも１つの通信チャンネルは、時間変動有限インパルス応答モデルとしてモデリングされ得る。この型のモデルは、種々の伝播システムをモデリングすることに適している。また、少なくとも１つの通信チャンネルは、有限インパルス応答モデルとしてモデリングされ得る。さらにまた、少なくとも１つの通信チャンネルは、非線形時間変動モデルとしてモデリングされ得る。
【００２６】
第１の状態空間システムは、確率モデルを用いてモデリングされる少なくとも１つのパラメータを有し得る。少なくとも１つの所望の信号は、ベイズの反転によって抽出され得る。また、少なくとも１つの所望の信号は、最尤法によって抽出され得る。さらにまた、少なくとも１つの所望の信号は、最小２乗法によって抽出され得る。信号抽出は、連続的なモンテカルロ方法によって実行され得る。また、信号抽出は、カルマンフィルタまたは拡張されたカルマンフィルタによって実行され得る。これらの技術は、複雑なモデルに特に有効であり、並列コンピュータに潜在的に実現可能である。
【００２７】
本発明の第２の局面によって、前述の請求項のいずれかに記載される方法を実行するために、コンピュータを制御するプログラムが提供される。
【００２８】
本発明の第３の局面によって、本発明の第２の局面によるプログラムを含むキャリアが提供される。
【００２９】
本発明の第４の局面によって、本発明の第２の局面によるプログラムによってプログラミングされたコンピュータが提供される。
【００３０】
例として、添付の図面を参照して、本発明をさらに説明する。
【００３１】
所望の信号を抽出するために、パラメトリックアプローチが用いられる。このパラメトリックアプローチにおいて、データは、時間ｔ、変数ｘｔによって表される潜在的な観測されていないプロセスによって生成されると想定される。変数ｘｔは、異なるソースの波形に関する情報を含む。所望の信号を抽出する際の問題は、問題の実際の構造がほとんど分かっていないため、良好な決定論的モデルを利用するには複雑すぎることである。あるいは、良好な決定論的モデルが利用可能である場合、扱い切れないほどの式が導かれ得ることである。
【００３２】
音を生成するためのモデルは、以前の時間ステップｘｔ−１における状態の値から、現在の「状態」ｘｔの尤度分布の式を提供する。この確率分布（状態遷移
0004402879.tif 000002
現在観測されているデータが現在の状態にどのように依存するかは、別の確率分
0004402879.tif 000003
る。
最終的には、最初の時刻において状態がどのように分布されるかは、ｐ（ｘ０）によって特定される。これらの３つの分布の特定の形態は、後で得られる。上記の問題を回避するための解決策は、確率分布によって式に不確定さを取り入れることを含む。より正確には、離散的な時間のセットアップにおいてｘｔ＋１が過去の値の決定論的な関数である（例えば、Ａが線形演算子である場合、ｘｔ＋１＝Ａｘｔ）と想定するのではなく、パラメータが存在する可能性の高い状態
0004402879.tif 000004
ｘｔ＋１の確率は、以前の値ｘｔからのχに等しい。これは、例えば、ｖｔがガウス分布によってゼロ付近に分布される不確定さまたはエラーである場合、ｘｔ＋１＝Ａｘｔ＋ｖｔと表され得る。分布の広がり方は、決定論的成分Ａｘｔの信用度を示す。広がりが狭いほど、信用は大きい。この型のモデリングは、非常に複雑なプロセスを表すために実際にしっかりしていることが実証される一方で、実際に使用されるのに十分なほど簡単である。
【００３３】
上記のように、プロセスの構造は既知であることが（可能性の高い物理的想定から）想定されるが、変数ｘｔは観測されず、ｙｔだけが観測可能である。概して、観測構造、すなわち、観測される前に対象のプロセスが経験する変形は、変数ｘｔによるが、現象の説明に幾分の不規則性が取り入れられる必要がある。例
0004402879.tif 000005
ｘｔ）によって行われる。潜在的なプロセスのパラメータがｘｔによって表されることから、ｙｔの確率はｙに等しい。上記の例において、可能な観測プロセスは、ｙｔ＝Ｃｘｔ＋ｗｔの形態であり、付加的な観測ノイズｗｔによって破壊されるプロセスの内部状態を表すパラメータｘｔが幾分変形する。
【００３４】
音響ソース分離の場合、パラメータｘｔは、時間ｔにおけるソースの所望の波形の値ｓｔを含み、ソース（ａｔ）および混合システム（ｈｔ）を導き出すためのパラメータを含む。図１は、単一の通信チャンネルに関して、これを図によって示す。遷移確率分布は、下記によって定義される。
ｐ（ａｔ＋１＝ａ│ａｔ）ｐ（ｈｔ＋１＝ｈ│ｈｔ）ａｔがｈｔ（突然導き出されることもあるが、これはめったに起こらない）に比べて素早く導き出されることを除いて、これの展開をどのように特徴付けるべきか分からないため、下記が与えられる。
【００３５】
【数１】
0004402879.tif 000006
【００３６】
ここで、
【００３７】
【数２】
0004402879.tif 000007
【００３８】
の分布はそれぞれ等方的であるが、異なる非定常時間尺度を反映するために異なる広がりを有する。パラメータａｔのセットは、ソースのスペクトル特性を導き出すと考えられ得、パラメータｈｔは、古典有限インパルス応答（ＦＩＲ）フィルタと考えられ得る。システムのパラメータが説明されたので、その波形が確率分布によって導き出されることが想定され得る。
ｐ（ｓｔ＋１＝ｓ│ｓｔ，ａｔ）これは、ソースのモデルを説明する。
【００３９】
音響信号の分離または抽出の一例として、ソースは音声ソースである。パラメータｘｔは、次いで、時間変動自己回帰プロセスとして、例えば、下記の式によって、モデリングされ得る。
【００４０】
【数３】
0004402879.tif 000008
【００４１】
ここで、パラメータａｉ，ｔは、人間の声の広がりを表す全共振フィルタシステムのフィルタ係数であり、ｅｔは、声帯によって生成されるノイズを表す。
【００４２】
図２は、会話をしている人などの２つの所望の音ソースが、ｓ（１）ｔおよびｓ（２）ｔによって表され、壁２〜５の形態で境界を有する部屋１の形態の密閉された空間に配置される本発明の典型的なアプリケーションを示す。マイクロフォン６〜８は、部屋１の中の固定されたそれぞれの位置に配置され、部屋の中の音領域をサンプリングして測定された汚染信号ｙ（１）ｔ、ｙ（２）ｔ、ｙ（３）ｔを生成する。測定された信号は、抽出方法を実行するためにプログラムによって制御されているコンピュータ９に供給される。このプログラムは、プログラムキャリア１０が保有し、図示の例において、プログラムキャリア１０は、コンピュータ９の操作を制御するメモリの形態である。このコンピュータは、個々の音声信号を抽出し、増幅器１１および１２ならびに拡声器１３および１４を含む個々のチャンネルを介してこれらの音声信号を供給するように示されている。その代わり、または、それに加えて、抽出された信号は格納され得るか、または、さらなる処理の対象となり得る。
【００４３】
図２は、音ソースから１つのマイクロフォン７へのいくつかの伝播経路を示す。それぞれの場合において、それぞれの音ソースからマイクロフォン７へは、直接経路１５および１６がある。それに加えて、それぞれの音ソースからマイクロフォン７への反射経路がある。例えば、１つの反射経路は、ソースから出て壁２で反射して反射線１８を形成する直接音線１７を含む。概して、それぞれの音ソースからそれぞれのマイクロフォンへの多くの反射経路があり、これらの経路は、異なる伝播の長さおよび異なる音減衰を有する。
【００４４】
観測ｙｔからソースｓｔを推定することが目的であり、実用的であるために、推定は、出来る限り遅延の少ない状態で、オンラインで実行されることが望ましい。モデリングプロセスは、効果（観測ｙｔ）をもたらす原因（ｘｔ）について説明するが、推定タスクは、原因および効果の反転を含む。確率モデリングのフレームワークにおいて、この反転は、ベイズの定理を用いて一貫して実行され得る。統計的フィルタリングに関して、統計的フィルタリングは時間ｔまでの全て
0004402879.tif 000009
１：ｔ）である。展開確率（ｅｖｏｌｕｔｉｏｎ  ｐｒｏｂａｂｉｌｉｔｙ）が
0004402879.tif 000010
イズの法則は下記をもたらす。
【００４５】
【数４】
0004402879.tif 000011
【００４６】
この式によって、時間ｔ−１におけるフィルタリング密度、すなわちｐ（ｘｔ
0004402879.tif 000012
非公開の形態でこの積分を計算することが出来ず、これらの数をオンラインのリアルタイムな方法で計算することが望ましいことである。
【００４７】
確率分布を用いると、機能空間の特定領域内の人口のメンバー密度が低く表されると考えられ得るが、モンテカルロ方法は、反対の方向に作用し、確率分布によって分布される意図的に生成されたサンプルのセットによって確率分布が表され得るという概念による。これは、機能空間の所定区域内のサンプル密度が、対象の分布において、この区域の確率を表すことが想定されることを必要とする。予想通り、人口が多いほど、表現はより正確になる。このアプローチは、多くの場合、同時に行われることが大いにあり得る計算を大幅に単純化するという利点を有する。
【００４８】
本発明に関与する技術の基本的な説明を以下に示す。その後、特定の実施形態の詳細な説明を行う。
【００４９】
0004402879.tif 000013
「粒子」
【００５０】
【数５】
0004402879.tif 000014
【００５１】
0004402879.tif 000015
−１）
または事前確率を用いて、粒子がどこで導き出されるかを推測することが可能である。展開分布または事前確率からのサンプリングは典型的に簡単である。従って、パラメータの展開の事前確率によって、尤もな領域に「スカウト」が
0004402879.tif 000016
れ得る新しい観測によってもたらされる情報を「スカウト」が考慮していないため、予測を訂正する必要がある。いくつかの粒子は対象の領域内にあるが、典型的にその数は不十分である。または、より少ないメンバーを含むべき領域に人口のメンバーが多すぎる場合もある。
【００５２】
人口を調整する方法は、人口不足領域内のメンバーを増やすこと、および、人口過密領域内のメンバーを抑制することからなる。このプロセスは、選択ステップと呼ばれる。「スカウト」の未来を決定する量が下記のインポータンス関数であることが、数学的に証明され得る。
【００５３】
【数６】
0004402879.tif 000017
【００５４】
有効なメカニズムは、
【００５５】
【数７】
0004402879.tif 000018
【００５６】
に最も近い整数を用いてスカウトの数ｉの内の「子供」の数を決定すること、または、確率
【００５７】
【数８】
0004402879.tif 000019
【００５８】
によって粒子
【００５９】
【数９】
0004402879.tif 000020
【００６０】
0004402879.tif 000021
って子供がおおよそ分布され、次のデータが処理され得る。これは、アルゴリズムの非常に一般的な説明であり、非常に特別な場合において、下記に示すように、多くの改善が可能である。
【００６１】
連続的なモンテカルロ推定を行う一般的な手法が与えられる。これは、最も基本的な形態の方法であり、Ｎ．Ｊ．Ｇｏｒｄｏｎ、Ｄ．Ｊ．ＳａｌｍｏｎｄおよびＡ．Ｆ．Ｍ．Ｓｍｉｔｈの「Ｎｏｖｅｌ  ａｐｐｒｏａｃｈ  ｔｏ  ｎｏｎｌｉｎｅａｒ／ｎｏｎ−Ｇａｕｓｓｉａｎ  Ｂａｙｅｓｉａｎ  ｓｔａｔｅ  ｅｓｔｉｍａｔｉｏｎ」、ＩＥＥ  Ｐｒｏｃｅｅｄｉｎｇｓ−Ｆ、ｖｏｌ．１４０、ｎｏ．２、１０７〜１１３ページ、１９９３に記載されており、本明細書中、同文献の内容を参考として援用する。
【００６２】
初期ステップ２０において、時間ｔはゼロに設定される。Ｎいわゆる「粒子」は、初期分布ｐ（ｘ０）から不規則に選択され、ｘ（ｋ）０と表示される。ここで、ｋは１からＮの整数である。粒子の数Ｎは、典型的に非常に大きく、例えば、ほぼ１，０００である。ステップ２０および２１は、手順の初期化を表す。
【００６３】
状態伝播がステップ２２において実行される。このステップにおいて、現在の
0004402879.tif 000022
現在の時間への更新が不規則に選択される。更新された粒子ｘ（ｋ）ｔ＋１は、状態遷移分布によって、空間の期待領域内にある尤もな値を有する。
【００６４】
ステップ２３において、次に測定される値ｙｔ＋１は、音領域を測定またはサンプリングすることによって得られる。ステップ２４において、粒子の重量が計算される。更新された粒子ｘ（ｋ）ｔ＋１が、新しく測定された値ｙｔ＋１を参照せずに生成されたので、新しい測定を考慮して各粒子が実際にどれほど「良い」かを表すために、各粒子の０と１との間の重量を計算する必要がある。各粒子の正確な重量の値は、粒子の観測分布の値に比例する。しかし、全ての重量の合計は、１に等しい必要がある。従って、各粒子ｘ（ｋ）ｔ＋１の重量は、下記によって得られる。
【００６５】
【数１０】
0004402879.tif 000023
【００６６】
ステップ２５は、次いで、粒子を再び選択またはサンプリングして、付加重量のない粒子のセットに戻る。ステップ２５は、以下により詳細に説明するように、その重量によって、Ｎ粒子を再び選択する。ステップ２６は、次いで、時間を増やして、ステップ２２に戻るように制御する。従って、ステップ２２〜２６は、この方法が用いられる限り繰り返される。
【００６７】
各サンプルの値は、例えば、事後平均推定法則によって、下記のように計算され得る。
【００６８】
【数１１】
0004402879.tif 000024
【００６９】
上記の技術は、比較的一般的であり、「汚染された」測定から音声または他の音響信号を分離または抽出する場合に特に限定されない。音響ソース信号の抽出をより明確にその目的とする方法について、より詳細に説明する。
【００７０】
ソースは、時間変動自己回帰プロセスによってモデリングされることが想定される。すなわち、時間ｔにおけるソースの数ｉは、同じソース（ｓｉ，ｔ−１、．．．、ｓｉ，ｔ−ｐｉ、または短いベクトルの形態ｓｉ，ｔ−１：ｔ−ｐｉ）の過去の値であるｐｉ（いわゆるモデルのオーダ）の一次結合であり、ノイズ（ガウス分布されていると想定される）
【００７１】
【数１２】
0004402879.tif 000025
【００７２】
によって摂動する。
【００７３】
【数１３】
0004402879.tif 000026
【００７４】
この種類のプロセスは、非常に順応性があり、共振がモデリングされることを可能にする。一次結合の係数ａｉ，ｔは、音声の共振の非定常を考慮するため、時間変動である。
【００７５】
【数１４】
0004402879.tif 000027
【００７６】
異なるソース（遅延された可能性のある）をマイクロフォンｊおよび時間ｔにおいて重ね合わせる観測は、下記のプロセスによって生成されることが想定される。
【００７７】
【数１５】
0004402879.tif 000028
【００７８】
すなわち、この観測は、観測ノイズｗｊ，ｔによって摂動するフィルタリングされた全てのソースの合計である（ソースｉからマイクロフォンｊへのフィルタの長さはｌｉｊであり、このフィルタは遅延されている）。ソースｉからマイクロフォンｊへの伝達関数は、下記の式によって、そのうち導き出される。
【００７９】
【数１６】
0004402879.tif 000029
【００８０】
ここで、
【００８１】
【数１７】
0004402879.tif 000030
【００８２】
およびｗｊ，ｔの特性が既知であることが想定されるが、アルゴリズムの全体のバージョンにおいて、これらのパラメータも推定される。
【００８３】
音響ソース分離のための連続的なモンテカルロアルゴリズムのステップは、基本的には、図３に示す通りであるが、下記のように変更または増強される。状態ベクトルｘｔは、全てのソース、自己回帰パラメータ、およびソースとマイクロフォンとの間の伝達関数を含む積重ねベクトルとして定義される。
【００８４】
初期値は、広がりの大きい通常分布から不規則に選択される。
【００８５】
ステップ２２において、ガウス分布から、不規則ノイズｖａ（ｋ）ｉ，ｔ＋１、ｖｓ（ｋ）ｉ，ｔ＋１およびｖｈ（ｋ）ｉ，ｊ，ｔ＋１が生成される。ガウス分布は、例えば、Ｂ．Ｄ．Ｒｉｐｌｅｙの「Ｓｔｏｃｈａｓｔｉｃ  Ｓｉｍｕｌａｔｉｏｎ」、Ｗｉｌｅｙ、Ｎ．Ｙ．１９８７に開示され、本明細書中、同文献の内容を参考として援用する。これらの量は、次いで、それぞれの状態変数に加算される。
【００８６】
【数１８】
0004402879.tif 000031
【００８７】
ステップ２４は、下記のように、非正常化重量を計算する。
【００８８】
【数１９】
0004402879.tif 000032
【００８９】
重量は、次いで、下記の式によって再び正常化される。
【００９０】
【数２０】
0004402879.tif 000033
【００９１】
励起ノイズおよび観測ノイズの分散の変動は、
【００９２】
【数２１】
0004402879.tif 000034
【００９３】
の展開式を下記のように定義することにより考慮され得る。
【００９４】
【数２２】
0004402879.tif 000035
【００９５】
ここで、
【００９６】
【数２３】
0004402879.tif 000036
【００９７】
は、ガウス分布によって分布されるｉ．ｉ．ｄ．シーケンスである。
【００９８】
次いで、それぞれの修正のいくつかまたは全ての性能を基本的なアルゴリズムに改善することが、下記のように可能である。
【００９９】
（ａ）時間変動ノイズ分散を推定する。不規則ノイズ
【０１００】
【数２４】
0004402879.tif 000037
【０１０１】
を通常分布から生成し、下記を計算する。
【０１０２】
【数２５】
0004402879.tif 000038
【０１０３】
上記の式は、指数関数によって変換されて、ノイズの分散を得る。すなわち、
【０１０４】
【数２６】
0004402879.tif 000039
【０１０５】
である。
【０１０６】
（ｂ）このモデルの構造を用いることにより、混合フィルタ係数および自己回帰係数が積算され得る。その含意は、重量の計算が下記のように修正されることである。
【０１０７】
【数２７】
0004402879.tif 000040
【０１０８】
この重量は、下記の式（１７）〜（１９）に示される自己回帰係数および混合フィルタ係数を含む状態の状態空間モデルに適用される標準的な手順であるカルマンフィルタ（式（５０）を参照）を用いて連続的に計算される。この方法の利点は、推定されるべきパラメータの数が著しく減少され、統計的な効率がとても改善されることである。
【０１０９】
（ｃ）ソースの値を伝播する分布は、フィルタリング密度
【０１１０】
【数２８】
0004402879.tif 000041
【０１１１】
の近似に修正される。このフィルタリング密度は、基本的なアルゴリズムとは対照的に、新しい観測ｙｔを考慮し、従って、改善された統計的な効率に再び導く。この密度は、状態空間モデルに適用されるカルマンフィルタの副産物（例えば、カルマンフィルタ＃２）である。この状態空間モデルにおいて、状態は、カルマンフィルタ＃１によって以前に求められた混合係数
【０１１２】
【数２９】
0004402879.tif 000042
【０１１３】
および自己回帰係数
【０１１４】
【数３０】
0004402879.tif 000043
【０１１５】
のフィルタリングされた平均推定によって、ソースｘｔおよび下記の式（１４）に示されるパラメータからなる。
【０１１６】
（ｄ）粒子間の多様性は、対象分布が
【０１１７】
【数３１】
0004402879.tif 000044
【０１１８】
である粒子にメトロポリス−ハスティング更新を用いて導入される。このステップは、ステップ２５の後に導入され、その詳細は以下に示される。
【０１１９】
（ｅ）改善のため、すなわち、時間ｔにおけるソースの値を計算していくつかのＬに関するいくつかの未来観測ｙｔ＋１、．．．、ｙｔ＋Ｙを待って考慮するため、推定が遅延される。このような技術は、固定遅延平滑化と呼ばれ、データｙｔ＋１を待つことのみが必要であるため、アルゴリズムは修正されない。
【０１２０】
これらの異なるステップの完全な説明を以下に示す。
【０１２１】
ソースおよびチャンネルを含むシステム全体は、状態空間システムとしてモデリングされる。これは、制御および信号処理の理論による標準的な概念である。状態空間システムの定義として、システムには、システムが推定または抽出を行うことが望ましい、ｘｔと表される、時間ｔにおけるいくつかの潜在的な状態がある。この場合、潜在的な所望のソース自体が含まれる。この状態は、以前の状態値の既知関数および不規則エラーまたは外乱状態として生成されると想定される。
【０１２２】
【数３２】
0004402879.tif 000045
【０１２３】
ここで、Ａｔ＋１（．．．）は既知関数であり、ある期間にわたる状態の想定された変遷パターンを表す。同様に、ｙｔと表される時間ｔにおける測定値は、現在の状態および別の不規則エラーまたはノイズ状態の既知関数であることが想定される。
【０１２４】
【数３３】
0004402879.tif 000046
【０１２５】
ここで、
【０１２６】
【数３４】
0004402879.tif 000047
【０１２７】
は、汚染プロセスを表す既知関数である。この場合、
【０１２８】
【数３５】
0004402879.tif 000048
【０１２９】
は、ソース（単数または複数）上のチャンネル（単数または複数）のフィルタリング効果を表し、
【０１３０】
【数３６】
0004402879.tif 000049
【０１３１】
は、システム内の測定ノイズを表す。
【０１３２】
本方法の重要な鍵は、ソース（単数または複数）とチャンネル（単数または複数）の両方が有する、時間変動特性が異なることであり、これは、さらに推定される必要がある。換言すると、関数Ａｔ＋１および
【０１３３】
【数３７】
0004402879.tif 000050
【０１３４】
は、それ自体、例えば、
【０１３５】
【数３８】
0004402879.tif 000051
【０１３６】
といった、いくつかの未知のパラメータに依存する。この場合、
【０１３７】
【数３９】
0004402879.tif 000052
【０１３８】
は、ソースの未知の時間変動自己回帰パラメータを表し、
【０１３９】
【数４０】
0004402879.tif 000053
【０１４０】
は、チャンネル（単数または複数）の未知の時間変動有限インパルス応答フィルタ（単数または複数）を表す。これらの時間変動特性は、後述のように、さらなる状態遷移関数を用いてモデリングされる。
【０１４１】
対処される問題は、ソース分離の問題、すなわち、ｎ個のソースは自己回帰（ＡＲ）プロセスとしてモデリングされ、このソースから、各時間ｔにおいて、ｎ個のソースの畳み込み混合（ｃｏｎｖｏｌｕｔｉｖｅ  ｍｉｘｔｕｒｅｓ）であるｍ回の観測が行なわれる。
ソースｉは、ｔ＝１，．．．の場合、
【０１４２】
【数４１】
0004402879.tif 000054
【０１４３】
としてモデリングされ得、ｐｉは、ｉ番目のＡＲモデルの次数である。ｉ＝１，．．．，ｎである場合、
【０１４４】
【数４２】
0004402879.tif 000055
【０１４５】
を想定する。
【０１４６】
【数４３】
0004402879.tif 000056
【０１４７】
は、ゼロ平均正規化ｉ．ｉ．ｄ．ガウスシーケンスである。ｉ＝１，．．．，ｎおよびｔ＝１，．．．である場合、ガウスシーケンスは、
【０１４８】
【数４４】
0004402879.tif 000057
【０１４９】
である。
【０１５０】
【数４５】
0004402879.tif 000058
【０１５１】
は、時間ｔにおける、各ソースのダイナミック雑音の分散である。導き出される自己回帰モデルは、線形ガウス状態空間表現
【０１５２】
【数４６】
0004402879.tif 000059
【０１５３】
に従うことが想定される。このモデルは、行列
【０１５４】
【数４７】
0004402879.tif 000060
【０１５５】
の異なったセット間の交換を含み得、例えば、沈黙を考慮に入れる。通常、
【０１５６】
【数４８】
0004402879.tif 000061
【０１５７】
である。
【０１５８】
混合モデルは、多次元の時間変動ＦＩＲフィルタであると想定される。ソースは以下の態様で混合され、追加ガウスｉ．ｉ．ｄ．雑音シーケンスにより破損されることが想定される。ｊ番目のセンサにおいて、およびｊ＝１，．．．，ｍの場合、
【０１５９】
【数４９】
0004402879.tif 000062
【０１６０】
であり、ここで、ｌｉｊは、ソースｉからセンサｊまでの長さである。級数（Ｗｊ，ｔ）ｔ＝１．．．は、ゼロ平均正規化ｉ．ｉ．ｄ．である。すなわち、ｊ＝１，．．．，ｎおよびｔ＝１，．．．の場合、ガウスシーケンスは
【０１６１】
【数５０】
0004402879.tif 000063
【０１６２】
である。
【０１６３】
量
【０１６４】
【数５１】
0004402879.tif 000064
【０１６５】
は、時間ｔにおける、各センサの観測雑音の分散である。Ｗｊ，ｔは、ＡＲモデルの励起から独立していることが想定される。拘束条件［ｈｉ，ｊ，ｔ］ｌ，ｌ＝１および［ｈｉ，ｊ，ｔ］ｋ，ｌ＝０（この場合、ｊ＝ｉ  ｍｏｄ  ｍおよびｉ＝１，．．．，ｎ）が課される。ｍ＝ｎの場合、この拘束条件は、［ｈｉ，ｊ，ｔ］ｌ，ｌ＝１および［ｈｉ，ｊ，ｔ］ｋ，ｌ＝０（この場合、ｊ＝ｉ）に対応する。ソースのモデルに関して、観測システムは、さらに、状態空間表現に従うことが想定される。
【０１６６】
【数５２】
0004402879.tif 000065
【０１６７】
ソースｉの各々に関して、信号は、以下の式で書き換えられる。
【０１６８】
【数５３】
0004402879.tif 000066
【０１６９】
ここで、
【０１７０】
【数５４】
0004402879.tif 000067
【０１７１】
はλｉ×λｉ行列であり、
【０１７２】
【数５５】
0004402879.tif 000068
【０１７３】
と定義される。
【０１７４】
動的システム式は、以下のように書き変えられ得、
【０１７５】
【数５６】
0004402879.tif 000069
【０１７６】
である。
「積み重ね（ｓｔａｃｋｅｄ）」パラメータベクトルを定義すると、
【０１７７】
【数５７】
0004402879.tif 000070
【０１７８】
であり、以下の状態空間表現を考慮することが可能である。
【０１７９】
【数５８】
0004402879.tif 000071
【０１８０】
である。
【０１８１】
これは、分散および状態Ｘｔを条件として、システム（１７）および（１９）が線形ガウス状態空間モデルであるため、混合フィルタおよび自己回帰フィルタが積分されることを可能にするので、実際に役立つ利点である。（１４）と共に、これらのシステムは、双線形ガウスプロセスを定義する。
【０１８２】
ソースの数をｍ、ｐｉおよびｌｉ，ｊとすると、ソース（Ｘｔ）ｔ＝１，．．．およびこれらのソースのパラメータ
【０１８３】
【数５９】
0004402879.tif 000072
【０１８４】
を、観測ｙｊ，ｔから逐次的に推定することが必要とされる。より厳密には、ベイズの推定の枠組みで、ｐ（ｄθｔ，ｄｘｔ｜ｙ１：ｔ＋Ｌ）型の事後分布（ｐｏｓｔｅｒｉｏｒ  ｄｉｓｔｒｉｂｕｔｉｏｎｓ）の再帰的、かつ時間どおりの推定に関心が持たれる。Ｌ＝１であるとき、これはフィルタ分布に対応し、Ｌ＞０であるとき、これは固定遅延平滑化分布（ｆｉｘｅｄ−ｌａｇ  ｓｍｏｏｔｈｉｎｇ  ｄｉｓｔｒｉｂｕｔｉｏｎ）に対応する。これはいかなる解析的解決も認めない、非常に複雑な問題であり、数値的方法（ｎｕｍｅｒｉｃａｌ  ｍｅｔｈｏｄｓ）を手段とすることを必要とする。このような数値的方法は、モンテカルロシミュレーションに基づく。次に、以下の表記法が用いられる。
【０１８５】
【数６０】
0004402879.tif 000073
【０１８６】
シミュレーションベースの最適フィルタ／固定遅延平滑子が用いられ、不観測ソースのフィルタリングされた／固定遅延平滑化推定値、および
【０１８７】
【数６１】
0004402879.tif 000074
【０１８８】
型のこれらの推定値のパラメータを取得する。
【０１８９】
標準的なベイズのインポータンスサンプリング法が最初に説明され、次に、カルマンフィルタに関連するアルゴリズムを用いて、高次元であり得るパラメータａｔおよびｈｔに積算することにより、この方法の解析構造を利用することはどのように可能であるかが示される。これは、明解で効率的なアルゴリズムに至る。このアルゴリズムにとって、唯一の軌道パラメータ（ｔｒａｃｋｅｄ  ｐａｒａｍｅｔｅｒ）は、ソースおよび雑音分散である。その後、最適フィルタリングのベイズインポータンスサンプリングの逐次バージョンが提示され、プロセスに選択およびダイバーシティを導入することが必要な理由が示される。最後に、モンテカルロフィルタ／固定遅延平滑子が説明される。
【０１９０】
0004402879.tif 000075
ばれるＮのｉ，ｉ，ｄ．サンプル
【０１９１】
【数６２】
0004402879.tif 000076
【０１９２】
をサンプリングすることが可能であると仮定すると、この分布の経験的推定は
【０１９３】
【数６３】
0004402879.tif 000077
【０１９４】
により与えられ、従って、周辺分布ｐ（ｄＸｔ，ｄθｔ｜ｙ１：ｔ＋Ｌ）のモンテカルロ近似
【０１９５】
【数６４】
0004402879.tif 000078
【０１９６】
になる。
【０１９７】
この分布を用いて、任意のｆｔに関するＩＬ（ｆｔ）の推定値
【０１９８】
【数６５】
0004402879.tif 000079
【０１９９】
が取得され得る。
【０２００】
この推定は不偏であり、大数の強法則
【０２０１】
【数６６】
0004402879.tif 000080
【０２０２】
による。さらなる想定のもとに、推定は、中心極限定理を満たす。任意のｆｔに関してＩＬ（ｆｔ）を推定することは容易であり、この推定の収束率（ｒａｔｅ  ｏｆ  ｃｏｎｖｅｒｇｅｎｃｅ）は、ｔまたは状態空間の特性に依存しないが、粒子Ｎの数および関数ｆｔの特性にのみ依存する。残念ながら、任意のｔにおける分布ｐ（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）から直接的にサンプリングすることは不可能であり、代替的な策が研究される必要がある。
【０２０３】
ｐ（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）およびＩＬ（ｆｔ）を推定する１つの解決策は、Ａ．Ｄｏｕｃｅｔ，Ｓ．Ｊ．ＧｏｄｓｉｌｌとＣ．Ａｎｄｒｉｅｕらによる「Ｏｎ  ｓｅｑｕｅｎｔｉａｌ  Ｍｏｎｔｅ  Ｃａｒｌｏ  ｓａｍｐｌｉｎｇ  ｍｅｔｈｏｄｓ  ｆｏｒ  Ｂａｙｅｓｉａｎ  ｆｉｌｔｅｒｉｎｇ」Ｓｔａｔｉｓｔｉｃｓ  ａｎｄ  Ｃｏｍｐｕｔｉｎｇ、２０００年、において開示された周知のベイズのインポータンスサンプリング法である。この方法は、任意のインポータンス分布π（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ０：ｔ＋Ｌ）の存在を想定する。この分布は容易にシミュレーションされ、この分布の支持（ｓｕｐｐｏｒｔ）は、ｐ（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）の支持を含む。この分布ＩＬ（ｆｔ）を用いて、
【０２０４】
【数６７】
0004402879.tif 000081
【０２０５】
として表現され得、ここで、インポータンス重みＷ（ｘ０：ｔ＋Ｌ，θ０：ｔ＋Ｌ）が
【０２０６】
【数６８】
0004402879.tif 000082
【０２０７】
によって与えられる。
【０２０８】
インポータンス重みは、通常、比例関係の定数までのみ評価され得る。なぜなら、ベイズの法則の結果として
【０２０９】
【数６９】
0004402879.tif 000083
【０２１０】
となるからであり、
ここで、正規化定数ｐ（ｙ１：ｔ＋Ｌ）は、通常、閉じた形式（ｃｌｏｓｅｄ−ｆｏｒｍ）で表され得ない。
【０２１１】
Ｎのｉ．ｉ．ｄ．サンプル
【０２１２】
【数７０】
0004402879.tif 000084
【０２１３】
が分布π（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）によりシミュレーションされ得る場合、（２５）におけるＩＬ（ｆｔ）のモンテカルロ推定は、
【０２１４】
【数７１】
0004402879.tif 000085
【０２１５】
として取得され得、
ここで、正規化インポータンス重みは、
【０２１６】
【数７２】
0004402879.tif 000086
【０２１７】
によって与えられる。
【０２１８】
この方法は、式
【０２１９】
【数７３】
0004402879.tif 000087
【０２２０】
の目標分布の点量近似（ｐｏｉｎｔ  ｍａｓｓ  ａｐｐｒｏｘｉｍａｔｉｏｎ）と等価である。
【０２２１】
完璧なシミュレーションケース、すなわち、π（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）＝Ｐ（ｄｘ０：ｔ＋Ｌ，ｄθ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）であるとき、
【０２２２】
【数７４】
0004402879.tif 000088
【０２２３】
に対応する。実際には、インポータンス分布は、所与の意味で目標分布に可能な限り近くなるように選択される。有限Ｎである場合、
【０２２４】
【数７５】
0004402879.tif 000089
【０２２５】
は偏りがある。なぜなら、大数
【０２２６】
【数７６】
0004402879.tif 000090
【０２２７】
の強法則によると、これは推定の比率を含むが、比対称であるからである。さらなる想定のもとに、Ｄｏｕｃｅｔ，ＧｏｄｓｉｌｌおよびＡｎｄｒｉｅｕらにより開示されるように、中心極限定理がさらに成り立つ。
【０２２８】
ｐ（ｄｘｔ，ｄθｔ｜ｙｔ：ｔ＋Ｌ）およびＩＬ（ｆｔ）の推定をｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）からのサンプリングの１つに制限することが可能であり、ここで、
【０２２９】
【数７７】
0004402879.tif 000091
【０２３０】
であることが回想される。実際、
ｐ（ｄαｔ，ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）＝ｐ（ｄαｔ｜γ０：ｔ＋Ｌ，ｙ１：ｔ＋Ｌ）×ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）
（３１）
であり、ここで、ｐ（ｄαｔ｜γ０：ｔ＋Ｌ，ｙ１：ｔ＋Ｌ）はガウス分布であり、このガウス分布のパラメータは、カルマンフィルタタイプの技術を用いて計算され得る。従って、ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）の近似と仮定して、ｐ（ｄｘｔ，ｄθｔ｜ｙ１：ｔ＋Ｌ）の近似は、容易に取得され得る。周辺インポータンス分布および関連のインポータンス重みを
【０２３１】
【数７８】
0004402879.tif 000092
【０２３２】
と定義し、π（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）により分布されたサンプル
【０２３３】
【数７９】
0004402879.tif 000093
【０２３４】
のセットが利用可能であると想定すると、ＩＬ（ｆｔ）の代替的なベイズインポータンスサンプリング推定は、ｐ（αｔ｜γ０：ｔ＋Ｌ，ｙ１：ｔ＋Ｌ）ｆｔ（ｘｔ，θｔ）が閉じた式（ｃｌｏｓｅｄ  ｆｏｒｍ  ｅｘｐｒｅｓｓｉｏｎ）で評価され得るとすれば、結果として
【０２３５】
【数８０】
0004402879.tif 000094
【０２３６】
になる。（３３）において、周辺正規化インポータンス重みは、
【０２３７】
【数８１】
0004402879.tif 000095
【０２３８】
により与えられる。
【０２３９】
直観的に、所与の精度に到達するために、
【０２４０】
【数８２】
0004402879.tif 000096
【０２４１】
は、
【０２４２】
【数８３】
0004402879.tif 000097
【０２４３】
にわたるサンプルの数を低減する必要がある。周辺分布
【０２４４】
【数８４】
0004402879.tif 000098
【０２４５】
からのサンプルのみを必要とするためである。Ａ．Ｄｏｕｃｅｔ，Ｊ．Ｆ．Ｇ．ｄｅ  ＦｒｅｉｔａｓおよびＮ．Ｊ．Ｇｏｒｄｏｎ（ｅｄｓ．）らによる「Ｓｅｑｕｅｎｔｉａｌ  Ｍｏｎｔｅ  Ｃａｒｌｏ  Ｍｅｔｈｏｄｓ  ｉｎ  Ｐｒａｃｔｉｃｅ」Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇ，２０００年６月、に開示されるように、推定の分散は低減されることが証明され得る。この開示の内容は、参考のため、本明細書中に援用される。この場合、これは重要である。なぜなら、パラメータの数は常時（すべての混合フィルタおよびＡＲプロセスが同じ長さを有すると想定する場合）、混合フィルタのパラメータ（複数）がｍ２Ｌ−ｍＬであり、ただし、Ｌは大きい可能性があり、
観測騒音のパラメータ（単数または複数）がｍまたは１であり、
自己回帰プロセスのパラメータ（複数）がｎｌ＋ｎであり、
ソースのパラメータ（複数）がｎである
からである。
【０２４６】
どの積分が分散を最も良く低減することを可能にするかは明確でないが、少なくともパラメータ空間における検索に関して、混合フィルタおよび自己回帰フィルタの積分が好適であるように思われる。
【０２４７】
これらの結果を仮定して、続く議論は、式
【０２４８】
【数８５】
0004402879.tif 000099
【０２４９】
からのインポータンス分布を用いて、
【０２５０】
【数８６】
0004402879.tif 000100
【０２５１】
およびＩＬ（ｆｔ）の近似を取得するベイズインポータンスサンプリング法に焦点を当てる。これまで説明された方法は、バッチ法である。逐次法がどのように取得され得るかは、以下に説明される。
【０２５２】
離散的な時間ｔにおけるインポータンス分布は因数分解されて
【０２５３】
【数８７】
0004402879.tif 000101
【０２５４】
になり得る。
【０２５５】
この目的は、任意の時間ｔにおいて、分布ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）の推定を取得し、この推定を、次に、過去にシミュレーションされた軌道
【０２５６】
【数８８】
0004402879.tif 000102
【０２５７】
を修正することなく、時間通りに伝搬できるようにすることである。これは、π（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）は、π（ｄγ０：ｔ−１＋Ｌ｜ｙ１：ｔ−１＋Ｌ）を周辺分散として認めることを意味する。これは、インポータンス分散が一般的な式
【０２５８】
【数８９】
0004402879.tif 000103
【０２５９】
のインポータンス分散であることに限定されるならば可能である。
【０２６０】
このようなインポータンス分散は、インポータンス重みの帰納的評価、すなわちｗ（γ０：ｔ＋Ｌ）＝ｗ（γ０：ｔ−１＋Ｌ）ｗｔ＋Ｌを可能にし、この特殊な場合は、
【０２６１】
【数９０】
0004402879.tif 000104
【０２６２】
である。
【０２６３】
量ｐ（ｄｘｔ＋Ｌ｜ｘｔ＋Ｌ，βｔ＋Ｌ）は、式（１９）およびｐ（ｙｔ＋Ｌ｜ｘ０：ｔ＋Ｌ，βｔ＋Ｌ）により与えられるシステムの、ワンステップアヘッドカルマンフィルタ（ｏｎｅ  ｓｔｅｐ  ａｈｅａｄ  Ｋａｌｍａｎ  ｆｉｌｔｅｒ）を用いて、正規化定数まで計算され得、ｐ（ｙｔ＋Ｌ｜ｘ０：ｔ＋Ｌ，βｔ＋Ｌ）は、式（１７）により与えられるシステムの、カルマンフィルタの１つ先の段階を用いて計算され得る。
【０２６４】
インポータンス分布π（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）に関する選択の数は無限にあり、唯一の制限は、その支持がｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）の支持を含むことである。２つの可能性は、次に考察される。可能な策は、時間ｔ＋Ｌにおいて、γ０：ｔ−１およびｙ１：ｔ＋Ｌと仮定して、インポータンス重みの分散を最小化するインポータンス分布を選択することである。この条件を満たすインポータンス分布は、ｐ（ｄγｔ＋Ｌ｜γ０：ｔ−１＋Ｌ，ｙ１：ｔ＋Ｌ）であり、ただし、関連の増分インポータンス重みは、
【０２６５】
【数９１】
0004402879.tif 000105
【０２６６】
により与えられる。
【０２６７】
最適インポータンス分布からの直接的なサンプリングは困難であり、インポータンス重みを評価することは、解析上取り扱いが困難である。従って、目的は、概して、取り扱い易い近似、典型的には局所的線形性ｐ（ｄγｔ｜γ０：ｔ−１，ｙ１：ｔ）を用いて最適分布を模倣することである。その代わりに、混合部分最適法（ｍｉｘｅｄ  ｓｕｂｏｐｔｉｍａｌ  ｍｅｔｈｏｄ）が説明される。αと１−αの比率を有する２つのインポータンス分布π１およびπ２により、時間ｔにおいて粒子をサンプリングし、従って、インポータンス重み
【０２６８】
【数９２】
0004402879.tif 000106
【０２６９】
が次に、式
【０２７０】
【数９３】
0004402879.tif 000107
【０２７１】
と推定される（パラメータαである、ベルヌーイ分布によりＮ１とＮ２を無秩序に導き出すことが可能であるが、これは、推定量分散を増加させることに留意されたい）ことが提唱される。
【０２７２】
インポータンス分散π１（ｄｘｔ＋Ｌ｜ｘ１：ｔＬ−１，β１：ｔ＋Ｌ，ｙ１：ｔ＋Ｌ）は、分散
【０２７３】
【数９４】
0004402879.tif 000108
【０２７４】
がゼロの周りに集中する正規分散と解され、π２（ｄｘ１＋Ｌ｜ｘ１：ｔ＋Ｌ−１，β１：ｔ＋Ｌ，ｙ１：ｔ＋Ｌ）は、
【０２７５】
【数９５】
0004402879.tif 000109
【０２７６】
と解される。これは、ａ（ｉ）およびｈ（ｉ）の値、ならびに初期分布Ｐａｔ＋Ｌ｜ｔ＋ＬおよびＰｈｔ＋Ｌ｜ｔ＋Ｌとしての
【０２７７】
【数９６】
0004402879.tif 000110
【０２７８】
を用いて（１４）において記載される状態空間モデルのワンステップアヘッドカルマンフィルタから取得されるガウス分布である。この分布は、これらの前の分布からサンプリングされ、数式（３７）は
【０２７９】
【数９７】
0004402879.tif 000111
【０２８０】
を計算するために用いられる。他のインポータンス分布が可能であるが、このアプローチは、良好な結果をもたらし、サンプルのダイバーシティを保持するように思われることに留意されたい。
【０２８１】
式（３６）により具体化された形式のインポータンス分布に関して、インポータンス重みの分散は、Ｄｏｕｃｅｔ、ＧｏｄｓｉｌｌおよびＡｎｄｒｉｅｕならびにその中の参考資料により開示されるように、時間の経過とともに増加する（確率的に）にすぎない。従って、縮重現象を回避することが可能である。実際、アルゴリズムがいくらか反復された後、正規化インポータンス重みの１つを除くすべてはゼロに非常に近く、最終推定への貢献がほぼゼロである軌道を更新することに多大な計算労力が充てられる。この理由で、選択およびダイバーシティを含むことは、非常に重要である。これは、本明細書中でより詳細に議論される。
【０２８２】
選択（またはリサンプリング）手順の目的は、低正規化インポータンス重みを有する粒子を廃棄し、これらを高正規化インポータンス重みに乗算し、アルゴリズムの退化を回避することである。選択手順は、各粒子、例えば
【０２８３】
【数９８】
0004402879.tif 000112
【０２８４】
子の数Ｎｉ∈Ｎと関連し、従って
【０２８５】
【数９９】
0004402879.tif 000113
【０２８６】
であり、Ｎ個の新しい粒子
【０２８７】
【数１００】
0004402879.tif 000114
【０２８８】
を取得する。Ｎｉ＝０の場合、
【０２８９】
【数１０１】
0004402879.tif 000115
【０２９０】
は廃棄され、Ｎｉ＝０でない場合、これは時間ｔ＋１においてＮｉの子を有する。選択工程の後、すべての粒子の正規化インポータンス重みはＮ−１にリセットされ、従って、過去のインポータンス重みに関するすべての情報を廃棄する。従って、次の時間工程における選択の前の正規化インポータンス重みは、（３７）に比例する。これらは、
【０２９１】
【数１０２】
0004402879.tif 000116
【０２９２】
で表される。なぜなら、これらは、正規化インポータンス重みの任意の過去の値には依存しないからである。選択手順が各時間工程で実行される場合、選択工程の前の近似化分布が
【０２９３】
【数１０３】
0004402879.tif 000117
【０２９４】
により与えられ、選択工程の後の近似分布は
【０２９５】
【数１０４】
0004402879.tif 000118
【０２９６】
という結果になる。Ｄｏｕｃｅｔ、ｄｅ  ＦｒｅｉｔｕｓおよびＧｏｒｄｏｎらにより開示されるシステム的サンプリングは、その良好な分散特性ゆえに選択される。
【０２９７】
しかしながら、選択は、別の問題を引き起こす。リサンプリングステージにおいて、高インポータンス重みを有する任意の特定の粒子が数回複製される。特に、Ｌ＞０のとき、軌道は時間ｔ＋１〜ｔ＋ＬにＬ回リサンプリングされ、非常にわずかな明確な軌道が時間ｔ＋Ｌにおいて残存する。これは、サンプルの減損（ｄｅｐｌｅｔｉｏｎ）という古典的な問題である。結果として、粒子のクラウドは最終的に個々の粒子に崩壊する。この退化は、目標の分布の不満足な近似に至る。いくつかの部分最適法は、この問題を克服し、粒子間でダイバーシティを導入するために提示される。これらのほとんどは、カーネル密度法（ｋｅｒｎｅｌ  ｄｅｎｓｉｔｙ  ｍｅｔｈｏｄ）に基づく（Ｄｏｕｃｅｔ、ＧｏｄｓｉｌｌおよびＡｎｄｒｉｅｕらにより、およびＮ．Ｊ．Ｇｏｒｄｏｎ、Ｄ．Ｊ．ＳａｌｍｏｎｄおよびＡ．Ｆ．Ｍ．Ｓｍｉｔｈらによる「Ｎｏｖｅｌ  ａｐｐｒｏａｃｈ  ｔｏ  ｎｏｎｌｉｎｅａｒ／ｎｏｎ−Ｇａｕｓｓｉａｎ  Ｂａｙｅｓｉａｎ  ｓｔａｔｅ  ｅｓｔｉｍａｔｉｏｎ」ＩＥＥ  Ｐｒｏｃｅｅｄｉｎｇｓ−Ｆ，ｖｏｌ．１４０、ｎｏ．２、１０７〜１１３ページ、１９９３年により開示される）（これらの内容は、参考のため、本明細書中に援用される）。このカーネル密度法は、粒子の現在のセットに基づいて、カーネル密度推定を用いて確率分布を近似化し、そこから新しいセットの別個の粒子をサンプリングする。しかしながら、特定のカーネルの選択および構成は、必ずしも容易であるとは限らない。さらに、これらの方法は、さらなるモンテカルロバージョンを導入する。以下において、ＭＣＭＣ法が逐次インポータンスサンプリングと組み合わされ得、モンテカルロバージョンを増加させることなく、ダイバーシティをサンプルの中に導入する方法が示される。
【０２９８】
サンプルの減損を制限する効率的な方法は、単に、ＭＣＭＣ工程をシミュレーションベースのフィルタ／固定遅延平滑子に（Ｄｏｕｃｅｔ、ＧｏｄｓｉｌｌおよびＡｎｄｒｉｅｕらにより、およびＣ．Ｐ．ＲｏｂｅｒｔおよびＧ．Ｃａｓｅｌｌａらにより引用された、ＢｅｒｚｕｉｎｉおよびＧｉｌｋｓらによる「Ｍｏｎｔｅ  Ｃａｒｌｏ  Ｓｔａｔｉｓｔｉｃａｌ  Ｍｅｔｈｏｄｓ」、Ｓｐｒｉｎｇｅｒ  Ｖｅｒｌａｇ、１９９９年、を参照）付加することからなる。これは、ダイバーシティをサンプルの中に導入し、従って、サンプルの減損の問題を極端に低減する。時間ｔ＋Ｌに、粒子
【０２９９】
【数１０５】
0004402879.tif 000119
【０３００】
は、ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）により、周辺に分布されることを想定されたい。不変量分布ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）を有する遷移カーネル（ｔｒａｎｓｉｔｉｏｎ  ｋｅｒｎｅｌ）Ｋ（γ０：ｔ＋Ｌ｜ｄγ’０：ｔ＋Ｌ）が粒子の各々に適用された場合、新しい粒子
【０３０１】
【数１０６】
0004402879.tif 000120
【０３０２】
は、依然として、目標の分布に従って分布される。Ｍｅｔｒｏｐｏｌｉｓ−Ｈａｓｔｉｎｇｓ（ＭＨ）アルゴリズムまたはＧｉｂｂｓサンプラといった、標準的ＭＣＭＣ法のいずれかが用いられ得る。しかしながら、従来のＭＣＭＣ法とは逆に、遷移カーネルは、エルゴード的である必要がない。この方法は、さらなるモンテカルロバージョンを導入しないだけでなく、この方法は、目標の分布に関して、粒子の現在の分布の全変動標準（ｔｏｔａｌ  ｖａｒｉａｔｉｏｎ  ｎｏｒｍ）のみを低減し得るように推定値を改善する。
【０３０３】
時間ｔ＋Ｌ−１においてＮ∈Ｎ＊粒子
【０３０４】
【数１０７】
0004402879.tif 000121
【０３０５】
がｐ（ｄγ０：ｔ＋Ｌ−１｜ｙ１：ｔ＋Ｌ−１）により近似的に分布されると仮定すれば、モンテカルロ固定遅延平滑子は、時間ｔ＋Ｌに以下のように進行する。
【０３０６】
（逐次インポータンスサンプリング工程）
ｉ＝１，．．．，Ｎ，の場合、
【０３０７】
【数１０８】
0004402879.tif 000122
【０３０８】
である。
【０３０９】
ｉ＝１，．．．，Ｎ，の場合、（３７）および（４０）を用いて正規化されたインポータンス重み
【０３１０】
【数１０９】
0004402879.tif 000123
【０３１１】
を計算する。
【０３１２】
（選択工程）
高／低正規化されたインポータンス重みに関して粒子
【０３１３】
【数１１０】
0004402879.tif 000124
【０３１４】
を乗算廃棄し、例えば、システム的なサンプリングを用いてＮ個の粒子
【０３１５】
【数１１１】
0004402879.tif 000125
【０３１６】
を取得する。
【０３１７】
（ＭＣＭＣ工程）
ｉ＝１，．．．，Ｎ，の場合、
【０３１８】
【数１１２】
0004402879.tif 000126
【０３１９】
に不変量分布（ｉｎｖａｒｉｅｎｔ  ｄｉｓｔｒｉｂｕｔｉｏｎ）ｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）を有するマルコフ遷移カーネル
【０３２０】
【数１１３】
0004402879.tif 000127
【０３２１】
を適用し、Ｎ個の粒子
【０３２２】
【数１１４】
0004402879.tif 000128
【０３２３】
を取得する。
【０３２４】
ＭＣＭＣ遷移カーネルについては無限数の選択がある。ここで、ｏｎｅ−ａｔ−ａ−ｔｉｍｅ  ＭＨアルゴリズムが採用される。これは、時間ｔ＋Ｌに、マルコフ過程の値を時間ｔからｔ＋Ｌに更新する。より具体的には、
【０３２５】
【数１１５】
0004402879.tif 000129
【０３２６】
は、
【０３２７】
【数１１６】
0004402879.tif 000130
【０３２８】
により
【０３２９】
【数１１７】
0004402879.tif 000131
【０３３０】
を用いてサンプリングされる。
このアルゴリズムがｐ（ｄγ０：ｔ＋Ｌ｜ｙ１：ｔ＋Ｌ）を不変量分布として認めることを証明することは容易である。
【０３３１】
【数１１８】
0004402879.tif 000132
【０３３２】
からのサンプリングは、Ａ．ＤｏｕｃｅｔおよびＣ．Ａｎｄｒｉｅｕらによる「Ｉｔｅｒａｔｉｖｅ  ａｌｇｏｒｉｔｈｍｓ  ｆｏｒ  ｏｐｔｉｍａｌ  ｓｔａｔｅ  ｅｓｔｉｍａｔｉｏｎ  ｏｆ  ｊｕｍｐ  Ｍａｒｋｏｖ  ｌｉｎｅａｒ  ｓｙｓｔｅｍｓ」Ｐｒｏｃ．Ｃｏｎｆ．ＩＥＥＥ  ＩＣＡＳＳＰ、１９９９年（この開示の内容は、参考のため、本明細書中に援用される）に開示されるように、Ｏ（Ｌ＋１）複雑性の後方（ｂａｃｋｗａｒｄ−ｆｏｒｗａｒｄ）アルゴリズムを介して効率的に行われ得る。時間ｔ＋Ｌにおいて、このサンプリングは、ｉ番目の粒子に関して以下にまとめられるように、進行する。
【０３３３】
ｋ＝ｔ＋Ｌ，．．．，ｔ，の場合、２つのシステム（１７）および（１９）に関して（５０）〜（５１）に定義される情報フィルタを動作することにより
【０３３４】
【数１１９】
0004402879.tif 000133
【０３３５】
を計算および格納する。
（前方ステップ）
ｋ＝ｔ，．．．，ｔ＋Ｌの場合、
（４３）における提案（ｐｒｏｐｏｓａｌ）分布を用いて、プロポーサルγｋ〜ｑ（ｄγｋ｜γ−ｋ，ｙ０：ｔ＋Ｌ）をサンプリングする。
【０３３６】
現在の値
【０３３７】
【数１２０】
0004402879.tif 000134
【０３３８】
および提案された値γｋに関して、（４８）〜（４９）におけるカルマンフィルタの１工程を実行し、（４１）を用いて、これらの事後確率を計算し、かつ２つのシステム（１７）および（１９）に関して計算する。
【０３３９】
（４４）に定義される、ＭＨ受容確率（ａｃｃｅｐｔａｎｃｅ  ｐｒｏｂａｂｉｌｉｔｙ）
【０３４０】
【数１２１】
0004402879.tif 000135
【０３４１】
を計算する。
【０３４２】
【数１２２】
0004402879.tif 000136
【０３４３】
である場合、セット
【０３４４】
【数１２３】
0004402879.tif 000137
【０３４５】
であり、他の場合は
【０３４６】
【数１２４】
0004402879.tif 000138
【０３４７】
である。
【０３４８】
ＭＨ工程の各々に関する目標の事後分布は、
【０３４９】
【数１２５】
0004402879.tif 000139
【０３５０】
により与えられる。
【０３５１】
ｙｋおよびｘｋという２つの表現の間には類似点がある。これらの２つの項は、システム（１７）および（１９）に関する２つのカルマンフィルタを最初に用いて、同じ方法で計算され得、
【０３５２】
【数１２６】
0004402879.tif 000140
【０３５３】
を取得し、次に、２つの情報フィルタを用いて、
【０３５４】
【数１２７】
0004402879.tif 000141
【０３５５】
および
【０３５６】
【数１２８】
0004402879.tif 000142
【０３５７】
に関する類似の表現を取得する。含まれる異なった行列は、以下のように定義される。
【０３５８】
【数１２９】
0004402879.tif 000143
【０３５９】
とすれば、ここで、
【０３６０】
【数１３０】
0004402879.tif 000144
【０３６１】
が
【０３６２】
【数１３１】
0004402879.tif 000145
【０３６３】
のゼロでない単一値
【０３６４】
【数１３２】
0004402879.tif 000146
【０３６５】
を含む対角行列であり、
【０３６６】
【数１３３】
0004402879.tif 000147
【０３６７】
は、ゼロでない単一の値に対応する
【０３６８】
【数１３４】
0004402879.tif 000148
【０３６９】
の列を含む行列であり、ここで、
【０３７０】
【数１３５】
0004402879.tif 000149
【０３７１】
は
【０３７２】
【数１３６】
0004402879.tif 000150
【０３７３】
の単一値の分解である。行列
【０３７４】
【数１３７】
0004402879.tif 000151
【０３７５】
は、
【０３７６】
【数１３８】
0004402879.tif 000152
【０３７７】
により与えられる。
【０３７８】
ＭＨ工程を用いる（４１）における分布からサンプリングするために、ここで、プロポーサル分布は、
【０３７９】
【数１３９】
0004402879.tif 000153
【０３８０】
および
ｑ（γｋ｜γ−ｋ，ｙ０：ｔ＋Ｌ）αｐ（βｋ｜βｋ−１，βｋ＋１）ｑ（ｘｋ｜ｘ−ｋ，β０：ｔ＋Ｌ，ｙ０：ｔ＋Ｌ）                                            （４３）
であると解され、
これは、システム（１４）に関するワンステップアヘッドカルマンフィルタを必要とする。両方の場合、
【０３８１】
【数１４０】
0004402879.tif 000154
【０３８２】
である。マルコフ連鎖の状態の現在の、および予定された新しい値がγｋおよびγ’ｋによりそれぞれ与えられた場合、ＭＨ受容確率は、
α（γｋ，γ’ｋ）＝ｍｉｎ｛１，ｒ（γｋ，γ’ｋ）｝          （４４）
となり、
ただし、受容比率は
【０３８３】
【数１４１】
0004402879.tif 000155
【０３８４】
により与えられる。
【０３８５】
各反復において、モンテカルロ固定遅延平滑子の計算上の複雑さは、Ｏ（（Ｌ＋１）Ｎ）であり、時間ｔからｔ＋Ｌのすべての軌道の経路、すなわち
【０３８６】
【数１４２】
0004402879.tif 000156
【０３８７】
および十分な統計学（ｓｕｆｆｉｃｉｅｎｔ  ｓｔａｔｉｓｔｉｃｓ）
【０３８８】
【数１４３】
0004402879.tif 000157
【０３８９】
を記憶する必要がある。
【０３９０】
各反復における、このアルゴリズムの計算の複雑さは、明らかにＯ（Ｎ）である。一見して、すべての軌道の経路
【０３９１】
【数１４４】
0004402879.tif 000158
【０３９２】
を記憶することが必要に思われ得るので、格納の要求は、時間とともに線形に増加する。実際、最適インポータンス分布と事前インポータンス分布の両方に関して、π（γｔ｜γ０：ｔ−１，ｙ１：ｔ）およびγ０：ｔ−１関連のインポータンス重みは、低次元の十分な統計
【０３９３】
【数１４５】
0004402879.tif 000159
【０３９４】
を介してのみγ０：ｔ−１に依存し、これらの値のみが記憶される必要がある。従って、格納の要求は、さらにＯ（Ｎ）であり、時間の経過とともに増加しない。
【０３９５】
（付録：カルマンフィルタ帰納）
システムは、
ｘｔ＋１＝Ａｔ＋１ｘ１＋Ｂｔ＋１υｔ＋１
（４６）
ｙｔ＝ＣｔＸｔ＋Ｄｔｗｔ                                        （４７）
が考察される。
【０３９６】
ここで公知であると想定されるシーケンス
【０３９７】
【数１４６】
0004402879.tif 000160
【０３９８】
すなわちカルマンフィルタの式は以下のとおりである。
【０３９９】
【数１４７】
0004402879.tif 000161
【０４００】
ここで、
【０４０１】
【数１４８】
0004402879.tif 000162
【０４０２】
と推定される。
【０４０３】
後方情報フィルタは、時間ｔ＋Ｌからｔへと進行する。
【０４０４】
【数１４９】
0004402879.tif 000163
【図面の簡単な説明】
【図１】    図１は、信号ソースおよび通信チャンネルを示す図である。
【図２】    図２は、部屋の中の音響音ソースおよび本発明の実施形態を構成する方法を実行する装置を示す図である。
【図３】    図３は、本発明の実施形態を構成する方法を示す流れ図である。
【図１】
0004402879.tif 000179
【図２】
0004402879.tif 000180
【図３】
0004402879.tif 000181