JP6539235B2 - 信号処理装置およびリスケーリング方法 - Google Patents

Description

本発明は、信号処理装置およびリスケーリング方法に関し、より詳細には、波長掃引光源（出力光の波長が時間とともに変動する光源）を用いた光干渉断層撮影法において、波長掃引光源が時間に対して出力光の波数が一次関数的に変化する（波数リニア）光源でない場合に、波数リニア光源を使用した時と同等の干渉波形信号に変換するリスケーリング方法およびこれを行う信号処理装置に関する。

光干渉断層撮影法（ＯＣＴ：Optical Coherence Tomography）は、赤外線の干渉を用いて断層画像を取得する方法であり、眼底検査などの医療分野で用いられている。赤外線光源として広帯域光源を用いて、多くの波長を含んだ光を対象の物体に照射し、反射した光をＣＣＤカメラ等の受光素子で取得する。反射した光は、波長によって位相のずれ方が変わるため、各々の波長における干渉強度が異なるので、分光器を通して波長分解することにより、各波長ごとの干渉信号を取得することができる。この干渉信号をフーリエ変換することにより、物体の反射面の情報を得ることができる。この方法は、ＳＤ−ＯＣＴ（Spectral Domain-OCT）と呼ばれている。また、全ての波長の光を同時に入射するのではなく、光源からの出射光の波長掃引を行って、波長ごとに干渉波形信号を時間的に取得することも行われている。この方法は、ＳＳ−ＯＣＴ（Swept Source-OCT）と呼ばれている。ＳＳ−ＯＣＴは、分光器が不要で、ＣＣＤカメラを用いる必要がないので、装置を小型化でき、ＳＤ−ＯＣＴよりも高速にデータを取得することができる点で優れている。

図１に、従来のＳＳ−ＯＣＴ装置の基本構成の一例を示す。ＳＳ−ＯＣＴ装置は、波長掃引光源２００１、干渉計２００２および信号処理装置２００４から構成され、サンプル２００３内の深さ方向の情報を得ることができる。波長掃引光源２００１は、時間とともに発振波長が連続変化する光源である。干渉計２００２は、ビームスプリッタ２００５、参照光ミラー２００６および受光素子２００７から構成される。

波長掃引光源２００１から出力された出射光２００９は、ビームスプリッタ２００５でサンプル光２０１０と参照光２０１２とに分けられる。サンプル光２０１０が通る光路をサンプルアーム、参照光２０１２が通る光路を参照アームと呼ぶ。サンプル２００３内に光が反射する面２００８がある場合、サンプル光２０１０は、反射面２００８にて反射され、ビームスプリッタを通って受光素子２００７に入射される。一方、参照光２０１２は参照光ミラー２００６にて反射され、ビームスプリッタ２００５を経て受光素子２００７に入射される。

受光素子２００７では、サンプル光２０１０および参照光２０１２の干渉した光２０１３を受光する。受光素子２００７は、干渉光２０１３を光電変換し、干渉信号２０１４を出力する。信号処理装置２００４は、干渉信号２０１４を信号処理する。干渉信号２０１４は、基準面２０１１から反射面２００８までの光学距離ｚに比例した周波数成分と、反射面２００８からの光強度に比例した強度とを有する。信号処理装置２００４における信号処理により、周波数解析を行うことによって、サンプル２００３の深さ方向（サンプル光の光軸方向）の反射面の情報を得ることができる。

信号処理装置２００４で実行する信号処理について、以下に説明する。ここで、図１に示す光学系において若干の定義を行う。サンプル光側において、ビームスプリッタ２００５と参照光ミラー２００６との間の光路長ｌと同じ距離に、便宜的に仮想的な面を定義することとし、これを基準面２０１１と呼ぶ。つまり、ビームスプリッタ２００５と基準面２０１１との間の光路長もｌである。基準面２０１１とサンプル内反射面２００８との光学距離をｚとする。このとき、干渉信号２０１４をｉ_m（ｔ）で表すと、ｉ_m（ｔ）は下記の式で表される。

ただし、ηは受光素子２００７の量子効率、Γ（・）は光源出射光２００９のコヒーレンス関数、Ｐ_r（・）とＰ_s（・）はそれぞれ参照光強度とサンプル光強度、ｋ（・）は出射光２００９の瞬時の波数である。ここで、波数ｋ（ｔ）が時間に対して一次関数的に変化すると仮定する。つまり、波数の時間変化が一定であるとする。このとき、波数の時間変化をｋ’とし、ｔ＝０の時の波数がｋ₀とすると、波数が下記の式で表される。

このような性質を持つ波長掃引光源を波数リニア光源と呼ぶ。この場合、干渉信号２０１４の強度ｉ_m（ｔ）の周波数ｆは下記の式で表される。

波数リニア光源を用いたＳＳ−ＯＣＴ装置では、サンプル内の位置ｚにある反射面からの反射光と参照ミラーからの反射光との干渉によって生じた干渉信号２０１４の強度ｉ_m（ｔ）は、その強度変調の周波数ｆがｚに比例する。従って、信号処理装置２００４で干渉信号２０１４を時間ｔに対してフーリエ変換すると、フーリエ変換の結果の信号には、ｚに比例する周波数に鋭いピークが現れる。

次に、波長掃引光源２００１が波数リニア光源でない場合、つまり、波数の時間変化が以下の式で表される場合を考える。

ここで、ｋ_Ω（・)は時間によって変化する関数である。このとき、干渉信号２０１４をｉ_m（ｔ）のフーリエ変換の結果Ｉ_m（ｆ）は、下記のように表される。

ここで、Ｆ［・］はフーリエ変換、＊は畳み込み積分、ｊは虚数単位を表す。

このようにＩ_m（ｆ）は、Ｆ［ｅｘｐ（±ｊ２ｚｋ_Ω（ｔ)）］で表され、±ｚｋ’／πだけずれた信号と、Ｆ[（Ｐｒ（ｔ）Ｐｓ（ｔ））^1/2]との畳み込み積分の結果となる。
波数リニア光源では、ｋ_Ω（・)が時間に対して変化が無いので、Ｆ［ｅｘｐ（±ｊ２ｚｋ_Ω（ｔ)）］はｆ＝０にピークを持つデルタ関数となる。このため、Ｉ_m（ｆ）は、Ｆ[（Ｐｒ（ｔ）Ｐｓ（ｔ））^1/2]を±ｚｋ’／πだけずらしたものとなる。Ｐｒ（ｔ）Ｐｓ（ｔ）がガウシアンの場合は、Ｉ_m（ｆ）は、±ｚｋ’／πにピークを持つ鋭いガウシアンとなる。Ｐｒ（ｔ）Ｐｓ（ｔ）が時間に対して広いほど、Ｉ_m（ｆ）の±ｚｋ’／πにピークを持つそれぞれのガウシアンの幅は狭くなり、画像の精細度が高くなる。

このような、ある一つの反射点に対するＩ_m（ｆ）をＰＳＦ（point spread function）と呼ぶ。ＳＳ−ＯＣＴ装置で得られる断層画像は、各反射点に対するＰＳＦの和であり、ＰＳＦの幅が狭いほど精細な画像が得られる。波数リニアではない光源では、干渉信号ｉ_m（ｔ）をフーリエ変換すると、上述の通り、ＰＳＦの幅が大きくなり断層画像の精細度が落ちてしまう。そこで、信号処理により、波数リニア光源を使った場合と同等の干渉信号になるように、ｉ_m（ｔ）を時間軸において伸縮する方法が取られる。この処理をリスケ―リングという。

リスケーリングの原理は、波数間隔が一定となる時間間隔で干渉信号２０１４であるｉ_m（ｔ）をサンプリングすることである。つまり、ｎ回目のサンプリング時間をｔ_nとしたときに、下記の式に示す通り、波数の間隔δｋが一定となるように、サンプリングする時刻ｔ_nを決めればよい。

ここで、δｋはｔ_nによらず一定である。

式（１）をみると、波数ｋはｃｏｓ関数の引数に組み込まれているため、ｉ_m（ｔ）の位相φ_m（ｔ）が観測できれば、下記の式で波数ｋが計算できる。

ここで、ｉ_m（ｔ）の位相φ_m（ｔ）の算出の方法について説明する。ｉ_m（ｔ）をフーリエ変換したＩ_m（ｆ）である式（５）の第２項を削除したものＩ_m ⁺（ｆ）を逆フーリエ変換した式ｉ_m’（ｔ）は、

となる。式（８）を見れば、ｉ_m’（ｔ）は複素数であるので、その偏角を求めれば、下記のようにφ_m（ｔ）が求められる。

ａｒｇ（・）は通常−π〜＋π（または、０〜２π）ラジアン等、２πラジアンの情報しか持たないが、位相が連続していると仮定すると、−πから＋πラジアン方向（０から２πラジアン方向）に算出した後に−π（０）ラジアン近辺に戻った時には、２π＋ａｒｇ（ｉ_m’（ｔ））というように、位相をつなげることができる。このような操作をアンラップ（unwrapping）という。式（９）（１０）を計算する場合は、位相をアンラップする必要がある。式（６）の関係となるような干渉信号のサンプリング時刻ｔ_nを得るためには、式（１０）でｋ（ｔ）を求めるときのｚは一定にしておく必要がある。実際にｋ（ｔ）を求める際には、サンプル光側のｚの位置に反射面（鏡など）を固定し、干渉信号ｉ_m（ｔ）を取得して信号処理する。

ｉ_m（ｔ）をフーリエ変換したＩ_m（ｆ）である式（５）の第２項を削除した後に逆フーリエ変換して、式（８）で示すｉ_m’（ｔ）を得るための正確な方法は複雑であるため、近似的な方法として、Ｉ_m（ｆ）の負の周波数成分を削除することが考えられる。Ｉ_m（ｆ）の負の周波数成分を削除したものをＩ_m ^+a（ｆ）とする。このようにすると、δ関数に畳み込まれた成分の内、正の周波数を持つδ関数に畳み込まれた成分の内、負の周波数へ漏れた部分が残らずに削除され、負の周波数成分を持つδ関数に畳み込まれた成分の内、正の周波数成分へ漏れた部分が削除されずに残るため、式（８）に示すような信号が正確に得られなくなる。その結果、正確なｋ（ｔ）が得られないことから、正確なｔ_nを得ることはできなくなる。しかし、Ｉ_m ⁺（ｆ）≒Ｉ_m ^+a（ｆ）と近似できる範囲では、うまくいくことが多いので、現状、この方法が多く取られる。

このような方法で時刻ｔ_nを求め、その時刻ｔ_nで干渉信号ｉ_m（ｔ）をサンプリングすることにより、波数リニア光源を使った場合と同等の干渉信号を得ることが可能となる。Ａ／Ｄ変換によって干渉信号ｉ_m（ｔ）をサンプリングした後に時刻ｔ_nで再度サンプリングすることから、リスケーリングにおけるサンプリングのことをリサンプリングと呼ぶ。なお、φ_m（ｔ）＝ａｒｇ（ｉ_m’（ｔ））を２ｚで割ってｋ（ｔ）を求めなくても、２ｚは一定と考えれば、下記の式によってｔ_nを求めることができる。

ただし、δφ_mは一定である。ｔ_nを求める際には、φ_m（ｔ）の逆関数ｔ（φ_m）を求め、任意のδφ_mを決めて、下記のように求めても良い。

ここで、φ₀は０番目のリサンプリング点（時刻）ｔ₀の位相である。

以上の原理を使って、実際にリスケーリングを行うための装置の構成と処理のフローを説明する。リスケーリング処理の動作は、２段階に分けられ、最初はリサンプリング点（時刻）を取得する段階であり、次に、観測対象のサンプルから得られた干渉信号をリスケーリングする段階である。図２に、従来のリスケーリング装置の機能ブロックを示し、図３に、従来のリスケーリング処理のフローチャートを示す。

