以下、実施形態を図面に基づいて説明する。
図1には、実施形態に係るスペクトル処理の流れが模式的に示されている。左上に示されたNMRスペクトル(以下、単に「スペクトル」という。)10が処理対象である。横軸は周波数fを示しており、縦軸はパワー又は強度を示している。後の説明ではスペクトル10はyで表現される。スペクトル10は、注目波形成分10aとベースライン成分10bとを含むものである。ベースライン成分は、注目波形成分10a以外の成分であって、概ね広帯域の成分であり、スペクトル10における底辺に相当するものである。スペクトル解析に際してはそのベースライン成分10bを事前に除去しておくことが望まれる。
符号12で示される処理において、スペクトル10の中から、ユーザーにより又は自動的に代表部分14が指定される。後の説明では代表部分14はyIで表現される。図示の例では、代表部分14は、複数の点として指定されている。それらは注目波形成分10aを避けて設定されており、換言すれば、ベースライン成分らしいところに集中的に設定されている。1又は複数の区間を指定することにより、代表部分14が指定されてもよい。
実施形態に係る方法では、ベースライン成分10bを模擬するベースラインモデル16が生成される。後の説明ではベースラインモデル16はAxで表現される。A及びxについては後述する。実際には、ベースライン成分の内の代表部分14に対して、ベースラインモデル中の対応部分がフィッティングされる。後の説明ではフィッティング部分はSIAxで表現される。SIについては後述する。
符号18で示された処理においては、第1条件及び第2条件が共に満たされるように、ベースラインモデル16が生成される。後に詳しく説明するように、第1条件は、代表部分14から対応部分を引いた残差のL2ノルムを小さくする条件である。第2条件は、以下に説明する重み列のLpノルムを小さくする条件である。ここで、0≦p≦1であり、解の安定性を求める場合には0<p≦1である。pが1以下の場合、Lpノルムを用いた解の探索過程において解のスパース性が増大することが知られている。実施形態においては、その性質が重み列の探索で利用されている。
実施形態においては、複数の基底関数からなる基底関数列20が用意されており、それに対して重み列22が乗算される。その乗算が符号24で示されている。1つの基底関数に対して1つの重みが与えられる。もっとも、重みの値が0であればそれが乗算される基底関数は事実上、機能しない。複数の基底関数に対して複数の重みを乗算した結果が加算され、これによってベースラインモデル16が生成される。通常、基底関数列20として、様々な次数をもった様々な種類の関数が用意されている。例えば、1から10次のN次対向式、及び、1から20次の余弦関数、等が用意されている。その他の関数が用意されてもよい。多数の基底関数を用意しておいても、過剰フィッティングのような問題が生じにくい。個々の関数は、インデックス(又は基底関数番号)sで識別される。
第1条件及び第2条件が共に満たされるように重み列22が最適化される。その過程において、特に第2条件に従って、重み列22のスパース性が徐々に高まっていく。図1に示した重み列22においては、一部の基底関数だけに対して有意な重みが与えられており、残りの多くの基底関数に対しては0又は0に近い重みが与えられている。基底関数を多くすると、フィッティング自由度が高まる一方、過剰フィッティングの問題が生じやすくなるが、Lpのノルムが発揮する解のスパース性を利用すれば、その問題を回避又は軽減すること可能となる。
重み列の探索結果が終了条件を満たした時点で探索が終了し、その時の重み列に基づいて最適ベースラインモデルが構成される。後の説明ではそれがx’で表現される。符号26で示された処理は、スペクトル10からの最適ベースラインモデルの減算を示している。これにより、ベースライン成分が取り除かれたスペクトル28が得られる。後の説明ではそれがy’で表現される。
図2には、NMR測定システムが示されている。NMR測定システムは、NMR測定装置30とスペクトル処理装置32とにより構成される。NMR測定装置30からスペクトル処理装置32へスペクトルデータが転送される。その転送は、例えば、ネットワークを介して又は記憶媒体を介して行われる。
NMR測定装置30は、図2においてはその詳細が示されていないものの、分光計及び測定部により構成される。測定部は、静磁場発生器、プローブ等により構成される。静磁場発生器は、垂直貫通孔としてのボアを有し、そのボアの内部にプローブの挿入部が挿入される。挿入部のヘッド内にはサンプルからのNMR信号を検出する検出回路が設けられている。分光計は、制御部、送信部、受信部等により構成される。送信部は、送信パルスシーケンスに従って送信信号を生成し、その送信信号がプローブに送られる。これにより、サンプルに対して電磁波が照射される。その後、プローブにおいて、サンプルからのNMR信号が検出される。その検出により生じた受信信号が受信部へ送られる。受信部においては、受信信号に対するFFT演算によりNMRスペクトルが生成される。そのNMRスペクトルが、必要に応じて、スペクトル処理装置32へ送られる。スペクトル処理装置32をNMR測定装置30内に組み込むようにしてもよい。
スペクトル処理装置32は、実施形態において、コンピュータ等の情報処理装置により構成される。図2においては、スペクトル処理装置32が有する複数の代表的機能が複数のブロックによって表現されている。
図2において、観測されたスペクトルyがスペクトル処理装置32に入力される。そのスペクトルyは選択手段として機能する選択部34に与えられる。選択部34は、スペクトル中において、ユーザーにより指定された部分を、ベースライン成分の中の代表部分yIとして、選択するものである。一般に、広帯域成分であるベースライン成分は、スペクトルyの全体に及んでいるところ、注目波形成分等のベースライン成分以外の成分がフィッティング対象とならないように、代表部分yIが選択される。なお、スペクトルyの内容次第では、代表部分yIの選択によらずに、スペクトルyの全体が以下に説明する探索部36に入力されてもよい。
探索部36は探索手段として機能するものである。探索部36は、以下の(1−2)式に示す前提条件(第1条件)、及び、以下の(1−1)式に示す基本条件(第2条件)に従って、重み列xの最適解を探索する。重み列xは複数の重みを要素としたベクトルである。
上記(1−2)式は、ベースライン成分をベースラインモデルでフィッティングすることを示しており、より詳しくは、ベースライン成分中の代表部分yIをベースラインモデルAx中の対応部分SIAxでフィッティングすることを示している。Aは基底関数列を示しており、具体的には、複数の基底関数である複数の基底信号ベクトルを列ベクトルとする行列である。Axは、複数の基底信号ベクトルに対して複数の重みを作用させることにより得られるベースラインモデル(計算上のベースライン成分)である。代表部分yIは、観測点Iに対応するスペクトル要素からなるベクトルであり(例えばIは1から256の値をとる)、SIは観測点Iに対応する標本ベクトルである。そのSIをベースラインモデルAxに作用させることにより、対応部分SIAxが抽出される。
上記(1−1)式は、重み列xのLpノルム(但しp≦1)を最小化する条件を示すものである。pが1以下の場合、Lpノルムは、その最小化問題を解く過程で、ノルム計算対象である解のスパース性を高める働きを発揮する。pの範囲は0≦p≦1であるが、pが0の場合に解が不安定になることもあるので、それを回避するためには、望ましくは0<p≦1である。スペクトル処理装置32は、最終的に求まる重み列がスパース性を有することを前提とし、Lpノルムが解のスパース性を高める性質を活用して、重み列の探索を行うものである。一般に、pは1であるが、重み列のスパース性が高いと予測される場合、pを例えば0.75又は0.5とするようにしてもよい。重み列を構成する要素の個数をNとし、代表部分yIを構成する要素の個数をMとした場合、代表部分yIはM行1列の行列、標本行列SIはM行N列の行列となる。基底関数の個数をKとした場合、基底関数列AはN行K列の行列である。
なお、ノルムに関しては、一般に、以下の(2−1)式で示す表現及び以下の(2−2)式で示す表現が認められる。本願において“Lpノルム”は基本的に以下の(2−2)式で定義されるノルムである。nはベクトルを構成する要素数である。
図2において、符号38は基底関数列Aを模式的に示しており、符号40は重み列xを模式的に示している。符号42は基底関数列A及び重み列xから生成されるベースラインモデルAxを示している。符号44は第1条件及び第2条件が満たされるように重み列xの最適解を探索する部分を示している。最適解の探索に際しては、Iterative Soft Thresholding (IST)法、Iterative Reweighted Least Square (IRLS)法、Spectroscopy by Integration of Frequency and Time domain (SIFT)法、等を利用し得る。p=1の場合には、IST法を利用するのが望ましい。Pが1よりも小さい場合には、IRLS法を用いるのが望ましい。以下の第1実施例はIST法に基づくものである。以下の第2実施例及び第3実施例はIRLS法に基づくものである。
減算部46は、重み列xの最適解x’が求まった段階で、スペクトルyから最適ベースラインモデルAx’を減算する減算手段である。これによりベースライン成分が除外又は低減されたスペクトルy’が得られる。それが解析部48における解析対象となる。
図3には、スペクトル処理装置として機能する情報処理装置の構成例が示されている。情報処理装置は、CPU50、メモリ56、入力器52、表示器54等を有する。メモリ56上には、スペクトル処理プログラム58、及び、スペクトル解析プログラム60が格納されている。スペクトル処理プログラム58は、図1及び図2を用いて説明したスペクトル処理を実施するためのプログラムである。スペクトル解析プログラム60は、前処理後のスペクトルを解析するためのプログラムである。それらのプログラム58,60は、CPU50において実行される。
CPU50は、探索手段、選択手段、減算手段等として機能する。入力器52はキーボード、ポインティングデバイス等によって構成され、入力器52を用いて上記代表部分がユーザーにより指定される。表示器54は例えばLCDによって構成され、そこには処理前のスペクトルが表示される。ユーザーにより、そのスペクトル上において代表部分を構成する範囲又は代表点群が指定されてもよい。
次に、第1実施例を説明する。この実施例ではp=1であり、以下の(3)式に示す評価値Jを最小化する条件に従って、最適解としての重み列が探索される。
上記(3)式中の第1項は、上記第1条件に対応するものであり、ベースライン成分の代表部分yIとベースラインモデルの対応部分SIAxとの間の残差(yI−SIAx)のL2ノルムを計算する部分である。上記(3)式中の第2項は、上記第2条件に対応するものであり、重み列xのL1ノルムを計算する部分である。λは正則化重みである。
一般に、評価値Jの最小化問題を解くためには、評価値Jの計算式の微分が必要となるところ、第1項については微分可能であるが、第2項については微分困難である。そこで、公知のIST(Iterative Soft Thresholding)法に従って上記問題を解くことにする。
図4には、第1実施例に係るアルゴリズムが示されている。S10では、スペクトルy及び代表点群を特定する座標列(観測点)Iが取り込まれる。S12では、標本行列SI及び基底関数列Aが構成される。S14に示される数式は、上記の(1)式及び(3)式に含まれる行列SIAを計算の便宜上、Bと置き換えるものである。S16では、スペクトルy及び座標列Iに基づいて、代表部分yIが構成される。S18では、リプテリッツ係数Lfが計算される。それは、Bについての複数の固有値の内の最大値つまり最大固有値である(リプテリッツ係数Lfを最大固有値以上としてもよい)。BHはBのエルミート転置である。S20では、初期重み列x0が代表部分yIから生成される。B+は行列Bの疑似逆行列である。また、S20ではkが初期化される。
S22では、第1終了条件が判定される。すなわち、演算回数を示すkが最大値kmax以下であるか否かが判断される。k≦kmaxであれば、S24及びS26が実行される。S24に示されている計算式及びS26に示されている計算式は、重み列xを二段階で更新するものである。S24に示されている計算式は上記(3)中の第1項の一階微分により求められるものであり、その計算式は、第1項のL2ノルムを最小化するように、重み列xを更新するものである。S26に示されている計算式は上記(3)式の第2項に相当するものである。第2項が微分困難であるために、ソフト閾値関数soft( )が利用されている。その関数は、S24の計算結果xk+1に基づいて、以下のように、当該xk+1を再構成するものである。
上記ソフト閾値関数soft( )の利用により、L1ノルムを最小化するように、重み列xが更新される。S28ではkがインクリメントされる。上記第1終了条件が満たされた時点での重み列xk+1が最適重み列として出力される。
上記アルゴリズム中に第2終了条件を加えてもよい。例えば、S26の直後において、以下の(5)式で定義される指標eを演算し、指標eが一定値以下になった場合にアルゴリズムを終了させてもよい。
上記(5)式の分母は、更新後の重み列xk+1のL2ノルムであり、上記(5)式の分子は、更新前後の重み列xk,xk+1の差分のL2ノルムである。これ以外の第2終了条件を定めてもよい。
図5には、第2実施例に係るアルゴリズムが示されている。このアルゴリズムはIRLS法に従うものである。評価値Jは(6)式に基づいて定義される。
例えば、pは0.75である。もちろん、pとして1以下の他の数値を与えることも可能である。IRLS法に従って、(6)式中における微分困難な第2項は重み行列Wを利用して以下の(7)式のように表現される。
上記の(7)式における第2項はL2ノルムの計算式となっており、それに対しては一階微分が可能である。
図5において、図4に示した工程と同様の工程には同一符号を付し、その説明を省略する。S10aでは、スペクトルy及び座標列Iが取り込まれ、またpの値が定義される。例えば、上記のようにpに0.75が与えられる。S30においては、0での割り算を防止するための係数ε、及び、正則化重みλに数値が与えられる。その後のS12〜S22の工程については既に説明した。
続くS32では、IRLS法に従って、複数の重み要素wiが定義される。S34では複数の重み要素を対角要素とし、対角要素以外をすべて0とした重み行列Wが定義される。S36に示されている計算式が重み列xの更新式である。それは上記(7)式を一階微分することにより求められるものである。S36に示されている計算式を繰り返し実行することによって、評価値Jを最小化する条件を満たす、重み列xの最適解が求められる。
図6乃至図8を用いて第2実施例の効果を説明する。図6は余弦形状でベースライン成分を模擬した場合における実験結果を示すものである。図7は直線形状でベースライン成分を模擬した場合における実験結果を示すものである。図8はブロードなローレンス信号(0.75Hz付近の成分)でベースライン信号を模擬した場合における実験結果を示すものである。
図6乃至図8において、(A)には人工的に加えられたベースライン成分を含むスペクトルが示されている。(B)には従来法(一般的な最小二乗フィッティング法)によるベースライン成分推定結果が示されている。実線がスペクトルを示しており、破線が推定されたベースライン成分を示している。(C)には従来法によるベースライン成分減算結果が示されている。(D)には従来法での重み列(クロス記号を参照)と第2実施例での重み列(丸記号を参照)が示されている。(E)には第2実施例によるベースライン成分推定結果が示されている。実線がスペクトルを示しており、破線が推定されたベースライン成分を示している。(E)には第2実施例によるベースライン成分減算結果が示されている。
なお、従来法及び第2実施例のいずれにおいても同じ基底関数列が利用されている。具体的には、基底関数列は、1〜10次のN次多項式、1から10次の余弦関数、及び、1から10次の正弦関数からなるものである。観測点の設定は手動で行っている。図6乃至図8において、観察の都合上、ベースライン成分の表示に際してはオフセットが加えられている。
図6乃至図8において、従来法では、(D)に示されるように、基底関数列の全体に対して概してまんべんなく重みが与えられており、その結果、(B)に示されている過剰フィッティングが生じている。特に観測点を設定していない区間に対する過剰なフィッティングが生じている。一方、第2実施例によれば、(D)に示されるように、重み列がスパース化されており、つまり、少数の基底関数によってベースライン成分が推定されている。この結果、(E)に示されるように、ベースライン成分が高い精度で推定されている。第1実施例及び以下に説明する第3実施例においても基本的に同様の良好な結果を得られる。
図9には、第3実施例に係るアルゴリズムが示されている。このアルゴリズムは、劣決定型のIRLS法に基づくものである。その方法は、データ総数Nよりも標本集合の要素数Mが小さい場合に利用可能なものである。例えば、pは0.75である。もちろん、pとして1以下の他の数値を与えることも可能である。
上記(1−1)式及び(1−2)式は、IRLS法に基づき、重み行列Wを用いて、以下の(8−1)式及び(8−2)式のように書き換えられる。(8−2)式は(1−2)式と同じものである。
上記の行列SIAが横長行列(つまり劣決定型)の場合、上記(8−1)式及び(8−2)式は、劣決定型のIRLS法に基づき、以下の等価式(9−1)式及び(9−2)式に書き換えることが可能である。
(9−1)式及び(9−2)式を満たす解は以下の(10)式のように計算される。
図9に示すアルゴリズムは以上のような考え方に基づくものである。なお、図9において、図4及び図5に示した工程と同様の工程には同一符号を付し、その説明を省略する。
S10aにおいてスペクトルy及び座標列Iが入力され、続いてS30aではεに対して所定の値が与える。S12〜S22、S32及びS34については既に説明したとおりである。S38に示されている計算式は重み列xの更新式である。それは上記の(10)式と同じである。
実施形態によれば、多くの基底関数を用意しておけるので、特殊なベースライン成分に対しても良好なフィッティング結果を得ることが可能となる。過剰フィッティングの問題も回避することが可能である。また、従来においては、ユーザーにおける関数選択や次数指定の適切さがフィッティング結果を左右していたが、上記構成によれば、ユーザーに格別な負担をかけることなく安定して良好なフィッティング結果を得られる。上記構成は従来の関数フィッティング法を発展させた手法と位置付けられ、その意味においてフィッティングのメカニズムをユーザーにおいて直感的に理解しやすいという面がある。上記実施形態ではNMRスペクトルが処理されていたが、X線スペクトル、マススペクトル等の他のスペクトルが処理されてもよい。