光断層画像表示システム
【課題】光断層画像表示システムにおいて、高精度で断層画像を表示システムの提供。
【解決手段】周波数領域の光断層画像表示システムにおいて、干渉光の周波数分析を行い、光断層画像を得る。このとき干渉信号を解析対象信号として解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるように周波数ｆ’、振幅Ａ’、初期位相φ’を求め、その解析結果により、断層画像を表示する。これにより高精度の光断層画像を得ることができる。
【解決手段】周波数領域の光断層画像表示システムにおいて、干渉光の周波数分析を行い、光断層画像を得る。このとき干渉信号を解析対象信号として解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるように周波数ｆ’、振幅Ａ’、初期位相φ’を求め、その解析結果により、断層画像を表示する。これにより高精度の光断層画像を得ることができる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理に特徴を有する光断層画像表示システムに関するものである。
【背景技術】
【０００２】
近年内視鏡治療などの医療技術の進歩に伴って、病理組織の診断を非深襲かつリアルタイムに行う診断方法が望まれている。従来例えばＣＣＤを用いた電子内視鏡や、ＣＴ、ＭＲＩ、超音波による画像化が診断方法として用いられている。電子内視鏡は生体の表面の観察に限定され、また後者の画像診断システムはミクロンオーダーの分解能で観察するには技術的な限界があった。このような方法を補完する技術として、光コヒーレンストモグラフィーシステム（ＯＣＴ）が注目されている。
【０００３】
ＯＣＴの中には、時間領域ＯＣＴ（ＴＤ−ＯＣＴ）と周波数領域ＯＣＴ（ＦＤ−ＯＣＴ）の２種類があり、またＦＤ−ＯＣＴの中にもスペクトロメータタイプ（ＳＤ−ＯＣＴ）と波長走査型光源タイプ（ＳＳ−ＯＣＴ）の２つがある。従来の時間領域ＯＣＴの場合には、広帯域な光を当て生体内部からの反射光の干渉成分を周波数分析していたが、この方法だと干渉光の中に異なる深さからの反射光も重なりあうために、ある特定の深さからの信号光だけを感度良く検出できなかった。
【０００４】
一方スペクトロメータタイプのＯＣＴは、特許文献１に示されているように、白色光源を生体に光を照射し、参照光と生体内の異なる深さから戻ってくる反射光とを干渉計で干渉させる。そして反射光を回折格子等で分光して１次元ＣＣＤで受光し、電気信号とする。波長走査型ＯＣＴは、特許文献２に記されているように、生体に光を照射し、照射光の波長を連続的に変化させ、参照光と生体内の異なる深さから戻ってくる反射光とを干渉計で干渉させる。これらはいずれも干渉信号の周波数成分を分析することによって、断層画像を得るシステムである。ＦＤ−ＯＣＴは得られた干渉信号光をフーリエ変換し、断層画像を得るようにしている。この技術は高分解能の断層画像を構築することができるため、高度なシステムとして期待されており、理論的に時間領域ＯＣＴの１００倍以上の感度が得られる。波長走査型ＯＣＴは測定感度も高く、動的ノイズに強いという点で内視鏡などの実使用に好適である。ここで照射する光の波長走査の帯域が広いほど周波数分析の帯域が上がるので、深さ方向の分解能が上がる。
【０００５】
しかしながらＦＤ−ＯＣＴはフーリエ変換時に窓の影響によりメインローブ以外にもサイドローブが発生するため、分解能に限界があるという欠点があった。以下この問題点について説明する。
【０００６】
従来から、フーリエ変換を用いた正弦波モデルによる信号解析方法としては、非特許文献１にも紹介されているように、様々な手法が提案されているが、その多くは、信号の周波数を事前に推定し、その推定した周波数における最適な振幅や初期位相を求めることによって行われている。
【０００７】
信号は複数の周波数の正弦波の重ね合わせであり、信号を解析するには、当該信号を構成する様々な周波数の正弦波毎に振幅や初期位相等の周波数パラメータを推定することになる。このような周波数解析手法としては、フーリエ変換が広く知られている。フーリエ変換における周波数パラメータの推定は、周期的な信号の場合には、次式（１）に示すように、有限長の窓によるフーリエ変換式を用いて行うことができ、この次式（１）に基づいて、特定の周波数ｆ（＝ｎ／Ｔ：ｎは整数、Ｔは解析窓長）のパラメータ（振幅Ａ及び初期位相φ）を推定することができる。ただし、振幅Ａと初期位相φとＸ（ｆ）との関係は次式（２）である。
【０００８】
【数１】
【数２】
【０００９】
このようなフーリエ変換によって推定することが可能な周波数パラメータは、解析窓長Ｔに依存し、１／Ｔの整数倍の周波数分しか計算することができない。したがって、上式（１）を用いてパラメータを推定している限り、周波数ｆは、原理的に離散的となり、非周期信号の周波数成分を特定することはできない。また、フーリエ変換を用いて周波数解析を行っている分野では、そのほとんどが、周期を持たない非周期信号を解析対象とするにもかかわらず、無限長の積分区間にわたる積分計算が極めて困難であることから、上式（１）を用いなければならないのが現状である。
【００１０】
この問題を解決するために、窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みがある。
【００１１】
例えば、非特許文献２に記載の短時間フーリエ変換（Short-Time Fourier Transform、Short-Term Fourier Transform；ＳＴＦＴ）は、音響信号を短時間の解析窓によって複数の信号に分割し、各信号に対して独立にフーリエ変換を行う手法である。
【００１２】
また、非特許文献３乃至非特許文献６に記載のＡＢＳ（Analysis by Synthesis）法やＧＨＡ（Generalized Harmonic Analysis）等においては、最適な周波数を定めるための周波数を複数用意して振幅及び初期位相を求め、その中から最適な周波数と振幅と初期位相との組み合わせを選択している。
【先行技術文献】
【特許文献】
【００１３】
【特許文献１】特開２００６−１３２９９６号
【特許文献２】特開２００７−０２４６７７号
【非特許文献】
【００１４】
【非特許文献１】東山三樹夫、小池恒彦、「高い周波数分析精度の信号分析手法」、日本音響学会論文誌、１９８８年、第５４巻、第８号
【非特許文献２】L. R. Rabiner and R. W. Schafer, "Digital Processing of Speech Signals", Prentice-Hall, Englewood Cliffs, NJ, 1978年
【非特許文献３】T. F. Quatieri and R. J. McAulay, "Speech transformations based on a sinesoidal representation", IEEE Trans. Acoust. Speech Signal Process. ASSP 34, 1986年, p.1449-1464
【非特許文献４】E. B. George and M. J. Smith, "Analysis-by-synthesis/overlap-add sinusoidal representation", J. Audio Eng. Soc. 40, 1992年, p.497-515
【非特許文献５】E. B. George and M. J. Smith, "Speech analysis/synthesis and modification using an analysis-by-synthesis/overlap-add sinusoidal model", IEEE Trans. Speech Audio Process. 5, 1997年, p.389-406
【非特許文献６】T. Terada, H. Nakajima, M. Tohyama and Y. Hirata, "Non-stationary waveform analysis and synthesis using generalized harmonic analysis", IEEE-SP, Int. Symp. TF/TS Analysis, 1994年, p.429-432
【発明の開示】
【発明が解決しようとする課題】
【００１５】
しかしながら、上述した非特許文献２に記載の短時間フーリエ変換においては、短時間の解析窓長の逆数の整数倍の周波数でしか、周波数パラメータを推定することができないという問題がなおも存在する。
【００１６】
また、上述した非特許文献３乃至非特許文献６に記載のＡＢＳ法やＧＨＡ等においては、元の信号に対して、予め用意された周波数と異なる周波数が存在する場合には、誤った組み合わせが複数検波されてしまい、解析精度が低下するという問題がある。
【００１７】
これら窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みは、
１）周波数分解能を上げるためには解析窓長を長くしなければならないが、時間的な周波数変動をともなう信号では周波数が平均化され、時間的な変動を捉えることが困難となること
２）一方、解析窓長を短くすると、周波数分解能が低下すること
という問題が指摘されている。
【００１８】
このように、周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にあり、信号解析上の大きな問題となっていた。
【００１９】
本発明はこのような実状に鑑みてなされたものであり、光断層画像表示システムにおいて周波数分解能が解析窓長に依存せずに高い周波数分解能を有し、高精度に周波数解析を行うことができ、高分解能、高感度、高速で画像表示することができる光断層画像表示システムを提供することを目的とする。
【課題を解決するための手段】
【００２０】
この課題を解決するために、本発明の光断層画像表示システムは、周期的に光の発振波長を走査する波長走査型光源と、前記波長走査型光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、前記干渉光学計より得られる干渉光を受光し、ビート信号を得る受光素子と、前記受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、前記信号解析部は、前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、前記演算手段より得られる断面画像を表示する表示部と、を具備するものである。
【００２１】
この課題を解決するために、本発明の光断層画像表示システムは、所定の範囲で一様な波長の光を発振する光源と、前記光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、前記干渉光学計から得られる干渉光を波長に基づいて分岐する波長分光手段と、前記波長分光手段より分光された範囲の干渉光を受光して電気信号に変換する一次元受光素子と、前記一次元受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、前記信号解析部は、前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、前記演算手段より得られる断面画像を表示する表示部と、を具備するものである。
【００２２】
ここで前記演算手段は、前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’のそれぞれについて適切な初期値を求め、求めた初期値から非線形方程式の最適解に収束させて前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を求めるようにしてもよい。
【００２３】
ここで前記演算手段は、前記正弦波モデル信号の位相を構成する前記周波数ｆ’及び前記初期位相φ’について最急降下法を適用して収束させ、当該周波数ｆ’及び当該初期位相φ’を求め、前記正弦波モデル信号の係数としての前記振幅Ａ’について最急降下法を適用して収束させ、当該振幅Ａ’を求めるようにしてもよい。
【００２４】
ここで前記演算手段は、前記最急降下法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させた後、さらに、ニュートン法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させるようにしてもよい。
【００２５】
ここで前記波長走査型光源の１走査の期間内に、前記波長走査型光源の光の等周波数間隔でのトリガ信号を発生するｋトリガ発生部を更に有し、前記信号解析部の記憶部は、前記ｋトリガ発生部により発生したｋトリガ信号のタイミングで前記干渉信号を保持するものとしてもよい。
【００２６】
ここで前記信号解析部の演算手段は、物体の測定位置に反射物を配置したときに、前記受光素子からの出力をデータ番号を付して等時間間隔で取り込み、前記ビート信号のピークとなるデータ番号を算出し、ピーク番号に対するデータ番号を示す近似曲線を算出し、前記ピーク番号変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出しデータ取り込み数列とし、測定位置に測定物体を配置したときの前記受光素子からの等時間間隔の出力から、前記データ取り込み数列の各要素でのタイミングの信号レベルを算出し、得られた等周波数間隔のデータを前記記憶手段に保持するものとしてもよい。
【発明の効果】
【００２７】
このような本発明による光断層画像表示システムによれば、干渉光信号を周波数分析する際に非周期信号を想定した周波数解析を目的として定式化を行い、非周期信号のフーリエ変換式のパラメータを求める問題を非線形方程式の最適解を求める問題に置き換えることにより、周波数分解能が解析窓長に依存しない。従って従来のフーリエ変換による周波数解析手法における精度と比べて１０万〜１００億倍以上という驚くべき高精度で周波数解析を行うことができ、高分解能で断層画像を表示することができるという効果が得られる。
【図面の簡単な説明】
【００２８】
【図１】図１は本発明の第１の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。
【図２】図２は本発明の実施の形態による光断層画像表示システムのｋトリガ発生部の構成を示すブロック図である。
【図３】図３は本発明の実施の形態による光断層画像表示システムの信号解析部の構成を示すブロック図である。
【図４】図４は走査角度と発振周波数の関係の一例を示すグラフである。
【図５】図５は本周波数解析手法とＤＦＴとの違いを説明するための図である。
【図６】図６は本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用した場合の具体例を示す図である。
【図７】図７は本周波数解析手法とＤＦＴとＧＨＡとの違いを説明するための図であり、各手法の誤差を求めた結果を示す図である。
【図８】図８は本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用して得られるスペクトルの解析精度の具体例を示す図である。
【図９】図９は本発明の第２の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。
【符号の説明】
【００２９】
１０ 波長走査型光源
１１，１４，１８ 光ファイバ
１５，１９，２２ レンズ
１６ 参照ミラー
２０ スキャニングミラー
２３ フォトダイオード
２４ 増幅器
２５ 信号解析部
２６ スキャントリガ発生部
２７ ｋトリガ発生部
３１ クロック発生回路
３２ トリガ生成回路
３３ メモリ
３４ ＲＷ制御部
４０ 光源
４１ 回折格子
４２ 一次元ＣＣＤ
４３ 増幅器
５０ 入力インターフェイス
５１ ＣＰＵ
５２ ＲＯＭ
５３ ＲＡＭ
５４ 記憶部
５５ 入力操作制御部
５６ 表示部
【発明を実施するための形態】
【００３０】
（第１の実施の形態）
［光断層画像表示システムの構成］
図１は本発明の第１の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。本実施の形態はＳＳ−ＯＣＴによる光断層画像表示システムであり、光源として波長走査型光源を用いる。波長走査型光源１０には一定の周波数範囲の光信号を発振する波長可変型のレーザを用いる。このレーザの出力は光ファイバ１１を介して分岐部１２に与えられる。光ファイバ１１の他端には光増幅器１３が設けられる。光増幅器１３は波長走査型光源１０のレーザ光をそのまま増幅するものであり、その出力は光ファイバ１４に与えられる。光ファイバ１４の他端にはコリメートレンズ１５及び参照ミラー１６を設ける。又この光ファイバ１４の中間部分には、他の光ファイバ１８を接近させて干渉させる結合部１７が設けられる。光ファイバ１８の一端には、波長走査型光源１０から結合部１７を介して得られた光信号を平行光とするコリメートレンズ１９、光をスキャニングするスキャニングミラー２０が設けられる。スキャニングミラー２０は紙面に垂直な軸を中心にして一定範囲で回動することによって平行光の反射角度を変化させるものである。集束レンズ２１はこの反射光を受光する位置に配置し、測定部位へ光を集束すると共に水平方向にスキャニング（走査）する。ここで結合部１７から参照ミラー１６までの光学距離Ｌ１と、結合部１７から測定部位の表面までの光学距離Ｌ２とを等しくしておく。さて光ファイバ１８の他端にはレンズ２２を介してフォトダイオード２３を接続する。フォトダイオード２３は、参照ミラー１６からの反射光と測定部位で反射された光の干渉光を受光することによって、そのビート信号を電気信号として得る受光素子である。ここで光ファイバ１４，１８と結合部１７、コリメートレンズ１５、参照ミラー１６、コリメートレンズ１９、スキャニングミラー２０、集束レンズ２１は干渉光学計を構成している。
【００３１】
さてフォトダイオード２３の出力は増幅器２４を介して信号解析部２５に入力される。信号解析部２５は後述するように干渉計から得られる受光信号を周波数領域の信号に変換することによって、断層画像信号を得るものである。
【００３２】
又波長走査型光源１０の出力の一部は分岐部１２によって分岐されてスキャントリガ発生部２６に与えられる。スキャントリガ発生部２６は波長走査型光源１０の光の１走査の開始を示すスキャントリガ信号を発生するものである。又ｋトリガ発生部２７はスキャントリガ信号に基づき１走査の範囲内で、等周波数間隔で多数のｋトリガ信号を発生させるものである。このｋトリガ信号は信号解析部２５に入力される。
【００３３】
次に図２はｋトリガ発生部２７の構成を示す図である。ｋトリガ発生部２７は特許文献２に示されているように、クロック発生回路３１及びトリガ生成回路３２及びメモリ３３を有している。クロック発生回路３１は一定のタイミングでクロック信号を発生するものである。メモリ３３は読み書きを制御するＲＷ制御部３４によってデータを任意に書き換えることができるものとする。このメモリ３３は、波長走査型光源１０の発振周波数ｆ（ｆ１≦ｆ≦ｆ２）を時間の関数として次式
ｆ＝ｆ（ｔ）
で表すものとし、周波数の変化範囲δｆが等間隔になるタイミングをｔ1，ｔ2・・・ｔm・・・ｔn（ｍ＝１〜ｎ）とすると、次式で示される発振周波数
ｆ（ｔm）＝ｆ１＋（ｍ−１）δｆ
となるタイミング
ｔm＝ｆ-1｛ｆ１＋（ｍ−１）δｆ｝
をテーブルとして保持するものである。トリガ生成回路３２はトリガ信号が入力される毎にメモリ３３のデータを読み出すことによって等周波数間隔のタイミングのｋトリガ信号を発生するものである。
【００３４】
［信号解析部の構成］
信号解析部２５は、非線形方程式を解いてフーリエ係数を推定することにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現するものである。信号解析部２５は、例えば図３に示すように、入力インターフェイス（ＩＦ）５０と、各部を統括的に制御するＣＰＵ（Central Processing Unit）５１と、各種プログラムを含む各種情報を格納する読み取り専用のＲＯＭ（Read Only Memory）５２と、ワークエリアとして機能するＲＡＭ（Random Access Memory）５３と、各種情報を読み出し及び／又は書き込み可能に記憶する記憶部５４と、ユーザインターフェイスとしての図示しない所定の操作デバイスを介した入力操作の処理及び制御を行う入力操作制御部５５と、各種情報を表示する表示部５６とを備える。
【００３５】
入力インターフェイス部５０は、解析対象信号をＣＰＵ５１によって読み出し可能に入力させる機能を有する。なお、入力インターフェイスは、ｋトリガ信号を入力した場合には、増幅器２４の出力のＡ／Ｄ変換を行ってデジタル信号に変換する機能を有する。このとき、入力インターフェイスは、必要に応じてアンチエイリアシングフィルタを含むようにしてもよい。
【００３６】
ＣＰＵ５１は、記憶部５４等に格納されている各種アプリケーションプログラムをはじめとする各種プログラムを実行し、各部を統括的に制御する。
【００３７】
ＲＯＭ５２は、各種プログラムをはじめとする各種情報を格納している。またＲＡＭ５３は、ＣＰＵ５１が各種プログラムを実行する際のワークエリアとして機能する。ＣＰＵ５１はＲＡＭ５３に各種情報を一時記憶するとともに、ＲＯＭ５２，ＲＡＭ５３に記憶している各種情報を読み出すものである。
【００３８】
記憶部５４は、本発明にかかる信号解析プログラム等のアプリケーションプログラムの他、周波数解析の対象となる解析対象信号のデータをはじめとする各種情報を記憶する。この記憶部５４としては、例えば、ハードディスクや不揮発性メモリ等を用いることができる。また、記憶部５４には、本体に対して着脱可能とされるフレキシブルディスクやメモリカード等の記憶媒体に対して、各種情報の読み出し及び／又は書き込みを行うドライブ装置も含まれる。この記憶部５４に記憶されている各種情報は、ＣＰＵ５１の制御のもとに読み出される。
【００３９】
入力操作制御部５５は、例えば、キーボード、マウス、キーパッド、赤外線リモートコントローラ、スティックキー、又はプッシュボタンなどのユーザインターフェイスとなる所定の操作デバイスを介した入力操作を受け付け、操作内容を示す制御信号をＣＰＵ５１に供給するものである。
【００４０】
表示部５６は、例えば、液晶ディスプレイなどの各種表示デバイスであり、ＣＰＵ５１の制御のもとに各種情報を表示する。例えば、表示部５６は、ＣＰＵ５１によって信号解析プログラムが起動されると、その画面を表示し、入力された解析対象信号や周波数解析結果等を表示する。
【００４１】
［光コヒーレントトモグラフィの原理］
次に、波長走査型光源を用いた光コヒーレントトモグラフィの原理について説明する。光源から光周波数が連続的にかつ周期的に変化するコヒーレント光を対象物体に照射させ、マイケルソン、あるいはマッハツェンダなどの干渉光学計を用いて物体内部、あるいは生体表皮下層で反射した後方散乱光と参照光とを干渉させる。この干渉光の強度分布を計測し、光周波数の変化に対応した強度分布の変化を測定することによって、深さ方向に沿った断層画像を構築できる。さらに物体上で１次元、２次元に空間ビームを走査することによって、夫々２次元、３次元の断層画像を構築することができる。
【００４２】
干渉計において結合部１７から２つの腕の光路、すなわち参照ミラー１６までの光路Ｌ１と、物体中の反射面までの光路Ｌ２とが等しいときには、干渉光のビート周波数はゼロとなる。次に、反射光が物体内部のある深さｚから反射するとき、光周波数が時間的に変化していると、その光路差の分、物体からの反射光と参照ミラー１６からの反射光の周波数に差が生じ、干渉光にビートが生じる。ここで、例えば光源の光周波数が時間的に線形に走査されているとする。干渉計の腕の長さが等しい位置に物体の表面があり、物体の反射面は表面から深さｚの位置にのみあるとする。結合部１７での参照光の周波数と物体からの反射光（物体光）の周波数の時間的変化は、夫々図４の直線Ａ，Ｂのようになる。ここで光周波数は走査レートα［Ｈｚ／ｓ］で、時間Ｔ［ｓ］内で周波数幅Δｆ＝αＴ［Ｈｚ］にわたって走査される。参照光に対する物体光の遅れ時間τは、物体の屈折率をｎとすると、
τ＝２ｎｚ／ｃ
となる。従ってフォトダイオード２３で受光される干渉光は、ビート周波数
ｆｂ＝ατ＝（Δｆ／Ｔ）（２ｎｚ／ｃ）
で変動することになる。
【００４３】
実際は反射光は物体内部の深さに沿って連続的に異なった位置から発生するので、反射光はそれぞれの深さに対応した異なったビート周波数成分をもつ。従って干渉光の強度変化を周波数分析することによって、ビート周波数に対応するある特定の深さからの反射光強度を検出することができる。この反射強度の空間分布をとることで、断層画像を構築できる。
【００４４】
Ｉdct＝（ηｑ／ｈν）｛Ｐｒ＋Ｐｏ∫ｒ（ｚ）ｄｚ＋２（ＰｒＰｏ）1/2
∫ｒ（ｚ）cos（２ｋ（ｔ）ｚ＋φ）ｄｚ｝
干渉光信号Ｉdctの第１，２項はそれぞれ参照ミラーと、物体からの反射光の直流成分であり、第３項が干渉信号光成分であり、式（２）においてＰｒは参照光強度、Ｐｏはプローブ光強度ｒ（ｚ）は深さ方向の反射率分布を示す。Ｉdctの干渉光信号を周波数分析することによって、物体中の任意の深さに対応する散乱光強度の関係を得ることができる。本発明の信号解析部２５は、この式で示される干渉光信号Ｉdctを解析対象信号Ｘ（ｎ）とし、ノンハーモニック分析によって周波数解析するようにしたものである。
【００４５】
［光断層画像表示システムの動作］
次に本実施の形態の動作について説明する。波長走査型光源１０を駆動する。これによって前述したように光ファイバ１４，１８を介して信号光が参照ミラー及び物体にまで照射され、その反射光が結合部１７を介して得られる。そのビート周波数を持つ信号がフォトダイオード２３に得られる。これを増幅した受光信号が信号解析部２５に入力される。又前述したようにｋトリガ発生部２７から等周波数間隔のｋトリガ信号が得られ、このｋトリガ信号に応じてＡ／Ｄ変換を行い、信号解析部２５で受光信号を周波数分析することにより深さ方向に１次元の断層画像が得られる。そしてスキャニングミラー２０を回動させることによって光の入射位置を変化させ、これによって2次元の断面画像を得ることができる。又、この干渉計自体又は測定対象をスキャニングミラー２０による光の走査方向と垂直に移動させることにより、３次元断面画像を得ることができる。
【００４６】
次に信号解析部２５の処理について更に詳細に説明する。ＣＰＵ５１の制御のもとに、信号解析プログラムを実行して入力された信号の周波数解析を行う。なお、周波数解析の対象となる解析対象信号Ｘ（ｎ）は、信号入力インターフェイス５０を介して入力され、ＣＰＵ５１によって読み出し可能に記憶部５４に記憶される。信号解析部２５は、ＣＰＵ５１の制御のもとに、このようにして入力された解析対象信号の周波数解析を行い、その解析結果等を表示部５６に表示出力する。
【００４７】
［周波数解析アルゴリズム］
以下、このような信号解析部２５における周波数解析アルゴリズムについて詳述する。本発明の信号解析部２５にて新たに提案する周波数解析手法（以下、本周波数解析手法という。）は、「周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にある」という従来の周波数解析手法の問題を本質的に解決するために、非線形方程式の解を求める問題に着目した。すなわち、本周波数解析手法においては、次式（３）に示す非周期信号のフーリエ変換式の周波数パラメータを求める問題を非線形方程式の最適解を求める問題に置き換えたものである。
【００４８】
【数３】
【００４９】
具体的には、本周波数解析手法においては、次式（４）に示すように、解析対象信号ｘ（ｎ）と正弦波モデル信号との差の二乗和で表される非線形方程式の最適解として、この非線形方程式の右辺が最小値になるような周波数ｆ’、振幅Ａ’、及び初期位相φ’を求める。なお、次式（４）において、Ｌはフレーム長（解析窓長）であり、ｆsはサンプリング周波数［Ｈｚ］である。このように、非線形方程式に周波数ｆ’をパラメータとして盛り込んだことは、いまだ事例がなく斬新である。換言すれば、本周波数解析手法においては、振幅Ａ’及び初期位相φ’のみならず、周波数ｆ’をも正確に求めることを目的とし、上式（３）に対して非線形方程式の解法を用いるものである。本周波数解析手法においては、このような最小二乗法によって非線形方程式の最適解を求める問題に帰着させることにより、解析窓の影響やエイリアシングの影響がなくなり、解析窓長が、１周期未満であってもよく、周期の整数倍でなくてもよく、さらには、不等間隔であってもよい等、柔軟な周波数解析処理を実現することが可能となる。
【００５０】
【数４】
【００５１】
なお、上式（４）は、解析対象信号ｘ（ｎ）が、受光信号、時系列信号等の１次元信号である場合を想定しているが、本周波数解析手法自体は、２次元の画像データ等の２次元信号や動画像データ等の３次元信号、ホログラムのような立体動画像データ等の４次元信号、さらにはそれ以上の次元の信号に対しても適用可能である。すなわち、本周波数解析手法は、１次元の解析対象信号をｓ（ｘn）、２次元の解析対象信号をｓ（ｘm，ｙn）、３次元の解析対象信号をｓ（ｘl，ｙm，ｚn）、４次元の解析対象信号をｓ（ｘk，ｙl，ｚm，ｗn）として一般化すると、解析対象信号が１次元、２次元、３次元、又は４次元の場合のそれぞれについて、次式（５）乃至次式（８）に示す非線形方程式の最適解を求めることになる。
【００５２】
【数５】
【数６】
【数７】
【数８】
【００５３】
換言すれば、本周波数解析手法は、任意のｎ次元（ｎは１以上の整数）の解析対象信号をｓ（ｘ１n1，ｘ２n2，ｘ３n3，・・・，ｘｎnn）として一般化すると、次式（９）に示す非線形方程式の最適解を求めることになる。
【数９】
【００５４】
なお、以下では、説明の便宜上、解析対象信号が１次元であり、上式（４）に示す非線形方程式の最適解を求めるものとする。
【００５５】
さて、上式（４）に示す非線形方程式の最適解を実際に求めるにあたっては、以下のような方法をとることができる。
【００５６】
本周波数解析手法においては、振幅Ａ’、周波数ｆ’、及び初期位相φ’のそれぞれについて適切な初期値を求め、これら初期値から非線形方程式の解法を用いて最適解に収束させる。この非線形問題では、上式（４）をコスト関数とする最小化問題とする。なお、適切な初期値は、離散フーリエ変換（Discrete Fourier Transform；以下、ＤＦＴという。）やウェーブレット変換等の任意の周波数変換を行ったり、フィルタリングを行うことによっておおよその見当をつけたりする等、既存の任意の方法を適用して求めることができる。
【００５７】
まず、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’について、いわゆる最急降下法を適用し、周波数パラメータｆm’，φm’を次式（１０）及び次式（１１）によって求める。
【００５８】
【数１０】
【数１１】
【００５９】
なお、上式（１０）及び上式（１１）においては、次式（１２）と略している。また、μmは、いわゆる減速法に基づく重み係数であり、各漸化式によって求められるコスト関数を単調減少数列にするために、適時０〜１の値をとる。
【００６０】
【数１２】
【００６１】
周波数パラメータｆm’，φm’を求めることができれば、上式（４）における正弦波モデル信号の係数としての周波数パラメータＡ’を一意に求めることができるため、本周波数解析手法においては、次式（１３）によって周波数パラメータＡm’を収束させる。
【００６２】
【数１３】
【００６３】
本周波数解析手法においては、これら一連の計算を反復して行うことにより、振幅Ａ’、周波数ｆ’、及び初期位相φ’を高精度に収束させることができる。特に、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’と、係数としての周波数パラメータＡ’とを別個に求めることにより、計算を簡便に行うことができる。
【００６４】
しかしながら、最急降下法は、比較的広い範囲から収束するものの、１回の反復では精度が低く、収束するまでに時間を要する。
【００６５】
そこで、本周波数解析手法においては、最急降下法を適用して周波数パラメータｆm’，φm’をある程度まで収束させた後、さらに、いわゆるニュートン法を適用して高精度に収束させるのが望ましい。具体的には、本周波数解析手法においては、ニュートン法として、次式（１４）及び次式（１５）に示す漸化式によって周波数パラメータｆm’，φm’を求める。
【００６６】
【数１４】
【数１５】
【００６７】
ただし、上式（１４）及び上式（１５）において、Ｊは次式（１６）とし、次式（１７）と略している。また、νmもμmと同様に減速法に基づく重み係数であり、適時０〜１の値をとる。
【００６８】
【数１６】
【数１７】
【００６９】
本周波数解析手法においては、上式（１４）及び上式（１５）によって周波数パラメータｆm’，φm’を求めた後、最急降下法と同様に、上式（１３）によって周波数パラメータＡm’を収束させ、この一連の計算をさらに反復して行う。
【００７０】
このように、本周波数解析手法においては、最急降下法とニュートン法とを組み合わせたハイブリッド型の解法を用いることにより、高速に且つ高精度に周波数パラメータＡ’，ｆ’，φ’を推定することができる。
【００７１】
また、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次減算処理することにより、近似的にスペクトルパラメータを導出することができる。ここで、解析対象信号ｘ（ｎ）が複数の正弦波の和であり、次式（１８）のように表されているとする。
【００７２】
【数１８】
【００７３】
パーセヴァル（Parseval）の定理より、解析対象信号ｘ（ｎ）の周波数ｆkと正弦波モデル信号の周波数パラメータｆ’とが全く一致しない場合、すなわち、次式（１９）である場合には、上式（４）に示す非線形方程式は次式（２０）となる。また、周波数パラメータｆ’，φ’の組が、周波数ｆk及び初期位相φkの組のいずれかに一致する場合には、上式（４）に示す非線形方程式は次式（２１）となる。さらに、振幅Ａjが周波数パラメータＡ’とも一致した場合には、解析対象信号から推定スペクトルに関する周波数成分を完全に消去することができる。そのため、最適解を求める問題は、周波数に対して独立であり、解析対象信号から順次個別に推定すれば、複数の正弦波で表される信号にも応用することができる。
【００７４】
【数１９】
【数２０】
【数２１】
【００７５】
すなわち、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次残差信号に対して同様に処理を行い、複数の正弦波を抽出することができる。
【００７６】
ただし、複数のスペクトルの場合には、上式（３）からもわかるように、無限長を仮定しているため、複数のスペクトルの周波数同士が近い場合には、合成スペクトルの周期が解析窓長よりも数段長くなり、上式（１７）が満たされなくなるため、誤差が発生する場合もある。
【００７７】
従来のＯＣＴでは、受光信号を複合正弦波によって表現するために多くのスペクトル数（正弦波の数）が必要であったが、本周波数解析手法においては、そのような信号であっても僅かなスペクトル数で誤差なく表現することができる。すなわち、信号をより少ないスペクトル数で表現可能であることは、情報圧縮の用途に有効であることは勿論のこと、それ以外にも、物理現象に関して波数表現をするようなアプローチに用いることにより、信号の本質的な特性も見抜くことができることを示している。
【００７８】
［本周波数解析手法の有効性］
以下、本周波数解析手法の有効性について具体的に説明する。
【００７９】
従来、センサによって計測された信号を、計算機を用いて解析するためには、一般に、例えば図５に示すように、その信号の一部を切り出し、ＤＦＴといった代表的な調和型の周波数解析手法によって解析していた。
【００８０】
しかしながら、解析対象信号を切り出して（窓掛けして）調和型の解析を行うことは、当該解析対象信号が、切り出した一部の信号が完全周期で繰り返される信号であると擬制して解析することに他ならず、当然ながらその周波数特性も、元のスペクトル特性とは完全に一致しないものとなる。また、この信号切り出しの影響を回避するために、窓関数によって影響を軽減することも試みられているが、切り出した一部の信号を完全周期として解析することには変わりなく、本質的な問題解決には至っていないのが現状である。
【００８１】
これに対して、本周波数解析手法は、非周期信号を想定した周波数解析を目的として定式化を行っていることから、解析窓長と周波数分解能の制約を受けることなく、正確に周波数とそのパラメータとを推定することができる。
【００８２】
また、本周波数解析手法は、周波数の推定についても、従来の離散的な探索から、動的且つ連続的な探索となるため、格段に精度を向上させることができ、革新的な手法である。
【００８３】
さらに、多次元化に関しても、従来の周波数解析手法においては、解析対象信号が１次元である場合と同様に窓関数の影響を受ける。具体的には、例えば図６に示すように、２次元の解析対象信号（原画像データ）のうち図６左図中の破線部を切り出してＤＦＴを施して解析したとしても、その切り出した区間を完全周期とした画像が延々と繰り返されるものを解析したことと等価であることから、右上図に示すように、復元した画像データが窓関数の影響を受けたものとなり、本質的な特性を見出すことは困難である。
【００８４】
これに対して、上式（６）を適用した本周波数解析手法においては、図６右下図に示すように、窓関数の影響が軽減されている。そのため本質的な周波数特性を見出すことができ、切り出した区間の周囲の画像も完全に再現することができる。
【００８５】
以上のように、本周波数解析手法は、非線形方程式の最適解を求めることにより、正弦波モデル信号の周波数ｆ’、振幅Ａ’、及び初期位相φ’を高速に且つ高精度に求めることができる。具体的な精度を立証するために、本願発明者は、ＤＦＴと、ＤＦＴの発展型のうち最も解析精度が高いといわれているＧＨＡ（Generalized Harmonic Analysis）とを比較対象として精度の検証を行った。
【００８６】
なお、ＤＦＴやＧＨＡは、１つの解析窓長に見かけ上複数の窓長を持たせていることから、周波数分解能が解析窓長に依存するが、その分解周波数が有限長であり、解析対象信号の周波数が分解周波数以外の周波数となった場合には解析することができず、解析対象信号が正確に解析できる周波数と異なる場合には、最も近い分解周波数の他に、その周辺に小さなスペクトルの周波数（側帯波成分）が現れ、複数の周波数が出現してしまう。
【００８７】
このような現象が本周波数解析手法においても生じるか否かについて、すなわち、本周波数解析手法の周波数分解能を検証するために、解析窓長を１秒（１０２４サンプル）とした１次元の非常に短い単一正弦波を解析し、各手法によって正弦波を１本抽出して元の信号との二乗誤差を調べた。その結果を図７に示す。
【００８８】
図７に示すように、ＤＦＴにおいては、基本周波数の整数倍以外の周波数における解析精度の悪化がみられた。また、ＧＨＡにおいては、１Ｈｚ以上の周波数ではＤＦＴと比べて２〜５桁程度の精度向上がみられた。これに対して、本周波数解析手法においては、１Ｈｚ以上の周波数ではＤＦＴと比べて１０桁以上、ＧＨＡと比べて５桁以上の精度向上がみられた。すなわち、本周波数解析手法は、既存の周波数解析手法と比べて１０万〜１００億倍以上の精度向上がみられた。
【００８９】
特に、１Ｈｚ以下の周波数を正確に推定することができるということは、解析窓長を超えた長い周期信号であっても解析可能であることを示しており、そのスペクトル構造をより正確に推定可能であることを示している。
【００９０】
このように、本周波数解析手法は、最も解析精度が高いといわれているＧＨＡと比べても驚くべき高精度に解析を行うことができる。国内外の様々な研究事例をみても、これほどの精度で信号を解析できる手法は存在せず、本周波数解析手法は、今後、正確な解析が必要とされる様々な分野への応用が期待できる手法であるといえる。
【００９１】
［画像の符号化］
本周波数解析手法においては、窓の影響が少ないため、スペクトルに側帯波成分が全く出ない。すなわち、本周波数解析手法は、上述したように、複雑な解析対象信号であっても、僅かなスペクトル数で効率的に誤差なく表現することができる。
【００９２】
図８に２次元信号をＤＦＴと本周波数解析手法とのそれぞれによって解析した具体例を示す。図８左図に示すように、本来であれば１本のスペクトルで表現される信号をＤＦＴによって変換した場合には、一般に、変換して得られるスペクトルも、解析窓の影響に起因して、図８右上図に示すように、複数の側帯波成分が群れて現れてしまう。これは、画像の圧縮符号化において広く用いられている離散コサイン変換（Discrete Cosine Transform；ＤＣＴ）等においても同様である。これに対して、本周波数解析手法は、図８右下図に示すように、元の信号のスペクトルを正確に推定することができる。正確に解析できると、画像データを表現するスペクトル数を格段に削減できる可能性があり、画質を落とさずに飛躍的な情報圧縮効果を期待することができる。
【００９３】
又この実施の形態では波長走査型光源の走査範囲をあまり広くすることなく高分解能の画像が得られる。そのため従来では使用できなかった安価な波長走査型光源を用いて必要な分解能の断層画像を得ることができる。又広範囲の波長を走査することができる光源を用いれば、より精密な断面画像とすることができる。
【００９４】
尚この実施の形態ではｋトリガ発生部により波長の１スキャニングの間に等周波数間隔となるタイミングで多数のｋトリガ信号を発生するようにしているが、波長の１走査の開示のみのスキャントリガ信号としてもよい。この場合に特開２００８−２６１７７８号に示されるように、物体の測定位置に反射物を配置したときに受光素子からの出力をデータ番号を付して等時間間隔で取り込み、ビート信号のピークとなるデータ番号を算出し、ピーク番号に対してデータ番号を示す近似曲線を算出し、ピーク番号の変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出し、データ取り込み数列とし、測定物体を配置したときの受光素子からの等時間間隔からの出力からデータ数列の各要素でのタイミングをトリガ信号として信号解析部に与えるようにしてもよい。
【００９５】
（第２の実施の形態）
図９は本発明の第２の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。本実施の形態はＳＤ−ＯＣＴによる光断層画像表示システムであり、光源として所定の波長範囲に一様な光を出射する光源４０を用いる。この光源４０からの光は光増幅器１３を介して光干渉計に導かれることは第１の実施の形態と同様である。この実施の形態ではレンズ２２に得られる光を回折格子４１に入射する。回折格子４１は入射光の波長に応じて光を分散させる波長分光手段であり、分散した光は一次元ＣＣＤ４２に入力される。一次元ＣＣＤ４２は所定の周期で各位置の光信号を一次元の受光信号に変換するものであり、その出力は増幅器４３を介して信号解析部２５に伝えられる。この場合には、ＣＣＤ４２の一次元の各画素からの出力のタイミングに周期してトリガ信号を入力インターフェイスに伝える。これ以降の構成については第１の実施の形態と同様である。このようにして１次元ＣＣＤで得られた受光信号に基づいて周波数を解析するようにしても、同様に断層画像を得ることができる。
【００９６】
又前述した実施の形態における信号解析装置ではソフトウェアによる周波数解析を行うものとして説明したが、本発明は、本周波数解析手法のアルゴリズムを実装したＤＳＰ（Digital Signal Processor）等、積和演算を行うことが可能であればハードウェアによっても実現することができる。
【００９７】
このように、本発明は、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【技術分野】
【０００１】
本発明は、信号処理に特徴を有する光断層画像表示システムに関するものである。
【背景技術】
【０００２】
近年内視鏡治療などの医療技術の進歩に伴って、病理組織の診断を非深襲かつリアルタイムに行う診断方法が望まれている。従来例えばＣＣＤを用いた電子内視鏡や、ＣＴ、ＭＲＩ、超音波による画像化が診断方法として用いられている。電子内視鏡は生体の表面の観察に限定され、また後者の画像診断システムはミクロンオーダーの分解能で観察するには技術的な限界があった。このような方法を補完する技術として、光コヒーレンストモグラフィーシステム（ＯＣＴ）が注目されている。
【０００３】
ＯＣＴの中には、時間領域ＯＣＴ（ＴＤ−ＯＣＴ）と周波数領域ＯＣＴ（ＦＤ−ＯＣＴ）の２種類があり、またＦＤ−ＯＣＴの中にもスペクトロメータタイプ（ＳＤ−ＯＣＴ）と波長走査型光源タイプ（ＳＳ−ＯＣＴ）の２つがある。従来の時間領域ＯＣＴの場合には、広帯域な光を当て生体内部からの反射光の干渉成分を周波数分析していたが、この方法だと干渉光の中に異なる深さからの反射光も重なりあうために、ある特定の深さからの信号光だけを感度良く検出できなかった。
【０００４】
一方スペクトロメータタイプのＯＣＴは、特許文献１に示されているように、白色光源を生体に光を照射し、参照光と生体内の異なる深さから戻ってくる反射光とを干渉計で干渉させる。そして反射光を回折格子等で分光して１次元ＣＣＤで受光し、電気信号とする。波長走査型ＯＣＴは、特許文献２に記されているように、生体に光を照射し、照射光の波長を連続的に変化させ、参照光と生体内の異なる深さから戻ってくる反射光とを干渉計で干渉させる。これらはいずれも干渉信号の周波数成分を分析することによって、断層画像を得るシステムである。ＦＤ−ＯＣＴは得られた干渉信号光をフーリエ変換し、断層画像を得るようにしている。この技術は高分解能の断層画像を構築することができるため、高度なシステムとして期待されており、理論的に時間領域ＯＣＴの１００倍以上の感度が得られる。波長走査型ＯＣＴは測定感度も高く、動的ノイズに強いという点で内視鏡などの実使用に好適である。ここで照射する光の波長走査の帯域が広いほど周波数分析の帯域が上がるので、深さ方向の分解能が上がる。
【０００５】
しかしながらＦＤ−ＯＣＴはフーリエ変換時に窓の影響によりメインローブ以外にもサイドローブが発生するため、分解能に限界があるという欠点があった。以下この問題点について説明する。
【０００６】
従来から、フーリエ変換を用いた正弦波モデルによる信号解析方法としては、非特許文献１にも紹介されているように、様々な手法が提案されているが、その多くは、信号の周波数を事前に推定し、その推定した周波数における最適な振幅や初期位相を求めることによって行われている。
【０００７】
信号は複数の周波数の正弦波の重ね合わせであり、信号を解析するには、当該信号を構成する様々な周波数の正弦波毎に振幅や初期位相等の周波数パラメータを推定することになる。このような周波数解析手法としては、フーリエ変換が広く知られている。フーリエ変換における周波数パラメータの推定は、周期的な信号の場合には、次式（１）に示すように、有限長の窓によるフーリエ変換式を用いて行うことができ、この次式（１）に基づいて、特定の周波数ｆ（＝ｎ／Ｔ：ｎは整数、Ｔは解析窓長）のパラメータ（振幅Ａ及び初期位相φ）を推定することができる。ただし、振幅Ａと初期位相φとＸ（ｆ）との関係は次式（２）である。
【０００８】
【数１】
【数２】
【０００９】
このようなフーリエ変換によって推定することが可能な周波数パラメータは、解析窓長Ｔに依存し、１／Ｔの整数倍の周波数分しか計算することができない。したがって、上式（１）を用いてパラメータを推定している限り、周波数ｆは、原理的に離散的となり、非周期信号の周波数成分を特定することはできない。また、フーリエ変換を用いて周波数解析を行っている分野では、そのほとんどが、周期を持たない非周期信号を解析対象とするにもかかわらず、無限長の積分区間にわたる積分計算が極めて困難であることから、上式（１）を用いなければならないのが現状である。
【００１０】
この問題を解決するために、窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みがある。
【００１１】
例えば、非特許文献２に記載の短時間フーリエ変換（Short-Time Fourier Transform、Short-Term Fourier Transform；ＳＴＦＴ）は、音響信号を短時間の解析窓によって複数の信号に分割し、各信号に対して独立にフーリエ変換を行う手法である。
【００１２】
また、非特許文献３乃至非特許文献６に記載のＡＢＳ（Analysis by Synthesis）法やＧＨＡ（Generalized Harmonic Analysis）等においては、最適な周波数を定めるための周波数を複数用意して振幅及び初期位相を求め、その中から最適な周波数と振幅と初期位相との組み合わせを選択している。
【先行技術文献】
【特許文献】
【００１３】
【特許文献１】特開２００６−１３２９９６号
【特許文献２】特開２００７−０２４６７７号
【非特許文献】
【００１４】
【非特許文献１】東山三樹夫、小池恒彦、「高い周波数分析精度の信号分析手法」、日本音響学会論文誌、１９８８年、第５４巻、第８号
【非特許文献２】L. R. Rabiner and R. W. Schafer, "Digital Processing of Speech Signals", Prentice-Hall, Englewood Cliffs, NJ, 1978年
【非特許文献３】T. F. Quatieri and R. J. McAulay, "Speech transformations based on a sinesoidal representation", IEEE Trans. Acoust. Speech Signal Process. ASSP 34, 1986年, p.1449-1464
【非特許文献４】E. B. George and M. J. Smith, "Analysis-by-synthesis/overlap-add sinusoidal representation", J. Audio Eng. Soc. 40, 1992年, p.497-515
【非特許文献５】E. B. George and M. J. Smith, "Speech analysis/synthesis and modification using an analysis-by-synthesis/overlap-add sinusoidal model", IEEE Trans. Speech Audio Process. 5, 1997年, p.389-406
【非特許文献６】T. Terada, H. Nakajima, M. Tohyama and Y. Hirata, "Non-stationary waveform analysis and synthesis using generalized harmonic analysis", IEEE-SP, Int. Symp. TF/TS Analysis, 1994年, p.429-432
【発明の開示】
【発明が解決しようとする課題】
【００１５】
しかしながら、上述した非特許文献２に記載の短時間フーリエ変換においては、短時間の解析窓長の逆数の整数倍の周波数でしか、周波数パラメータを推定することができないという問題がなおも存在する。
【００１６】
また、上述した非特許文献３乃至非特許文献６に記載のＡＢＳ法やＧＨＡ等においては、元の信号に対して、予め用意された周波数と異なる周波数が存在する場合には、誤った組み合わせが複数検波されてしまい、解析精度が低下するという問題がある。
【００１７】
これら窓関数の調整や解析区間を複数利用して見かけ上の周波数分解能を向上させる試みは、
１）周波数分解能を上げるためには解析窓長を長くしなければならないが、時間的な周波数変動をともなう信号では周波数が平均化され、時間的な変動を捉えることが困難となること
２）一方、解析窓長を短くすると、周波数分解能が低下すること
という問題が指摘されている。
【００１８】
このように、周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にあり、信号解析上の大きな問題となっていた。
【００１９】
本発明はこのような実状に鑑みてなされたものであり、光断層画像表示システムにおいて周波数分解能が解析窓長に依存せずに高い周波数分解能を有し、高精度に周波数解析を行うことができ、高分解能、高感度、高速で画像表示することができる光断層画像表示システムを提供することを目的とする。
【課題を解決するための手段】
【００２０】
この課題を解決するために、本発明の光断層画像表示システムは、周期的に光の発振波長を走査する波長走査型光源と、前記波長走査型光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、前記干渉光学計より得られる干渉光を受光し、ビート信号を得る受光素子と、前記受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、前記信号解析部は、前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、前記演算手段より得られる断面画像を表示する表示部と、を具備するものである。
【００２１】
この課題を解決するために、本発明の光断層画像表示システムは、所定の範囲で一様な波長の光を発振する光源と、前記光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、前記干渉光学計から得られる干渉光を波長に基づいて分岐する波長分光手段と、前記波長分光手段より分光された範囲の干渉光を受光して電気信号に変換する一次元受光素子と、前記一次元受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、前記信号解析部は、前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、前記演算手段より得られる断面画像を表示する表示部と、を具備するものである。
【００２２】
ここで前記演算手段は、前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’のそれぞれについて適切な初期値を求め、求めた初期値から非線形方程式の最適解に収束させて前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を求めるようにしてもよい。
【００２３】
ここで前記演算手段は、前記正弦波モデル信号の位相を構成する前記周波数ｆ’及び前記初期位相φ’について最急降下法を適用して収束させ、当該周波数ｆ’及び当該初期位相φ’を求め、前記正弦波モデル信号の係数としての前記振幅Ａ’について最急降下法を適用して収束させ、当該振幅Ａ’を求めるようにしてもよい。
【００２４】
ここで前記演算手段は、前記最急降下法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させた後、さらに、ニュートン法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させるようにしてもよい。
【００２５】
ここで前記波長走査型光源の１走査の期間内に、前記波長走査型光源の光の等周波数間隔でのトリガ信号を発生するｋトリガ発生部を更に有し、前記信号解析部の記憶部は、前記ｋトリガ発生部により発生したｋトリガ信号のタイミングで前記干渉信号を保持するものとしてもよい。
【００２６】
ここで前記信号解析部の演算手段は、物体の測定位置に反射物を配置したときに、前記受光素子からの出力をデータ番号を付して等時間間隔で取り込み、前記ビート信号のピークとなるデータ番号を算出し、ピーク番号に対するデータ番号を示す近似曲線を算出し、前記ピーク番号変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出しデータ取り込み数列とし、測定位置に測定物体を配置したときの前記受光素子からの等時間間隔の出力から、前記データ取り込み数列の各要素でのタイミングの信号レベルを算出し、得られた等周波数間隔のデータを前記記憶手段に保持するものとしてもよい。
【発明の効果】
【００２７】
このような本発明による光断層画像表示システムによれば、干渉光信号を周波数分析する際に非周期信号を想定した周波数解析を目的として定式化を行い、非周期信号のフーリエ変換式のパラメータを求める問題を非線形方程式の最適解を求める問題に置き換えることにより、周波数分解能が解析窓長に依存しない。従って従来のフーリエ変換による周波数解析手法における精度と比べて１０万〜１００億倍以上という驚くべき高精度で周波数解析を行うことができ、高分解能で断層画像を表示することができるという効果が得られる。
【図面の簡単な説明】
【００２８】
【図１】図１は本発明の第１の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。
【図２】図２は本発明の実施の形態による光断層画像表示システムのｋトリガ発生部の構成を示すブロック図である。
【図３】図３は本発明の実施の形態による光断層画像表示システムの信号解析部の構成を示すブロック図である。
【図４】図４は走査角度と発振周波数の関係の一例を示すグラフである。
【図５】図５は本周波数解析手法とＤＦＴとの違いを説明するための図である。
【図６】図６は本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用した場合の具体例を示す図である。
【図７】図７は本周波数解析手法とＤＦＴとＧＨＡとの違いを説明するための図であり、各手法の誤差を求めた結果を示す図である。
【図８】図８は本周波数解析手法とＤＦＴとの違いを説明するための図であり、２次元の解析対象信号に各手法を適用して得られるスペクトルの解析精度の具体例を示す図である。
【図９】図９は本発明の第２の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。
【符号の説明】
【００２９】
１０ 波長走査型光源
１１，１４，１８ 光ファイバ
１５，１９，２２ レンズ
１６ 参照ミラー
２０ スキャニングミラー
２３ フォトダイオード
２４ 増幅器
２５ 信号解析部
２６ スキャントリガ発生部
２７ ｋトリガ発生部
３１ クロック発生回路
３２ トリガ生成回路
３３ メモリ
３４ ＲＷ制御部
４０ 光源
４１ 回折格子
４２ 一次元ＣＣＤ
４３ 増幅器
５０ 入力インターフェイス
５１ ＣＰＵ
５２ ＲＯＭ
５３ ＲＡＭ
５４ 記憶部
５５ 入力操作制御部
５６ 表示部
【発明を実施するための形態】
【００３０】
（第１の実施の形態）
［光断層画像表示システムの構成］
図１は本発明の第１の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。本実施の形態はＳＳ−ＯＣＴによる光断層画像表示システムであり、光源として波長走査型光源を用いる。波長走査型光源１０には一定の周波数範囲の光信号を発振する波長可変型のレーザを用いる。このレーザの出力は光ファイバ１１を介して分岐部１２に与えられる。光ファイバ１１の他端には光増幅器１３が設けられる。光増幅器１３は波長走査型光源１０のレーザ光をそのまま増幅するものであり、その出力は光ファイバ１４に与えられる。光ファイバ１４の他端にはコリメートレンズ１５及び参照ミラー１６を設ける。又この光ファイバ１４の中間部分には、他の光ファイバ１８を接近させて干渉させる結合部１７が設けられる。光ファイバ１８の一端には、波長走査型光源１０から結合部１７を介して得られた光信号を平行光とするコリメートレンズ１９、光をスキャニングするスキャニングミラー２０が設けられる。スキャニングミラー２０は紙面に垂直な軸を中心にして一定範囲で回動することによって平行光の反射角度を変化させるものである。集束レンズ２１はこの反射光を受光する位置に配置し、測定部位へ光を集束すると共に水平方向にスキャニング（走査）する。ここで結合部１７から参照ミラー１６までの光学距離Ｌ１と、結合部１７から測定部位の表面までの光学距離Ｌ２とを等しくしておく。さて光ファイバ１８の他端にはレンズ２２を介してフォトダイオード２３を接続する。フォトダイオード２３は、参照ミラー１６からの反射光と測定部位で反射された光の干渉光を受光することによって、そのビート信号を電気信号として得る受光素子である。ここで光ファイバ１４，１８と結合部１７、コリメートレンズ１５、参照ミラー１６、コリメートレンズ１９、スキャニングミラー２０、集束レンズ２１は干渉光学計を構成している。
【００３１】
さてフォトダイオード２３の出力は増幅器２４を介して信号解析部２５に入力される。信号解析部２５は後述するように干渉計から得られる受光信号を周波数領域の信号に変換することによって、断層画像信号を得るものである。
【００３２】
又波長走査型光源１０の出力の一部は分岐部１２によって分岐されてスキャントリガ発生部２６に与えられる。スキャントリガ発生部２６は波長走査型光源１０の光の１走査の開始を示すスキャントリガ信号を発生するものである。又ｋトリガ発生部２７はスキャントリガ信号に基づき１走査の範囲内で、等周波数間隔で多数のｋトリガ信号を発生させるものである。このｋトリガ信号は信号解析部２５に入力される。
【００３３】
次に図２はｋトリガ発生部２７の構成を示す図である。ｋトリガ発生部２７は特許文献２に示されているように、クロック発生回路３１及びトリガ生成回路３２及びメモリ３３を有している。クロック発生回路３１は一定のタイミングでクロック信号を発生するものである。メモリ３３は読み書きを制御するＲＷ制御部３４によってデータを任意に書き換えることができるものとする。このメモリ３３は、波長走査型光源１０の発振周波数ｆ（ｆ１≦ｆ≦ｆ２）を時間の関数として次式
ｆ＝ｆ（ｔ）
で表すものとし、周波数の変化範囲δｆが等間隔になるタイミングをｔ1，ｔ2・・・ｔm・・・ｔn（ｍ＝１〜ｎ）とすると、次式で示される発振周波数
ｆ（ｔm）＝ｆ１＋（ｍ−１）δｆ
となるタイミング
ｔm＝ｆ-1｛ｆ１＋（ｍ−１）δｆ｝
をテーブルとして保持するものである。トリガ生成回路３２はトリガ信号が入力される毎にメモリ３３のデータを読み出すことによって等周波数間隔のタイミングのｋトリガ信号を発生するものである。
【００３４】
［信号解析部の構成］
信号解析部２５は、非線形方程式を解いてフーリエ係数を推定することにより、周波数分解能が解析窓長に依存しない新たな周波数解析手法を実現するものである。信号解析部２５は、例えば図３に示すように、入力インターフェイス（ＩＦ）５０と、各部を統括的に制御するＣＰＵ（Central Processing Unit）５１と、各種プログラムを含む各種情報を格納する読み取り専用のＲＯＭ（Read Only Memory）５２と、ワークエリアとして機能するＲＡＭ（Random Access Memory）５３と、各種情報を読み出し及び／又は書き込み可能に記憶する記憶部５４と、ユーザインターフェイスとしての図示しない所定の操作デバイスを介した入力操作の処理及び制御を行う入力操作制御部５５と、各種情報を表示する表示部５６とを備える。
【００３５】
入力インターフェイス部５０は、解析対象信号をＣＰＵ５１によって読み出し可能に入力させる機能を有する。なお、入力インターフェイスは、ｋトリガ信号を入力した場合には、増幅器２４の出力のＡ／Ｄ変換を行ってデジタル信号に変換する機能を有する。このとき、入力インターフェイスは、必要に応じてアンチエイリアシングフィルタを含むようにしてもよい。
【００３６】
ＣＰＵ５１は、記憶部５４等に格納されている各種アプリケーションプログラムをはじめとする各種プログラムを実行し、各部を統括的に制御する。
【００３７】
ＲＯＭ５２は、各種プログラムをはじめとする各種情報を格納している。またＲＡＭ５３は、ＣＰＵ５１が各種プログラムを実行する際のワークエリアとして機能する。ＣＰＵ５１はＲＡＭ５３に各種情報を一時記憶するとともに、ＲＯＭ５２，ＲＡＭ５３に記憶している各種情報を読み出すものである。
【００３８】
記憶部５４は、本発明にかかる信号解析プログラム等のアプリケーションプログラムの他、周波数解析の対象となる解析対象信号のデータをはじめとする各種情報を記憶する。この記憶部５４としては、例えば、ハードディスクや不揮発性メモリ等を用いることができる。また、記憶部５４には、本体に対して着脱可能とされるフレキシブルディスクやメモリカード等の記憶媒体に対して、各種情報の読み出し及び／又は書き込みを行うドライブ装置も含まれる。この記憶部５４に記憶されている各種情報は、ＣＰＵ５１の制御のもとに読み出される。
【００３９】
入力操作制御部５５は、例えば、キーボード、マウス、キーパッド、赤外線リモートコントローラ、スティックキー、又はプッシュボタンなどのユーザインターフェイスとなる所定の操作デバイスを介した入力操作を受け付け、操作内容を示す制御信号をＣＰＵ５１に供給するものである。
【００４０】
表示部５６は、例えば、液晶ディスプレイなどの各種表示デバイスであり、ＣＰＵ５１の制御のもとに各種情報を表示する。例えば、表示部５６は、ＣＰＵ５１によって信号解析プログラムが起動されると、その画面を表示し、入力された解析対象信号や周波数解析結果等を表示する。
【００４１】
［光コヒーレントトモグラフィの原理］
次に、波長走査型光源を用いた光コヒーレントトモグラフィの原理について説明する。光源から光周波数が連続的にかつ周期的に変化するコヒーレント光を対象物体に照射させ、マイケルソン、あるいはマッハツェンダなどの干渉光学計を用いて物体内部、あるいは生体表皮下層で反射した後方散乱光と参照光とを干渉させる。この干渉光の強度分布を計測し、光周波数の変化に対応した強度分布の変化を測定することによって、深さ方向に沿った断層画像を構築できる。さらに物体上で１次元、２次元に空間ビームを走査することによって、夫々２次元、３次元の断層画像を構築することができる。
【００４２】
干渉計において結合部１７から２つの腕の光路、すなわち参照ミラー１６までの光路Ｌ１と、物体中の反射面までの光路Ｌ２とが等しいときには、干渉光のビート周波数はゼロとなる。次に、反射光が物体内部のある深さｚから反射するとき、光周波数が時間的に変化していると、その光路差の分、物体からの反射光と参照ミラー１６からの反射光の周波数に差が生じ、干渉光にビートが生じる。ここで、例えば光源の光周波数が時間的に線形に走査されているとする。干渉計の腕の長さが等しい位置に物体の表面があり、物体の反射面は表面から深さｚの位置にのみあるとする。結合部１７での参照光の周波数と物体からの反射光（物体光）の周波数の時間的変化は、夫々図４の直線Ａ，Ｂのようになる。ここで光周波数は走査レートα［Ｈｚ／ｓ］で、時間Ｔ［ｓ］内で周波数幅Δｆ＝αＴ［Ｈｚ］にわたって走査される。参照光に対する物体光の遅れ時間τは、物体の屈折率をｎとすると、
τ＝２ｎｚ／ｃ
となる。従ってフォトダイオード２３で受光される干渉光は、ビート周波数
ｆｂ＝ατ＝（Δｆ／Ｔ）（２ｎｚ／ｃ）
で変動することになる。
【００４３】
実際は反射光は物体内部の深さに沿って連続的に異なった位置から発生するので、反射光はそれぞれの深さに対応した異なったビート周波数成分をもつ。従って干渉光の強度変化を周波数分析することによって、ビート周波数に対応するある特定の深さからの反射光強度を検出することができる。この反射強度の空間分布をとることで、断層画像を構築できる。
【００４４】
Ｉdct＝（ηｑ／ｈν）｛Ｐｒ＋Ｐｏ∫ｒ（ｚ）ｄｚ＋２（ＰｒＰｏ）1/2
∫ｒ（ｚ）cos（２ｋ（ｔ）ｚ＋φ）ｄｚ｝
干渉光信号Ｉdctの第１，２項はそれぞれ参照ミラーと、物体からの反射光の直流成分であり、第３項が干渉信号光成分であり、式（２）においてＰｒは参照光強度、Ｐｏはプローブ光強度ｒ（ｚ）は深さ方向の反射率分布を示す。Ｉdctの干渉光信号を周波数分析することによって、物体中の任意の深さに対応する散乱光強度の関係を得ることができる。本発明の信号解析部２５は、この式で示される干渉光信号Ｉdctを解析対象信号Ｘ（ｎ）とし、ノンハーモニック分析によって周波数解析するようにしたものである。
【００４５】
［光断層画像表示システムの動作］
次に本実施の形態の動作について説明する。波長走査型光源１０を駆動する。これによって前述したように光ファイバ１４，１８を介して信号光が参照ミラー及び物体にまで照射され、その反射光が結合部１７を介して得られる。そのビート周波数を持つ信号がフォトダイオード２３に得られる。これを増幅した受光信号が信号解析部２５に入力される。又前述したようにｋトリガ発生部２７から等周波数間隔のｋトリガ信号が得られ、このｋトリガ信号に応じてＡ／Ｄ変換を行い、信号解析部２５で受光信号を周波数分析することにより深さ方向に１次元の断層画像が得られる。そしてスキャニングミラー２０を回動させることによって光の入射位置を変化させ、これによって2次元の断面画像を得ることができる。又、この干渉計自体又は測定対象をスキャニングミラー２０による光の走査方向と垂直に移動させることにより、３次元断面画像を得ることができる。
【００４６】
次に信号解析部２５の処理について更に詳細に説明する。ＣＰＵ５１の制御のもとに、信号解析プログラムを実行して入力された信号の周波数解析を行う。なお、周波数解析の対象となる解析対象信号Ｘ（ｎ）は、信号入力インターフェイス５０を介して入力され、ＣＰＵ５１によって読み出し可能に記憶部５４に記憶される。信号解析部２５は、ＣＰＵ５１の制御のもとに、このようにして入力された解析対象信号の周波数解析を行い、その解析結果等を表示部５６に表示出力する。
【００４７】
［周波数解析アルゴリズム］
以下、このような信号解析部２５における周波数解析アルゴリズムについて詳述する。本発明の信号解析部２５にて新たに提案する周波数解析手法（以下、本周波数解析手法という。）は、「周波数の変化を細かく捉えることと、周波数分解能を上げることとは、互いに相反する関係にある」という従来の周波数解析手法の問題を本質的に解決するために、非線形方程式の解を求める問題に着目した。すなわち、本周波数解析手法においては、次式（３）に示す非周期信号のフーリエ変換式の周波数パラメータを求める問題を非線形方程式の最適解を求める問題に置き換えたものである。
【００４８】
【数３】
【００４９】
具体的には、本周波数解析手法においては、次式（４）に示すように、解析対象信号ｘ（ｎ）と正弦波モデル信号との差の二乗和で表される非線形方程式の最適解として、この非線形方程式の右辺が最小値になるような周波数ｆ’、振幅Ａ’、及び初期位相φ’を求める。なお、次式（４）において、Ｌはフレーム長（解析窓長）であり、ｆsはサンプリング周波数［Ｈｚ］である。このように、非線形方程式に周波数ｆ’をパラメータとして盛り込んだことは、いまだ事例がなく斬新である。換言すれば、本周波数解析手法においては、振幅Ａ’及び初期位相φ’のみならず、周波数ｆ’をも正確に求めることを目的とし、上式（３）に対して非線形方程式の解法を用いるものである。本周波数解析手法においては、このような最小二乗法によって非線形方程式の最適解を求める問題に帰着させることにより、解析窓の影響やエイリアシングの影響がなくなり、解析窓長が、１周期未満であってもよく、周期の整数倍でなくてもよく、さらには、不等間隔であってもよい等、柔軟な周波数解析処理を実現することが可能となる。
【００５０】
【数４】
【００５１】
なお、上式（４）は、解析対象信号ｘ（ｎ）が、受光信号、時系列信号等の１次元信号である場合を想定しているが、本周波数解析手法自体は、２次元の画像データ等の２次元信号や動画像データ等の３次元信号、ホログラムのような立体動画像データ等の４次元信号、さらにはそれ以上の次元の信号に対しても適用可能である。すなわち、本周波数解析手法は、１次元の解析対象信号をｓ（ｘn）、２次元の解析対象信号をｓ（ｘm，ｙn）、３次元の解析対象信号をｓ（ｘl，ｙm，ｚn）、４次元の解析対象信号をｓ（ｘk，ｙl，ｚm，ｗn）として一般化すると、解析対象信号が１次元、２次元、３次元、又は４次元の場合のそれぞれについて、次式（５）乃至次式（８）に示す非線形方程式の最適解を求めることになる。
【００５２】
【数５】
【数６】
【数７】
【数８】
【００５３】
換言すれば、本周波数解析手法は、任意のｎ次元（ｎは１以上の整数）の解析対象信号をｓ（ｘ１n1，ｘ２n2，ｘ３n3，・・・，ｘｎnn）として一般化すると、次式（９）に示す非線形方程式の最適解を求めることになる。
【数９】
【００５４】
なお、以下では、説明の便宜上、解析対象信号が１次元であり、上式（４）に示す非線形方程式の最適解を求めるものとする。
【００５５】
さて、上式（４）に示す非線形方程式の最適解を実際に求めるにあたっては、以下のような方法をとることができる。
【００５６】
本周波数解析手法においては、振幅Ａ’、周波数ｆ’、及び初期位相φ’のそれぞれについて適切な初期値を求め、これら初期値から非線形方程式の解法を用いて最適解に収束させる。この非線形問題では、上式（４）をコスト関数とする最小化問題とする。なお、適切な初期値は、離散フーリエ変換（Discrete Fourier Transform；以下、ＤＦＴという。）やウェーブレット変換等の任意の周波数変換を行ったり、フィルタリングを行うことによっておおよその見当をつけたりする等、既存の任意の方法を適用して求めることができる。
【００５７】
まず、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’について、いわゆる最急降下法を適用し、周波数パラメータｆm’，φm’を次式（１０）及び次式（１１）によって求める。
【００５８】
【数１０】
【数１１】
【００５９】
なお、上式（１０）及び上式（１１）においては、次式（１２）と略している。また、μmは、いわゆる減速法に基づく重み係数であり、各漸化式によって求められるコスト関数を単調減少数列にするために、適時０〜１の値をとる。
【００６０】
【数１２】
【００６１】
周波数パラメータｆm’，φm’を求めることができれば、上式（４）における正弦波モデル信号の係数としての周波数パラメータＡ’を一意に求めることができるため、本周波数解析手法においては、次式（１３）によって周波数パラメータＡm’を収束させる。
【００６２】
【数１３】
【００６３】
本周波数解析手法においては、これら一連の計算を反復して行うことにより、振幅Ａ’、周波数ｆ’、及び初期位相φ’を高精度に収束させることができる。特に、本周波数解析手法においては、上式（４）における正弦波モデル信号の位相を構成する周波数パラメータｆ’，φ’と、係数としての周波数パラメータＡ’とを別個に求めることにより、計算を簡便に行うことができる。
【００６４】
しかしながら、最急降下法は、比較的広い範囲から収束するものの、１回の反復では精度が低く、収束するまでに時間を要する。
【００６５】
そこで、本周波数解析手法においては、最急降下法を適用して周波数パラメータｆm’，φm’をある程度まで収束させた後、さらに、いわゆるニュートン法を適用して高精度に収束させるのが望ましい。具体的には、本周波数解析手法においては、ニュートン法として、次式（１４）及び次式（１５）に示す漸化式によって周波数パラメータｆm’，φm’を求める。
【００６６】
【数１４】
【数１５】
【００６７】
ただし、上式（１４）及び上式（１５）において、Ｊは次式（１６）とし、次式（１７）と略している。また、νmもμmと同様に減速法に基づく重み係数であり、適時０〜１の値をとる。
【００６８】
【数１６】
【数１７】
【００６９】
本周波数解析手法においては、上式（１４）及び上式（１５）によって周波数パラメータｆm’，φm’を求めた後、最急降下法と同様に、上式（１３）によって周波数パラメータＡm’を収束させ、この一連の計算をさらに反復して行う。
【００７０】
このように、本周波数解析手法においては、最急降下法とニュートン法とを組み合わせたハイブリッド型の解法を用いることにより、高速に且つ高精度に周波数パラメータＡ’，ｆ’，φ’を推定することができる。
【００７１】
また、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次減算処理することにより、近似的にスペクトルパラメータを導出することができる。ここで、解析対象信号ｘ（ｎ）が複数の正弦波の和であり、次式（１８）のように表されているとする。
【００７２】
【数１８】
【００７３】
パーセヴァル（Parseval）の定理より、解析対象信号ｘ（ｎ）の周波数ｆkと正弦波モデル信号の周波数パラメータｆ’とが全く一致しない場合、すなわち、次式（１９）である場合には、上式（４）に示す非線形方程式は次式（２０）となる。また、周波数パラメータｆ’，φ’の組が、周波数ｆk及び初期位相φkの組のいずれかに一致する場合には、上式（４）に示す非線形方程式は次式（２１）となる。さらに、振幅Ａjが周波数パラメータＡ’とも一致した場合には、解析対象信号から推定スペクトルに関する周波数成分を完全に消去することができる。そのため、最適解を求める問題は、周波数に対して独立であり、解析対象信号から順次個別に推定すれば、複数の正弦波で表される信号にも応用することができる。
【００７４】
【数１９】
【数２０】
【数２１】
【００７５】
すなわち、本周波数解析手法においては、解析対象信号ｘ（ｎ）が複合正弦波の場合であっても、逐次残差信号に対して同様に処理を行い、複数の正弦波を抽出することができる。
【００７６】
ただし、複数のスペクトルの場合には、上式（３）からもわかるように、無限長を仮定しているため、複数のスペクトルの周波数同士が近い場合には、合成スペクトルの周期が解析窓長よりも数段長くなり、上式（１７）が満たされなくなるため、誤差が発生する場合もある。
【００７７】
従来のＯＣＴでは、受光信号を複合正弦波によって表現するために多くのスペクトル数（正弦波の数）が必要であったが、本周波数解析手法においては、そのような信号であっても僅かなスペクトル数で誤差なく表現することができる。すなわち、信号をより少ないスペクトル数で表現可能であることは、情報圧縮の用途に有効であることは勿論のこと、それ以外にも、物理現象に関して波数表現をするようなアプローチに用いることにより、信号の本質的な特性も見抜くことができることを示している。
【００７８】
［本周波数解析手法の有効性］
以下、本周波数解析手法の有効性について具体的に説明する。
【００７９】
従来、センサによって計測された信号を、計算機を用いて解析するためには、一般に、例えば図５に示すように、その信号の一部を切り出し、ＤＦＴといった代表的な調和型の周波数解析手法によって解析していた。
【００８０】
しかしながら、解析対象信号を切り出して（窓掛けして）調和型の解析を行うことは、当該解析対象信号が、切り出した一部の信号が完全周期で繰り返される信号であると擬制して解析することに他ならず、当然ながらその周波数特性も、元のスペクトル特性とは完全に一致しないものとなる。また、この信号切り出しの影響を回避するために、窓関数によって影響を軽減することも試みられているが、切り出した一部の信号を完全周期として解析することには変わりなく、本質的な問題解決には至っていないのが現状である。
【００８１】
これに対して、本周波数解析手法は、非周期信号を想定した周波数解析を目的として定式化を行っていることから、解析窓長と周波数分解能の制約を受けることなく、正確に周波数とそのパラメータとを推定することができる。
【００８２】
また、本周波数解析手法は、周波数の推定についても、従来の離散的な探索から、動的且つ連続的な探索となるため、格段に精度を向上させることができ、革新的な手法である。
【００８３】
さらに、多次元化に関しても、従来の周波数解析手法においては、解析対象信号が１次元である場合と同様に窓関数の影響を受ける。具体的には、例えば図６に示すように、２次元の解析対象信号（原画像データ）のうち図６左図中の破線部を切り出してＤＦＴを施して解析したとしても、その切り出した区間を完全周期とした画像が延々と繰り返されるものを解析したことと等価であることから、右上図に示すように、復元した画像データが窓関数の影響を受けたものとなり、本質的な特性を見出すことは困難である。
【００８４】
これに対して、上式（６）を適用した本周波数解析手法においては、図６右下図に示すように、窓関数の影響が軽減されている。そのため本質的な周波数特性を見出すことができ、切り出した区間の周囲の画像も完全に再現することができる。
【００８５】
以上のように、本周波数解析手法は、非線形方程式の最適解を求めることにより、正弦波モデル信号の周波数ｆ’、振幅Ａ’、及び初期位相φ’を高速に且つ高精度に求めることができる。具体的な精度を立証するために、本願発明者は、ＤＦＴと、ＤＦＴの発展型のうち最も解析精度が高いといわれているＧＨＡ（Generalized Harmonic Analysis）とを比較対象として精度の検証を行った。
【００８６】
なお、ＤＦＴやＧＨＡは、１つの解析窓長に見かけ上複数の窓長を持たせていることから、周波数分解能が解析窓長に依存するが、その分解周波数が有限長であり、解析対象信号の周波数が分解周波数以外の周波数となった場合には解析することができず、解析対象信号が正確に解析できる周波数と異なる場合には、最も近い分解周波数の他に、その周辺に小さなスペクトルの周波数（側帯波成分）が現れ、複数の周波数が出現してしまう。
【００８７】
このような現象が本周波数解析手法においても生じるか否かについて、すなわち、本周波数解析手法の周波数分解能を検証するために、解析窓長を１秒（１０２４サンプル）とした１次元の非常に短い単一正弦波を解析し、各手法によって正弦波を１本抽出して元の信号との二乗誤差を調べた。その結果を図７に示す。
【００８８】
図７に示すように、ＤＦＴにおいては、基本周波数の整数倍以外の周波数における解析精度の悪化がみられた。また、ＧＨＡにおいては、１Ｈｚ以上の周波数ではＤＦＴと比べて２〜５桁程度の精度向上がみられた。これに対して、本周波数解析手法においては、１Ｈｚ以上の周波数ではＤＦＴと比べて１０桁以上、ＧＨＡと比べて５桁以上の精度向上がみられた。すなわち、本周波数解析手法は、既存の周波数解析手法と比べて１０万〜１００億倍以上の精度向上がみられた。
【００８９】
特に、１Ｈｚ以下の周波数を正確に推定することができるということは、解析窓長を超えた長い周期信号であっても解析可能であることを示しており、そのスペクトル構造をより正確に推定可能であることを示している。
【００９０】
このように、本周波数解析手法は、最も解析精度が高いといわれているＧＨＡと比べても驚くべき高精度に解析を行うことができる。国内外の様々な研究事例をみても、これほどの精度で信号を解析できる手法は存在せず、本周波数解析手法は、今後、正確な解析が必要とされる様々な分野への応用が期待できる手法であるといえる。
【００９１】
［画像の符号化］
本周波数解析手法においては、窓の影響が少ないため、スペクトルに側帯波成分が全く出ない。すなわち、本周波数解析手法は、上述したように、複雑な解析対象信号であっても、僅かなスペクトル数で効率的に誤差なく表現することができる。
【００９２】
図８に２次元信号をＤＦＴと本周波数解析手法とのそれぞれによって解析した具体例を示す。図８左図に示すように、本来であれば１本のスペクトルで表現される信号をＤＦＴによって変換した場合には、一般に、変換して得られるスペクトルも、解析窓の影響に起因して、図８右上図に示すように、複数の側帯波成分が群れて現れてしまう。これは、画像の圧縮符号化において広く用いられている離散コサイン変換（Discrete Cosine Transform；ＤＣＴ）等においても同様である。これに対して、本周波数解析手法は、図８右下図に示すように、元の信号のスペクトルを正確に推定することができる。正確に解析できると、画像データを表現するスペクトル数を格段に削減できる可能性があり、画質を落とさずに飛躍的な情報圧縮効果を期待することができる。
【００９３】
又この実施の形態では波長走査型光源の走査範囲をあまり広くすることなく高分解能の画像が得られる。そのため従来では使用できなかった安価な波長走査型光源を用いて必要な分解能の断層画像を得ることができる。又広範囲の波長を走査することができる光源を用いれば、より精密な断面画像とすることができる。
【００９４】
尚この実施の形態ではｋトリガ発生部により波長の１スキャニングの間に等周波数間隔となるタイミングで多数のｋトリガ信号を発生するようにしているが、波長の１走査の開示のみのスキャントリガ信号としてもよい。この場合に特開２００８−２６１７７８号に示されるように、物体の測定位置に反射物を配置したときに受光素子からの出力をデータ番号を付して等時間間隔で取り込み、ビート信号のピークとなるデータ番号を算出し、ピーク番号に対してデータ番号を示す近似曲線を算出し、ピーク番号の変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出し、データ取り込み数列とし、測定物体を配置したときの受光素子からの等時間間隔からの出力からデータ数列の各要素でのタイミングをトリガ信号として信号解析部に与えるようにしてもよい。
【００９５】
（第２の実施の形態）
図９は本発明の第２の実施の形態による光断層画像表示システムの全体構成を示すブロック図である。本実施の形態はＳＤ−ＯＣＴによる光断層画像表示システムであり、光源として所定の波長範囲に一様な光を出射する光源４０を用いる。この光源４０からの光は光増幅器１３を介して光干渉計に導かれることは第１の実施の形態と同様である。この実施の形態ではレンズ２２に得られる光を回折格子４１に入射する。回折格子４１は入射光の波長に応じて光を分散させる波長分光手段であり、分散した光は一次元ＣＣＤ４２に入力される。一次元ＣＣＤ４２は所定の周期で各位置の光信号を一次元の受光信号に変換するものであり、その出力は増幅器４３を介して信号解析部２５に伝えられる。この場合には、ＣＣＤ４２の一次元の各画素からの出力のタイミングに周期してトリガ信号を入力インターフェイスに伝える。これ以降の構成については第１の実施の形態と同様である。このようにして１次元ＣＣＤで得られた受光信号に基づいて周波数を解析するようにしても、同様に断層画像を得ることができる。
【００９６】
又前述した実施の形態における信号解析装置ではソフトウェアによる周波数解析を行うものとして説明したが、本発明は、本周波数解析手法のアルゴリズムを実装したＤＳＰ（Digital Signal Processor）等、積和演算を行うことが可能であればハードウェアによっても実現することができる。
【００９７】
このように、本発明は、その趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【特許請求の範囲】
【請求項１】
周期的に光の発振波長を走査する波長走査型光源と、
前記波長走査型光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、
前記干渉光学計より得られる干渉光を受光し、ビート信号を得る受光素子と、
前記受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、
前記信号解析部は、
前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、
前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、
前記演算手段より得られる断面画像を表示する表示部と、を具備する光断層画像表示システム。
【請求項２】
所定の範囲で一様な波長の光を発振する光源と、
前記光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、
前記干渉光学計から得られる干渉光を波長に基づいて分岐する波長分光手段と、
前記波長分光手段より分光された範囲の干渉光を受光して電気信号に変換する一次元受光素子と、
前記一次元受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、
前記信号解析部は、
前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、
前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、
前記演算手段より得られる断面画像を表示する表示部と、を具備する光断層画像表示システム。
【請求項３】
前記演算手段は、前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’のそれぞれについて適切な初期値を求め、求めた初期値から非線形方程式の最適解に収束させて前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を求めることを特徴とする請求項１又は２記載の光断層画像表示システム。
【請求項４】
前記演算手段は、前記正弦波モデル信号の位相を構成する前記周波数ｆ’及び前記初期位相φ’について最急降下法を適用して収束させ、当該周波数ｆ’及び当該初期位相φ’を求め、前記正弦波モデル信号の係数としての前記振幅Ａ’について最急降下法を適用して収束させ、当該振幅Ａ’を求めることを特徴とする請求項３記載の光断層画像表示システム。
【請求項５】
前記演算手段は、前記最急降下法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させた後、さらに、ニュートン法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させることを特徴とする請求項４記載の光断層画像表示システム。
【請求項６】
前記波長走査型光源の１走査の期間内に、前記波長走査型光源の光の等周波数間隔でのトリガ信号を発生するｋトリガ発生部を更に有し、
前記信号解析部の記憶部は、前記ｋトリガ発生部により発生したｋトリガ信号のタイミングで前記干渉信号を保持するものである請求項１記載光断層画像表示システム。
【請求項７】
前記信号解析部の演算手段は、
物体の測定位置に反射物を配置したときに、前記受光素子からの出力をデータ番号を付して等時間間隔で取り込み、前記ビート信号のピークとなるデータ番号を算出し、ピーク番号に対するデータ番号を示す近似曲線を算出し、前記ピーク番号変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出しデータ取り込み数列とし、測定位置に測定物体を配置したときの前記受光素子からの等時間間隔の出力から、前記データ取り込み数列の各要素でのタイミングの信号レベルを算出し、得られた等周波数間隔のデータを前記記憶手段に保持するものである請求項１記載の光断層画像表示システム。
【請求項１】
周期的に光の発振波長を走査する波長走査型光源と、
前記波長走査型光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、
前記干渉光学計より得られる干渉光を受光し、ビート信号を得る受光素子と、
前記受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、
前記信号解析部は、
前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、
前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、
前記演算手段より得られる断面画像を表示する表示部と、を具備する光断層画像表示システム。
【請求項２】
所定の範囲で一様な波長の光を発振する光源と、
前記光源からの光を参照光と物体への照射光とに分岐し、物体からの反射光と参照光との干渉光を発生する干渉光学計と、
前記干渉光学計から得られる干渉光を波長に基づいて分岐する波長分光手段と、
前記波長分光手段より分光された範囲の干渉光を受光して電気信号に変換する一次元受光素子と、
前記一次元受光素子から得られる干渉信号に対して周波数解析を行い、物体の断層画像を形成する信号解析部と、を具備し、
前記信号解析部は、
前記受光素子から得られる干渉信号を解析対象信号として保持する記憶部と、
前記記憶部に記憶された前記解析対象信号を読み出し、前記解析対象信号と、周波数ｆ’及び初期位相φ’を用いた位相と振幅Ａ’とによって表される正弦波モデル信号との差の二乗和が最小値になるような前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を算出し、求めた前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’に基づき断面画像を出力する演算手段と、
前記演算手段より得られる断面画像を表示する表示部と、を具備する光断層画像表示システム。
【請求項３】
前記演算手段は、前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’のそれぞれについて適切な初期値を求め、求めた初期値から非線形方程式の最適解に収束させて前記周波数ｆ’、前記振幅Ａ’、及び前記初期位相φ’を求めることを特徴とする請求項１又は２記載の光断層画像表示システム。
【請求項４】
前記演算手段は、前記正弦波モデル信号の位相を構成する前記周波数ｆ’及び前記初期位相φ’について最急降下法を適用して収束させ、当該周波数ｆ’及び当該初期位相φ’を求め、前記正弦波モデル信号の係数としての前記振幅Ａ’について最急降下法を適用して収束させ、当該振幅Ａ’を求めることを特徴とする請求項３記載の光断層画像表示システム。
【請求項５】
前記演算手段は、前記最急降下法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させた後、さらに、ニュートン法を適用して前記周波数ｆ’及び前記初期位相φ’を収束させることを特徴とする請求項４記載の光断層画像表示システム。
【請求項６】
前記波長走査型光源の１走査の期間内に、前記波長走査型光源の光の等周波数間隔でのトリガ信号を発生するｋトリガ発生部を更に有し、
前記信号解析部の記憶部は、前記ｋトリガ発生部により発生したｋトリガ信号のタイミングで前記干渉信号を保持するものである請求項１記載光断層画像表示システム。
【請求項７】
前記信号解析部の演算手段は、
物体の測定位置に反射物を配置したときに、前記受光素子からの出力をデータ番号を付して等時間間隔で取り込み、前記ビート信号のピークとなるデータ番号を算出し、ピーク番号に対するデータ番号を示す近似曲線を算出し、前記ピーク番号変化分を任意の数に分割して、分割したピーク番号に対応するデータ番号を近似曲線により算出しデータ取り込み数列とし、測定位置に測定物体を配置したときの前記受光素子からの等時間間隔の出力から、前記データ取り込み数列の各要素でのタイミングの信号レベルを算出し、得られた等周波数間隔のデータを前記記憶手段に保持するものである請求項１記載の光断層画像表示システム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１０−２２３６７０（Ｐ２０１０−２２３６７０Ａ）
【公開日】平成２２年１０月７日（２０１０．１０．７）
【国際特許分類】
【出願番号】特願２００９−６９６３２（Ｐ２００９−６９６３２）
【出願日】平成２１年３月２３日（２００９．３．２３）
【出願人】（３０５０６０５６７）国立大学法人富山大学 (194)
【出願人】（５９１１０２６９３）サンテック株式会社 (57)
【Ｆターム（参考）】
【公開日】平成２２年１０月７日（２０１０．１０．７）
【国際特許分類】
【出願日】平成２１年３月２３日（２００９．３．２３）
【出願人】（３０５０６０５６７）国立大学法人富山大学 (194)
【出願人】（５９１１０２６９３）サンテック株式会社 (57)
【Ｆターム（参考）】
[ Back to top ]