(57)【特許請求の範囲】
【請求項１】
撮影装置の手ぶれ予測装置において、
手ぶれ振動を検出し電気信号に変換する検出手段と、
前記検出手段からの時系列信号のうちの所定の数の信号を抽出し、それぞれを多次元空間の座標の各成分とすることによって、時系列信号を多次元空間内の点に対応させる軌跡生成手段と、
前記軌跡生成手段で生成された多次元空間内の複数の点の時系列に基づいて、多次元空間内で時系列的に次に配置される点を予測する予測手段と、
を有することを特徴とする手ぶれ予測装置。
【請求項２】
前記検出手段からの信号を記憶する記憶手段を更に有し、
前記多次元空間内の複数の点と、前記多次元空間内の複数の点よりも時系列的に前の前記記憶手段に格納された信号から生成される多次元空間内の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする、請求項１に記載の手ぶれ予測装置。
【請求項３】
あらかじめ測定された時系列の手ぶれ信号を記憶した記憶手段を更に有し、
前記予測手段は、前記軌跡生成手段によって生成された多次元空間内の複数の点と、前記記憶手段からの出力から生成した多次元空間内の複数の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする、請求項１に記載の手ぶれ予測装置。
【請求項４】
前記検出手段からの出力を平滑化する平滑手段を更に有し、
前記軌跡生成手段は、前記平滑手段からの出力に基づいて多次元空間内の点を生成することを特徴とする、請求項１から請求項３のいずれか１項に記載の手ぶれ予測装置。
【請求項５】
前記記憶手段からの出力を平滑化する平滑手段を更に有し、
前記予測手段は、前記軌跡生成手段によって生成された多次元空間内の複数の点と、前記平滑手段からの出力から生成した多次元空間内の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする、請求項１に記載の手ぶれ予測装置。
【請求項６】
前記予測手段からの出力と当該出力に対する検出手段からの出力とを比較する比較手段と、
前記比較手段からの出力が最小となるように前記所定の数を変化させる次元制御手段と、を更に有することを特徴とする、請求項１から請求項５のいずれか１項に記載の手ぶれ予測装置。
【請求項７】
前記予測手段からの出力と当該出力に対する検出手段からの出力とを比較する比較手段と、
前記比較手段からの出力が最小となるように前記軌跡生成手段から出力されて前記予測手段に入力される多次元空間内の複数の点の個数を制御する手段と、
を更に有することを特徴とする、請求項１から請求項６のいずれか１項に記載の手ぶれ予測装置。
【請求項８】
前記予測手段は、Ｊａｃｏｂｉ法のアルゴリズムに基づいて予測することを特徴とする、請求項１から請求項７のいずれか１項に記載の手ぶれ予測装置。
【請求項９】
前記予測手段は、直交規格化のアルゴリズムに基づいて予測することを特徴とする、請求項１から請求項６のいずれか１項に記載の手ぶれ予測装置。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は、撮影装置の手ぶれによる撮像画質の劣化を防止する装置に関し、特に手ぶれ振動等による撮影装置の移動を補正する手ぶれ予測装置に関する。
【０００２】
【従来技術】
従来よりスチルカメラやビデオカメラの様な撮影装置の手ぶれによる画像の劣化の防止やその画像を改善する機能を有する装置が提案されている。例えば、機械的な振動検出器を用い検出した信号に応じて、撮影装置の光学系をぶれないように補正駆動したり、ＣＣＤなどの電子撮像素子においては、読出しを行うアドレスを補正する等の対策が採れていた。また、画像の撮影速度を向上させたり、ぶれの大きい場合にあっては、その撮影を禁止する等の方法が採用されていた。その他には、手ぶれ状態で撮影された画像を所定の画像処理により精細な画像に再現する方法もある。
【０００３】
【発明が解決しようとする課題】
前者のようなリアルタイムに手ぶれ信号を補正する方式においては、手ぶれ振動に対応して光学系、電子撮像装置に所定の補正処理を施すのであるが、実際には、手ぶれ振動の検出素子や、光学系補正のための機械的な駆動装置の反応時間（  レスポンス）  や、撮像素子までに光路を導く装置の動作時間などがかかるために、手ぶれ検出の実時刻から、撮影装置が撮影動作を行う時の手ぶれ振動を予測して、光学系などの補正を行わなければならない。予測の方法としては、手ぶれ信号を周期的変動として近似して補正を行う方法や、フーリエ変換による特徴的な周波数の「揺らぎ」として近似する方法がある。しかしながら、実際の手ぶれ信号は非周期的な変動である故に、前者の近似による方法は適切とは言えない。
【０００４】
また、フーリエ変換による方式では、「揺らぎ」の時間的情報が失われてしまうので関連する計算が複雑であるにも関わらず、その予測精度は望ましくはない。
【０００５】
これらの補正の方法に対し、手ぶれの検出時刻から予測時刻までの変動を線形な関数と考えて、過去の手ぶれ振動の波形から線形な時間的関数の係数を逐次決定し予測を行う手法が、特願平４−１１２２５  号公報、特願平４−１１３７２  号公報、特願平４−１１３７３  号公報、特願平４−１１２７８０号公報および特願平３−３４２４４１号公報に提案されている。
【０００６】
しかしながら、これらの従来技術の手法では、予測可能な時間が短いとともに、手ぶれ信号そのものに起因するいわゆる「ホワイトノイズ」的な小振幅の揺らぎ成分と手ぶれ振動のセンサの感度や測定系の量子化に伴う測定誤差成分の予測値に対する影響を考慮していない。
【０００７】
以上のように、従来の手ぶれ防止装置にはまだ実用上におけるいろいろな不具合を解決することが依然としてある。本発明の課題は、従来よりも効果的な手ぶれ防止を実現し、カメラ等の撮影装置の手ぶれ振動の検出器の応答時間、光学系の補正装置の動作時間および撮像系の動作時間を十分に考慮することである。
【０００８】
本発明の目的は、不規則な手ぶれ振動の長時間の予測が可能な不規則信号予測機能を有する手ぶれ予測装置を提供することにある。
【０００９】
【課題を解決するための手段】
そこで本発明は上記の課題に鑑みてなされた手ぶれ防止装置である処の手ぶれ予測装置であり、その課題を解決し上述の目的を達成するために次のような手段を講じている。
【００１０】
本発明の予測装置は、撮影装置の手ぶれ予測装置であって、手ぶれ振動を検出し電気信号に変換する検出手段と、前記検出手段からの時系列信号のうちの所定の数の信号を抽出し、それぞれを多次元空間の座標の各成分とすることによって、時系列信号を多次元空間内の点に対応させる軌跡生成手段と、前記軌跡生成手段で生成された多次元空間内の複数の点の時系列に基づいて、多次元空間内で時系列的に次に配置される点を予測する予測手段とを有することを特徴とする。
【００１１】
【作用】
以上の手段を有する手ぶれ予測装置は次のような作用を奏する。
【００１２】
撮影者固有の手ぶれ信号は検出手段により検知され電気信号に変換される。時系列的に抽出された所定数の電気信号は、軌跡生成手段により多次元空間の座標の各成分としてその多次元空間内の点に対応づけられる。そして、予測手段によってこの多次元空間内において時系列的に次に配置されるべき点が予測される。その結果、検知された手ぶれ信号を基にして局所線形化の手法が適用されることにより、所定時間後の点の位置、即ちぶれの状態が予測される。
【００１３】
【実施例】
まず最初に、本発明に適用する不規則時系列の予測のための主な２つの手法について説明する。
【００１４】
（１）「Ｊａｃｏｂｉ法」：
局所線形化予測の原理について解説する。すなわち、
手ぶれ信号検出部によって検出された信号は、直交座標系    ｘ，  ｙ，  ｚ    の１つの座標軸方向の変異（変位）に関してみれば、１次元の時系列    ｓｉ  ，  ｉ＝１，２，３，．．．Ｎより以下のような時系列の組によるベクトルを構成する。
【００１５】
【数１】
0003548268.tif 000002
【００１６】
図５（ａ）  には時系列における手ぶれ信号のｘ軸方向の変位が示されている。上式のｄは「埋め込み次元」と呼ばれるものであり、また、図５（ｂ）  が示すｘｉ  の描く軌道を「アトラクタ」と呼ぶ。これらの図には１次元の時系列から３次元のアトラクタを構成する様子が示されている。このアトラクタ上のｉ番目の点からある経過時間ｍ後のｉ＋ｍ  番目の点までの「写像」を、形式的にｘｉ＋ｍ  ＝ｆｍ  （Ｘｉ  ）  と表現する。予測問題としては、現在の状態から所定時間後の状態への「写像」を表すこの関数の近似解を求めることである。ここでは、１次近似について次に述べる。測定によって記録された過去のデータから、ｘｊ  とｘｊ＋ｍ  が既知であると仮定すると、ｘｊ  の近傍の点ｘｊ  ＋δｘのｍステップ後の状態は級数展開を用いて次式のように表される。
【００１７】
【数２】
0003548268.tif 000003
【００１８】
ここで、Ｄｆｍ（ｘｊ  ）は、ｘｊ  のまわりでの「Ｊａｃｏｂｉ行列」と呼ばれるが、例えばｘが３次元においてｉ，  ｊ，  ｋの３成分をもつものとして写像がこれら各成分ごとに次式のような３つの関数ｆ１  ，ｆ２  ，ｆ３  で表現できると仮定する。
【００１９】
【数３】
0003548268.tif 000004
【００２０】
このとき、Ｊａｃｏｂｉ行列は次のような行列式で表される。
【００２１】
【数４】
0003548268.tif 000005
【００２２】
Ｏ（２）  は、δｘに関して２次以上の高次の項を表している。１次近似ではδｘの値は十分に小さいと仮定して、２次以上の高次項を０と考える。従って、Ｊａｃｏｂｉ行列を近似的に求めることができれば、式（  ２）  を近似的に表すことができる。
【００２３】
ここで図６を用いて説明すると、与えられた時系列で基準となる点ｘｊ  とその点の所定時間後の状態およびｘｊ  の近傍の点とその時間後の状態が既知のものとして与えられている場合に、Ｊａｃｏｂｉ行列を近似的に求める方法を説明する。現在予測したい点をｘＮＯＷ  とする。次にアトラクタ上のｘＮＯＷ  からの最近接点ｘｊ  を探し、これを写像の基準点とする。次にｘｊ  の周り半径εｍａｘ  以内にアトラクタ上の点の集合｛ｘｋｉ｜‖ｘｋｉ−ｘｊ  ‖    ＜    εｍａｘ  ｝を    ｘｋｉ，  ｉ  ＝１，２，・・，Ｎとする。次にこれらの点のｘｊ  に対する偏差ベクトルをｙｉ  ＝ｘｋｉ−ｘｊ  とする。また、ｘｊ  およびｘｋｉをτ＝ｍΔｔ（  但し、Δｔは測定の時間分解能、離散時系列ではτ＝ｍ）だけ時間発展した後の偏差ベクトルをｚｉ  ＝ｘｋｉ＋ｍ−ｘｊ＋ｍ  とする。これらの間の写像を形式的に次式のように表す。
【００２４】
【数５】
0003548268.tif 000006
【００２５】
ここで、Ａｊ  は、以下のように表される行列である。
【００２６】
【数６】
0003548268.tif 000007
【００２７】
ただし、Ｃ，Ｖは「共変行列」と呼ばれるもので、次の二式で表されるｋ行ｌ列の要素を持っている。すなわち、
【数７】
0003548268.tif 000008
【００２８】
【数８】
0003548268.tif 000009
【００２９】
ここで、ｙｉ  ｌ  ，ｚｉ  ｋ  はｙｉ  ，ｚｉ  のｌ  番目およびｋ  番目の成分を表している。
【００３０】
次に、ｘＮＯＷ  とｘｊ  の偏差ベクトルをγとすると、ｘＮＯＷ  のｍ点後の予測ベクトルは次式のように表せる。
【００３１】
【数９】
0003548268.tif 000010
【００３２】
前述の式１の定義によれば、ｘＮＯＷ＋ｍ  の１番目の要素は現在からｍステップ後のスカラ変数の予測値である。
【００３３】
以上のように、Ｊａｃｏｂｉ行列を求めるいわゆる「Ｊａｃｏｂｉ法」によるアルゴリズムを適用することにより上記の予測値を求めることができ、この値を基にして手ぶれの補正制御にも使用できる。
【００３４】
次に、本発明に適用するもう１つの予測手法を説明する。
（２）「正規直交基底法」：
前述の計算手法では、逆行列Ｖ−１の計算が含まれる。この計算は複雑であり、行列Ｖが正則行列であってもしばしば誤差を生じることがある。そこで次に図７（ａ）  〜（ｄ）  を参照しながら、前述の手法を簡素化して逆行列の計算を含まずに予測値を与える当手法について説明する。
【００３５】
まず前述の式（  １）  と同様に、埋め込み次元ｄでアトラクタを構成する。ここで、ｘＮＯＷ  に接近している点をその近い順にｘｐ１，ｘｐ２，ｘｐ３，・・・，ｘｐｋ，ｋ  ＝ｄ  とする。（参照：図７（ａ）  ）
実際にアトラクタ上の接近点を高速に計算する方法としては、先行技術である特願平５−２１２９７６号公報に詳しく教示されている。
【００３６】
次に、これらの点の重心Ｒとする。（参照：図７（ｂ）  ）
重心とｘｐｉ    ｉ＝１，２，３，・・ｄ  の偏差ベクトルをξｉ  、ｘＮＯＷ  の偏差ベクトルをξとする。（参照：図７（ｃ）  ）
ξｉ    ｉ＝  １，２，３，  ・・ｄ  から次の正規直交系を求める。（参照：図７（ｄ）  ）
【数１０】
0003548268.tif 000011
【００３７】
ただし、‖‖はベクトルのノルムであり、〈，  〉は内積である。
【００３８】
次に、直交系とξとの内積を次式のように定義する。
【００３９】
【数１１】
0003548268.tif 000012
【００４０】
この定義により、ξは直交基底Ｖｍ  を用いて次式のように表される。
【００４１】
【数１２】
0003548268.tif 000013
【００４２】
次に、重心および偏差ベクトルのｍステップ後の状態をそれぞれ次のように定義する。
【００４３】
【数１３】
0003548268.tif 000014
【００４４】
ここで、ｘＮＯＷ＋ｍ  だけは未知の点であり、これを求めることが予測問題である。ここで再び正規直交系を求め、次式のように表す。
【００４５】
【数１４】
0003548268.tif 000015
【００４６】
【数１８】
0003548268.tif 000016
【００４７】
【数１５】
0003548268.tif 000017
【００４８】
ｍステップ後を予測したベクトルは次式のように表される。
【００４９】
【数１６】
0003548268.tif 000018
【００５０】
この方法は図８が示すように、ｍステップ後の写像によって正規直交基底がどの様になっているかということを表わして新たな直交系Ｖｉ  と以前の係数の組Ｌｉ  を用いて未知の点の予測値を与えるものである。
【００５１】
以上のように、「正規直交基底法」によるアルゴリズムを適用することにより上記のような予測値を求めることができ、この値を手ぶれの補正制御にも使用できる。
【００５２】
次にその他の手法として、測定値の微小な「揺らぎ」を考慮して、「平滑化」による上述の写像に関する予測値の補正を行う手法について説明する。
【００５３】
測定値の微小な揺らぎをキャンセルする方法には、数値的に移動平均をとる方法、電気的に帯域制限を行うフィルタを利用する方法、等があげられる。ここでは、上記の「局所線形化の手法」に付加できるような方法について述べる。
【００５４】
アトラクタ上でｘｉ  からｘｉ＋ｍ  までの運動ベクトルを平均化する方法として次式に与えられるようなｘ’  ｉ＋ｍ  をｘｉ＋ｍ  の代わりの値として用いる。
【００５５】
【数１７】
0003548268.tif 000019
【００５６】
あるいは、ｘｉ  からｘｉ＋ｍ  までの運動を各成分ごとに追跡してそれをローパスフィルタを用いて周波数の高い成分の揺らぎを遮断した軌跡でのｍ点後の値、ｘ’  ｉ＋ｍ  をｘｉ＋ｍ  の代わりの値として用いる方法が考えられる。この様な平均化操作を備えたＪａｃｏｂｉ法による局所線形化によって行う予測について、図９（  ａ）  〜（  ｅ）  を参照しながら説明する。
【００５７】
まず、図９（  ａ）  に示すように、予測したい点ｘｉ  の近傍の点ｘｏ  を探し出し、そのｘｏ  からｘｏ＋ｍ  までの軌跡を追跡して平均化するような補正値ｘ’  ｏ＋ｍ  を求める。
【００５８】
図９（  ｂ）  に示すように、上記のｘｏ  を基準点としてそのε近傍の点Ｘｏ  （ｋ）  ｋ  ＝１，２，・・，Ｎを探索する。また、図９（  ｃ）  に示すように、ｘｏ  の近傍の点についても前述と同様にｍステップ後までに軌跡を追跡して平均化のための操作を行う。
【００５９】
次に、図９（  ｄ）  に示されるように、上述の手順で構成したｘｏ  とｘｏ  （ｋ）  を用いて、前述の図６の説明の手法によってＪａｃｏｂｉ行列Ａの計算を行う。そして、図９（  ｅ）  の示すように、ｘｉ  とｘｏ  の偏差ベクトルにＪａｃｏｂｉ行列Ａを作用させて求める予測値ｘｉ＋ｍ  を得ることができる。
【００６０】
以下に、図面を参照して本発明に係わる実施例ついて説明する。
【００６１】
（第１実施例）
図１には本発明の第１実施例に係る手ぶれ予測装置の構成が機能ブロック図で示されている。画像撮影装置の手ぶれ防止装置において、第１実施例としての本発明の装置は、例えばジャイロ式角速度センサからなる手ぶれ信号検出部１０１と、検出された信号を記憶する手ぶれ信号記録部１０２と、記録した手ぶれ信号の時系列から指定時間後の手ぶれ信号の状態を予測する予測演算部１０３と、状態予測を基にしてモーターまたはボイスコイルから成る画像撮影装置の光学系の光軸を機械的に補正するアクチュエータ１０５へ動作信号を生成する補正信号演算部１０４から構成されている。
【００６２】
不規則な手ぶれ信号に関わる小信号の「揺らぎ成分」および計測に関わる誤差成分を統計的に処理して平滑化を図ることが考えられるが、この平滑化の操作は、前述の図１が示す第１実施例の構成においては、信号記録部１０２に記録された信号から移動平均などの数値演算によって平滑化を行うこともよい。
【００６３】
この第１実施例によれば、信号記録部１０２が手ぶれ信号検出部１０１によって検出された手ぶれ振動の時系列信号の現在までの状態を記録する。予測演算部１０３に対してあらかじめ、予測を行う時間と局所線形化を行う際の線形化行列の次元と線形化に必要な過去のデータのサンプル数を設定しておく。この予測演算部１０３では記録された信号から現在の状態に近い過去の履歴を選び出して「局所線形化」を行う。（但し、実際の「局所線形化」の詳細な手法については後述する。）
局所線形化の手法により、現在から所定時間後のぶれの状態を予測する。補正信号演算部１０４は予測した値を電気信号に変換して画像撮影装置の撮像光学系１０６の光軸を機械的に補正するアクチュエータ１０５への動作信号を生成する。以上の作用によりこのアクチュエータは光軸を手ぶれによる移動を補正するように調整し、手ぶれの補正制御を行う。
【００６４】
図２に示された構成は前述の第１実施例の１つの変形例に係るものである。この構成の特徴は、信号の「揺らぎ成分」および計測に関わる誤差成分を平滑化するために、角速度センサ１０１と信号記録部１０２との間にアナログ信号の低域ろ過のための平滑化処理部１０７（フィルタ手段）を配設したものである。
【００６５】
図示のように、手ぶれ信号検出部１０１で検出された信号は、この平滑化処理部（例えば高周波フィルタ）１０７を通過することで高周波の成分が除かれ、その信号が信号記録部１０２に送られる。
【００６６】
変形例によれば、平滑化処理部１０７で処理された誤差信号を予測演算部１０３に読み込ませ、後述する所定のアルゴリズムを利用する方法による予測手段のパラメータを変化させて誤差を最小になるように「学習」させる。これにより、撮影者固有の手ぶれ信号の特徴が学習され、信号記録部に一時記憶される。
【００６７】
次に、図３に示された構成は、前述の第１実施例のその他の変形例に係るものである。この構成の特徴は、誤差演算部１０８を前述の角速度センサ１０１と予測演算部１０３との間に配設したことである。予測演算部１０３で予測された数値と実際に計測された手ぶれ信号をこの誤差演算部１０８に送り込み、予測値と実際の計測値との誤差を計算する。そしてこの算出された誤差を予測演算部１０３にフィードバックすることを特徴とする。
【００６８】
この変形例によれば、予測値と実際の計測値との誤差が計算によって求められるので、その誤差値に基づいて制御および予測の補正処理を行う。その結果、補正制御等の精度の向上が図れる。
【００６９】
（第２実施例）
次に、図４には本発明の手ぶれ予測装置の第２実施例に係る構成が機能ブロック図で示されている。
【００７０】
撮像装置が前述の第１実施例のような光学撮像装置ではない電子撮像素子２０１から成る場合に、第１実施例の構成要素であったアクチュエータ１０５の代わりに、電子撮像素子２０１の画素のアドレスを補正する信号を生成する画素補正信号演算部２０５を備えることを特徴とする手ぶれ予測装置である。また、手ぶれ信号検出部を構成する要素として、この電子撮像素子２０１からの信号を基にぶれを検出するセンサとしての動きベクトル検出部２０２を備えている。このような構成の装置において、予測演算部２０４が現在の状態から所定時間後のぶれの状態を予測するとき、過去の手ぶれ振動データの履歴から局所的な「線形射像関数」を生成してその予測を行う。
【００７１】
また、この予測演算部２０４においては、不規則な手ぶれ信号に関わる小信号の「揺らぎ成分」および計測に関わる誤差成分を統計的に処理して平滑化を行い、局所的な線形射像関数を生成して予測を行うことを特徴とする。
【００７２】
本第２実施例によれば、画素補正信号演算部２０５が、画像を取り込む電子撮像素子２０１の画素のアドレスを補正する信号を生成する。また、動きベクトル検出部２０２は、電子撮像素子２０１からの信号を基に手ぶれの検出センサとして機能する。そして予測演算部２０４が撮影者固有の不規則な信号記録部２０３に蓄積された過去からの手ぶれ振動データ（  即ち、履歴）  を基にして、局所的な線形射像関数の生成および、信号の「揺らぎ成分」や計測誤差成分を統計処理して平滑化を図り予測を行う。
【００７３】
（実施例の動作１）
次に、本実施例の装置の動作について図１０〜図１３を参照しながら説明する。図１０には、図１〜図４に示した装置の予測演算部１０３，２０４が行う処理手順が前述の「Ｊａｃｏｂｉ法」に基づいて流れ図で示されている。
【００７４】
まず手ぶれの予測を行うに際して、あらかじめ手ぶれ検出部の応答速度および、撮像系の動作時間を考慮して予測する時間ｍを計算しておき、予測演算部にこの値を設定する（Ｓ１）。
【００７５】
線形化を行う次元を設定する（Ｓ２）。
【００７６】
線形化を行うに当たって必要な過去の履歴の測定点数ｎを設定する（Ｓ３）。
【００７７】
設定した点数の手ぶれの測定を図１〜図４の要素１０１と１０２を含む手ぶれ信号検出部で検出し、図１〜図４の信号記録部１０２，２０２で記録しておく。予測は記録された最も新しいデータｓ（  ｎ）  について行う。式（  １）  の方法によってこのｓ（  ｎ）  について再編成ベクトルを構成する（Ｓ５）。
【００７８】
次に再編成の状態空間でｘｎ  の最近接点を求めるが、近接点を求める具体的方法については、先行技術を開示した特願平５−２１２９７６号公報にも示されるように、スカラ列ｓ（  ｉ）  ，ｉ＝１，２，・・，  ｎ  からソート列を生成して、ベクトルの各要素毎に近接点の探索を行う方法がある。
【００７９】
図１１には、図１０の流れ図中のステップＳ６の詳細手続きを表すサブルーチンが示されている。
【００８０】
まずＳ６１  では、時系列データからソート列を計算する。
【００８１】
Ｓ６２  では、ソート列を参照してｘｎ  の１番目の要素ｓ（  ｎ）  から順次、ソート順位が近い点からの差分を計算してｍ次元のノルムの計算を行う。
【００８２】
Ｓ６３  では、ノルムが最小であるか否かの判断を行い、この最小値とされた場合には処理を終了する。
【００８３】
次に、最近接点ｘｏ  の近接点の集合を求める（Ｓ７）。
図１２には、図１０の流れ図中のステップＳ７の詳細手続きを表すサブルーチンが示されている。
【００８４】
Ｓ７１  では、接近点の基準をεを設定する。
【００８５】
Ｓ７２  以降では、前述のＳ６１  で計算したソート列を参照してｘｏ  の接近点を計算する。
【００８６】
Ｓ７３  では、１番目の要素の差分Ｄｓ  を計算する。前述のＳ６２  と同様にソート順位の近い順番に探索していき、このときの差分Ｄｓ  が基準εに達した時に検索を終了する（Ｓ７４  ）。
【００８７】
またこの手順以降のＳ７５  〜７７  においては、ｄ次元のノルムを計算してその値が基準値εを越えていなければ、その点をｘ０  の近傍の点のリストに加える。
【００８８】
なお、図１０のＳ８では、前述のＳ７で探索された近傍点を用いてＪａｃｏｂｉ行列を計算する。この計算は前記の式（  ５）  〜（  ９）  に従って行う。
【００８９】
図１３には、図１０の流れ図中のステップＳ８の詳細手続きを表すサブルーチンが示されている。
【００９０】
まずＳ８１  では、前述のＳ７で探索した近接点のリストを参照する。
【００９１】
Ｓ８２  では、近接点に対してｘｏ  からの偏差ベクトルｙｉ  ＝ｘｋｉ−ｘｊ  を計算する。
【００９２】
Ｓ８３  では、信号記録部１０２，２０２に記録された過去の信号の履歴を参照して、ｘｏ＋ｍ  と近接点のｍステップ後の時間発展後の偏差ベクトルｚｉ  ＝ｘｋｉ＋ｍ−ｘｊ＋ｍ  を計算する。
【００９３】
Ｓ８４  では、式（  ７）  に従って共変行列Ｃ，Ｖを計算する。
【００９４】
Ｓ８５  では、Ｊａｃｏｂｉ行列    Ａ＝Ｃ・Ｖｉｎｖ  を計算する。ここで、Ｖｉｎｖ  とは逆行列であることを意味する。
【００９５】
Ｓ９では、前述のＳ８で得られたＪａｃｏｂｉ行列をｘｎ  −ｘｏ  に作用させ、式（  ９）  に従って予測ベクトルｘｎ＋ｍ  を得る。
【００９６】
なお、このときのベクトルの１番目の要素は現在のスカラ観測値ｓｎ  のｍステップ後の予測値である。
【００９７】
（実施例の動作２）
本実施例について、図１４〜図１６を参照しながら本装置の動作を説明する。
【００９８】
図１４には、図１〜図４に示した装置の予測演算部１０３，２０４が行う処理手順が前述の「正規直交基底法」に基づいて流れ図で示されている。
【００９９】
まず手ぶれの予測を行うに際して、あらかじめ、手ぶれ検出部の応答速度や撮像系の動作時間を考慮して予測時間ｍを計算しておき、予測演算部にこの値を設定する（Ｓ１）。
【０１００】
また、線形化を行うための次元を設定する（Ｓ２）。
【０１０１】
線形化の処理を行うにあたって必要となる過去からの履歴の測定点数ｎを設定する（Ｓ３）。
【０１０２】
この設定した測定点数の手ぶれの測定を、図１〜４の手ぶれ信号検出部１０１，２０１で検出し、図１の信号記録部１０２，２０２で記録しておく。予測は記録された最も新しいデータｓ（  ｎ）  について式（  １）  の方法により再編成ベクトルを構成する（Ｓ５）。
【０１０３】
また、ｓ（  ｎ）  の最近接点の探索は、図１０のステップＳ６と同じ手続きによって行う。
【０１０４】
なお、このＳ６では、ｄ次元空間でｘｎ  に対して最も近接しているｄ個の点を探索している。
【０１０５】
図１５には、図１４の流れ図中のステップＳ６の詳細手続きを表すサブルーチンが示されている。この例では、探索処理を高速化するために図１０のＳ６と同様にソーティング列を用いている。
【０１０６】
まずＳ６１  では、探索点の候補の数ｑを設定する。このｑの値は正規直交系の次元ｄよりも大きな値が選ばれる。
【０１０７】
Ｓ６２  では、ソーティング列を参照する。
【０１０８】
Ｓ６３  では、接近点の候補の探索を終了する判断手続きである。
【０１０９】
また、Ｓ６４  とＳ６５  では、ｘｎ  の１番目の要素ｓ（  ｎ）  のソート順位を参照してソート順位の近い順に各要素の差分を計算してｄ次元のノルムを計算する。
【０１１０】
Ｓ６６  では、ｘｎ  近接点の候補のインデックスとノルムを対応させたリストを作成する。探索点数がｑを超過したところで検索を終了し、ノルムのソーティングを行い（Ｓ６７  ）、ノルムの小さい順にｄ個の点でリストを作成する。これをｘｐｉ｜ｉ  ＝１，２，，，ｄ  とする。
【０１１１】
次に、図１４のステップＳ７では、正規直交系を用いてｘｎ  のｍステップ後の予測値を計算している。
【０１１２】
図１６には、この図１４の流れ図中のステップＳ７の詳細手続きを表すサブルーチンが示されている。
【０１１３】
Ｓ７１  では、Ｓ６で作成したｘｎ  のｄ個の近傍の点を参照してその重心Ｒを計算している。
【０１１４】
Ｓ７２  では、近傍の点と重心Ｒとの偏差ベクトルξｉ  を計算する。
【０１１５】
Ｓ７３  では、式（  １０）  の基づきξｉ  からｄ次元の正規直交系Ｖを構成している。
【０１１６】
Ｓ７４  では、ξとＶｉ  の内積Ｌｌ  を計算する。
【０１１７】
またＳ７５  では、信号記録部１０２，２０２に記録された過去の信号の履歴を参照してｘｎ  の近傍の点のｍステップ後の点ｘ＊ ｐｉの重心Ｒ＊ を計算する。
【０１１８】
Ｓ７６  では、重心Ｒ＊  とｘ＊ ｐｉの偏差ベクトルξｉ ＊  を計算する。
【０１１９】
Ｓ７７  では、式（  １４）  に従ってξｉ  ＊  の正規直交系Ｖｉ  ＊  を計算する。
【０１２０】
Ｓ７８  では、前述の式（  １５）  に従いＶｉ  ＊  とＬｌ  からξのｍステップ後の状態ξ＊  を予測している。
【０１２１】
Ｓ７９  では、前述の式（  １６）  に従いｘｎ  のｍステップ後の点ｘｎ＋ｍ  を与える。
【０１２２】
なお、この状態ベクトルの１番目の要素は現在のスカラ観測値ｓｎ  のｍステップ後の予測値である。
【０１２３】
（実施例の動作３）
本実施例について、図１７〜図２０を参照しながら装置の動作を説明する。
【０１２４】
図１７には、図１〜図４に示した装置の予測演算部１０３、２０４が行う処理手順を示す図１０の手順に、微少な小信号の「揺らぎ」による予測誤差を低減するために移動平均による履歴データの平滑化手続き処理動作を備えた実施例が示されている。
【０１２５】
ただし、Ｓ１〜Ｓ４の手続きは前述の図１０の説明内容と同じである。
【０１２６】
Ｓ５では、計測データｓ（  ｉ）  を移動平均によって平滑化しｓａ（  ｉ）  を得る。
【０１２７】
Ｓ６では、式（  １）  に従って最新計測データｓ（  ｎ）  のベクトル化を行う。
【０１２８】
Ｓ７では、次に平滑化したデータｓａ（  ｉ）  から再編成した状態空間で、ｘｎ  の最近接点を求める。
【０１２９】
なお図１８には、図１７の流れ図中のステップＳ７の詳細手続きを表すサブルーチンが示されている。
【０１３０】
Ｓ７１  では、平滑化した時系列データｓａ（  ｉ）  とｘｎ の要素とからソート列を計算する。
【０１３１】
Ｓ７２  では、ソート列を参照してｘｎ  の１番目の要素ｓ（  ｎ）  から順次、ソート順位が近い点からの差分を計算してｍ次元のノルムの計算を行う。
【０１３２】
Ｓ７３  では、ノルムが最小か否かの判断を行い、この最小値と判断されたときに処理を終了する。
【０１３３】
次に、平滑化データで構成したアトラクタでｘｎ  に最も近い点ｘｏ  の近接点の集合を求める（Ｓ８）。
図１９には、図１７の流れ図中のこのステップＳ８の詳細手続きを表すサブルーチンが示されている。
【０１３４】
まずＳ８１  では、接近点の基準εを設定する。
【０１３５】
Ｓ８２  以降では、Ｓ７１  で計算されたソート列を参照してｘｏ  の近接点を計算する。
【０１３６】
Ｓ８３  では、１番目の要素の差分Ｄｓ  を計算する。
【０１３７】
また、前述のＳ７２  と同様にソート順位の近い順番に探索していき、このときの差分Ｄｓ  が基準εに到達したとき、一連の探索処理を終了する（Ｓ８４  ）。
【０１３８】
以降の    Ｓ８５  〜８７  では、ｄ次元のノルムを計算して、その値が基準εを超過していなければ、その点をｘ０  の近傍の点のリストに加える。
【０１３９】
次に、図１７のＳ９では、ステップＳ８で探索された近傍点を用いてＪａｃｏｂｉ行列を計算する。この計算は式（  ５）  〜（  ９）  に従って行われる。
【０１４０】
図２０には、図１７の流れ図中のステップＳ９の詳細手続きを表すサブルーチンが示されている。
【０１４１】
まずＳ９１  では、前述のＳ８で探索した近傍点のリストを参照する。
【０１４２】
Ｓ９２  では、近接点に対してｘｏ  からの偏差ベクトルｙｉ  ＝ｘｋｉ−ｘｊ  を計算する。
【０１４３】
Ｓ９３  では、平滑化データ列ｓａ（  ｉ）  を参照する。
【０１４４】
Ｓ９４  では、ｘｏ＋ｍ  と近傍点のｍステップ時間発展後のとの間に生じる偏差ベクトルｚｉ  ＝ｘｋｉ＋ｍ−ｘｊ＋ｍ  を計算する。
【０１４５】
また、Ｓ９５  では定義された式（  ７）  に従い、共変行列Ｃ，Ｖを計算している。
【０１４６】
Ｓ９６  では、Ｊａｃｏｂｉ行列    Ａ＝Ｃ・Ｖｉｎｖ  を計算する。ただし、Ｖｉｎｖ  は逆行列であることを表す。
【０１４７】
次にＳ１０においては、前述のＳ９で得られたＪａｃｏｂｉ行列をｘｎ  −ｘｏ  に作用させ、式（  ９）  に従って予測ベクトルｘｎ＋ｍ  を得る。
【０１４８】
なおこのとき、このベクトルの１番目の要素は、現在のスカラ観測値ｓ（  ｎ）  のｍステップ時間後の予測値である。
【０１４９】
（実施例の動作４）
本実施例について、図２１〜図２２を参照しながら装置の動作を説明する。
【０１５０】
図２１には、図１〜図４に示した装置の予測演算部１０３，２０４が行う処理手順を示す図１０の手順に、微少な小信号の「揺らぎ」による予測誤差を低減するために、「時間発展」の補正を行った局所線形化により予測を行った実施例が示されている。なお、Ｓ１〜Ｓ７の手続きは、前述の図１０の各ステップと同じである。
【０１５１】
Ｓ８では、ｘｏ  およびその近接点の時間発展を、式（  １７）  または、図９に示した手法を用いて補正したＪａｃｏｂｉ行列を式（  ５）  〜（  ９）  に従って計算を行う。
【０１５２】
図２２には、図２１の流れ図中のステップＳ８の詳細手続きを表すサブルーチンが示されている。
【０１５３】
まずＳ８１  では、前述のＳ７で探索された近傍点のリストを参照する。
【０１５４】
Ｓ８２  では、近傍点に対してｘｏ  からの偏差ベクトルｙｉ  ＝ｘｋｉ−ｘｊ  を計算する。
【０１５５】
Ｓ８３  では、信号記録部１０２，２０２に記録された過去の信号の履歴を参照してｘｏ  の経過時間ｍステップの時間発展およびｘｏ  と近傍点ｘｏ  （  ｉ）  の時間発展を式（  １７）  を用いて補正を行い、ｘ’  ｏ  ，ｘ’  ｏ  （  ｉ）  を得る。
【０１５６】
Ｓ８４  では、補正点の偏差ベクトルｚｉ  ＝ｘ’  ｏ  −ｘ’  ｏ  （  ｉ）  を計算する。
【０１５７】
Ｓ８５  では、定義された式（  ７）  に従って共変行列Ｃ，Ｖを計算する。
【０１５８】
Ｓ８６  では、Ｊａｃｏｂｉ行列Ａ＝Ｃ・Ｖｉｎｖ  を計算している。但し、このＶｉｎｖ  は逆行列を表している。
【０１５９】
Ｓ９では、前述のＳ８で得られた補正Ｊａｃｏｂｉ行列をｘｎ  −ｘｏ  に作用させ、式（  ９）  に従って予測ベクトルｘｎ＋ｍ  を得る。なおこのとき、このベクトルの１番目の要素は現在のスカラ観測値ｓｎ  のｍステップ後の予測値である。
【０１６０】
（実施例の動作５）
図２３には、本発明の第１および第２実施例に係わる予測方法に準じた変形として装置の動作が流れ図で示されている。
【０１６１】
前述のＳ６において、経過時間ｍステップ後の予測が求められる。
【０１６２】
Ｓ７では、実際の経過時間ｍステップ後の状態を読み出す。
【０１６３】
Ｓ８では、実測値と予測値の誤差を計算する。ここで、今回の実測値を記憶しておき、Ｓ１０において線形化の次元ｄを変更するＳ２に戻り、再度予測値を計算を行う。
【０１６４】
この様に、線形化の次元ｄを変化させ逐次予測値と実測値の差分が計算される。その結果、使用者の手ぶれの特徴を線形化の次元として記憶することができる。
【０１６５】
また、Ｓ１０において、次元ｄを変更する代わりに、予測に用いる過去のデータのサンプル数ｎの値を変更してもよい。このようにすることで、ｎを変更しながら誤差を最小とする予測値を求めることが可能となると共に、予測に必要なサンプルの最適な数を求めることができる。
【０１６６】
（その他の変形実施例）
なお、手ぶれ振動に関するデータの保存は、多人数の平均データを統計的に収集しＲＯＭ化して記録しておくことも可能である。それにより、統計学上の確率度が向上するために、アトラクタの予測精度が大きくなる効果も生まれる。
【０１６７】
本発明装置に関する実施例としては、以上に列挙した例以外にも本発明の主旨を逸脱しない範囲で種々の変形実施が可能である。
【０１６８】
以上に例示した各実施例のそれぞれには、以下の発明が含まれている。
【０１６９】
（  １）    撮影装置の手ぶれ予測装置において、手ぶれ振動を検出し電気信号に変換する検出手段と、前記検出手段からの時系列信号のうちの所定の数の信号を抽出し、それぞれを多次元空間の座標の各成分とすることによって、時系列信号を多次元空間内の点に対応させる軌跡生成手段と、前記軌跡生成手段で生成された多次元空間内の複数の点の時系列に基づいて、多次元空間内で時系列的に次に配置される点を予測する予測手段とを有することを特徴とする手ぶれ予測装置。
【０１７０】
作用１  ：  検出手段は検出センサとして動き、検出した信号を手ぶれのベクトル検出のための信号に変換する。軌跡生成手段は検出信号を基に時系列的に空間座標化して連続する軌跡を生成し、予測手段はこの軌跡を基に局所線形化を行って現在より所定時間後の状態を予測する。
【０１７１】
効果１  ：  手ぶれ振動の長時間の予測を可能とする手ぶれ予測装置が提供でき、その予測情報に基づき撮影装置の手ぶれ振動に関し補正制御を行って実質的に手ぶれを防止する装置が実現できる。
【０１７２】
（  ２）    前記検出手段からの信号を記憶する記憶手段を更に有し、前記多次元空間内の複数の点と、前記多次元空間内の複数の点よりも時系列的に前の前記記憶手段に格納された信号から生成される多次元空間内の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする（  １）  に記載の手ぶれ予測装置。
【０１７３】
作用２  ：  検出手段で検出された手ぶれ振動の時系列信号の現在までの状態を、記憶手段は撮影者固有の不規則な信号としてあらかじめ履歴データとして蓄積して記録する。予測手段はその履歴データに基づき、局所的な線形射像関数の生成および信号の揺らぎ成分や計測誤差成分を統計処理して予測を行う。つまり予測手段は、生成された手ぶれ信号を基に「線形化」および「ベクトル化」処理を行いアトラクタとしての写像を求めることにより「局所線形予測」を行い手ぶれ信号の予測を行う。
【０１７４】
効果２  ：  検出データの蓄積により、予測の信頼性が向上した手ぶれ予測装置を提供できる。
【０１７５】
（  ３）    あらかじめ測定された時系列の手ぶれ信号を記憶した記憶手段を更に有し、前記予測手段は、前記軌跡生成手段によって生成された多次元空間内の複数の点と、前記記憶手段からの出力から生成した多次元空間内の複数の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする（  １）  に記載の手ぶれ予測装置。
【０１７６】
作用３  ：  記憶手段が撮影者固有の不規則な信号を履歴データとして蓄積して記録する。予測手段は、この検出信号と蓄積されている履歴データを利用して多次元空間内での時系列的な点に基づき、次に配置されるべき点を予測する。
【０１７７】
効果３  ：  記憶手段と予測手段との組合せによって、撮影者の個人的な手ぶれ信号の特徴を学習する手ぶれ予測装置を提供できる。
【０１７８】
（  ４）    前記検出手段からの出力を平滑化する平滑手段を更に有し、前記軌跡生成手段は、前記平滑手段からの出力に基づいて多次元空間内の点を生成することを特徴とする（  １）  から（  ３）  のいずれかに記載の手ぶれ予測装置。
【０１７９】
作用４  ：  手ぶれ信号の移動平均による検出データの平滑化を行うフィルタを用いて高周波成分が除去される。
【０１８０】
効果４  ：  微小な信号の揺らぎによる予測誤差が軽減される。
【０１８１】
（  ５）    前記記憶手段からの出力を平滑化する平滑手段を更に有し、前記予測手段は、前記軌跡生成手段によって生成された多次元空間内の複数の点と、前記平滑手段からの出力から生成した多次元空間内の点とに基づいて、多次元空間内で時系列的に次に配置される点を予測することを特徴とする（  １）  に記載の手ぶれ予測装置。
【０１８２】
作用５  ：  手ぶれ信号の移動平均による履歴データの平滑化を行うフィルタを用いて高周波成分が除去される。
【０１８３】
効果５  ：  微小な信号の揺らぎによる予測誤差が軽減される。
【０１８４】
（  ６）    前記予測手段からの出力と当該出力に対する検出手段からの出力とを比較する比較手段と、前記比較手段からの出力が最小となるように前記所定の数を変化させる次元制御手段とを更に有することを特徴とする（  １）  から（  ５）  のいずれかに記載の手ぶれ予測装置。
【０１８５】
作用６  ：  次元制御手段は、比較手段が出力する差分を所定のパラメータを変化させて最小となるように制御する。つまり、線形化の次元を変化させて逐次予測値と実測値を一致させる。
【０１８６】
効果６  ：  予測値が実測値に一致するように制御することで、実質的に手ぶれ振動によるぶれは抑制できる。
【０１８７】
（  ７）    前記予測手段からの出力と当該出力に対する検出手段からの出力とを比較する比較手段と、前記比較手段からの出力が最小となるように前記軌跡生成手段から出力されて前記予測手段に入力される多次元空間内の複数の点の個数を制御する手段とを更に有することを特徴とする（  １）  から（  ６）  のいずれかに記載の手ぶれ予測装置。
【０１８８】
作用７  ：  点の個数を制御する手段は、比較手段から出力される差分が最小になるような点の個数にその多次元空間内の複数の点を調整し、予測値と実測値を一致させる。
【０１８９】
効果７  ：  予測値が実測値に一致するように制御することで、実質的に手ぶれ振動によるぶれが抑制できる。
【０１９０】
（  ８）    前記予測手段は、Ｊａｃｏｂｉ法のアルゴリズムに基づいて予測することを特徴とする（  １）  から（  ７）  のいずれかに記載の手ぶれ予測装置。
【０１９１】
作用８  ：  Ｊａｃｏｂｉ法のアルゴリズムを利用する手ぶれ振動の時系列に関する局所線形化の手法によって、所定時間後のぶれの状態が導き出せる。
【０１９２】
効果８  ：  手ぶれ振動の長時間の予測が、この周知なＪａｃｏｂｉ法のアルゴリズムの適用によって容易となる。
【０１９３】
（  ９）    前記予測手段は、直交規格化のアルゴリズムに基づいて予測することを特徴とする（  １）  から（  ６）  のいずれかに記載の手ぶれ予測装置。
【０１９４】
作用９  ：  直交規格化のアルゴリズムを利用する手ぶれ振動の時系列に関する局所線形化の手法によって、所定時間後のぶれの状態が導き出せる。
【０１９５】
効果９  ：  手ぶれ振動の長時間の予測が、この周知な直交規格化のアルゴリズムの適用によって容易となる。
【０１９６】
【発明の効果】
以上、本発明全体としては、画像撮影装置のユーザ固有の手ぶれ振動に関するデータの蓄積機能を有し、更に従来に確立された予測のための所定のアルゴリズムとを不規則な振動予測のために適用することにより、手ぶれ振動の長期間にわたる予測が容易に可能となり、しかも予測精度が一層向上された信頼性の高い手ぶれ予測装置を提供できる。
【０１９７】
また、この予測された値を電気信号に変換し画像撮影装置の光軸を機械的に補正するアクチュエータ等への動作信号を基にして処理すれば、アクチュエータは光軸を手ぶれによる移動分補正するように調整して手ぶれ補正を行うこともでき、よって、従来に増して信頼性の高い手ぶれ防止機構をも実現できる。
【図面の簡単な説明】
【図１】第１実施例に係わる手ぶれ予測装置の構成を示す機能ブロック図。
【図２】第１実施例の変形例に係わる構成を示す機能ブロック図。
【図３】第１実施例のその他の変形例に係わる構成を示す機能ブロック図。
【図４】第２実施例に係わる手ぶれ予測装置の構成を示す機能ブロック図。
【図５】（  ａ）  は時系列における手ぶれ信号のｘ軸方向の変位を示すグラフ、（  ｂ）  はｘｉ  の軌道を描くアトラクタ。
【図６】Ｊａｃｏｂｉ行列を近似的に求める方法を示すＪａｃｏｂｉ法の概念図。
【図７】（  ａ）  〜（  ｄ）  はＪａｃｏｂｉ法を簡素化して予測する方法を説明する概念図。
【図８】写像による正規直交基底の変化を表す説明図。
【図９】（  ａ）  〜（  ｅ）  はＪａｃｏｂｉ法による局所線形化による予測方法を説明する概念図。
【図１０】実施例の動作手順１に関するメインルーチンの流れ図。
【図１１】動作手順１のＳ６のサブルーチンの流れ図。
【図１２】動作手順１のＳ７のサブルーチンの流れ図。
【図１３】動作手順１のＳ８のサブルーチンの流れ図。
【図１４】実施例の動作手順２に関するメインルーチンの流れ図。
【図１５】動作手順２のＳ６のサブルーチンの流れ図。
【図１６】動作手順２のＳ７のサブルーチンの流れ図。
【図１７】実施例の動作手順３に関するメインルーチンの流れ図。
【図１８】動作手順３のＳ７のサブルーチンの流れ図。
【図１９】動作手順３のＳ８のサブルーチンの流れ図。
【図２０】動作手順３のＳ９のサブルーチンの流れ図。
【図２１】実施例の動作手順４に関するメインルーチンの流れ図。
【図２２】動作手順のＳ８のサブルーチンの流れ図。
【図２３】実施例の動作手順５に関するメインルーチンの流れ図。
【符号の説明】
１０１…角速度センサ，１０２…信号記録部，１０３…予測演算部，
１０４…補正信号演算部，１０５…アクチュエータ，１０６…撮像光学系，
１０７…平滑化処理部，１０８…誤差演算部，
２０１…電子撮像素子，２０２…動きベクトル検出部，２０３…信号記録部，
２０４…予測演算部，２０５…画素補正信号演算部，２０６…電子撮像素子，
Ｓ１〜Ｓ１０…メイン処理ルーチン，Ｓ６１〜Ｓ６８…Ｓ６のサブルーチン，
Ｓ７１〜Ｓ７９…Ｓ７のサブルーチン，Ｓ８１〜Ｓ８７…Ｓ８のサブルーチン，
Ｓ９１〜Ｓ９６…Ｓ９のサブルーチン。
【図１】
0003548268.tif 000020
【図２】
0003548268.tif 000021
【図３】
0003548268.tif 000022
【図４】
0003548268.tif 000023
【図５】
0003548268.tif 000024
【図６】
0003548268.tif 000025
【図７】
0003548268.tif 000026
【図８】
0003548268.tif 000027
【図９】
0003548268.tif 000028
【図１０】
0003548268.tif 000029
【図１１】
0003548268.tif 000030
【図１２】
0003548268.tif 000031
【図１３】
0003548268.tif 000032
【図１４】
0003548268.tif 000033
【図１５】
0003548268.tif 000034
【図１６】
0003548268.tif 000035
【図１７】
0003548268.tif 000036
【図１８】
0003548268.tif 000037
【図１９】
0003548268.tif 000038
【図２０】
0003548268.tif 000039
【図２１】
0003548268.tif 000040
【図２２】
0003548268.tif 000041
【図２３】
0003548268.tif 000042