最初に、干渉信号をリスケールするためのリサンプリング点（時刻）を取得する段階について説明する。

［ステップＳ−２２０１（干渉波形信号取得１）］サンプル光側のｚの位置に鏡などの反射面を置き、干渉計２００２で生成した干渉信号を、干渉波形信号取得手段２１０１が取得して出力する。受光素子２００７がアナログ出力であり、そのアナログ出力の信号を干渉波形信号取得手段２１０１が受ける場合は、Ａ／Ｄ変換を行い、ディジタル信号として出力する。

干渉波形信号取得手段２１０１が干渉信号を取得するタイミングは、波長掃引光源が短波長から長波長へ、または、逆に長波長から短波長へ掃引を開始するタイミングとする。例えば、波長掃引光源から掃引開始時刻を示すトリガ信号を、干渉波形信号取得手段２１０１が得るようにしても良い。または、トリガ信号を出力するためのトリガ信号出力手段を別途用意し、波長掃引光源にトリガ信号を入力して、トリガ信号のタイミングで（あるいは同期して）波長掃引光源が掃引を開始しても良い。

［ステップＳ−２２０２（離散フーリエ変換）］離散フーリエ変換手段２１０２は、干渉波形取得手段から得た干渉信号を離散フーリエ変換して出力する。

［ステップＳ−２２０３（マイナス周波数領域除去）］マイナス周波数領域除去手段２１０３は、離散フーリエ変換した干渉信号のうち、負の周波数を除去し、除去した信号を出力する。

［ステップＳ−２２０４（離散フーリエ逆変換）］離散フーリエ逆変換手段２１０４は、負の周波数を除去した信号を逆フーリエ変換して出力する。

［ステップＳ−２２０５（位相変化曲線算出）］位相変化曲線算出手段２１０５は、逆フーリエ変換した信号の偏角、つまり位相を求め、アンラップした位相信号φ_m（ｔ）を出力する。

［ステップＳ−２２０９（位相−時間関数生成）］アンラップされた位相φ_m（ｔ）は離散的であるが、式（１１）に示すような、δφ_m間隔のリサンプリング点（時刻）に、位相算出手段２１０５で算出した時間−位相の点があるとは限らないので、式（１１）によってリサンプリング点（時刻）を取得するには不都合である。そこで、位相−時間関数生成手段２１０９は、φ_m（ｔ）からその逆関数ｔ（φ_m）の近似式ｔ’（φ_m）を求めて、どの位相に対しても時間が対応するようにする。非特許文献１では、ｔ（φ_m）を４次関数に近似している。

［ステップＳ−２２１０（リサンプリング点算出）］リサンプリング点算出手段２１１０は、式（１２）のｔ（φ_m）を、その近似式ｔ’（φ_m）に置き換えて、定数δφ_m間隔となるリサンプリング点（時刻）ｔ_n（＝ｔ’（ｎδφ_m＋φ₀）を得て出力する。

次に、サンプル光側の鏡の代わりに観測対象のサンプルを置き、測定を行う段階について説明する。

［ステップＳ−２２１１（干渉波形信号取得２）］観測対象のサンプルから反射されたサンプル光と参照光ミラーからの参照光とから得られる干渉信号を、干渉波形信号取得手段２１０１を用いて取得して出力する。

［ステップＳ−２２１２（リサンプリング）］ステップＳ−２２１１で得た干渉信号を、ステップＳ−２２１０で得たリサンプリング点（時刻）ｔ_nでサンプリング（リサンプリング）する。

上記、離散フーリエ変換手段２１０２、マイナス周波数領域除去手段２１０３、離散フーリエ逆変換手段２１０４、および位相変化曲線算出手段２１０５は、干渉信号ｉ_m（ｔ）から干渉信号位相変化曲線φ_m（ｔ）を取得する手段であり、図２ではこれらの手段をまとめて位相変化曲線取得手段２１１２としている。位相変化曲線取得手段２１１２の動作は、ステップＳ−２２０２〜Ｓ−２２０５で示され、図３ではそれらをまとめてステップＳ−２２１３（位相変化曲線取得）としている。

Yoshiaki Yasuno, et al., "Three-dimensional and high-speed swept-source optical coherence tomography for in vivo investigation of human anterior eye segments" OPTICS EXPRESS, Vol.13, No.26, pp.10652-10664, 2005.

干渉信号ｉ_m（ｔ）の位相φ_m（ｔ）は、以下の式で表される。

従来は、反射面の位置ｚが時間的に変動が無い事を前提としていた。しかし、サンプル光の光路と参照光の光路のいずれかに、光の波長による屈折率の違い（波長分散）がある物体がある挿入されていた場合には、ＳＳ−ＯＣＴ光源が時間的に波長が変化する光源であることから、ｚがｋ（ｔ）によって時間変動する可能性がある。例えば、ｚが次のように表せるとする。

ただし、ｚ₀はｋ＝ｋ₀時のｚであり定数である。ｋ₀は、例えば波長掃引光源の掃引する波数の中心とする。ｒ（ｋ）はｋ＝ｋ₀を基準とした屈折率の増加率を表す関数（屈折率増加関数）である。なお、ｒ（ｋ）は、屈折率の波長分散に係り、波数ｋの時の屈折率をｎ（ｋ）と表すと、ｒ（ｋ）＝ｎ（ｋ）／ｎ（ｋ₀）となる。このとき、式（１３）は次のようになる。

図４に、リファレンス光側には共通の物体（空気等）をおき、サンプル光側には、波長分散の異なる２つの物体aまたは物体ｂを置いた場合の干渉信号の位相の時間推移（位相変化曲線）の例を示す。φ_ma（ｔ）は、物体ａを置いた場合の位相変化曲線であり、φ_mb（ｔ）は、物体ｂを置いた場合の位相変化曲線である。ｔ_anとｔ_bn（図４ではｎ＝０〜１０）は、それぞれの位相変化曲線から得られたリサンプリング点（時刻）である。この図から、波長分散が異なる場合は、リサンプリング点（時刻）が異なるため、適切にリスケーリングできないことが分かる。

リサンプリング点（時刻）を取得した時と同じ波長分散を持つ物体を、観測対象にする場合には、図４に示す２つの位相変化曲線は一致する。従って、リサンプリング点（時刻）も一致するため、適正なリスケーリングができるため問題は生じない。通常のＳＳ−ＯＣＴでは、リサンプリング点（時刻）の取得においては、空気を媒体とする遅延を使用したり、ファイバによる遅延を利用することが多い。しかし、リサンプリング点（時刻）を取得した時の物体（空気、ファイバ）の波長分散と観測対象の屈折率の波長分散とが異なることが多く、図４に示すようなリサンプリング点（時刻）のずれが生じ、適正なリスケーリングができないという問題があった。

その結果、リスケーリングされた干渉信号は、波数リニア光源から得られた干渉信号からずれるため、断層画像の点広がり関数（ＰＳＦ）の幅（半値全幅）が広がることから、ＳＳ−ＯＣＴ画像（断層画像）の精細度が低くなるという問題があった。また、ＰＳＦのピーク値が低くなることによって、ＳＳ−ＯＣＴ画像（断層画像）のＳＮＲが悪化するという問題もあった。

本発明の目的は、観測対象についての干渉信号に対して適正なリサンプリング（リスケーリング）を行うことにより、ＳＳ−ＯＣＴ画像（断層画像）の精細度を高くし、ＳＮＲを向上させるための信号処理装置およびリスケーリング方法を提供することにある。

本発明は、このような目的を達成するために、一実施態様は、波数が時間的に一方向に連続的に変化する光を出力する波長掃引光源と、該波長掃引光源から出力された光を入力し、サンプルアームと参照アームと間の光路長差に応じた干渉信号を出力する干渉計とを含み、該干渉計から出力された干渉信号を周波数解析して断層信号を出力する光断層撮像装置においてリスケーリング処理を行うための信号処理装置であって、サンプルａを用いて単一の光路長差に係る位相変化曲線φ_m(ｔ)を出力する位相変化曲線取得手段と、前記サンプルａの屈折率分散データｎ_a（ｋ）、観測対象の屈折率分散データｎ_b（ｋ）、および前記波長掃引光源から出力された光の時間−波数対応データｋ（ｔ）から、波数ｋ₀における屈折率を基準としたときの、サンプルａの時間−屈折率変動データｒ_a(ｔ)（＝ｎ_a（ｋ）／ｎ_a（ｋ₀））に対する観測対象の時間−屈折率変動データｒ_b(ｔ)（＝ｎ_b（ｋ）／ｎ_b（ｋ₀））の比ｒ_b(ｔ)／ｒ_a(ｔ)を時間−補正係数対応データｒ_a→b(ｔ)として出力する波長分散補正係数算出手段と、前記時間−補正係数対応データｒ_a→b(ｔ)と前記位相変化曲線φ_m(ｔ)との積を補正された位相変化曲線φ_m’(ｔ)として出力する位相変化曲線補正手段と、前記補正された位相変化曲線φ_m’(ｔ)からリサンプリング点ｔ_nを算出するリサンプリング点算出手段とを備えたことを特徴とする。

前記信号処理装置は、前記時間−波数対応データｋ（ｔ）を出力する時間−波数対応データ取得手段をさらに備えることができる。

本発明によれば、観測対象の屈折率の波長分散に応じたリサンプリング点（時刻）を使って、観測対象についての干渉信号をリサンプリング（リスケーリング）するので、ＰＳＦの幅（半値全幅）が狭くなり、ＳＳ−ＯＣＴ画像（断層画像）の精細度を高くすることができる。また、ＰＳＦのピークが高くなることにより、ＳＳ−ＯＣＴ画像のＳＮＲを向上させることができる。

従来のＳＳ−ＯＣＴ装置の基本構成の一例を示す図である。 従来のリスケーリング装置の機能ブロックを示す図である。 従来のリスケーリング処理を示すフローチャートである。 干渉信号の位相の時間推移（位相変化曲線）の例を示す図である。 本発明の実施例１にかかるリスケーリング装置の機能ブロックを示す図である。 本発明の実施例１にかかるリスケーリング処理を示すフローチャートである。 本発明の実施例１にかかるＳＳ−ＯＣＴ装置の基本構成を示す図である。 実施例２にかかる時間−波数対応データ取得手段の機能ブロックの第１例を示す図である。 実施例２にかかる時間−波数対応データ取得手段における処理を示すフローチャートである。 時間−波数対応データ取得手段の機能ブロックの第２例を示す図である。 リットマン型波長掃引光源の光学系の基本構成を示す図である。 ＫＴＮ光偏向器の構成を示す図である。 実施例３にかかる時間−波数対応データ取得手段の機能ブロックを示す図である。 実施例３にかかる時間−波数対応データ取得手段における処理を示すフローチャートである。 疎な時間−波数対応データの外挿または内挿の例を示す図である。 時間−波数対応データ取得手段におけるマッハツェンダ干渉計の構成を示す図である。 実施例４にかかる時間−波数対応データ取得手段における機能ブロックを示す図である。 実施例４にかかる時間−波数対応データ取得手段における処理を示すフローチャートである。 位相変化曲線補正手段の機能ブロックを示す図である。 位相変化曲線補正手段における処理を示すフローチャートである。 時間−波数対応データ生成手段の機能ブロックを示す図である。 時間−波数対応データ生成手段における処理を示すフローチャートである。 屈折率分散データの温度補正手段の機能ブロックを示す図である。 屈折率分散データの温度補正手段における処理を示すフローチャートである。

