JP6160366B2

JP6160366B2 - Ｎｍｒ信号処理システム

Info

Publication number: JP6160366B2
Application number: JP2013174931A
Authority: JP
Inventors: 孝子末松; 博明内海; 朋喜中尾; 利博古川
Original assignee: Jeol Ltd
Current assignee: Jeol Ltd
Priority date: 2013-08-26
Filing date: 2013-08-26
Publication date: 2017-07-12
Anticipated expiration: 2033-08-26
Also published as: US20150057979A1; US10451695B2; JP2015042964A

Description

本発明はＮＭＲ（核磁気共鳴）信号処理システムに関し、特に、数学的なモデルを用いてＮＭＲ信号に含まれる複数の信号成分を推定するシステムに関する。

核磁気共鳴分光法は、一般に、化合物中の観測核が生じさせるＮＭＲ信号としてのＦＩＤ（Free Induction Decay）信号を検出し、その検出されたＦＩＤ信号をフーリエ変換し、ＮＭＲスペクトルを生成するものである。そのスペクトルに基づいて化合物の分子構造等が解析される。近時、核磁気共鳴分光法を利用して定量分析を行うシステムが実用化されており、それはｑＮＭＲ（quantitative NMR）とも言われている。ｑＮＭＲにおいては、スペクトル中において解析対象となる注目ピークが有する強度（信号強度）から、サンプルに含まれる特定物質の含量や濃度が求められる。よって、ｑＮＭＲにおいては強度を精度良く求めることが非常に大切である。従来、強度を求める方法として以下の方法が知られている。なお、一般に、qNMRでは、定性的な性質がある程度分かっているものを解析対象とするので、解析に使用するNMRスペクトルにおけるピーク数、周波数などの情報は既知である。

第１の方法は積分法と言われている方法である。この方法では、スペクトル上において注目ピークを包含する周波数範囲が設定され、その周波数範囲内で積分処理を行うことにより、当該注目ピークについての強度が求められる。その強度は注目ピーク波形が有する面積に相当する。この方法では、どこまでを積分範囲として決定するのかの判断が難しく、積分範囲如何によって演算される強度が変化してしまうという問題がある。また、複数のピークが部分的に重なり合っているような場合に正確な強度を求めることが困難であるという問題がある。

第２の方法は波形分離法（特にカーブフィッティング法）である。この方法は、スペクトル上において注目ピークを波形で近似（模擬）するに当たり、その波形を定義するパラメータ（例えば強度、線幅、周波数）の最適値を探索するものである。しかしながら、この方法では、スペクトル形態次第では、演算結果がパラメータの初期値に依存するという傾向が認められ、つまり、それらの初期値が異なると、演算結果も異なってしまうという問題がある。

第３の方法は、実時間軸上のＦＩＤ信号を表現する数学的モデルを導入し、その数学的モデルに含まれるパラメータ（パラメータ値としての数値又は数値列）を推定する方法である。この方法の１つとしてProny法が知られている（非特許文献１）。このProny法では、ＦＩＤ信号が複数の信号成分で構成されていることが仮定され、個々の信号成分を定義するパラメータの内容が推定される。個々の信号成分はスペクトル上の個々のピークに相当するものである。もっとも、実際には、スペクトル上には様々な成分が含まれ、またスペクトル波形はＦＩＤ信号に対して行われた各処理の影響を受ける。そこで、スペクトル上の注目ピーク数よりも多い信号成分を想定してProny法を適用することが提案されている（高次推定法）（非特許文献２）。この手法によれば、上記第１の手法について指摘した積分範囲の設定上の問題を回避できる。また、上記第２の手法について指摘したような初期値依存性の問題も回避できる。よって、この第３の手法は、上記第１及び第２の手法に比べて、より客観的な手法であると言いうる。

F. Malz, H. Jancke: J. Pharm. Biomed. Anal., 38, 813(2005). Santosh Kumar Bharti, Raja Roy: TrAC, 35, 5(2012). William H. Press, Saul A. Teukolsky, Willian Vetterlin, Bruan P. Flannery: "Numerical Recipes in C [日本語版]", p.505 (2009), （技術評論社）. P.Koehl: Prog. nucl. magn. reson. spectrosc, 34, 257(1998). Amab K. Shaw, Krishna Naishadham: "ARMA-Based Time- Signature Estimator for Analyzing Resonant Structures by the FDTD Method", IEEE Antennas and Propagation Society, 49, 327(2001). Parks, T.W., C.S. Burrus.: Digital Filter Design (Topics and Digital Signal Processing), p.226(1987), (Wiley-Interscience). S. Haykin, 鈴木博: "適応フィルタ理論", p.912 (2001), (科学技術出版). Y. T. Chan: "Spectrum Estimation via the High-Order Yule- Walker Equations" IEEE Transactions on Acoustics, 30, 682(1982). Steven M. Kay: "Noise Compensation for Autoregeressive Spectral Estimates" IEEE Transactions on Acoustics, 28, 292(1980). S. Braun, Y. M. Ram; "Determination of structural modes via the Prony model: System order and noise induced poles" J. Acoust. Soc. Arm., 81, 1447(1987) 中溝高好;"信号解析とシステム同定"p.198(1988), (コロナ社). 村田恵介;"AR過程に基づくピーク推定によるＦＩＤの変換：周波数と減衰率の同時推定"日本核磁気共鳴学会,第47回ＮＭＲ討論会(2008). 村田恵介,久保田一;"ＮＭＲにおける極推定"信学技報p.35(2007-08),電子情報通信学会.

しかしながら、上記の第３の手法では、想定する信号成分の個数（次数）に推定結果が依存するという面が認められ、想定する信号成分の個数を正しく定めることができない場合には推定精度が劣化するという問題を指摘できる。

本発明の目的は、ＮＭＲスペクトル上の注目ピークについて強度を精度良く求めることにある。あるいは、本発明の目的は、ＮＭＲスペクトルの元になるＦＩＤ信号を表す数学的モデルに対して適切な推定次数を与えて、注目ピークに対応するパラメータを精度良く求めることにある。

本発明に係るＮＭＲ信号処理システムは、ＮＭＲ測定で得られる実信号をｑ個の信号成分の和として表現する数学的モデルに基づき、当該数学的モデルにおける次数ｑを可変しながら、次数ｑごとに、前記実信号に基づいて前記ｑ個の信号成分を定義する複数のパラメータを推定する推定手段と、前記推定手段が推定した次数ｑごとの推定結果を評価する評価手段と、前記推定手段が評価した次数ｑごとの評価結果に基づいて最適次数を決定することにより、前記推定手段の推定結果の中から前記最適次数に対応する複数のパラメータを特定する決定手段と、を含むことを特徴とするものである。

上記構成によれば、数学的モデルにおける次数ｑ（推定次数、信号成分数）を試行的に可変させてみて、次数ｑごとに実信号に基づく推定結果を得て、それらを評価することにより、最適次数（最適推定次数）を決定することが可能である。よって、推定次数に依存した推定誤差増大という問題を解消することができる。換言すれば、実信号に応じて実際に設定、利用する推定次数を最適化できる。この点に着目して、本手法を推定次数最適化法と称することが可能である。

上記構成においては、最適次数の探索過程において、次数ｑごとに複数のパラメータが推定されるから、最適次数が決定された時点で再度の推定演算は基本的に不要であり、既に演算されている推定結果の中から、最適次数に対応する複数のパラメータを特定すればよい。もっとも、再度の推定演算を行うことも可能である。数式モデルにおいて、次数ｑは一般にＮＭＲ信号をｑ個の信号成分で模擬することを意味する。数学的モデルにおいて、個々の信号成分はパラメータによって定義される。パラメータは１又は複数の変数により構成されるものであり、後述する実施形態において、パラメータは極及び複素強度の組である。通常、ｑは1,2,3,…である。

