みなとみらい夕日T
2010.12．23みなとみらい夕日、遠くには富士山が見えます。

OCT理論学習, Optical Coherence Tomography

　参考文献をもとにOCTにおける光干渉波形をシミュレーションしながら、原理の理解を進めていきたいと思います。
　2002年に"Handbook of Optical Coherence Tomography "の初版が出ていますので、その頃がOCT開発のピークだと推測されます。2000年ごろと言えば、光通信の開発もバブル状態でした。半導体レーザーや波長可変レーザ、光ファイバ増幅器など光部品の開発が進み、OCTにも同様な光部品が使われることから、OCT開発にも弾みがついたものと思います。

　そして、2008年の春から、眼科において3次元画像解析装置を使った検査は保険がきくようになったといいます。市民権を得たということでしょう。3次元画像を見れば、なくてはならない装置であることは自明ですが、ここまで育てた苦労話はいろんなwebサイトに紹介されています。（これを書いているのは2010年12月です）

　開発当初は低コヒーレント光源（SLDやASE光源などの広帯域光源）を使っていて、参照ミラーを動かしていました。最近の3次元画像解析装置のOCTは、干渉し易いレーザ光が使われています。当初からレーザ光を使えばいいじゃないかと思うかも知れませんが、OCTは低コヒーレント光源から出発しています。コヒーレント長が短い低コヒーレント光源を使ったところに特長があったのです。

　その後眼科での応用などがきっかけでOCTの改造が進んでいきます。そしてOCTの光源は、干渉に欠かせない波長可変レーザ光を使うようになりました、さらに解像度を上げるため波長可変レーザの帯域を広げる（可変波長範囲が広いほうが解像度（分解能）が上がる）などの開発が進みました。測定原理はよく知られたマイケルソン干渉光学系なのですが、光源に何を使ったかで進化してきました。また光ファイバで干渉系を構成しているのも特徴かもしれません。

　 代表的なOCTの種類
・Time domain OCT-- 参照ミラーを動かすタイプ
・Fourier domain OCT (FD-OCT)--参照ミラーは固定、広帯域光源、検出光を分光器で周波数領域に変換
・Swept Source OCT(SS-OCT)--参照ミラー固定波長可変レーザ
・Full-field OCTもあります
などが提案されています。光源、光源の走査やデータ処理、３D表示などが重要要素ですが、ここでは触れていません。応用分野は、光と相性の良い眼科です。

波長可変レーザを使ったマイケルソン干渉計

　まずマイケルソン干渉計の復習から始めてみます。図はこちらを参照してください。次のように書いてあります。
「2つの経路の長さが波長の整数倍（0を含む）の場合、2つの光線は互いに強め合うように干渉し、検出器は強い信号を検出する。経路長の差が波長の整数倍と2分の1の場合（例えば、0.5波長、1.5波長、2.5波長など）、2つの光線は互いに弱め合うように干渉し、検出器は弱い信号を検出する。」

　ここで、光路差を固定した場合を考てみます。光源が単一波長の場合、干渉光の光強度は明暗のどこか一定値になるでしょう。ここで波長を変化させれば位相差が変わるので明暗の強度も変わります。
　すなわち、光路差一定でも光源を波長変レーザにすれば、干渉光は明暗を繰り返します。そして、この明暗の繰り返しは位相差を意味します。光路差が大きくなれば、レーザの波数変化に対して位相差は早く変化します、すなわち対象物の位置が離れているほど明暗の繰り返し周波数は高くなります。対象物(サンプル)の位置というのは、OCTで言えば、生体内部で測定光を反射する組織の位置ということになるのでしょう。このようにして出来上がったのがSS-OCTです。
　歴史から言うと、SS-OCTの前にFD-OCTが開発されています（たぶん）のでFD-OCTから調べてみます。

FD-OCT(Fourier domain OCT)

SS-OCT(Swept Source OCT)

　FD-OCTとSS-OCTは広帯域光源を使うのか、波長可変光源を使うのか異なりますが、光干渉の考え方は同じですのでSS-OCTで説明します。
　干渉系の図、式は、参考文献（２）、あるいは参考文献（６）があります。参考文献（６）は光波の実数部だけを扱っていますのでわかりやすいのでこちらでシミュレーションしてみます。

【光波を実数部だけで表わす計算例】
　参考文献（６）で
　検出器に入る２つの光の光路差は
　s2-s1 = 2(Δz) -----------------------------(1)
　　*注）参考文献（６）の光路差L2-L1をここではΔzとした。
　検出器に入る光強度は
　I=(A1cos(wt-ks1)+A2cos(wt-ks2))^2
　=(1/2)*A1^2+(1/2)*A2^2+A1*A2*cos2k(Δz)-----(2)
　　計算が簡単に見えますね。
　(2)をFFTすることで、サンプルの内部情報を求めることができます。
　FFT(I)--------------------------------------(3)

