(57)【特許請求の範囲】
【請求項１】
装置により収集された信号であって、当該装置による当該信号の収集により歪みが生じた当該信号のサンプルの系列から、当該信号の周波数領域表現を推定するための信号処理方法であって、
前記サンプルの系列を取得する工程と、
信号領域における装置関数のシフトされたバージョンによってそれぞれが定義される複数の基底関数を含む、基底関数の集合を、信号領域の振動成分を含む複素指数関数にフィッティングする係数の集合を取得する工程と、
前記係数の集合を使用して前記サンプルの系列をスケーリングし、スケーリングされたサンプルの総和を求めることによって、前記信号の前記周波数領域表現を推定する工程と、を含み、
前記装置関数は、前記装置による信号の前記収集が当該信号に対して有する、前記歪みの影響を表し、かつ、前記複素指数関数に含まれる前記振動成分の周波数であって当該複素指数関数に関連付けられた周波数における、前記周波数領域の前記歪みの影響の特性に基づいている
ことを特徴とする方法。
【請求項２】
前記サンプルは、前記信号領域においてサンプル間間隔により分離されており、
各基底関数は、前記信号領域において整数のサンプル間間隔によりシフトされ、
前記フィッティングは、前記複数の基底関数により定義される空間への射影を含む
ことを特徴とする請求項１に記載の方法。
【請求項３】
前記サンプルは、前記信号領域においてサンプル間間隔により分離されており、
各基底関数は、前記信号領域において非整数のサンプル間間隔によりシフトされ、
前記フィッティングは、メリット関数を減少させることにより実行される
ことを特徴とする請求項１に記載の方法。
【請求項４】
前記係数は、前記複数の基底関数のうちの、それぞれ対応する１つの基底関数の周波数領域表現を含む
ことを特徴とする請求項１から３のいずれか１項に記載の方法。
【請求項５】
前記複素指数関数は、信号領域の振動成分を含み、前記装置関数が特徴付けられる周波数は、前記振動成分の周波数を含む
ことを特徴とする請求項１から４のいずれか１項に記載の方法。
【請求項６】
（ａ）前記周波数は、前記周波数領域において等間隔に位置している、
（ｂ）前記周波数は、周波数間の間隔を増加させながら、前記周波数における前記装置関数の大きさの変動を減少させるように選択される、及び／又は
（ｃ）前記信号は、イノベーション率有限（ＦＲＩ）信号の種類に属し、前記周波数は、前記種類の信号のイノベーションパラメータについての最尤推定器の性能に基づいて選択される
ことを特徴とする請求項５に記載の方法。
【請求項７】
前記周波数領域表現に基づいて、前記信号の推定値を提供する工程を更に含む
ことを特徴とする請求項１から６のいずれか１項に記載の方法。
【請求項８】
前記複素指数関数の前記係数及び前記複数の基底関数によって提供される前記周波数領域表現の誤差を取得する工程と、
前記誤差に基づいて、前記信号の前記周波数領域表現を再推定するかどうかを判定する工程と、
を含むことを特徴とする請求項１から７のいずれか１項に記載の方法。
【請求項９】
前記信号の前記周波数領域表現を再推定するかどうかを判定する工程は、
原信号における雑音の推定値を取得する工程と、
前記推定された雑音と前記複素指数関数の前記周波数領域表現における前記誤差とに基づく比較を実行する工程と、
を含むことを特徴とする請求項８に記載の方法。
【請求項１０】
前記信号の前記周波数領域表現の再推定は、
前記信号の前記周波数領域の推定表現のエネルギーを最小化するように構成されるフィルタを判定し、前記フィルタに基づいて前記周波数領域表現を調節すること、及び
行列ペンシル問題から固有値を導出することによって、前記信号領域における前記信号の位置を表すパラメータの推定値を取得すること、
のうちの１つを含むことを特徴とする請求項８又は９に記載の方法。
【請求項１１】
前記フィルタに基づいて調整することは、前記フィルタの根を使用して、前記信号領域における前記信号の位置を表すパラメータの推定値を取得することを含む
ことを特徴とする請求項１０に記載の方法。
【請求項１２】
前記信号は画像を含み、
前記方法は、
前記複素指数関数に関連付けられた前記周波数に基づいて、前記画像の前記周波数領域表現をスケーリングする工程と、
前記画像の前記周波数領域表現のエネルギーを最小化するように構成されるフィルタ係数ｂkjの集合を判定する工程と、
前記フィルタ係数から前記画像内の境界を判定する工程と、
を更に含むことを特徴とする請求項１から１１のいずれか１項に記載の方法。
【請求項１３】
前記境界と第２の画像について取得される境界とに基づいて、前記画像と前記第２の画像とを共同位置合わせするための変換を決定する工程を更に含む
ことを特徴とする請求項１２に記載の方法。
【請求項１４】
画像解像度を向上させる方法であって、
請求項１２に記載の方法に従って、複数の画像の各々について少なくとも１つの境界を特定する工程と、
境界毎に、各境界と選択された基準境界とを共同位置合わせするための変換を決定する工程と、
前記複数の画像及び前記変換に基づいて、合成画像を決定する工程と、
を含むことを特徴とする方法。
【請求項１５】
前記合成画像を決定する工程は、前記複数の画像のうちの少なくとも１つの画像からのサンプル及び当該画像に関連付けられた前記変換に基づいて、前記合成画像における画素の値を決定する工程を含む
ことを特徴とする請求項１４に記載の方法。
【請求項１６】
前記装置関数は、前記装置のサンプリングカーネルを含む
ことを特徴とする請求項１から１５のいずれか１項に記載の方法。
【請求項１７】
前記装置は光学系を含み、前記サンプリングカーネルは点広がり関数を含む
ことを特徴とする請求項１６に記載の方法。
【請求項１８】
装置により収集された信号のサンプルの系列から、当該信号の周波数領域表現を提供するデジタル信号処理装置であって、前記信号には前記装置による当該信号の収集により歪みが生じ、前記デジタル信号処理装置は、
前記信号の前記サンプルの系列を取得するように構成されたデータ取得器と、
信号領域におけるサンプル数だけシフトされた装置関数によって定義される基底関数の集合から、前記信号領域の振動成分を含む複素指数関数を近似するように選択される係数の集合を使用して、前記サンプルの系列をスケーリングするように構成されたデータスケーラと、
前記スケーリングされたサンプルの総和を求めることによって、前記信号の前記周波数領域表現を提供するデータ提供器と、を備え、
前記装置関数は、前記装置による信号の収集が当該信号に対して有する、前記歪みの影響を表し、かつ、前記複素指数関数に含まれる前記振動成分の周波数であって当該複素指数関数に関連付けられた周波数における、前記周波数領域の前記歪みの影響の特性に基づいている
ことを特徴とするデジタル信号処理装置。
【請求項１９】
請求項１から１７のいずれか１項に記載の方法を実行するよう構成された、請求項１８に記載のデジタル信号処理装置。
【請求項２０】
信号収集装置を備え、かつ、較正方法に従って構成されたデジタル信号処理装置であって、前記較正方法は、
前記信号収集装置の伝達関数の推定値を取得する工程と、
前記周波数領域において等間隔に位置する周波数の集合を選択する工程と、
前記選択された周波数における前記伝達関数のサンプルに基づく係数を、前記信号収集装置のメモリに格納する工程と、を含む
ことを特徴とする請求項１８又は１９に記載のデジタル信号処理装置。
【請求項２１】
前記周波数は、周波数間の間隔を増加しながら、前記周波数における前記装置関数の大きさの変動を減少させるように選択される
ことを特徴とする請求項２０に記載のデジタル信号処理装置。
【請求項２２】
請求項１から１７のいずれか１項に記載の方法を実行するようにプロセッサをプログラムするよう動作可能なコンピュータプログラム命令を含むコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号のサンプリング及び再構成に関し、特に、イノベーション率有限（ＦＲＩ：finite rate of innovation）信号等の帯域制限のない信号のサンプリング及び再構成に関する。
【背景技術】
【０００２】
全サンプリングの基本的な問題は、原信号のサンプルの離散集合から当該原信号を推定することである。通常、信号には、例えば実際の信号と測定値を収集する計器のサンプリングカーネルとの畳み込みにより、当該計器による歪みが生じる。殆どの信号処理において、この歪んだ信号が離散サンプルに変換される。サンプルの信号領域の間隔により、測定の周波数領域の帯域幅が制限されることは既知である。
【０００３】
何らかの特定の信号収集計器は、当該計器により実行される測定に影響を及ぼす物理的特徴を有する。それらの特徴は、サンプリング前に原信号に歪みを生じさせる場合がある。この歪みの特性は、装置関数により特徴付けることができる。そのような装置関数の一例は、光学系の点広がり関数である。この装置関数は、サンプリング前の原信号に対して計器が有する歪み効果を表現する。
【０００４】
一般に、信号の測定のサンプルに基づく信号の推定精度を向上させる２つのアプローチが存在する。第１のアプローチは、サンプリングの帯域幅を増加させる（サンプリング周波数を増加させる）ことである。第２のアプローチは、装置関数の歪み効果を減少させることであり、例えば光学系では、これは点広がり関数の空間範囲を減少させることにより行われうる。
【０００５】
イノベーション率有限（ＦＲＩ）信号等の特定の種類の信号に対して、信号復元問題に対する１つの理論的アプローチは、装置関数が指数関数再生カーネルとして既知である関数の種類に属することを想定している。この関数の種類は、そのようなカーネルのシフトされた表現の線形結合から任意の複素指数関数を再構成する、そのようなカーネルの能力によって定められる。
【０００６】
現実の実際の計器において、装置関数は、計器の実際の物理的制限により決定され、例えば光学系では、開口及び光学列の他の構成要素によって装置関数が定義される。従って、一般に、現実の計器の装置関数は指数関数再生カーネルではない。更に、指数関数再生カーネルを使用するＦＲＩ再構成から開発された理論的アプローチは、雑音（例えばモデル挙動からの原信号のずれ）が存在する場合に、安定性が非常に低いという欠点を有する。当然、雑音は実際のあらゆる物理的計器に伴う。従って、これらの理論的アプローチは、当該技術分野において長年既知であったが、現実世界の実際のシステムには適用できない。
【発明の概要】
【０００７】
本開示の態様は、実際の現実世界のシステムにおける、及び測定雑音の存在下における、信号のサンプリング、周波数空間推定及び復元を対象としている。信号は、ＦＲＩ信号及びその他の種類の信号を含みうる。
【０００８】
本開示の実施形態は、信号の周波数領域推定値を提供するために装置関数の知識に基づいて当該信号の離散サンプルを修正することにより、この問題を解決する。特に、サンプルは、係数Ｃｍｎの集合によりスケーリングされる。それらの係数は、基底関数ｈｎの集合を複素指数関数にフィッティングするように選択される。それらの基底関数の各々は、信号領域においてシフトされた、装置関数の信号領域表現により定義される。
【０００９】
装置関数のバージョンは、選択された周波数αｍにおける、周波数領域の装置関数の特性（characterisation）に基づいている。これらの周波数は更に、係数Ｃｍｎを選択するために使用される複素指数関数に関連付けられる。係数Ｃｍｎを選択するために使用される複素指数関数は、ｅｘｐ｛−αｍｔ｝を含んでもよく、ｔは信号領域の変位（displacement）を示す。例えば時間領域関数の場合、パラメータｔは時間遅延を示し、空間関数の場合、ｔは時間遅延ではなく空間変位を示す。
【００１０】
周波数αｍは、ラプラス変換の意味で複素周波数であってもよく、その場合、複素周波数の虚数部は振動成分の周波数を示し、複素周波数の実数部は信号の減衰又は増大を示す。
【００１１】
本開示の実施形態は、上記の周波数を選択し、且つ、装置関数の推定値に基づいて計器についての対応する係数Ｃｍｎを生成することにより、信号処理装置１２を構成する方法に関連している。
【００１２】
本発明の理解を助けるために、発明の特定の実現例をいくつか選択して以下に説明する。それらのうちの第１の実現例は、パルス系列の飛行時間に基づく距離計測（range finding）に関する。パルス系列を含む信号の推定に関連するため、この例が選択されている。当業者に理解されるように、殆どの実際の信号及び特にパルス系列は、信号点の系列（例えば、ディラックのデルタ関数の系列）のタイミング及び振幅により表されうる。従って、この例は本開示を理解するための１つの方法をして提供されるが、本開示の装置及び方法はあらゆる信号の信号処理に適用可能であることが理解されるだろう。
【００１３】
本開示は、特にＦＲＩ信号に適用されるが、他の種類の信号に適用されてもよい。ＦＲＩ信号は、名称が示すように、単位時間毎に有限な数の自由度を有し、それらは周期的であってもよいし、あるいは非周期的であってもよい。ＦＲＩ信号の例には、ディラックのデルタ関数のストリーム、不均一なスプライン、正弦波及び区分的多項式が含まれる。
【００１４】
それらの種類の信号の各々は、「イノベーションパラメータ」により特徴付けられうる。例えばディラックのデルタ関数の系列のイノベーションパラメータは、デルタ関数の信号領域（例えば、タイミング又は位置）における位置ｔｋ及び振幅ａｋである。ほぼ全ての信号がデルタ関数の系列としての表現に変形可能であるため、この種のＦＲＩ信号を以下の説明の大部分の基礎として使用する。しかし、本開示の教示はあらゆる種類の信号に適用可能であることが理解されるだろう。
【００１５】
一態様において、装置関数（ｈ）により歪みが生じたサンプル（ｙｎ）の系列から信号（ｘ）の周波数領域表現（Ｘ）を推定する信号処理方法が提供される。本方法は、サンプル（ｙｎ）の系列を取得する工程と、信号領域における装置関数のシフトされたバージョンによってそれぞれが定義される複数の基底関数を含む、基底関数（ｈｎ）の集合を、複素指数関数にフィッティングする係数（Ｃｍｎ）の集合を取得する工程と、サンプル（ｙ）の系列及び係数（Ｃｍｎ）に基づいて、信号（ｘ）の周波数領域表現（Ｘ）を推定する工程と、を含み、装置関数の推定値は、複素指数関数に関連付けられた周波数における、周波数領域の装置関数の特性に基づいている。
【００１６】
一態様において、装置関数ｈにより歪みが生じた信号ｘに対応する、信号ｙ（ｔ）のサンプルｙｎの系列を取得するように構成されたデータ取得器と、サンプルの番号ｎだけシフトされた装置関数によって定義される基底関数の集合ｈｎから、複素指数関数を近似するように選択される係数Ｃｍｎを使用して、サンプルｙｎの系列をスケーリングするように構成されたデータスケーラと、当該スケーリングされたサンプルに基づいて、信号ｘの周波数領域表現Ｘを提供するデータ提供器と、を備えるデジタル信号処理装置が更に提供され、装置関数は、複素指数関数に関連付けられた周波数における、周波数領域の装置関数の特性に基づいている。
【００１７】
信号のサンプルは、信号領域においてサンプル間間隔又はサンプル周期により分離されてもよい。各基底関数は、信号領域において整数ｎのサンプル間間隔により他の基底関数に対してシフトされてもよく、それらの例において、フィッティングは、複数の基底関数により定義される空間への射影を含みうる。いくつかの例において、各基底関数は、信号領域において非整数のサンプル間間隔によりシフトされてもよく、それらの例において、フィッティングは、メリット関数を減少させることにより実行されうる。
【００１８】
係数（ｃｍｎ）は、複数の基底関数（ｈｎ）のうちの、それぞれ対応する１つの基底関数の周波数領域表現（Ｈｎ）を含んでもよい。複素指数関数は、信号領域の振動成分を含んでもよく、装置関数（ｈ）が特徴付けられる周波数は、振動成分の周波数を含む。複素指数関数は、信号領域の振動成分のみで構成されてもよい。これは、特に雑音に対してロバストな信号推定を提供することがわかっている。
【００１９】
いくつかの可能な例において、係数は、選択された周波数における基底関数ｈｎの周波数領域表現Ｈｎの逆数を含む。いくつかの可能な例において、係数は、基底関数ｈｎの周波数領域表現Ｈｎの反転バージョンを含む。係数は、信号領域における装置関数の自己相関のｚ変換によりスケーリングされてもよい。
【００２０】
選択される周波数は、周波数領域において等間隔に位置してもよく、いくつかの例において、周波数は、周波数間の間隔を増加させながら、当該周波数における装置関数の大きさの変動を減少させるように選択される。これにより、雑音に対する復元力（resilience）が改善することがわかっている。いくつかの可能な例において、信号はＫの自由度を有し、Ｋは整数であり、Ｋが整数である場合、ｎ＝０，１,..., Ｐ及びＰ≦２Ｋ−１である。
【００２１】
いくつかの可能な例において、推定される信号は、イノベーション率有限（ＦＲＩ）信号という種類に属し、周波数は、この種類の信号のイノベーションパラメータの最尤推定器の性能に基づいて選択される。これは、クラメール・ラオ限界を減少させることにより達成されうる。
【００２２】
いくつかの可能な例において、本方法は、複素指数関数の係数（ｃｍｎ）及び基底関数（ｈｎ）により提供される表現の誤差を取得する工程と、当該誤差に基づいて、信号の周波数領域表現（Ｘ）を再推定するかどうかを判定する工程と、を含む。信号の周波数領域表現（Ｘ）を再推定するかどうかの判定は、原信号における雑音の推定値を取得する工程と、推定された雑音と複素指数関数の表現における誤差とに基づく比較を実行する工程と、を含んでもよい。
【００２３】
いくつかの可能な例において、信号の周波数領域表現の再推定は、信号（ｘ）の周波数領域の推定表現（Ｘ）のエネルギーを最小化するように構成されるフィルタを判定する工程と、当該フィルタに基づいて周波数領域表現を調節する工程と、を含む。これは、消滅フィルタ（annihilating filter）への近似により提供されうる。いくつかの可能な例において、フィルタに基づいて調節する工程は、当該フィルタの根を使用して、信号領域における信号の位置を表すパラメータの推定値を取得する工程を含む。
【００２４】
消滅フィルタを使用する方法が使用される必要はなく、例えばいくつかの可能な例において、信号の周波数領域表現の再推定は、行列ペンシル問題から固有値を導出することにより、信号領域における信号の位置を表すパラメータの推定を取得する工程を含む。
【００２５】
いくつかの可能な例において、信号は画像を含み、本方法は、複素指数関数に関連付けられた周波数に基づいて、画像の周波数領域表現をスケーリングする工程と、画像の周波数領域表現のエネルギーを最小化するように構成されるフィルタ係数ｂｋｊの集合を判定する工程と、当該フィルタ係数から画像内の境界を判定する工程と、を更に含む。理解されるように、画像データは、各行の周波数領域表現を推定し、その後に各列の周波数領域表現を推定することにより処理されてもよい。
【００２６】
いくつかの可能な例において、画像処理方法は、画像の共同位置合わせ（co-registration）のために使用されてもよく、それらの可能な例において、本方法は、境界と第２の画像について取得された境界とに基づいて、画像と第２の画像とを共同位置合わせするための変換を決定する工程を含む。画像解像度を向上させる方法であって、複数の画像の各々について少なくとも１つの境界を特定する工程と、境界毎に、各境界と選択された基準境界とを共同位置合わせするための変換を決定する工程と、複数の画像及び変換に基づいて、合成画像を決定する工程と、を含む方法が更に提供される。合成画像を決定する工程は、複数の画像のうちの少なくとも１つの画像からのサンプル及び当該画像に関連付けられた変換に基づいて、合成画像における画素の値を決定する工程を含んでもよい。
【００２７】
一態様において、信号収集計器を備える信号収集デバイスを較正する方法であって、計器の伝達関数の推定値を取得する工程と、周波数領域において等間隔に位置する周波数の集合を選択する工程と、選択された周波数における伝達関数のサンプルに基づく係数（Ｃｍｎ）を、デバイスのメモリに格納する工程と、を含む方法が提供される。これらの周波数は、周波数間の間隔を増加させながら当該周波数における装置関数の大きさの変動を減少させるように選択されてもよい。いくつかの可能な例において、信号収集デバイスは、イノベーション率有限（ＦＲＩ）信号を収集するように適合しており、この種類の信号のイノベーションパラメータの最尤推定器の性能に基づいて、周波数が選択される。
【００２８】
付録１において説明するように、周波数は、それらの周波数に基づいて判定されるクラメール・ラオ限界を減少させるように選択されてもよい。
【図面の簡単な説明】
【００２９】
添付の図面を参照して、ほんの一例として本開示の実施形態を以下で説明する。
【図１】図１は、信号処理装置を備える距離計測装置を概略的に示す図である。
【図２】図２は、信号をサンプリングする処理を示す図である。
【図３】図３は、図１の信号処理装置１２を更に詳細に示す距離計測装置を概略的に示す図である。
【図４】図４は、図３に示すような信号処理装置１２の動作の１つの方法を示す図である。
【図５】図５は、図３の信号処理装置１２及び／又は図４に示す方法を使用して実現可能な画像処理方法及び装置を示す図である。
【図６】図６は、画像の解像度を向上させる方法を示す図である。
【図７】図７は、図６に示す方法の例を更に示す図である。
【発明を実施するための形態】
【００３０】
図１は、信号の飛行時間（time of flight）に基づいて距離を判定する距離計測装置２を示す。
【００３１】
距離計測装置２は、信号４を物体６に向け、当該物体から反射した信号８を受信し、装置２における送出信号４の送信と反射信号８の受信との間の時間遅延に基づいて、装置２と物体６との間の距離を判定するように構成される。
【００３２】
距離計測装置２は、送出信号４を物体６へ送信するように構成される信号提供器１０を備える。信号提供器１０は更に、当該送出信号をタイミング推定器に提供する。
【００３３】
距離計測装置２は、反射（到来）信号を検出する計器を更に備える。ＡＤＣ１８は、計器１６により受信された信号をデジタルサンプルに変換するために、検出計器１６と接続される。
【００３４】
ＡＤＣ１８は、信号処理装置１２（ＤＳＰ）に受信信号のサンプルを提供するために接続される。ＤＳＰ１２は、サンプルを処理し、処理したサンプルをタイミング推定器１４に提供するように構成される。タイミング推定器１４は、到来した反射信号のタイミングを推定し、送出信号と反射信号との間の遅延を判定するように構成される。これにより、距離計測装置は、信号の飛行時間に基づいて物体６までの距離を判定できる。
【００３５】
図２に示すように、計器により受信される反射信号８は、時間領域信号ｘ（ｔ）を含む。検出計器１６は、例えば遅延及び減衰を生じさせることで、この到来した信号に歪みを生じさせる。この歪みは、例えば計器１６のインパルス応答である装置関数によって表現されうる。検出計器によって信号ｘ（ｔ）に歪みが生じた信号ｙ（ｔ）は、歪んだ信号の離散的表現ｙ（ｎ）を提供するためにＡＤＣ１８によってサンプリングされる。
【００３６】
距離計測装置２と物体６との間の距離を判定するために、距離計測装置２は、信号ｙ（ｎ）の量子化表現に基づいて、送出信号４の送信と反射信号８の受信との間の遅延を判定する。この測定の分解能は、装置関数によって生じる歪み及びＡＤＣの分解能によって制限される。従って、本例における距離推定の問題は、例えば原信号のサンプルの離散集合から原信号を復元する方法である、上述した全サンプリングの基本的な問題に帰着する。
【００３７】
パルス系列に基づく距離計測の例は、本開示を理解するのに有用であり、パルス列を含む多くの実際の信号がデルタ関数の系列として表されうることが、当業者には理解されるだろう。従って、本例は距離推定との関連で説明されるが、原理はあらゆる信号に適用可能であり、特に、デルタ関数の系列として表されうる信号に適用可能である。
【００３８】
図３は、図１及び図２を参照して上述した距離計測装置の更に詳細な例を示す。図３に示す装置は、信号の測定値を取得する計器、測定信号のサンプルの系列を提供するＡＤＣ１８、信号処理装置１２及びタイミング推定器を備える。
【００３９】
信号処理装置１２は、ＡＤＣからサンプルを受信し、計器により測定された信号の周波数領域表現の推定値をタイミング推定器に提供するように構成される。当該信号の信号領域表現は、この周波数領域表現から導出可能である。従って、タイミング推定器１４は、信号処理装置からの出力に基づいて、送信信号と当該信号領域表現との間の時間領域遅延を推定できる。
【００４０】
信号処理装置１２は、図３の距離計測装置の残りの部分とは別に提供されてもよい。信号処理装置１２及びその動作方法は、他の状況に適用されてもよい。この信号処理装置１２及びその動作を以下で更に詳細に説明する。
【００４１】
信号処理装置１２は、ＡＤＣから信号のサンプルｙ（ｎ）を取得するデータ取得器２２を備える。データ取得器２２は、データスケーラ３０と接続され、データスケーラ３０は、係数提供器２０と接続される。データスケーラ３０は、信号処理装置からの出力を提供するデータ提供器３４と更に接続される。信号処理装置１２は、データ提供器３４、係数提供器２０及びデータ取得器２２と接続されるデータ調整器３２を更に備える。
【００４２】
データ取得器２２は、信号のサンプルを格納するためのメモリ２４と、ＡＤＣ１８等のデジタルサンプルの供給元からメモリ２４にデータが提供されるようにする入力インタフェース２６とを備える。データ取得器２２は、データスケーラ３０及びデータ調整器３２がメモリ２４から信号のサンプルを読み出せるようにする連結器２８を更に備える。
【００４３】
係数提供器２０は、係数Ｃｍｎを格納し、且つ、データスケーラに当該係数を提供するように構成されたメモリ２１を備える。係数Ｃｍｎは、例えば次式のように、基底関数の集合の線形結合を定義することによって複素指数関数を近似するように選択される。
0006397001.tif 000002
【００４４】
基底関数ｈｎ（ｔ）は、信号領域における装置関数のシフトされたバージョンｈｎにより定義される。ここで、ｎは、例えばサンプル周期の数ｎによる信号領域におけるシフトを示し、周波数αｍは、虚数部が振動成分の周波数を示し、実数部が信号の減衰又は増大（減少又は増加）を示す複素周波数であってもよい。
【００４５】
図３に示す例において、係数提供器２０は、以下の関係式に従って選択される係数Ｃｍｎを備える。
0006397001.tif 000003
ここで、Ｈ（αｍ）は、周波数αｍにおける装置関数の周波数領域バージョンｈｎを含む。周波数αｍは、それらの周波数における装置関数の大きさ|Ｈ（αｍ）|の変動を減少させる一方で周波数間の間隔を増加させるように選択される。上述のように、パラメータｎは、信号領域におけるシフトに関連する。
【００４６】
データスケーラ３０は、データ取得器２２を介して信号のサンプルｙ（ｎ）を取得するように構成される。データスケーラ３０は、係数Ｃｍｎを使用してそれらのサンプルをスケーリングし、以下の関係式に従って、スケーリングされたサンプルの総和を求めることで信号の周波数領域表現Ｘ（αｍ）を推定するように構成される。
0006397001.tif 000004
【００４７】
データスケーラ３０は、信号の周波数領域表現の当該推定値Ｘ（αｍ）をデータ提供器に提供するように構成される。本段階での完全性のために、当該推定値Ｘ（αｍ）は以下のように記述されうることを示すことができる（付録１を参照）。
0006397001.tif 000005
ここで、
0006397001.tif 000006
であり、ａｋ及びｔｋは、信号のイノベーションパラメータである（例えば、信号を表すデルタ関数系列のタイミングｔｋ及び振幅ａｋ）。ζｍは、以下に定義する推定値の誤差を表す。パラメータα０及びλは、以下の関係式に従って複素平面における周波数αｍを定義する複素パラメータである。
αｍ＝α０＋ｍλ
ここで、ｍは０≦ｍ≦２Ｋ−１の整数であり、Ｋは信号の自由度の数を定義する。
【００４８】
データ提供器３４は、データを格納するためのメモリ３５と、データスケーラ３０及びデータ調整器３２がメモリ３５との間でデータの読み出し及び書き込みを行えるようにする連結器３６とを備える。データ提供器３４は、データが出力として例えばタイミング推定器に提供されるようにする出力インタフェース３８を更に備える。
【００４９】
データ調整器３２は、データ取得器から測定雑音の推定値Ｅを取得し、データ提供器から推定値Ｘ（αｍ）を取得するように構成される。データ調整器３２は、以下の関係式に従って信号の周波数領域表現の推定値Ｘ（αｍ）における誤差ζｍを推定するように更に構成される。
0006397001.tif 000007
ここで、推定される信号は、Ｋ個のデルタ関数から成る系列と考えられてもよく、ｔｋは、信号領域におけるそれらのデルタ関数の位置（例えば、それらのタイミング）を表し、ａｋは、それらの振幅を表す。振幅ａｋは、上記の式４において定義された関係式に従って、周波数領域表現Ｘ（αｍ）から（例えば、ζｍ＝０を設定し、且つ、Ｘ（αｍ）の現在の推定値に基づいて上記の式４を解くことによって）推定可能である。
【００５０】
式５におけるパラメータεｍは、指数関数ｅｘｐ｛αｍ｝と、係数Ｃｍｎによりスケーリングされた、信号領域における装置関数のシフトされたバージョンｈｎ（ｔ）から組み立てられたモデルとの間のミスマッチを表す。すなわち、次式で表される。
0006397001.tif 000008
【００５１】
データ調整器３２は、誤差の推定値ζｍと測定雑音の推定値Ｅとを比較し、当該比較に基づいて、信号を再推定するかどうかを判定するように構成される。
【００５２】
信号が再推定される場合、データ調整器３２は、例えば畳み込み積分Ｆ＊Ｘを０に向けて調節するフィルタＦの係数を選択することにより、周波数領域推定表現Ｘ（αｍ）のエネルギーを減少するフィルタＦを判定することで信号を再推定するように動作可能である。これは、Y.Hua及びT.K.Sakar、「Matrix Pencil Method for Estimating Parameters of Exponentially Damped Undamped Sinusoids in Noise」、IEEE Transactions on Acoustics、Speech and Signal Processing、vol.38、814〜824ページ、1990年5月において定義される行列ペンシル法を使用して、あるいは、M.Vetterli、P.Marliliano.及びT.Blu.、「Sampling signals with finite rate of innovation」、IEEE Transactions on Signal Processing、vol.50、1417〜1428ページ、2002年6月において定義される消滅フィルタ法を使用して実行されうる。
【００５３】
データ調整器３２は、フィルタＦのｚ変換の根ｕｋ（例えば、Ｆのｚ変換が０の値になる複素平面内の点）を取得し、以下の関係式に基づいてＫ個のデルタ関数から成る系列のタイミングｔｋの新しい推定値を判定するように構成される。
0006397001.tif 000009
ここで、Ｔは、サンプリング周期であり、λは、複素平面における周波数αｍの間隔を表す。データ調整器３２は、タイミングｔｋの新しい推定値に基づいて、振幅ａｋの新しい推定値を判定するように更に構成される。
【００５４】
データ調整器３２は、以下の関係式に従って、既存の推定値Ｘ（αｍ）を調節することにより、信号の周波数領域表現を再推定するように更に構成される。
0006397001.tif 000010
ここで、上付き文字は反復回数を表し、例えばＸ（αｍ）ｉ＋１は、（ｉ＋１）番目の反復における再推定された信号の周波数領域表現を表し、Ｘ（αｍ）ｉは、ｉ番目の反復における既存の推定を表す。データ調整器３２は、パラメータａｋ、ｔｋが収束するまで、フィルタＦ及び信号の周波数領域表現の再推定値Ｘを判定する処理を反復するように構成される。収束基準は、連続する反復の間に推定されたパラメータａｋ、ｔｋにおける絶対変動若しくは相対変動に基づいて、及び／又は、原信号における測定雑音の推定値（例えば、反復の間の調節が、推定された測定雑音を下回るかの判定）に基づいて、選択可能である。
【００５５】
図４を参照して、信号処理装置１２の動作を以下で説明する。
動作において、データスケーラ３０は、データ取得器２２からサンプルｙｎの系列を取得し（１００）、係数提供器２０から係数Ｃｍｎを取得する（１０２）。データスケーラ３０は、当該係数を使用してサンプルをスケーリングすることで（１０４）、上記の式３のように信号の周波数領域表現の推定値を提供する。
【００５６】
データ調整器３２は、データ取得器から原信号における測定雑音の推定値を取得する（１０６）。その後、データ調整器３２は、式５で定義された関係に従ってＸ（αｍ）から取得された推定振幅ａｋを使用して、モデルミスマッチをスケーリングすることで（１１０）、Ｘ（αｍ）の推定値における誤差の推定値ζｍを判定する。
【００５７】
データ調整器３２は、推定測定雑音と上記の式５で定義された誤差の推定値ζｍとを比較する（１１２）。測定雑音がモデルミスマッチに関連する誤差を上回る場合、データ調整器３２は推定値を調整せず、データ提供器３４は推定値Ｘ（αｍ）を提供する（１１４）。
【００５８】
データ調整器３２は、上記の式６で定義されたように、係数Ｃｍｎ及び基底関数ｈｎにより提供される、複素指数関数の表現における誤差に基づいて、「モデルミスマッチ（model mismatch）」の推定値を取得する（１０８）。
【００５９】
データ調整器３２は、例えば畳み込み積分Ｆ＊Ｘを０に向けて調節するフィルタＦの係数を選択することにより、周波数領域推定表現Ｘのエネルギーを減少するフィルタＦを判定する（１１６）ことで、振幅ａｋを取得する。データ調整器３２は、当該フィルタのｚ変換の根と、タイミングｔｋを生じさせるｕｋ（式７を参照）と振幅ａｋとを判定する（１１８）ことで、Ｘを再推定する。これは、付録１において説明され、且つ、Y.Hua及びT.K.Sakar、「Matrix Pencil Method for Estimating Parameters of Exponentially Damped Undamped Sinusoids in Noise」、IEEE Transactions on Acoustics、Speech and Signal Processing、vol.38、814〜824ページ、1990年5月において定義されるような行列ペンシル法を使用して実行されてもよい。
【００６０】
その後、ｔｋ及びａｋの値は、図３を参照して上述したように、収束が試験される（１２２）。値が収束した場合、データ提供器は、信号の周波数領域表現の推定値Ｘを提供する。値が収束していない場合、ステップ１１６、１１８、１２０及び１２２が繰り返される。
【００６１】
信号ｘ（ｔ）の周波数領域推定値Ｘ（αｍ）は、信号領域における信号ｘ（ｔ）を推定するために使用されてもよく、それにより、タイミング推定器１４は送信信号と反射信号との間の時間遅延を判定できる。この推定は、イノベーションパラメータを特定するために行列ペンシル問題の固定値を判定することにより、「消滅」フィルタ法若しくはそれに対する近似を使用することにより、又は、付録１において定義されるProny法の他の何らかの適用により、実行されうる。
【００６２】
理解されるように、これまでに説明した例は、画像が同様の方法で処理されうるような、時間領域信号、空間信号に関連している。パルス列をディラックのデルタ関数のタイミング及び振幅により特徴付けることができるのと同様に、画像は、画素の位置及び強度、又は、画像の境界の位置、振幅及び傾斜により定義されうるため、画像がＦＲＩ信号として扱われてもよい。図３及び図４を参照して説明した信号処理装置１２及び方法は、画像のイノベーションパラメータの判定に同様に適用されてもよい。
【００６３】
２次元画像における領域間の境界の定義に関連する１つの特定の例を図５に示す。
図５は、画像処理装置の一例を示す。図５に更に示されているのは、原画像と、装置関数により及び画素へのサンプリングにより修正された画像の第２のバージョンとである。図５の例において、画像の第２のバージョンの画素値は、原信号のサンプルｙ（ｎ）に対応する。
【００６４】
図５の装置は、図３及び図４を参照して上述した信号処理装置１２を備える。信号処理装置１２は境界判定器と接続され、境界判定器１４４は補間器と接続される。
【００６５】
信号処理装置１２は、画素値を受信し、図４を参照して説明した方法に従って当該画素値を処理して、画像の周波数領域表現を提供するように構成される。図３及び図４において説明した例と同様に、画像の周波数領域表現は、撮像（imaging）の装置関数が特徴付けられた周波数αｍに基づいている。当然、これらは水平成分及び垂直成分を含みうるため、周波数領域表現はＸ（αｍｘ，αｍｙ）で表記することができる。
【００６６】
境界判定器１４４は、画像における局所周波数により値Ｘ（αｍｘ，αｍｙ）をスケーリングすることで、スケーリングされた周波数領域画像データＩ（ｘ，ｙ）を次式のように提供するように構成される。
0006397001.tif 000011
ここで、
0006397001.tif 000012
であり、ｊは、−１の平方根を示し、Ｍ及びＮは、ｘ方向及びｙ方向における画像の大きさをそれぞれ表し、ｐ及びｑは、インデックスである。
【００６７】
境界判定器１４４は、以下により定義される式を解くか又は（例えば、最良のフィッティングの意味で）その近似解を見つけることにより、量Ｉを０に向けて減少させる係数ｂの集合を判定するように更に構成される。
0006397001.tif 000013
ここで、＊は畳み込みを示す。これは、量Ｉを消滅させるフィルタを見つけることに類似し、上述したような及び付録１で説明するような消滅フィルタ又は行列ペンシル法を使用して対処されてもよい。境界判定器１４４は、ｂｐｑが判定されると、ｆ（ｘ，ｙ）＝０である場合の画像内の位置（ｘ，ｙ）として境界を特定するために、以下の連立方程式を解くように構成される。
0006397001.tif 000014
ここで、
0006397001.tif 000015
である。
【００６８】
この連立方程式において、Ｍ及びＮは、原画像の大きさであり、ｊは、−１の平方根を示す。この結果の検証は、H.Pan、T.Blu及びP.L.Dragotti、「Sampling Curves with Finite Rate of Innovation」、in Proc. of Sampling Theory and Application Conference、Singapore、2011年5月において見つけられうる。これは、画像における境界を特定する１つの特定の方法を提供する。しかし、理解されるように、画像がＦＲＩ信号を含む場合、当該画像の構造（境界等）の位置を定義するために画像のイノベーションパラメータを使用でき、このため、推定値Ｘ（αｍｘ，αｍｙ）を使用して境界を復元する他の方法が、例えば画像のイノベーションパラメータの判定に基づいて更に使用されてもよい。
【００６９】
補間器１４６は、境界判定器から境界データを取得し、画像のサンプルｙ（ｎ）に基づいて複数の境界の間を補間するように構成される。その後、補間された画像は、境界から取得された情報と初期画像データ値とを組み込むことで、画像の強調（enhanced）バージョンを提供する。
【００７０】
画像強調を提供するために境界間を補間することに加えて、本開示の実施形態は、画像を合成して解像度を向上させる方法を含む。図６及び図７を参照して、１つのそのような方法を説明する。
【００７１】
図６及び図７に示すように、本方法において、第１の画像８１２が取得され（７０２）、図３、図４及び図５を参照して説明したように、画像の周波数空間表現Ｘ１（αｍｘ，αｍｙ）が推定される。第１の境界が、上記の式９〜１１に基づいて当該第１の画像８１２において特定される（７０４）。
【００７２】
第２の画像８１０が更に取得され（７０６）、当該画像の周波数空間表現Ｘ２（αｍｘ，αｍｙ）が推定され、第２の境界が当該第２の画像において識別される（７０８）。これは第１の画像と同様の方法で実行可能であるが、第１の境界及び第２の境界は同様の方法で取得される必要がなく、同時に取得される必要もないことが理解されるだろう。
【００７３】
第１の境界と第２の境界とを共同位置合わせするために、アフィン変換等の変換が判定される（７１０）。これは、例えば２つの境界間の差に基づくメリット関数を減少させる（例えば最小化する）変換を見つけることにより実行されうる。
【００７４】
当該変換から、画像の画素位置の集合８０６が判定されることで、第１の画像８１２及び第２の画像８１０よりも高い解像度を有する画素位置の集合８０６が提供される（７１４）。その後、当該変換が第１の画像データ８１２に適用されることで、変換された第１の画像データ８０２が提供される。
【００７５】
変換された第１の画像データ８０２及び第２の画像データ８１０は、変換された第１の画像データ８０２及び第２の画像データ８１０から補間することによって合成されて（７１６）、新しい画素位置の集合８０６における合成画像データ８２０が判定される。
【００７６】
２つの画像の共同位置合わせに関連して本方法を説明したが、当然、これはあらゆる数の画像に適用されうる。
【００７７】
本開示の理解を助けるために、図２の信号処理装置１２の機能を機能要素に名目上細分化している。しかし、それらの機能要素は物理的に別個の構造として提供される必要はなく、機能のいくつか又は全ての態様は共通の物理的構造において提供されてもよく、機能のいくつか又は全ての態様は更に細分化されてもよいことが理解されるだろう。１つの特定の例として、信号処理装置１２は単一のプロセッサにより提供されてもよく、データ提供器３４及びデータ取得器２２は、チップにより読み出し可能及び／又は書き込み可能であるレジスタ又はバッファ等のメモリ要素を備えてもよい。例えばデータ提供器３４及びデータ取得器２２は、プロセッサの内蔵キャッシュ等のメモリキャッシュであってもよい。更に、上記機能は、特定用途向け集積回路（ＡＳＩＣ）、フィールドプログラマブルゲートアレイ（ＦＰＧＡ）又はデジタル論理ゲートのアセンブリ等の何らかの適切なプロセッサにより提供されてもよい。いくつかの実施形態において、上記機能はコンピュータプログラムにより提供される。コンピュータプログラムは、信号処理装置の機能を提供するようにプロセッサをプログラムするように動作可能なプログラム命令を含む。このコンピュータプログラムは、一時的又は非一時的な形態で格納可能である。
【００７８】
図１は、本開示の一実施形態が使用されうる方法の一例を紹介している。この特定の例は距離計測に関するものであり、パルスの系列を含む信号の推定に関するものであるためにこの例が選択されている。当業者には理解されるように、あらゆる信号は、単一の点の系列（例えば、ディラックのデルタ関数の系列）により表されうる。従って、本例は本開示を理解する１つの方法として提供されるが、本開示の装置及び方法は他の多くの状況における信号処理に適用可能であることが理解されるだろう。パルス系列の推定は一例にすぎない。装置関数の例には、フィルタのインパルス応答及び光学系の点広がり関数が含まれる。
【００７９】
図１の距離計測の例では、光信号又は音響（例えば、超音波）信号等の種々の信号を使用できる。しかし、本開示との関連で理解されるように、ＡＤＣの時間分解能、並びに送出信号及び反射信号の速度（例えば、光信号の場合の光の速度）によって、測定の精度が制限される。
【００８０】
装置関数は、単位インパルスに対する応答により特徴付けられてもよく、別の方法は、例えば既知の周波数の正弦波を使用して、既知の周波数及び振幅の信号に対する計器の応答を測定することである。検出計器を特徴付ける他の方法は、当業者には理解されるだろう。
【００８１】
上記の説明において、パラメータｎは信号領域におけるサンプルにインデックスを付けるために使用され、信号領域の例におけるシフトを示すために使用される。信号が時間領域信号である場合、パラメータｎは時間遅延に関連し、信号が空間的に変動する信号である場合、パラメータｎは空間変位に関連する。
【００８２】
遅延ｎは、整数のサンプルを含んでもよく、その場合、係数Ｃｍｎは、シフトされた基底関数上への指数関数ｅｘｐ｛αｍ｝の直交射影から取得されてもよい。しかし、いくつかの実施形態において、遅延ｎは整数及び非整数のサンプルを含んでもよい。これらの可能な例において、係数Ｃｍｎは、次式のようなモデルミスマッチを表すメリット関数の減少に基づくフィッティング手順から取得可能である。
0006397001.tif 000016
【００８３】
これらの係数Ｃｍｎを選択する１つの方法は、上記の式２で定義される関係に従うことである。しかし、他の数値解及び解析解が使用されてもよい。例えば次式が使用されてもよい。
0006397001.tif 000017
ここで、
0006397001.tif 000018
は、ｚ＝ｅｘｐ｛αｍ｝において評価された装置関数
0006397001.tif 000019
の信号領域表現の自己相関のｚ変換である。
【００８４】
係数Ｃｍｎは、他の方法で選択されてもよい。例えば、特定の位置において（例えば、選択された許容範囲内で０になる）特定の値を有するようにモデルミスマッチを制約することに基づいて、係数が選択されてもよい。本アプローチの１つの例は、ｔが整数である指数関数ｅｘｐ｛αｍｔ｝をモデルΣＣｍｎｈｎが正確に補間することを保証することである。本アプローチを用いる場合、係数を次式のように定義できる。
0006397001.tif 000020
モデルミスマッチにより判定されるメリット関数を減少させる他の係数が更に使用されてもよい。一般に、それらの係数の例は装置関数の周波数領域表現を含む。
【００８５】
周波数αｍは、それらの周波数における装置関数の大きさＨ（αｍ）の変動を減少させる一方で周波数間の間隔を増加させるように選択されると説明されている。これは、種々の方法で達成されうる。いくつかの実施形態において、信号は、イノベーション率有限（ＦＲＩ）信号という種類に属し、周波数αｍは、この種類の信号のイノベーションパラメータの最尤推定器の性能に基づいて選択される。これを達成する１つの方法は、クラメール・ラオ限界を減少させる（例えば、最小化する）周波数αｍを選択することであり、これは付録１の式１６により定義される。これらの周波数を選択するために、他の方法が使用されてもよい。係数Ｃｍｎの１つ以上の集合を、係数提供器２０に格納できる。更に、例えばデータ調整器がアクセスするために、モデルミスマッチεｍ（ｔ）もメモリに格納されてもよい。
【００８６】
距離計測の例では、原信号は時間領域信号であるため、周波数領域及び信号領域への言及はその文脈で理解される必要がある。他の例は画像処理に関連しており、その場合、信号は空間領域信号であるため、信号領域及び周波数領域への言及はその文脈で理解される必要がある。当業者には理解されるように、周波数領域と信号領域との間のリンクは、フーリエ級数、フーリエ変換、ラプラス変換、ｚ変換又は他の何らかの周波数変換等の変換の選択肢のうちのいずれか１つにより提供されうる。
【００８７】
特に、本開示の信号処理方法及び装置は、ニューロン信号の推定、地震学又は画像処理において使用されてもよい。本開示の例は、各チャネルが図３を参照して説明したような信号処理装置１２を含むマルチチャネルデータ収集装置を提供する。
【００８８】
計器により収集され、且つ、当該計器による信号の収集に関連付けられた装置関数により歪みが生じた信号のサンプル系列から、信号の周波数領域表現を推定する信号処理方法が開示されたことが理解されるだろう。本方法は、本明細書中で説明される種々の方法のうちのいずれか１つで適用されうることが更に理解されるだろう。例えば、本方法は、サンプルの系列を取得する工程と、信号領域における装置関数のシフトされたバージョンによってそれぞれが定義される複数の基底関数を含む、基底関数の集合を、複素指数関数にフィッティングする係数の集合を取得する工程と、サンプルの系列及び係数に基づいて、信号の周波数領域表現を推定する工程とを含みうる。装置関数は、複素指数関数に関連付けられた周波数における、周波数領域における装置関数の特性に基づいている。本方法のこれらの特徴の各々は、本明細書中で説明される方法のうちのいずれか１つ以上で実現され、一般化され又は更に改良されてもよい。
【００８９】
従って、説明される実施形態は例示として理解されるべきであることが理解されるだろう。更なる実施形態が考えられる。いずれか１つの実施形態に関して説明されるあらゆる特徴は、単独で使用されてもよく、説明される他の特徴と組み合わせて使用されてもよく、何らかの他の実施形態又は何らかの他の実施形態の何らかの組み合わせの１つ以上の特徴と組み合わせて使用されてもよいことが理解されるべきである。更に、添付の特許請求の範囲において定義される本発明の範囲から逸脱することなく、上述されなかった均等物及び変形が使用されてもよい。
【００９０】
図面全体を参照すると、概略的な機能ブロック図は本明細書中で説明するシステム及び装置の機能を示すために使用されることが理解されるだろう。しかし、機能はそのように分割される必要がなく、説明され且つ以下でクレームされる以外のハードウェアの何らかの特定の構造を示すと理解されるべきでないことが理解されるだろう。図面に示す要素のうちの１つ以上の要素の機能は、更に細分化されてもよいし、及び／又は本開示の装置全体に分散されてもよい。いくつかの実施形態において、図面に示す１つ以上の要素の機能は、単一の機能ユニットに一体化されてもよい。
【００９１】
いくつかの例において、１つ以上のメモリ要素は、本明細書中で説明される動作を実現するために使用されるデータ及び／又はプログラム命令を格納できる。本開示の実施形態は、本明細書中で説明され、及び／又はクレームされる方法のうちのいずれか１つ以上を実行するようにプロセッサをプログラムするよう動作可能なプログラム命令と、本明細書中で説明され、及び／又はクレームされるデータ処理装置を提供するように動作可能なプログラム命令との少なくともいずれかを含む有形の非一時的な格納媒体を提供する。
【００９２】
本明細書中で概略を示す動作及び装置は、論理ゲートのアセンブリ等の固定論理、あるいはプロセッサにより実行されるソフトウェア及び／又はコンピュータプログラム命令等のプログラム可能論理により提供されてもよい制御器及び／又はプロセッサを使用して実現されてもよい。他の種類のプログラム可能論理は、プログラム可能プロセッサ、プログラム可能デジタル論理（例えば、フィールドプログラマブルゲートアレイ（ＦＰＧＡ）、消去可能プログラム可能読み出し専用メモリ（ＥＰＲＯＭ）、電気的消去可能プログラム可能読み出し専用メモリ（ＥＥＰＲＯＭ）、特定用途向け集積回路ＡＳＩＣ、あるいは他の何らかの種類のデジタル論理、ソフトウェア、コード、電子的命令、フラッシュメモリ、光ディスク、ＣＤ−ＲＯＭ、ＤＶＤ  ＲＯＭ、磁気カード又は光カード、電子的命令を格納するのに適した他の種類の機械読取可能媒体、あるいはそれらの何らかの適切な組み合わせを含む。
【００９３】
［付録１］
以下は、本発明が考察されてもよい方法の詳細な説明である。種々の特徴が必須又は必要であると説明されうるが、これは、例えば対象となっている特定の数学的状況により求められる他の条件に起因して、説明される特定の例に対してのみ当てはまる場合がある。従って、それらの記載は、本発明を何らかの方法で限定するものであると解釈されてはならない。
【００９４】
Starang-Fixの近似：
任意のカーネルを用いるＦＲＩのサンプリング
Jose Antonio Uriguen、Thierry BIu、及びPier Luigi Dragotti
【００９５】
＜要約＞
近年、イノベーション率有限（ＦＲＩ）信号として知られる特定の種類の帯域制限のない信号のサンプリング及び正確な再構成のための複数の新しい方法が開発された。[1]において使用されるサンプリングカーネルは指数関数再生カーネルと呼ばれ、一般化Strang-Fix条件を満たす。これらは、カーネルとそのシフトされたバージョンとの適切な線形結合が多項式又は指数関数を正確に再生することを保証する。
【００９６】
本稿の第１の目的は、雑音を含むサンプルからＦＲＩ信号を推定する場合に正しい指数関数再生カーネルを選択する方法に対する明確な指針を提供することである。次に、サンプリングカーネルを選択できる状況から離れ、あらゆるカーネルを扱うことができるという意味で汎用である新しい方法を開発する。そのために、正確な指数関数再生条件を満たすことは非常に困難な制約であることに注目し、指数関数再生の考えをあらゆるカーネルと共に使用し、且つ、雑音に対するロバスト性の高い再構成方法を開発するために、再生式における制御誤差を許容する。
【００９７】
数値的結果は、本稿の種々の目的を実証し、特に、近似指数関数再生方法が従来の完全復元方法を使用する場合に取得される結果より安定性の高い再構成結果を生じることを示す。
【００９８】
＜索引用語＞
サンプリング、イノベーション率有限、雑音、ＭＯＭＳ、行列ペンシル
【００９９】
＜ＥＤＩＣＳカテゴリ：ＤＳＰ−ＳＡＭＰ＞
I．序論
大部分のデジタル収集システムは、アナログ信号からデジタル信号への変換を含み、サンプリング定理は、連続時間の世界と離散時間の世界との間を橋渡しする。通常、収集処理はFigure 1のようにモデル化され、平滑化関数φ（ｔ）はサンプリングカーネルと呼ばれ、一般に、収集デバイスによる歪みをモデル化する。フィルタリングされた連続時間信号ｙ（ｔ）＝ｘ（ｔ）＊φ（−ｔ／Ｔ）は、ｆs＝１／Ｔのレートで均一にサンプリングされる。この設定に従うと、測定は次式により与えられる。
0006397001.tif 000021
【０１００】
サンプリングの基本的な問題は、サンプルｙｎの集合を使用して元の連続時間波形ｘ（ｔ）を復元することである。信号が帯域制限される場合、シャノンによる解は周知である。最近、特定の種類の帯域制限のない信号をサンプリングして完全に再構成できることが[2],[1],[3]において示されている。そのような信号は、有限の個数の自由パラメータにより完全に説明されるため、イノベーション率有限（ＦＲＩ）信号と呼ばれる。完全再構成は、消滅フィルタ法[4]としても既知であるProny法の変形を使用して達成される。当該フレームワーク内でサンプリング可能な信号は、ディラック等のパルスのストリーム[2],[1],[3],[5]、区分的多項式信号、区分的正弦波信号[6]及び２次元信号の種類[7],[8],[9],[10]を含む。雑音が存在する場合、ＦＲＩ再構成技術は不安定になり、雑音に対する復元力を向上させる方法が[11],[12],[13],[14],[15],[10]において提示されている。
【０１０１】
ＦＲＩ信号を完全再構成するために、ＦＲＩに関するオリジナル論文において最初に提案されたｓｉｎｃ関数及びガウス関数[2]、ＳｏＳ（Sum of Sines）のファミリー等のコンパクトサポートカーネル[3]、並びに多項式再生カーネル及び指数関数再生カーネル[1],[16]等の種々のサンプリングカーネルを使用できることに注目する。これらは全て、雑音を含まない設定において完全再構成が可能であり、雑音が存在する場合は挙動が変化する。従って、当然、性能の低下の原因となる要素を理解し、雑音に対する復元力が最も高いカーネルを設計する方法に関する正確な指針を提供することが試みられる。
【０１０２】
本稿では、指数関数再生カーネルのファミリー[1]を中心として、安定性が最も高いカーネルを設計する方法に関する指針を提供する。指数関数再生カーネルを中心とするのには２つの理由がある。第１に、それらは、雑音を含む測定値を処理する際に適切な特性であるコンパクトサポートを有することができる。第２の更に重要な理由は、Figure 1の設定を使用してＦＲＩサンプリングにおいて従来使用されてきたコンパクトサポートカーネルが指数関数再生カーネルのファミリーの特定の例であることである。これを付録Ｂにおいて証明する。
【０１０３】
我々の貢献は２つに分けられる。最初に、雑音を含む環境においてＦＲＩ信号をサンプリングし、且つ、再構成する場合に最も効果的な指数関数再生カーネルを設計する方法を説明する。ＦＲＩ復元は雑音内のパラメータの集合を推定することに等しいため、このＦＲＩ推定問題のクラメール・ラオ限界を最小化する指数関数再生カーネルのファミリーを見つけることにより、この最適な設計を達成する。第２の貢献において、要求に最も適合する方法でカーネルを作成できる状況から離れ、収集デバイスを制御できない場合について考える。この新しいシナリオにおいて、あらゆる可能なカーネルを扱うことができるという意味で汎用であるＦＲＩサンプリングに対する新しい方法を開発する。入力の正確なパラメータを見つけようとする既存の技術と異なり、雑音が存在する場合も高精度が可能な近似復元を実行する方法を説明する。新しい方法の利点は、あらゆるサンプリングカーネルに適用可能であることである。更に、新しい方法は、完全再構成を理論上達成できるが実際は非常に不安定である多項式スプライン等の、カーネルに対する従来の完全復元フレームワークの代わりに使用できる。この新しい近似フレームワークを完全復元方法の代わりに使用することにより、「不安定な」カーネルを使用する場合でも、雑音を含む設定におけるＦＲＩ信号の信頼性の高い再構成が可能であることを示す。
【０１０４】
本稿の概要は次のようである。第II節において、雑音を含まない例を紹介し、その例において、原型的なＦＲＩ信号であるディラックの系列をサンプリングして完全再構成する。指数関数再生カーネル[1]及び一般化Strang-Fix条件[17]についても説明する。第III節において、収集処理に雑音が存在する場合の更に現実的な設定を扱う。当問題に対してディラックの系列を回復し、且つ、クラメール・ラオ限界を計算する実用的な技術を説明する。第IV節において、雑音に対する復元力が最も高い指数関数再生カーネルのファミリーを設計する。第V節において、近似ＦＲＩフレームワークを提示し、あらゆるカーネルを用いてＦＲＩ信号をサンプリングする基本的な考えを開発する。再構成段階を適切に定義することにより、この新しいフレームワークを用いることで既存のカーネルの再構成精度を向上できることを更に示す点が興味深い。第VI節において、本稿の種々の目的を証明するシミュレーション結果を提示する。最後に、第VII節において結論を示す。
【０１０５】
II．イノベーション率有限信号のサンプリング
本節において、ＦＲＩ信号のサンプリングの理論の概要を説明する。特に、サンプルからディラックのストリームを完全再構成する方法を説明する。なお、更に高度な信号の再構成は常に本例に帰着させることができる。
【０１０６】
Ａ．ディラックのストリームの完全再構成
理解を容易にするために、入力ｘ（ｔ）は、ｔｋ∈［０，τ）及びｋ＝０,..., Ｋ−１である異なる瞬間ｔｋに位置し、振幅ａｋを有するＫ個のディラックのストリームであると考える。信号は全体で２Ｋの自由度を有し、次式として示すことができる。
0006397001.tif 000022
【０１０７】
次に、Figure 1の収集モデルに基づいて、カーネルφ（ｔ）を用いて入力をフィルタリングし、以下の測定値を取得する。
0006397001.tif 000023
ここで、ｎ＝０，１,..., Ｎ−１である。ここで、サンプリング周期Ｔがτ＝ＮＴを満たすと仮定する。更に、φ（ｔ）はコンパクトサポートであり、ｍ＝０,..., Ｐ、且つ、αｍ∈Ｃである適切な係数ｃｍ,ｎに対し次式を満たすと更に仮定する。
0006397001.tif 000024
（３）を満たすカーネルを次数Ｐ＋１の指数関数再生カーネルと呼ぶ。第II−Ｂ節において、これを詳細に説明する。本稿の残りの部分において、α０, λ∈Ｃであるαｍ＝α０＋ｍλの形式のパラメータに解析を限定する。これを選択した理由は、本節において後で明らかになる。
【０１０８】
入力をサンプリングした後、入力を測定値ｙｎの集合から一義的に回復できることを示す必要がある。そのために、最初に、サンプルｙｎと（３）の係数ｃｍ,ｎとを線形結合して、ｍ＝０,..., Ｐに対する新しい測定値を取得する。
0006397001.tif 000025
その後、信号ｘ（ｔ）はディラックのストリーム（１）であるとし、（４）と（２）とを合成することにより、次式を得られる[1]。
0006397001.tif 000026
ここで、
0006397001.tif 000027
である。なお、値ｓｍは信号ｘ（ｔ）の（指数関数）モーメントであり[6]、指数関数
0006397001.tif 000028
により張られる部分空間上にｘ（ｔ）を射影することに等しい。ｍ＝０,..., Ｐであるαｍ＝α０＋ｍλの形式のパラメータを使用することは、値ｓｍが冪和（power sum）の級数の形式を有するために必要であり、これは復元段階において重要な点である。
【０１０９】
ｋ＝０,..., Ｋ−１とする未知数（ｕｋ，ｘｋ）の新しい対は、スペクトル推定の分野[4]ではProny法としても知られる消滅フィルタ法[2],[1],[12]を使用して測定値ｓｍから回復可能である。ｍ＝０,..., Ｋであるｈｍがｚ変換
0006397001.tif 000029
を有するフィルタを示すとすると、すなわち、その根は探される値ｕｋに対応する。その場合、ｈｍは観察された級数ｓｍを消滅させることになる。
0006397001.tif 000030
【０１１０】
更に、位置ｔｋは別個であるとすると、当該フィルタの０は値ｕｋを一意に定義する。恒等式（６）は行列とベクトルの形式で以下のように書くことができる。
Ｓｈ＝０      （７）
これは、テプリッツ行列Ｓがランク不足であることを示す。上記の式を解くことにより、フィルタ係数ｈｍを見つけ、
0006397001.tif 000031
の根を計算することによりｕｋを回復する。ｕｋが与えられる場合、
0006397001.tif 000032
であるため、位置ｔｋが取得される。アルゴリズム１における消滅フィルタ法の主要なステップを詳述する。
【０１１１】
なお、λが純虚数である場合、ｕｋは周期
0006397001.tif 000033
の繰り返しであり、位置ｔｋは、一義的に回復されるために、ｋ＝０,..., Ｋ−１とする
0006397001.tif 000034
を満たす必要がある。最後に、例えば（５）において最初のＫ個の連続する式を解くことにより、重みａｋを判定する。なお、少なくとも未知数と同数の式が存在する場合のみ問題を解くことができ、これは、Ｐ≧２Ｋ−１を意味する。これは、指数関数再生カーネルの次数Ｐ＋１が入力信号ｘ（ｔ）の自由度の数に従って選択される必要があることを示す。
0006397001.tif 000035
【０１１２】
以上で提案した全てのＦＲＩ再構成方法（[2],[1],[12],[3]）は、この同一の方法を使用する。特に、形式（２）のＮ個のサンプルは、ベクトルｙ＝（ｙ０，ｙ１,..., ｙＮ−１）Ｔを形成すると考えられる。その場合、それらは（４）のように係数ｃｍ,ｎと線形結合される。これは、ベクトルｙにｓ＝Ｃｙ＝（ｓ０，ｓ１,..., ｓＰ）Ｔとなる適切な重み行列Ｃを乗算することに等しい。ｓはＰ＋１個の新しい測定値を有するベクトルであり、Ｃは位置（ｍ，ｎ）に係数ｃｍ,ｎを有する（Ｐ＋１）×Ｎ行列である。値ｓは、ｋ＝０,..., Ｋ−１とするイノベーションパラメータ（ａｋ，ｔｋ）を含む冪和の級数の形式で常に表すことができる。その場合、それらは消滅フィルタ法を使用してｓから回復可能である。Ｃの選択はサンプリングカーネルに依存し、雑音を含まない設定では完全再構成が常に達成されるが、本稿において説明される雑音を含む非理想的な状況において、Ｃの役割は非常に重要になる。
【０１１３】
Ｂ．指数関数再生カーネル
指数関数再生カーネルは、そのシフトされたバージョンの線形結合と共に、複素パラメータαｍを用いる形式ｅｘｐ｛αｍｔ｝の関数を再生できるあらゆる関数φ（ｔ）である。これは、適切に選択された係数ｃｍ,ｎ∈Ｃ対して、以下のように数学的に表現できる。
0006397001.tif 000036
ここで、ｍ＝０,..., Ｐであり、αｍ∈Ｃである。（８）が当てはまる指数関数再生カーネルは、いわゆる一般化Strang-Fix条件[17]を満たす（付録Ａを参照）。特に、式（８）は、以下のとき、及び、そのときに限り成立する。
0006397001.tif 000037
ここで、ｍ＝０,..., Ｐ、
0006397001.tif 000038
であり、
0006397001.tif 000039
は、ｓ＝αｍにおけるφ（ｔ）の両側ラプラス変換を表す。
【０１１４】
更に、（８）の係数ｃｍ,ｎは次式により与えられる。
0006397001.tif 000040
ここで、関数
0006397001.tif 000041
は、φ（ｔ）を有する正規双直交集合を形成し、すなわち、
0006397001.tif 000042
が
0006397001.tif 000043
を満たし、
0006397001.tif 000044
である。
【０１１５】
あらゆる指数関数再生カーネルは、
0006397001.tif 000045
[1],[18],[19]として書くことができ、ここで、γ（ｔ）は任意の関数であり、分布であってもよく、
0006397001.tif 000046
はＥスプラインである。フーリエ変換
0006397001.tif 000047
を有する関数βα（ｔ）は、α∈Ｃである一次基数指数スプラインと呼ばれる[18]。そのような関数の時間領域表現は、ｔ∈（０，１）の場合はβα（ｔ）＝ｅαｔであり、それ以外の場合は０である。従って、関数βα（ｔ）は、コンパクトサポートであり、そのシフトされたバージョンの線形結合は、指数関数ｅαｔを再生する。高次の基数指数スプライン（Ｅスプライン）は、一次基数指数スプラインの畳み込みを介して取得可能であり、例えば
0006397001.tif 000048
であり、ここで、
0006397001.tif 000049
は次数Ｐ＋１のＥスプラインである。これはフーリエ領域において以下のように書くこともできる。
0006397001.tif 000050
【０１１６】
最高次のＥスプラインの長さはＰ＋１であり、その正則性はＰと共に増加する（すなわち、それらはＰ−１個の連続する導関数を有する）。これらの関数は、スパン
0006397001.tif 000051
において、あらゆる指数関数を再生できる[1],[18]。ｍ＝０,..., Ｐに対してαｍ＝０である場合は常に、関数
0006397001.tif 000052
はＢスプラインに帰着し、指数関数を再生しないがＰ次までの多項式を再生することにも注目する。
【０１１７】
本稿では、実数値のサンプリングカーネルが扱われる。従って、
0006397001.tif 000053
及び
0006397001.tif 000054
は実数値である必要がある。この第２の条件（すなわち、
0006397001.tif 000055
が実数値であること）は、指数関数αｍが実数であるか又は複素共役対である場合に満たされることに注目する。更に、指数関数パラメータをｍ＝０,..., Ｐであるαｍ＝α０＋ｍλの形式に制限する。これらの条件は、λ及びα０がいずれも実数であるか、あるいはλ＝ｊλ０が純虚数であることを示す。後者の場合、Ｉｍ｛α０｝＝−Ｐλ０／２である。
【０１１８】
III．雑音が存在する場合のＦＲＩ信号のサンプリング
「雑音」、又は更に一般的にはモデルミスマッチが、データ収集において存在し、第II−Ａ節において提示した解決策を非現実的にする。従って、次式で表す測定値を使用できると仮定する。
0006397001.tif 000056
ここで、ｎ＝０,..., Ｎ−１であり、εｎは、平均０及び標準偏差σのi.i.d.ガウス確率変数である。サンプルが雑音により破損する場合、式（４）の測定値ｓｍの集合が変化し、完全再構成がもはやできなくなる。その場合、ｍ＝０,..., Ｐとして次式で示す雑音を含むモーメントを有する。
0006397001.tif 000057
ここで、ｋ＝０,..., Ｋ−１であり、
0006397001.tif 000058
である。
【０１１９】
次に、雑音を含むサンプルからＦＲＩ信号を復元するために本稿において使用したアルゴリズムのうちのいくつかを再確認する。また、雑音が存在する場合にどのＦＲＩ信号を回復できるかの精度を解析するために、クラメール・ラオの下限（ＣＲＢ）を導入する。これは、不偏推定量に適用する平均二乗誤差（ＭＳＥ）に対する下限である[20]。
【０１２０】
Ａ．雑音を含むＦＲＩを再構成するための実用的な方法
第II−Ａ節で説明したように、ＦＲＩ信号の回復処理は、新しい測定値の集合の計算（４）に基づいている。サンプル（２）が雑音により破損する場合、式（４）の測定値ｓｍは（１３）になる。その結果、
0006397001.tif 000059
であり、Ｂが（１３）からのエントリｂｍを有するテプリッツ行列であるため、（７）は満たされない。[12]において使用される手順の背景にある考えは、
制約||ｈ||２＝１の下で
0006397001.tif 000060
を最小化する解を見つけることである。これは、特異値分解（ＳＶＤ）を使用して解決できる従来の総最小二乗（ＴＬＳ）問題である。アルゴリズムは、Cadzowの反復アルゴリズム[12],[21]を使用することにより、ＴＬＳを適用前に
0006397001.tif 000061
の雑音を除去することで更に改善されてもよい。同様の精度を達成する他の方法が存在する。１つのそのような方法は、行列ペンシル問題[22],[23]を解くことに基づき、[11]においてＦＲＩに対して導入された。これは、[3],[16]等の他のＦＲＩに関する出版物において使用されており、本稿のシミュレーションにおいても使用される。
【０１２１】
乱れが白色である場合、すなわち、雑音Ｂの共分散行列が、αを定数係数とし、且つ、Ｉを単位行列として、ＲＢ＝Ｅ｛ＢＨＢ｝＝αＩを満たす場合、上記の全ての方法は効果的に動作する。そのような場合、雑音は、
0006397001.tif 000062
のＳＶＤ分解の主値に対して同様の影響を有する[24],[11]。これは、これまでに再確認した回復手順が最も適切である理由を説明する。しかし、多くのＦＲＩカーネルに対して、白色ガウス雑音の仮定は成立しない。その結果、テプリッツ雑音行列Ｂのエントリは不均一な変動を有するか又は相関性を有するようにさえなる場合がある。ＳＶＤが適切に動作するために、雑音を「事前に白色化」することが必要である。これは、スペクトル推定の分野において様々な著者により（例えば、[24]においてDE Moorにより）提案される周知の方法である。本稿のシミュレーションにおいて、ＲＡ＝Ｅ｛ＡＨＡ｝＝Ｉ、且つ、Ａ＝ＢＷであるように、重み行列Ｗ＝ＲＢ１／２[25]を使用することからなる一般的な方法を使用する。従って、白色雑音によって特徴付けられた
0006397001.tif 000063
のＳＶＤ分解を実行する。
【０１２２】
Ｂ．性能の測定
ＦＲＩ信号は、各自のイノベーションパラメータにより完全に特徴付けられる。例えばＫ個のディラックのストリームは、位置ｔｋ及び振幅ａｋから完全に判定可能である。従って、ベクトルΘ＝（ｔ０,..., ｔＫ−１，ａ０,..., ａＫ−１）Ｔを導入し、（１２）により与えられる雑音を含むＮ個のサンプル
0006397001.tif 000064
のベクトルからΘを推定することを目的とする。簡略にするために、サンプリング周期はＴ＝１であると仮定する。この推定問題のＣＲＢを判定する方法は、εｎが共分散行列Ｒ＝Ｅ｛ｅｅＨ｝を有する平均０のガウス雑音であり、ｅが値εｎを有する長さＮのベクトルであると仮定して、[12]において与えられた。この設定において、未知パラメータ
0006397001.tif 000065
のあらゆる不偏推定値は、次式により下限を示される共分散行列を有する。
0006397001.tif 000066
ここで、行列Φｙは（１５）により与えられる（１）。
0006397001.tif 000067
（（注釈１）行列は、Θにおける各パラメータに関して
0006397001.tif 000068
の導関数を計算することにより取得可能である。すなわち、｜の左側のΦｙの列は
0006397001.tif 000069
であり、｜の右側のΦｙの列は
0006397001.tif 000070
である。）
【０１２３】
以上は、雑音を含むサンプル
0006397001.tif 000071
を与えられた場合に種々のＦＲＩ復元技術[12],[11]の性能を測定する正しい方法であるが、ＦＲＩ再構成方法が実際はモーメント
0006397001.tif 000072
のシーケンスに対して作用することに注目する。従って、測定値
0006397001.tif 000073
（サンプル数より少ない場合がある）に関連するＣＲＢが実際に達成可能な最良の性能を示すため、このＣＲＢを見つけることが有用である。この新しい例において、限界は[26]により与えられる。
0006397001.tif 000074
ここで、ΦＳは以下の形式をとる。
0006397001.tif 000075
また、Ｒｂ＝Ｅ｛ｂｂＨ｝である。ここで（・）Ｈはエルミート転置であり、ｂは雑音を含むＰ＋１個の値ｂｍのベクトルである。
【０１２４】
更に、サンプルに加えられる雑音ｅｎは、分散σ２を有する加法性白色ガウス雑音であると仮定したため、フィルタリングされた雑音の共分散行列は、Ｒｂ＝Ｅ｛ｂｂＨ｝＝σ２ＣＣＨである。そのため、公式（１４）及び（１６）の間には、行列Ｃを介して表されうる直接的な関係が存在する。これを確認するために、
0006397001.tif 000076
であること、及びΦＳ＝ＣΦｙであることに注目する。更に、雑音ｅｎがＡＷＧＮであるため、次式が成立することが既知である。
0006397001.tif 000077
そのため、次式が成立する。
0006397001.tif 000078
ここで、
0006397001.tif 000079
はＣの左側擬似逆行列である。モーメントの数Ｐ＋１がサンプル数Ｎと等しい場合、双方の公式は等価である。これは、
0006397001.tif 000080
から容易にわかる。係数の行列Ｃが正方行列であり可逆である場合、各公式を使用して予測されるＣＲＢは一致する。サンプル自体に基づいて推定を向上させるサンプルの線形結合が存在しないため、これは直感的に最適な構成である。Ｐ＋１＜Ｎの場合にモーメントを介してイノベーションパラメータを回復することは、次に最適な構成である。
【０１２５】
Ｃの条件が十分に適切である場合、ＦＲＩアルゴリズムは限界（１６）に到達することを実験的に確認した。従って、次の目的は、Ｐのあらゆる選択に対して（１６）を最小化する適切な条件のＣを生じるカーネルを設計することである。（２）
（（注釈２）本節の範囲を超える過剰に多くの技術的詳細を説明しなくても、条件Ｐ−１＝Ｎがブロック型のサンプリングのみに適用され、例えばＮ個のサンプルを使用して周期信号をサンプリングする場合のみに適用されることがわかる。この条件は、サンプル数が無限である無限長信号には適用できず、順次再構成アルゴリズムは、サンプル数が変化する可能性のあるブロックに対して作用する。）
【０１２６】
IV．最適な指数関数再生カーネル
上述したように、指数関数再生カーネルはあらゆる関数
0006397001.tif 000081
[1],[18]であり、ここで、γ（ｔ）は任意であり、
0006397001.tif 000082
はＥスプラインである。本節において、これらのカーネルを用いるＦＲＩ復元技術を可能な限り安定化するために、ｍ＝０,..., Ｐである指数関数パラメータαｍ＝α＋ｊ（ｍ―Ｐ／２）λ０及び関数γ（ｔ）を選択する方法のルールを見つけたい。最適なパラメータを探すことは、係数の行列Ｃの最適化につながる。従って、最初に、Ｃが満たす必要のある特性を判定し、その後、正しいＣをもたらすカーネルを設計する。
【０１２７】
Ａ．Ｃの正しい選択
ＦＲＩ再構成段階の第１のステップは、サンプルｙのベクトルをモーメントｓ＝Ｃｙのベクトルに変換することである。従って、最初の目的は、適切な条件のＣを取得することである。なお、本ステップの妥当性は[3]において既に強調されている。更に、行列Ｃは、ｎ＝０,..., Ｎ−１及びｍ＝０,..., Ｐである場合に、位置（ｍ,ｎ）の要素ｃｍ,ｎ＝ｃｍ，０ｅｘｐ｛αｍｎ｝によって次式のように構成される。
0006397001.tif 000083
ここで、Ｄは対角行列であり、Ｖはヴァンデルモンド行列である。従って、安定したＣを有するために、Ｄの対角要素の絶対値が同一であること、すなわち｜ｃｍ，０｜＝１であることが求められる。更に、Ｖの要素が次式のように純虚数であることが求められる。
ｍ＝０,..., Ｐに対して
0006397001.tif 000084
すなわち、α＝０及びλ０＝２π／Ｌ。      (18)
【０１２８】
明らかに、純虚数αｍはヴァンデルモンド行列Ｖの条件を改善する[27]。従って、（１８）において最良のＬを見つけるという問題のみが残る。Ｃの条件が適切である場合にＦＲＩアルゴリズムがＣＲＢ（１６）に達することができることが実験からわかっているため、Ｌを決定する１つの方法は、単一のディラックの位置に対して（１６）を最小化する値を選択することである。Ｐ及びＬの種々の選択に対して、Figure 2に示すように、Ｌ＝Ｐ＋１の場合に最小化が常に達成されることがわかる。
0006397001.tif 000085
【０１２９】
この選択は、指数関数が単位円全体を範囲に含むことを保証し、これは冪和の級数のパラメータを復元する場合に最良の構成であることが既知であるため[28]、このことはある程度まで当然である。最後に、（１６）を最小化することに加えて、Ｐ＋１＝Ｎ及びＬ＝Ｐ＋１とする場合、（１６）におけるモーメントに基づくＣＲＢが（１４）におけるサンプルに基づく限界と一致することが更に保証され、最良の可能な性能が達成される。この状況において、行列Ｃは正方のユニタリ行列になる。これは、条件数が１であるため、最も安定する数値変換である。
【０１３０】
要約すると、最良の指数関数再生カーネルは、ｍ＝０,..., Ｐとする場合、
0006397001.tif 000086
の形式の指数を用いて指数関数を再生し、且つ、｜ｃｍ，０｜＝１を有する必要がある。最後に、カーネルの次数は可能な限りＰ＋１＝Ｎであるべきである。次節において、そのようなカーネルを取得する方法を示す。
【０１３１】
Ｂ．指数関数ＭＯＭＳ
前節の解析を用いて、最高次及び最小サポート（ｅＭＯＭＳ）の最適な指数関数再生カーネルを設計する。ｍ＝０,..., Ｐに対して｜ｃｍ，０｜＝１であること、及び、指数関数パラメータは以下の形式であることを必要とする。
0006397001.tif 000087
しかし、以下の方法はあらゆる値のＰに対して使用できるため、Ｐに関する仮定を行わない。
【０１３２】
あらゆる指数関数再生カーネルφ（ｔ）を
0006397001.tif 000088
と書くことができることを考慮することにより、｜ｃｍ，０｜＝１が満たされるようにγ（ｔ）を設計する。（３）を使用することで次式が得られることに注目する。
0006397001.tif 000089
その結果、
0006397001.tif 000090
となり、ここで、（ａ）はポアソンの和公式（３）から得られ、（ｂ）は一般化Strang-Fix条件（９）の適用から得られる。
（（注釈３）ポアソンの和公式：
0006397001.tif 000091
）
従って、あらゆる指数関数再生カーネルに対して、
0006397001.tif 000092
であることがわかる。｜ｃｍ，０｜＝１とすることは、
0006397001.tif 000093
を必要とすることと等価であることが認識される。最後に、
0006397001.tif 000094
であることを使用し、且つ、αｍ＝ｊωｍにおいてラプラス変換を評価することにより、
0006397001.tif 000095
に対する以下の条件に到達する。
0006397001.tif 000096
この場合、各関数のフーリエ変換を扱うことになる。
【０１３３】
（２０）を満たす全ての許容できるカーネルの中で、最短サポートＰ＋１を有するカーネルに関心を有する。従って、元のＥスプライン
0006397001.tif 000097
の種々の導関数の線形結合により与えられるカーネルについて考える。すなわち、
0006397001.tif 000098
ここで、
0006397001.tif 000099
は
0006397001.tif 000100
の第ｌ番目の導関数であり、
0006397001.tif 000101
であり、且つ、ｄｌは係数の集合である。これは、γ（ｔ）がディラックのデルタとその次数Ｐまでの導関数との線形結合であると述べるのと同様である[19]。これらのカーネルは、依然として指数関数ｅｘｐ｛αｍｔ｝を再生でき、[29]において紹介された最高次最小サポート（ＭＯＭＳ）カーネルの変形である。そのため、それらを指数関数ＭＯＭＳ（又はｅＭＯＭＳ）と呼ぶ。それらは、[30]において提示された一般化Ｅスプラインの限界ファミリーの特定の例でもある。この公式の利点は２つある。第１に、変更されたカーネルφ（ｔ）は、
0006397001.tif 000102
と同様に最小サポートＰ＋１のカーネルである。第２に、｜ｃｍ，０｜＝１を達成するためには、制約（２０）を満たす係数ｄｌを見つけるだけでよい。次式により与えられる、（２１）のフーリエ変換
0006397001.tif 000103
を使用する場合、結果として得られる多項式
0006397001.tif 000104
が、ｍ＝０，１,..., Ｐに対して点
0006397001.tif 000105
の集合を補間するように係数ｄｌを選択することにより、（２０）を満たしうることがわかる。
【０１３４】
ｃｍ，０が全てのｍに対して絶対値１を有するということを満たすカーネルを設計すると、位相の曖昧さが残る。位相を決定することは、更なる自由度を導入して、（２１）におけるＥスプラインに対して時間シフトΔを使用することと同等である。指数関数ＭＯＭＳが｜ｃｍ，０｜＝１であり、且つ、パラメータ（１９）が連続時間関数であるためには、ｍ＝０，１,..., Ｐに対して、Δが１以上、且つ、Ｐ以下の整数である場合に、｜ｃｍ，０｜＝｜ｃｍ，０｜ｅｘｐ｛ｊωｍΔ｝であることを示すことが可能である。
【０１３５】
Figure 3には、上述した手順を実現することにより取得されるカーネルのうちのいくつかを提示している。興味深いことには、付録Ｂに示すように、これらの特定の関数は、ディリクレ（Dirichlet）カーネルの１周期と常に等しい。また、Ｐ＋１＝Ｎの場合、この指数関数再生カーネルのファミリーを使用して派生するシナリオは、入力、すなわち、等価的にはサンプリングカーネルを周期化する場合に[2]の元のＦＲＩ公式に収束することを指摘する。
0006397001.tif 000106
【０１３６】
V．イノベーション率有限信号の汎用サンプリング
前節において、雑音を含むＦＲＩサンプリングに対して最適な指数関数再生カーネルを設計する方法を示した。しかし、多くの実際の状況において、サンプリングカーネルφ（ｔ）を選択する自由はない。
【０１３７】
ＦＲＩ設定における重要な点は、φ（ｔ）が指数関数を再生する能力であり、これにより、信号再構成問題をスペクトル線推定理論におけるProny法に帰着できる。本節では、この条件を緩め、指数関数再生特性（３）が必ずしも成立しなくてもよいあらゆる関数φ（ｔ）について考える。それらの関数に対して、指数関数の再生が正確なものではなく近似であるように係数ｃｍ,ｎを見つけることは依然として可能である。この近似再生及び対応する係数ｃｍ,ｎを使用することで、それらのカーネルを使用して取得されたサンプルからＦＲＩ信号を回復することを提案する。
【０１３８】
この新しいアプローチは、複数の利点を有する。第１に、あらゆるカーネルφ（ｔ）と共に使用できるという点で汎用である。実際、以下の節において示すように、この新しい定式化はカーネルの正確な知識すら必要としない。第２の、この新しい方法を用いるＦＲＩ信号の再構成は正確なものにはならないが、多くの場合、適切な反復アルゴリズムにより再構成誤差を任意に小さくできることを示す。最後に、ＦＲＩに関する文献において提案されるいくつかの不安定なカーネルの雑音に対する復元力を増加するために本方法を使用できる。例えば多項式スプライン又はガウス関数等のカーネルは、条件が非常に悪い再構成手順の原因となる。指数関数の近似再生に基づいて適切に選択された係数ｃｍ,ｎから形成された行列に元のＣを置換することにより、同一のカーネルを用いて安定性がはるかに高い再構成が達成されることを示す。
【０１３９】
Ａ．指数関数の近似再構成
指数関数ｅαｎを近似するために、関数φ（ｔ）及びその整数シフトを使用したいと仮定する。特に、次式のように係数ｃｎを見つけたい。
0006397001.tif 000107
【０１４０】
この近似は、φ（ｔ）が一般化Strang-Fix条件（９）を満たす場合のみ正確である。他の何らかの関数の場合、（２２）に最良にフィッティングする係数ｃｎを見つけることが特に有用である。そのような係数を見つけるために、ｃｎ＝ｃ０ｅαｎを直接使用し（４）、１周期関数
0006397001.tif 000108
を使用する。
（（注釈４）正確な指数関数再生係数はｃｎ＝ｃ０ｅａｎを常に満たす。以下、本節を通して導出する異なる近似係数の集合も、同一の形式を有すると予想する。）
【０１４１】
この場合、φ（ｔ）の整数シフトを用いて指数関数ｅαｔを近似することは、定数１によりｇα（ｔ）を近似することに変換可能であることがわかる。これは、ｃｎ＝ｃ０ｅαｎを減算し、且つ、ｅαｔを左側に移動させることにより、（２３）の右側の形式に（２２）を書き直すことができるためである。
【０１４２】
ポアソンの和公式の結果、ｇα（ｔ）のフーリエ級数展開が次式により得られることがわかる。
0006397001.tif 000109
また、本稿の近似問題が次式に帰着することがわかる。
0006397001.tif 000110
【０１４３】
これは、一般化Strang-Fix条件（９）と指数関数の近似との間の関係を更に深く示している。φ（ｔ）が一般化Strang-Fix条件（９）を満たす場合、
0006397001.tif 000111
に対して
0006397001.tif 000112
であり、（２４）は、
0006397001.tif 000113
の場合に正確に成立する。それ以外の場合、
0006397001.tif 000114
に対して項
0006397001.tif 000115
は０にならず、
0006397001.tif 000116
になる係数ｃ０のみを見つけられる。しかし、値
0006397001.tif 000117
が０に近いほど、（２２）における近似が改善する。
【０１４４】
一般に、φ（ｔ）は任意の関数であってもよく、（２２）が成立するために異なる係数ｃｎの集合を見つけることができる。使用する係数に関係なく、ｇα（ｔ）のフーリエ級数展開を使用することにより、本開示の近似の精度を判定できる。実際、係数ｃｎ＝ｃ０ｅαｎを用いる関数
0006397001.tif 000118
によるｆ（ｔ）＝ｅαｔの近似の誤差は、次式に等しい。
0006397001.tif 000119
なお、φ（ｔ）のラプラス変換が十分速く減衰する場合、フーリエ級数展開の殆どの項は誤差に対する正確な限界を有するために必要ない。
【０１４５】
係数ｃｎ＝ｃ０ｅαｎの自然な選択は、最小二乗近似により与えられる係数である。ｆ（ｔ）は二乗可積分でないが、φ（ｔ−ｎ）により張られる部分空間へのｆ（ｔ）の正射影を計算することにより、係数を依然として取得できる[31]。それらは次式の形式をとる。
0006397001.tif 000120
ここで、
0006397001.tif 000121
は、ｚ＝ｅαにおいて評価された
0006397001.tif 000122
のｚ変換である。
【０１４６】
最小二乗近似は、φ（ｔ）の正確な知識を必要とするという欠点を有する。しかし、上述のように、φ（ｔ）のラプラス変換が十分速く減衰する場合、項
0006397001.tif 000123
は
0006397001.tif 000124
の場合に０に近いと仮定できる。この場合、（２５）における誤差は
0006397001.tif 000125
を選択することにより容易に最小化されることがわかる。この第２の種類の近似を定数最小二乗（constant least-squares）と呼ぶ。その簡潔さに加えて、
0006397001.tif 000126
を選択することの第２の利点は、αにおけるφ（ｔ）のラプラス変換の知識のみが必要であることである。ｍ＝０,..., Ｐである指数関数ｅαｍｔの近似再生を必要とするＦＲＩ設定の場合、この簡略化された公式では、ｍ＝０,..., Ｐであるαｍにおけるφ（ｔ）のラプラス変換の知識のみが必要である。
【０１４７】
最後に、第３の興味深い係数選択方法は、ｓ（ｔ）が
0006397001.tif 000127
においてｆ（ｔ）を正確に補間することを保証する係数を選択することである[32],[18]。これらの係数は以下の通りである。
0006397001.tif 000128
なお、補間係数を使用するためには、整数の瞬間におけるφ（ｔ）に関する情報のみが必要である。以上の結果を表１に要約する。
0006397001.tif 000129
【０１４８】
経験上、近似するための指数関数ｅαｔが、純虚数でない指数αを有する場合、最小二乗係数は、より少ない誤差（２５）を提供する。それ以外の場合、定数最小二乗係数も同様に適切である。整数の時点におけるカーネルの値を与えられた場合、補間係数は計算が非常に容易である。しかし、それらは低い近似品質を常に提供する。従って、本稿の残りの部分では、定数最小二乗係数を使用する。
【０１４９】
上記の解析の一例をFigure 4に示す。ここでは、線形スプラインの線形結合を使用して指数関数を近似したい。まず、当該スプラインは次数０及び１の多項式を正確に再生する。これをFigure 4の（ａ）及び（ｂ）に示す。同一の関数を用いて、ｍ＝３,..., ０に対する４つの複素指数関数
0006397001.tif 000130
を定数最小二乗係数
0006397001.tif 000131
で近似再生する。それらの実数部の近似をFigure 4の（ｃ）〜（ｆ）に示す。いくつかの指数関数は他の指数関数より適切に近似されることがわかる。本例において、それらは低周波数を有するものである。更に高次のスプラインを使用した場合、近似品質は向上するだろう。しかし、例示するために線形スプラインを選択している。また、近似可能な指数関数の数は、スプラインの次数に依存しない。
0006397001.tif 000132
【０１５０】
Ｂ．ＦＲＩの復元の近似
ディラックのストリーム
0006397001.tif 000133
及び以下のサンプルについて再度考える。
0006397001.tif 000134
【０１５１】
サンプル（２６）からディラックの位置及び振幅を回復したいが、本例ではサンプリングカーネルについて何も仮定しない。ｍ＝０,..., Ｐ、αｍ＝α０＋ｍλ且つα０，λ∈Ｃである指数関数ｅｘｐ｛αｍｔ｝を近似するφ（ｔ）に対する適切な係数を探す。前節から、
0006397001.tif 000135
であるｃｍ,ｎ＝ｃｍ，０ｅｘｐ｛αｍｎ｝を選択する場合、適切な近似が達成されることがわかっている。従って、αｍ、ｍ＝０,..., Ｐにおけるφ（ｔ）のラプラス変換がわかるだけでよい。また、Ｐはφ（ｔ）のサポートに関連付けられる必要がないが、Ｐ≧２Ｋ−１を条件とするあらゆる値を使用できる。
【０１５２】
イノベーションパラメータ（ｔｋ，ａｋ）を回復するために、指数関数の正確な再生の場合と同一の方法を使用するが、本例において、モーメントは次式であることがわかっている。
0006397001.tif 000136
ここで、
0006397001.tif 000137
である。（２５）の近似誤差εｍ（ｔ）によるモデルミスマッチが存在し、これはζｍに等しい。これを雑音として扱い、第III節の方法を使用して信号のパラメータを回復する。モデルミスマッチは、係数ｃｍ,ｎ、値αｍ及びＰ、並びにカーネルφ（ｔ）により決定される近似の品質に依存する。他の形態の雑音と比較してζｍが無視できるほど小さい場合、手順は十分に適切である。しかし、雑音を殆ど含まない設定において、ディラックの推定は、アルゴリズム２の欄に示す反復アルゴリズムを使用して改善可能である。本アルゴリズムの基本的な考えは、ディラックの位置の推定値を与えられた場合、ζｍの近似を計算でき、それを使用してモーメントｓｍの計算を改善できることである。
0006397001.tif 000138
【０１５３】
Ｃ．パラメータαｍを選択する方法
第IV節において、サンプリングカーネルの設計を完全に制御できる場合、サンプルと同数のモーメントを使用する必要があり、すなわちＰ＋１＝Ｎであり、指数関数パラメータは純虚数であり、且つ、
0006397001.tif 000139
の形式であり、係数ｃｍ,ｎはｍ＝０,..., Ｐに対して｜ｃｍ，０｜＝１である必要があると判定した。しかし、近似されたＦＲＩのシナリオでは、サンプリングカーネルは固定され、モーメントの数Ｐ＋１及び値αｍのみを選択でき、｜ｃｍ，０｜＝１とすることができない。
【０１５４】
このことは、αｍの選択におけるトレードオフにつながる。一方で、αｍは純虚数であること、すなわちαｍ＝ｊωｍであること及び単位円全体に及ぶことが求められる。他方で、値
0006397001.tif 000140
は可能な限り１に近い絶対数を有することが求められる。
0006397001.tif 000141
は通常はローパスフィルタであるため、ｍ＝０,..., Ｐの場合の条件
0006397001.tif 000142
は全てのωｍが０に非常に近い場合に満たされ、これは単位円に及ぶこととは対照的である。従って、次式の形式となるように指数関数パラメータを選択し、上記のトレードオフを最適化するＰ及びＬを判定する。
0006397001.tif 000143
【０１５５】
この場合も、単一のディラックの位置を回復する場合にＣＲＢ（１６）を最小化するＰ及びＬの値を選択するという基準に従う。解析された殆どの例において、最良のＰは、通常はサンプリングカーネルφ（ｔ）のサポート以上であり、Ｌは、Ｐ＋１≦Ｌ≦４（Ｐ＋１）の範囲内であることがわかった。φ（ｔ）が６次のＢスプラインである場合のＬの選択の一例をFigure 5に示す。本例では、最初にＬに対する値を判定し、パラメータ（２８）を計算する。その後、ｍ＝０,..., Ｐに対して係数
0006397001.tif 000144
を計算する。ここで、
0006397001.tif 000145
はφ（ｔ）のフーリエ変換である。最後に、行列Ｃを作成し、単一のディラックについて（１６）を計算する。最小値は常にＬ＝１．５（Ｐ＋１）の近辺である。
0006397001.tif 000146
【０１５６】
VI．シミュレーション
次に、本稿の主要な貢献を実証するためのシミュレーション結果を提示する。特に、サンプリングカーネルが多項式スプラインである場合の、第IV−Ｂ節で紹介したｅＭＯＭＳカーネルの性能及び第V節で紹介した近似ＦＲＩ復元方法の性能を示す。
【０１５７】
Ａ．実験の設定
Ｋ個のディラックから成る系列を入力として有するため、ｎ＝０,..., Ｎ−１に対して
0006397001.tif 000147
を直接計算することにより、Figure 1の方法に従ってＮ個のサンプルを取得する。
0006397001.tif 000148
として定義される目標信号対雑音比に従って選択される分散σ２の加法的白色ガウス雑音を用いてサンプルを破損させる。最後に、雑音を含むＰ＋１個のモーメントを計算し、[11]と同様に行列ペンシル法を使用して、入力のイノベーションパラメータ
0006397001.tif 000149
を回復する。
【０１５８】
時間位置の推定における誤差は回復が最も困難なパラメータであるため、そのような誤差に主に関心がある。各ディラックについて、この誤差の標準偏差を次式に示す。
0006397001.tif 000150
ここで、
0006397001.tif 000151
は反復ｉにおける推定される時間位置であり、Ｉは反復の総数である。固定された信号対雑音比の範囲に対して（２９）を計算し、ＳＮＲ毎にＩ＝１０００の雑音の実現を使用して結果を平均する。性能（２９）と、第III節の２つの異なるクラメール・ラオ限界（ＣＲＢ）、すなわちサンプルに基づくＣＲＢ（１４）及びモーメントに基づくＣＲＢ（１６）により予測した分散の平方根とを比較する。
【０１５９】
Ｂ．指数関数ＭＯＭＳ
Figure 6の（ａ）及び（ｂ）において、標準的なＥスプライン及び第IV−Ｂ節の指数関数ＭＯＭＳカーネルを使用して、Ｎ＝３１個のサンプルからＫ＝２個のディラックを回復する場合のシミュレーション結果を示す。前者は、ｍ＝０,..., Ｐに対して純虚数指数
0006397001.tif 000152
により特徴付けられる。
【０１６０】
考察されるあらゆる次数Ｐ＋１の場合において、ｅＭＯＭＳは、Ｅスプラインより性能が優れていることがわかる。更に、ｅＭＯＭＳは、モーメントに基づくＣＲＢ（赤色で示され、凡例ではｓ−ＣＲＢと表記）を常に達成する。この限界は、Ｐ＋１の値が増加するにつれてサンプルに基づくＣＲＢ（黒色で示され、凡例ではｙ−ＣＲＢと表記）に近付き、Ｐ＋１＝Ｎの場合、予想通りＣＲＢと一致する。ｅＭＯＭＳの安定性を更に示すために、Figure 6の（ｃ）に、τ＝ＮＴ＝１にわたりランダムな間隔で位置し、且つ、任意の振幅を有するＫ＝２０個のディラックの回復を示す。信号対雑音比は１５ｄＢであり、Ｎ＝６１個のサンプルを使用し、Ｐ＋１＝Ｎである。
0006397001.tif 000153
【０１６１】
Ｃ．ＦＲＩの近似
次に、次数Ｐ＋１のＢスプラインカーネルを使用してサンプルを取得し、
0006397001.tif 000154
及びｍ＝０,..., Ｐ２である指数関数ｅｘｐ｛αｍｔ｝を近似する場合、の第V−Ｂ節の近似回復手順の結果を示す。定数最小二乗公式
0006397001.tif 000155
を使用して係数を取得し、
0006397001.tif 000156
はＢスプラインのラプラス変換を表す。本例では、多項式の正確な再生に基づく従来の復元方法と性能を比較する。
【０１６２】
Figure 7の（ａ）及び（ｂ）において、ＳＮＲ＝２０ｄＢを使用し、雑音を含むＮ＝３１個のサンプルからの、τ＝ＮＴにわたってランダムに位置するＫ＝６個のディラックの推定を比較する。（ａ）では、次数Ｐ＋１＝２６のＢスプラインを用いてサンプリング及び再構成している。全てのディラックを見つけられるわけではないことがわかる。（ｂ）では、次数Ｐ＋１＝６のＢスプラインを用いてサンプリングしているが、近似再構成の精度を向上させるためにＰ２＋１＝２６個の指数関数モーメントを作成している。６個のディラックの位置の推定の平均誤差は、（ａ）において適切に復元されたディラックの場合は約１０−２であり、（ｂ）の近似回復の場合は１０−３である。
0006397001.tif 000157
【０１６３】
近似方法を使用して、次数Ｐ＋１＝６のＢスプラインカーネルにより取得されるＮ＝３１個のサンプルからＫ＝２個のディラックを回復する場合の、更なる結果を示す。カーネルの次数を固定する場合でも、より多くのＰ２＋１個のモーメントを再構成することで性能を改善できる。Figure 8の（ａ）〜（ｄ）は、ｍ＝０,..., Ｐ２、且つ、Ｌ＝３／２（Ｐ２＋１）であるパラメータ
0006397001.tif 000158
の場合の結果を示す。モーメントの数Ｐ２＋１が増加するほど性能は改善し、最終的にサンプルに基づくＣＲＢに到達する。
0006397001.tif 000159
【０１６４】
VII．結論
本稿では、雑音の存在下におけるＦＲＩ再構成問題を考察してきた。最初に、雑音を含まない状況における既存の結果、及び最新の文献における最も効果的な雑音の処理を再検討した。次に、前述の処理に対するロバストな代替方法を検討した。
【０１６５】
更に詳細には、本稿の貢献は２つに分けられる。安定性が最も高い信号再構成をもたらすという意味で、最適な指数関数再生カーネルを設計する方法を見つけ出した。更に、サンプリングカーネルを完全に制御できる理想的な状況から離れ、任意の収集デバイスを用いて取得される破損したサンプルを与えられる場合について考察した。この状況において、あらゆるカーネルを扱う汎用的なＦＲＩ再構成方法を開発した。入力信号の正確なパラメータを見つけようとする元のＦＲＩフレームワークと異なり、指数関数の近似再構成に基づく入力の近似復元を提案している。我々の新しいアプローチの利点は、あらゆるサンプリングカーネルに適用可能であるとともに、過去に使用された特定の種類のカーネルよりも安定性が高く正確な再構成を提供する点である。
【０１６６】
＜付録Ａ＞
一般化Strang-Fix条件
指数関数再生カーネルは、そのシフトされたバージョンの線形結合と共に、
0006397001.tif 000160
の形式の多項式指数関数を生成できる任意の関数φ（ｔ）である。パラメータαｍは複素数であり、ｍ＝０,..., Ｐ及びｒ＝０,..., Ｒである。指数関数再生カーネルは、多項式再生カーネルのファミリーの一般化である。後者は、カーネルψ（ｔ）が多項式を再生できることを記述する既知のStrang-Fix条件[1]、[33]、すなわち、
0006397001.tif 000161
の場合、及びその場合に限り、
0006397001.tif 000162
という条件によって特徴付けられる。ここで、
0006397001.tif 000163
はψ（ｔ）のフーリエ変換であり、
0006397001.tif 000164
はそのｒ次導関数を表す。
【０１６７】
多項式再生カーネルと同様の方法で、指数関数再生カーネルについてのStrang-Fix条件の更に一般的な形式を導出できる。より具体的には、カーネルφ（ｔ）は指数関数多項式、すなわち、
0006397001.tif 000165
の場合、及びその場合に限り、
0006397001.tif 000166
を再生できる。ここで、
0006397001.tif 000167
はφ（ｔ）の両側ラプラス変換のｒ次導関数を表す。証明は、ｔｒを明らかに再生する関数ψ（ｔ）＝ｅｘｐ｛−αｍｔ｝φ（ｔ）を考察することにより、多項式再生カーネルについてのStrang-Fix条件から取得される。従って、ψ（ｔ）は、多項式再生カーネルであるため、Strang-Fix条件を満たす。更に、
0006397001.tif 000168
であるため、ψ（ｔ）のフーリエ変換はφ（ｔ）のラプラス変換に関する。
0006397001.tif 000169
にStrang-Fix条件を単に適用することにより、
0006397001.tif 000170
についての一般化Strang-Fix条件に達し、これで証明が完成する。
【０１６８】
＜付録Ｂ＞
ｅＭＯＭＳのファミリーはＳｏＳカーネル及びディリクレカーネルを含む
０を中心とした、サポートＰ＋１の指数関数再生カーネル
0006397001.tif 000171
について考える。ここで、
0006397001.tif 000172
とし、
0006397001.tif 000173
は（１１）により与えられるフーリエ変換を有する。Ｐが偶数である場合に解析を制限し、以下の指数関数パラメータを使用する。
0006397001.tif 000174
ここで、ｍ＝０,..., Ｐである。次に、φ０（ｔ）の周期的拡張を使用する。すなわち、
0006397001.tif 000175
であり、これは、以下のように書くことができる。
0006397001.tif 000176
ここで、最後の項はポアソンの和公式の適用から得られる。Ｐが奇数の場合は、２（Ｐ＋１）にわたり周期化することにより、同様に導出可能である。また、シフトされたカーネルφ０（ｔ）のフーリエ変換は、以下のように書くことができる。
0006397001.tif 000177
【０１６９】
以下のように表記する場合、第IV−Ｂ節のカーネルとは異なるが依然として[30]の一般化指数関数再生カーネルの特定のサブファミリーに対応する最高次及び最小サポート（ｅＭＯＭＳ）の指数関数再生カーネルが設計される。
0006397001.tif 000178
【０１７０】
（３１）において、フーリエ変換
0006397001.tif 000179
はωｋ＝２πｋ／（Ｐ＋１）で評価される。（３３）を考慮すると、範囲ｋ＝−Ｐ／２,..., Ｐ／２に対して
0006397001.tif 000180
であることがわかる。また、
0006397001.tif 000181
に等しい積（３２）における項を見つけられないため、他のあらゆるｋに対して
0006397001.tif 000182
であることがわかる。従って、式（３１）を次式に帰着できる。
0006397001.tif 000183
【０１７１】
なお、全てのｋに対して値ηｋ＝１である場合、式（３４）は、周期Ｐ＋１のディリクレカーネルの１周期、
0006397001.tif 000184
に帰着する。これはまさに第IV−Ｂ節のｅＭＯＭＳカーネルのファミリーである。
【０１７２】
最後に、（３４）の１周期について考え、
0006397001.tif 000185
を使用する。更に、サンプル数Ｎ＝Ｐ＋１を制限し、τ＝ＮＴ＝（Ｐ＋１）Ｔと表記する。その場合、ＳｏＳカーネル[3]の時間領域定義、
0006397001.tif 000186
を得る。この場合、サンプル数Ｎは奇数である必要があり、インデックスの集合は、
0006397001.tif 000187
である。
【０１７３】
0006397001.tif 000188
0006397001.tif 000189
【図１】
0006397001.tif 000190
【図２】
0006397001.tif 000191
【図３】
0006397001.tif 000192
【図４】
0006397001.tif 000193
【図５】
0006397001.tif 000194
【図６】
0006397001.tif 000195
【図７】
0006397001.tif 000196