上記の実信号はＮＭＲ信号つまりＦＩＤ信号である。もっとも、推定演算に先立ってＦＩＤ信号に対して前処理を施すようにしてもよい。数学的モデルは、望ましくは、時間軸上のＦＩＤ信号を複数の指数減衰正弦波の和として表現する数式である。次数ｑの可変範囲はユーザー指定されあるいは自動的に設定される。最適次数を効率的に探知できるように、次数ｑの可変順序を定めるようにしてもよい。その場合、推定と評価とを同時進行で行うことも可能である。基本的に次数ｑごとにｑ個のパラメータが推定される。但し、例外的に一部のパラメータについての推定を省略してもよい。推定結果の評価に際しては望ましくは実信号が参照される。

望ましくは、周波数軸上におけるｐ個の注目ピークに応じて定められた選抜条件に従って、次数ｑごとに、前記推定手段が推定したｑ個のパラメータの中から評価対象とする複数のパラメータを選抜する選抜手段を含み、前記評価手段は前記選抜手段の選抜結果に基づいて次数ｑごとに評価値を演算し、前記決定手段は前記評価手段が演算した次数ｑごとの評価値の中から最良評価値を特定することにより前記最適次数を決定する。

上記構成によれば、次数ｑごとに推定されたｑ個のパラメータ全部を評価対象とするのではなく、その中で選抜条件を満たす複数のパラメータだけを評価対象にすることが可能である。つまり、ｐ個のピークに適合する可能性の高いものだけを評価対象にすることが可能である。これにより評価精度を高められる。この点に着目して、本手法を選抜評価法と称することが可能である。ここで、ｐは1以上の整数であり、ｑの可変範囲は望ましくはｐとｑ_maxとの間である。なお、次数ｑの最大数をｑ_maxとする。

望ましくは、決定手段は、最適次数が決定された場合に、その最適次数に相当する個数のパラメータ全部を特定するのではなく、その内でｐ個のピークに適合する（つまり選抜された）ｐ個のパラメータだけを特定する。これによれば、最適次数の下で推定されたパラメータ群の中で、例えばノイズ波形に相当するパラメータを棄却して、真のピークに相応するパラメータだけを取り出すことが可能である。よって、目的に合致した推定を実現でき、推定精度を高められる。この点に着目して本手法を選抜推定法と称することが可能である。

望ましくは、前記選抜手段は、次数ｑごとに、前記推定手段が推定したｑ個のパラメータの中から評価対象とするｐ個のパラメータを選抜する機能を有し、前記評価手段は、次数ｑごとに、前記選抜手段が選抜したｐ個のパラメータに基づいて評価値を演算する機能を有する。すなわち、基本的に、高次推定法の下で、ｐ以上の次数ｑごとにｑ個のパラメータが推定され、その内で、ｐ個のパラメータが選抜される。選抜の結果、例外的にｐ未満の個数のパラメータが選抜されてもよい。その場合に、それを評価対象に含めてもよいし、評価対象から除外してもよい。

望ましくは、前記選抜条件には、前記ｐ個の注目ピークを包含する選抜領域を利用して次数ｑごとに選抜を行うための第１選抜条件が含まれ、前記選抜手段は、前記推定手段が推定したｑ個のパラメータの中から前記第１選抜条件を満たす複数のパラメータを選抜する第１選抜手段を含む。この構成によれば、周波数の観点から個々のパラメータの良否を判断して、注目ピークに該当する可能性の高い優良なパラメータだけを残すことが可能である。

望ましくは、前記選抜条件には、真のピークである可能性が高い順で上位ｐ個のパラメータを選抜するための第２選抜条件が含まれ、前記選抜手段は、前記第１選抜手段が選抜した前記第１選抜条件を満たす複数のパラメータの中から前記第２選抜条件に従って上位ｐ個のパラメータを選抜する第２選抜手段を含む。この構成によれば、第１選抜処理によって残された複数のパラメータの中から、相対評価によって、注目ピークの可能性がより高いものを必要数だけ厳選することが可能である。第２選抜条件は望ましくは周波数とは別の物理量に関する条件であるのが望ましい。例えば、線幅、強度、形状等に関する条件であってもよい。

望ましくは、前記次数ｑの可変範囲の下限はｐであり、前記次数ｑの可変範囲の上限はｐに基づいて設定される。次数ｑの可変範囲内に最適次数が含まれるようにその可変範囲を定めるのが望ましい。最適次数は注目ピーク数ｐ以上であると言えるから、可変範囲の下限としてｐが設定される。可変範囲の上限は例えばｐに対して係数ｍを乗算することにより定めるようにしてもよい。

望ましくは、前記評価手段は、前記推定された複数のパラメータを前記数学的モデルに与えることにより生成される合成信号と、前記実信号と、を比較することにより、評価値を演算する。この構成によれば、推定信号の波形と、実信号の波形とが比較されて、一致度又は相違度をもって評価値が算出される。時間軸上の評価に代えて周波数軸上の評価（例えばスペクトル波形比較）を行うことも可能である。後者において、周波数軸上に評価対象となる周波数範囲を定めてもよい。

望ましくは、前記実信号に対して前処理としての帯域通過処理を施すフィルタ手段を含み、前記推定手段は前記帯域通過処理後の実信号に基づいて前記複数のパラメータの推定を行う。この構成によれば信号処理で不要となる周波数成分を事前に除外又は抑圧して、演算量を軽減でき、あるいは推定精度を高められる。

本発明に係るプログラムは、ＮＭＲ測定で得られる実信号を入力し、当該実信号をｑ個の信号成分の和として表現する数学的モデルに基づき、当該数学的モデルにおける次数ｑを可変しながら、次数ｑごとに、前記実信号に基づいて前記ｑ個の信号成分を定義する複数のパラメータを推定し、それらの推定結果を第１記憶部に格納する機能と、前記記憶部に書き込まれた次数ｑごとの推定結果を評価し、次数ｑごとに評価結果を第２記憶部に格納する機能と、前記第２記憶部に格納された次数ｑごとの評価結果に基づいて最適次数を決定することにより、前記第１記憶部に格納された推定結果の中から前記最適次数に対応する複数のパラメータを特定する機能と、を含むことを特徴とするものである。このプログラムは、望ましくはＣＤ−ＲＯＭ等の記憶媒体に格納され、記憶媒体を介してＮＭＲ信号処理システムにインストールされる。もちろん、ネットワークを経由してインストールされてもよい。

本発明によれば、ＮＭＲスペクトル上の注目ピークについて強度を精度良く求めることができる。あるいは、本発明によれば、ＮＭＲスペクトルの元になるＦＩＤ信号を表す数学的モデルに対して適切な推定次数を与えて、注目ピークに対応するパラメータを精度良く求めることが可能である。

従来のProny法による強度推定方法を示す概念図である。従来の高次Prony法の問題点を説明するための図である。本発明に係るパラメータ推定方法の好適な実施形態を示すフローチャートである。図３に示したパラメータ推定方法の実行に伴って生じる推定結果、選抜結果及び評価結果を示す概念図である。推定次数の可変に伴う評価値の変化及び最適推定次数の頻度分布を示す図である。模擬ＦＩＤ信号に基づくＮＭＲスペクトルを示す図である。実測ＦＩＤ信号に基づくＮＭＲスペクトルを示す図である。本発明に係るＮＭＲ信号処理システムの好適な実施形態を示すブロック図である。選抜領域の第１設定例を示す図である。選抜領域の第２設定例を示す図である。評価範囲を示す図である。

