(57)【要約】【課題】  事前情報（事前測定）を行うことなく、１入力多出力線形系１１の共通極を比較的簡単に求めることを可能とする。【解決手段】  まず白色雑音を系１１に入力し、系１１の出力端のセンサ１８1 ，１８2 の出力を１サンプルずつ遅延し（Ｓ３−１）、これら遅延信号の自己相関行列ＲA を計算し（Ｓ３−２）、その擬似逆行列ＲA-＝（ＲA ＋δ2 Ｉ）-1を計算し（Ｓ３−３）、また各センサ出力と遅延出力との相互相関行列ＲC を計算し（Ｓ３−４）、ＲA-とＲC の積Ｑ＝ＲA-ＲC を求め、そのδ2  を０に最小化して積行列Ｑを求め（Ｓ３−５）、Ｑの特性多項式を計算し（Ｓ３−６）、そのＡＲ（自己回帰）過程における分母多項式と同等のものを求め、その逆数を求めて共通極とする（Ｓ３−７）、この共通極を用いて、入力信号又はセンサ出力信号に対してフィルタ処理して等化する。【選択図】  図４
2005086774.tif 000002
【特許請求の範囲】
【請求項１】
１入力多出力線形系のＮ個の出力端（Ｎは２以上の整数）に備えられた第１，第２，…，第Ｎのセンサよりの各出力信号にそれぞれ１サンプル遅延を付与する過程と、
上記Ｎ個の遅延された出力信号から自己相関行列を求める過程と、
上記自己相関行列の擬似逆行列を求める過程と、
上記Ｎ個のセンサ出力信号と上記Ｎ個の遅延された出力信号とから相互相関行列を求める過程と、
上記自己相関行列の擬似逆行列と上記相互相関行列との積行列を求める過程と、
上記積行列の特性多項式を求める過程と、
上記特性多項式の逆元を求めて共通極伝達関数を求める過程と、
上記共通極伝達関数により、上記１入力多出力線形系の入力信号又は上記Ｎ個のセンサ出力信号の少くとも１つをフィルタ処理する過程と
を有する伝達特性等化方法。
【請求項２】
１入力多出力線形系のＮ個の出力端（Ｎは２以上の整数）に備えられた第１，第２，…，第Ｎセンサと、
これら第１，第２，…，第Ｎセンサの出力信号にそれぞれ１サンプル遅延を付与する第１，第２，…，第Ｎ遅延器と、
これら第１，第２，…，第Ｎ遅延器の出力信号が入力され、自己相関行列を計算する自己相関計算部と、
上記自己相関行列が入力され、その擬似逆行列が計算される逆行列計算部と、
上記第１，第２，…，第Ｎセンサの各出力信号と上記第１，第２，…，第Ｎ遅延器の出力信号とが入力され、相互相関行列を計算する相互相関計算部と、
上記擬似逆行列と上記相互相関行列とが入力され、これら積行列を計算する乗算部と、
上記積行列が入力され、その特性多項式を計算する多項式計算部と、
上記特性多項式が入力され、その逆元を計算する逆元計算部と、
入力された信号を上記逆元でフィルタ処理する少くとも１つのフィルタと、
上記１入力多出力線形系の入力端に直列に又は上記Ｎ個のセンサの少くとも１つの出力側に直列に上記フィルタを挿入する手段と
を具備する伝達特性等化装置。
【請求項３】
請求項１に記載した伝達特性等化方法の各過程をコンピュータに実行させるための伝達特性等化プログラム。
【請求項４】
請求項３に記載した伝達特性等化プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は例えば、オーディオ機器用、前段増幅器用のイコライザや、ハウリング防止装置、移動無線受信時、ビル、高速道路などに起因するマルチパス歪の領域技術などに用いられ、１入力多出力線形系において、出力信号が系の伝達関数に起因する歪みを伴う場合に、伝達関数固有の事前情報を用いず、伝達関数に共通に存在する極を高精度にブラインド推定し、出力信号を平坦化する、つまり信号伝達経路の特性を等化する方法、及び装置、そのプログラム、記録媒体に関する。
【背景技術】
【０００２】
信号伝達特性等化装置とは、信号伝達経路への入力信号又は出力信号を処理するフィルタを有し、経路の伝達特性を平坦化することにより、前記出力信号の周波数特性を平坦に調節する装置である。例えば、音響等化装置では、複数マイクロホンを用いて音源信号を観測し、それらに共通する極の影響を取り除く方法が提案されている。これは、ある音響系における複数の伝達特性はその音響系固有の極を含むという音響学的な知見に基づいている。
【０００３】
以下、共通極推定方法について説明する。特許文献１に開示されているように、複数の伝達特性Ｈj(ｚ)は、全ての音響伝達特性に共通する共通極Ａ（ｚ）と、音源と受音点の配置が異なるとその伝達特性の異なるＢj(ｚ)とによって次式のように表される。
Ｈj(ｚ)＝Ｂj(ｚ)／Ａ（ｚ）（ｊ＝１，２，…，Ｍ）
この表現は、共通極・零モデルと呼ばれる音響伝達特性の表現である。この共通極Ａ（ｚ）の特性を用いることにより、各伝達特性に共通な振幅周波数特性のピークが抑圧された信号となる。なお、共通極は、特許文献１に示すようにそれと等価な表現である共通ＡＲ（Autoregressive：自己回帰）係数の形で推定される。
【０００４】
従来においては、計測された各伝達特性を用いて、直接的に最小二乗誤差で共通ＡＲ係数を推定している。そのためには、計測されたＭ個のインパルス応答ｈj(ｋ)（ｊ＝１，２，…，Ｍ）と、これに対応する共通極・零モデルのインパルス応答との誤差の二乗和を最小にすればよい。ただし、この誤差を最小にする係数の推定は困難なため、式誤差ｅｒｒj(ｎ)が用いられている。この式誤差の全てのインパルス応答に渡る和は以下の様に表される。
ｅ＝ｈ−Ａｘ
ｅ＝［ｅ1 ，ｅ2 ，…，ｅM ］T ，ｅj ＝［ｅｒｒj(０)，ｅｒｒj(１)，…，ｅｒｒj(Ｎ＋Ｐ)］T 
ｈ＝［ｈ1 ，ｈ2 ，…，ｈM ］T ，ｈj ＝［ｈj(0)，ｈj(1)，…，ｈj(Ｎ)，０，…，０］T 
【数１】
2005086774.tif 000003
［・］T は転置を表わす。
【０００５】
誤差ｅの二乗ノルムを最小化する解はｘ＝（ＡT Ａ）-1ＡT ｈで与えられ、これを計算することで、共通ＡＲ係数ａi(ｉ＝１，２，…，Ｐ)が最小二乗解として求められる（非特許文献１）。
【特許文献１】特公平７−３９９６８号公報、段落［００３１］
【非特許文献１】Y.Haneda,et al,“Multiple-point equalization of room transfer functions by using common acoustical poles”，IEEE Trans.SAP,325-333(1997)。
【発明の開示】
【発明が解決しようとする課題】
【０００６】
非特許文献１に示す方法において複数の伝達特性から共通極を推定するためには、複数の伝達特性が既知である必要がある。そのため、第一の課題として、入力信号と複数の出力信号の同期を取る必要があり、演算やシステム構成が複雑となる点が挙げられる。また、第二の課題として、複数の伝達特性に含まれる共通極の推定は、受音点の位置の違いなどにより影響を受けるという点が挙げられる。
【課題を解決するための手段】
【０００７】
図１Ａに示す１入力２出力線形系１１を参照して、この発明の原理を説明する。この系１１の入力端１２と、２つの出力端１３1 ，１３2 との間の各伝達関数をｈ1(ｚ)，ｈ2(ｚ)とする。これら伝達関数ｈ1(ｚ)，ｈ2(ｚ)は共通な零点を持たないと仮定する。伝達関数ｈ1(ｚ)，ｈ2(ｚ)に共通な極が存在する時、伝達関数ｈ1(ｚ)，ｈ2(ｚ)は以下のように書ける。
ｈ1(ｚ)＝ｇ1(ｚ)／（１−ａ（ｚ）），ｈ2(ｚ)＝ｇ2(ｚ)／（１−ａ（ｚ））
（１）
ここで、１／（１−ａ（ｚ））は共通極に対応するＡＲ多項式を表している。従って白色信号が入力された系１１は図１Ｂに示すように白色信号により駆動されたＡＲ過程信号が１入力２出力線形系１４で伝達関数ｇ1(ｚ)，ｇ2(ｚ)を通過するものと表わせる。すなわち、白色信号を用いて伝達関数ｈ1(ｚ)，ｈ2(ｚ)の共通極をブラインド推定する問題は、ＡＲ過程信号ｕn が伝達関数ｇj(ｚ)（ｊ＝１，２）を持つ伝達系に入力された時に、入力信号のＡＲ過程をブラインド推定する問題に帰着される。この問題は、非特許文献２に記載の方法と同様な手法で解くことができる。
【０００８】
ＡＲ過程で生成された入力信号ｕ（ｎ）のｐ個分の時刻よりなるベクトルｕn は以下の様に表される。
【数２】
2005086774.tif 000004
但し、ｅ（ｎ）：線形予測誤差
出力端１３1 ，１３2 に設けられた各センサの出力をｘ1(ｎ)，ｘ2(ｎ)とすると、
ｘn ＝［ｘ1(ｎ)，…，ｘ1(ｎ−ｍ)，ｘ2(ｎ)，…，ｘ2(ｎ−ｍ)］T        （３）
ｍ：ｍ≧Ｎ−１，ｍ＝ｐ−Ｎを満たす整数
は両出力端１３1 ，１３2 の各センサ出力信号ｘ1(ｎ)，ｘ2(ｎ)とこれらに対し１サンプル遅延を与えた信号ｘ1(ｎ−１)，ｘ2(ｎ−１)との相互相関行列ＲC を計算し、
ＲC ＝Ｅ＜ｘn-1 ｘnT ＞                                          （４）
また前記１サンプル遅延を与えた信号の自己相関行列ＲA を計算する。
ＲA ＝Ｅ＜ｘn-1 ｘn-1T ＞                                          （５）
これら相関行列ＲC とＲA を用いて次の行列Ｑを求める。
Ｑ＝lim（ＲA ＋δ2 Ｉ）-1ＲC                                       （６）
ｌｉｍはδ2 を最小化することを表わし、
δ2 は小さな正数、Ｉは単位行列を表わす。
【０００９】
行列Ｑの特性多項式ｆc(Ｑ)は、式（３）のコンパニオン行列Ｃの特性多項式ｆc(Ｃ)と等しい（M.Miyoshi,“Estimating AR parameter-sets for linear-recurrent signals in convolutive mixtures”，4thInternational Symposium on ICA,585-589(2003)）。コンパニオン行列の性質（例えば非特許文献２：児玉，須田「システム制御のためのマトリクス理論」コロナ社、１０３頁（１９７８）参照）により、特性多項式ｆc(Ｃ)は多項式（１−ａ（ｚ））に等価である。よって次の関数が成立する。
ｆc(Ｑ)＝ｆc(Ｃ)＝１−ａ（ｚ）                                      （７）
つまり共通極（１−ａ（ｚ））がブラインド計算される。
【００１０】
以上の原理に基づき、この発明の方法によれば、その対象とする１入力多出力線形系にまず白色雑音を入力し、その１入力多出力線形系のＮ個の出力端（Ｎは２以上の整数）に備えられた第１，第２，…，第Ｎセンサの各出力信号にそれぞれ１サンプル遅延を付与し、
これら遅延された出力信号から自己相関行列を求め、この自己相関行列の擬似逆行列を求め、前記第１，第２，…，第Ｎセンサの各出力信号と、その各１サンプル遅延信号とから相互相関行列を求め、前記自己相関行列の擬似逆行列と前記相互相関行列との積行列を求め、この積行列の特性多項式を求め、その特性多項式の一部の関数の逆元を求める。次にその逆元により前記１入力多出力線形系の入力信号又はＮ個のセンサ出力信号の少くとも１つに対してフィルタ処理をする。
【発明の効果】
【００１１】
この発明によれば、出力信号のみを用いる構成とすることによって従来の複雑な構成を簡略化し、系の出力端に備える複数センサの各出力を１サンプル遅延させ、これらセンサ出力と遅延信号から求めた多項式の逆元に相当する伝達関数をセットしたフィルタで入力信号を処理することにより、受信点の位置によらない高精度な等化が出来る。
【発明を実施するための最良の形態】
【００１２】
以下、図面を用いてこの発明の実施形態を説明する。図２はこの発明を１入力２出力線形系１１に適用した場合の機能構成を示す。この系１１の入力端１２に、フィルタ１５の出力端と、白色雑音源１０とが切替えスイッチ１６により切替え接続される。フィルタ１５の入力端が信号印加端子１７に接続される。つまり１入力２出力線形系１１の入力信号ｕ（ｎ）を信号印加端子１７よりフィルタ１５を介して印加したり、白色雑音を系１１の入力端１２に印加したりすることができるようにされる。
系１１の出力端１３1 ，１３2 において、入力信号ｕ（ｎ）を検出する２個のセンサ１８1 ，１８2 が設けられ、これらセンサ１８1 ，１８2 の出力信号は遅延器１９1 ，１９2 でそれぞれ１サンプル遅延されるセンサ１８1 ，１８2 の各出力信号と遅延器１９1 ，１９2 の各出力信号はフィルタ計算部２３へ供給される。フィルタ計算部２３は系１１の入力端１２に白色雑音源１０を接続した状態で、系１１の共通極を次のように求める。
【００１３】
フィルタ計算部２３の機能構成例を図３に示す。
遅延器１９1 ，１９2 の各出力信号より自己相関計算部２４で式（５）の自己相関行列ＲA が計算され、この自己相関行列ＲA が擬似逆行列計算部２５に入力されてＲA擬似逆行列、つまり式（６）中の行列（ＲA ＋δ2 Ｉ）-1が計算される。遅延器１９1 ，１９2 の各出力信号とセンサ１８1 ，１８2 の各出力信号とが相互相関計算部２６に入力されてこれらの相互相関行列ＲC が式（４）により計算される。
上記行列（ＲA ＋δ2 Ｉ）-1と相互相関行列ＲC との積行列Ｑ＝ (ＲA ＋δ2 Ｉ)-1ＲCが乗算部２７で計算され、そのδ2 が最小化され、つまり式（６）の行列Ｑが計算される。この積行列Ｑの特性多項式ｆc(Ｑ)が多項式計算部２８で計算され、その特性多項式中の、ＡＲ生成過程の分母多項式と同等な多項式（１−ａ（ｚ））の逆元が逆元計算部２９で計算される。これが系１１の共通極１／（１−ａ（ｚ））として用いられる。
【００１４】
この共通極１／（１−ａ（ｚ））、つまり伝達関数がフィルタ１５に設定される。この状態でスイッチ１６がフィルタ１５側に切替えられる。従って、入力端１７からの入力信号はフィルタ１５により、フィルタ処理されて、系１１における共通極の影響が除去され、つまり等化される。なお、図２中に破線で示すように、フィルタ１５を省略して、スイッチにより白色雑音源１０と信号印加端子１７とを入力端子１２に切替え接続できるようにし、センサ１８1 ，１８2 に必要に応じてスイッチ１６1 ，１６2 を通じてフィルタ１５1 ，１５2 をそれぞれ接続し、これらフィルタ１５1 ，１５2 に、フィルタ計算部２３で求めた共通極１／（１−ａ（ｚ））を設定し、センサ１８1 ，１８2 の出力信号をそれぞれフィルタ処理して系１１の共通極による影響を等化してもよい。なお、必ずしも全てのセンサ出力信号に対し、共通極の等化を行う必要がない場合は、一部のセンサ出力信号、図２ではセンサ１８1 ，１８2 の２つのみであるからその一方の出力信号のみ、一般には１乃至複数のセンサ出力信号に対して共通極の等化を行うようにしてもよい。
【００１５】
前記図３に示す処理により共通極が得られることは次のように説明することもできる。線形系１１における入力端１２と出力端１３j(ｊ＝１，２)間のインパルス応答を式（８）で表わし、これを用いて下記のように行列を定義する。
【数３】
2005086774.tif 000005
遅延器１９1 ，１９2 の出力信号の自己相関行列ＲA を求めると、次式（１０）が求まる。
ＲA ＝ＨＥ＜ｕn-1ｕn-1T ＞ＨT                                     （１０）
自己相関計算部２４で自己相関行列ＲA を求め、その擬似逆行列ＲA-を擬似逆行列計算部２５で、擬似逆行列の定義式により求めると次式（１１）を求めたことになる。
ＲA-＝（ＲA ＋δ2 Ｉ）-1＝（ＨＥ＜ｕn-1ｕn-1T ＞ＨT ＋δ2 Ｉ）-1
（１１）
【００１６】
また遅延器１９1 ，１９2 及びセンサ１８1 ，１８2 の各出力信号より、相互相関行列ＲC を相互相関計算部２６で計算すると、次式（１２）を求めたことになる。
ＲC ＝ＨＥ＜ｕn-1ｕnT ＞ＨT                                     （１２）
よって、乗算部２７より得られる積行列Ｑは次式（１３）を求めたことになる。
Ｑ＝ (ＲA ＋δ2 Ｉ)-1ＲC＝（ＨＥ＜ｕn-1ｕn-1T＞ＨT ＋δ2 Ｉ）-1ＨＥ＜ｕn-1ｕnT＞ＨT                                                  （１３）
δ2 を十分小さな正数とすることにより、式（１３）の関係は非特許文献２の３３９頁の記載より次のようになる。
Ｑ＝lim (ＨＥ＜ｕn-1ｕn-1T＞ＨT ＋δ2 Ｉ)-1ＨＥ＜ｕn-1ｕnT＞ＨT＝Ｈ（ＨT Ｈ）-1ＣTＨT                                                  （１４）
lim はδ2 を０に最小化する。
【００１７】
積行列Ｑの特性方程式ｆc(Ｑ)よりＡＲ過程の分母多項式１−ａ（ｚ）を求める。ここで、積行列Ｑの零でない固有値λi ＜Ｑ＞は次の関係を満足する（例えば非特許文献２の１６６〜１６７頁参照）。
λi ＜Ｑ＞＝λi ＜Ｈ（ＨT Ｈ）-1  ＣTＨT ＞＝λi ＜ＨT Ｈ（ＨTＨ）-1ＣT＞＝λi ＜Ｃ＞                                                    （１５）
ｉ＝１，２，…，ｐ＜ｍ＋１について
行列Ｃの特性多項式ｆc(Ｃ)は、式（２）より、次の様に表される。
ｆc(Ｃ)＝ｚp −（ａ1 ｚp-1 ＋ａ2 ｚp-2 ＋…＋ａp ）＝ｚp ｛１−（ａ1 ｚ-1＋ａ2 ｚ-2＋…＋ａpｚ-p ）｝＝ｚp ｛１−ａ（ｚ）｝                  （１６）
式（１５）の関係より、積行列Ｑの特性多項式ｆc(Ｑ)は次の様に表される。
ｆc(Ｑ)＝±ｚ2(L+1)-pｚp｛１−ａ（ｚ）｝＝±ｚ2(L+1)｛１−ａ（ｚ）｝
（１７）
上記特性多項式ｆc(Ｑ)にｚ-2(L+1)を乗ずる事より、ＡＲ過程の分母多項式に同等な多項式±｛１−ａ（ｚ）｝を得る。即ち、伝達系の共通極は高精度に推定される。
尚、この実施形態は１入力２出力線形系を対象としているが、１入力Ｎ（Ｎは２以上の整数）出力線形系を対象とする際には、センサ及び遅延器の数をそれぞれＮ個として、同様の処理を行えばよい。
【００１８】
次にこの発明の方法を図２に示した１入力２出力線形系１１に適用した実施形態を図４を参照して簡単に説明する。この方法の各手順を、図５に簡略に示すコンピュータにより実行させることを想定する。
まず系１１の入力端子１２に白色雑音を印加し、センサ１８1 ，１８2 の出力信号を図に示してないが通常はそれぞれ一定周期で標本化してディジタル信号系列としてインタフェース３１を通して記憶部３２に一旦取り込む（Ｓ１）。これらセンサ出力信号を用いてまず共通極の推定を行う（Ｓ２）。このため各センサ出力信号に対し、それぞれ１サンプルの遅延を付与し（Ｓ３−１）、例えば図３に示した処理と同様に両遅延信号の自己相関行列ＲA を計算し（Ｓ３−２）、更にその擬似逆行列ＲA- ＝（ＲA ＋δ2 Ｉ）-1を計算して一時記憶する（Ｓ３−３）。また両センサ出力信号と両遅延信号との相互相関行列ＲC を計算し（Ｓ３−４）、相互相関行列ＲC と擬似逆行列ＲA -との積行列を計算し、そのδ2 を最小化して行列Ｑを求める（Ｓ３−５）。
積行列Ｑの特性多項式ｆc(Ｑ)＝（１−ａ（ｚ））を計算し（Ｓ３−６）、その逆元（１−ａ（ｚ））を計算し、これを系１１の共通極として記憶する（Ｓ３−７）。
【００１９】
このようにして共通極を推定した後、入力信号又は少なくとも１つのセンサ出力信号に対し共通極１／（１−ａ（ｚ））によりフィルタ処理し（Ｓ４）、その結果を出力部３３を通じて出力する（Ｓ５）。
この図４に示した各過程をコンピュータに実行させるための伝達特性等化プログラムをＣＤ−ＲＯＭ、磁気ディスクなどの記録媒体から、あるいは通信回線を介して図５中のプログラムメモリ３４にダウンロードして、中央演算処理部（ＣＰＵ）により、プログラムメモリ３４に格納されているプログラムを実行させればよい。
【００２０】
次の様な伝達経路について、共通極ブラインド計算シミュレーションを行った。
ｈ1(ｚ)＝（１−２ｚ-1＋２ｚ-2）／（１−ａ（ｚ））                （１８ａ）
ｈ2(ｚ)＝（ｚ-1−３ｚ-2＋２．２５ｚ-3）／（１−ａ（ｚ））        （１８ｂ）
ただし、共通極がない場合（１−ａ（ｚ））＝１、共通極がある場合（１−ａ（ｚ））＝１−ｚ-1＋０．５ｚ-2（共通極配置はｚ＝０．５√２ｅ±ｊπ／４）を用いた。又、［背景技術］の項で述べた従来方法と比較するため、線形系１１への入力信号は分散１の定常白色雑音とした。
共通極が存在しない場合の２つの伝達系の振幅周波数特性を図６Ａに、抽出されたＡＲ多項式（共通極）の振幅周波数特性を図６Ｂにそれぞれ示す。横軸は０〜πに正規化した角周波数ωである。共通極が存在しない場合、従来の共通極計算方法では二乗誤差最小に基づいて共通特性を求めているため、零点の影響を受けた特性が得られている。この発明方法による特性は平坦であり、共通極がないことを示している。
【００２１】
共通極が存在する場合の２つの伝達系の振幅周波数特性を図７Ａに示す。抽出されたＡＲ多項式（共通極）の振幅周波数特性を図７Ｂに示す。伝達系ｈ1(ｚ)の特性は零点の影響によりほぼ平坦になっている。抽出された特性においては、およそπ／４付近にピークが見られ、極が抽出されていることが分かる。図８に抽出されたＡＲ多項式（共通極）の各係数の値を示す。この発明方法の場合、精度良くＡＲ係数ａi(ｉ＝１，２，…)の値が求まっていることが分かる。
この発明は室内の音声収録、エコー消去などに限らず、移動無線におけるマルチパスによる影響の除去にも適用でき、この場合はセンサとしてアンテナが用いられ、そのアンテナ出力信号をベースバンド信号に変換して前述と同様の処理を行えばよい。その他同様の問題にもこの発明を適用することができる。
【図面の簡単な説明】
【００２２】
【図１】Ａは共通極をもつ１入力２出力線形系１１を示す図、Ｂはその線形系１１を、共通極と、共通極をもたない１入力２出力線形系とに分離して示した図である。
【図２】この発明装置の実施形態の機能構成を示す図。
【図３】図２中のフィルタ計算部における処理の機能構成例を示す図。
【図４】この発明方法の実施形態を示す流れ図。
【図５】この発明装置をコンピュータにより機能させるための簡略構成を示す図。
【図６】共通極をもたない伝達関数ｈ1(ｚ)，ｈ2(ｚ)の振幅周波数特性、その抽出されたＡＲ多項式（共通極）の振幅周波数特性を示す図。
【図７】共通極をもつ伝達関数ｈ1(ｚ)，ｈ2(ｚ)の振幅周波数特性、その抽出されたＡＲ多項式（共通極）の振幅周波数特性を示す図。
【図８】ＡＲ係数の真値、従来法により抽出された値、この発明方法により抽出された値を示す図。
【図１】
2005086774.tif 000006
【図２】
2005086774.tif 000007
【図３】
2005086774.tif 000008
【図４】
2005086774.tif 000009
【図５】
2005086774.tif 000010
【図６】
2005086774.tif 000011
【図７】
2005086774.tif 000012
【図８】
2005086774.tif 000013