λ=c/freq
波数k=2*pi/λ=2pi*freq/cなのでサンプルの屈折率n、光速cを考慮すると式（２）において
位相差2k(Δz)=2*2*pi*freq/c*(n*Δz)となる。------(4)

周波数freqを変化(sweep)させた時に式（２）の出力はどうなるのだろうか
対象物（サンプル）位置の計算
式（４）の位相差が波長の整数倍の時、検出される干渉光の強度は大きくなるので
強度が大きくなる隣り合う波長をλ1、λ2とすれば
　2*n*Δz=N1*λ1 (N1は整数,2は往復)
　2*n*Δz=N2*λ2 (N2は整数,2は往復)
　N2-N1=2*n*Δz(1/λ2-1/λ1)=1　（隣り合う波長なのでN2-N1=1）
　周波数で表わすと
　Δz=(1/(2*n))/((freq2-freq1)/c)----------------(5)
　freq1=干渉光強度がピークになる周波数
　freq2=freq1の隣で干渉光強度がピークになる周波数
　
式（５）からサンプルの位置は求められますが、実際には式（３）のように式（２）をFFTして求めることになります。

シミュレーション
　式（２）を計算してみます。
　波数　k=2*pi/λ
　シミュレーションは波長1.31um付近とします。波長と周波数の関係は下記のようになります。
　　　周波数と　　　　,波長の関係　　,波数は
　2.2880*10^14(Hz)　,　1.31(um)　,　4.7963(1/um)
　2.1257*10^14(Hz)　,　1.41(um)　,　4.4562(1/um)
　1.9849*10^14(Hz)　,　1.51(um)　,　4.1610(1/um)
　1.9982*10^14(Hz)　,　1.50(um)　,　4.1888(1/um)

ここで、波長可変レーザーの発振周波数を
freq=1.9982*10^14（波長1.5um）～2.2542*10^14Hz（波長1.3296um）まで0.001*10^14Hz刻みで、256段階とします。

　Fig.1とFig.2はサンプル位置がそれぞれ50um、200umの干渉光強度とそれをFFT変換したものです。サンプル位置Δzが大きくなると干渉光ピークの数が増えます。オシロスコープで見れば繰り返し周波数が大きくなったように見えるのでしょう。

SS-OCT

Fig.1-1はサンプルの位置Δz=50umの例です。 Fig.1-1において
　freq1=1.9992*10^14
　freq2=2.0222*10^14
　式（５）よりΔz=48.8um
　サンプル位置50umに対して、1.2umの測定誤差で、サンプル位置を測定することが出来ます。
Fig.1-2はFFTしたものです。最初のピークの位置は０から数えて11pointsのところです。
式（５）Δz=(1/(2*n))/((freq1-freq2)/c)=(1/(2*n))/((終了周波数-開始周波数)/FFTピーク位置）/c)
　　　　　=(1/n)/(((2.2542*10^14-1.9982*10^14)/11)/c)=48.27um
　　　　　誤差1.8umです。
FFTするとFFTpoints/2を中心にして折り返すような形になりますが、これは実数をFFTする場合の性質ですので、左半分で考えます。

SS-OCT

Fig.2-1はサンプルの位置Δz=200umの例です。


以下計算プログラムです。ソフトウェアはOctaveを使っています。
clear;
n=1.334;%屈折率水
c=2.9972458*10^14;%光速 u/s;
z1=50;%サンプルの位置um
z2=200;%サンプルの位置um

m=0;
wave=1.31;
%freq=c/wave;
freq=1.9982*10^14;

freq=1.9982*10^14:0.001*10^14:2.2542*10^14;#周波数Hz 256step
octi1=(cos(2*2*pi*freq/c*n*z1)); # z=50um
octi2=(cos(2*2*pi*freq/c*n*z2)); # z=200um   

octf1=fft(octi1);
octf2=fft(octi2);

figure(1)
subplot(2,1,1)
x=freq;
plot(x,real(octi1),x,imag(octi1))
grid 'on'
title('Fig.1-1 SS-OCT freq=199.82THz to 225.36THz, z=50um Blue=real part') 
xlabel('Frequency(Hz)') 


subplot(2,1,2)
x=1:1:256;
plot(x,real(octf1(x)))
grid 'on'
title('Fig.1-2 SS-OCT Fourier transformation Blue=real part') 
xlabel('FFT points') 


figure(2)
subplot(2,1,1)
x=freq;
plot(x,real(octi2),x,imag(octi2))
grid 'on'
title('Fig.2-1 SS-OCT freq=199.82THz to 225.36THz, z=200um Blue=real part,Green=imag') 
xlabel('Frequency(Hz)') 