以下、本発明の好適な実施形態を説明する。最初にＦＩＤ信号モデルに関する基礎的な説明を行い、続いて本実施形態に係る方法及びシステムを具体的に説明する。

（１）ＦＩＤ信号モデル
ＮＭＲ測定で得られた時間軸上のＦＩＤ信号は、以下の式（1）に示すように、複数の指数減衰正弦波（複数の信号成分）の線形結合として表され得る。

但し、Ｔｓはサンプリング周期であり、ｈ（ｎ）はＦＩＤ信号における時刻ｔ（=ｎＴｓ）での標本値であり（但しｎ＝0,1,2,3,…）、ｐはＦＩＤ信号を構成する信号成分の個数を示しており、以下においてｐを「構成次数」と呼ぶ。また、

はそれぞれＦＩＤ信号の強度、位相、線幅（半値半幅）、共鳴周波数である。それらはいずれも信号成分ごとに特定されるものである。添え字である

は個々の信号成分を特定する変数であり、それは1,2,…,pの値をとる。ここで、

をそれぞれＦＩＤ信号についての「複素強度」及び「極」と呼ぶ。「極」も複素数である。信号成分ごとに複素強度及び極の組（後において「パラメータ」と称する。）が存在する。最終的に求めたい値は、式（1）中の強度

であり、つまり信号成分ごとの強度である。これはＮＭＲスペクトルにおけるピークごとの強度に相当する。
さて、計算上の都合から、上記の式（1）に対して、離散関数を対象とした数学的操作であるｚ変換（ｚtransform）を適用すると、以下の式（2）が得られる（非特許文献５）。

式（2）中のz‐¹は遅延演算子であり、

はそれぞれAR係数列、MA係数列である。式（2）における右辺の表現はAuto Regressive Moving Average (ARMA)モデルと呼ばれている。なお、

は（2）式の関数Ｈ（ｚ）における極に相当する。

式（1）及び（2）に基づき、ＦＩＤ信号を構成する信号成分ごとの極及び複素強度を推定した上で、信号成分ごとの強度を求める手法として、以下に説明するProny法が知られている（非特許文献６）。

（２）Prony法
ＦＩＤ信号を表す（1）式のｈ（ｎ）のｚ変換を以下の式（3）のように表現することもできる。

但し、Ｎはデータポイント数である。以下において、ｈ_ｎ=ｈ（ｎ）とする。
式（2）の右辺と式（3）の右辺から、以下の式（4）

が成り立つ。式（4）の左辺を展開し、両辺間でｚ‐ⁿの係数を比較すると、

が導かれる。但し、（5）式中の

はそれぞれ以下に示すとおりである。

但し、Ｔは転置を意味しており、Ｈはエルミート転置を意味している。式（5）に基づいて、式（2）の分母多項式の係数列

を求めたいのであるが、式（6）〜（8）からわかるように、式（5）から係数列ａを求めるために必要となる、式（2）の分子多項式中の係数列

は未知である。そこで、

を満足するａが、式（5）を満足するａの必要条件であることに基づき、式（9）を解くことにする。但し、Ｈａ、ｈａは以下のとおりである。

Ｈａの列正則性が成立するものと仮定する。この仮定は、必ずしも一般的に成立するものではないが、後に説明する実験結果から、その仮定の成立が肯定され得る。非特許文献７においても同様の仮定をおいて議論を行っている。上記仮定によれば、式（9）は次の式（12）

と等価である。よって、係数列ａは

で与えられる。上記の式（13）によって得られた係数列ａから、式（2）中の分母多項式を特定した上で、その分母多項式の因数分解を行うことにより、信号成分ごとの根

が求められる。これから、式（2）に基づき、信号成分ごとの極が以下の式（14）のように求められる。

一方、上記の式（14）で得られた、信号成分ごとの極

を式（1）に代入した上で、n=0、1、…、N‐1に対応する連立多項式を解くことにより、ＦＩＤ信号を構成する信号成分ごとの複素強度

を推定できる。その複素強度計算について行列表記を行うと、以下の式（15）ように表される。

但し、Ｅ、ｃ、ｈはそれぞれ以下のとおりである。

また、極

は、多重極を含まないものとし、つまり

が成立するものとする。よって、式（16）のＥは列正則であり、式（15）のｃは

で与えられる。強度

は複素強度

の絶対値である。すなわち、以下のとおり信号成分ごとの強度が推定される。

図１には、以上説明したProny法に基づく強度推定方法が概念図として示されている。Ｓ１０では式（13）が計算され、その計算結果としてＳ１４に示す係数列ａが特定される。式（13）中のＨａ及びｈａは、Ｓ１２に示すように、式（10），（11）によって特定される。すなわち、それらの行列はＦＩＤ信号としての時間軸上の振幅値列により構成し得るものである。後述する高次Prony法に基づく強度推定方法では、次数が異なるものの、基本的に、図１に示した手法と同じ手法が適用される。その場合において、式（10），（11）で表される各行列には実測ＦＩＤ信号が与えられる。

Ｓ１６では、係数列ａが式（2）中の右辺の分母多項式に与えられ、その上で、Ｓ１８において、その分母多項式の因数分解を行えば、関数Ｈ（ｚ）について、ｐ個の根を特定できる。Ｓ２０では、式（14）に基づいてｐ個の根からｐ個の極が特定される。Ｓ２４では、式（19）が計算される。そこに含まれるＥ及びｈは、Ｓ２２で示すように、式（16）,（18）で特定される。ここで、ＥはＳ２０で算出された極によって特定され、ｈはＦＩＤ信号によって定義される。式（19）を計算することにより、Ｓ２６で示すように、信号成分ごとの複素強度ｃが特定される。そして、信号成分ごとの複素強度ｃの絶対値として、信号成分ごとの強度Ｃが求まる。各強度ＣはＮＭＲスペクトルにおける各ピークの強度である。それらの強度の推定をもって各ピークを定量化することが可能である。

（３）高次Prony法
実際には、雑音や測定誤差の影響を受けて、実測ＦＩＤ信号は、式（1）を修正した以下の式（20）のように表される。

ここでｖ（ｎ）は加法性白色雑音である。式（20）によって表されている

は、解析対象つまり実測ＦＩＤ信号であり、これは理想的なあるいは純粋な数学的表現形式としての

と区別される。上記の式（20）を前提とした場合、式(1)、（2)と実測ＦＩＤ信号とが不一致となるために、精確な強度

を推定することはできない。
そこで、高次Prony法が提案されている（非特許文献１１）。この手法は、構成次数ｐよりも多い信号成分数ｑをもって（つまりｑ＞ｐの条件で）ＦＩＤ信号モデルを定義するものである。具体的には実測ＦＩＤ信号が以下のように表現される。

この式は、Prony法における式（2）に対応するものである。上記の式（21）によって表される

を「推定モデル」と呼び、ｑを「推定次数」と呼ぶ。なお、構成次数ｐはＮＭＲスペクトル上の注目ピーク数に対応し、推定次数はＮＭＲスペクトル上のノイズを含む諸ピーク数に対応する。上記の式（21）は、実時間上のＦＩＤ信号を

と表現することを前提とするものである。この式（22）は上記の式（1）に対応する。

