以下、添付図面を参照しながら本発明によるデータ作成装置、光制御装置、データ作成方法、及びデータ作成プログラムの実施の形態を詳細に説明する。なお、図面の説明において同一の要素には同一の符号を付し、重複する説明を省略する。
(第1実施形態)
図1は、本発明の第1実施形態に係る光制御装置1Aの構成を概略的に示す図である。図2は、光制御装置1Aが備える光学系10の構成を示す図である。本実施形態の光制御装置1Aは、入力光Laから、該入力光Laとは異なる任意の時間強度波形を有する出力光Ldを生成する。図1に示されるように、光制御装置1Aは、光源2、光学系10、及び変調パターン算出装置(データ作成装置)20を備える。
光源2は、光学系10に入力される入力光Laを出力する。光源2は例えば固体レーザ光源やファイバーレーザ光源等のレーザ光源であり、入力光Laは例えばコヒーレントなパルス光である。光学系10は、SLM14を有しており、変調パターン算出装置20からSLM14の各画素を制御するための制御信号SCをSLM14に受ける。光学系10は、光源2からの入力光Laを、任意の時間強度波形を有する出力光Ldに変換する。変調パターンは、SLM14を制御するためのデータであり、複素振幅分布の強度あるいは位相分布の強度をファイルに出力されたデータである。変調パターンは、例えば、計算機合成ホログラム(Computer-Generated Holograms(CGH))である。
図2に示されるように、光学系10は、回折格子12、レンズ13、SLM14、レンズ15、及び回折格子16を有する。回折格子12は本実施形態における分光素子であり、光源2と光学的に結合されている。SLM14はレンズ13を介して回折格子12と光学的に結合されている。回折格子12は、入力光Laを波長成分毎に分光する。なお、分光素子として、回折格子12に代えてプリズム等の他の光学部品を用いてもよい。また、分光素子は反射型であってもよく、透過型であってもよい。入力光Laは、回折格子12に対して斜めに入射し、複数の波長成分に分光される。この複数の波長成分を含む光Lbは、レンズ13によって各波長成分毎に集光され、SLM14の変調面に結像される。レンズ13は、光透過部材からなる凸レンズであってもよく、凹状の光反射面を有する凹面鏡であってもよい。
SLM14は、入力光Laとは異なる任意の時間強度波形を有する出力光Ldを生成するために、光Lbの位相変調と強度変調とを同時に行う。SLM14は、強度変調のみを行ってもよい。SLM14は、例えば位相変調型である。一実施例では、SLM14はLCOS(Liquid crystal on silicon)型である。或いは、SLM14はデジタルマイクロミラーデバイス(DMD)などの強度変調型SLMであってもよい。また、SLM14は反射型であってもよく、透過型であってもよい。図3は、SLM14の変調面17を示す図である。図3に示されるように、変調面17には、複数の変調領域17aが或る方向Aに沿って並んでおり、各変調領域17aは方向Aと交差する方向Bに延びている。方向Aは、回折格子12による分光方向である。この変調面17はフーリエ変換面として働き、複数の変調領域17aのそれぞれには、分光後の対応する各波長成分が入射する。SLM14は、各変調領域17aにおいて、入射した各波長成分の位相及び強度を他の波長成分から独立して変調する。なお、本実施形態のSLM14は位相変調型であるため、強度変調は、変調面17に呈示される位相パターン(位相画像)によって実現される。
SLM14によって変調された変調光Lcの各波長成分は、レンズ15によって回折格子16上の一点に集められる。このときのレンズ15は、変調光Lcを集光する集光光学系として機能する。レンズ15は、光透過部材からなる凸レンズであってもよく、凹状の光反射面を有する凹面鏡であってもよい。また、回折格子16は合波光学系として機能し、変調後の各波長成分を合波する。すなわち、これらのレンズ15及び回折格子16により、変調光Lcの複数の波長成分は互いに集光・合波されて出力光Ldとなる。
レンズ15よりも前の領域(スペクトル領域)と、回折格子16よりも後ろの領域(時間領域)とは、互いにフーリエ変換の関係にあり、スペクトル領域における位相変調及び強度変調は、時間領域における時間強度波形に影響する。従って、出力光Ldは、SLM14の変調パターンに応じた、入力光Laとは異なる所望の時間強度波形を有することとなる。ここで、図4の(a)は、一例として、単パルス状の入力光Laのスペクトル波形(スペクトル位相G11及びスペクトル強度G12)を示し、図4の(b)は、該入力光Laの時間強度波形を示す。また、図5の(a)は、一例として、SLM14において矩形波状の位相スペクトル変調を与えたときの出力光Ldのスペクトル波形(スペクトル位相G21及びスペクトル強度G22)を示し、図5の(b)は、該出力光Ldの時間強度波形を示す。図4の(a)及び図5の(a)において、横軸は波長(nm)を示し、左の縦軸は強度スペクトルの強度値(任意単位)を示し、右の縦軸は位相スペクトルの位相値(rad)を示す。また、図4の(b)及び図5の(b)において、横軸は時間(フェムト秒)を表し、縦軸は光強度(任意単位)を表す。この例では、矩形波状の位相スペクトル波形を出力光Ldに与えることにより、入力光Laのシングルパルスが、出力光Ldとして高次光を伴うダブルパルスに変換されている。なお、図4及び図5に示されるスペクトル及び波形は一つの例であって、様々な位相スペクトル及び強度スペクトルの組み合わせにより、出力光Ldの時間強度波形を様々な形状に整形することができる。
再び図1を参照する。変調パターン算出装置20は、例えば、パーソナルコンピュータ、スマートフォン及びタブレット端末といったスマートデバイス、あるいはクラウドサーバといった、プロセッサを有するコンピュータである。変調パターン算出装置20は、SLM14と電気的に接続されており、出力光Ldの時間強度波形を所望の波形に近づけるための位相変調パターンを算出し、該位相変調パターンを含む制御信号SCをSLM14に提供する。本実施形態の変調パターン算出装置20は、所望の波形を得る為の位相スペクトルを出力光Ldに与える位相変調用の位相パターンと、所望の波形を得る為の強度スペクトルを出力光Ldに与える強度変調用の位相パターンとを含む位相パターンをSLM14に呈示させる。そのために、変調パターン算出装置20は、任意波形入力部21と、位相スペクトル設計部22と、強度スペクトル設計部23と、変調パターン生成部(データ生成部)24とを有する。すなわち、変調パターン算出装置20に設けられたコンピュータのプロセッサは、任意波形入力部21の機能と、位相スペクトル設計部22の機能と、強度スペクトル設計部23の機能と、変調パターン生成部24の機能とを実現する。それぞれの機能は、同じプロセッサにより実現されてもよいし、異なるプロセッサにより実現されてもよい。
図6は、変調パターン算出装置20のハードウェアの構成例を概略的に示す図である。図6に示されるように、変調パターン算出装置20は、物理的には、プロセッサ(CPU)201、ROM202及びRAM203等の主記憶装置、キーボード、マウス及びタッチスクリーン等の入力デバイス204、ディスプレイ(タッチスクリーン含む)等の出力デバイス205、他の装置との間でデータの送受信を行うためのネットワークカード等の通信モジュール206、ハードディスク等の補助記憶装置207などを含む、通常のコンピュータとして構成され得る。
コンピュータのプロセッサ201は、変調パターン算出プログラム(データ作成プログラム)によって、上記の各機能(任意波形入力部21、位相スペクトル設計部22、強度スペクトル設計部23、及び変調パターン生成部24)を実現することができる。故に、変調パターン算出プログラムは、コンピュータのプロセッサ201を、変調パターン算出装置20における任意波形入力部21、位相スペクトル設計部22、強度スペクトル設計部23、及び変調パターン生成部24として動作させる。変調パターン算出プログラムは、例えば補助記憶装置207といった、コンピュータの内部または外部の記憶装置(記憶媒体)に記憶される。記憶装置は、非一時的記録媒体であってもよい。記録媒体としては、フレキシブルディスク、CD、DVD等の記録媒体、ROM等の記録媒体、半導体メモリ、クラウドサーバ等が例示される。
任意波形入力部21は、操作者からの所望の時間強度波形の入力を受け付ける。操作者は、所望の時間強度波形に関する情報(例えばパルス幅、パルス数など)を任意波形入力部21に入力する。所望の時間強度波形に関する情報は、位相スペクトル設計部22及び強度スペクトル設計部23に与えられる。位相スペクトル設計部22は、与えられた所望の時間強度波形の実現に適した、出力光Ldの位相スペクトルを算出する。強度スペクトル設計部23は、与えられた所望の時間強度波形の実現に適した、出力光Ldの強度スペクトルを算出する。変調パターン生成部24は、位相スペクトル設計部22において求められた位相スペクトルと、強度スペクトル設計部23において求められた強度スペクトルとを出力光Ldに与えるための位相変調パターン(例えば、計算機合成ホログラム)を算出する。そして、算出された位相変調パターンを含む制御信号SCが、SLM14に提供され、SLM14は、制御信号SCに基づいて制御される。
図7は、強度スペクトル設計部23の内部構成を示すブロック図である。図7に示されるように、強度スペクトル設計部23は、初期値設定部25、評価値算出部26、個体選定部27、及び次世代生成部28を含む。初期値設定部25は初期個体生成部25aを含む。また、図8は、変調パターン算出装置20による強度スペクトル設計方法(データ作成方法)を示すフローチャートである。以下、図7及び図8を参照しながら、本実施形態の変調パターン算出装置20の動作、すなわち強度スペクトル設計方法(データ作成方法)について説明する。
まず、強度スペクトル設計部23が、任意波形入力部21から入力された所望の時間強度波形に適した強度スペクトル関数A(ω)を生成する(強度スペクトル関数生成ステップS1)。詳細には、強度スペクトル関数生成ステップS1は、初期値設定ステップS11、評価値算出ステップS12、個体選定ステップS13、および次世代生成ステップS14を含んで構成されている。
初期値設定ステップS11では、初期値設定部25が、強度スペクトル関数A(ω)に関する第1世代のM個(但しMは2以上の整数)の個体(遺伝情報)A
1(ω)~A
M(ω)、及び位相スペクトル関数Ψ(ω)を設定する。個体A
1(ω)~A
M(ω)及び位相スペクトル関数Ψ(ω)は、周波数ωの関数である。位相スペクトル関数Ψ(ω)は、操作者により入力されたものであってもよく、或いは、位相スペクトル設計部22により計算されたものであってもよい。この初期値設定ステップS11により、強度スペクトル関数A(ω)の第1世代の個体A
1(ω)~A
M(ω)のそれぞれと、位相スペクトル関数Ψ
0(ω)とを含む周波数領域のM個の波形関数(1)が定義される。この波形関数(1)は、本実施形態における第1波形関数である。但し、iは虚数である。
本実施形態の初期値設定ステップS11は、初期個体生成ステップS11aを含んでいる。初期個体生成ステップS11aでは、初期個体生成部25aが、強度スペクトル関数A
IFTA(ω)を反復フーリエ法により生成し、この強度スペクトル関数A
IFTA(ω)を変化させることにより第1世代の個体A
1(ω)~A
M(ω)を生成する。図9は、初期個体生成部25aにおける強度スペクトル関数A
IFTA(ω)の算出手法を概念的に示す図である。図9に示されるように、まず、初期個体生成部25aは、初期の強度スペクトル関数A
k=0(ω)及び位相スペクトル関数Ψ
0(ω)を用意する(図中の処理番号(1))。一例では、初期の強度スペクトル関数A
k=0(ω)及び位相スペクトル関数Ψ
0(ω)は、入力光Laのスペクトル強度及びスペクトル位相に基づいて定められる。次に、初期個体生成部25aは、強度スペクトル関数A
k(ω)及び位相スペクトル関数Ψ
0(ω)を含む周波数領域の波形関数(2)を用意する(図中の処理番号(2))。この波形関数(2)は、本実施形態における第3波形関数である。
添え字kは、第k回目のフーリエ変換処理後を表す。最初(第1回目)のフーリエ変換処理の前においては、強度スペクトル関数A
k(ω)として上記の初期強度スペクトル関数A
k=0(ω)が用いられる。iは虚数である。
続いて、初期個体生成部25aは、上記関数(2)に対して周波数領域から時間領域へのフーリエ変換を行う(図中の矢印A1)。これにより、時間強度波形関数b
k(t)及び時間位相関数Θ
k(t)を含む時間領域の波形関数(3)が得られる(図中の処理番号(3))。
続いて、初期個体生成部25aは、フーリエ変換後の波形関数b
k(t)と関数Target
0(t)に係数αを乗じたもの(α×Target
0(t))との差が、波形関数b
k(t)と関数Target
0(t)との差よりも小さくなるような係数αを求める(図中の処理番号(4))。一例では、次の数式(4)で示されるように、フーリエ変換後の波形関数b
k(t)に対する、α×Target
0(t)の標準偏差σが最小(σ
min)となる係数αを探査的に導出する。なお、数式(4)において、Dはデータ点数を表し、t
s、t
eはそれぞれ時間軸の始点及び終点を表す。
続いて、初期個体生成部25aは、フーリエ変換後の関数(3)に含まれる時間強度波形関数b
k(t)に対して所望の波形に基づく置き換えを行う(第1の置き換え)。このとき、初期個体生成部25aは、所望の波形を表す関数Target
0(t)に係数αを乗じたもの(α×Target
0(t))を使用して置き換えを行う。一例では、数式(5)により算出されるTarget
k(t)に置き換える(図中の処理番号(5)、(6))。
続いて、初期個体生成部25aは、上記関数(6)に対して時間領域から周波数領域への逆フーリエ変換を行う(図中の矢印A2)。これにより、強度スペクトル関数C
k(ω)及び位相スペクトル関数Ψ
k(ω)を含む周波数領域の波形関数(7)が得られる(図中の処理番号(7))。
続いて、初期個体生成部25aは、上記関数(7)に含まれる位相スペクトル関数Ψ
k(ω)を拘束するため、初期の位相スペクトル関数Ψ
0(ω)に置き換える(第2の置き換え、図中の処理番号(8))。
また、初期個体生成部25aは、逆フーリエ変換後の周波数領域における強度スペクトル関数C
k(ω)に対し、入力光Laの強度スペクトルに基づくフィルタ処理を行う。具体的には、強度スペクトル関数C
k(ω)により表される強度スペクトルのうち、入力光Laの強度スペクトルに基づいて定められる各波長毎のカットオフ強度を超える部分をカットする。一例では、波長毎のカットオフ強度は、入力光Laの強度スペクトル(本実施形態では初期強度スペクトル関数A
k=0(ω))と一致するように設定される。その場合、次の数式(9)に示されるように、強度スペクトル関数C
k(ω)が初期強度スペクトル関数A
k=0(ω)よりも大きい周波数では、強度スペクトル関数A
k(ω)の値として初期強度スペクトル関数A
k=0(ω)の値が取り入れられる。また、強度スペクトル関数C
k(ω)が初期強度スペクトル関数A
k=0(ω)以下である周波数では、強度スペクトル関数A
k(ω)の値として強度スペクトル関数C
k(ω)の値が取り入れられる。
初期個体生成部25aは、上記関数(7)に含まれる強度スペクトル関数C
k(ω)を、上記数式(9)によるフィルタ処理後の強度スペクトル関数A
k(ω)に置き換える。また、C
k(ω)に任意の係数を乗じた関数C’
k(ω)を定義して、カットオフ強度を相対的に変化させる方法を用いても良い(図中の処理番号(9))。
以降、初期個体生成部25aが上記の処理(1)~(9)を複数回繰り返し行うことにより、波形関数中の強度スペクトル関数Ak(ω)を、所望の時間強度波形に対応する強度スペクトル形状に近づけることができる。最終的に、強度スペクトル関数AIFTA(ω)が得られる。
初期個体生成部25aは、強度スペクトル関数A
IFTA(ω)を変化させることにより第1世代の個体A
1(ω)~A
M(ω)を生成する。具体的には、初期個体生成部25aは、下記の数式(10)を用いて第1世代の個体A
1(ω)~A
M(ω)を生成する。
ここで、mは個体番号を表し、m=1,2,・・・,Mである。B
m(ω)は強度スペクトル関数A
IFTA(ω)に変化を与える確率関数である。この確率関数B
m(ω)を適切に設定することにより、強度スペクトル関数A
IFTA(ω)に準じた適切な第1世代の個体A
1(ω)~A
M(ω)を生成することができる。
例えば、被変調光である入力光Laの強度スペクトル関数A
pulse(ω)をA
IFTA(ω)×B
m(ω)が超えることはできない。確率関数B
m(ω)が取り得る実数範囲の上限を関数b
m(ω)、下限を関数a
m(ω)として表現するとき、これらの関数b
m(ω)、a
m(ω)は下記のように強度スペクトル関数A
IFTA(ω)及びA
pulse(ω)を用いて表すことができる。なお、式中のsは、分母が0にならないように便宜上挿入された微小値である。
図10は、入力光Laの強度スペクトル関数Apulse(ω)、及び強度スペクトル関数AIFTA(ω)の一例を示すグラフである。図10において、グラフG31は強度スペクトル関数Apulse(ω)を表し、グラフG32は強度スペクトル関数AIFTA(ω)を表す。また、横軸は波長(単位:nm)を表し、縦軸は強度(任意単位)を表す。横軸に関しては、波長を周波数ωに変換し、本文中の式や図で表現されているように周波数ωとして取り扱うことができる。図10に示されるように、この強度スペクトル関数Apulse(ω)はガウス分布に従う。また、強度スペクトル関数AIFTA(ω)は、極大及び極小を繰り返しながらも、前述した数式(9)の作用により全周波数域において強度スペクトル関数Apulse(ω)を超えない現実的な強度スペクトル関数となっている。
図11は、図10に示された強度スペクトル関数Apulse(ω)及びAIFTA(ω)を用いる場合の関数bm(ω)と、関数bm(ω)を上限、関数am(ω)(=0)を下限とする範囲内においてランダムに(全く無秩序に、且つ出現確率が同じになるように)生成された確率関数Bm(ω)とを例示するグラフである。図11において、グラフG41は関数bm(ω)を表し、グラフG42は確率関数Bm(ω)を表す。また、横軸は波長(単位:nm)を表し、縦軸は関数bm(ω)及び確率関数Bm(ω)の値(実数)を表す。図11に示されるように確率関数Bm(ω)の上限及び下限を適切に設定することにより、前述した数式(10)を用いて、現実的に取り得る実数値の範囲内でランダムに第1世代の個体A1(ω)~AM(ω)を生成することができる。また、上限と下限との差、すなわち関数bm(ω)と関数am(ω)との差を任意に狭めることにより、第1世代の個体A1(ω)~AM(ω)の分散を小さくすることができる。なお、確率関数Bm(ω)は完全な意味でのランダムでなくてもよく、例えば出現確率が正規分布に従ってもよい。
再び図7及び図8を参照する。次に、評価値算出ステップS12において、評価値算出部26は、第n世代(nは1以上の整数)の個体A
1(ω)~A
M(ω)のそれぞれと、共通の位相スペクトル関数Ψ(ω)とを含む周波数領域のM個の第1波形関数(12)
を、時間強度波形関数I
1(t)~I
M(t)それぞれと、時間位相波形関数Φ
1(t)~Φ
M(t)それぞれとを含む時間領域のM個の波形関数(13)
に変換する。これらの波形関数(13)は、本実施形態における第2波形関数である。そして、評価値算出部26は、時間強度波形関数I
1(t)~I
M(t)それぞれと所望の時間強度波形T(t)(=Target
0(t))との相違の度合いを示すM個の評価値を算出する。例えば、評価値算出部26は、所望の時間強度波形T(t)に対する時間強度波形関数I
1(t)~I
M(t)それぞれの標準偏差を評価値として算出する。このとき、所望の時間強度波形T(t)と時間強度波形関数I
1(t)~I
M(t)との間にエネルギー差が存在すると、このエネルギー差に起因して評価値が変動してしまう。本実施形態では、このエネルギー差を補償するために、探査型評価関数を導入する。具体的には、評価値算出部26は、次の数式(14)で表されるように、時間強度波形関数I
1(t)~I
M(t)それぞれと、所望の時間位相波形を表す関数T(t)に係数α
1~α
Mそれぞれを乗じたものとの相違の度合いを示すM個の評価値を算出する。
係数α
1~α
Mは、係数α
1~α
Mの乗算前と比較して、各評価値が良好になる値を有する。数式(14)は、評価値の一例として、所望の時間位相波形を表す関数T(t)に係数α
1~α
Mを乗じたものに対する、時間強度波形関数I
1(t)~I
M(t)の標準偏差を示している。この例では、各標準偏差が最小値をとるように、α
1~α
Mを変化させる。そして、標準偏差の最小値σ
min1~σ
minMそれぞれを、時間強度波形関数I
1(t)~I
M(t)それぞれの評価値とする。
続いて、評価値算出ステップS12において算出したM個の評価値(具体的には標準偏差の最小値σmin1~σminM)に基づいて、個体選定部27は、第n世代の複数の個体A1(ω)~AM(ω)の中から、第(n+1)世代の複数の個体A1(ω)~AM(ω)の生成に用いられる二以上の個体を選定する(個体選定ステップS13)。この個体選定ステップS13では、M個の評価値の優良さに基づいて、二以上の個体を選定する。ここで、「優良さに基づいて」とは、例えば、第n世代のM個の個体A1(ω)~AM(ω)から選定された少なくとも1つの個体からなる個体群G1(第1の個体群)が、M個の個体A1(ω)~AM(ω)のうちその個体群G1に含まれない他の全ての個体よりも評価値が優れていることを意味する。或いは、第n世代のM個の個体A1(ω)~AM(ω)から選定された1つ以上の個体からなる個体群G1における評価値の平均が、M個の個体A1(ω)~AM(ω)の評価値の平均よりも優れていることを意味してもよい。以下、この個体群G1を「エリート個体群」と称することがある。
また、本実施形態では、個体選定ステップS13において個体選定部27が選定する二以上の個体は、エリート個体群G1に加えて、少なくとも1つの別の個体からなる個体群G2(第2の個体群)を含んでもよい。この場合、個体群G2の評価値の平均は、第n世代のM個の個体A1(ω)~AM(ω)の評価値σmin1~σminMの平均よりも劣っている。以下、この個体群G2を「非エリート個体群」と称することがある。M個の個体A1(ω)~AM(ω)の評価値を優良な順に並べた場合、非エリート個体群G2の評価値はエリート個体群G1の評価値から連続していない。すなわち、M個の個体A1(ω)~AM(ω)の中には、エリート個体群G1の中で最も劣る評価値よりも劣り、且つ、非エリート個体群G2の中で最も優良な評価値よりも優良な評価値を有する1つ以上の個体が存在する。
続いて、次世代生成ステップS14において、次世代生成部28は、個体選定ステップS13において個体選定部27が選定した二以上の個体に基づいて、第(n+1)世代のM個の個体A1(ω)~AM(ω)を生成する。ここで、「選定された二以上の個体に基づいて第(n+1)世代の複数の個体を生成する」とは、例えば交叉、突然変異、増殖といった処理を意味し、第(n+1)世代のM個の個体A1(ω)~AM(ω)それぞれが、第n世代のいずれかの個体の少なくとも一部の成分を含んでいることを意味する。なお、選定した二以上の個体の一部(例えば最も優良な個体群)を、そのまま第(n+1)世代の個体A1(ω)~AM(ω)のいずれかとしてもよい。
強度スペクトル関数生成ステップS1においては、上述した評価値算出ステップS12、個体選定ステップS13、及び次世代生成ステップS14が、所定の条件が満たされるまでnを1ずつ加算しながら繰り返される(ステップS15)。言い換えると、評価値算出部26、個体選定部27、及び次世代生成部28は、所定の条件が満たされるまでnを1ずつ加算しながら処理を繰り返す。そして、強度スペクトル設計部23は(強度スペクトル関数生成ステップS1においては)、所定の条件が満たされた場合の第n世代のM個の個体A1(ω)~AM(ω)から、所望の時間強度波形に適した強度スペクトル関数A(ω)を生成する。例えば、M個の個体A1(ω)~AM(ω)のうち1つの個体Am(ω)を抜き出して、強度スペクトル関数A(ω)としてもよい。なお、所定の条件とは、例えば任意に設定した繰り返し試行回数を終了した場合や、任意に定めた評価値を満たした場合である。
上記の処理の後、データ生成ステップS2において、変調パターン生成部24は、位相スペクトル関数Ψ(ω)と、強度スペクトル関数生成ステップS1において生成された強度スペクトル関数A(ω)とに基づいて、SLM14に呈示させる変調パターンに関するデータを生成する。変調パターン生成部24は、生成したデータを、制御信号SCとしてSLM14に提供する。
以上に説明した本実施形態の光制御装置1A、変調パターン算出装置20、変調パターン算出方法、及び変調パターン算出プログラムによって得られる効果について説明する。従来、所望の時間波形を有する光をSLMを用いて実現する際、所望の時間波形に対応するスペクトル強度を精度良く算出するために、反復フーリエ法若しくは反復フーリエ法を修正した方法(例えば特許文献1,2を参照)が用いられている。しかしながら、同方法を用いてマルチパルスなどの生成を試みると、波形制御精度が大きく向上するものの、詳細に波形形状の分析を行うと、各パルスのピーク値やパルス幅に分散(ばらつき)が確認された。このことは、波形制御パターンの設計手法として改善の余地があることを意味する。特に、パルス光の顕微鏡応用や加工応用を考えた場合、パルス幅の変化やピーク値の変化は、信号のS/N比や加工状態の変化に大きく影響を及ぼす可能性がある。従って、波形制御パターンをより高精度に設計できる手法が望まれる。
このような課題に対し、本実施形態の変調パターン算出装置20、変調パターン算出方法、及び変調パターン算出プログラムにおいては、時間強度波形関数I1(t)~IM(t)と、所望の時間強度波形T(t)との相違の度合いを示す評価値の優良さに基づいて、次世代のM個の個体A1(ω)~AM(ω)の生成に用いられる二以上の個体を選定する。そして、選定された二以上の個体に基づいて、第(n+1)世代のM個の個体A1(ω)~AM(ω)を生成する。このような処理を、所定の条件が満たされるまでnを1ずつ加算しながら繰り返し、所定の条件が満たされた場合の第n世代のM個の個体A1(ω)~AM(ω)から、所望の時間強度波形T(t)に適した強度スペクトル関数A(ω)を生成する。本発明者は研究の末、このような方式(遺伝的アルゴリズム)により、反復フーリエ法及び反復フーリエ法を修正した方法と比較して、強度スペクトル関数A(ω)が局所解に導かれてしまう割合を低減し、最適解をより正確に探索することができることを見出した。
図12の(a)は、一例として、本実施形態の変調パターン算出装置20及び変調パターン算出方法によって算出された、8パルスの出力光Ldを生成するための強度スペクトル関数A(ω)を示すグラフである。グラフG61は、その強度スペクトル関数A(ω)を示す。また、グラフG62は、この例で用いられた位相スペクトル関数Ψ(ω)を示す。この位相スペクトル関数Ψ(ω)は、反復フーリエ法により算出されたものである。図12の(b)は、図12の(a)に示される強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)に基づく変調パターンをSLM14に呈示させて得られる出力光Ldの時間強度波形を示すグラフである。また、図13の(a)は、比較例として、反復フーリエ法のみを用いて算出された、8パルスの出力光Ldを生成するための強度スペクトル関数A(ω)を示すグラフである。グラフG71は、その強度スペクトル関数A(ω)を示す。また、グラフG72は、この比較例で用いられた位相スペクトル関数Ψ(ω)を示す。なお、この位相スペクトル関数Ψ(ω)の形状は図12の(a)のグラフG62と同一である。図13の(b)は、図13の(a)に示される強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)に基づく変調パターンをSLM14に呈示させて得られる出力光Ldの時間強度波形を示すグラフである。
図13の(b)を参照すると、反復フーリエ法のみを用いて強度スペクトル関数A(ω)を算出した場合、生成される出力光Ldの時間強度波形において、8つのパルスのピーク強度に大きなばらつきが見られる。これに対し、図12の(b)を参照すると、本実施形態の変調パターン算出装置20及び変調パターン算出方法によって強度スペクトル関数A(ω)を算出した場合、強度損失量は39%となり反復フーリエ法の31%よりも高くなるものの、生成される出力光Ldの時間強度波形において、8つのパルスのピーク強度のばらつきは小さく抑えられている。評価値に関しては、本実施形態による値が反復フーリエ法による値の1/8.7となり、大きく改善された。
このように、本実施形態の変調パターン算出装置20及び変調パターン算出方法によれば、反復フーリエ法のみを用いる従来の装置及び方法と比較して、局所解に導かれてしまう割合を低減し、最適解をより正確に探索することができる。すなわち、本実施形態によれば、出力光Ldの時間波形を所望の波形T(t)に近づけるためのスペクトル強度を精度良く算出して、所望の時間波形を精度良く得ることができる。なお、本実施形態では各世代の個体の数を全てM個に統一しているが、各世代の個体の数は変化してもよい。
ここで、本実施形態の有効性を確認するために、マルチパルスを含む時間強度波形を有する出力光Ldを生成するための変調パターンを、パルス数を変えながら複数計算した。各パルスはTLパルス(時間幅135fsのシングルパルス)であり、パルス間隔は1psで等間隔であった。初期位相スペクトルΨ0(ω)としては、反復フーリエ法を用いて算出したものを使用した。図14は、その際のパルス数と評価値(数式(14)に示される標準偏差の最小値)との関係をプロットしたグラフである。図14において、プロットP11は、反復フーリエ法のみを用いる従来の方法により強度スペクトル関数A(ω)を算出した場合を示し、プロットP12は、本実施形態の方法により強度スペクトル関数A(ω)を算出した場合を示す。なお、本実施形態の方法による強度スペクトル関数A(ω)の算出では、評価値が収束するまでステップS12~S14を十分に繰り返した。その繰り返し数(世代数)は1700世代であった。図14に示されるように、いずれのパルス数においても、本実施形態の方法によれば、従来の方法と比較して波形制御精度(評価値)が大幅に改善される。パルス数によって差はあるが、概ね3倍から24倍程度の改善が見込めることが確認された。
また、図15は、上記の実施例におけるパルス数と平均パルス幅(半値全幅)との関係をプロットしたグラフである。図16は、上記の実施例におけるパルス数とピーク値の分散との関係をプロットしたグラフである。図15及び図16において、プロットP21,P31は、反復フーリエ法のみを用いる従来の方法により強度スペクトル関数A(ω)を算出した場合を示し、プロットP22,P32は、本実施形態の方法により強度スペクトル関数A(ω)を算出した場合を示す。これらの図から、図14に示された評価値の改善は、パルス幅が狭くなったこと(TLパルスの半値幅に近くなったこと)、及び、ピーク値のばらつきが小さくなったことに因ると考えられる。
本実施形態のように、初期値設定部25(初期値設定ステップS11)は、第1世代の個体A1(ω)~AM(ω)を生成する初期個体生成部25a(初期個体生成ステップS11a)を含んでもよい。そして、初期個体生成部25a(初期個体生成ステップS11a)は、反復フーリエ変換によって強度スペクトル関数AIFTA(ω)を生成し、該強度スペクトル関数AIFTA(ω)を変化させることにより第1世代の個体A1(ω)~AM(ω)を生成してもよい。本発明者の知見によれば、本実施形態の変調パターン算出装置20において最適解をより正確に探索するためには、第1世代の個体A1(ω)~AM(ω)の設定が極めて重要である。そして、反復フーリエ法には、短時間で評価値の優れた解を算出できるという特徴がある。その上、その解近傍にさらに評価値の高い解が、存在することがある。従って、反復フーリエ法を用いて第1世代の個体A1(ω)~AM(ω)の基となる強度スペクトル関数AIFTA(ω)を生成することにより、第1世代の個体A1(ω)~AM(ω)を適切に設定することができる。
また、本実施形態のように、評価値算出部26(評価値算出ステップS12)は、時間強度波形関数I1(t)~IM(t)それぞれと、所望の時間位相波形を表す関数T(t)に係数α1~αMそれぞれを乗じたものとの相違の度合いを示すM個の評価値を算出し、係数α1~αMは、該係数α1~αMの乗算前と比較して、乗算後の評価値が良好になる値を有してもよい。これにより、所望の時間強度波形T(t)と時間強度波形関数I1(t)~IM(t)との総エネルギーの違いが評価値の算出に影響することを抑制し、所望の時間強度波形T(t)と時間強度波形関数I1(t)~IM(t)との形状の違いに主に基づいて評価値を算出することができる。
また、本実施形態のように、個体選定部27(個体選定ステップS13)において選定される二以上の個体は、少なくとも一つの個体からなる個体群G1と、少なくとも一つの別の個体からなる個体群G2とを含んでもよい。そして、個体群G1の評価値の平均は、第n世代の個体A1(ω)~AM(ω)の評価値の平均よりも優れており、個体群G2の評価値の平均は、第n世代の個体A1(ω)~AM(ω)の評価値の平均よりも劣っていてもよい。
本実施形態のような遺伝的アルゴリズムにおいては、世代数が進んでいくと、複数の個体A1(ω)~AM(ω)が次第に均一に近づく。故に、評価値が収束し、それ以上の改善が見られなくなるか、または、改善の度合いが著しく低下してしまう。そこで、本実施形態では、個体選定部27(個体選定ステップS13)において選定される二以上の個体に、非エリート個体からなる個体群G2を含める。図17は、8パルスの出力光Ldを生成するための強度スペクトル関数A(ω)(図12の(a)を参照)を算出するときに、世代が進むに従い評価値が変化する様子を示すグラフである。図17の横軸は世代を表し、縦軸は評価値(標準偏差の最小値σmin)を表す。図17に示される例では、約150世代の辺りで評価値が一旦収束しているが、約250世代の辺りで非エリート個体を導入(図中に矢印A3で示す)したことにより、評価値が不連続的に低下(改善)している。そして、その後に非エリート個体を導入する毎に、更に評価値が不連続的に低下(改善)している。このように、個体選定部27(個体選定ステップS13)において選定される二以上の個体に非エリート個体を導入することによって、一旦収束した値から抜け出し、評価値を更に良好にすることが可能となる。
本発明者は、更に、どのような非エリート個体が効果的かを調べた。図18は、非エリート個体の導入前後での評価値の変化量を図示したものである。この実施例では、個体群G2を構成する非エリート個体を評価値の範囲別に4つのグループ(0.0025<σmin<0.003、0.0035<σmin<0.004、0.0045<σmin<0.005、及び0.0055<σmin<0.006)に分けて、それぞれ個別に計算を行った。非エリート個体の導入は、評価値が概ね0.0015から0.0025で収束値を迎えた演算に対して行った。図18において、グラフG51は各範囲毎における計算を複数回試行したときの評価値の変化量の平均を示し、グラフG52はそのときの各範囲における評価値の変化量の存在範囲を示す。図18を参照すると、非エリート個体の評価値にかかわらず、評価値が収束値から更に優良な方向へ有意に変化した。また、この計算条件においては、非エリート個体の評価値が優れている(標準偏差が小さい)ほど、評価値が収束値から優良な方向へ大きく変化し、顕著な効果が得られることがわかった。
このように、本発明者の試行によれば、評価値が比較的劣る非エリート個体を選定後の個体の一部に含めることにより、最終的なスペクトル強度の算出精度をより高めることができる。
また、本実施形態の光制御装置1Aによれば、変調パターン算出装置20を備えることにより、局所解に導かれる割合を低減してスペクトル強度を精度良く算出し、出力光Ldの時間波形を所望の波形T(t)に近づけることができる。
なお、上記の説明では主に強度スペクトル設計部23の構成及びスペクトル強度の算出方法について説明したが、位相スペクトル設計部22の構成及びスペクトル位相の算出方法は、従来の構成及び方法(例えば反復フーリエ法若しくはその改良方法)を用いてもよく、或いは、本実施形態の強度スペクトル設計部23の構成及びスペクトル強度の算出方法と同様の構成及び方法を用いてもよい。
(第2実施形態)
図19は、本発明の第2実施形態に係る光制御装置1Bの構成を概略的に示す図である。本実施形態の光制御装置1Bは、入力光Laから、時間間隔をあけて並ぶ複数の光パルスを含む光パルス列を出力光Ldとして生成する。光パルス列は、複数の光パルスLP1~LPN(Nは2以上の整数)を含む。各光パルスLP1~LPNは、例えば100fs以下といった極めて短い時間幅を有する。各光パルスLP1~LPNの出現タイミングは任意であり、隣接する直前の光パルスとの時間間隔は、各光パルスLP2~LPNにおいて互いに等しくてもよく、異なってもよい。また、各光パルスLP2~LPNのピーク強度は互いに等しくてもよく、異なってもよい。
光制御装置1Bは、光源2、光学系10、及び変調パターン算出装置(データ作成装置)30を備える。光源2及び光学系10の構成は、第1実施形態(図2を参照)と同様である。変調パターン算出装置30のハードウェア構成は、第1実施形態の変調パターン算出装置20(図4を参照)と同様である。変調パターン算出装置30は、SLM14と電気的に接続されており、出力光Ldの時間強度波形を所望のピーク強度及び時間間隔を有する光パルス列に近づけるための位相変調パターンを算出し、該位相変調パターンを含む制御信号SC2をSLM14に提供する。本実施形態の変調パターン算出装置30は、所望の光パルス列を得る為の位相スペクトルを出力光Ldに与える位相変調用の位相パターンと、所望の光パルス列を得る為の強度スペクトルを出力光Ldに与える強度変調用の位相パターンとを含む位相パターンをSLM14に呈示させる。そのために、変調パターン算出装置30は、任意波形入力部31と、スペクトル設計部32と、変調パターン生成部(データ生成部)33とを有する。すなわち、変調パターン算出装置30に設けられたコンピュータのプロセッサ201(図4を参照)は、任意波形入力部31の機能と、スペクトル設計部32の機能と、変調パターン生成部33の機能とを実現する。
任意波形入力部31は、操作者からの所望の光パルス列に関する情報の入力を受け付ける。操作者は、所望の光パルス列に関する情報(光パルス数、各光パルスのピーク強度、時間間隔など)を任意波形入力部31に入力する。所望の光パルス列に関する情報は、スペクトル設計部32に与えられる。スペクトル設計部32は、与えられた所望の光パルス列の実現に適した、出力光Ldの位相スペクトル及び強度スペクトルを算出する。変調パターン生成部33は、スペクトル設計部32において求められた位相スペクトル及び強度スペクトルを出力光Ldに与えるための位相変調パターン(例えば、計算機合成ホログラム)を算出する。そして、算出された位相変調パターンを含む制御信号SC2が、SLM14に提供され、SLM14は、制御信号SC2に基づいて制御される。
図20は、スペクトル設計部32の内部構成を示すブロック図である。図20に示されるように、スペクトル設計部32は、初期値設定部35、評価値算出部36、個体選定部37、及び次世代生成部38を含む。初期値設定部35は初期個体生成部39を含む。また、図21は、変調パターン算出装置30による強度スペクトル及び位相スペクトルの設計方法(データ作成方法)を示すフローチャートである。以下、図20及び図21を参照しながら、本実施形態の変調パターン算出装置30の動作、すなわち強度スペクトル及び位相スペクトルの設計方法(データ作成方法)について説明する。
まず、スペクトル設計部32が、任意波形入力部31から入力された所望の光パルス列に関する情報に基づいて、該光パルス列に適した強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を生成する(スペクトル関数生成ステップS3)。詳細には、スペクトル関数生成ステップS3は、初期値設定ステップS31、評価値算出ステップS32、個体選定ステップS33、および次世代生成ステップS34を含んで構成されている。
初期値設定ステップS31では、初期値設定部35が、強度スペクトル関数A(ω)に関する第1世代のM個(但しMは2以上の整数)の個体(遺伝情報)A
1(ω)~A
M(ω)、及び位相スペクトル関数Ψ(ω)に関する第1世代のM個の個体(遺伝情報)Ψ
1(ω)~Ψ
M(ω)を設定する。各個体A
m(ω)と各個体Ψ
m(ω)(但し、m=1,2,・・・,M)とは、第m番目の個体対(ペア)を構成する。すなわち、個体A
1(ω)~A
M(ω)と個体Ψ
1(ω)~Ψ
M(ω)とは、M個の個体対を構成する。個体A
1(ω)~A
M(ω)及び個体Ψ
1(ω)~Ψ
M(ω)は、周波数ωの関数である。この初期値設定ステップS31により、強度スペクトル関数A(ω)の第1世代の個体A
1(ω)~A
M(ω)のそれぞれと、位相スペクトル関数Ψ
0(ω)の第1世代の個体Ψ
1(ω)~Ψ
M(ω)のそれぞれとを含む周波数領域のM個の波形関数(15)が定義される。但し、iは虚数である。
本実施形態の初期値設定ステップS31は、初期個体生成ステップS36を含んでいる。初期個体生成ステップS36では、初期個体生成部39が、所望の光パルス列に含まれる各光パルスの振幅及びタイミングと同じ振幅及びタイミングを有し、時間位相がそれぞれ異なるM個のデルタ関数群をそれぞれフーリエ変換することにより、第1世代のM個の個体対(Am(ω),Ψm(ω))(但しm=1,2,・・・,M)を生成する。
図22は、初期個体生成部25aにおける、M個のデルタ関数群から第1世代のM個の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)を算出する方法の一例を概念的に示す図である。図22の(a)部は所望の光パルス列の時間波形を示す。図22の(b)部ないし(d)部はそれぞれ第1番目、第2番目および第M番目のデルタ関数群を示す。なお、図22各部の上部には複数の光パルスの番号(1,2,・・・,N)が示されており、図中のt1,t2,・・・,tNは各光パルスの出現タイミングを示す時間であり、A1,A2,・・・,ANは各光パルスのピーク強度である。また、図23は、初期個体生成ステップS36の詳細を示すフローチャートである。図23に示されるように、初期個体生成ステップS36は、ステップS36a~S36cを含む。
図22の(a)に示される所望の光パルス列の時間波形を実現するために、初期個体生成部25aは、まず、該所望の光パルス列のパラメータの入力を受ける(ステップS36a)。光パルス列のパラメータとは、例えば当該光パルス列に含まれる光パルスの本数N、各光パルスのピーク強度(A1,A2,・・・,AN)、各光パルスの出現タイミング(t1,t2,・・・,tN)等である。
次に、初期個体生成ステップS36は、光パルス列に関する上記パラメータを用いて、図22の(b)~(d)に示されるようなM個のデルタ関数群を用意する(ステップS36b)。各デルタ関数群のN個の光パルスは、所望の光パルス列に含まれるN個の光パルスの振幅及びタイミングと同じ振幅及びタイミングを有する。また、各デルタ関数群のN個の光パルスの時間位相は任意であり、それぞれ異なる。時間位相とは、前述した数式(3)のΘk(t)に相当する値である。具体的には、図22の(b)~(d)に示されるように、デルタ関数群に含まれる各デルタ関数の位相Φ1~ΦNの大きさを設定する。
続いて、初期個体生成部25aは、これらのデルタ関数群のそれぞれに対してフーリエ変換を行う(ステップS36c)。これにより、数式(15)に示された周波数領域のM個の波形関数が求まり、その結果、第1世代の個体A1(ω)~AM(ω),Ψ1(ω)~ΨM(ω)が算出される。
再び図20及び図21を参照する。次に、評価値算出ステップS32において、評価値算出部36は、第n世代(nは1以上の整数)の個体A
1(ω)~A
M(ω)それぞれを、最大値が均等となるように規格化する。一例では、評価値算出部36は、第n世代の個体A
1(ω)~A
M(ω)それぞれを、最大値が1となるように規格化する。そして、評価値算出部36は、個体A
1(ω)~A
M(ω)それぞれの積分値に基づいて評価値を算出する。例えば、評価値算出部36は、個体A
1(ω)~A
M(ω)に基づく強度スペクトル変調によって生じる損失量を示すM個の評価値を算出してもよい。その場合、評価値算出部26は、規格化後の個体A
1(ω)~A
M(ω)の積分値(グラフの面積)Area
1~Area
Mを、全ての周波数で上記の均等値h(例えばh=1)をとる強度スペクトル関数A(ω)=hの面積Areaから差し引き、その値をAreaで除算した値(すなわち下記の数式(16)に示される値)を、評価値Loss
1~Loss
Mとして算出してもよい。
続いて、評価値算出ステップS32において算出したM個の評価値(例えばLoss1~LossM)に基づいて、個体選定部37は、第n世代のM個の個体対(Am(ω),Ψm(ω))(但しm=1,2,・・・,M)の中から、第(n+1)世代のM個の個体対(Am(ω),Ψm(ω))の生成に用いられる二以上の個体対(Am(ω),Ψm(ω))を選定する(個体選定ステップS33)。この個体選定ステップS33では、M個の評価値の優良さに基づいて、二以上の個体対を選定する。ここで、「優良さに基づいて」とは、例えば、第n世代のM個の個体対(Am(ω),Ψm(ω))から選定された少なくとも1つの個体対(Am(ω),Ψm(ω))からなる個体対群G3(第1の個体対群)が、M個の個体対(Am(ω),Ψm(ω))のうちその個体対群G3に含まれない他の全ての個体対よりも評価値が優れていることを意味する。或いは、第n世代のM個の個体対(Am(ω),Ψm(ω))から選定された1つ以上の個体対からなる個体対群G3における評価値の平均が、M個の個体対(Am(ω),Ψm(ω))の評価値の平均よりも優れていることを意味してもよい。以下、この個体対群G3を「エリート個体対群」と称することがある。
また、本実施形態では、個体選定ステップS33において個体選定部37が選定する二以上の個体対(Am(ω),Ψm(ω))は、エリート個体対群G3に加えて、少なくとも1つの別の個体対からなる個体対群G4(第2の個体対群)を含んでもよい。この場合、個体対群G4の評価値の平均は、第n世代のM個の個体対(Am(ω),Ψm(ω))の評価値の平均(例えばLoss1~LossMの平均)よりも劣っている。以下、この個体対群G4を「非エリート個体対群」と称することがある。M個の個体対(Am(ω),Ψm(ω))の評価値を優良な順に並べた場合、非エリート個体対群G4の評価値はエリート個体対群G3の評価値から連続していない。すなわち、M個の個体対(Am(ω),Ψm(ω))の中には、エリート個体対群G3の中で最も劣る評価値よりも劣り、且つ、非エリート個体対群G4の中で最も優良な評価値よりも優良な評価値を有する1つ以上の個体対が存在する。
続いて、次世代生成ステップS34において、次世代生成部38は、個体選定ステップS33において個体選定部37が選定した二以上の個体対(Am(ω),Ψm(ω))に基づいて、第(n+1)世代のM個の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)を生成する。ここで、「選定された二以上の個体対に基づいて第(n+1)世代の複数の個体対を生成する」とは、例えば交叉、突然変異、増殖といった処理を意味し、第(n+1)世代のM個の個体対(Am(ω),Ψm(ω))それぞれが、第n世代のいずれかの個体対(Am(ω),Ψm(ω))の少なくとも一部の成分を含んでいることを意味する。なお、選定した二以上の個体対(Am(ω),Ψm(ω))の一部(例えば最も優良な個体対群)を、そのまま第(n+1)世代の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)のいずれかとしてもよい。
スペクトル関数生成ステップS3においては、上述した評価値算出ステップS32、個体選定ステップS33、及び次世代生成ステップS34を、所定の条件が満たされるまでnを1ずつ加算しながら繰り返す(ステップS35)。言い換えると、評価値算出部36、個体選定部37、及び次世代生成部38は、所定の条件が満たされるまでnを1ずつ加算しながら処理を繰り返す。そして、スペクトル設計部32は(スペクトル関数生成ステップS3においては)、所定の条件が満たされた場合の第n世代のM個の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)から、所望の光パルス列の形態(光パルスの数、各光パルスのピーク強度、各光パルスの出現タイミング等)に適した強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を生成する。例えば、M個の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)のうち1つの個体対(Am(ω),Ψm(ω))(mは任意の1つの整数)を抜き出して、強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)としてもよい。なお、所定の条件とは、例えば任意に設定した繰り返し試行回数を終了した場合や、任意に定めた評価値を満たした場合である。
上記の処理の後、データ生成ステップS2において、変調パターン生成部33は、スペクトル関数生成ステップS3において生成された強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)に基づいて、SLM14に呈示させる変調パターンに関するデータを生成する。変調パターン生成部33は、生成したデータを、制御信号SC2としてSLM14に提供する。
以上に説明した本実施形態の光制御装置1B、変調パターン算出装置30、変調パターン算出方法、及び変調パターン算出プログラムによって得られる効果について説明する。本実施形態においては、強度スペクトル変調により生じる損失量を示す評価値の優良さに基づいて、次世代のM個の個体対(Am(ω),Ψm(ω))(m=1,2,・・・,M)の生成に用いられる二以上の個体対(Am(ω),Ψm(ω))を選定する。そして、選定された二以上の個体対(Am(ω),Ψm(ω))に基づいて、第(n+1)世代の複数の個体対(Am(ω),Ψm(ω))を生成する。このような処理を、所定の条件が満たされるまでnを1ずつ加算しながら繰り返し、所定の条件が満たされた場合の第n世代のM個の個体対(Am(ω),Ψm(ω))から、所望の光パルス列に適した強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を生成する。第1実施形態と同様に、このような方式(遺伝的アルゴリズム)により、反復フーリエ法及び反復フーリエ法を修正した方法と比較して、強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)が局所解に導かれてしまう割合を低減し、最適解をより正確に探索することができる。すなわち、本実施形態によれば、局所解に導かれてしまう割合を低減しつつ、光パルス列の時間波形を所望の波形に近づけるためのスペクトル強度A(ω)及びスペクトル位相Ψ(ω)を算出することができる。更に、本実施形態によれば、光パルス列の生成時に生じる損失量を最小化するためのスペクトル強度A(ω)及びスペクトル位相Ψ(ω)を精度良く算出することができる。
図24の(a)は、一例として、本実施形態の変調パターン算出装置30及び変調パターン算出方法によって算出された、50本の光パルスからなる光パルス列を生成するための強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を示すグラフである。グラフG81は、その強度スペクトル関数A(ω)を示す。グラフG82は、その位相スペクトル関数Ψ(ω)を示す。図24の(b)は、図24の(a)に示される強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)に基づく変調パターンをSLM14に呈示させて得られる出力光Ldの時間強度波形を示すグラフである。図24の(b)を参照すると、本実施形態の変調パターン算出装置30及び変調パターン算出方法によって強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を算出した場合、生成される出力光Ldの時間強度波形において、50本の光パルスのピーク強度のばらつきが小さく抑えられていることがわかる。
また、従来より、所望の形態(光パルスの数、各光パルスのピーク強度、各光パルスの出現タイミング等)を有する光パルス列をSLMを用いて実現する際、損失の少ない強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を算出するために、光パルス列の時間強度波形関数及び時間位相関数の最適な組み合わせを網羅的に探索する方法が用いられている。しかしながら、このような方法を用いて光パルス列の生成を試みると、光パルスの数の増加に伴って計算量が指数関数的に増大し、コンピュータの計算能力及び計算時間の制約等から、光パルスの数が実質的に制限されることが確認された。これに対し、本実施形態の方式によれば、各世代の個体対(Am(ω),Ψm(ω))の数がM個に限定されており、また、光パルスの数の増加に伴う計算量の増加は僅かである。従って、本実施形態によれば、所望の形態を有する光パルス列を実現するための強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)を、光パルスの数によらず設計することができる。
図25は、光パルス列の時間強度波形関数及び時間位相関数の最適な組み合わせを網羅的に探索する従来の方式において必要となる高速フーリエ変換(FFT)の演算回数と、本実施形態において必要となるFFTの演算回数を比較するグラフである。縦軸はFFTの演算回数を表し、横軸は光パルスの数を表す。また、円形のプロットP41は本実施形態の方式を示し、菱形のプロットP42は従来の方式を示す。図25に示されるように、従来の方式では、光パルス数の増加に伴いFFT回数が指数関数的に増加する。これに対し、本実施形態の方式では、FFT回数がパルス数に殆ど依存していない。このことは、本実施形態の方式が、所望の光パルス列のための変調パターンの算出を、光パルス数の制限なしに可能とすることを示している。
図26は、従来の方式により算出した強度スペクトル変調パターンによる損失量と、本実施形態の方式により算出した強度スペクトル変調パターンによる損失量とを比較するグラフである。縦軸は損失量(任意単位)を表し、横軸は光パルスの数を表す。また、円形のプロットP51は本実施形態の方式を示し、菱形のプロットP52は従来の方式を示す。従来の方式では、10パルスを超えると計算量が膨大になり、計算ができなかったため、10パルス以下の結果を示している。図26を参照すると、本実施形態の方式では、従来の方式と比較して、いずれのパルス数においても損失量が小さいことがわかる。すなわち、本実施形態によれば、計算量を少なく抑えつつ、損失が小さい変調パターンを作成することができる。
また、本実施形態のように、個体選定部37(個体選定ステップS33)において選定される二以上の個体対(Am(ω),Ψm(ω))は、少なくとも一つの個体対(Am(ω),Ψm(ω))からなる個体対群G3と、少なくとも一つの別の個体対(Am(ω),Ψm(ω))からなる個体対群G4とを含んでもよい。そして、個体対群G3の評価値の平均は、第n世代のM個の個体対(Am(ω),Ψm(ω))の評価値の平均よりも優れており、個体対群G4の評価値の平均は、第n世代のM個の個体対(Am(ω),Ψm(ω))の評価値の平均よりも劣っていてもよい。
第1実施形態と同様に、本実施形態においても、個体選定部37(個体選定ステップS33)において選定される二以上の個体対(Am(ω),Ψm(ω))に、非エリート個体対からなる個体対群G4を含める。図27は、50パルスの光パルス列を生成するための強度スペクトル関数A(ω)及び位相スペクトル関数Ψ(ω)(図23の(a)を参照)を算出するときに、世代が進むに従い評価値が変化する様子を示すグラフである。図27の横軸は世代を表し、縦軸は評価値(強度スペクトル変調による損失量)を表す。図27に示される例では、約40世代の辺りで評価値が一旦収束しているが、約100世代の辺りで非エリート個体対を導入(図中に矢印A4で示す)したことにより、評価値が不連続的に低下(改善)している。そして、その後に非エリート個体対を導入する毎に、更に評価値が不連続的に低下(改善)している。このように、個体選定部37(個体選定ステップS33)において選定される二以上の個体対(Am(ω),Ψm(ω))に非エリート個体対を導入することによって、一旦収束した値から抜け出し、評価値を更に良好にすることが可能となる。
また、本実施形態の光制御装置1Bによれば、変調パターン算出装置30を備えることにより、局所解に導かれる割合を低減しつつスペクトル強度A(ω)及びスペクトル位相Ψ(ω)を算出し、出力光Ldの時間波形を所望の光パルス列の時間波形に近づけることができる。更に、本実施形態の光制御装置1Bによれば、光パルス列の生成時に生じる損失量を最小化するためのスペクトル強度A(ω)及びスペクトル位相Ψ(ω)を精度良く算出することができる。
なお、本実施形態においては各世代の個体対の数を全てM個に統一しているが、各世代の個体対の数は変化してもよい。
(第3実施形態)
上記実施形態の変調パターン算出装置20、変調パターン算出方法、及び変調パターン算出プログラムは、時間パルス整形に代表される強度スペクトル変調パターン(1次元パターン)の設計に限らず、例えばビーム強度分布整形に代表される、2次元強度変調パターンの設計にも用いられ得る。言い換えると、例えばホログラムといった、所望の強度パターンと光学的フーリエ変換の関係にあるような領域にあるパターンの強度分布の設計にも用いられ得る。
図28は、本発明の第3実施形態に係る、2次元強度変調パターンを有効利用する際の光制御装置1Cの構成を概略的に示す図である。なお、図28では、光制御装置1Cが備える変調パターン算出装置20(または30)の図示を省略している。光源は、上記各実施形態の光源2のようにパルス光源であってもよく、或いはCW(Continuous Wave)レーザ光源であってもよい。光制御装置1Cは、上記各実施形態の応用として、所望の光強度分布をスクリーン45に表示する。光制御装置1Cは、2つのSLM41,46と、一対のレンズ42,43と、フーリエ変換レンズ44とを備えている。2つのSLM41,46は、一対のレンズ42,43を介して光学的に結合されている。SLM46とレンズ42との光学距離はレンズ42の焦点距離f1であり、SLM41とレンズ43との光学距離はレンズ43の焦点距離f2である。一例では、焦点距離f1と焦点距離f2とは互いに等しい。レンズ42とレンズ43との間の光学距離は、焦点距離f1と焦点距離f2との和である。フーリエ変換レンズ44は、SLM41と光学的に結合されており、その間の光学距離はフーリエ変換レンズ44の焦点距離f3である。この光制御装置1Cは、フーリエ変換レンズ44からSLM41とは反対側に焦点距離f3だけ離れたスクリーン45上に出力光像を結像する。SLM46は2次元の強度変調用のSLMであって、変調パターン算出装置20(または30)から提供された強度変調用の変調パターンを呈示する。SLM41は2次元の位相変調用のSLMであって、変調パターン算出装置20(または30)から提供された位相変調用の変調パターンを呈示する。強度変調用のホログラムパターンを表示するSLM46と、位相変調用のホログラムパターンを表示するSLM41とは、互いに入れ替えてもよい。
本発明によるデータ作成装置、光制御装置、データ作成方法、及びデータ作成プログラムは、上述した実施形態に限られるものではなく、他に様々な変形が可能である。例えば、上述した第1実施形態では初期値設定部が初期個体生成部を含み、初期個体生成部が反復フーリエ法を用いて第1世代の複数の個体を生成しているが、第1世代の複数の個体の決定方式はこれに限られず、例えば任意の複数の個体を入力してもよい。また、上述した第1実施形態では、評価値算出部が、第2波形関数の時間強度波形関数と、所望の時間位相波形を表す関数に係数を乗じたものとの相違の度合いを示す評価値を算出しているが(数式(14))、評価値の算出式はこれに限られず、第2波形関数の時間強度波形関数と所望の時間強度波形との相違の度合いを表すものであれば任意の算出式を用いることができる。