subplot(2,1,2)
x=1:1:256;

plot(x,real(octf2(x)))
grid 'on'
title('Fig.2-2 SS-OCT Fourier transformation Blue=real part') 
xlabel('FFT points')

【光電界を複素表示で解析する計算例】
　参考文献（２）は光電界を複素表示で扱っています。
　式(2.2)=参照ミラーからの反射光
　式(2.4)=サンプル（被測定物）からの反射光
　式(2.5)=式(2.2)+式(2.4)=参照ミラーの反射光とサンプルからの反射光の合成波（干渉光）

　マイケルソ干渉系の出力はPDなどで受光する場合、光電界の絶対値の2乗が受光されるから
　式(2.7)=｜式(2.5｜^2 = 式(2.5)*(式(2.5)の共役複素）
　簡単にしてまとめると次の形になる。
　I(ω)=A1*|H(ω)|^2+A1+1*A1*H(ω)--------------(2.10)
　
　ここではサンプル側にはサンプルにおける光の入出力応答関数H(ω)が使われているが、どのような形か定義されていない。一般的に考えれば、位相遅れ関数と考えればよさそうな気がする。そこで
　H(ω)=exp[-jφ(Δz)]で表わし、φ(Δz)=2ωnΔz/c とする------式(A）

式(2.10)I(ω)をシミュレーションしてみます。条件は【光波を実数部だけで表わす計算例】と同じです。結果も同じになります。複素計算をしていますが、実際にPDに受光されるのは実部だけですので実数部を表示しています。

SS-OCT

Fig.3-1はサンプルの位置Δz=50umの例

SS-OCT

Fig.4-1はサンプルの位置Δz=200umの例

SS-OCT

Fig.3-1-AはFig.3-1と同じ条件で複素表示したもの。Fig.3-2-AでFFTしたときに左半分のスペクトラムが消えている？なぜか？


複素計算の場合、表示グラフは実数部のみ
clear;
n=1.334;%屈折率水
c=2.9972458*10^14;%光速 u/s;
z1=50;%サンプルの位置um
z2=200;%サンプルの位置um

m=0;
wave=1.31;
%freq=c/wave;
freq=1.9982*10^14;

freq=1.9982*10^14:0.001*10^14:2.2542*10^14;#周波数Hz 256step
octi1=exp(-j*(4*pi*freq*n*z1/c)); # z
octi1=octi1*octi1.'+octi1;

octi2=exp(-j*(4*pi*freq*n*z2/c)); # z   
octi2=octi2*octi2.'+octi2;

octf1=fft(real(octi1));
octf2=fft(real(octi2));

figure(3)
subplot(2,1,1)
x=freq;
plot(x,real(octi1))

grid 'on'
title('Fig.3-1 SS-OCT freq=199.82THz to 225.36THz, z=50um Blue=real part') 
xlabel('Frequency(Hz)') 

subplot(2,1,2)
x=1:1:256;
plot(x,real(octf1(x)))
grid 'on'
title('Fig.3-2 SS-OCT Fourier transformation Blue=real part') 
xlabel('FFT points') 

figure(4)
subplot(2,1,1)
x=freq;
plot(x,real(octi2))
grid 'on'
title('Fig.4-1 SS-OCT freq=199.82THz to 225.36THz, z=200um Blue=real part') 
xlabel('Frequency(Hz)')  

subplot(2,1,2)
x=1:1:256;

plot(x,octf2(x))
axis([0 256 -100 50])
grid 'on'
title('Fig.4-2 SS-OCT Fourier transformation Blue=real part') 
xlabel('FFT points')

広帯域光源と波長可変ファイバレーザは何時も競合する

　広帯域光源を使った測定には、測定機として光スペクトラムアナライザ（分光器）が必要になります。また波長可変レーザを使えば、光パワーメータが必要になります。これはいろいろなパターン当てはまります。これらは何時も競合しているのです。
　OCTの例ですと
（１）FD-OCT、TD-OCT
（２）SS-OCT
OCT以外では
（３）FBGセンシング装置
（４）光部品の波長損失特性の測定
（５）光ファイバの偏波分散特性測定
などがあり、比較して光源を選択する必要があります。
参考文献
（１）Wikipedia, マイケルソン干渉計
（２）草刈　修「光駆動による分散チューニングを用いた高速・広帯域波長可変ファイバレーザ」東京大学修士論文、平成22年2月9日
（３）OCT news
（４）Wikipedia,Optical coherence tomography
（５）edited by Brett E.Bouma,Guillermo J.Tearney "Handbook of Optical Coherence Tomography "Marcel Dekker,Inc 2002.
（６）熊谷、「量子エレクトロニクス」2008