この高次Prony法は、式（21）、（22）に基づいて、ＦＩＤ信号の極と強度を推定するものである。次数を高めた点を除き、高次Prony法に基づく推定方法はProny法における推定方法と同じである。すなわち、基本的に、図１に示したプロセスによってＦＩＤ信号を構成する信号成分ごとに極と複素強度が推定される。

（４）高次Prony法の問題点
図２は、構成次数ｐを４とした模擬ＦＩＤ信号に対して高次Prony法を適用した結果を示している。上段の（A）には、模擬ＦＩＤ信号のスペクトルが示されている。そのスペクトルには４つの注目ピークが含まれている。中段には推定結果としての各ピークの極が示されており（縦軸は線幅（半値半幅））、下段には推定結果としての各ピークの強度が示されている。中段及び下段において、左側に示す(B)は推定次数ｑが４（＝ｐ）の場合の推定結果を示しており、中央に示す(C)は推定次数ｑが１６の場合の推定結果を示しており、右側に示す(D)は推定次数ｑが３２の場合の推定結果を示している。いずれも横軸は周波数を示している。記号“○”は真値を表しており、記号“＋”は推定値を示している。なお、真値が存在する範囲から大きく外れる一部の推定値については図示省略されている。

左側の(B)に示されているように、推定次数ｑが構成次数ｐと同一の場合、観測信号中の雑音成分の影響により、満足できる推定精度が得られていない。これに対して、中央の(C)に示されているように、ｑが１６、つまり推定次数ｑを構成次数ｐの４倍にした場合、４つの真値の上又は近傍に４つの推定値が存在しており、その部分だけを着目するならば、ｑが４の場合の推定結果に比べて、推定精度がかなり向上している。一方、右側の(D)に示されているように、ｑが３２、つまり推定次数ｑを構成次数ｐの８倍にした場合、逆に、推定精度が低下してしまっている。つまり、ｑが４の場合に比べれば推定精度が良いものの、ｑが１６の場合に比べると推定精度が劣化している。ＦＩＤ信号を過剰な信号成分数で近似しようとしたために、逆に推定精度が低下したものと推認される。本例では、ｑが１６の場合において最良の推定結果が得られている。これと同様の認識が非特許文献１０、１１においても示されている。特に非特許文献１１ではパラメータの推定精度が推定次数ｑの大きさに依存することが報告されている。

以上のように、高次Prony 法によれば、推定次数ｑを的確に設定できる限りにおいて非常に精度良く極や強度を推定できる。しかし、推定次数ｑを的確に設定できない場合、推定精度が大きく低下してしまうおそれがある。解析対象となった実測ＦＩＤ信号の態様如何によって最適な推定次数ｑが変わり得ると考えられるから、事前に推定次数ｑを一義的に特定しておくことは困難である。そこで、真値に最も近い結果を与える推定次数、つまり最適な推定次数、を探索して特定する必要がある。

（５）本実施形態に係る方法
本実施形態に係る方法（本手法）は、高次Prony法に対して選抜再構成評価法(SRE; Select Reconstitution Evaluation method)を適用したものである。本手法に含まれる幾つかの特徴事項に着目し、本手法を推定次数最適化法、選抜評価法及び選抜推定法と称することも可能である。

以下、図３のフローチャートに基づいて、本手法を具体的に説明する。図３における左側にはメインルーチンが示されており、同図の右側にはサブルーチンが示されている。なお、以下においては、便宜上、推定モデル（式（21）、（22））を構成する信号成分ごとの極と複素強度の組を

と表記し、また、その組を「パラメータ」と呼ぶ。推定次数ｑを初期値から最大値q_maxまで変化させていく過程において、個々の推定次数ｑごとにｑ個のパラメータが推定される。

図３において、Ｓ３０は、初期条件の入力工程であり、すなわち、構成次数ｐ及び選抜領域Ｒ（ω1,ω2）が設定される。構成次数ｐは、ＮＭＲスペクトル上の注目ピーク群を構成するピーク数に相当する。注目ピーク群は例えば注目原子団（例えば注目塩基）に対応する例えば互いに密集した複数のピークである。ただし、１つの注目ピークだけが対象となってもよい。選抜領域Ｒは、注目ピーク群を含む周波数領域である。ｑＮＭＲではω1及びω2は選抜領域Ｒの下限と上限に相当する。注目ピーク又は選抜領域Ｒについては後述する。qNMRでは、定性的な性質がある程度分かっているものを解析対象とするので、解析に使用するNMRスペクトルにおけるピーク数、周波数などの情報は既知であることから、ｐ、ω1及びω2を、定量分析の対象となる化合物の分子構造から与えることが可能である。

Ｓ３２は初期条件の設定工程であり、すなわち、推定次数ｑについての可変範囲につき、初期値及び最大値ｑ_maxが指定される。図示の例では、初期値ｑ_startとして構成次数ｐが指定されている。注目ピーク数以上の信号成分数が想定されるからである。最大値ｑ_maxとして例えばｐの１０倍である１０ｐが指定される。倍率ｍがユーザーによって指定されてもよいし、又は、倍率ｍが諸条件に基づき自動的に指定されてもよい。いずれにしても、探索対象となる最適推定次数が確実に含まれるように可変範囲を定めるのが望ましい。

Ｓ３４では、個々の推定次数ｑにおいて、上述した高次Prony法に基づき、ｑ個のパラメータからなるパラメータ列が推定される。なお、高次Prony法でのパラメータ推定プロセスは、次数の点を除き、Prony法に基づくパラメータ推定プロセスと基本的に同じである（図１参照）。図３において、推定されたｑ個のパラメータからなるパラメータ列が以下のように表現されている。なお、ｋは、推定次数ｑごとに、１からｑまでの間の整数値をとる。

Ｓ３６は評価工程である。この工程では、基本的に、推定されたｑ個のパラメータの中から、ｐ個の注目ピーク（真のピーク）に対応する可能性の高いｐ個のパラメータ

が選抜される。それらを選抜パラメータ列と呼ぶ。続いて、選抜パラメータ列が評価されて、評価値Ｊが算出される。但し、選抜されたパラメータ数ｐ´がｐ未満の場合、例外的にｐ´個のパラメータが評価対象となる。もっとも、そのような場合には評価値Ｊが悪くなる可能性が高いので、評価値演算を省略してもよい。評価工程の具体的内容についてはサブルーチンに示されており、それについては後に詳述する。

Ｓ３８では、ｑがｑ_max未満であるか否かが判断され、ＹＥＳであれば、Ｓ３９においてｑに１を加算したものがあらたにｑと定義されて、Ｓ３４以降の工程が繰り返し実行される。Ｓ３８の判断結果がＮＯであればＳ４０が実行される。

Ｓ４０では、初期値ｑ_startから最大値ｑ_maxまでの範囲内における推定次数ｑごとの評価値Ｊが相互に比較され、最良評価値（後述する一致度）が判定される。その最良評価値に対応する推定次数が最適推定次数である。すなわち、本手法では、推定次数の試行的な可変により最適推定次数を探索することが可能である。このように最適推定次数が特定されると、それに対応する選抜パラメータ列が特定される。その選抜パラメータ列に基づいて、ｐ個の注目ピークに対応するｐ個の強度が推定される。なお、全パラメータの推定が完了した後に選抜処理を行うのではなく、個々のパラメータが推定された時点で、その都度、選抜処理（例えば後述する選抜領域Ｒに基づく選抜処理）が実行されるようにしてもよい。

次に、上記の評価工程Ｓ３６の具体例について図３の右側に示されるサブルーチンに基づいて説明する。