以下、図面を参照しながら本発明の実施形態について詳細に説明する。
本実施形態においては、観測対象となるサンプルの波長分散に合わせて位相変化曲線に波数ｋ毎の係数をかけることにより、観測対象の波長分散に合った位相変化曲線を作成することによって、観測対象に合ったリサンプリング点（時刻）を得る。リサンプリング点（時刻）を求めるとき、サンプルａをサンプルとして用いることとし、光の波数ｋ＝ｋ₀時の屈折率ｎ_a（ｋ₀）を基準としたとき、波数ｋに対する屈折率の増加率を表す屈折率増加関数をｒ_a（ｋ）（＝ｎ_a（ｋ）／ｎ_a（ｋ₀））とする。ｋ₀は、例えば波長掃引光源の掃引する波数の中心としてもよい。このとき、サンプルａを用いた時の干渉信号の位相変化曲線φ_ma（ｔ）は下記のようになる。

ただし、ｚ_a（ｋ）は、サンプルａのｋに対する光路長であり、ｚ_a0は、ｋ＝ｋ₀時のｚ_a（ｋ）である。この式から、波長分散による光路長の変化を補償した干渉信号の位相は、下記のように示される。

つまり、サンプルａを使って観測された干渉信号の位相φ_ma（ｔ）をｒ_a（ｋ（ｔ））で割れば、ｋ＝ｋ₀時の光路長ｚ_a0に対応した位相φ_ma0（ｔ）が得られる。

ｋ（ｔ）が時間ｔの関数であるので、ｒ_a（ｋ）は時間ｔの関数となるが、このためには、波長掃引光源からの出射光の波数の時間変化の関数ｋ（ｔ）を知っておく必要がある。このｋ（ｔ）を知る方法については、実施例２に記す。

ｋ（ｔ）が分かるとｒ_a（ｋ）を時間の関数とすることができる。ここでは説明を簡単にするため、以下の説明ではｒ_a（ｋ（ｔ））をｒ_a（ｔ）と置き換えることとする。このｒ_a（ｔ）を使って式（１−２）を改めて書き直すと下記となる。

式（１−３）に倣うと、観測対象としてサンプルｂを使ったときの、ｋ＝ｋ₀時の光路長ｚ_b0に対応した干渉信号の位相は、以下のように示されることになる。

ただし、φ_mb（ｔ）はサンプルｂを使った時に観測される干渉信号であり、ｒ_b（ｔ）は、ｋ＝ｋ₀を基準としたときの時間による屈折率変化を表す関数ｒ_b（ｋ）（＝ｎ_b（ｋ）／ｎ_b（ｋ₀））であり、ｎ_b（ｋ）はサンプルｂの屈折率、ｒ_b（ｔ）＝ｒ_b（ｋ（ｔ））である。また、ｚ_b0はｋ＝ｋ₀時のｚ_b（ｋ）である。

式（１−３）、式（１−４）から、φ_ma（ｔ）とφ_mb（ｔ）の関係は下記となる。

つまり、φ_ma（ｔ）に（ｚ_b0／ｚ_a0）とｒ_a→b（ｔ）（＝ｒ_b（ｔ）／ｒ_a（ｔ））をかけるとφ_mb（ｔ）となる。ここで、ｒ_a→b（ｔ）＝（ｒ_b（ｔ）／ｒ_a（ｔ））である。

φ_mb（ｔ）は、（ｚ_b0／ｚ_a0）・ｒ_a→b（ｔ）・φ_ma（ｔ）と等しいので、式（１１）に示すように、定数δφ_mbを決めて、

となるｔ_a,nとを同じにすること、つまりｔ_b,n＝ｔ_a,nとできる。式（１−８）右辺において（ｚ_b0／ｚ_a0）で括ることができる理由は、（ｚ_b0／ｚ_a0）が時間に対して不変であるからである。ところで、

とすることができる。つまり、位相間隔δφ_mbを（ｚ_a0／ｚ_b0）倍すれば、ｔ_b,n＝ｔ_a,nとなるリサンプリング点（時刻）が得られる。このことは、サンプルａの位相変化曲線にｒ_a→b（ｔ）をかけて補正した位相変化曲線を用いて、式（１１）または式（１２）からサンプルｂ用のリサンプリング点（時刻）を取得できることを意味している。

ところで、（ｚ_b0／ｚ_a0）は時間に対して不変であるので、式（１−６）からφ_mb（ｔ）とｒ_a→b（ｔ）・φ_ma（ｔ）には下記のような関係があることが分かる。

このようにｒ_a→b（ｔ）φ_ma（ｔ）は、φ_mb（ｔ）と比例関係にあることが分かる。このように比例関係にある波数変化曲線では、波数リニア光源を用いた干渉信号を得るためのリサンプリング点（時刻）は、両者とも同一にすることが可能であるという言い方もできる。

以上の説明では、サンプルａを使って得られる位相変化曲線を補正することによって得られるリサンプリング点（時刻）を、サンプルｂ用のリサンプリング点（時刻）と同一にできる、という言い方をしたが、リサンプリングの性質上、同一であることは重要ではない。例えば、式（１−９）では、δφ_ma=（ｚ_a0／ｚ_b0）δφ_mbとしたが、補正したサンプルａの位相変化曲線ｒ_a→b（ｔ）・φ_ma（ｔ）は、適宜決めたδφ_maでリサンプリング点（時刻）を得ても、サンプルｂに適したリサンプリング点（時刻）を得ることはできる。

以上の原理を使って、実際にリスケーリングを行うための装置の構成と処理のフローを説明する。リスケーリング処理の動作は、従来と同様に２段階に分けられ、最初はリサンプリング点（時刻）を取得する段階であり、次に、観測対象のサンプルから得られた干渉信号をリスケーリングする段階である。リサンプリング点（時刻）を取得する段階では、従来の位相変化曲線を取得するフローと位相変化曲線の補正係数を算出するフローとを含む段階と、その補正係数を使って位相変化曲線を補正してリサンプリング点（時刻）を取得する段階とがある。

図５に、本発明の実施例１にかかるリスケーリング装置の機能ブロックを示し、図６に、本発明の実施例１にかかるリスケーリング処理のフローチャートを示す。最初に、干渉信号をリスケールするためのリサンプリング点（時刻）を取得する段階について説明する。この段階における位相変化曲線を算出するフローは、従来と同様であり、下記の通りとなる。

［ステップＳ−２０１（干渉波形信号取得１）］サンプル光側のｚの位置に鏡などの反射面を置き、干渉計２００２で生成した干渉信号を、干渉波形信号取得手段１０１が取得して出力する。受光素子２００７がアナログ出力であり、そのアナログ出力の信号を干渉波形信号取得手段１０１が受ける場合は、Ａ／Ｄ変換を行い、ディジタル信号として出力する。

干渉波形信号取得手段１０１が干渉信号を取得するタイミングは、波長掃引光源が短波長から長波長へ、または、逆に長波長から短波長へ掃引を開始するタイミングとする。例えば、波長掃引光源から掃引開始時刻を示すトリガ信号を、干渉波形信号取得手段１０１が得るようにしても良い。または、トリガ信号を出力するためのトリガ信号出力手段を別途用意し、波長掃引光源にトリガ信号を入力して、トリガ信号のタイミングで（あるいは同期して）波長掃引光源が掃引を開始しても良い。

［ステップＳ−２０２（離散フーリエ変換）］離散フーリエ変換手段１０２は干渉波形取得手段から得た干渉信号を離散フーリエ変換して出力する。

［ステップＳ−２０３（マイナス周波数領域除去）］マイナス周波数領域除去手段１０３は、離散フーリエ変換した干渉信号のうち、負の周波数を除去し、除去した信号を出力する。

［ステップＳ−２０４（離散フーリエ逆変換）］離散フーリエ逆変換手段１０４は、負の周波数を除去した信号を逆フーリエ変換して出力する。

［ステップＳ−２０５（位相変化曲線算出）］位相変化曲線算出手段１０５は、逆フーリエ変換した信号の偏角、つまり位相を求め、アンラップした位相信号φ_m（ｔ）を出力する。

次に、位相変化曲線の補正係数を算出するフローについて説明する。
［ステップＳ−２０６（時間−波数対応データ取得）］時間−波数対応データ取得手段１０６は、波長掃引光源から出力される光の波数の時間変化データ（時間−波数対応データ）を取得し出力する。具体的な方法は実施例２に示す。

［ステップＳ−２０７（波長分散補正係数算出）］波長分散補正係数算出手段１０７は、ステップＳ−２０６において時間−波数対応データ取得手段１０６から得た時間−波数対応データｋ（ｔ）と、屈折率分散データｎ_a（ｋ）、ｎ_b（ｋ）と、基準波数ｋ₀における屈折率を基準としたときの位相変化曲線の補正係数ｒ_a→b（ｔ）を計算して出力する。ここで、ｎ_a（ｋ）は、リサンプリング点（時刻）を得るために使用するサンプルａの屈折率の波数分散（波長分散）であり、ｎ_b（ｋ）は、観測対象のサンプルｂの屈折率の波数分散（波長分散）である。

計算の一例としては、最初に、波数ｋ毎のｒ_a→b（ｋ）＝（ｎ_b（ｋ）／ｎ_b（ｋ₀））／（ｎ_a（ｋ₀）／ｎ_a（ｋ））のデータを作り、次に、ｒ_a→b（ｋ）のデータとｋ（ｔ）のデータを使って、波数ｋを時間ｔに置き換えて（つまりｒ_a→b（ｔ）＝ｒ（ｋ（ｔ））として）、波数ｋ毎の補正係数ｒ_a→b（ｋ）のデータを、時間ｔ毎の補正係数ｒ_a→b（ｔ）のデータにする。

ｋ（ｔ）データとｒ_a→b（ｋ）データとは、両方とも離散的であるので、あるｔに対するｋ（ｔ）が、ｒ_a→b（ｋ）データの中にあるとは限らない。その場合は、ｒ_a→b（ｋ）データの近似関数（ｋの多項式関数等）を算出するか、または、ｒ_a→b（ｋ）データを補間して（ｋの多項式による補間等）、波数ｋに対するｒ_a→b（ｋ）データの近似値を算出する。

基準波数ｋ₀を波長分散補正係数算出手段１０７に明示的に入力するのではなく、掃引波数の最小値ｋ_minと最大値ｋ_maxを入力し、ｋ₀＝（ｋ_max＋ｋ_min）／２としてもよい。また、波長から換算し、最小波長λ_minと最大波長λ_maxを入力し、ｋ₀＝２π（１／λ_min＋１／λ_max）／２としても良い。

上記ｒ_a→b（ｋ）のデータは、ステップＳ−２０７を実行する前に作っておいても良い。この場合は、その計算時間が省けるため、高速化ができると共に、ステップＳ−２０１〜Ｓ−２０５と並列に動作することが無い。従って、計算リソースが省け、本実施形態のリスケーリング装置の小型化、低価格化を実現することができる。

次に、補正係数を使って位相変化曲線を補正した後、リサンプリング点（時刻）を取得する段階について説明する。

［ステップＳ−２０８（位相変化曲線補正）］位相変化曲線補正手段１０８は、ステップＳ−２０７において波長分散補正係数算出手段から出力された時間−補正係数対応データｒ_a→b（ｔ）を用いて、ステップＳ−２０５において位相変化曲線算出手段１０５から出力された位相変化曲線φ_m（ｔ）を補正し、観測対象であるサンプルｂに合った位相変化曲線φ_m’（ｔ）にする。理想的には式（１−１１）に示すφ_mb（ｔ）がφ_m’（ｔ）と一致する。

具体的な動作としては、時間ｔ毎にφ_m（ｔ）ｒ_a→b（ｔ）（＝φ_m’（ｔ））を計算し、ｔに対する位相φ_m’（ｔ）のデータを作成し、出力する。

