(57)【特許請求の範囲】
【請求項１】
Ｎチャネル（Ｎは２以上の整数）の受話信号をそれぞれ対応するチャネル疑似反響路に印加してチャネル疑似反響信号を生成し、
これらチャネル疑似反響信号を加算して総合疑似反響信号を生成し、
その総合疑似反響信号を、上記各受話信号の再生信号に対する集音信号から差し引いて残留反響信号を生成し、
上記各チャネルの受話信号の再生信号と上記集音信号との間の各インパルス応答を入力して、各チャネルのインパルス応答の級数の係数の二乗和の相対比と各チャネルの更新ステップサイズの相対比が等しくなるように各チャネルの更新ステップサイズを算出し、
上記各チャネルごとに更新ステップサイズと、各チャネルの受話信号と、上記残留反響信号とからそのチャネルの上記再生信号の反響路のインパルス応答を逐次推定し、
その逐次推定したインパルス応答により対応するチャネルの疑似反響路を生成することを特徴とする多チャネルエコーキャンセル方法。
【請求項２】
Ｎチャネル（Ｎは２以上の整数）の受話信号がそれぞれ印加されて、チャネル疑似反響信号を生成するＮ個のチャネル疑似反響路と、
上記Ｎ個のチャネル疑似反響信号を加算して総合疑似信号を生成する合成回路と、
上記Ｎ個の受話信号の再生信号に対する集音信号から上記総合疑似信号を差し引いて残留反響信号を生成する減算回路と、
上記Ｎ個の受話信号の再生信号と上記集音信号との間の各インパルス応答を入力して、各チャネルのインパルス応答の級数の係数の二乗和の相対比と各チャネルの更新ステップサイズの相対比が等しくなるようにＮ個の更新ステップサイズを求めるステップサイズ生成部と、
上記各チャネルごとに上記更新ステップサイズと、各チャネルの受話信号と、上記残留反響信号とが入力され、そのチャネルの上記再生信号の反響路のインパルス応答を逐次推定し、その推定したインパルス応答を対応する上記チャネル疑似反響路に設定するＮ個の反響路推定部と、
を具備する多チャネルエコーキャンセル装置。
【請求項３】
請求項１記載の多チャネルエコーキャンセル方法をコンピュータに実行させるプログラム。
【請求項４】
請求項３記載のプログラムを記録したコンピュータ読み出し可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
この発明は、多チャネル再生系を有する例えば通信会議システムにおいて、ハウリングの原因および聴覚上の障害となる室内反響信号を消去する多チャネル反響消去方法、その装置、そのプログラム及びその記録媒体に関するものである。
【０００２】
【従来の技術】
同時通話性能に優れ反響感の少ない拡声通話システムを提供するために、反響消去装置がある。
まず、１チャネル用の反響消去装置について、その反響消去方法および装置構成を図４を参照して説明する。拡声通話において、受話端子１１に相手の発話等で得られる受話信号は、そのままスピーカ１２から再生される場合と、スピーカ１２へ送る前に、受話信号の振幅やパワー等の大きさに応じて自動的に利得を調節するなど、受話信号に何らかの加工が施された後に、スピーカ１２から再生される場合とがある。このため、この明細書で受話信号ｘ１（ｋ）とは、相手からの受話信号そのものとは限らず、受話信号に対し加工が施された場合は、その加工された後の受話信号を指すものとする。ｋは離散時間を表す。反響消去装置１４は、受話信号ｘ１（ｋ）がスピーカ１２から反響路１５を経て、マイクロホン１６に集音されて得られる反響信号ｙ１（ｋ）を消去する。ここで、反響信号ｙ１（ｋ）は、時刻ｋにおける反響路１５のインパルス応答をｈ１１（ｋ，ｎ）として、
0003616341.tif 000002
Σはｎ＝０からＬ−１まで
のような畳み込み演算で得られるものとモデル化できる。Ｌはタップ数で、反響路１５の残響時間に対応させて、あらかじめ設定しておく定数である。まず、受話信号蓄積・ベクトル生成部１７において、受話信号ｘ１（ｋ）をＬ−１時刻過去のものまで蓄積しておく。蓄積された信号は、受話信号ベクトルｘ１（ｋ）、すなわち
0003616341.tif 000003
として出力される。但し、＊Ｔ  はベクトルの転置を表す。疑似反響信号生成部１８では、式（２）の受話信号ベクトルｘ１（ｋ）と、反響路推定部１９から得られる疑似反響ベクトルｈ＾１１（ｋ）との内積演算
0003616341.tif 000004
を行ない、その結果として、疑似反響信号ｙ＾１（ｋ）を生成する。この内積演算は、式（１）のような畳み込み演算と等価である。反響路推定部１９では、疑似反響信号生成部１８で用いる疑似反響路ベクトルｈ＾１１（ｋ）を生成する。この反響路推定に用いる最も一般的なアルゴリズムは、ＮＬＭＳアルゴリズム（学習同定法）である。ＮＬＭＳアルゴリズムでは、時刻ｋにおける受話信号ベクトルｘ１（ｋ）と、残留反響信号ｅ１（ｋ）、すなわち差回路２１でマイクロホン１６の出力ｙ１（ｋ）から疑似反響信号ｙ＾１（ｋ）を差し引いた信号
0003616341.tif 000005
とから、時刻ｋ＋１において用いる疑似反響路ベクトルｈ＾１１（ｋ＋１）を次式のように求める。
0003616341.tif 000006
但し、μはステップサイズパラメータと呼ばれ０＜μ＜２の範囲で適応動作の調整に用いる。以上のような処理を繰り返すことにより、反響路推定部１９では、次第に疑似反響路ベクトルｈ＾１１（ｋ）を、真の反響路１５のインパルス応答ｈ１１（ｋ，ｎ）の時系列を各要素として持つ反響路ベクトルｈ１１（ｋ）、すなわち
0003616341.tif 000007
と一致させることが可能となり、その結果式（４）の残留反響信号ｅ１  （ｋ）を小さくすることができる。
【０００３】
一般にＮ（＞２）チャネルの再生系とＭ（＞１）チャネルの集音系とで構成される通信会議システムの場合の反響の消去は、図５に示すような構成により行われる。すなわち再生側の全Ｎチャネルと集音側の各１チャネルとの間にＮ入力１出力時系列信号を処理するＮチャネル反響消去装置２２１，２２２，…，２２Ｍ  をそれぞれ接続した反響消去システム２３として実現される。この場合システム全体でＮ×Ｍ個の反響路１５ｎｍ（１＜ｎ＜Ｎ，１＜ｍ＜Ｍ）が存在する。このシステムの構成単位である再生側の全Ｎチャネルと集音側の各１チャネルとの間に接続されるＮチャネル反響消去装置２２１，２２２，…，２２Ｍ  については、図４に示した反響消去装置１４の構成を拡張して、図６に示すように構成される。これは例えば電子情報通信学会論文誌  ’８６／１０  Ｖｏｌ．Ｊ６９−Ａ  Ｎｏ．１０「多チャネル適応ディジタルフィルタ」に詳しく述べられている。ここで、集音側が第ｍ集音チャネル（１＜ｍ＜Ｍ）に接続されているＮチャネル反響消去装置２２ｍ  を考える。第ｍチャネルのマイクロホン１６ｍ  で集音される反響信号は、各再生チャネルの受話信号がそれぞれの反響路１５１ｍ〜１５Ｎｍを経て集音側で全て加算されることにより得られるために、反響路推定をどの再生チャネルについても、統一的に同じ１つの残留反響信号ｅｍ  （ｋ）のみを評価して行なうための工夫が必要となる。まず、各再生チャネルの受話信号について、受話信号蓄積・ベクトル生成部（１７１，１７２，…，１７Ｎ  ）により、受話信号ベクトル
0003616341.tif 000008
を生成する。但し、Ｌ１，Ｌ２，…，ＬＮ  はタップ数で、各反響路１５１ｍ，１５２ｍ，…，１５Ｎｍの残響時間に対応させて、あらかじめ設定する定数である。これらのベクトルをベクトル結合部２４によって、
0003616341.tif 000009
と結合する。また、反響路推定部１９ｍ  においても、各再生チャネルと第ｍ集音チャネルとの間のＮ個の反響路を模疑するための、各疑似反響路ベクトルｈ＾１ｍ（ｋ），ｈ＾２ｍ（ｋ），…，ｈ＾Ｎｍ（ｋ）  を結合して
0003616341.tif 000010
として扱う。疑似反響路結合ベクトルｈ＾ｍ  （ｋ）  の更新は、ＮＬＭＳアルゴリズムを用いた場合、
0003616341.tif 000011
のように行なわれる。疑似反響信号生成部１８ｍ  では、内積演算
0003616341.tif 000012
により、第ｍ集音チャネルで集音された反響信号ｙｍ  （ｋ）に対する疑似反響信号ｙ＾ｍ  （ｋ）を生成する。このように、各チャネル毎のベクトルを結合して１つのベクトルとして扱うことにより、基本的な処理の流れは、図４に示した１チャネル反響消去装置と同様となる。
【０００４】
【発明が解決しようとする課題】
従来の多チャネルエコー消去方法はその疑似反響路が真の反響路にほぼ収束するのが遅い。この点を改善すべく研究した結果、各チャネルの受話信号の再生信号ごとの反響路についての疑似反響路において、その反響路の逐次推定における更新ステップサイズμが同一値である点に問題があると考えるに至った。
説明を簡単にするため、図７に示すように受話チャネルが２つ、集音チャネルが１つの場合を例とする。各受話信号ｘ１  （ｋ），ｘ２  （ｋ）はそれぞれチャネル疑似反響路１８１  ，１８２  へ印加され、チャネル疑似反響路１８１  ，１８２  よりの各チャネル疑似反響信号ｙ＾１  （ｋ），ｙ＾２  （ｋ）は合成回路２６で合成されて総合疑似反響信号ｙ＾（ｋ）とされて差回路２１へ供給され、マイクロホン１６１  よりの集音信号ｙ（ｋ）からｙ＾（ｋ）が引算され、残留反響信号ｅ（ｋ）が求められる。残留反響信号ｅ（ｋ）と、受話信号ｘ１  （ｋ），ｘ２  （ｋ）がチャネル反響路推定部１９１  ，１９２  へ供給され、チャネル反響路推定部１９１  ，１９２  はそれぞれ、受話信号ｘ１  （ｋ），ｘ２  （ｋ）の各再生信号がマイクロホン１６１  に至る反響路１５１  ，１５２  のインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）を逐次推定して、これら推定したインパルス応答ｈ＾１  （ｋ），ｈ＾２  （ｋ）をチャネル疑似反響路１８１  ，１８２  にそれぞれ設定する。
【０００５】
チャネル反響路推定部１９１  ，１９２  における各推定チャネルインパルス応答ｈ＾１  （ｋ），ｈ＾２  （ｋ）は式（１２）から以下のようになる。式（１２）は受話チャネルが２つであるから式（１４）となる。
【０００６】
【数１】
0003616341.tif 000013
【０００７】
反響路１５1 ，１５2 の残響時間により決るタップ数Ｌ1 とＬ2 を等しくＬとした。
このようにＸT（ｋ）Ｘ（ｋ）はスカラーであるから、式（１４）は行列の各行ごとに次のように書ける。
0003616341.tif 000014
つまりチャネル反響路推定部１９1 ，１９2 はそれぞれ式（１５），（１６）によりチャネル反響路１５1 ，１５2 のインパルス応答ｈ1（ｋ），ｈ2（ｋ）を逐次推定することになる。Ｘ1（ｋ）／（ＸT（ｋ）Ｘ（ｋ）），Ｘ2 （ｋ）／（ＸT（ｋ）Ｘ（ｋ））はそれぞれ全チャネルの受話信号のパワーで正規化された受話信号（ベクトル）と言える。ＸT（ｋ）Ｘ（ｋ）はパワー計算部２７で計算されて各推定部１９1 ，１９2 へ供給される。また図７中には図６中の受話信号蓄積ベクトル生成部は必要とするが省略してある。式（１５），（１６）より、その更新ステップサイズは共にμと等しいものとなっている。
【０００８】
しかし、例えば図７に示すように、スピーカ１２１  ，１２２  とマイクロホン１６１  ，１６２  を互いに対向させて、方形の各コーナに位置させた場合、マイクロホン１６１  に対し正面のスピーカ１２１  からの反響路１５１  のインパルス応答ｈ１  （ｋ）の方が、対角位置のスピーカ１２２  からの反響路１５２  のインパルス応答ｈ２  （ｋ）よりエネルギーが大きい。このことは、スピーカとマイクロホンの方形コーナの配置に限らず、他の配置においても反響路によりインパルス応答のエネルギーに差があり、かつ３チャネル以上の多チャネルエコー消去方法についても同様のことが言える。
【０００９】
そこでこの発明の目的はこのインパルス応答のエネルギーの差を考慮して従来よりも収束が速い多チャネルエコーキャンセル方法、その装置、そのプログラム及びその記録媒体を提供することにある。
【００１０】
【課題を解決するための手段】
この発明によれば各チャネルの受話信号をチャネル疑似反響路に印加してチャネル疑似反響信号をそれぞれ生成し、これらチャネル疑似反響信号を加算して総合疑似反響信号を生成し、この総合疑似反響信号を集音信号から差し引いて残留反響信号を生成し、特に各チャネルの受話信号の再生信号と集音信号との各インパルス応答と対応する信号の大きさを入力して、各チャネルのインパルス応答の級数の係数の二乗和の相対比と各チャネルの更新ステップサイズの相対比が等しくなるように各チャネルの更新ステップサイズを生成し、各チャネルごとにその更新ステップサイズと、各受話信号と、上記残留反響信号とからそのチャネルの反響路のインパルス応答を逐次推定し、その推定したインパルス応答によりそのチャネルの疑似反響路を生成する。
【００１１】
【発明の実施の形態】
図１にこの発明の実施形態を示す。説明を簡単にするため、受話２チャネル、集音１チャネルを例とし、図７と対応する部分に同一参照符号を付けてある。この実施形態では、各チャネルの反響路１５１  ，１５２  のインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）と対応する信号をステップサイズ生成部３１に入力して、各チャネルごとのステップサイズμ１  ，μ２  をそのインパルス応答に対応する信号の大きさに応じて決定する。インパルス応答ｈ１  （ｋ），ｈ２  （ｋ）と対応する信号としては受話信号ｘ１  （ｋ），ｘ２  （ｋ）の各再生信号と集音信号ｙ（ｋ）との間の各インパルス応答、つまり反響路１５１  ，１５２  の各インパルス応答ｈ１  （ｋ），ｈ２  （ｋ）それ自体、あるいはスピーカ１２１  ，１２２  、マイクロホン１６１  の配置に応じて経験的に推測される大略の値、またはチャネル疑似反響路１８１  ，１８２  に設定される推定インパルス応答などが用いられる。インパルス応答ｈ１  （ｋ），ｈ２  （ｋ）を用いる場合は、予め、各１つのスピーカのみからインパルス音を放射して、マイクロホン１６の出力信号を測定して求める。
【００１２】
更新ステップサイズμ１  ，μ２  は各チャネルのインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）と対応する信号（以下ｈ１  （ｋ），ｈ２  （ｋ）で代表する）の大きさ、つまり振幅又はパワーに応じ、大きければ大きな値とする。特に各チャネルのインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）のパワーの、これらパワーの和に対する比に応じて更新ステップサイズμ１  ，μ２  を決定する。この場合、これら各チャネルの更新ステップサイズμ１２とμ２２の和が２となるようにする。例えば次のようにして求める。いまｈ１  （ｋ），ｈ２  （ｋ）がそれぞれ１要素のみとして、これらをｈ平面上に書くと図２に示すように、式（１４）の更新ベクトルｈ＾（ｋ）は、ｈ１  （ｋ）とｈ２  （ｋ）のパワーが等しい場合は、横軸に対し４５度のベクトル４１上にあることになる。ベクトルｈ１  （ｋ），ｈ２  （ｋ）の各長さをａ＝√（ｈ１  Ｔ  （ｋ）ｈ１  （ｋ）），ｂ＝√（ｈ２  Ｔ  （ｋ）ｈ２  （ｋ））とすると、更新ベクトルｈ＾（ｋ）は図２中のベクトル４２上のように設定することで、その疑似反響路が真の反響路にほぼ収束するのを速くすることができる。この更新ベクトルｈ＾（ｋ）に応じてステップサイズμ１  ，μ２  を、√（ａ２  ）と√（ｂ２  ）の比と等しくする、つまり
μ１  ：μ２  ＝√（ａ２  ）：√（ｂ２  ）
とすればよい。従って下記の連立方程式が得られる。
【００１３】
μ1 √（ｂ2 ）＝μ2 √（ａ2 ）
μ12＋μ22＝２
この連立方程式を解くと以下のようになる。
μ1 ＝√（２ａ2 ）／（√（ａ2＋ｂ2 ））
μ2 ＝√（２ｂ2 ）／（√（ａ2＋ｂ2 ））
つまり
μ1 ＝√（２ｈ1T（ｋ）ｈ1（ｋ））／（√（ｈ1T（ｋ）ｈ1（ｋ））＋√（ｈ2T（ｋ）ｈ2 （ｋ）））        （17）
μ2 ＝√（２ｈ2T（ｋ）ｈ2（ｋ））／（√（ｈ1T（ｋ）ｈ1（ｋ））＋√（ｈ2T（ｋ）ｈ2 （ｋ）））          （18）
をステップサイズ生成部３１で計算すればよい。
【００１４】
従ってインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）がそれぞれ２乗和回路３２１  ，３２２  に入力されて、各ベクトルの要素の２乗和
ｈ１  Ｔ  （ｋ）ｈ１  （ｋ）＝Σｉ＝０  Ｌ−１  ｈ１２（ｋ，ｉ）
ｈ２  Ｔ  （ｋ）ｈ２  （ｋ）＝Σｉ＝０  Ｌ−１  ｈ２２（ｋ，ｉ）
がそれぞれ計算され、これら２乗和は開平回路３３１  ，３３２  で平方根がとられ、ａ＝√（ｈ１  Ｔ  （ｋ）ｈ１  （ｋ），ｂ＝√（ｈ２  Ｔ  （ｋ）ｈ２  （ｋ））が求められる。これらａとｂが加算回路３４で加算され、式（１７），（１８）の分母が求められる。乗除算回路３５１  で加算回路３４の加算値でａを除算した値に２を乗算して式（１７）の計算結果μ１  が求められる。同様に乗除算回路３５２  で加算回路３４の加算値でｂを除算した値に２を乗算して式（１８）の計算結果μ２  が求められる。
【００１５】
これら更新ステップサイズμ１  ，μ２  がそれぞれ対応するチャネル反響推定部１９１  ，１９２  へ供給される。チャネル反響路推定部１９１  ，１９２  には図６に示した場合と同様にそれぞれ受話信号ｘ１  （ｋ），ｘ２  （ｋ）が供給されると共にパワー計算部２７よりの所定時間の全チャネルの受話信号パワーの和と残留反響信号ｅ（ｋ）が共通に供給される。チャネル反響路推定部１９１  ，１９２  ではそれぞれ以下の式によりチャネル反響路のインパルス応答の推定ｈ＾１  （ｋ＋１），ｈ＾２  （ｋ＋１）が逐次行われる。
【００１６】
0003616341.tif 000015
これら推定されたインパルス応答ｈ＾１  （ｋ＋１），ｈ＾２  （ｋ＋１）がそれぞれチャネル疑似反響路１８１  ，１８２  に設定される。その他は図６に示した装置と同様な処理を行う。
【００１７】
このように更新ステップサイズμ１  ，μ２  が反響路１５１  ，１５２  のインパルス応答ｈ１  （ｋ），ｈ２  （ｋ）のエネルギーの相対比に応じて、大きいエネルギーのチャネルに対応する更新ステップサイズが大きくされ、それだけ全体としての疑似反響路１８１  ，１８２  の収束が速くなる。
ステップサイズ生成部３１で各時刻ｋごとにμ１  （ｋ），μ２  （ｋ）を計算することにより、反響路１５１  ，１５２  が変動してｈ１  （ｋ），ｈ２  （ｋ）が変動してもこれに追従して、適切なμ１  （ｋ），μ２  （ｋ）が設定されることになる。また最初にステップサイズ生成部３１に、インパルス応答ｈ１  （ｋ），ｈ２  （ｋ）の経験的推定値を入力しても、その後、チャネル疑似反響路１８１  ，１８２  の各インパルス応答ｈ＾１  （ｋ），ｈ＾２  （ｋ）をステップサイズ生成部３１に入力することにより、チャネル疑似反響路１８１  ，１８２  が反響路１５１  ，１５２  に収束するに従って、更新ステップサイズμ１  ，μ２  が適切な値μ１  （ｋ），μ２  （ｋ）に更新される。あるいはμ１  ，μ２  の初期値はμ１  ＝μ２  とし従来において使用していた適当な値を設定し、その後は、ｈ＾１  （ｋ），ｈ＾２  （ｋ）をステップサイズ生成部３１に入力してμ１  ，μ２  を更新してもよい。
【００１８】
この発明は受話チャネルが３チャネル以上の場合にも適用できる。例えばＮチャネルの場合、ステップサイズ生成部３１で各チャネルｎ（ｎ＝１，２，・・・，Ｎ）に対し、次の計算をすればよい。
μn＝√Ｎ（ｈnT（ｋ）ｈn（ｋ））／√（ｈ1T（ｋ）ｈ1（ｋ）＋(ｈ2T（ｋ）ｈ2（ｋ）・・・＋ｈNT（ｋ）ｈN（ｋ））
この式から理解されるように、各チャネルのインパルス応答の大きさ（振幅又はパワー）の、チャネル間の相対的大きさ比により、インパルス応答の大きさが大きいチャネルの疑似反響路の収束速度が大きくなるようにステップサイズを求めればよい。また上記実施形態ではチャネル反響路推定部において、全チャネルの受話信号のパワーの和を求めたものを用いたが、必ずしもその必要はない。例えば電子情報通信学会論文誌’８６／１０  Ｖｏｌ．Ｊ６９−Ａ  Ｎｏ．１０「多チャネル適応ディジタルフィルタ」で提案しているように、チャネルごとにその受話信号をその電力で正規化したものを用いてもよい。
【００１９】
またチャネル反響路推定のアルゴリズムとしては学習同定法に限らず、射影法その他のアルゴリズムを用いてもよい。更にこの多チャネルエコーキャンセル装置を、コンピュータによりプログラムを実行させて機能させることもできる。この場合はこの多チャネルエコーキャンセルプログラムを、ＣＤ−ＲＯＭ、フロッピーディスク、磁気ディスクなどからコンピュータのプログラムメモリにインストールし、あるいは通信回線を通じてダウンロードして、そのプログラムを実行することになる。
【００２０】
【発明の効果】
図１に示した実施形態において、マイクロホン１６１  と正面スピーカ１２１  間の距離が約１．５ｍ、スピーカ１２１  と１２２  間の距離が約１．２ｍの場合、インパルス応答ｈ１  （ｋ）とｈ２  （ｋ）のパワー比は壁の反射などにより変わり、１：０．７５〜０．８５程度となるが、１：０．８とした。この場合の計算機シミュレーションの結果を図３に示す。図３中の実線はこの発明方法による疑似反響路の収束状態を示し、破線は従来法による収束状態を示す。反響路を変化させてから係数誤差（インパルス応答推定値誤差）が−１０ｄＢになる時間が、従来法では８．３秒掛ったが、この発明方法によれば５．９秒となり、従来法より約１．５倍の収束速度が得られ、この発明方法が優れていることが確認できた。
【００２１】
以上述べたようにこの発明によれば、各チャネルの反響路のインパルス応答に応じてそのチャネルの反響路推定における更新ステップサイズを決定することにより、多チャネルエコーキャンセルの収束速度を向上させることができる。
【図面の簡単な説明】
【図１】この発明の実施形態の機能構成を示す図。
【図２】インパルス応答に応じてステップサイズを決めるための説明図。
【図３】この発明の効果を確認するための、収束状態の計算機シミュレーションを示す図。
【図４】従来の１チャネルエコーキャンセル装置の機能構成を示す図。
【図５】従来の多チャネルエコーキャンセル装置の構成図。
【図６】従来の多チャネルエコーキャンセル装置の具体的機能構成を示す図。
【図７】従来の２チャネルエコーキャンセル装置の機能構成を示す図。
【図１】
0003616341.tif 000016
【図２】
0003616341.tif 000017
【図３】
0003616341.tif 000018
【図４】
0003616341.tif 000019
【図５】
0003616341.tif 000020
【図６】
0003616341.tif 000021
【図７】
0003616341.tif 000022