Ｓ４２では、ｑ個のパラメータからなるパラメータ列の中から、選抜領域Ｒ内に属するｐ´個のパラメータが選抜される（但しｐ´≦ｑ）。その場合、各パラメータを構成する極に含まれる共鳴周波数ωが参照される。この第１選抜によって、注目している周波数領域以外に存在するピークに相当する信号成分が評価対象から外されるので、逆に言えば、ｐ個の真のピークに近い共鳴周波数を有する１又は複数のパラメータだけが評価対象に選ばれるので、評価の精度を高められる。Ｓ４４では、ｐ´がｐより大きいか否か（第２選抜で除外可能な余剰パラメータが含まれるか否か）が判断され、ＹＥＳであれば、第２選抜を行うためにＳ４６が実行される。

Ｓ４６では、選抜領域Ｒ内に属すると判断されたｐ´個のパラメータの中から、例えば線幅（半値半幅）が小さい順で上位ｐ個のパラメータが選抜される。つまり、選抜領域Ｒ内に属すると判断されたものの中から、より真のピークである可能性の高いものがｐ個選ばれる。この第２選抜により、相互比較をもって真のピークにより近いパラメータを選び出すことが可能となる。この場合、第１選抜処理で考慮した周波数以外の情報に基づいて選抜を行えば、選抜効果をより高められる。強度等が選抜基準として参照されてもよい。また、線幅、周波数、強度といった複数の情報が選抜基準として参照されてもよい。なお、Ｓ４４において、ｐ´がｐ以下であれば第２選抜を行う余地がないので、そのままＳ４８が実行される。Ｓ４８を経由せずにサブルーチン処理を終了させてもよい。

Ｓ４８においては、以上のように選抜された、原則ｐ個（例外的にｐ´個）のパラメータからなる選抜パラメータ列

が評価される。具体的には、まず、選抜パラメータ列を用いて以下の評価信号

が構成される。但し、式（23）は例示である。式（23）はｐ個のパラメータが選抜された場合を前提としている。

評価値Ｊは、以下の式（24）に基づいて計算される。

式（24）において、

は、時間軸上の実測ＦＩＤ信号であり、それは式（20）のように表現されるものである。式（24）によって、サンプルポイントごとの振幅差の絶対値についての総和（誤差エネルギーに相当）として、評価値Ｊが求められる。評価量Jは実測FID信号と評価信号の一致度を評価する量である。ある推定次数qにおいて、FID信号の真の極及び強度を推定でき、かつこれらを前述の選抜基準で選び出すことができれば、Jは雑音成分のエネルギーとなる。なお、式（24）に対して規格化のための分母項を追加してもよいが、推定次数ｑにかかわらず実測ＦＩＤ信号は不変であるので、相互比較を行う上ではそのような規格化は必ずしも必要ではない。評価値Ｊつまり誤差エネルギーを最小とする推定次数ｑをもって最適推定次数と言い得る。

図４には本手法の実行に伴って段階的に生成される成果物が模式的に示されている。なお、上記で説明したように、本例では構成次数ｐとして４を想定している。

（Ａ）には、上記Ｓ３４の推定工程の実行結果つまりパラメータ推定結果が示されている。横軸である推定次数ｑは、本実施形態において、ｐからｑ_maxまで１つ刻みで変化している。演算量を削減するために、刻みをより大きくすることも可能であるし、粗い刻みと細かい刻みとを併用してもよい。個々の推定次数ｑごとに、ｑ個のパラメータ１０Ａが推定される。それらはパラメータ列１０を構成する。よって、推定次数が順次上がると、パラメータ推定数も順次増加する。なお、パラメータ１０Ａは上述したように極及び複素強度の組である。

（Ｂ）には、上記Ｓ４２及びＳ４６での選抜工程の実行結果つまり選抜結果が示されている。本実施形態では、基本的に、２段階の選抜処理が適用されている。それらは、選抜領域Ｒを用いた第１選抜処理、及び、線幅が小さい順において上位ｐ個を選抜する第２選抜処理、により構成される。符号１２は、推定次数ｑごとに推定されたパラメータ列１０に対する第１選抜処理を示している。第１選抜処理の結果、5個以上のパラメータが選抜された場合には更に第2選抜処理が適用される。（Ｂ）において、図示の例では、ｑがｐの場合、及び、ｑがｐ＋１の場合において、ｐ＝４未満のｐ´個のパラメータが選抜されている。それらの場合には第２選抜処理は適用されない。一方、ｐ＋２以上の各推定次数ｑにおいて、第２選抜処理が適用され、その適用後に、いずれも４つのパラメータ（選抜パラメータ列）が選抜されている。すなわち、真のピークである可能性の高い順で筆頭から４番目までのパラメータが選抜されている。それ以外は棄却される。第2選抜処理においてパラメータ列の並び換えまでは基本的に不要である。

（Ｃ）には、上記Ｓ４８での評価値演算工程の実行結果つまり評価結果が示されている。符号１４で示すように、推定次数ｑごとに選抜パラメータ列に基づいて評価値Ｊが計算されている。その上で、複数の評価値の中から、最良評価値１６が特定される。これは最適推定次数の決定に相当する。続いて、符号１８で示すように、最良評価値１６に対応する選抜パラメータ列２０が特定される。そして、符号２２で示すように、その特定された選抜パラメータ列が出力され、あるいは、それに基づいて計算される構成次数ｐ分の強度が出力される。なお、ｐ´個のパラメータ又は強度が特定された場合、その事実を報知した上で、それらを出力してもよいし、エラーを報知してもよい。後者の場合、上記Ｓ４４の判断がＮＯであればエラーを記録しておき、評価値の演算を省略してもよい。

以上のように本手法は、第１に、推定次数ｑを事前に一義的に定めるものではなく、それを試行的に変化させながら、推定次数ｑごとに推定結果を評価し、それらの評価結果に基づいて最適推定次数を特定するものである。よって、ＦＩＤ信号に応じて計算上使用する推定次数を適応的に選択できるから、不適切な推定次数の設定に起因する推定精度の低下という問題は生じない。すなわち、高次Prony法をそのまま適用した場合に生じる問題を克服できる。第２に、本手法は、評価に際して推定結果の全部を常に評価するのではなく、推定結果の中から選抜された優良部分だけが評価されるようにし、また、最終的な評価結果も選抜された優良部分によって構成されるようにしたものである。よって、真値から大きく外れたパラメータが推定されている場合であっても同時に真値に近いパラメータが推定されているならば、その後者の良い部分のみを評価しまた採用することができる。すなわち、本手法は、高次Prony法を基礎としつつも、同法とは異なり、構成次数つまり注目ピーク数と同数の推定結果を選抜して評価及び採用するものである。

（６）本実施形態に係る方法についての検証
本実施形態に係る方法の有効性を検証するため、２つの計算機シミュレーション及び実測信号を利用した実験を行った。

第１のシミュレーションでは、図２に示したものと同様の模擬ＦＩＤ信号に対して本手法を適用し、推定次数に対する評価量Ｊの推移を求めると共に、本手法の繰り返し適用により最適推定次数の分布を求めた。第２のシミュレーションでは、エタノール中の¹HのＮＭＲ測定を前提とする模擬FID信号に対して本手法を適用し、原子団（置換基等）ごとの強度を推定した。

まず第１のシミュレーションについて説明する。構成次数を４とし、選抜領域Ｒを注目ピーク波形の中心周波数を基準としてその上下に0．1Hzの範囲とした。そのような前提で推定モデルを利用して、推定次数ごとに、模擬ＦＩＤ信号を構成する信号成分ごとの強度を求め、それらの合計値を算出した。そのような処理を外乱を変えつつ１００回試行した。その結果、真値の2に対して、推定した強度の平均値は1．9997となり、標準偏差は0．0156であった。