［ステップＳ−２０９（位相−時間関数生成）］アンラップされた位相φ_m（ｔ）は、離散的であるが、式（１１）に示すような、δφ_m間隔のリサンプリング点（時刻）に、位相変化曲線算出手段１０５で算出した時間−位相の点があるとは限らないので、式（１１）によってリサンプリング点を取得するには不都合である。そこで、位相−時間関数生成手段１０９は、φ_m（ｔ）からその逆関数ｔ（φ_m）の近似式ｔ’（φ_m）を求めて、どの位相に対しても時間が対応するようにする。非特許文献１では、ｔ（φ_m）を４次関数に近似している。

［ステップＳ−２１０（リサンプリング点算出）］リサンプリング点算出手段１１０は、式（１２）のｔ（φ_m）を、その近似式ｔ’（φ_m）に置き換えて、定数δφ_m間隔となるリサンプリング点（時刻）ｔ_n（＝ｔ’（ｎδφ_m＋φ₀）を得て出力する。

最後に、サンプル光側の鏡の代わりに観測対象のサンプルを置き、測定を行う段階について説明する。この段階は、従来と同様であり、下記のように処理を行う。

［ステップＳ−２１１（干渉波形信号取得２）］観測対象のサンプルから反射されたサンプル光と参照光ミラーからの参照光とから得られる干渉信号を、干渉波形信号取得手段１０１を用いて取得して出力する。

［ステップＳ−２１２（リサンプリング）］ステップＳ−２１１で得た干渉信号を、ステップＳ−２１０で得たリサンプリング点（時刻）ｔ_nでサンプリング（リサンプリング）する。

以上の説明では、ステップＳ−２１１で観測対象のサンプルについての干渉波形を得る前に、サンプル光側の鏡の代わりに観測対象のサンプルを置くこととした。干渉計を２つ用いて、一方はリサンプリング点（時刻）を得るための干渉信号を取得し、他方は観測対象についての干渉信号を取得するための干渉計としても良い。２つの干渉計を使う場合、波長掃引光源の１掃引ごとに２つの干渉計から同時に干渉信号を取得し、１掃引ごとにリサンプリング点（時刻）のデータを算出し、同じ掃引で取得した観測対象の干渉波形をリサンプリングしても良い。この方法を、セルフクロッキング法（self clocking method)と呼ぶ。

２つの干渉計を使用する利点は、リサンプリング点（時刻）を取得するための鏡と観測対象のサンプルを置き換える必要が無いことである。また、セルフクロッキング法を適用する場合は、波長掃引光源の時間変動等で、波長掃引光源からの出射光の波長（または波数）の時間変化が１掃引ごとに異なる場合であっても、１掃引ごとにリサンプリング点（時刻）を取得するので、常に最適なリサンプリング点（時刻）が得られる利点がある。

以上の説明では、ステップＳ−２０１〜Ｓ−２０５とステップＳ−２０６〜Ｓ−２０７をシーケンシャルに動作するように記した。ステップＳ−２０１〜Ｓ−２０５とステップＳ−２０６〜Ｓ−２０７とは、図６のフローチャートに示す通り、並列して動作しても良い。また、先にステップＳ−２０６〜Ｓ−２０７を動作させてからステップＳ−２０１〜Ｓ−２０５を動作させても良い。並列して動作させる場合は、動作時間を短縮できる。また、シーケンシャルに動作させる場合は、計算リソース（ハードウェア）を共用できるため、本実施形態のリスケーリング装置の小型化、低価格化を実現することができる。

図７に、本発明の実施例１にかかるＳＳ−ＯＣＴ装置の基本構成を示す図である。干渉計は、光ファイバ２４１５、２４１７、２４１９とサーキュレータ２４１６、ファイバコリメータ２４１８で構成され、２面が平行なサンプル２４０３を用いている。サンプル２４０３の入射側の面は基準面となり、対向する面は反射面となる。入射面と反射面とは平行になっている。サンプル２４０３の例としてはエタロンなどが挙げられる。サンプル２４０３の屈折率がサンプル２４０３の周りと異なれば、基準面２４１１とサンプル内の反射面２４０８での光の反射はフレネル反射となる。

波数が時間的に一方向に連続的に変化する光を出力する波長掃引光源２４０１から出力された出射光２４０９は、光ファイバ２４１５を使って出力され、サーキュレータ２４１６、光ファイバ２４１７、ファイバコリメータ２４１８を介してサンプル２４０３の入射面（基準面）に垂直入射する。サンプル２４０３に入射した光は、反射面２４０８で反射される。サンプル２４０３を透過する光２４１０は、サンプル光であり、入射面からは、サンプル２４０３内から出射された光と入射面で反射された光（参照光２４１２）とが混合され干渉光２４１３となる。すなわち、サンプルアームと参照アームと間の光路長差に応じた干渉光が得られる。

干渉光２４１３は、ファイバコリメータ２４１８、光ファイバ２４１７、サーキュレータ２４１６を介して光ファイバ２４１９に入り、受光素子２４０７で光電変換される。受光素子２４０７からは、干渉光が光電変換された干渉信号２４１４が出力され、信号処理装置２４０４に入って、リスケーリング処理がなされ、ＰＳＦが生成される。

図１に示した干渉計の代わりに図７に示した干渉計を用いて、観測対象をサンプル２４０３の位置に配置しても良い。この場合は、観測対象の厚みを測定するＳＳ−ＯＣＴ装置として働く。観測対象が平行平板であり、その厚みを測定するのであれば、図１の構成よりも図７の構成とした方が、ＳＳ−ＯＣＴ装置の構成が単純であり、安価に作製できるという利点がある。

上記、離散フーリエ変換手段１０２、マイナス周波数領域除去手段１０３、離散フーリエ逆変換手段１０４、および位相変化曲線算出手段１０５は、干渉信号ｉ_m（ｔ）から干渉信号位相変化曲線φ_m（ｔ）を取得する手段であり、図５ではこれらの手段をまとめて位相変化曲線取得手段１１２としている。位相変化曲線取得手段１１２の動作は、ステップＳ−２０２〜Ｓ−２０５で示され、図６ではそれらをまとめてステップＳ−２１３（位相変化曲線取得）としている。

以上の説明では、Ｓ−２１１〜Ｓ−２１２（観測対象の干渉波形を取得して、リサンプリングする動作）を１回行う前に、Ｓ−２０１〜Ｓ−２１０を１回行う例を示した。光源の波長の時間変化ｋ（ｔ）が各掃引において同一であるか、または、各掃引のｋ（ｔ）の変動が無視できる場合は、次の２段階の動作としても良い。
（ａ）Ｓ−２０１〜Ｓ−２１０を行い、その後はｋ（ｔ）の掃引ごとの変動が無視できる掃引回数だけ、光源の時間変動Ｓ−２１１〜Ｓ−２１２のみを繰り返し行う。
（ｂ）変動が無視できない掃引回数に達した場合に、再度Ｓ−２０１〜Ｓ−２１０を行い、位相変化曲線と時間−補正係数対応データとを初期化する。

このような構成にすると、Ｓ−２１１〜Ｓ−２１２（観測対象の干渉波形を取得して、リサンプリングする動作）を共通のリサンプリング点（時刻）を用いて行うことが可能となるため、リサンプリング点を生成する動作を少なくでき、リサンプリング点を生成するハードウェアリソースを削減できる利点がある。

（時間−波数対応データ取得手段）
実施例１における時間−波数対応データｋ（ｔ）の取得手段の具体的な構成について述べる。ここでは、光の狭帯域バンドパスフィルタの中心波長と時間の関係を取得する方法について述べる。図８に、時間−波数対応データ取得手段の機能ブロックの第１例を示し、図９にフローチャートを示す。

［ステップＳ−４０１（測定波数範囲および初期波数設定）］波数指示手段３０４に波長掃引光源の光の波数の測定範囲ｋ_min〜ｋ_maxと、波数のステップδｋを設定する。ここで、ｋ_min＜ｋ_max、δｋ＞０である。また、測定のスタート時の波数ｋを設定する。例えば、ｋ＝ｋ_minとする。

［ステップＳ−４０２（終了判定）］コントローラ３０８は、測定波数ｋがｋ_min≦ｋ≦ｋ_maxか否かを判定する。ｋ_min≦ｋ≦ｋ_maxが真（図９中ではＹｅｓ）の場合、ステップＳ−４０３へ進み、そうでない（図９中ではＮｏ）場合は、終了する。

［ステップＳ−４０３（トリガ信号発信）］波長掃引光源３０１は、掃引開始時にトリガ信号を出力し、光を出力する。

［ステップＳ−４０４（波数指示）］波数指示手段３０４は、波長掃引光源３０１からトリガ信号を受け、時間を測る対象となる光源の波数を出力する。ここでは波数としたが、光の波長、光の周波数としても良い。

［ステップＳ−４０５（フィルタリング）］狭帯域バンドパスフィルタ３０２は、波数指示手段３０４から波数を受け、その波数を中心周波数とする透過帯域として、波長掃引光源３０１から得た光を透過させ、その透過光を出力する。

［ステップＳ−４０６（光検知）］狭帯域バンドパスフィルタ３０２から出力されたフィルタリングされた光は、光検知器３０３で受光され、光の強度信号を出力する。

［ステップＳ−４０７（時間計測）］時間計測手段３０６は、光検知器３０３から出力された光強度信号と、波長掃引光源３０１からトリガ信号とを受け、トリガ信号と光強度信号との時間間隔を得る。光強度信号やトリガ信号の受信時刻を得る方法としては、光強度信号のピークを検出して、その検出時刻を光強度の信号の受信時刻とする方法、または、それぞれあるスレッショルドをあらかじめ定めておいて、そのスレッショルドを超える時刻をそれぞれの受信時刻とする方法がある。

［ステップＳ−４０８（時間−波数対応データ記録）］波数指示手段３０４が狭帯域バンドパスフィルタ３０２に設定させた波数と、波長掃引光源３０１のトリガ信号から光検知器３０３が光を検知するまでの時間とを対応させた時間−波数対応データｋ（ｔ）を記録手段３０５に記録する。

［ステップＳ−４０９（波数更新）］波数ｋを更新する。例えば、ｋをｋ＋δｋと置き換える。そして、ステップＳ−４０２に戻る。

上記の動作手順では、測定スタート時の光の波数ｋをｋ_minとしたが、ｋ_maxとして波数更新ステップＳ−４０９において、ｋをｋ−δｋとしても良い。上記の動作手順では、光の波数ｋをパラメータとしたが、光の波長λや周波数ｆをパラメータとしても良い。その場合は、ｋ＝２π／λ、ｋ＝２πｆ／ｃとして波数ｋに置き換え、上記の計算を行う。

図１０に、時間−波数対応データ取得手段の機能ブロックの第２例を示す。第１例では、波長掃引光源３０１がトリガ信号を出力して、波数指示手段３０４と時間計測手段３０６が受ける構成であった。第２例では、トリガ信号出力手段３０７を付加し、波長掃引光源３０１と波数指示手段３０４と時間計測手段３０６に波長掃引の開始時を指定するトリガ信号をトリガ信号出力手段３０７が出すようにしても良い。その場合は、動作手順は下記のようになる。

［ステップＳ−４０１（測定波数範囲および初期波数設定）］波数指示手段３０４に波長掃引光源の光の波数の測定範囲ｋ_min〜ｋ_maxと、波数のステップδｋを設定する。ここで、ｋ_min＜ｋ_max、δｋ＞０である。また、測定のスタート時の波数ｋを設定する。例えば、ｋ＝ｋ_minとする。

［ステップＳ−４０２（終了判定）］コントローラ３０８は、測定波数ｋがｋ_min≦ｋ≦ｋ_maxか否かを判定する。ｋ_min≦ｋ≦ｋ_maxが真（図９中ではＹｅｓ）の場合、ステップＳ−４０３へ進み、そうでない（図９中ではＮｏ）場合終了する。

［ステップＳ−４０３（トリガ信号発信）］トリガ信号出力手段３０７から出力されたトリガ信号を波長掃引光源３０１が受け取り、波長掃引光源３０１は掃引を開始し、光を出力する。

