(57)【特許請求の範囲】
【請求項１】
Ｎチャネル（Ｎは２以上の整数）の受話信号に対して、それぞれ再生信号のチャネル間の相互相関が絶えず変動するような付加信号を加算して上記再生信号を生成し、
上記再生信号を、Ｎ個の反響経路を模擬した疑似反響経路に印加して疑似エコー信号を生成し、
上記Ｎ個の反響経路から得られたエコー信号から疑似エコー信号を差し引いてエコー信号を消去して誤差信号を求め、
再生信号よりも、付加信号をより多く含む修正用基本ベクトルを、付加信号ベクトルと受話信号ベクトルから生成し、
その修正用基本ベクトルと、上記誤差信号とから修正ベクトルを求め、
その修正ベクトルを用いて疑似エコー経路のインパルス応答を逐次更新する
各ステップを含む反響消去方法。
【請求項２】
請求項１記載の方法において、
受話信号ベクトルをａ倍（０＜ａ＜１）して付加信号ベクトルに加算して上記修正用基本ベクトルとすることを特徴とする反響消去方法。
【請求項３】
請求項１記載の方法において、
上記修正用基本ベクトルの生成は、
再生信号ベクトルを、受話信号ベクトルの線形和のベクトルと受話信号ベクトルに直交するベクトルとに分解し、
再生信号ベクトルから、受話信号ベクトルの線形和ベクトルのｂ倍（０＜ｂ＜１）を差し引いて上記修正用基本ベクトルとすることを特徴とする反響消去方法。
【請求項４】
請求項１乃至３の何れかに記載の方法において、
上記修正ベクトルを求めるステップは、
各修正用基本ベクトルの係数を、上記誤差信号、Ｎチャネルの再生信号および修正用基本ベクトルから決定し、
その決定した係数を対応する修正用基本ベクトルに与えて修正用基本ベクトルの線形和を求めて上記修正ベクトルとすることを特徴とする反響消去方法。
【請求項５】
請求項１乃至２の何れかに記載の方法において、
上記疑似エコー信号生成ステップは上記再生信号を周波数領域に変換し、その周波数領域の再生信号に対し、周波数領域で上記疑似反響経路によるフィルタ処理を行い、その処理結果を時間領域に変換して上記疑似エコー信号を生成し、
上記修正ベクトルを求めるステップは上記誤差信号及び上記修正用基本ベクトルをそれぞれ周波数領域に変換して、周波数領域の上記修正ベクトルを求め、上記インパルス応答の逐次更新ステップは周波数領域で行うことを特徴とする反響消去方法。
【請求項６】
請求項５記載の方法において、
上記周波数領域での上記疑似反響経路の更新ステップは、
Ａ．周波数領域で対応スペクトルごとに再生信号と修正用基本信号のクロススペクトルの全チャネルの総和を求め、
Ｂ．周波数領域で各クロススペクトルの総和の逆数をそれぞれ修正ベクトルにかけて修正ベクトルを補正し、
Ｃ．その補正された修正ベクトルを用いて周波数領域で疑似反響経路のインパルス応答を更新する
ことを特徴とする反響消去方法。
【請求項７】
請求項６記載の方法において、
前回求めたスペクトルごとの再生信号と修正用信号のクロススペクトルの短時間平均の総和と、今回求めた対応スペクトルの再生信号と修正用基本信号のクロススペクトルの総和とを重み付け加算して今回のクロススペクトルの短時間平均の総和を求め、この短時間平均の総和を上記ステップＡで求める総和とすることを特徴とする反響消去方法。
【請求項８】
Ｎチャネル（Ｎは２以上の整数）の受話信号を入力し、それぞれ再生信号のチャネル間の相互相関が絶えず変動するような付加信号を加算した上記再生信号を生成する手段と、
上記Ｎ個の再生信号を入力し、Ｎ個の反響経路を模擬した疑似反響経路を備え、疑似エコー信号を生成出力する疑似エコー信号生成手段と、
上記Ｎ個の反響経路から得られたエコー信号から上記疑似エコー信号を差し引いてエコー信号を消去して誤差信号を求める消去手段と、
再生信号よりも、付加信号をより多く含む修正用基本ベクトルを、付加信号ベクトルと受話信号ベクトルとから生成する手段と、
上記修正用基本ベクトルと、上記誤差信号とを入力して修正ベクトルを求める手段と、
上記修正ベクトルを用いて、上記疑似エコー信号生成手段の疑似反響経路のインパルス応答を逐次更新する逐次更新手段と
を具備する反響消去装置。
【請求項９】
請求項８記載の装置において、
上記修正用基本ベクトルを求める手段は受話信号ベクトルをａ倍（０＜ａ＜１）して付加信号ベクトルと加算する手段であることを特徴とする反響消去装置。
【請求項１０】
請求項８記載の装置において、
上記修正用基本ベクトルを求める手段は
再生信号ベクトルを、受話信号ベクトルの線形和のベクトルと、受話信号ベクトルと直交するベクトルとに分解する手段と、
上記再生信号ベクトルから、受話信号ベクトルの線形和ベクトルのｂ倍（０＜ｂ＜１）を差し引いて上記修正用基本ベクトルを出力する手段とよりなることを特徴とする反響消去装置。
【請求項１１】
請求項８乃至１０のいずれかに記載の装置において、
上記修正ベクトル生成手段は、上記誤差信号と、Ｎチャネルの再生信号および修正用基本ベクトルから上記各修正用基本ベクトルの係数を求める手段と、
上記求めた係数を対応する修正用基本ベクトルに与えて、これらの線形和を求めて修正ベクトルとして出力する手段とを備えることを特徴とする反響消去装置。
【請求項１２】
請求項８乃至９の何れかに記載の装置において、
上記疑似エコー信号生成手段は、上記再生信号を周波数領域に変換する手段と、その変換された周波数領域の再生信号に対し、上記疑似反響経路によるフィルタ処理を周波数領域で行う手段と、そのフィルタ処理された結果を時間領域に変換して上記疑似エコー信号を出力する手段とを備え、
上記修正ベクトルを求める手段は、上記誤差信号を周波数領域に変換する手段と、上記修正用基本ベクトルを周波数領域に変換する手段と、これら周波数領域に変換された誤差信号と修正用基本ベクトルにより周波数領域で上記修正ベクトルを求める手段とを備え、
上記逐次更新手段は上記周波数領域で求められた上記修正ベクトルが入力され、上記疑似反響経路のインパルス応答の逐次更新を周波数領域で行う手段であることを特徴とする反響消去装置。
【請求項１３】
請求項１２記載の装置において、
上記逐次更新手段は対応スペクトルごとに、周波数領域の再生信号と周波数領域の修正用基本信号とが入力され、これらのクロススペクトルの全チャネルの総和を求めて出力する第１手段と、
上記スペクトルごとのクロススペクトルの総和と上記周波数領域の修正ベクトルが入力され、各クロススペクトルの総和の逆数を修正ベクトルにかけて修正ベクトルを補正する第２手段と、
上記周波数領域の補正された修正ベクトルが入力され、上記疑似反響経路のインパルス応答を周波数領域で更新する第３手段とを備えることを特徴とする反響消去装置。
【請求項１４】
請求項１３記載の装置において、
上記第１手段は、前回求めたスペクトルごとのクロススペクトルの短時間平均の総和と、今回求めた対応スペクトルの再生信号と修正用基本信号のクロススペクトルの全チャネルの総和とを重み付け加算して、今回のクロススペクトルの短時間平均の総和を求めて上記出力する総和とする手段であることを特徴とする反響消去装置。
【請求項１５】
請求項１乃至７の何れかに記載の反響消去方法をコンピュータにより実行する反響消去プログラム。
【請求項１６】
請求項１５記載の反響消去プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
この発明は、例えば多チャネル音響再生系を有する通信会議システムに適用され、ハウリングの原因及び聴覚上の障害となる音響エコーを消去する多チャネル反響消去方法、その装置、そのプログラム及びその記録媒体に関するものである。
【０００２】
【従来の技術】
近年のデジタルネットワークと音声、画像の高能率符号化技術の進展により、複数の人が容易に参加でき、より自然な通話環境を提供できる多チャネルの拡声通話方式が研究されはじめている。その実現のためには、複数のスピーカからマイクロホンへの音響的回り込みを消去する多チャネル音響エコー消去の技術的課題と解決策の検討が必要となる。
Ｎ（≧２）チャネルの再生系とＭ（≧１）チャネルの収音系とで構成される通信会議システムは、図１に示すような構成により音響エコーの消去を行う。即ち各受話端子１1 〜１N からの受話信号は各スピーカ２1 〜２N で音響信号として再生され、各Ｎ個の音響エコー経路１０1 〜１０N を経て各マイクロホン３m （ｍ＝１，…，Ｍ）に回り込む。受話側の全Ｎチャネルの受話端子１1 〜１N と、Ｍチャネル送話側の送話端子５1 〜５M それぞれとの間にＮチャネルエコーキャンセル部４1 〜４M を接続して音響エコーを消去する。
【０００３】
上記Ｎチャネルエコーキャンセル部４m は、各収音チャネル毎に再生側の全Ｎチャネルと収音側の１チャネルとの間のＮ入力１出力時系列信号を処理する構成をとる。このＮチャネルエコーキャンセル部４m （ｍ＝１，…，Ｍ）の構成を図２に示す。Ｎチャネルの各受話信号ｘ1 （ｋ）…ｘN （ｋ）は疑似エコー信号生成部４１に入力されて疑似エコー信号が生成され、減算器４２により疑似エコー信号とマイクロホン３m からの収音信号ｙ（ｋ）との差である残留信号（誤差信号）が取り出され、この残留信号がエコー経路推定部４３に帰還され、推定エコー経路が逐次修正される。疑似エコー信号生成部４１は一般にフィルタで構成され、そのフィルタの係数がエコー経路推定部４３により逐次修正される。これら疑似エコー信号生成部４１、およびエコー経路推定部４３は適応フィルタを構成しており、以後、これら全体を適応フィルタと記すこともある。
【０００４】
実際の通信会議では、多くの場合１人の話者音声が対地から多チャネルで送出されて多チャネル受話信号となる。この受話信号のチャネル間相互相関は非常に高いために、エコーが消去されている状態であっても、推定されたエコー伝達特性と真のエコー伝達特性は必ずしも一致しないことが知られており、文献M.M. Sondhi,D.R.Morgan,andJ.L.Hall,“Stereo-phonic Acoustic Echo Cancellation - An Overview of the Fundamental Problem,”IEEE Signal Processing Letters, vol.2, no.8,pp.148-151(1995)に詳細に解析されている。推定されたエコー伝達特性と真のエコー伝達特性が一致していないと、対地で話者が交代して受話信号のチャネル間相互相関が変化すると突然音響エコーが消去されなくなり、送話信号として対地に送出される現象が生じる。
【０００５】
このことを、図１の第ｍ収音チャネルに接続されているＮチャネルエコーキャンセル部４m について見てみる。Ｎチャネル入力信号をｘ1 （ｋ）…ｘN （ｋ）、収音された信号をｙ（ｋ）、第ｎチャネル（ｎ＝１，…，Ｎ）の再生器２n から収音器３m までの音響エコー経路１０n のインパルス応答をｈn （ｋ）、その長さをＬとする。Ｎチャネル入力信号と収音信号の間には次の関係がある。
ｙ(k)＝Σi=0L-1ｈ1(i)ｘ1（ｋ−ｉ）＋…＋Σi=0L-1ｈN(i)ｘN（ｋ−ｉ）
各チャネルのインパルス応答と入力信号を
ｈn＝［ｈn(０）…ｈn(Ｌ−１）］T
ｘn(ｋ） ＝［ｘn(ｋ）…ｘn(ｋ−Ｌ＋１）］T
のようにベクトル化し、さらに全Ｎチャネルのインパルス応答と入力信号を
ｈ ＝［ｈ1T…ｈNT］T
ｘ（ｋ)＝［ｘ1T(ｋ)…ｘNT(ｋ)］T
のように１つのベクトルにまとめると、Ｎチャネル入力信号と収音信号の関係は次のように記述される。
【０００６】
ｙ（ｋ）＝ｈT ｘ（ｋ）＝ｈ1Tｘ1(ｋ）＋…＋ｈNTｘN(ｋ）
第ｍ収音チャネルに接続されているＮチャネル・エコーキャンセル部４m は、図２に示すように構成されており、収音される信号ｙ（ｋ）をＮチャネル入力信号ｘnT(ｋ） から疑似エコー信号生成部４１により予測する。実際に収音された信号と予測された信号の差ｅおよび過去のＮチャネル入力信号に基づいて、収音信号と予測信号の差が小さくなるようにエコー経路推定部４３で、疑似エコー信号生成部４１と構成するフィルタの係数ｗ（ｋ）が逐次修正される。
【０００７】
過去のＮチャネル入力信号ベクトルをどこまで考慮するかにより、ＮＬＭＳ法、射影法、ＲＬＳ法などの適応アルゴリズムがある。射影法では、
【０００８】
【数１】
0003673727.tif 000002
【０００９】
のように修正ベクトルｄｗ（ｋ）が過去ｐ個の入力信号ベクトルの線形和である、という制約条件のもとで、過去ｐ個の入力信号の関係
ｙ（ｋ）＝ｗT（ｋ＋１)ｘ（ｋ）
：
ｙ（ｋ−ｐ＋１）＝ｗT（ｋ＋１)ｘ（ｋ−ｐ＋１）
を満たす適応フィルタ係数ｗ（ｋ＋１）＝ｗ（ｋ）＋ｄｗ（ｋ）を求める。この修正ベクトルｄｗ（ｋ）は
Ｘ（ｋ）＝［ｘ（ｋ）…ｘ（ｋ−ｐ＋１）]
ｅT（ｋ)＝［ｙ（ｋ）…ｙ（ｋ−ｐ＋１）］−ｗT（ｋ)Ｘ（ｋ）
ｃ＝（ＸT（ｋ)Ｘ（ｋ))-1ｅ（ｋ）
ｄｗ（ｋ）＝Ｘ（ｋ）ｃ
なる計算により得られる。Ｘ（ｋ）は入力信号ベクトルからなる入力信号行列であり、ｅ（ｋ）は収音信号と疑似エコー信号との誤差からなるベクトル、ｃは修正ベクトルを構成するための修正係数である。エコーキャンセル後の残留信号ｙ（ｋ）−ｗT（ｋ)ｘ（ｋ）を用いて、図２及び図３に示すように多チャネルエコーキャンセラ４m を構成できる。実際には推定を安定にするために０〜２の値をとるステップサイズμを用いて
ｗ（ｋ＋１）＝ｗ（ｋ）＋μＸ（ｋ）ｃ
により、適応フィルタの係数を更新する。
【００１０】
以上の適応信号処理が、図２中の音響エコー経路推定部４３で行われる。音響エコー経路推定部４３内では、図３に示すように、入力信号行列生成部４３１にて入力信号ｘ1 （ｋ）…ｘN （ｋ）から入力信号行列Ｘ（ｋ）が生成される。誤差ベクトル生成部４３４では、これまでの残留信号から誤差ベクトルを生成し、修正係数算出部４３２では誤差信号ベクトルｅと入力信号行列Ｘ（ｋ）から修正係数ｃを算出する。フィルタ係数更新部４３３では、修正係数ｃと入力信号行列Ｘ（ｋ）とから修正ベクトルｄｗ（ｋ）を求め、適応フィルタ係数ｗ（ｋ）を更新する。なお３次以上の射影アルゴリズムを用いる場合は誤差ベクトル生成部４３４にこれまでの入力信号と修正係数も入力する必要がある。
ところで適応フィルタ係数の修正法としてのＮＬＭＳ法（学習同程法）は射影法をｐ＝１とした時と一致する。実際に収音された信号ｙ（ｋ）と適応フィルタにより予測された信号との差ｅ（ｋ）は、
ｅ（ｋ）＝ｙ（ｋ）−Σn=1NｗnT（ｋ）ｘn(ｋ)
により計算される。この誤差をもちいて修正ベクトル
ｄｗn(ｋ)＝ｅ（ｋ）ｘn(ｋ)／Σn=1NｘnT（ｋ）ｘn（ｋ）（ｎ＝１，…,Ｎ）
を求め、各チャネルの適応フィルタを
ｗn(ｋ＋１)＝ｗn（ｋ）＋μｄｗn（ｋ）（ｎ＝１，…,Ｎ）
により更新する。ただしｗn(ｋ)は要素数Ｌのベクトルであり、第ｎチャネルの適応フィルタ係数のベクトルである。またμは推定を安定にするために設定されるステップサイズである。
ＮＬＭＳ法では、疑似エコーを生成するための畳み込み演算と適応フィルタの修正を逐次行うために、演算量が非常に大きくなる。文献E.R.Ferrara,“Fast Implementation of LMS adaptive filters,”IEEE Trans.Acoust,Speech,Signal Processing,vol.ASSP-28,pp.474-475(1980)で提案されている適応アルゴリズム（以下ＦＬＭＳ法と記す）は、適応フィルタの更新を逐次処理からＬサンプル毎のブロック処理に変更し、ＦＦＴをもちいてブロック信号処理を行うことで演算量を大幅に削減している。このアルゴリズムは、時刻ｋで適応フィルタが更新されるとき
ｄｗn(ｋ)＝Σi=0L-1ｅ（ｋ−ｉ）ｘn(ｋ−ｉ）（ｎ＝１，…,Ｎ）
のような畳込み演算により修正ベクトルｄｗn(ｋ)を計算する。この部分と疑似エコー生成部分の畳込み演算は、チャネル毎にＦＦＴ（高速離散フーリエ変換）をもちいて効率よく実行できるので、演算量が大幅に減少する。
このＦＬＭＳ法に、さらに文献D.Mansour and A.H.Gray,“Unconstrained Frequency-Domain Adaptive Filter,”IEEE Trans.on Acoust,Speech,Signal Processing,vol.ASSP-30,No.5,pp.726-734(1982)で提案されている白色化処理を組み合わせることによって、音声信号のようにスペクトルに偏りのある信号が入力されても、適応フィルタの収束特性は劣化しなくなる。
ここでは、多入力１出力適応フィルタに白色化処理付きのＦＬＭＳ法を適用する従来方法について説明する。このアルゴリズムでは、適応フィルタ長がＬのとき、Overlap-save方式をもちいてＬサンプル毎に長さ２Ｌの信号ベクトルをＦＦＴして処理することで、効率の高い畳込み演算処理を実現している。このアルゴリズムは、以下のステップからなる。
【００１１】
ステップ１
各チャネルの入力信号ｘn(ｋ)（ｎ＝１，…,Ｎ）を、Ｌサンプル毎に長さ２Ｌの入力信号ベクトルにブロック化してＦＦＴにより周波数領域に変換し、ベクトルの要素を対角成分に持つ行列Ｘnf（ｋ）を算出する。数式を用いると、
Ｘnf（ｋ）＝diag（ＦＦＴ（［ｘn(ｋ−２Ｌ＋１),…,ｘn(ｋ)］T））（ｎ＝１，…,Ｎ）
と記述される。ただし関数ＦＦＴ（ｘ）はベクトルｘをＦＦＴ変換する関数である。また関数diag（ｘ）によって、ベクトルｘはその要素を対角成分とする行列に変換される。すなわち  ｘ＝［ｘ（１）…ｘ（２Ｌ）］T のとき
【数２】
0003673727.tif 000003
である。
【００１２】
ステップ２
周波数領域でＸnf（ｋ）とｗnf（ｋ）を掛けることで、入力信号ベクトルをチャネルごとにフィルタ処理する。そして計算結果を逆ＦＦＴ（ＩＦＦＴ）処理し、時間領域での信号ベクトルｙ＾n（ｋ）（ｎ＝１，…，Ｎ）を得る。
ｙ＾n(ｋ）＝［ＩL ０L］ＩＦＦＴ（Ｘnf（ｋ）ｗnf（ｋ））
ただしｗnf（ｋ）（ｎ＝１，…，Ｎ）は要素数２Ｌの複素数ベクトルであり、
逆ＦＦＴ変換して前半Ｌ個を取り出すと、第ｎチャネル適応フィルタのインパルス応答になる。また０L はＬ×Ｌの零行列、ＩL はＬ×Ｌの単位行列である。
ステップ３
信号ベクトルｙ＾n(ｋ）（ｎ＝１，…，Ｎ）を加算して、疑似エコー信号のベクトルｙ＾（ｋ）を得る。
ｙ＾（ｋ）＝Σn=1Nｙ＾n（ｋ）
ステップ４
時間領域にて収音信号ベクトルｙ（ｋ）と疑似エコーの信号ベクトルｙ＾（ｋ）との差から誤差信号ベクトルを求め、ＦＦＴにより周波数領域に変換する。
ｅf(ｋ)＝ＦＦＴ（［０，…,０，ｙT(ｋ）−ｙ＾T(ｋ）］T )
ただし  ｙ（ｋ）＝［ｙ（ｋ−Ｌ＋１）…ｙ（ｋ）］T であり、ＦＦＴ［  ］内の０の数はＬ個であり、ｅf(ｋ)の要素数を２Ｌ個にするためである。
【００１３】
ステップ５
誤差信号と入力信号を周波数領域で処理し、修正ベクトルｄｗnf（ｋ）（ｎ＝１，…，Ｎ）を求める。
先ず以下のようにＸ*nf(ｋ)とｅf(ｋ)の積を逆ＦＦＴし、その結果の前半のＬ個を取出しｖnf(ｋ)を求める。
ｖnf(ｋ)＝［ＩL ０L］ＩＦＦＴ（Ｘ*nf（ｋ）ｅf(ｋ)）
ただし行列Ｘ*nf(ｋ)の各成分は行列Ｘnf（ｋ)各成分の複素共役である。
次にｖnf(ｋ)Tの後にＬ個の０を埋めてＦＦＴを行う。
ｄｗnf(ｋ)＝ＦＦＴ（［ｖnfT(ｋ),０，…，０］T ）
ステップ６
各チャネルの適応フィルタを次式で更新する。
ｗnf（ｋ＋Ｌ）＝ｗnf（ｋ）＋Ｐ（ｋ）ｄｗnf（ｋ）
行列Ｐ（ｋ）は、修正ベクトルｄｗnf（ｋ）を補正しており、
【数３】
0003673727.tif 000004
により計算される。μは０〜１の値をとるステップサイズである。関数Ｔ（Ｘnf(ｋ),ｉ）は行列Ｘnf(ｋ)の（ｉ，ｉ）要素を引き出す。行列Ｐ（ｋ）の対角要素の分母に含まれるｐ（ｋ，ｉ）は、周波数成分ごとに第１〜Ｎチャネルの入力信号パワーの短時間平均の総和を求めたものである。δは分母が０になることを防止するための微小な正定数である。βは前回の短時間平均パワーの総和ｐ（ｋ−Ｌ，ｉ）と今回の短時間パワーとの短時間平均をとるための平滑化定数であり、０〜１の値をとる。入力信号が音声のように有色性信号のとき、ｄｗnf（ｋ）に行列Ｐ（ｋ）をかけることは入力信号の白色化処理に対応し、有色信号が入力されたときの適応フィルタの収束速度を向上させることが知られている。
エコー経路の特性は、周波数領域でｗnf（ｋ）（ｎ＝１，…,Ｎ）として推定される。このベクトルを逆フーリエ変換することで、各エコー経路インパルス応答の推定値が得られる。
Ｎ入力１出力適応フィルタについてチャネル当りの適応フィルタ長をＬとするとき、Ｌサンプル分の疑似エコー信号を算出するのに必要となる積算の演算量は、ＮＬＭＳ法では、ＮＬ（２Ｌ＋４）である。一方、ＦＬＭＳ法で必要となる積算の演算量はＮＬ（１０ｌｏｇＬ＋８）である。チャネル当りの適応フィルタ長をＬ＝１０２４とするとき、ＦＬＭＳ法の演算量はＮＬＭＳ法の約５．３％になり、演算処理が非常に効率的になる。
【００１４】
図１中のＮチャネルエコーキャンセル部４mは、ＦＬＭＳ法では図４に示す構成で実現される。第ｎチャネルの入力信号ｘn(ｋ)（ｎ＝１，…，Ｎ）は、ＴＦ変換部４４ｎにてステップ１のようにブロック化され周波数領域に変換される。ステップ２のように入力信号がフィルタ係数により周波数領域でフィルタ処理部（疑似反響経路）４５ｎによりフィルタ処理され、その処理結果がＦＴ変換部４６ｎ（ｎ＝１，…,Ｎ）により時間領域に変換されて時間領域の信号ベクトルｙ＾n(ｋ)が得られる。信号ベクトル加算部４７では各信号ベクトルｙ＾n(ｋ)がステップ３のように加算されて時間領域での疑似エコーｙ＾(ｋ)が算出される。収音信号ｙ(ｋ)は、ブロック化部４８にてＬ個のサンプル（要素）にブロック化される。ＴＦ変換部４４ｎおよび収音信号のブロック化部４８は、各入力信号と収音信号の間に時間のズレが発生しないように信号をブロック化して、それぞれ信号ベクトルを生成する。
信号ベクトル減算部４９では、ステップ４のように収音信号ベクトルｙ(ｋ)から疑似エコーの信号ベクトルｙ＾(ｋ)が引かれ、誤差信号ベクトルｅ（ｋ）が求められ、これはＴＦ変換部５１にて周波数領域の誤差信号ベクトルｅf(ｋ)へ変換される。フィルタ係数更新部５２ｎ（ｎ＝１，…,Ｎ）では、ＴＦ変換部４４ｎからのＸnf（ｋ）とＦＴ変換部５１からのｅf(ｋ)を用いて、ステップ５及び６にしたがって周波数領域でフィルタ（疑似反響経路）を更新する。更新されたフィルタはフィルタ処理部４５ｎ（ｎ＝１，…，Ｎ）に反映される。なおステップ６での行列Ｐ（ｋ）の計算には全チャネル分のＸnf（ｋ）（ｎ＝１，…,Ｎ）を必要とするが、見やすくするために図４ではこの信号流れを省略している。
【００１５】
ところで入力信号のチャネル間の相互相関が一定で大きい場合には入出力信号の関係ｙ（ｋ）＝ｗT（ｋ)ｘ（ｋ）を満たすｗ（ｋ）が複数存在することが知られている。このため上記適応アルゴリズムにより推定されたインパルス応答が、対応する音響エコー経路のインパルス応答と一致するとは限らない。
このようなエコー伝達特性の誤推定を防ぐために、図５に示すように相関変動処理部６1 ，…，６Nを設けて、チャネル毎に受話信号を乱数で振幅変調して元の受話信号に付加して相互相関が絶えず変動している信号を生成し、各スピーカから再生すると同時に多チャネル・エコーキャンセラへの入力信号とする手法が特願平７−５０００２，文献S.Shimauchi and S.Makino,“Stereo Projection Echo Canceller with True Echo Path Estimation,”Proc.ICASSP95,vol.5,pp.3059-3062(1995)にて提案されている。その後、より効率的に相互相関が変動する信号を生成する手法として、文献J.Benesty,D.R.Morgan,and M.M.Sondhi,“A Better Understanding and an Improved Solution to the Problems of Stereophonic Acoustic Echo  Cancellation, ”Proc.ICASSP97,vol.1,pp.303-306(1997)では、受話信号を非線形関数で処理して元の受話信号に付加する方法が提案されている。
【００１６】
【発明が解決しようとする課題】
しかし受話信号に付加信号を加えてスピーカから再生したとき、元の受話信号と比較して聴感上違和感のない範囲におさめなければならないため、付加信号の信号パワーは制限され、受話信号のチャネル間の相互相関は依然高い。そのため真のエコー伝達特性を推定するにはＲＬＳ法のように計算量が大きくてノイズに敏感な適応アルゴリズムを用いる必要があると考えられており、ＮＬＭＳ法や射影法、ＦＬＭＳ法のように低演算量の適応アルゴリズムを用いた場合には、チャネル間相互相関の高い受話信号から修正ベクトルが生成されるための相互相関変動処理によるエコー経路インパルス応答推定性能の改善幅は小さい。
【００１７】
実際に数値シミュレーションを行った結果を図６に示す。この数値シミュレーションでは、サンプリング周波数を８ｋＨｚに設定し、音響エコー経路として残響時間２００ｍｓの部屋で実測した室内伝達関数を７００タップに打ち切って音響エコーを生成した。相互相関一定の２チャネル受話信号ｕ1 （ｋ），ｕ2 （ｋ）は、２本のマイクロホンで単一話者の音声を収音している状況を模擬することで生成した。適応フィルタのタップ数は１チャネル当り６００タップに設定し、適応アルゴリズムとして２次射影（ｐ＝２）をステップサイズμ＝０．５で適用した。
【００１８】
相関変動処理には、文献P.Eneroth,T.Gaensler,S.Gay and J.Benesty,“Studies of a wideband stereophonic acoustic echo canceller,”Proc.1999 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics,pp.207-210(1999)で用いられている半波整流方式
ｇ1（ｕ（ｋ））＝ｄ（ｕ(ｋ）＋｜ｕ(ｋ）｜)／２
ｇ2（ｕ（ｋ））＝ｄ（ｕ(ｋ）−｜ｕ(ｋ）｜)／２
を、聴感上違和感のほとんどないｄ＝０．２６で適用した。２チャネルエコーキャンセル部への入力は
ｘ1(ｋ）＝ｕ1(ｋ）＋ｇ1（ｕ1（ｋ））
ｘ2(ｋ）＝ｕ2(ｋ）＋ｇ2（ｕ2（ｋ））
になる。ｘ1(ｋ），ｘ2(ｋ）を以後再生信号と呼ぶ。またこれ以降、受話信号と付加信号をそれぞれ
ｕ（ｋ）＝［ｕ1(ｋ）…ｕ1(ｋ−Ｌ＋１）ｕ2(ｋ）…ｕ2(ｋ−Ｌ＋１）］
ｇ(ｋ）＝［ｇ1（ｕ1（ｋ））…ｇ1（ｕ1(ｋ−Ｌ＋１）ｇ2（ｕ2（ｋ）…ｇ2（ｕ2(ｋ−Ｌ＋１）］
とベクトル化して取り扱う。
【００１９】
相関変動処理を適用した場合（Ｂ）と適用しなかった場合（Ａ）の適応フィルタの推定性能を図６に示す。推定性能は、音響エコー経路のインパルス応答からなるベクトルｈと、適応フィルタの各インパルス応答の後ろに０詰めしてｈとサイズをそろえたベクトルｗ′(ｋ）との相対誤差
｜ｈ−ｗ′(ｋ）｜／｜ｈ｜
で評価した。図６のグラフによれば、相関変動処理を適用しない場合、はじめのｌｓ間は係数推定誤差がすばやく減少しているが、すぐに飽和し約−４．５ｄＢにとどまる。一方相関変動処理を用いた場合、係数推定誤差は飽和しないが減少は緩やかであり、１０ｓ後でも−７ｄＢ程度にとどまる。
【００２０】
この発明の第１の目的は従来よりも係数推定誤差を速く小さくすることができ、エコー消去性能を向上させた反響消去方法、その装置、そのプログラム及びその記録媒体を提供することにある。
この発明の第２の目的は第１の目的を達成しかつ演算量を大幅に減らすことができる反響消去方法、その装置、そのプログラム及びその記録媒体を提供することにある。
【００２１】
【課題を解決するための手段】
まずこの発明に至る考え方を説明する。
所で先の数値シミュレーションにおける条件で元の受話信号に由来する誤差ｅ0 （ｋ）と相互相関変動のための付加信号に由来する誤差ｅa(ｋ）、つまり
ｅ0(ｋ）＝（ｈ−ｗ′(ｋ))Tｕ（ｋ）
ｅa(ｋ）＝（ｈ−ｗ′(ｋ))Tｇ(ｋ）
の信号パワーをプロットすると図７のようになっている。点線がｅ0 (ｋ）、実線がｅa (ｋ）である。付加信号ｇn（ｕn (ｋ））（ｎ＝１，２）の信号パワーは元の受話信号ｕn（ｋ）（ｎ＝１，２）から約−１８ｄＢと小さなものであるが、このグラフによれば、付加信号に由来する誤差信号のパワーｅa (ｋ）は受話信号に由来する誤差信号のパワーｅ0 (ｋ）とほぼ同等である。すなわち、誤差ｙ（ｋ）−ｗT（ｋ)ｘ（ｋ）への付加信号ベクトルの寄与は受話信号ベクトルとほぼ同等である。
【００２２】
しかし、射影法をｐ＝１で適用したとき、すなわちＮＬＭＳ法では適応フィルタの係数は
ｗ(k+1)＝ｗ(k)＋μｅ(k)［（ｕ(k)＋ｇ(k)）／｜ｕ(k)＋ｇ(k)｜2 ］
のように更新されている。この更新式によれば、付加信号の修正ベクトルへの寄与は受話信号の約−１８ｄＢに過ぎず、付加信号ベクトルの情報は適応フィルタ係数の更新に対して過小評価されていることになる。
【００２３】
そこでこの発明では、付加信号と受話信号の修正ベクトルへの寄与が、各信号の誤差信号への寄与を反映するように、付加信号の比率が再生信号よりも多い修正用基本ベクトルｚ（ｋ）を受話信号と付加信号から生成する。そしてこのベクトルから適応フィルタの修正ベクトルを構成する。付加信号の比率が再生信号よりも多い修正用基本ベクトルｚ（ｋ）の一例としては、
ｚ（ｋ）＝ａｕ（ｋ）＋ｇ（ｋ），  ０＜ａ＜１
のようにすることが考えられる。
このように付加信号を強調した修正用基本ベクトルｚ（ｋ）を疑似反響経路のインパルス応答の修正ベクトルに反映させればチャネル間相互相関の小さい修正ベクトルを生成できる。つまり、受話信号ｕn(ｋ)に付加信号ｇn（ｕn（ｋ））が付加された信号を
ｘn(ｋ)＝ｕn(ｋ)＋ｇn（ｕn（ｋ））
とし、前記付加信号が強調された修正用信号を
ｚn(ｋ)＝ａｕn(ｋ)＋ｇn（ｕn（ｋ））
とし、これらを下記のようにベクトル化する。
ｘn(ｋ）＝［ｘn(k)…ｘn(ｋ−Ｌ＋１)］T （ｎ＝１，…，Ｎ）
ｚn(ｋ）＝［ｚn(k)…ｚn(ｋ−Ｌ＋１)］T （ｎ＝１，…，Ｎ）
この時、疑似反響経路により予測された信号と収音信号ｙ（ｋ）との誤差信号ｅ（ｋ）は次式で求められる。
ｅ（ｋ）＝ｙ（ｋ）−Σn=1NｗnT（ｋ）ｘn（ｋ）
この誤差信号を修正用基本ベクトルから修正ベクトルを次式により求められる。
ｄｗn(ｋ)＝ｅ（ｋ）ｚn(ｋ）    （ｎ＝１，…，Ｎ）
この修正ベクトルにより各チャネルの疑似反響経路のインパルス応答を次式により更新すればよい。
ｗn(k＋１)＝ｗn(k)＋μｄｗn(ｋ)（ｎ＝１，…，Ｎ）
ここでステップサイズμは毎回の繰り返しにおける補正の大きさを制御するパラメータである。
【００２４】
また再生信号より付加信号情報の比率が大きいベクトルｚ（ｋ）により、次のように修正ベクトルｄｗ（ｋ）を求めてもよい。つまりｄｗ（ｋ）がベクトルｚ（ｋ）…ｚ（ｋ−ｐ＋１）の線形和という制約条件のもとで、過去ｐ個の入出力信号の関係
ｙ（ｋ）＝ｗT（ｋ＋１)ｘ（ｋ）
：
ｙ（ｋ−ｐ＋１）＝ｗT（ｋ＋１)ｘ（ｋ−ｐ＋１）
を満たす修正ベクトルは、
Ｘ（ｋ）＝［ｘ（ｋ）…ｘ（ｋ−ｐ＋１）］
Ｚ（ｋ）＝［ｚ（ｋ）…ｚ（ｋ−ｐ＋１）］
ｅT（ｋ)＝［ｙ（ｋ）…ｙ（ｋ−ｐ＋１）］−ｗT（ｋ)Ｘ（ｋ）
ｃ＝（ＸT（ｋ)Ｚ（ｋ))-1ｅ（ｋ）
ｄｗ（ｋ）＝Ｚ（ｋ）ｃ
より求められる。実際にはステップサイズμを用いて
ｗ（ｋ＋１）＝ｗ（ｋ）＋μＺ（ｋ）ｃ
により、適応フィルタ係数を更新する。
【００２５】
つまりこの発明によれば
（Ａ）Ｎチャネルにおける受話信号に対して、それぞれ付加信号が付加された再生信号をそれぞれ生成し、
（Ｂ）この再生信号をＮ個のエコー経路を模擬した疑似反響経路に印加して疑似エコーを生成し、
（Ｃ）Ｎチャネルの再生信号がエコー経路を介して収音された音響エコーから疑似エコーを差し引くことで音響エコー消去を行い、
（Ｄ）音響エコーと疑似エコーの差、Ｎチャネルの受話信号および付加信号から修正ベクトルを求め、
（Ｅ）その修正ベクトルを用いて疑似反響経路のインパルス応答を逐次修正する
というステップにより多チャネル音響エコー消去を行い、
特にこの発明の１形態では上記ステップ（Ｄ）は
（Ｄ１）付加信号ベクトルと受話信号ベクトルから、再生信号よりも付加信号情報をより多く含む修正用基本ベクトルを生成し、
（Ｄ２）その修正用基本ベクトルの線形和を修正ベクトルとし、
（Ｄ３）その線形和に用いる各修正用基本ベクトルの係数を、音響エコーと疑似エコーの差、再生信号および修正用信号から決定する、
というステップを含むことがよい。さらにステップ（Ｄ１）において、受話信号ベクトルをａ倍（ａは０〜１の値）して付加信号ベクトルに加算して修正用基本ベクトルとする処理、もしくは
（Ｄ１−ａ）付加信号ベクトルと受話信号ベクトルから生成された再生信号ベクトルを、受話信号ベクトルの線形和とその受話信号ベクトルに直交するベクトルに分解して、再生信号ベクトルから受話信号ベクトルの線形和ベクトルのｂ倍（ｂは０〜１の値）を差し引いたベクトルを修正用基本ベクトルとする処理を行うとよい。
【００２６】
この発明の他の実施形態によれば前記付加信号を強調した信号ｚn(ｋ)を用いる考えをＦＬＭＳ法に導入する。この場合は
再生信号ｘn(ｋ）＝ｕn(ｋ）＋ｇn（ｕn(ｋ））を短時間区間ごとに周波数領域に変換し、周波数領域でＭ×Ｎ個の疑似反響経路によるフィルタ処理を行ない、時間領域に再変換してＭ個の疑似エコーを生成し（Ｎは受話チャネル数、Ｍは収音チャネル数）、
音響エコー信号と疑似エコー信号の誤差信号を短時間区間ごとに周波数領域に変換し、
修正用信号ｚn(ｋ)を短時間区間ごとに周波数領域に変換し、
周波数領域において変換された誤差信号と変換された修正用信号を処理して修正ベクトルを求め、
その修正ベクトルをもちいて周波数領域で疑似反響経路を更新する。
ここで短時間区間は疑似反響経路のタップ数Ｌと対応した時間又はこれより短かい時間である。
【００２７】
【発明の実施の形態】
実施例１
Ｎ（≧２）チャネルの再生系とＭ（≧１）チャネルの収音系とで構成される通信会議システムは、収音チャネル毎に図８のような再生側の全Ｎチャネルと収音側１チャネルとの間のＮ入力１出力時系列信号を処理するＮチャネルエコーキャンセル部７m を備える。
Ｎチャネルエコーキャンセル部７m には、受話信号と、その相関変動処理を経た受話信号が図９に示すように別々に入力され、これらから生成される再生信号
ｘn(ｋ）＝ｕn(ｋ）＋ｇn（ｕn(ｋ））（ｎ＝１,…,Ｎ）
が疑似エコー信号生成部（疑似反響経路）７１に入力されて疑似エコー信号が生成され、減算器７２により疑似エコー信号とマイクロホン３m からの収音信号ｙ（ｋ）との差である誤差信号ｅ（ｋ）が求められ、この誤差信号ｅ（ｋ）がエコー経路推定部７３に帰還される。
【００２８】
エコー経路推定部７３内は、図１０のようになっている。Ｚ（ｋ）生成部７３１では、受話信号ｕn(ｋ）と付加信号ｇn（ｕn(ｋ））から、各ｕn(ｋ）に対しａ（０＜ａ＜１）を乗算し、修正基本ベクトルとして
ｚ（ｋ）＝ａｕ（ｋ）＋ｇ(ｋ)
Ｚ（ｋ）＝［ｚ（ｋ）…ｚ（ｋ−ｐ＋１）］
のように付加信号情報に対する受話信号情報の比率が再生信号よりも小さい信号ベクトルｚ（ｋ）を生成し、更に修正用信号行列Ｚ（ｋ）を生成する。Ｘ（ｋ）生成部７３２では
ｘ（ｋ）＝ｕ（ｋ）＋ｇ(ｋ)
Ｘ（ｋ）＝［ｘ（ｋ）…ｘ（ｋ−ｐ＋１）］
のように受話信号ベクトルと付加信号ベクトルから再生信号行列Ｘ（ｋ）を生成する。ただし、ａはあらかじめ設定された０より大きく１より小さい値であり、実験により良い値を決めておく。
【００２９】
誤差ベクトル生成部７３５では、これまでの残留信号から誤差ベクトル
ｅT（ｋ）＝［ｙ（ｋ）…ｙ（ｋ−ｐ＋１）］−ｗT（ｋ）Ｘ（ｋ）
を生成し、修正係数算出部７３３ではＺ（ｋ），Ｘ（ｋ）と誤差ベクトルから修正用の係数からなるベクトル
ｃ＝（ＸT（ｋ）Ｚ（ｋ））-1ｅ（ｋ）
を算出する。フィルタ係数更新部７３４では、修正係数ｃとこれまでの修正用信号行列Ｚ（ｋ）から修正ベクトルＺ（ｋ）ｃを求め
ｗ（ｋ＋１）＝ｗ（ｋ）＋μＺ（ｋ）ｃ
により適応フィルタの係数を更新する。ただしμはステップサイズである。このときの計算量は、通常の射影アルゴリズムとほとんど変わらない。なお誤差ベクトル生成部７３５では、３次以上の射影アルゴリズムの場合は、これまでの再生信号及びこれまでの修正係数も用いて誤差ベクトルを生成する。
実施例２
Ｎ（≧２）チャネルの再生系とＭ（≧１）チャネルの収音系とで構成される通信会議システムは、収音チャネル毎に図８に示すように再生側の全Ｎチャネルと収音側１チャネルとの間のＮ入力１出力時系列信号を処理するＮチャネルエコーキャンセル部７m を備える。Ｎチャネルエコーキャンセル部の内部は図９、図１０に示したようになっている。
【００３０】
図９のＮチャネルエコーキャンセル部７mにおいて、疑似エコー信号生成部７１への入力信号のベクトルは、
ｘ（ｋ）＝ｕ（ｋ）＋ｇ(ｋ)
のように生成される。
この入力信号ベクトルは、受話信号成分および受話信号と直交する成分に分離できる。一例として受話信号成分として２時点の受話信号ベクトルｕ（ｋ），ｕ（ｋ−１）を考慮に入れた場合に、
【００３１】
【数４】
0003673727.tif 000005
【００３２】
のように再生信号ベクトルに含まれｕ（ｋ），ｕ（ｋ−１）のなす平面に直交するベクトルとしてｖ（ｋ）が求められる。Ｕ（ｋ）＝［ｕ（ｋ）  ｕ（ｋ−１）］とおき、上式に左からＵT（ｋ）をかけると
【００３３】
【数５】
0003673727.tif 000006
の関係から、ｓ0 ，ｓ1 は
【００３４】
【数６】
0003673727.tif 000007
【００３５】
により求まる。各ベクトルｕ（ｋ），ｕ（ｋ−１），ｖ（ｋ）の関係は図１１に示すようになる。つまり再生信号ベクトルｘ（ｋ）は受話信号ベクトルの線形和ｓ0 ｕ（ｋ）＋ｓ1 ｕ（ｋ−１）と、これに直交するベクトルｖ（ｋ）に分解できる。
このとき受話信号ベクトルの線形和からなる成分を１−ｂ倍することで、付加情報に対する受話信号情報の比率が小さい修正基本ベクトル
【００３６】
【数７】
0003673727.tif 000008
【００３７】
が生成できる。この式の右辺の第１項は受話信号ベクトルであり、第２項以下の項は受話信号ベクトル線形和のｂ倍の信号である。以上の式は２時点の受話信号ベクトルを用いる場合であるが、Ｕ（ｋ）をｒ時点の受話信号ベクトルから構成すれば、受話信号情報の付加信号情報に対する比率が再生信号よりも小さい修正基本ベクトルｚ（ｋ）は式（３）の右辺から求められる。ただし、ｂはあらかじめ設定された０〜１の範囲の値であり、実験により良い値を求めておく。
この処理がＺ（ｋ）生成部７３１にて行われたのち、
Ｘ（ｋ）＝［ｘ（ｋ）…ｘ（ｋ−ｐ＋１）］
Ｚ（ｋ）＝［ｚ（ｋ）…ｚ（ｋ−ｐ＋１）］
ｅT（ｋ）＝［ｙ（ｋ）…ｙ（ｋ−ｐ＋１）］−ｗT（ｋ）Ｘ（ｋ）
ｃ＝（ＸT（ｋ）Ｚ（ｋ））-1ｅ（ｋ）
ｗ（ｋ＋１）＝ｗ（ｋ）＋μＺ（ｋ）ｃ
により適応フィルタの係数を更新する。ただしμはステップサイズである。
【００３８】
実施例３
図８中のＮチャネルエコーキャンセル部７mの処理に、ＦＦＴを用いるブロック信号処理を適用する疑似反響経路のインパルス応答の更新処理の手順の例を以下に示す。
ステップ１
各チャネルの受話信号ｕn(ｋ)と相関変動処理のための付加信号ｇn(ｕn(ｋ))（ｎ＝１，…，Ｎ）から、再生信号ｘn(ｋ)と修正用信号ｚn(ｋ)を
ｘn(ｋ)＝ｕn(ｋ)＋ｇn(ｕn(ｋ))
ｚn(ｋ)＝ａｕn(ｋ)＋ｇn(ｕn(ｋ))    （ｎ＝１，…，Ｎ）
により生成する。ただしａは０より大きく１以下の値である。これら信号ｘn(ｋ),ｚn(ｋ)を、Ｌサンプル毎に長さ２Ｌの信号ベクトルにブロック化し、ＦＦＴをもちいて
Ｘnf（ｋ）＝diag（ＦＦＴ（［ｘn(ｋ−２Ｌ＋１），…，ｘn(ｋ)］T ））
Ｚnf（ｋ）＝diag（ＦＦＴ（［ｚn(ｋ−２Ｌ＋１），…，ｚn(ｋ)］T ））
（ｎ＝１，…,Ｎ）
のように周波数領域に変換する。
ステップ２
周波数領域でＸnf（ｋ）とｗnf（ｋ）を掛けることで、チャネルごとに入力信号ベクトルを、疑似反響経路でフィルタ処理する。このフィルタ処理結果を逆ＦＦＴ処理し、時間領域での信号ベクトルｙ＾n（ｋ）（ｎ＝１，…,Ｎ）を得る。
ｙ＾n(ｋ）＝［ＩL ０L］ＩＦＦＴ（Ｘnf（ｋ）ｗnf（ｋ））
ただし、０L はＬ×Ｌの零行列、ＩL はＬ×Ｌの単位行列である。
ステップ３
信号ベクトルｙ＾n(ｋ）（ｎ＝１，…，Ｎ）を加算して、疑似エコー信号のベクトルｙ＾（ｋ）を得る。
ｙ＾（ｋ）＝Σn=1Nｙ＾n（ｋ）
ステップ４
時間領域にて収音信号ベクトルｙ（ｋ）と疑似エコーのベクトルｙ＾（ｋ）から誤差信号ベクトルを求め、その誤差信号ベクトルをＦＦＴにより周波数領域に変換する。
ｅf(ｋ)＝ＦＦＴ（［０，…,０，ｙT(ｋ）−ｙ＾T(ｋ）］T )
ただし
ｙ（ｋ）＝［ｙ（ｋ−Ｌ＋１）…ｙ（ｋ）］T 
であり、ＦＦＴ［  ］中の０の数はＬ個である。
【００３９】
ステップ５
誤差信号ｅf(ｋ)と修正用信号ｚn(ｋ)を周波数領域で処理し、修正ベクトルｄｗnf（ｋ）を求める。
周波数領域でＺ*nf (ｋ)とｅf(ｋ)を乗算し、その結果を逆ＦＦＴして時間領域に変換し、その前半のＬ個を取出してｖ*nf (ｋ)とする。
ｖnf(ｋ)＝［ＩL ０L］ＩＦＦＴ（Ｚ*nf（ｋ）ｅf(ｋ)）
更にこのｖnf(ｋ)にＬ個の０を後詰めして、ＦＦＴにより周波数領域に変換する。
ｄｗnf(ｋ)＝ＦＦＴ（［ｖnfT(ｋ),０，…，０］T ）
ただし行列Ｚ*nf(ｋ)の各成分は修正用信号ｚn(ｋ)から生成された行列
Ｚnf(ｋ)各成分の複素共役である。
ステップ６
各チャネルの適応フィルタを次式で更新する。
ｗnf（ｋ＋Ｌ）＝ｗnf（ｋ）＋Ｐ（ｋ）ｄｗnf（ｋ）
ただし行列Ｐ（ｋ）は、
【数８】
0003673727.tif 000009
により計算される。μは０〜１の値をとるステップサイズである。関数Ｔ（Ｘnf(ｋ),ｉ）は行列Ｘnf(ｋ)の（ｉ，ｉ）番目の要素を引き出している。δは分母が０になることを防止するための微小な正定数である。行列Ｐ（ｋ）中のｐ（ｋ，ｉ）は、入力信号スペクトルＸnf（ｋ）と修正信号スペクトルＺnf（ｋ）のクロススペクトル短時間平均になっている。つまり前回の修正用信号と再生信号のクロススペクトルの短時間平均の全チャネル分の総和と、今回の修正用信号と再生信号のクロススペクトルの短時間の全チャネル分の総和とをβで重み付け加算して、今回の短時間平均総和を求める。
【００４０】
Ｎチャネルエコーキャンセル部７mの機能構成は、図１２に示すようになる。受話信号および付加信号をＴＦ変換する８１ｎは、図４中のＴＦ変換部４４ｎに対応している。受話信号ｕn(ｋ)には加算器８１１ｎにより付加信号ｇn（ｕn(ｋ)）が加算されて再生信号ｘn(ｋ)が生成され、再生信号ｘn(ｋ)はＴＦ変換部８１２ｎによってＸnf（ｋ）に変換される。また受話信号ｕn(ｋ)は減衰器８１３ｎによりａ倍（ただしａは０から１の値）され、加算器８１４ｎにより付加信号ｇn（ｕn(ｋ)）が加算されて修正用信号ｚn(ｋ)が生成される。修正用信号ｚn(ｋ)はＴＦ変換部８１５ｎによりＺnf（ｋ）に変換される。Ｘnf（ｋ）はフィルタ処理部（疑似エコー信号生成部）８２ｎに、Ｚnf（ｋ）はフィルタ更新部８８ｎにそれぞれ渡される。フィルタ処理部８２ｎ、ＦＴ変換部８３ｎ、ベクトル加算部８４では、ステップ２,３の処理を行い疑似エコー信号が生成される。マイクロホン３mからの収音信号ｙ（ｋ）は、ブロック化部８５でＬサンプルごとにブロック化され、ステップ４にしたがってベクトル減算部８６にて疑似エコー信号ベクトルとの誤差がとられ、その誤差ベクトルはＴＦ変換部８７で周波数領域へ変換される。フィルタ更新部８８ｎ（ｎ＝１，…，Ｎ）は、ステップ５，６にしたがって周波数領域でｗnf（ｋ）を更新することで、適応フィルタを更新する。
フィルタ更新部８８ｎは、図１３Ａに示すように、誤差信号ｅf(ｋ)と修正用基本ベクトルＺnf（ｋ）が修正ベクトル生成部８８１ｎに入力されて周波数領域で処理されて修正ベクトルｄｗnf（ｋ）が出力され、この修正ベクトルｄｗnf（ｋ）により逐次更新部８８２ｎにおいて周波数領域で疑似反響経路のインパルス応答がｗnf（ｋ）からｗnf（ｋ＋Ｌ）に更新される。
【００４１】
この際に、入力信号の白色化処理を行う場合は、ｄｗnf（ｋ）に対し、行列Ｐ（ｋ）により補正部８８３ｎで補正して、逐次更新部８８２ｎへ供給する。行列Ｐ（ｋ）の生成は、図１３Ｂに示すようにＸnf（ｋ）,Ｚnf（ｋ）の各ｉ番目の要素（スペクトル）（ｉ＝１，…，２Ｌ）ごとに乗算部８８４ｎで乗算してクロススペクトルを求め、これらクロススペクトルを加算部８８５で全チャネル分を加算し、この加算値ｐ′（ｋ，ｉ）と、前回の対応する（ｉ番目の）クロススペクトルの短時間平均値ｐ（ｋ−Ｌ，ｉ）とが平均化部８８６で荷重平均され、今回のｉ番目のクロススペクトル短時間平均ｐ（ｋ，ｉ）とする。この荷重平均は例えばβｐ（ｋ−Ｌ，ｉ）＋（１−β）ｐ′（ｋ，ｉ）＝ｐ（ｋ，ｉ）とする。βは０〜１の値であり、ｐ（ｋ，ｉ）の値が平滑化される。更にこれら各ｉ番目のクロススペクトル短時間平均の逆数にステップサイズμが乗算された各値を要素とする対角行列Ｐ（ｋ）が補正行列生成部８８９で生成される。
【００４２】
実施例４
ＦＬＭＳ法および実施例３の手法は、適応フィルタ（疑似反響経路）長がＬのとき、Ｌサンプル毎に過去２Ｌサンプル分の信号ブロックをもちいて、計算効率よく適応信号処理を行う手法である。この手法では、信号がＬサンプル分蓄積してから１フレーム分の適応信号処理が開始されるために、信号処理にＬサンプルの処理遅延が生じる。会議室用エコーキャンセラでは適応フィルタ長が部屋の残響時間と同等の例えば３００ｍｓ以上になるため、処理遅延が無視できない影響を持つ。またフィルタの更新頻度も低くなるために、例えばマイクロホンが動くなどしてエコー経路の特性が変動すると、エコーがすぐには消去されない問題が生じる。
文献J.S.Soo and K.K.Pang,“Multidelay Block Frequency Domain Adaptive Filter,”IEEE Trans.on ASSP,vol.ASSP 38,no.2,pp.373-376(1990)では、マルチディレイ・フィルタ（以下ＭＤＦと略す）をもちいて、処理遅延が大きく更新レートが低いというＦＬＭＳ法の問題を解決している。
周波数領域の信号処理では、畳み込み処理はオーバーラップセーブ法により実現されている。ＭＤＦは、この畳み込み処理がより小さいブロック同士のオーバーラップセーブ処理に分割できることを利用する。適応フィルタのタップ長をＬ、分割数をＤ（ただしＬはＤで割り切れる）、Ｌ′＝Ｌ／Ｄとすると、ＭＤＦ法ではＬ′サンプル毎に畳み込み処理が可能なため、Ｌ′サンプル毎に適応信号処理を適用することが可能になる。
実施例３の手法も、以下のステップのようにＭＤＦ法と組合わせることで、処理遅延と低更新レートの問題が解決される。
【００４３】
ステップ１
各チャネルの受話信号ｕn(ｋ)と相関変動処理のための付加信号ｇn(ｕn(ｋ))（ｎ＝１，…，Ｎ）から、再生信号ｘn(ｋ)と修正用信号ｚn(ｋ)を
ｘn(ｋ)＝ｕn(ｋ)＋ｇn(ｕn(ｋ))
ｚn(ｋ)＝ａｕn(ｋ)＋ｇn(ｕn(ｋ))    （ｎ＝１，…，Ｎ）
により生成する。ただしａは０より大きく１以下の値である。これらｘn(ｋ),ｚn(ｋ)をＬ′サンプル毎に長さ２Ｌ′の信号ベクトルにブロック化し、ＦＦＴをもちいて
Ｘnf（ｋ,Ｄ）＝diag（ＦＦＴ（［ｘn(ｋ−２Ｌ′＋１），…，ｘn(ｋ)］T ））
Ｚnf（ｋ,Ｄ）＝diag（ＦＦＴ（［ｚn(ｋ−２Ｌ′＋１），…，ｚn(ｋ)］T ））    （ｎ＝１，…,Ｎ）
のように周波数領域に変換する。また、疑似反響経路（適応フィルタ）長はＬであり、Ｄ−１個前まで計算結果を用いて、各Ｌ′についてフィルタ処理する必要があるから、
Ｘnf（ｋ,ｄ）＝Ｘnf（ｋ−Ｌ′,ｄ＋１）（ｄ＝１，…,Ｄ−１）
Ｚnf（ｋ,ｄ）＝Ｚnf（ｋ−Ｌ′,ｄ＋１）（ｄ＝１，…,Ｄ−１）
とする。
ステップ２
チャネルごとに周波数領域で掛け算処理を行うことで、入力信号ベクトルをフィルタ処理する。計算結果を逆ＦＦＴ処理し、時間領域での信号ベクトルｙ＾n(ｋ）を得る。
ｙ＾n(ｋ）＝［ＩL' ０L'］ＩＦＦＴ（Σd=1DＸnf（ｋ,ｄ）ｗnf（ｋ,ｄ））
ただし、０L'はＬ′×Ｌ′の零行列、ＩL'はＬ′×Ｌ′の単位行列である。
ステップ３
信号ベクトルｙ＾n(ｋ）（ｎ＝１，…，Ｎ）を加算して、疑似エコー信号のベクトルｙ＾（ｋ）を得る。
ｙ＾（ｋ）＝Σn=1Nｙ＾n（ｋ）
ステップ４
収音信号と疑似エコーの誤差信号のベクトルを
ｅf(ｋ)＝ＦＦＴ（［０，…,０，ｙT(ｋ）−ｙ＾T(ｋ）］T )
で算出する。ただし
ｙ（ｋ）＝［ｙ（ｋ−Ｌ′＋１）…ｙ（ｋ）］T 
であり、ＦＦＴ（［…］）内の０の数はＬ′個である。
【００４４】
ステップ５
誤差信号と修正用信号を周波数領域で処理し、修正ベクトルｄｗnf（ｋ）を求める。
ｖnf(ｋ,ｄ)＝［ＩL' ０L'］ＩＦＦＴ（Ｚ*nf（ｋ,ｄ）ｅf(ｋ)）
ｄｗnf(ｋ,ｄ)＝ＦＦＴ（［ｖnfT(ｋ,ｄ),０，…，０］T ）
（ｄ＝１，…，Ｄ）
ただし行列Ｚ*nf(ｋ)の各成分は行列Ｚnf(ｋ)各成分の複素共役であり、
ＦＦＴ（［  ］）内の０の数はＬ′個である。
ステップ６
各チャネルの適応フィルタを次式で更新する。
ｗnf（ｋ＋Ｌ′，ｄ）＝ｗnf（ｋ，ｄ）＋Ｐ（ｋ）ｄｗnf（ｋ,ｄ）
（ｄ＝１，…,Ｄ）
ただし、行列Ｐ（ｋ）は、
【数９】
0003673727.tif 000010
により計算され、μは０〜１の値をとるステップサイズである。またδは分母が０になることを防止するための微小な正定数である。
【００４５】
実施例４のＮチャネルエコーキャンセル部７m内部は、実施例３と同様に図１２に示したような機能構成をとる。受話信号ｕn(ｋ)には加算器８１１ｎにより付加信号ｇn（ｕn(ｋ)）が加算されて、再生信号ｘn(ｋ)が生成され、更にＴＦ変換部８１２ｎによってＸnf（ｋ）に変換される。また受話信号ｕn(ｋ)は減衰器８１３ｎによりａ倍（ただしａは０から１の値）され、加算器８１４ｎにより付加信号ｇn（ｕn(ｋ)）が加算されて、修正用信号ｚn(ｋ)が生成される。ｚn(ｋ)はＴＦ変換部８１５ｎによりＺnf（ｋ）に変換される。Ｘnf（ｋ）はフィルタ処理部８２ｎへ、Ｚnf（ｋ）はフィルタ更新部８８ｎに渡される。
フィルタ処理部８２ｎ、ＦＴ変換部８３ｎ、ベクトル加算部８４では、ステップ２，３の処理を経て疑似エコー信号が生成される。マイクロホン３mからの収音信号ｙ（ｋ）は、ブロック化部８５でブロック化され、ステップ４にしたがってベクトル加算部８６にて疑似エコー信号ベクトルとの誤差がとられ、ＴＦ変換部８７で周波数領域へ変換される。フィルタ更新部８８ｎではステップ５，６にしたがって適応フィルタが更新される。Ｎ入力１出力適応フィルタについて、チャネル当りの適応フィルタ長をＬとすると、Ｌサンプル分の疑似エコー信号を算出するのに必要となる積算の演算量は、ＮＬＭＳ法では、ＮＬ（２Ｌ＋４）である。一方、実施例４の方法で必要となる積算の演算量はＮＬ（（４Ｄ＋８）ｌｏｇ2 （Ｌ／Ｄ）＋１５Ｄ＋５）である。チャネル当りの適応フィルタタップ数をＬ＝１０２４とするとき、実施例４の方法で適応フィルタをＬ／４タップ毎に更新する場合の演算量はＮＬＭＳ法の約１２．５％であり、Ｌ／８タップ毎に更新する場合の演算量は約２０％である。このように演算量をＮＬＭＳ法と比較して低く抑えたまま、ＦＬＭＳ法と比較して、処理遅延を大幅に小さくすることができる。
【００４６】
以上述べたようにこの発明は再生信号ｘn(ｋ）と比較して付加信号ｇn(ｕn（ｋ））の比率が大きい修正用信号を用いて、疑似反響経路のインパルス応答を逐次更新するための修正ベクトルｄｗ（ｋ）を作る点に特徴がある。よってこの基本構成を図１４に示すと共に以下にその処理手順を説明する。
第１〜第Ｎチャネルの各受話信号を
ｕ1（ｋ）…ｕN（ｋ）
第１〜第Ｎチャネルの各付加信号を
ｇ1(ｕ1（ｋ））…ｇN(ｕN（ｋ））
第１〜第Ｎチャネルの適応フィルタ（疑似反響経路）のフィルタ係数（インパルス応答）を
ｗn ＝［ｗn(０)…ｗn(Ｌ−１)］T （ｎ＝１，…，Ｎ）
とする。ただし、Ｌは適応フィルタのチャネル当りのタップ数である。
第１〜第Ｎチャネルの受話信号に付加信号を付加して再生信号
ｘn(ｋ)＝ｕn(ｋ)＋ｇn(ｕn(ｋ)）（ｎ＝１，…，Ｎ）
とし、第１〜第Ｎチャネルの修正用信号を
ｚn(ｋ)＝ａｕn(ｋ)＋ｇn(ｕn（ｋ)）（ｎ＝１，…，Ｎ）
とし、それぞれｘn(ｋ)生成部９１、ｚn(ｋ)生成部９２で
ｘn(ｋ)＝［ｘn(k)…ｘn(ｋ−Ｌ＋１)］T （ｎ＝１，…，Ｎ）
ｚn(ｋ)＝［ｚn(k)…ｚn(ｋ−Ｌ＋１)］T （ｎ＝１，…，Ｎ）
のようにベクトル化する。
【００４７】
実際に収音された信号ｙ（ｋ）と適応フィルタ（疑似エコー信号生成部）９３により予測された信号ｙ＾（ｋ）との差ｅ（ｋ）を、減算部９４により
ｅ（ｋ）＝ｙ（ｋ）−Σn=1NｗnT（ｋ）ｘn(ｋ)
と求める。この誤差信号ｅ（ｋ）と修正用基本ベクトルｚn(ｎ＝１，…，Ｎ)とを用いて修正ベクトル生成部９５で
ｄｗn(ｋ)＝ｅ（ｋ)ｚn(ｋ)（ｎ＝１，…，Ｎ）
を求める。各チャネルの適応フィルタ９３の係数を逐次更新部９６により
ｗn （ｋ＋１）＝ｗn(k)＋μｄｗn(ｋ)（ｎ＝１，…，Ｎ）
と更新する。μは毎回の繰り返しにおける補正の大きさを制御するパラメータであり、ステップサイズと呼ばれる。なお修正用信号の生成はｚn（ｋ）＝ｕn（ｋ）＋ｂｇn（ｕn(ｋ））（ｎ＝1,・・・，Ｎ），ｂ＞１としてもよい。
【００４８】
効果の実証例（１）
再生チャネル数Ｎ＝２、収音チャネル数Ｍ＝１の音響系と多チャネル・エコーキャンセラに実施例１の手法を適用して数値シミュレーションを行った。サンプリング周波数を８ｋＨｚに設定し、音響エコー経路として残響時間２００ｍｓの部屋で実測した室内伝達関数を７００タップに打ち切って使用した。相互相関一定の２チャネル受話信号は、２本の４０ｄＢＳＮＲのマイクロホンで単一話者の音声を収音している状況をシミュレートして生成した。適応フィルタのタップ数は１チャネル当り６００タップに設定し、適応アルゴリズムに２次射影アルゴリズムを用いた。
【００４９】
相関変動処理として、半波整流方式
ｇ1（ｕ（ｋ））＝ｄ（ｕ(ｋ)＋｜ｕ(ｋ)｜）／２
ｇ2（ｕ（ｋ））＝ｄ（ｕ(ｋ)−｜ｕ(ｋ)｜）／２
を、ｄ＝０．２６として用いた。
推定性能は、音響エコー経路のインパルス応答からなるベクトルｈと適応フィルタの各インパルス応答後部に０詰めしてｈとサイズをそろえたベクトルｗ′(ｋ）との相対誤差
｜ｈ−ｗ′(ｋ）｜／｜ｈ｜
で評価した。
【００５０】
付加信号なしで従来の２次射影アルゴリズムをμ＝０．５で適用した場合（Ａ）、付加信号を加えて従来の２次射影アルゴリズムをμ＝０．５で適用した場合（Ｂ）、この発明の実施例１の手法をｐ＝２，ａ＝０．１，μ＝０．５で適用した場合（Ｃ）の適応フィルタの推定性能を図１５に示す。
このグラフによれば、従来の２次射影アルゴリズムでは、係数誤差は飽和しないものの減少は緩やかで、１０ｓ後の係数誤差は−７．０ｄＢ程度にとどまる。しかしこの発明法によれば、１０ｓ後の係数誤差は−１３．６ｄＢまで減少し、この発明が優れていることがわかる。
【００５１】
効果の実証例（２）
実際に数値シミュレーションを行った結果を図１６に示す。この数値シミュレーションでは、サンプリング周波数を８ｋＨｚに設定し、音響エコー経路として残響時間２００ｍｓの部屋で実測した室内伝達関数を７００タップに打ち切って音響エコーを生成した。相互相関一定の２チャネル受話信号ｕ1(ｋ),ｕ2(ｋ)は、２本のマイクロホンで単一話者の音声を収音している状況を模擬することで生成した。適応フィルタのタップ数は１チャネル当り５１２タップに設定し、従来適応アルゴリズムとしてＮＬＭＳ法とＦＬＭＳ法を適用した場合と、この発明の実施例４の方法とを比較した。
相関変動処理には、文献P.Eneroth,T.Gaensler,S.Gay and J.Benesty,“Studies of a wideband stereophonic acoustic echo canceler，”Proc.1999 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics,pp.207-210(1999)でもちいられている半波整流方式
ｇ1(ｕ（ｋ）)＝ｄ（ｕ（ｋ）＋｜ｕ（ｋ）｜）／２
ｇ2(ｕ（ｋ）)＝ｄ（ｕ（ｋ）−｜ｕ（ｋ）｜）／２
を、聴感上違和感のほとんどないｄ＝０．２６で適用した。２チャネルエコーキャンセル部への入力は
ｘ1(ｋ)＝ｕ1(ｋ)＋ｇ1(ｕ1(ｋ）)
ｘ2(ｋ)＝ｕ2(ｋ)＋ｇ2(ｕ2(ｋ）)
とした。
推定性能は、音響エコー経路のインパルス応答からなるベクトルｈと適応フィルタの各インパルス応答後部に０詰めしてｈとサイズをそろえたベクトルｗ′（ｋ）との相対誤差
｜ｈ−ｗ′（ｋ）｜／｜ｈ｜
で評価した。
付加信号を加えて従来のＮＬＭＳアルゴリズムをμ＝０．５で適用した場合（Ａ）、付加信号を加えてＦＬＭＳアルゴリズムをμ＝０．５で適用した場合（Ｂ）とこの発明の実施例４の手法を分割数Ｄ＝４，ａ＝０．１，μ＝０．５で適用した場合（Ｃ）の適応フィルタの推定性能を図１５に示す。
このグラフによれば、従来のＮＬＭＳアルゴリズムでは、係数誤差は飽和しないものの減少は緩やかで、１０ｓ後の係数誤差は−６．０ｄＢ程度にとどまり、ＦＬＭＳ法をもちいると白色化処理により係数誤差は約−１２ｄＢまで減少する。この発明実施例２の方法によれば、１０ｓ後の係数誤差はさらに低下し約−１８ｄＢにまで減少する。
【００５２】
上述においては受話信号に付加信号を付加して再生信号としたが、受話信号を処理して再生信号を得てもよい。この場合は、再生信号から受話信号を引算して付加信号を求めて、前述したこの発明の方法を行えばよい。また付加信号は受話信号を処理したものに限らず、受話信号とは独立に生成したものでもよい。
上述したこの発明による多チャネルエコー消去はコンピュータにより機能させることもできる。つまり例えば図１７に示すように、受話信号ｕ1(ｋ)，…，ｕN(ｋ)は入力部２１より入力され、音響エコー信号ｙ（ｋ）は入力部２２より入力され、これら入力信号はデータ記憶部２４に一時格納され、記憶部２４から読み出されて、付加信号の生成、再生信号行列Ｘ（ｋ）の生成、修正用基本ベクトル行列Ｚ（ｋ）の生成、疑似反響経路の生成、疑似エコー信号の生成、音響エコー信号から疑似エコー信号の除去、その誤差信号、修正用基本ベクトルから修正ベクトルの算出、修正ベクトルにより疑似反響経路のインパルス応答の逐次修正などを、ワーク用メモリ２５を必要に応じて用いて、プロセッサ２６によりメモリ２７に格納されているプログラムを実行させることにより行わせる。エコー消去された信号は出力部２３から出力される。この場合プロセッサを複数用いて、それぞれに処理を分担させると共に１つのプロセッサにより統括的処理を行うように、それぞれプロセッサに対応したプログラムを各別のメモリに格納してもよい。このプログラムはＣＤ−ＲＯＭ、磁気ディスクあるいは通信回線からインストールされて用いられる。
【００５３】
【発明の効果】
以上述べたようにこの発明によれば、付加信号情報に対する受話信号情報の比率を小さくした信号から適応フィルタの修正ベクトルを求める新しい適応アルゴリズムにより、多チャネル・エコー消去方法の推定性能を向上させている。特に適応フィルタ更新処理を周波数領域で行う場合は演算量を大幅に減少できる。これにより、対地で話者が交代し受話信号の相互相関が変化しても、エコーの増加を抑えることができる。
【図面の簡単な説明】
【図１】多チャネルエコー消去装置の一般的構成を示す図。
【図２】図１中のＮチャネルエコーキャンセル部４m の機能構成を示す図。
【図３】図２中のエコー経路推定部４３の機能構成を示す図。
【図４】周波数領域で適応フィルタ更新処理を行う従来機能を示す機能構成図。
【図５】受話信号に付加信号を加えた、多チャネルエコー消去装置の構成を示す図。
【図６】従来の方法による疑似エコー経路のインパルス応答係数推定誤差の時間経過を示す図。
【図７】受話信号による誤差信号パワー（点線）と、付加信号による誤差信号パワー（実線）の時間変化を示す図。
【図８】この発明が適用された多チャネルエコー消去装置の構成例を示す図。
【図９】図７中のこの発明によるＮチャネルエコーキャンセル部７m の機能構成例を示す図。
【図１０】図９中のエコー経路推定部７３の機能構成例を示す図。
【図１１】再生信号ベクトルを受話信号ベクトルの線形和と、これに直交するベクトルとに分解した様子を示す図。
【図１２】適用フィルタの更新処理を周波数領域で行うこの発明の実施例の機能構成を示す図。
【図１３】Ａは図１２中のフィルタ更新部を更に具体化した例を示す図、Ｂは図１２中のフィルタ更新部における白色化処理のための機能構成を示す図である。
【図１４】この発明の基本的な機能構成を示す図。
【図１５】従来法とこの発明方法（実施例１）による疑似エコー経路のインパルス応答係数推定誤差の時間経過を示す図。
【図１６】従来法とこの発明方法（実施例４）による疑似エコー経路のインパルス応答係数推定誤差の時間経過を示す図。
【図１７】この発明装置をコンピュータにより実行させる場合の構成例を示す図。
【図１】
0003673727.tif 000011
【図２】
0003673727.tif 000012
【図３】
0003673727.tif 000013
【図４】
0003673727.tif 000014
【図５】
0003673727.tif 000015
【図６】
0003673727.tif 000016
【図７】
0003673727.tif 000017
【図８】
0003673727.tif 000018
【図９】
0003673727.tif 000019
【図１０】
0003673727.tif 000020
【図１１】
0003673727.tif 000021
【図１２】
0003673727.tif 000022
【図１３】
0003673727.tif 000023
【図１４】
0003673727.tif 000024
【図１５】
0003673727.tif 000025
【図１６】
0003673727.tif 000026
【図１７】
0003673727.tif 000027