図５には第１のシミュレーションの結果が示されている。図５の上段には、推定次数の変化に伴う評価値Ｊの変化が示されており、図５の下段には最適推定次数の頻度分布が示されている。図５から、推定精度が悪いと思われる低次範囲及び高次範囲においてＪの値が大きくなっていること、及び、推定次数ｑ＝１９周辺において最適な推定次数が決定される頻度が多いこと、を理解できる。最適推定次数の頻度分布についての平均値及び標準偏差から理解できるように、本手法を用いることで、適切でない推定次数を設定した場合に生じるパラメータ推定精度の悪化を効果的に回避でき、つまり強度推定精度を向上できる。

次に第２のシミュレーションについて説明する。サンプリング周波数を10kHzとし、観測核の共鳴周波数を400MHzとし、観測中心周波数を5ppmとし、データポイント数を2¹⁵とした。図６に示すように、模擬ＦＩＤ信号に基づくスペクトルには、注目する３つのピーク波形（注目ピーク又は注目ピーク群）であるＯＨ（Ａ）（図６中のピーク波形（Ａ））、ＣＨ_２（Ｂ）（図６中のピーク波形（Ｂ））、ＣＨ_３（Ｃ）（図６中のピーク波形（Ｃ））が含まれ、それらを生じさせたパラメータセットは以下の表１に示すとおりである。

化合物の定量解析においては、個々のピークごとの強度よりも原子団（例えば置換基）ごとの強度Cgroupが重要である。よって、以下においては原子団の強度Cgroupの推定精度について検討する。

ＯＨ（Ａ）、ＣＨ_２（Ｂ）及びＣＨ_３（Ｃ）に対しては、それぞれ構成次数ｐとして１，４，３を設定した。また、それらに対しては、ピーク波形の中心周波数として、それぞれ4.664、3.516、1.045ppmを設定した。それらに対する選抜領域Rとして、各ピーク波形の中心周波数を基準としてその上下に0.1ppmの範囲とした。但し、計算時間を削減するために、前処理として、帯域通過フィルタによる信号抽出を行った。フィルタの通過帯域は、各原子団についての中心周波数を基準としてその上下に約0.25ppmの範囲とした。なお、フィルタ処理は時間領域で行った。

以上の解析条件の下、模擬FID信号に本手法を適用し、構成次数分のピークの強度を推定すると共に、それらの総和としての原子団ごとの強度を算出した。以下の表２に外乱を変え、強度推定を１００回試行した結果を示す。表２の上段は、推定強度についての平均値及び標準偏差を示しており、表２の下段は、推定誤差ついての平均値及び標準偏差を示している。

表２から、Ｓ／Ｎが１００以上の場合において、原子団ごとの強度が誤差１％以内で推定できていることを確認できる。すなわち、本手法によれば、実用的に問題のない定量解析精度を得られる。

続いて、アセトアミノフェン中の¹Hに対するＮＭＲ測定で得た実測ＦＩＤ信号に対して本手法を適用し、分子内積分比（強度比）を推定する実験を行った。併せて、同じ信号に対して従来法としてのカーブフィッティング法及び積分法も適用し、それらによる推定結果との比較も行った。但し、本手法はＦＩＤ信号である時系列データそれ自体に対して適用し、従来法はＦＩＤ信号の吸収スペクトルに対して適用した。サンプリング周波数を10.016kHzとし、観測核の共鳴周波数を399.78MHzとし、観測中心周波数を5ppmとし、データポイント数を2¹⁵とし、積算回数を8回とした。図７は、上記の実測ＦＩＤ信号のＦＦＴスペクトルを示すものである。

解析条件は次のとおりである。本手法では、ＣＨ（Ａ）、ＣＨ（Ｂ）、ＣＨ_３（Ｃ）に対してそれぞれ構成次数pとして6、6、1を適用し、また、それらの原子団についての中心周波数をそれぞれ7.295、6.63、1.936ppmとし、それらの原子団に対する選抜領域として、それぞれ中心周波数を基準としてその上下に0.025、0.05、0.1 ppmの範囲を設定した。なお、原子団ごとに、前処理として、中心周波数から上下に約0.25ppmの範囲で信号抽出を行った。また、データポイントの序盤21点分はダウンサンプリング処理のためのディジタルフィルタの遅延部分に相当する。このため、22点目以降のデータポイントを解析で用いた。

参考までに、カーブフィッティング法の解析条件について説明すると、初期パラメータは以下の表３に示すとおりである。

左側の条件１は、まず吸収スペクトルからピークごとの周波数を求め、続いて線幅、強度の順でそれらを求めて、これらを初期パラメータとしたものである。右側の条件２は、周波数については条件１と同じ値を適用し、線幅及び強度については、条件１における最下段に示されている条件つまり1.43E-03及び2.44E-03を全てのピークに対して適用し、これらを初期パラメータとしたものである。なお、パラメータの最適化を行う範囲として、原子団ごとに中心周波数を基準としてその上下に約0.05ppmの範囲とした。

積分法においては、原子団ごとに、そのピーク波形の両端（線幅）を基準とし、そのピーク波形の中心周波数を中心として、線幅の３２倍及び６４倍となるように積分範囲を設定した。但し、線幅や周波数の値は表３における条件1中の値を用いた。

以下の表４には、本手法（提案法）、カーブフィッティング法及び積分法による推定結果が示されている。但し、それらの方法の適用に当たり、ＣＨとＣＨ_３の強度比は２：３であると仮定し、ＣＨ_３（Ｃ）（図７参照）の強度を3として、それに対するＣＨ（Ａ）（図７参照）、ＣＨ（Ｂ）（図７参照）の強度を推定した。

表４の上段には、推定強度が示されており、表４の下段には、推定誤差が示されている。表４に示されているように、カーブフィッティング法の推定誤差が一番大きく、条件１、２でそれぞれ１０％、４％以上の誤差が生じた。さらに、同方法では条件１の結果と条件２の結果との間に１０％程度の差が見られる。このことからカーブフィッティング法が初期値に大きく依存する手法であることを確認できる。積分法は、積分範囲が線幅の６４倍の場合には１％以内で強度が推定されたが、積分範囲が線幅の３２倍の場合には１％を超える誤差が生じた。このことから積分法の推定精度は、積分範囲に大きく依存することがわかる。

これに対して、本手法では、誤差１％以内で強度を推定できており、本手法は定量解析に対して十分に適用可能である。さらに本手法は、従来法において必要であった位相補正処理を行わずに、ＦＩＤ信号の強度を推定できるという利点を得られる。このような点から、本手法は、従来の２つの方法に比べて、客観性の高い手法と言える。また、表４より、選抜範囲の大きさの違いが推定精度に影響を与えないことを理解できる。他の実験データについても同様の結果が得られるとまでは言えないが、選抜範囲はパラメータ列の適切な選抜を行うためのものであるため、適切な選抜ができる程度に選抜範囲の設定を行えれば十分である。よって、選抜範囲は厳密に決める必要性は少なく、それは注目しているピーク群を含むようにある程度広くとればよい。

（７）実施形態に係るＮＭＲ信号処理システム
図８には、本実施形態に係るＮＭＲ信号処理システムがブロック図として示されている。このシステムは、ＮＭＲ測定システムの一部に組み込まれ、あるいは、単体の情報処理装置として構成されるものである。図８に示されている各機能は基本的にソフトウエア処理により実現されている。但し、その全部又は一部が専用ハードウエアモジュールによって構成されてもよい。