［ステップＳ−４０４（波数指示）］波数指示手段３０４は、トリガ信号出力手段３０７からトリガ信号を受け、時間を測る対象となる光源の波数を出力する。ここでは波数としたが、光の波長、光の周波数としても良い。

［ステップＳ−４０５（フィルタリング）］狭帯域バンドパスフィルタ３０２は、波数指示手段３０４から波数を受け、その波数を中心周波数とするバンドパスとして、波長掃引光源３０１から得た光を透過させ、その透過光を出力する。

［ステップＳ−４０６（光検知）］狭帯域バンドパスフィルタ３０２から出力されたフィルタリングされた光は、光検知器３０３で受光され、光の強度信号を出力する。

［ステップＳ−４０７（時間計測）］時間計測手段３０６は光検知器３０３から出力された光強度信号と、トリガ信号出力手段３０７からトリガ信号とを受け、トリガ信号と光強度信号の時間間隔を得る。光強度信号やトリガ信号の受信時刻を得る方法としては、光強度信号のピークを検出してその検出時刻を光強度の信号の受信時刻とする方法、または、それぞれあるスレッショルドをあらかじめ定めておいて、そのスレッショルドを超える時刻をそれぞれの受信時刻とする方法がある。

［ステップＳ−４０８（時間−波数対応データ記録）］波数指示手段３０４が狭帯域バンドパスフィルタ３０２に設定させた波数と、トリガ信号出力手段３０７のトリガ信号から光検知器３０３が光を検知するまでの時間とを対応させた時間−波数対応データｋ（ｔ）を記録手段３０５に記録する。

［ステップＳ−４０９（波数更新）］波数ｋを更新する。例えば、ｋをｋ＋δｋと置き換える。そして、ステップＳ−４０２に戻る。

狭帯域バンドパスフィルタ３０２の帯域は、波長掃引光源３０１の掃引幅の１％以下であることが望ましい。例えば、波長掃引幅が１００ｎｍである場合は、狭帯域バンドパスフィルタのバンド幅は中心波長±０．５ｎｍ以下であることが望ましい。

図８、１０において、狭帯域バンドパスフィルタ３０２、光検知器３０３、波数指示手段３０４、記録手段３０５、時間計測手段３０６、コントローラ３０８をまとめたものを、時間−波数対応データ実測手段３０９と記載した。これについては、後述の別の時間−波数対応データ取得手段の説明で使用する。

（時間−波数対応データ補間生成手段）
上記の時間−波数対応データ取得手段においては、波数のステップδｋを小さくし、波数を密に取得する構成であった。δｋを大きくして波数を疎に取得し、時間−波数対応データ取得後、その疎な時間−波数対応データを、時間−波数対応データ補間生成手段によって補間して、密な時間−波数対応データを作成することが考えられる。

時間−波数対応データ補間生成手段で計算する例としては、ｔの多項式でｋを表し、その多項式から必要な数の時刻を入力して波数を得るということが考えられる。ｎ次関数で表す場合は、ｋ（ｔ）＝ａ_nｔⁿ＋ａ_n-1ｔ^n-1＋…＋ａ₀のような多項式が考えられる。ここで、ａ_n、ａ_n-1、…ａ₀は未知パラメータであり、すでに取得した時間−波数対応データ取得を使って、最小二乗法によりａ_n、ａ_n-1、…ａ₀を推定する。このように生成した多項式のｔに必要なだけの時刻を入力し、それに対する波数を得る。この方法で時間−波数対応データを作成する場合は、時間−波数対応データの測定数が少ないため、短時間で時間−波数対応データを作成することが期待できる。

また、時間−波数対応データ補間生成手段で計算する別の例としては、光源の構成が分かっている場合は、その構成から得られる時間−波数対応データの式を、上記の疎な時間−波数対応データに当てはめて、疎な時間−波数対応データを補間して、密な時間−波数対応データを作成することが考えられる。

図１１に、リットマン型波長掃引光源の光学系の基本構成を示す。半導体光増幅器（ＳＯＡ）８０４から出射した光は、コリメータレンズ８０６でコリメートされ、電源８１０からの印加電圧によって光偏向器８０７により偏向される。偏向光は、回折格子８１１で回折された後、ミラー８１２に達する。ミラー８１２に垂直入射した光のみが選択的に逆の光路をたどりＳＯＡ８０４に戻る。戻り光は、ＳＯＡ８０４内部で増幅された後、ＳＯＡ反射面８０５で一部反射し、再度ＳＯＡ８０４内部で増幅された後にコリメータレンズ８０６側から出射する。このようにしてＳＯＡ反射面８０５とミラー８１２との間で光が繰り返し往復することによって、回折格子８１１とミラー８１２とで構成される光のバンドパスフィルタによる発振線幅の狭窄化が生じてレーザ発振する。ＳＯＡ反射面８０５はハーフミラーとなっており、ＳＯＡ８０４からの光が一部通過し、集光レンズ８０３によってアイソレータ８０２にカップリングし、光ファイバ８０１を通してレーザ光が出力される。

このような構成の光源では、レーザ光の波数ｋは以下のように表される。

ここで、ｍは回折格子８１１の回折光次数、Ｎは回折格子８１１の刻線数、Ψは光偏向器８０７の光偏向角、λ₀はΨ＝０時の発振波長であり、下記の式で表される。

ここで、φは光偏向器の出射光軸（Ψ＝０時の光の出射光軸）とミラーの法線との成す角である。Ψの時間変動を表す関数が既知であり、その関数をΨ（ｔ）と表すと、式（２−１）はｋの時間関数となり、下記のようになる。

例えば、光偏向器８０７として、ＫＴＮ光偏向器を使う場合を考える。

図１２に、ＫＴＮ光偏向器の構成を示す。直方体の形状をしているＫＴＮ（タンタル酸ニオブ酸カリウム）単結晶９０１の対向する面に電極９０２、９０３がオーミック接合されている。電極間の結晶厚みをｄ、ＫＴＮ単結晶の光軸９０４方向の結晶長さをＬで表している。ＫＴＮ偏向器は、電極９０２、９０３間に直流電圧を印加してＫＴＮ結晶９０１内部に電子がトラップした状態で、使用する。電子がトラップされると、ＫＴＮ結晶９０１は、ガウスの法則とカー効果に従って凸レンズとなる。この状態で電極９０２、９０３にさらにＡＣ電圧を印加すると、ＡＣ電圧強度に従って、凸レンズが光軸９０４とは垂直な方向に移動するため、ＫＴＮ結晶９０１に入射した光は、凸レンズの移動とともに偏向されて、出射光９０５となってＫＴＮ結晶９０１から出射される。

ＫＴＮ偏向器の光の偏向角Ψは下記のように表される。

ただし、ｎ₀はＫＴＮ結晶９０１に電子をトラップさせないときの屈折率、ｇ₁₁は電気光学定数、εはＡＣ電圧の周波数以下でのＫＴＮ結晶９０１の誘電率、Ｖ（ｔ）は電極９０２、９０３間の電位差、ρはＫＴＮ結晶９０１内の電荷密度である。式（２−４）を式（２−３）に代入すると、ｋは以下のようになる。

ただし、ι₀＝−（ｎ₀）²√ｇ₁₁ε・ｓｉｎ（ｎ₀√ｇ₁₁ρＬ）／ｄである。もし、Ｖ（ｔ）＝Ｖ₀ｓｉｎ（２πｆｔ）であり、Ｖ₀は電圧振幅、ｆは電圧変動の周波数とすると、ｋは以下のようになる。

式（２−６）の中で、λ₀、ｍ、Ｎ、α、ι₀Ｖ₀、ｆの内、分かっていないものについては、実測した疎の時間−波数対応データを式（２−６）で当てはめることにより、未知パラメータを推定できる。例えば、Levenberg-Marquardt法を使って推定することが考えられる。得られた時間−波数の関数から計算する方法としては、上記で時間−波数の関数を多項式で表した場合と同様に、必要な数の時刻を入力して波数を得るということが考えられる。

このような方法で時間−波数対応データを作成する場合は、時間−波数対応データの測定数が少ないため、短時間で時間−波数対応データを作成することが可能となる。また、光源の構成に整合した式を当てはめるため、光源の構成とは関係のない多項式を使用するよりは、高精度な時間−波数対応データを生成できる。

上記は、ｋ（ｔ）が正確に予測できないことを前提に、時間−波数対応データを実測するか、疎の時間−波数対応データを実測して補間するという構成について述べた。もし、時間−波数対応データが計算で正確に算出できるのであればその計算式を使用して、時間−波数対応データを作成しても良い。例えば、図１１に示した波長掃引光源を用いた場合、式（２−３）、（２−５）、（２−６）でｋ（ｔ）を算出、つまり、時間−波数対応データを算出できる。このように計算によって時間−波数対応データを作成できれば、時間−波数対応データの作成時間は、極めて短縮できる利点がある。

上記の式（２−１）〜（２−６）は、図１１に示したリットマン型波長掃引光源用の式であるが、φ＝０とすれば、リトロー型波長掃引光源でも同様に適用することができる。

実施例２に引き続き、時間−波数対応データｋ（ｔ）の取得手段の具体的な方法について述べる。実施例２では、時間−波数対応データを直接的に求める方法について述べたが、この方法では、多くの時間−波数対応データを取得する必要がある。実施例３では、直接的に求める時間−波数対応データを削減し、処理の軽減や時間短縮を行う。最初に、光の狭帯域バンドパスフィルタの中心波長と時間との関係（時間−波数対応データ）をいくつか取得する。次に、それら時間と波数の対応データの間を、図５に示す位相変化曲線取得手段１１２と同様の方法で取得した位相変化曲線を使って外挿または内挿をする。そのために、波長分散が無い（真空）場合、または、観測対象に対して波長分散が無視できる光路長を使い、干渉信号を得る。

図１３に、実施例３にかかる時間−波数対応データ取得手段の機能ブロックを示し、図１４にフローチャートを示す。波長掃引光源１００１の出力は、光ファイバ１０１５を通ってビームスプリッタ１０１６で分岐され、一方は光ファイバ１０１７を通して干渉計１００２に入り、他方は光ファイバ１０１９を通して時間−波数対応データ実測手段１００９に入る。干渉計１００２では、光ファイバ１０１７を通して入力された光を使い、波長分散が無い光路長差、または、波長分散が観測対象と比較して無視できる光路長差の光を干渉させ、受光素子１００７で干渉光を受光して干渉信号１０１４を得る。この時の光路長差は予め設定された既知のものであり、その光路長差をｚとする。

位相変化曲線実測手段１０１０では、波長掃引光源１００１から出力される掃引開始を示すトリガ信号に同期して、干渉波形信号取得手段１０１１が、干渉信号ｉ_m（ｔ）を生成し、位相変化曲線取得手段１０１２が、干渉信号ｉ_m（ｔ）から位相変化曲線φ_m（ｔ）を生成し出力する。

時間−波数対応データ実測手段１００９では、波長掃引光源１００１のトリガ信号を使って、光ファイバ１０１９を通して入力された光から疎な時間−波数対応データｋ_s（ｔ）を生成し出力する。時間−波数対応データ外内挿手段１０１３では、時間−波数対応データの時刻ごとに位相変化曲線φ_m（ｔ）を分割し、分割された位相変化曲線を使って、時間−波数対応データｋ_s（ｔ）を補間する（離散的なｋ_s（ｔ）間を内挿、または、ｋ_s（ｔ）が存在する時間の外側を外挿）。補間されたデータ、つまり、密な時間−波数対応データｋ（ｔ）を出力する。

時間−波数対応データ実測手段１００９は、実施例２で説明した時間−波数対応データ実測手段３０９を使用しても良い。また、１または複数の固定の帯域通過中心波長を有する狭帯域バンドパスフィルタを使って、トリガ信号の発生時刻からこれらのバンドパスフィルタを通った光の時刻（例えば、フィルタ通過後の光強度がピークとなる時刻）までの差の時間と、それぞれのバンドパスフィルタの中心波数の対応したものを時間−波数対応データとしてもよい。位相変化曲線取得手段１０１２は、実施例１で説明した位相変化曲線取得手段１１２を使用しても良い。

以下、動作手順を説明する。ただし、波長掃引光源１００１は動作をしており、トリガ信号を出し続けているものとする。また、干渉計１００２の受光素子１００７からも干渉信号（アナログ）１０１４が出続けているものとする。

［ステップＳ−１１０１（時間−波数対応データ実測）］時間−波数対応データ実測手段１００９は、光ファイバ１０１９を通ってきた掃引光とトリガ信号から時間−波数対応データｋ_s（ｔ）を生成する。

［ステップＳ−１１０２（干渉波形信号取得）］サンプル光側のｚの位置に鏡などの反射面を置き、干渉計１００２で生成した干渉信号１０１４を、干渉波形信号取得手段１０１１が取得して干渉信号ｉ_m（ｔ）を出力する。受光素子１００７がアナログ出力であり、アナログ出力の信号を干渉波形信号取得手段１０１１が受ける場合は、Ａ／Ｄ変換を行い、ディジタル信号として出力する。

干渉波形信号取得手段１０１１が干渉信号を取得するタイミングは、波長掃引光源が短波長から長波長へ、または、逆に長波長から短波長へ掃引を開始するタイミングとする。例えば、波長掃引光源から掃引開始時刻を示すトリガ信号を、干渉波形信号取得手段１０１１が得るようにしても良い。または、トリガ信号を出力するためのトリガ信号出力手段を別途用意し、波長掃引光源にトリガ信号を入力して、トリガ信号のタイミングで（あるいは同期して）波長掃引光源が掃引を開始しても良い。

［ステップＳ−１１０３（位相変化曲線取得）］位相変化曲線取得手段１０１２は、干渉波形信号取得手段１０１１が出力した干渉信号ｉ_m（ｔ）から位相変化曲線φ_m（ｔ）を生成し出力する。

［ステップＳ−１１０４（時間−波数対応データ外内挿）］時間−波数対応データ実測手段１００９から得た疎な時間−波数対応データｋ_s（ｔ）は離散的なデータであり、（時刻、波数）という組み合わせデータｋ_s（ｔ_i）の集合である。ここで、ｉは整数であり、ｋ_s（ｔ_i）は時刻ｔ_iと波数の対応した離散的なデータであることを示している。

時間−波数対応データ外内挿手段１０１３は、時間−波数対応データ実測手段１００９から得た疎な時間−波数対応データｋ_s（ｔ_i）の時刻ｔ_iごとに、位相変化曲線取得手段１０１２が出力した位相変化曲線φｍ（ｔ）を分割する。説明のため、時刻ｔ_iのｉが０〜ｎであるとし、時刻ｔ_iから時刻ｔ_i+1区間の位相変化曲線をφｍｔ_i，ｔ_i+1（ｔ）とする。また時刻ｔ₀以前のものをφｍｓ，ｔ₀（ｔ）、時刻ｔ_n以降のものをφｍｔ_n，e（ｔ）とする。それら分割した位相変化曲線φｍｔ_i，ｔ_i+1（ｔ）から分割した波数変化曲線ｋｔ_s，ｔ₀（ｔ）、ｋｔ_i，ｔ_i+1（ｔ）、ｋｔ_n，ｔ_e（ｔ）を算出する。ただし、iは０〜ｎである。この方法としては、例えば、分割した位相変化曲線を２ｚで割ることが好例である。

図１５に、疎な時間−波数対応データの外挿または内挿の例を示す。時間−波数対応データ外内挿手段１０１３は、ｋｔ_i，ｔ_i+1（ｔ）の時刻ｔ_i、ｔ_i+1の波数をそれぞれｋ_s（ｔ_i）、ｋ_s（ｔ_i+1）とし、ｋ_s（ｔ_i）を内挿、および、外挿する。

分割された波数変化曲線ｋｔ_s，ｔ₀（ｔ）、ｋｔ_i，ｔ_i+1（ｔ）、ｋｔ_n，ｔ_e（ｔ）（ただし、iは０〜ｎ）は、干渉信号に重畳するノイズ等の影響で、ｔ_i〜ｔ_i+1の変化量が想定しているδｋと厳密に一致しない場合がある。この場合は、時間−波数対応データ外内挿手段１０１３は、干渉信号に時間的に均等にノイズが重畳して均等に波数変化曲線が伸縮していると仮定して、分割された波数変化曲線ｋｔ_i，ｔ_i+1（ｔ）を補正して整合をとる。例えば、分割された波数変化曲線のｔ_i〜ｔ_i+1の変化量をδｋｔ_i，ｔ_i+1＝ｋｔ_i，ｔ_i+1（ｔ_i+1）−ｋｔ_i，ｔ_i+1（ｔ_i）とすると、補正された（分割された）波数変化曲線ｋ’ｔ_i，ｔ_i+1（ｔ）は下記の式とする。

時間−波数対応データ外内挿手段１０１３は、すべての（疎の）時間−波数対応データｋ_s（ｔ_i）について、内挿・外挿が済めば、内挿・外挿した時間−波数対応データを出力して終了する。

上記は内挿と外挿の両方を行う方法について述べたが、必要に応じて内挿のみを行ってもよい。例えば、（疎の）時間−波数対応データｋ_s（ｔ_i）を２つだけ取得し、その間を内挿してもよい。また、逆に、（疎の）時間−波数対応データｋ_s（ｔ_i）を１つだけ取得し、ｋ_s（ｔ_i）の前後を外挿してもよい。

図１６に、時間−波数対応データ取得手段におけるマッハツェンダ干渉計の構成を示す。図１３に示した干渉計１００２には、図１と同様に、マイケルソン干渉計を示したが、マッハツェンダ干渉計を用いても良い。

サンプル１３０３は、ブロック１（１３２０）とブロック２（１３２１）から構成される。ブロック１の面１（１３２２）と面２（１３２３）とは平行である。ブロック２の面１（１３２４）と面２（１３２５）も平行である。ブロック１の面２（１３２３）とブロック２の面１（１３２４）は平行に配置されている。ブロック１（１３２０）の面２（１３２３）とブロック２（１３２１）の面１（１３２４）の間は、真空、または、観測対象の屈折率の波長分散と比較して波長分散の小さいもの（例えば、観測対象がシリコンの場合は、空気や石英等）を充填する。ブロック１の面２（１３２３）を基準面１３１１とし、ブロック２の面１（１３２４）をサンプル内反射面１３０８とする。基準面１３１１とサンプル内反射面１３０８の間隔をｚとする。このｚは既知であるとする。ブロック１の面１（１３２２）やブロック２の面２（１３２５）でも反射が生じ、それらの反射光と、基準面１３１１やサンプル内反射面１３０８からの反射光の４つの反射光がお互いに干渉し合うことになる。しかし、波長掃引光源のコヒーレンス長Ｌｃに対してブロック１、２内の光路長が十分に長くなるようにし、ｚを十分小さくすることにより、基準面１３１１の反射光とサンプル内反射面１３０８の反射光との干渉光以外の干渉強度を十分小さくすることができる。例えば、ｚ＝Ｌｃ／８とし、ブロック１、２内の光路長をＬｃとすると基準面１３１１の反射光とサンプル内反射面１３０８の反射光の干渉光よりも、それ以外の干渉光の強度を約４０ｄＢ小さくすることが可能である。

また、ブロック１の面１（１３２２）とブロック２の面２（１３２５）に反射防止コート（ＡＲコート）をすれば、基準面１３１１の反射光とサンプル内反射面１３０８の反射光強度をブロック１の面１（１３２２）とブロック２の面２（１３２５）より桁違いに大きくすることが可能となる。このため、基準面１３１１の反射光とサンプル内反射面１３０８の反射光による干渉信号（光路長ｚの干渉信号）をノイズなく受光素子１３０７から取り出すことが可能となる。

実施例２、３に引き続き、時間−波数対応データｋ（ｔ）の取得手段の具体的な方法について述べる。ここでは、干渉計の光路長差を生じる媒質が波長分散媒質で構成され、その屈折率増加関数ｒ（ｋ）は既知であり、基準波数ｋ₀の時の光路長ｚ₀は既知であり、干渉信号ｉ_m（ｔ）から算出した位相φ_m（ｔ）が与えられた時に、ｋ（ｔ）を求める方法について述べる。

ｋ（ｔ）を得る原理を以下に示す。式（１５）を再掲すると以下のようになる。

ある時刻ｔ＝ｔ_aにおいて式（４−３）をｋについて解くことを考える。もし、ｒ（ｋ）がｋの０次関数、つまり、ｒ（ｋ）＝ｒ₀であるか、それに近似できるとすると、θ（ｋ）＝２・ｚ₀・ｒ₀・ｋとなるので、式（４−３）からｋは以下のように算出される。

もし、ｒ（ｋ）がｋの１次関数、つまり、ｒ（ｋ）＝ｒ₀＋ｒ₁ｋであるか、それに近似できるとすると、θ（ｋ）＝２・ｚ₀・（ｒ₀＋ｒ₁ｋ）・ｋとなるので、式（４−３）からｋは以下のように算出される。

式（４−５）は２つの解が存在するが、実際の値に応じて、物理的に意味のあるものを選択する。ｒ（ｋ）の近似式中に現れる波数ｋは０以上を前提としているので、式（４−５）の右辺が０以上となる場合に意味がある。ここで、φ_m（ｔ）、ｚ₀、ｒ₀は０以上であることから、（細かい議論は割愛するが）結果的にはｒ₁≦０、かつ下記の式が物理的に意味のある式となる。

もし、ｒ（ｋ）がｋの２次関数や３次関数であるか、それに近似できるとすると、θ（ｋ）はｋの３次関数や４次関数であるか、それに近似できるので、それぞれカルダノの方法やフェラーリの方法を使ってｋ（ｔ）を得る。

ｒ（ｋ）がｋの４次以上の関数や、ｋの多項式ではない関数であるか、それに近似できるとすると（あるいは、ｒ（ｋ）がｋの１、２、３次の関数であっても）、

となる式のｋを求めればよい。このような非線型方程式の解法としては、例えば、ニュートン・ラフソン（Newton-Raphson）法、割線法（secant method）、二分法、ＤＫＡ（Durand-Kerner-Aberth）法等の反復解法がある。

式（４−６）からｋ（ｔ_α）を求める例としてニュートン法を使う場合について説明する。ただし、波数ｋの移動量Δｋについては、φ_m（ｔ_α）をｋに対しては定数とみなして、Δｋ＝−（θ（ｋ（ｔ_α））―φ_m（ｔ_α））／（（∂θ／∂ｋ）｜ｔ＝ｔ_α）とする。ここで、（∂θ／∂ｋ）｜ｔ＝ｔ_αは、ｔ＝ｔ_α時の（∂θ／∂ｋ）の値である。反復処理の終了を判断するためのΔｋの閾値をδとする。

処理１）処理を始める前に終了条件として使用するパラメータδを入力する。
処理２）φｍ（ｔ_α）／（２ｚ₀）を計算し、初期値の波数ｋ（ｔ_α）とする。
処理３）Δｋを計算する。
処理４）Δｋ≦δであれば終了（そのときのｋ（ｔ_α）が求めるｋ（ｔ_α））。そうでなければ、処理５を行う。
処理５）ｋ（ｔ_α）＋Δｋをｋ（ｔ_α）に置き換え、処理３に戻る。