図８において、ＮＭＲ信号処理システム４０は、ＮＭＲ測定によって取得された受信信号としてのＦＩＤ信号を処理するシステムであり、それは例えばパーソナルコンピュータによって構成される。そのＦＩＤ信号は一般に直交検波後の複素信号（デジタル信号）である。ＦＦＴ演算部は、従来同様に、ＦＩＤ信号に対するＦＦＴ演算によってスペクトルを生成する。スペクトルを表すデータは表示処理部４６を経由して図示されていない表示部へ送られる。

ＮＭＲ信号処理システム４０は、上述した本実施形態に係る方法を実行するパラメータ推定モジュール４８を含む。パラメータ推定モジュール４８は、既知の化合物を前提として、換言すれば、注目ピークの個数が既知であることを前提として、選抜再構成評価法が組み込まれた高次Prony法に基づいて、個々の注目ピークごとの強度を演算するものである。

パラメータ推定モジュール４８について具体的に説明する。ＢＰＦ（バンドパスフィルタ）５６は、ＦＩＤ信号中における、注目ピーク波形を包含する一定帯域内の周波数成分を抽出するものである。その帯域は、後述する選抜領域を包含するより大きな範囲である。このＢＰＦ５６は演算量削減等の目的から必要に応じて前処理部として設けられる。ＢＰＦ５６における周波数特性、つまり通過帯域の中心周波数及び帯域幅はフィルタ特性設定部５２によって定められる。そのフィルタ特性設定部５２は、図示の例において、選抜領域Ｒの両端を特定する周波数ω1，ω2に基づいてフィルタ特性を定めている。

推定部５８は、上記の高次Prony法を実行するユニットである。すなわち、信号モデルに基づいて、推定次数ｑごとに、ｑ個のパラメータ（極及び複素強度の組）を演算するものである。推定次数ｑは、図示の例において、初期値ｐから最大値ｑ_maxまで１つ刻みで順次変更される。個々の推定次数でパラメータ列の推定演算が実行される。そのような推定演算は上述した数学的な信号モデル６０を基礎とするものである。演算部６２は、この例において、ユーザー指定で又は自動的に決定された係数ｍを構成次数ｐに乗算することにより、上限値としてｍ×ｐを設定する。係数ｍとして例えば１０が設定される。

推定結果記憶部６４はメモリによって構成される。すなわち、そこには、個々の推定次数ｑごとに演算されたパラメータ列が記憶される（図４の（Ａ）を参照）。各パラメータ列はｑ個のパラメータからなるものである。

選抜部６６は、推定結果の評価に際して、上述した２種類の選抜処理を行うものである。具体的には、まず、推定次数ｑごとに推定されたｑ個のパラメータの中から、共鳴周波数が選抜領域Ｒ内に属する条件を満たすパラメータを選抜する。それらは選抜パラメータ列を構成する。次に、選抜部６６は、選抜パラメータ列を構成するパラメータ数がｐよりも大きい場合に、選抜パラメータ列の中から真のピークである可能性の高い順に上位ｐ個を取り出して選抜パラメータ列を構成する。これにより、推定次数ｑごとに基本的にｐ個のパラメータが選抜されることになり、それらが評価対象となる（図４の（Ｂ）を参照）。但し、例外的にｐ個未満のパラメータが選抜されることもある。なお、選抜結果については図示されていない記憶部に格納され、あるいは、記憶部６４上に選抜結果を表す情報が格納される。

評価値演算部６８は、推定次数ｑごとに、選抜パラメータ列を上記式（24）に与えて評価値を演算する（図４の（Ｃ）を参照）。この評価値は、図示されていない記憶部に格納され、推定次数ｑごとに演算された選抜パラメータ列が真のピーク群を表している可能性の度合いを示すものである。よって、そのような評価値を演算できる限りにおいて各種の評価値演算式を用いることができる。決定部７０は、推定次数ｑを初期値から最大値まで変化させることにより得られた複数の評価値の中で最良評価値を決定する。この決定により同時に最適推定次数が決定される。決定部７０は、最良評価値つまり最適推定次数に対応する最良選抜パラメータ列を推定結果記憶部６４から取り出して表示処理部へ出力する。その最良選抜パラメータ列は、原則として構成次数ｐと同数のパラメータからなるものである。それらはいずれも極と複素強度の組であるから、これによって原則としてｐ個の複素強度が特定され、それらからｐ個の強度が特定される。例外として、最良選抜パラメータ列がｐ未満の個数のパラメータにより構成される場合、それと同数の強度が求められる。その場合でも各極から共鳴周波数を特定できるので、求められた個々の強度と注目ピークとの関係を対応付けることは可能である。

表示処理部４６は、以上のように求められたパラメータ列、複数の強度等を表す推定結果画像を生成し、それを図示されていない表示部へ送る。その場合、注目ピークを特定する情報とそれについて推定された強度等とが対応付けられつつ表示されるのが望ましい。その場合、ＦＦＴ演算部４２で演算されたスペクトル（実測スペクトル）が同時に表示されてもよい。推定されたパラメータ列をＦＩＤ信号モデルに代入することにより実測ＦＩＤ信号を模擬する推定ＦＩＤ信号を構成し、その波形を表示したり、そのスペクトル（推定スペクトル）を表示したりしてもよい。後者の場合において、実測スペクトルと推定スペクトルとを同時に表示すれば推定演算の正しさを視覚的に確認することができる。

設定用テーブル５０は必要に応じて設けられるものであり、それは記憶部により構成されている。設定用テーブル５０には、化合物、注目原子団等の推定対象と、それに対応する初期値セットと、が対応付けられており、いずれかの推定対象が指定されると、設定用テーブル５０から、その推定対象に対応した初期値セットが出力される。その初期値セットは、構成次数ｐ、選抜領域Ｒの下限及び上限を特定する情報ω1、ω2を含む。更に係数ｍが特定されてもよい。もっとも、それらの初期値をユーザーが指定してもよく、所定の関数式に基づいて演算してもよい。ｑＮＭＲを実行する上で使い勝手のよいシステムが構成されるのが望ましい。

上記の実施形態では、注目ピーク数つまり構成次数が既知であったが、注目ピーク数つまり構成次数の特定を求めない変形例も考えられる。例えば、注目ピークの個数を試行的に変更しながら注目ピーク数ごとに本手法を繰り返し適用してもよい。また、上記の実施形態では、個々の注目ピークについての中心周波数の情報が積極的に利用されていなかったが、事前に特定されている中心周波数と推定された中心周波数とを比較する等の変形例も考えられる。そのような比較を評価値の演算に取り込むようにしてもよい。

図９及び図１０にはＦＩＤ信号のスペクトルが示されており、それらには選抜領域の設定例が示されている。図９に示す第１例において、選抜領域Ｒ１は、注目ピーク波形を含むように設定されている。このような選抜範囲Ｒ１を用いた選抜処理により、関心のないピーク波形あるいは妨害ピーク波形を推定ないし評価の対象から除外して、注目ピーク波形についての推定精度を高められる。選抜領域Ｒ１の他、２つの選抜領域Ｒ２、Ｒ３が設定されており、注目ピーク波形の順次選択に伴ってそれらの選抜領域が順次選択される。ＢＰＦの通過帯域１００には３つの選抜領域Ｒ１，Ｒ２，Ｒ３の全体が含まれている。