なお、ニュートン・ラフソン法を使用すると、ｒ（ｋ）がｋ（ｔ_α）の１次関数の場合はθ（ｋ）がｋ（ｔ_α）の２次関数となるので、１回の反復でｋ（ｔ_α）を算出できる。

ところで、式（４−１）の左辺φ_m（ｔ）は、観測した干渉信号ｉ_m（ｔ）から、Ｓ−２２０２〜Ｓ−２２０５によって得られるが、ｉ_m（ｔ）の観測誤差やｉ_m（ｔ）に含まれる雑音によって、または、Ｓ−２２０５で偏角を求める際のアンラップ処理において、φ_m（ｔ）に誤差が含まれる可能性がある。これを考慮して、1つ以上の時間−波数対応データを、実施例２、３に記載した方法で実測し、φ_m（ｔ）を補正する。補正方法としては、光路長差２ｚ₀と掃引スタート時の波数ｋ₀に誤差が含まれていると仮定し、２ｚ₀に関してはφ_m（ｔ）の係数、ｋ₀に関してはφ_m（ｔ）のバイアス値を、実測のｋ（ｔ）に近づける方法を取る。

具体的には、下記Ｊが最小となる係数ｐ、ｑを求める（最小２乗法）。

ここで、ｔ_iは時間−波数対応データｋ（ｔ_i）の取得時刻であり、ｉ＝０〜Ｎ−１のＮ個の時間−波数対応データｋ（ｔ_i）を取得し、φ_m（ｔ）の補正に使用していることを示している。また、φ_m（ｔ_i）はその時刻ｔ_i時のｉ_m（ｔ）の位相である。

ｐ、ｑを求める方法としては、上記Ｊのｐ、ｑそれぞれで偏微分した式が０となる連立方程式を解くことによりｐ、ｑを求める方法がある。具体的には下記の正規方程式を解くことによってｐ、ｑを求める。

このようにして求めたｐ、ｑを使い、下記のように位相を補正し、補正された位相変化曲線φ_m’（ｔ）を得る。

以上の原理にを使った、時間−波数対応データを取得する方法について説明する。図１７に、時間−波数対応データ取得手段の機能ブロックを示し、図１８にフローチャートを示す。

図１７の１４０１，１４０９〜１４１２は、それぞれ図１３の１００１，１００９〜１０１２と同じであり、また、図１８のＳ−１５０１〜Ｓ−１５０３は、それぞれ図１４のＳ−１１０１〜Ｓ−１１０３と同じであるため、符号１４０１，１４０９〜１４１２で表される構成要素とＳ−１５０１〜Ｓ−１５０３の手順の説明と割愛し、図１３にはない、屈折率増加関数生成手段１４２１、位相変化曲線補正手段１４２２、時間−波数対応データ生成手段１４２３、および、図１４にはない、Ｓ−１５０４、Ｓ−１５０５、Ｓ−１５０６について述べる。以下、Ｓ−１５０４、Ｓ−１５０５、Ｓ−１５０６の手順について述べる。

［ステップＳ−１５０４（屈折率増加関数生成）］屈折率増加関数生成手段１４２１は、屈折率分散データｎ（ｋ）と基準波数ｋ₀が入力され、屈折率増加関数ｒ（ｋ）＝ｎ（ｋ）／ｎ（ｋ₀）を計算し、ｒ（ｋ）を出力する。

［ステップＳ−１５０５（位相変化曲線補正）］位相変化曲線補正手段１４２２は、時間−波数対応データ実測手段１４０９から得た時間−波数対応データｋ_s（ｔ）と、位相変化曲線実測手段１４１０から得た位相変化曲線φ_m（ｔ）と、屈折率増加関数生成手段１４２１から出力された屈折率増加関数ｒ（ｋ）と、光路長差２ｚ₀とが入力され、式（４−７）のＪが最少となるｐ、ｑを算出し、ｐ、ｑを使って式（４−１１）によって補正された位相変化曲線φ_m’（ｔ）を出力する。位相変化曲線補正手段１４２２とステップＳ−１５０５については、後で詳述する。

［ステップＳ−１５０６（時間−波数対応データ生成）］時間−波数対応データ生成手段１４２３は、屈折率増加関数生成手段１４２１から出力された屈折率増加関数ｒ（ｋ）と、位相変化曲線補正手段１４２２から出力された位相変化曲線φ_m’（ｔ）と、光路長差２ｚ₀とが入力され、ｒ（ｋ）の近似式を求め、φ_m’（ｔ）の時間範囲において、各ｔにおけるｋ（ｔ）を、φ_m’（ｔ）と２ｚ₀とｒ（ｋ）を算出し、算出したｋ（ｔ）を出力する。時間−波数対応データ生成手段１４２３とステップＳ−１５０６については、後述する。

図１９に、位相変化曲線補正手段の機能ブロックを示し、図２０に、フローチャートを示す。以下、動作手順を示す。

［ステップＳ−１７０１（補正係数算出）］補正係数算出手段１６０１は、時間−波数対応データｋ_s（ｔ）と位相変化曲線φ_m（ｔ）と光路長差２ｚ₀が入力され、式（４−８）をｐ、ｑについて解くか、または、式（４−９）や式（４−１０）の計算をしてｐ、ｑを算出し、出力する。ここで、ｋ_s（ｔ）のデータはＮ個あり、それぞれの時刻はｔ_i（ｉは０〜Ｎ−１）で表される。また、上式（４−７）〜（４−１０）にあるｋ（ｔ_i）は、ｋ_s（ｔ_i）を代入して計算する。式（４−８）をｐ、ｑについて解く方法としては、ガウスの消去法等の直接法や、ガウス・ザイデル法、ＳＯＲ法等の反復法や、共役勾配法等がある。

［ステップＳ−１７０２（補正係数演算手段）］補正係数演算手段１６０２は、式（４−１１）の計算をして位相変化曲線φ_m（ｔ）を補正し、補正された位相変化曲線φ_m’（ｔ）を出力する。

図２１に、時間−波数対応データ生成手段の機能ブロックを示し、図２２に、フローチャートを示す。以下、動作手順を示す。

［ステップＳ−１９０１（近次関数パラメータ算出）］近似関数パラメータ算出手段１８０１は、屈折率増加関数ｒ（ｋ）が入力され、ｒ（ｋ）の近似関数ｒ_App（ｋ）のパラメータを算出し、出力する。

近似関数は、予め所定の関数が決められている。例えば、ｋの０次関数や１次関数といった具合に、ｋのｍ次関数（ｍは０以上の整数）とすることが考えられる。この場合、近似関数のパラメータは、ｍ次関数の場合は、ｒ_mｋ^m＋ｒ_m-1ｋ^m-1＋…＋ｒ₀のｒ_m〜ｒ₀のことである。これらのパラメータは、ｒ（ｋ）を使って最小二乗法で求めることができる。

また、ｋのｍ次関数以外の関数の線形結合、つまり、ｒ_mΨ_m（ｋ）＋ｒ_m-1Ψ_m-1（ｋ）＋…＋ｒ₀Ψ₀（ｋ）の場合も同様である。この場合、Ψ_m（ｋ）〜Ψ₀（ｋ）は予め決めておく。（なお、Ψ_i（ｋ）＝ｋⁱとすると、ｋのｍ次の関数となる）この場合、次の正規方程式を係数ｒ_m〜ｒ₀について解くことで係数ｒ_m〜ｒ₀を求めることができる。

ただし、上式は波数ｋがＭ個の離散的に与えられ、ｋ_i（ｉ＝０〜Ｍ）としている。式（４−１２）を係数ｒ_m〜ｒ₀について解く方法としては、ガウスの消去法等の直接法や、ガウス・ザイデル法、ＳＯＲ法等の反復法や、共役勾配法等がある。

上記は、近似関数を与えた方法であったが、最適な近似関数を自動的に得る方法が考えられる。例えば、Ψ_m（ｋ）〜Ψ₀（ｋ）を予め決めておき（例えば、Ψ_i（ｋ）＝ｋⁱ等）、Ψ₀（ｋ）からｒ_mΨ_m（ｋ）＋ｒ_m-1Ψ_m-1（ｋ）＋…＋ｒ₀Ψ₀（ｋ）までをｒ（ｋ）に当てはめた近似式のそれぞれのパラメータ（係数）ｒ_iとＲ²を計算し、Ｒ²が最も小さい近似式を採用しても良い。Ｒ²の代わりに、近似式と実データとの差の２乗和（Σ｛ｒ（ｋ_i）−ｒ_mΨ_m（ｋ_i）＋ｒ_m-1Ψ_m-1（ｋ_i）＋…＋ｒ₀Ψ₀（ｋ_i）｝²（残差ともいう）を使用しても良い。なお、Ｒ²は決定係数と呼ばれるものであり、その定義は多数存在する。以下にその一例を示す。

ただし、ｆ_i＝ｒ_mΨ_m（ｋ_i）＋ｒ_m-1Ψ_m-1（ｋ_i）＋…＋ｒ₀Ψ₀（ｋ_i）であり、ｋ_meanはｋ_i（ただし、ｉ＝０〜Ｍ）の平均値である。また、式（４−１３）の分子、および、分母を自由度で調整した、下記のような例もある。

式（４−１４）は自由度調整済みの決定係数（adjusted R2）と呼ばれる。

［ステップＳ−１９０２（波数算出）］波数算出手段１８０２は、補正された位相変化曲線φｍ’（ｔ）と、近似関数算出手段１８０１から得た屈折率増加関数近似関数ｒ_App（ｋ）と、光路長差２ｚ₀が入力され、ｒ_App（ｋ）に応じて時間−波数対応データｋ（ｔ）を算出し、出力する。
・ｒ_App（ｋ）がｋの０次関数の場合は、式（４−４）に示す式によりｋ（ｔ）を算出する。
・ｒ_App（ｋ）がｋの１次関数の場合は、式（４−５）’に示す式によりｋ（ｔ）を算出する。
・ｒ_App（ｋ）がｋの２次関数の場合は、式（４−１）の右辺がｋの３次関数となるが、φ_m（ｔ）を定数として、カルダノの方法によりｋ（ｔ）を算出する。
・ｒ_App（ｋ）がｋの３次関数の場合は、式（４−１）の右辺がｋの４次関数となるが、φ_m（ｔ）を定数として、フェラーリの方法によりｋ（ｔ）を算出する。
・ｒ_App（ｋ）がｋの４次関数以上、または、ｋの多項式ではない関数であるとすると、式（４−６）が成り立つｋ（ｔ）を算出する。このような非線型方程式の解法としては、例えば、ニュートン・ラフソン（Newton-Raphson）法、割線法（secant method）、二分法、ＤＫＡ（Durand-Kerner-Aberth）法等の反復解法を用いる。

上記のｋ（ｔ）の計算は、時間ｔを少しずつ変えながら、φｍ’（ｔ）のデータの全時間領域内で行う。例えば、φｍ’（ｔ）が、時間ｔ_i（ｉ＝０〜Ｌ）それぞれに位相データを持つ離散データであり、ｔ_i+1＝ｔ_i＋Δｔ（ただし、Δｔ＞０）という場合は、ｔ₀から順番にｔ_Lまで、それぞれのｔ_iに対して上記のｋ（ｔ_i）を算出し、ｔ_iとｋ（ｔ_i）を組みにしたデータを時間−波数対応データとして出力しても良い。

屈折率は、サンプルの温度によって変化する。ここでは、実施例１における、ステップＳ−２０７波長分散補正係数算出において、波長分散補正係数算出手段１０７に入力する屈折率分散データｎ_a（ｋ）、ｎ_b（ｋ）を、サンプルａ、ｂの温度に応じて補正する装置について述べる。

図２３に、屈折率分散データの温度補正手段の機能ブロックを示し、図２４に、フローチャートを示す。以下の説明ではサンプルａの屈折率分散データｎ_a（ｋ）について述べるが、サンプルｂの屈折率分散データｎ_b（ｋ）についても同様である。