図１０に示す第２例においては、２つのサブ領域Ｒ−１、Ｒ−２が設定されている。それらはそれぞれ１又は複数の注目ピークを含んでおり、それら全体として単一の選抜領域Ｒを構成している。ＢＰＦの通過帯域１０２は、２つのサブ領域Ｒ−１、Ｒ−２を包含するように設定されている。この第２例では、離間した２つの２つのサブ領域Ｒ−１、Ｒ−２内の複数の注目ピークについて１回の推定演算でそれぞれの強度を推定することが可能である。換言すれば、一定範囲にわたる選抜領域中において推定又は評価の対象としたくない部分を演算対象から除外することができる。例えば、妨害ピークが存在しているような場合にそれを推定、評価の対象から除外できる。

以上説明したように、本実施形態によれば、特にｑＮＭＲにおいて、注目ピークの強度を高精度に推定可能である。従来の積分法では、積分範囲の設定に困難が伴っていたが、本実施形態に手法によればそのような問題は生じない。また、従来のカーブフィッティング法で見られるような初期パラメータ依存性という問題も生じない。よって、より客観的な手法を提供できる。

(８）他の評価方法
上記の式(24)では時間領域上の実測ＦＩＤ信号及び再構成信号に基づいて評価値が演算されていたが、周波数領域上の実測ＦＩＤ信号及び再構成信号に基づいて評価値を演算することも可能である。すなわち、時間領域上の実測ＦＩＤ信号及び再構成信号をそれぞれフーリエ変換すれば、周波数領域上の実測ＦＩＤ信号及び再構成信号を得ることができ、それらの残差から評価値を求めることが可能である。フーリエ変換後の周波数領域上の実測ＦＩＤ信号及び再構成信号をそれぞれ以下のように表現する。但し、n=0,1,…,N‐1である。

fsをサンプリング周波数とした場合、周波数領域での標本点w_nは以下のように表現される。

以上に基づき周波数領域での評価値Ｊが例えば以下のように定義される。

２つの信号間の残差エネルギーは、時間領域であれ周波数領域であれ本質的には同じであるため、式（24）で定義される評価値と式（25）で定義される評価値は等価である。
但し、周波数領域における評価値が有益に機能する場合がある。それは全帯域ではなく局所的な帯域内のみで評価を行いたい場合である。そのような評価は時間領域で行うことは難しい。周波数領域において、評価する周波数範囲を指定した上で、その範囲内のスペクトルから周波数領域での評価量が以下のように定義される。

但し、ω_a及びω_bは、それぞれ評価範囲の下限と上限である。上記式(26)を用いることで、局所的な帯域に限ったところでの評価が可能となる。そのような評価法は、特に、観測周波数領域内に他の信号成分が含まれるような場合に有効となる。

これについて図１１を用いて説明する。図１１には、８個の信号成分（つまりp=8）からなるＦＩＤ信号のスペクトルが示されている。図示されてはいないが、これと同様の実測ＦＩＤ信号のスペクトルが別途生成される。符号１０４は注目ピークであり、それが含まれるようにフィルタ通過帯域１０６が設定されている。符号１０８は注目ピーク１０４を中心とした選抜範囲を示している。符号１１０は注目ピークを中心とした評価範囲を示している。実測ＦＩＤ信号及び再構成信号とも、評価範囲１１０に含まれる周波数成分だけが評価値の演算で利用され、それ以外の周波数成分は評価値の演算では利用されない。

このように上記式(26)による評価値の演算によれば、評価したい周波数範囲に限って評価を行えるから、他の信号成分による影響を効果的に除外又は軽減でき、評価精度を高めることが可能である。フィルタ通過範囲１０６の設定によって上記同様の範囲限定（特に注目ピークに近接する他の信号成分の除外）を行うことは非常に難しいが、周波数領域での処理であればそれを容易に実現できる。図１１に示す例では、選抜範囲１０８と評価範囲１１０とは別々に設定されていたが、両者を一致させることも可能である。

４８強度推定モジュール、５６バンドパスフィルタ、５８推定部、６４推定結果記憶部、６６選抜部、６８評価値演算部、７０決定部。

Claims

ＮＭＲ測定で得られる実信号をｑ個の信号成分の和として表現する数学的モデルに基づき、当該数学的モデルにおける次数ｑを可変しながら、次数ｑごとに、前記実信号に基づいて前記ｑ個の信号成分を定義するｑ個のパラメータを推定する推定手段と、
周波数軸上におけるｐ個の注目ピークに応じて定められた選抜条件に従って、次数ｑごとに、前記推定手段が推定したｑ個のパラメータの中から評価対象とする上位ｐ個のパラメータを選抜する機能を有する選抜手段と、
前記選抜手段の選抜結果に基づいて、次数ｑごとに評価値を演算する評価手段と、
前記評価手段が演算した次数ｑごとの評価値の中から最良評価値を特定することにより最適次数を決定した上で前記最適次数に対応する上位ｐ個のパラメータを特定する機能を有する決定手段と、
を含むことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記実信号は時間軸上のＦＩＤ信号であり、
前記数学的モデルは前記ＦＩＤ信号をｑ個の指数減衰正弦波の和として表現するモデルである、
ことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記各パラメータは極及び複素強度の組である、
ことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記選抜条件には、前記ｐ個の注目ピークを包含する選抜領域を利用して次数ｑごとに選抜を行うための第１選抜条件が含まれ、
前記選抜手段は、前記推定手段が推定したｑ個のパラメータの中から前記第１選抜条件を満たす複数のパラメータを選抜する第１選抜手段を含む、
ことを特徴とするＮＭＲ信号処理システム。
請求項４記載のシステムにおいて、
前記選抜条件には、真のピークである可能性が高い順で上位ｐ個のパラメータを選抜するための第２選抜条件が含まれ、
前記選抜手段は、前記第１選抜手段が選抜した前記第１選抜条件を満たす複数のパラメータの中から前記第２選抜条件に従って上位ｐ個のパラメータを選抜する第２選抜手段を含む、
ことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記次数ｑの可変範囲の下限はｐであり、
前記次数ｑの可変範囲の上限はｐに基づいて設定される、
ことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記評価手段は、前記推定された複数のパラメータを前記数学的モデルに与えることにより生成される合成信号と、前記実信号と、を比較することにより、評価値を演算する、
ことを特徴とするＮＭＲ信号処理システム。
請求項１記載のシステムにおいて、
前記実信号に対して前処理としての帯域通過処理を施すフィルタ手段を含み、
前記推定手段は前記帯域通過処理後の実信号に基づいて前記複数のパラメータの推定を行う、
ことを特徴とするＮＭＲ信号処理システム。
ＮＭＲ信号処理システムにおいて実行されるプログラムであって、
ＮＭＲ測定で得られる実信号を入力し、当該実信号をｑ個の信号成分の和として表現する数学的モデルに基づき、当該数学的モデルにおける次数ｑを可変しながら、次数ｑごとに、前記実信号に基づいて前記ｑ個の信号成分を定義するｑ個のパラメータを推定し、それらの推定結果を第１記憶部に格納する機能と、
周波数軸上におけるｐ個の注目ピークに応じて定められた選抜条件に従って、次数ｑごとに、前記推定されたｑ個のパラメータの中から評価対象とする上位ｐ個のパラメータを選抜する機能と、
次数ｑごとに選抜された上位ｐ個のパラメータを評価し、次数ｑごとに評価値を第２記憶部に格納する機能と、
次数ｑごとの評価値の中から最良評価値を特定することにより最適次数を決定した上で前記最適次数に対応する上位ｐ個のパラメータを特定する機能と、
を含むことを特徴とするプログラム。