［Ｓ−７０１（サンプル温度取得）］サンプル温度取得手段６０２がサンプルａの温度Ｔ_aを取得し、屈折率分散温度補正手段６０１へ送る。サンプルａの温度の取得は、サーミスタ等の温度センサをサンプルａに接着して、その温度センサの電気的な特性を検知する装置を用いて温度変化を取得することが考えられる。また、サーモビューワ等、赤外線等の光の分布と強度を電気的に取得して温度に換算する方法でも良い。

［Ｓ−７０２（屈折率分散データ、および、想定温度取得）］屈折率分散データ温度補正手段６０１が、屈折率分散データｎ_a（ｋ）と、該屈折率分散データｎ_a（ｋ）が想定している（屈折率分散データｎ_a（ｋ）を取得した時のサンプルの）温度Ｔ₀を取得する。また、Ｓ−７０１（サンプル温度取得）にてサンプル温度取得手段６０２が取得したサンプル温度Ｔ_aを取得する。

［Ｓ−７０３（屈折率分散データ温度補正）］屈折率分散データ温度補正手段６０１が、サンプルａの温度係数Δｎ（ｋ）を使い、補正された屈折率分散データｎ_a’（ｋ）を算出し、出力する。ｎ_a’（ｋ）は例えばｎ_a’（ｋ）＝ｎ_a（ｋ）＋Δｎ（ｋ）（Ｔ_a−Ｔ₀）のように計算される。

上記は、Ｓ−７０１、Ｓ−７０２、Ｓ−７０３の順での動作を示したが、図７に示すように、Ｓ−７０１とＳ−７０２は並列動作させても良いし、Ｓ−７０２、Ｓ−７０１、Ｓ−７０３の順で動作させても良い。

ＳＳ−ＯＣＴ装置の基本構成は、波長掃引光源、干渉計および信号処理装置から構成され、本実施形態のリスケーリング方法を実行する信号処理装置は、コンピュータとプログラムによっても実現することができ、リスケーリング方法を実行するためのプログラムを、記録媒体に記録することも、ネットワークを通して提供することもできる。

８０１，１０１５，１０１７，１０１９，１３１５，１３１７，１３１９，１４１５，１４１７，１４１９，２４１５，２４１７，２４１９ 光ファイバ
８０２ アイソレータ
８０３ 集光レンズ
８０４ 半導体光増幅器（ＳＯＡ）
８０５ ＳＯＡ反射面
８０６ コリメータレンズ
８０７ 光偏向器
８１０ 電源
８１１ 回折格子
８１２ ミラー
９０１ ＫＴＮ結晶
９０２，９０３ 電極
９０４ 光軸
９０５，１３０９，２００９，２４０９ 出射光
１００１，１４０１，２００１，２４０１ 波長掃引光源
１００２，１３０２，１４０２，２００２，２４０２ 干渉計
１００７，１３０７，１４０７，２００７，２４０７ 受光素子
１０１４，１３１４，１４１４，２０１４，２４１４ 干渉信号
１０１６，１４１６，２００５ ビームスプリッタ
１３０３，２００３，２４０３ サンプル
１３０８，２００８，２４０８ 反射面
１３１０，２０１０，２４１０ サンプル光
１３１１，２０１１，２４１１ 基準面
１３１２，２０１２，２４１２ 参照光
１３１３，２０１３，２４１３ 干渉光
１３１６，２４１６ サーキュレータ
１３２０ ブロック１
１３２１ ブロック２
１３２２，１３２４ 面１
１３２３，１３２５ 面２
２００４，２４０４ 信号処理装置
２００６ 参照光ミラー
２４１８ ファイバコリメータ

Claims (7)

1. 波数が時間的に一方向に連続的に変化する光を出力する波長掃引光源と、該波長掃引光源から出力された光を入力し、サンプルアームと参照アームと間の光路長差に応じた干渉信号を出力する干渉計とを含み、該干渉計から出力された干渉信号を周波数解析して断層信号を出力する光断層撮像装置において、リスケーリング処理を行うための信号処理装置であって、
   サンプルａを用いて単一の光路長差に係る位相変化曲線φ_m(ｔ)を出力する位相変化曲線取得手段と、
   前記サンプルａの屈折率分散データｎ_a（ｋ）、観測対象の屈折率分散データｎ_b（ｋ）、および前記波長掃引光源から出力された光の時間−波数対応データｋ（ｔ）から、波数ｋ₀における屈折率を基準としたときの、サンプルａの時間−屈折率変動データｒ_a(ｔ)（＝ｎ_a（ｋ）／ｎ_a（ｋ₀））に対する観測対象の時間−屈折率変動データｒ_b(ｔ)（＝ｎ_b（ｋ）／ｎ_b（ｋ₀））の比ｒ_b(ｔ)／ｒ_a(ｔ)を時間−補正係数対応データｒ_a→b(ｔ)として出力する波長分散補正係数算出手段と、
   前記時間−補正係数対応データｒ_a→b(ｔ)と前記位相変化曲線φ_m(ｔ)との積を補正された位相変化曲線φ_m’(ｔ)として出力する位相変化曲線補正手段と、
   前記補正された位相変化曲線φ_m’(ｔ)からリサンプリング点ｔ_nを算出するリサンプリング点算出手段と
   を備えたことを特徴とする信号処理装置。
2. 前記波長掃引光源から出力された光に対して、指示された波数ｋを中心としたバンドパスフィルタリングを行う狭帯域バンドパスフィルタと、
   該狭帯域バンドパスフィルタでフィルタリングされた光を検知する光検知器と、
   前記波長掃引光源のトリガ信号から前記光検知器が検知した光のピークまでの時間を計測し、前記指示された波数と対応させた、前記時間−波数対応データｋ（ｔ）を出力する時間計測手段と
   を含む時間−波数対応データ取得手段をさらに備えたことを特徴とする請求項１に記載の信号処理装置。
3. 前記波長掃引光源から出力された光とトリガ信号から疎な時間−波数対応データｋ_s（ｔ）を出力する時間−波数対応データ実測手段と、
   前記波長掃引光源から出力された光を光路長差のある干渉計を通して得た干渉信号ｉ_m（ｔ）から位相変化曲線φｍ（ｔ）を出力する位相変化曲線実測手段と、
   前記疎な時間−波数対応データｋ_s（ｔ）の時刻ｔ_iごとに、前記位相変化曲線実測手段が出力した位相変化曲線φｍ（ｔ）を分割し、分割された位相変化曲線により、前記疎な時間―波数対応データを外挿または内挿して前記時間−波数対応データｋ（ｔ）を出力する時間−波数対応データ外内挿手段と
   を含む時間−波数対応データ取得手段をさらに備えたことを特徴とする請求項１に記載の信号処理装置。
4. 前記波長掃引光源から出力された光とトリガ信号から疎な時間−波数対応データｋ_s（ｔ）を出力する時間−波数対応データ実測手段と、
   前記波長掃引光源から出力された光を光路長差のある干渉計を通して得た干渉信号ｉ_m（ｔ）から位相変化曲線φｍ（ｔ）を出力する位相変化曲線実測手段と、
   屈折率分散データｎ（ｋ）と所定の基準波数ｋ₀から屈折率増加関数ｒ（ｋ）を出力する屈折率増加関数生成手段と、
   前記疎な時間−波数対応データと前記位相変化曲線φｍ（ｔ）と前記屈折率増加関数ｒ（ｋ）と前記光路長差とから補正された位相変化曲線φ_m’（ｔ）を出力する位相変化曲線実測手段と、
   前記屈折率増加関数ｒ（ｋ）の近似関数ｒ_App（ｋ）を生成し、前記補正された位相変化曲線φ_m’（ｔ）と前記光路長差とから前記近似関数ｒ_App（ｋ）に応じた時間−波数対応データｋ（ｔ）を出力する時間−波数対応データ生成手段と
   を含む時間−波数対応データ取得手段をさらに備えたことを特徴とする請求項１に記載の信号処理装置。
5. 前記サンプルａの温度Ｔ_aを取得するサンプル温度取得手段と、
   前記屈折率分散データｎ_a（ｋ）の想定温度Ｔ₀との差に従って、前記サンプルａの温度係数Δｎ（ｋ）を用いて、前記屈折率分散データを補正する屈折率分散データ温度補正手段と
   をさらに備えたことを特徴とする請求項１ないし４のいずれかに記載の信号処理装置。
6. 波数が時間的に一方向に連続的に変化する光を出力する波長掃引光源と、該波長掃引光源から出力された光を入力し、サンプルアームと参照アームと間の光路長差に応じた干渉信号を出力する干渉計とを含み、該干渉計から出力された干渉信号を周波数解析して断層信号を出力する光断層撮像装置における信号処理装置のリスケーリング方法であって、
   前記信号処理装置の位相変化曲線取得手段が、サンプルａを用いて単一の光路長差に係る位相変化曲線φ_m(ｔ)を出力する第１ステップと、
   前記信号処理装置の波長分散補正係数算出手段が、前記サンプルａの屈折率分散データｎ_a（ｋ）、観測対象の屈折率分散データｎ_b（ｋ）、および前記波長掃引光源から出力された光の時間−波数対応データｋ（ｔ）から、波数ｋ₀における屈折率を基準としたときの、サンプルａの時間−屈折率変動データｒ_a(ｔ)（＝ｎ_a（ｋ）／ｎ_a（ｋ₀））に対する観測対象の時間−屈折率変動データｒ_b(ｔ)（＝ｎ_b（ｋ）／ｎ_b（ｋ₀））の比ｒ_b(ｔ)／ｒ_a(ｔ)を時間−補正係数対応データｒ_a→b(ｔ)として出力する第２ステップと、
   前記信号処理装置の位相変化曲線補正手段が、前記時間−補正係数対応データｒ_a→b(ｔ)と前記位相変化曲線φ_m(ｔ)との積を補正された位相変化曲線φ_m’(ｔ)として出力する第３ステップと、
   前記信号処理装置のリサンプリング点算出手段が、前記補正された位相変化曲線φ_m’(ｔ)からリサンプリング点ｔ_nを算出する第４ステップと
   を備えたことを特徴とするリスケーリング方法。
7. 波数が時間的に一方向に連続的に変化する光を出力する波長掃引光源と、該波長掃引光源から出力された光を入力し、サンプルアームと参照アームと間の光路長差に応じた干渉信号を出力する干渉計とを含み、該干渉計から出力された干渉信号を周波数解析して断層信号を出力する光断層撮像装置におけるリスケーリング方法をコンピュータに実行させるためのプログラムであって、前記コンピュータを、
   サンプルａを用いて単一の光路長差に係る位相変化曲線φ_m(ｔ)を出力する位相変化曲線取得手段と、
   前記サンプルａの屈折率分散データｎ_a（ｋ）、観測対象の屈折率分散データｎ_b（ｋ）、および前記波長掃引光源から出力された光の時間−波数対応データｋ（ｔ）から、波数ｋ₀における屈折率を基準としたときの、サンプルａの時間−屈折率変動データｒ_a(ｔ)（＝ｎ_a（ｋ）／ｎ_a（ｋ₀））に対する観測対象の時間−屈折率変動データｒ_b(ｔ)（＝ｎ_b（ｋ）／ｎ_b（ｋ₀））の比ｒ_b(ｔ)／ｒ_a(ｔ)を時間−補正係数対応データｒ_a→b(ｔ)として出力する波長分散補正係数算出手段と、
   前記時間−補正係数対応データｒ_a→b(ｔ)と前記位相変化曲線φ_m(ｔ)との積を補正された位相変化曲線φ_m’(ｔ)として出力する位相変化曲線補正手段と、
   前記補正された位相変化曲線φ_m’(ｔ)からリサンプリング点ｔ_nを算出するリサンプリング点算出手段
   として機能させるためのプログラム。