特許 6041325 信号処理装置、信号処理方法およびプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6041325

(24)【登録日】2016年11月18日

(45)【発行日】2016年12月7日

(54)【発明の名称】信号処理装置、信号処理方法およびプログラム

(51)【国際特許分類】

   G01S 13/90 20060101AFI20161128BHJP

   A61B 5/055 20060101ALI20161128BHJP

【ＦＩ】

   G01S13/90 127

   A61B5/05 374

【請求項の数】7

【全頁数】24

(21)【出願番号】特願2014-511131(P2014-511131)

(86)(22)【出願日】2013年2月22日

(86)【国際出願番号】JP2013054596

(87)【国際公開番号】WO2013157299

(87)【国際公開日】20131024

【審査請求日】2016年2月8日

(31)【優先権主張番号】特願2012-95589(P2012-95589)

(32)【優先日】2012年4月19日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】304021417

【氏名又は名称】国立大学法人東京工業大学

(74)【代理人】

【識別番号】100099759

【弁理士】

【氏名又は名称】青木 篤

(74)【代理人】

【識別番号】100092624

【弁理士】

【氏名又は名称】鶴田 準一

(74)【代理人】

【識別番号】100114018

【弁理士】

【氏名又は名称】南山 知広

(74)【代理人】

【識別番号】100119987

【弁理士】

【氏名又は名称】伊坪 公一

(74)【代理人】

【識別番号】100161089

【弁理士】

【氏名又は名称】萩原 良一

(72)【発明者】

【氏名】山田 功

(72)【発明者】

【氏名】北原 大地

【審査官】 山下 雅人

(56)【参考文献】

【文献】 特開２００７−３３３５８３（ＪＰ，Ａ）

【文献】 特開２００３−７０７６３（ＪＰ，Ａ）

【文献】 特開平１０−９０１１２（ＪＰ，Ａ）

【文献】 特開２０１０−２５２３３１（ＪＰ，Ａ）

【文献】 特開平５−７５４７６（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｓ １３／９０

Ａ６１Ｂ ５／０５５

(57)【特許請求の範囲】

【請求項1】

１組の信号列を取得する信号取得部と、
前記信号取得部により取得された前記１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出する第１の算出部と、
前記複数の区間のそれぞれにおいて、前記第１の算出部が算出した前記第１の関数と前記第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に当該区間内のいずれかの点を代入した値の数列を算出する第２の算出部と、
前記第２の算出部により算出された前記数列の各項の符号に基づいて、前記いずれかの点における前記１組の信号列の位相を決定する位相決定部と、
前記複数の区間のそれぞれにおいて前記位相決定部により決定された位相の信号列を出力する信号出力部と
を備える信号処理装置。

【請求項2】

前記位相決定部が、前記数列内の隣接する２項間で当該数列の値の符号が変化した回数に基づいて、前記いずれかの点における前記１組の信号列の位相がもつπの整数倍の不定部分の値を決定する、請求項１に記載の信号処理装置。

【請求項3】

前記第１の算出部が、前記１組の信号列のそれぞれについて、当該信号列の各標本点を通過し、かつ隣り合う区間の多項式同士が滑らかに連続する区分的多項式を算出する、請求項１または２に記載の信号処理装置。

【請求項4】

前記第２の算出部が、前記複数の区間のそれぞれにおいて、前記第１の関数と前記第２の関数に関する終結式行列の小行列である複数の部分終結式行列の行列式に基づいて、前記多項式剰余列から得られる前記数列に代えて、当該数列の各項を正の定数倍した数列を算出する、請求項１から３のいずれか１項に記載の信号処理装置。

【請求項5】

前記位相決定部が位相を決定した点の中に位相差がπより大きい隣接する２点が含まれる位相の信号列を、前記信号出力部が出力する、請求項１から４のいずれか１項に記載の信号処理装置。

【請求項6】

１組の信号列を取得するステップと、
取得された前記１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出するステップと、
前記複数の区間のそれぞれにおいて、算出された前記第１の関数と前記第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に当該区間内のいずれかの点を代入した値の数列を算出するステップと、
算出された前記数列の各項の符号に基づいて、前記いずれかの点における前記１組の信号列の位相を決定するステップと、
前記複数の区間のそれぞれにおいて決定された位相の信号列を出力するステップと
を含む信号処理方法。

【請求項7】

コンピュータに、
１組の信号列を取得する機能と、
取得された前記１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出する機能と、
前記複数の区間のそれぞれにおいて、算出された前記第１の関数と前記第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に当該区間内のいずれかの点を代入した値の数列を算出する機能と、
算出された前記数列の各項の符号に基づいて、前記いずれかの点における前記１組の信号列の位相を決定する機能と、
前記複数の区間のそれぞれにおいて決定された位相の信号列を出力する機能と
を実現させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、１組の信号列の位相を決定する信号処理に関する。

【背景技術】

【0002】

リモートセンシングや磁気共鳴映像法（ｍａｇｎｅｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｍａｇｉｎｇ：ＭＲＩ）などの信号処理では位相情報を必要とし、位相情報を精度よく求めることが重要である。そうした信号処理では、２つの信号間の位相（干渉信号など）を連続値として求めることが必要であり、この問題は「位相アンラップ（Ｐｈａｓｅ Ｕｎｗｒａｐｐｉｎｇ）問題」と呼ばれている。

【0003】

特許文献１には、移動する目標に対して電波を送信し、目標からの反射波を受信して目標の画像を得るため、目標の移動に伴う反射波の位相の変化を補償する位相補償回路であって、反射波の受信信号列から特定のレンジビンのデータ列を取り出すレンジビン決定手段と、取り出されたデータ列から位相の時間変化を算出し、この算出した位相の時間変化の内、位相の折り返しを除去する位相アンラップ手段と、位相の折り返しが除去された位相の時間変化を示すデータに最小二乗法を用いたフィッティングを行い、位相の時間変化の二次以上の変化成分を検出し、この検出結果から二次以上の変化成分を低減させる位相補償量を算出する位相補償量算出手段と、得られた位相補償量に応じて反射波の受信信号列を補償する位相補償手段とを有する位相補償回路が記載されている。

【0004】

特許文献２には、３次元温度分布情報を含む３次元ＭＲ撮影シーケンスを実行し、得られる３次元複素ＭＲ画像から３次元位相分布を計算し、次に３次元位相アンラップ処理を行い、アンラップ処理後の３次元位相分布から３次元温度分布を計算し、この計算結果に対し、ボリュームレンダリング処理を行って、３次元温度画像を作り表示する、ＭＲＩ装置を用いた３次元温度計測方法が記載されている。

【0005】

特許文献３には、受信器において狭帯域干渉信号を相殺する方法であって、受信入力信号から参照信号を減算するステップと、アークタンジェント関数に基づいて減算の結果の位相を計算するステップと、アークタンジェント関数により生じられるモジュロ２π制限を除外することによって、アークタンジェント関数からの出力信号にアンラップ関数を実行するステップであって、これにより、絶対値位相表現を生成するステップと、所定の時間だけずらされる位相表現値を比較することによって周波数オフセットを決定するステップと、決定された周波数オフセットの結果に基づいて狭帯域干渉物を相殺するステップと、を有する方法が記載されている。

【0006】

非特許文献１〜３には、２つの実多項式関数Ｂ_（０）（ｒ），Ｂ_（１）（ｒ）について、ユークリッドの互除法を拡張した計算法により、連続関数である位相θ（ｒ）＝ｔａｎ^−１（Ｂ_（１）（ｒ）／Ｂ_（０）（ｒ））を求める代数的な解法が記載されている。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】特許第３３９５６０６号公報

【特許文献2】特開２０００−３００５３６号公報

【特許文献3】特開２００７−５２１７２９号公報

【非特許文献】

【0008】

【非特許文献1】I. Yamada, K. Kurosawa, H. Hasegawa, K. Sakaniwa, “Algebraic Multidimensional Phase Unwrapping and Zero Distribution of Complex Polynomials − Characterization of Multivariate Stable Polynomials”, IEEE Transactions on Signal Processing, 1998, vol. 46, no. 6, p.1639−1664

【非特許文献2】I. Yamada, N. K. Bose, “Algebraic Phase Unwrapping and Zero Distribution of Polynomial for Continuous−Time Systems”, IEEE Transactions on Circuits and Systems I： Fundamental Theory and Applications, 2002, vol. 49, no. 3, p.298−304

【非特許文献3】I. Yamada, K. Oguchi, “High−Resolution Estimation of the Directions−of−Arrival Distribution by Algebraic Phase Unwrapping Algorithms”, Multidimensional Systems and Signal Processing, Springer, 2011, vol. 22, no. 1−3, p.191−211

【発明の概要】

【発明が解決しようとする課題】

【0009】

実際の信号処理での信号は必ずしも実多項式関数で表されず、多項式で近似しようとしても次数が高くなると数値計算が不安定になる。このため、非特許文献１〜３の代数的な解法を実際の信号処理に適用して連続な位相を求めることは困難である。本発明の目的は、１組の信号列の連続した位相を求める数値計算の安定性を向上させることにある。

【課題を解決するための手段】

【0010】

かかる目的のもと、本発明は、１組の信号列を取得する信号取得部と、信号取得部により取得された１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出する第１の算出部と、複数の区間のそれぞれにおいて、第１の算出部が算出した第１の関数と第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に区間内のいずれかの点を代入した値の数列を算出する第２の算出部と、第２の算出部により算出された数列の各項の符号に基づいて、いずれかの点における１組の信号列の位相を決定する位相決定部と、複数の区間のそれぞれにおいて位相決定部により決定された位相の信号列を出力する信号出力部とを備える信号処理装置を提供する。
位相決定部が、数列内の隣接する２項間で数列の値の符号が変化した回数に基づいて、いずれかの点における１組の信号列の位相がもつπの整数倍の不定部分の値を決定するものであってよい。
第１の算出部が、１組の信号列のそれぞれについて、信号列の各標本点を通過し、かつ隣り合う区間の多項式同士が滑らかに連続する区分的多項式を算出するものであってよい。
第２の算出部が、複数の区間のそれぞれにおいて、第１の関数と第２の関数に関する終結式行列の小行列である複数の部分終結式行列の行列式に基づいて、多項式剰余列から得られる数列に代えて、その数列の各項を正の定数倍した数列を算出するものであってよい。
位相決定部が位相を決定した点の中に位相差がπより大きい隣接する２点が含まれる位相の信号列を、信号出力部が出力するものであってよい。
また、本発明は、１組の信号列を取得するステップと、取得された１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出するステップと、複数の区間のそれぞれにおいて、算出された第１の関数と第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に区間内のいずれかの点を代入した値の数列を算出するステップと、算出された数列の各項の符号に基づいて、いずれかの点における１組の信号列の位相を決定するステップと、複数の区間のそれぞれにおいて決定された位相の信号列を出力するステップとを含む信号処理方法を提供する。
また、本発明は、コンピュータに、１組の信号列を取得する機能と、取得された１組の信号列のそれぞれを複数の区間ごとに多項式で近似した区分的多項式である第１の関数および第２の関数を算出する機能と、複数の区間のそれぞれにおいて、算出された第１の関数と第２の関数にユークリッドの互除法を適用することにより得られる多項式剰余列に区間内のいずれかの点を代入した値の数列を算出する機能と、算出された数列の各項の符号に基づいて、いずれかの点における１組の信号列の位相を決定する機能と、複数の区間のそれぞれにおいて決定された位相の信号列を出力する機能とを実現させるためのプログラムを提供する。

【発明の効果】

【0011】

本発明によれば、１組の信号列の連続した位相を求める数値計算の安定性を、本発明の構成を備えない場合と比べて向上させることができる。

【図面の簡単な説明】

【0012】

【図1】干渉合成開口レーダにより地表高度を計測する方法の原理を説明するための図である。

【図2】位相アンラップを行う位相データの例を示した画像である。

【図3】本実施形態の信号処理装置の機能構成例を示したブロック図である。

【図4】本実施形態の信号処理装置の動作例の概略を示したフローチャートである。

【図5】本実施形態の信号処理装置におけるスプライン算出部の動作例を示したフローチャートである。

【図6】第１の実施形態の信号処理装置における多項式列算出部の動作例を示したフローチャートである。

【図7】第２の実施形態の信号処理装置における多項式列算出部の動作例を示したフローチャートである。

【図8】位相アンラップの数値実験に用いた一方の入力信号のグラフである。

【図9】位相アンラップの数値実験に用いた他方の入力信号のグラフである。

【図10】図８および図９の入力信号のラップ位相のグラフである。

【図11】図８および図９の入力信号のアンラップ位相のグラフである。

【図12】図８および図９の入力信号について位相アンラップを行った数値実験の結果を示したグラフである。

【図13】本実施形態を実現可能なコンピュータのハードウェア構成図である。

【発明を実施するための形態】

【0013】

以下、添付図面を参照して、本発明の実施の形態について詳細に説明する。
本実施形態の信号処理装置は、取得した１組の信号列Ｆ，Ｇから位相Θ＝ｔａｎ^−１（Ｇ／Ｆ）を計算し、出力する。この位相Θはｍπ（ｍは整数）の不定性があり一意に求まらないが、Θが連続になる適切な整数値ｍを定めることを「位相アンラップ」という。本実施形態の信号処理装置は、この位相アンラップを行う。

【0014】

まず、本実施形態の信号処理装置が適用されるシステムの例について説明する。
合成開口レーダ（Ｓｙｎｔｈｅｔｉｃ Ａｐｅｒｔｕｒｅ Ｒａｄａｒ：ＳＡＲ）を用いて地表の情報を計測するシステムがある。ＳＡＲは、飛行機や衛星に搭載されたレーダから地上に向けて電波を発射し、対象物から反射された電波を自身のアンテナで受信しており、アンテナ自身が飛行機や衛星によって移動することで仮想的に大きな開口面を実現したレーダである。合成開口レーダの１つである干渉合成開口レーダ（干渉ＳＡＲ）は、同一地域について得られた２セットの複素画像を干渉させて生成される干渉縞から、地表高度や、地殻変動による地表高度の変化を計測する。

【0015】

その２セットの複素画像をそれぞれＡ_１＝Ａ_０ｅｘｐ（ｉφ_１），Ａ_２＝Ａ_０ｅｘｐ（ｉφ_２）とすると、干渉縞はＡ_１Ａ_２^＊＝|Ａ_０|^２ｅｘｐ（ｉΔφ）である。ここで、Δφ＝φ_１−φ_２は位相差であり、次式で表される。ｒｅａｌは実部、ｉｍａｇは虚部、＊は複素共役を意味する。

【0016】

【数1】

【0017】

地表高度を求めるためにはΔφを求める必要があるが、正接関数ｔａｎには周期πの周期性があるため、Δφの値は一意に決まらない。そこで本実施形態の信号処理装置は、位相アンラップを行って、位相Δφが連続になる適切な整数値ｍを定める。

【0018】

位相アンラップにより連続な位相Δφが求まれば、地表高度は例えば以下のように求められる。図１は、干渉ＳＡＲにより地表高度を計測する方法の原理を説明するための図である。Ｒ_１，Ｒ_２はそれぞれの軌道１，２のレーダから地表の対象物までの距離、Ｂは軌道１，２のレーダ間の距離、Ｈはレーダの高度、ｈは測定対象の高度であり、図１のように角θ，α，δθを定める。未知数はＲ_１，Ｒ_２，ｈ，θである。すると、Ｒ_１はＲ_１＝Ｂｓｉｎ（θ−α）＋Ｒ_２ｃｏｓδθとなる。実際にはδθ≒０なので、ｃｏｓδθ≒１であり、Ｒ_１＝Ｂｓｉｎ（θ−α）＋Ｒ_２となる。すると、Ｒ_１とＲ_２の差ΔＲは、ΔＲ＝Ｒ_１−Ｒ_２＝Ｂｓｉｎ（θ−α）である。したがって、位相差Δφは次式で表される。

【0019】

【数2】

【0020】

ここで、λは電波の波長である。また、ｈ＝Ｈ−Ｒ_１ｃｏｓθなので、Δφからθを求め、さらにｈを求めることで、地表高度ｈが計測される。

【0021】

本実施形態の信号処理装置が適用されるシステムの別の例として、ＭＲＩを用いた核磁気共鳴画像システムがある。ＭＲＩは、体内のプロトン（水素原子の原子核（陽子））の密度と状態を画像化する技術である。人体は９割ほどが水と脂肪によって構成されており、それらはどちらも水素原子を多く含むため、プロトンの密度と状態からおおよその体の構造がわかる。ＭＲＩの撮像法の１つであるグラジエントフィールドエコー（ｇｒａｄｉｅｎｔ ｆｉｅｌｄ ｅｃｈｏ：ＧＲＥ）法では、水の信号と脂肪の信号が混合した状態から両者の成分を分離して撮像する、Ｄｉｘｏｎ法という手法が用いられる。

【0022】

Ｄｉｘｏｎ法では、水と脂肪の分子構造の違いから核磁気共鳴スペクトルにシフトが生じる（化学シフト現象）ことを利用する。外部から適切にパルスを加えると、次式のように、水と脂肪の信号が同位相になっている信号と、水と脂肪の信号が逆位相になっている信号が得られる。
同位相（ｉｎ−ｐｈａｓｅ） ｓ_０＝（Ｗ＋Ｆ）ｅｘｐ（ｉθ_０）
逆位相（ｏｕｔ−ｐｈａｓｅ） ｓ_１＝（Ｗ−Ｆ）ｅｘｐ（ｉ（θ_０＋φ））
ここで、Ｗ，Ｆ，θ_０は定数である。すると、ｓ_０^＊ｓ_１／ｓ_０ｓ_１^＊＝ｅｘｐ（ｉ２φ）であるから、位相差φは以下のように求められる。

【0023】

【数3】

【0024】

本実施形態の信号処理装置は、ここで適切な整数値ｍを定めて連続な位相２φを求めるために用いられる。連続な位相２φが求まれば、次式により水画像と脂肪画像が得られる。

【0025】

【数4】

【0026】

このように、信号処理や画像処理では、２つの実数値連続関数の組ｆ（ｘ），ｇ（ｘ）が与えられ、ｔａｎΘ（ｘ）＝ｇ（ｘ）／ｆ（ｘ）を満足する連続関数Θ（ｘ）を決定することが要請される。一般に独立変数は１つとは限らず、２変数の場合の応用が特に重要である。しかし、２変数の位相アンラップも、１変数の位相アンラップを各々の変数について計算することで実現できる。このため、本実施形態では簡単化のため１つの実数変数ｘの関数を考える。

【0027】

一般に、２つの実数の組（ｘ，ｙ）（ただしｘ≠０）が与えられるとき、ｔａｎθ＝ｙ／ｘを満たす実数θを（ｘ，ｙ）がなす位相という。しかし、上記の通り、正接関数ｔａｎは周期πの周期関数であるため、このθは一意に決まらない。すなわち、（ｘ，ｙ）がなす位相は、ｔａｎθ_０＝ｙ／ｘを満足する唯一のθ_０∈（−π／２，π／２）（主値という）を用いたθ＝θ_０＋ｍπ（ｍは任意の整数）のいずれでもよい。

【0028】

上記の連続関数Θ（ｘ）を決定する位相アンラップ問題は、ｆ（ｘ）≠０となるｘでは主値にπの適当な整数倍を加える補正を施し、かつｆ（ｘ）＝０となるｘでは適当な実数値をΘ（ｘ）に割り当てて連続関数Θを構成するという問題である。従来の位相アンラップの方法は試行錯誤的なものであり、系統的かつ安定な位相アンラップを実現する方法は未だ確立されていない。多くの最先端の計測システムの性能は位相アンラップの成否に大きく左右されるため、安定した位相アンラップを実現する技術が確立できれば大きな波及効果が期待される。

【0029】

図２は、位相アンラップを行う位相データの例を示した画像である。これらの画像は、法政大学情報科学部花泉弘研究室のＷｅｂページ［ｏｎｌｉｎｅ］［平成２４年４月２日検索］（ＵＲＬ：http：／／cis.k.hosei.ac.jp／research／laboratory／digital／hanaizumi／remote＿sensing.html）から引用した。

【0030】

図２（ａ）は、位相アンラップを行う前の２次元の位相データの画像である。図２の各画像では白黒の濃淡により位相の値の大小を示しており、図２（ａ）では位相が−πからπの間に制限されている。このため、値が−πからπまたはπから−πに不連続に変わる点が存在することから、濃淡の変化が不連続な、白と黒にはっきり分かれた部分が見られる。このように、値（値域）が−πからπの間などに制限された位相のことを、以下では「ラップ位相（Ｗｒａｐｐｅｄ Ｐｈａｓｅ）」という。

【0031】

図２（ｂ）は、図２（ａ）の位相データに対して適切に位相アンラップを行って得られる画像である。図２（ｂ）の場合、位相の値は−πからπの間に制限されず、値が不連続に変わる点が存在しない。位相がすべての点で連続的に変化していることを反映して、濃淡が連続的に変化していることがわかる。このように、位相アンラップを行って値域の制限を取り除いた位相のことを、以下では「アンラップ位相（Ｕｎｗｒａｐｐｅｄ Ｐｈａｓｅ）」という。

【0032】

一方、図２（ｃ）は、図２（ａ）の位相データに対する位相アンラップに失敗した場合の画像である。すなわち、図２（ｃ）では、位相アンラップを行った後でも位相が２πだけ不連続に変わる点が存在するため、濃淡の変化が不連続な部分が見られる。従来の試行錯誤的な位相アンラップの方法では、隣接する標本点間のアンラップ位相の変動が±π以下であるという仮定に基づいているため、例えば位相が急激に変化する点がある場合、図２（ｃ）のように位相アンラップが失敗する。

【0033】

本実施形態の信号処理装置は、必ずしも実多項式関数で表されない１組の信号列について、連続なアンラップ位相を系統的な方法で求める。そのために、本実施形態の信号処理装置は、取得した１組の信号列をそれぞれスプライン関数で近似し、得られたスプライン関数の各区間の多項式を用いて位相アンラップを行う。位相アンラップは、スプライン関数の多項式にユークリッドの互除法を適用することにより多項式列を算出し、各区間の位相を求める点におけるその多項式列の値を並べた数列の符号の変化の回数を調べ、その変化の回数に基づいてπの整数倍の不定部分を決定するという手順で行う。

【0034】

なお、本実施形態において、「連続」という用語は、関数に対しては数学的な意味で用いるが、実際のシステムで扱われる離散的な信号の位相に対しては、図２（ａ）や図２（ｃ）で見られるような値の変化の不自然な「飛び」がないという意味で用いる。

【0035】

＜第１の実施形態＞
図３は、本実施形態の信号処理装置１０の機能構成例を示したブロック図である。図示するように、信号処理装置１０は、信号取得部１１と、区間決定部１２と、スプライン算出部１３と、関数格納部１４と、多項式列算出部１５と、数列生成部１６と、符号カウント部１７と、アンラップ処理部１８と、信号出力部１９とを備える。

【0036】

信号取得部１１は、処理対象となる１組の信号列を取得する。例えば、上記の干渉ＳＡＲの場合はＡ_１Ａ_２^＊の実部と虚部を、ＭＲＩの場合はｓ_０^＊ｓ_１／ｓ_０ｓ_１^＊の実部と虚部をそれぞれ１組の信号列として取得する。以下では、信号取得部１１が取得する１組の信号列のそれぞれをＦ（ｘ），Ｇ（ｘ）と表す。例えばこれらの信号列が時系列信号である場合には、ｘは時刻を表す。

【0037】

区間決定部１２は、信号取得部１１が取得した信号列からスプライン算出部１３がスプライン関数を求める際の区間を決定する。標本点ｘ_０，ｘ_１，…，ｘ_Ｎでの信号Ｆ（ｘ），Ｇ（ｘ）を信号取得部１１が取得したとすると、区間決定部１２は、スプライン算出部１３がスプライン関数を求める際の区間を例えば［ｘ_０，ｘ_１］，…，［ｘ_Ｎ−１，ｘ_Ｎ］と定める。なお、区間［ａ，ｂ］という表記は、端点ａおよびｂを含む区間を意味する。

【0038】

スプライン算出部１３は、信号取得部１１が取得した信号列Ｆ（ｘ），Ｇ（ｘ）をそれぞれ近似する区分的多項式の一例として、スプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）を算出する。スプライン算出部１３がスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）を求める処理の詳細については、図５を用いて後述する。本実施形態では、第１の算出部の一例として、スプライン算出部１３を設けている。

【0039】

関数格納部１４は、スプライン算出部１３が算出したスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）や、多項式列算出部１５が算出した多項式列を格納する。

【0040】

多項式列算出部１５は、スプライン算出部１３が算出したスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）がそれぞれ多項式で表される各区間［ｘ_０，ｘ_１］，…，［ｘ_Ｎ−１，ｘ_Ｎ］について、その区間でのＳ_Ｆ（ｘ）とＳ_Ｇ（ｘ）にユークリッドの互除法を適用して、多項式列を算出する。多項式列算出部１５が多項式列を求める処理は、非特許文献３に示されている手順に従って行う。その詳細については、図６を用いて後述する。

【0041】

数列生成部１６は、多項式列算出部１５が多項式列を算出した各区間内の点ｘにおけるその多項式列の値を並べた数列を生成する。例えば、区間［ｘ_０，ｘ_１］については、その区間内のｘ_０における各多項式の値（すなわち、ｘ_０を代入したときの値）を計算し、その計算結果を順に並べて数列とする。本実施形態では、第２の算出部の一例として、多項式列算出部１５と数列生成部１６を設けている。

【0042】

符号カウント部１７は、数列生成部１６が生成した数列内の隣接する２項間でその数列の値の符号が変化した回数を数える。

【0043】

アンラップ処理部１８は、数列生成部１６が数列を生成した点ｘ＝ｔにおける信号列Ｆ（ｔ），Ｇ（ｔ）の位相がもつπの整数倍の不定部分の値を、符号カウント部１７が得た符号変化の回数に基づき決定する。アンラップ処理部１８がπの整数倍の値を決定する処理の詳細についても後述する。本実施形態では、位相決定部の一例として、符号カウント部１７とアンラップ処理部１８を設けている。

【0044】

信号出力部１９は、アンラップ処理部１８が各点で求めた位相の値を信号列として出力する。本実施形態のアンラップ処理部１８が決定する位相は各点で連続的に変化するものであるため、信号出力部１９が出力する位相の信号列も連続的なものになる。

【0045】

以下では、本実施形態の信号処理装置１０の動作について説明する。図４は、信号処理装置１０の動作例の概略を示したフローチャートである。

【0046】

まず、標本点ｘ_０，ｘ_１，…，ｘ_Ｎでの１組の信号列Ｆ（ｘ），Ｇ（ｘ）を信号取得部１１が取得する（ステップ１０）。そして、ｘの区間［ｘ_０，ｘ_１］，…，［ｘ_Ｎ−１，ｘ_Ｎ］を区間決定部１２が決定した上で、それぞれの信号列を近似するスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）をスプライン算出部１３が算出する（ステップ２０）。算出されたスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）は、関数格納部１４が格納する。

【0047】

そして、ｘの最初の区間［ｘ_０，ｘ_１］について、関数格納部１４に格納されているスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）から、多項式列算出部１５が多項式列を算出する（ステップ３０）。区間［ｘ_０，ｘ_１］においてＳ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）がそれぞれ多項式ｆ_０（ｘ），ｇ_０（ｘ）で表されるとすると、多項式列算出部１５はこのｆ_０（ｘ），ｇ_０（ｘ）にユークリッドの互除法を適用して、多項式列｛ψ_０（ｘ），…，ψ_ｑ（ｘ）｝を算出する。ここで、ｑはψ_ｑ（ｘ）が０次（定数）となる番号である。例えばψ_０（ｘ）とψ_１（ｘ）を３次とすると、この多項式列は、例えばψ_２（ｘ）が２次、ψ_３（ｘ）が１次、ψ_４（ｘ）が０次（定数）のように、次数が順に下がる有限個の列となる（この例ではｑ＝４である）。本実施形態では、この多項式列のことを「多項式剰余列」とも呼ぶ。算出された多項式列｛ψ_０（ｘ），…，ψ_ｑ（ｘ）｝は、関数格納部１４が格納する。

【0048】

次に、区間［ｘ_０，ｘ_１］内の点ｘ＝ｔについて、数列生成部１６が多項式列｛ψ_０（ｘ），…，ψ_ｑ（ｘ）｝の値を計算し、それらの値を並べた数列｛ψ_０（ｔ），…，ψ_ｑ（ｔ）｝を作る（ステップ４０）。そして符号カウント部１７が、その数列｛ψ_０（ｔ），…，ψ_ｑ（ｔ）｝について、隣接する２項間で符号が変化する回数を数える（ステップ５０）。この符号変化の回数をＶ｛ψ（ｔ）｝と書く。

【0049】

符号変化の回数Ｖ｛ψ（ｔ）｝の求め方について、具体的に説明する。例えば、数列の項数が６であり、各項の符号が｛＋，＋，−，０，０，＋｝であったとする。第１項から第２項は符号が正のまま変わらないが、第２項から第３項にかけて正から負に変化するので、これを１回と数える。第４項と第５項は０であるが、このように０の項は飛ばして次に０でない第６項を見る。すると、第３項から第６項にかけて負から正に変化するので、これを１回と数える。よって、この数列についてはＶ｛ψ（ｔ）｝＝２である。

【0050】

非特許文献１〜３によると、ｘ＝ｔでの多項式ｆ（ｘ），ｇ（ｘ）のアンラップ位相θ（ｔ）は、次式で与えられることが証明されている。

【0051】

【数5】

【0052】

ここで、ｔ_０はアンラップ位相θ（ｔ）を求めるための基準点である。右辺第２項と第３項は（−π／２，π／２）の範囲の値（主値）であり、Ｖ｛ψ（ｔ）｝はステップ５０で求めたｘ＝ｔでの符号変化の回数である。また、Ｖ｛ψ（ｔ_０）｝は、基準点ｘ＝ｔ_０での多項式列｛ψ_０（ｘ），…，ψ_ｑ（ｘ）｝の値を並べた数列の符号変化の回数である。すなわち、主値にπの何倍を加えればアンラップ位相が得られるのかということは、ステップ５０で求めた符号変化の回数Ｖ｛ψ（ｔ）｝から上式により厳密に求められる。

【0053】

スプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）は区間ごとに異なる多項式で表されることから、本実施形態では上記の基準点ｔ_０も区間ごとに定める。例えば、「区間［ｘ_０，ｘ_１］ではｘ_０，…，区間［ｘ_Ｎ−１，ｘ_Ｎ］ではｘ_Ｎ−１」というように、区間の開始点（左端の点）を基準点ｔ_０とすればよい。

【0054】

したがって、アンラップ処理部１８は、区間［ｘ_０，ｘ_１］については基準点をｔ_０＝ｘ_０として、関数格納部１４に格納されている区間［ｘ_０，ｘ_１］での多項式ｆ_０（ｘ），ｇ_０（ｘ）とステップ５０で求めたＶ｛ψ（ｔ）｝を用いて、数５を計算する。すなわち、アンラップ処理部１８は、次式を計算することにより、区間［ｘ_０，ｘ_１］内のｘ＝ｔでのアンラップ位相θ（ｔ）を求める（ステップ６０）。

【0055】

【数6】

【0056】

そして、区間［ｘ_０，ｘ_１］内の別の点についてもアンラップ位相を求める場合はステップ４０に戻り、次の区間［ｘ_１，ｘ_２］の処理に移る場合はステップ８０に進む（ステップ７０）。なお、各区間内のどの点でアンラップ位相を計算してもよいが、本実施形態では、少なくとも各区間の端点ｘ_０，ｘ_１，…，ｘ_Ｎでのアンラップ位相は計算しておく。

【0057】

次の区間がある場合、処理はステップ３０に戻る（ステップ８０でＹｅｓ）。区間［ｘ_１，ｘ_２］でアンラップ位相を求めるときは、アンラップ処理部１８は、基準点をｔ_０＝ｘ_１として、区間［ｘ_１，ｘ_２］での多項式ｆ_１（ｘ），ｇ_１（ｘ）とステップ５０で求めたＶ｛ψ（ｔ）｝を用いて、数６と同様の計算を行う。区間［ｘ_２，ｘ_３］以降についても同様である。

【0058】

一方、次の区間がない場合、すなわち、区間［ｘ_Ｎ−１，ｘ_Ｎ］までアンラップ位相を求めた場合、処理はステップ９０に進む（ステップ８０でＮｏ）。そして最後に、今までに求めた各アンラップ位相の値を信号出力部１９が信号列として出力する（ステップ９０）。以上で信号処理装置１０の動作は終了する。

【0059】

次に、図４のステップ２０でスプライン算出部１３がスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）を算出する処理について説明する。
本実施形態では、信号列Ｆ（ｘ），Ｇ（ｘ）の標本値を用いて、スプライン算出部１３がＦ（ｘ），Ｇ（ｘ）のそれぞれをスプライン関数で近似することにより、各標本点を通る２つの区分的多項式を求める。ｎ次のスプライン関数Ｓ（ｘ）は、「各区間においてＳ（ｘ）がｎ次以下の多項式であり、かつ定義域全体でＳ（ｘ）とその（ｎ−１）次以下の導関数が連続な関数」と定義される。すなわち、スプライン関数は、複数の多項式を互いに接続した区分的多項式であって、多項式同士のつなぎ目（節点という）も含めて数学的な意味で滑らかな関数である。

【0060】

一般に、関数ｆ（ｘ）の有限個の標本点（ｘ_ｋ，ｆ（ｘ_ｋ））（ｋ＝１，２，３，…，ｐ）が与えられたときに、これらの標本点を通過する何らかの関数ｈ（ｘ）を求めるには、ｈ（ｘ）を（ｐ―１）次多項式として、その多項式の係数に関する連立方程式を解けばよい。しかし、標本点の個数ｐが大きくなると、多項式ｈ（ｘ）は標本点以外のｘで大きく振動するため、１つの多項式で補間することは適切とはいえない。一方、スプライン関数を用いれば、１つの多項式で補間する場合より低次なものにすることができ、振動の小さな関数になる。そこで本実施形態では、信号列Ｆ（ｘ），Ｇ（ｘ）のそれぞれに対してスプライン関数の補間アルゴリズムを利用する。

【0061】

次に、標本点を３次のスプライン関数で補間する場合の計算方法について説明する。ここでは、標本点とスプライン関数の節点とを一致させる。つまり、標本（ｘ_ｋ，ｙ_ｋ）（ｋ＝０，１，…，Ｎかつｘ_０＜ｘ_１＜…＜ｘ_Ｎ）が与えられたときに、区間［ｘ_０，ｘ_１］ではある３次関数、区間［ｘ_１，ｘ_２］では別のある３次関数というように、標本点間をそれぞれ別の３次関数でつなぐものとする。求めるスプライン関数をＳ（ｘ）とすると、Ｓ（ｘ）は、標本点（ｘ_ｋ，ｙ_ｋ）（ｋ＝０，１，…，Ｎ）をすべて通過し、かつ区間（ｘ_０，ｘ_Ｎ）上で２階微分までが連続になる。

【0062】

さて、標本点ｘ_ｋでの２階微分の値をＭ_ｋ（ｋ＝０，１，…，Ｎ）と書くと、区間［ｘ_ｋ−１，ｘ_ｋ］（ｋ＝１，２，…，Ｎ）でのＳ（ｘ）の２階微分Ｓ’’（ｘ）は次式で表される。

【0063】

【数7】

【0064】

ただし、ｈ_ｋ＝ｘ_ｋ−ｘ_ｋ−１である。これを２回積分し、Ｓ（ｘ_ｋ−１）＝ｙ_ｋ−１，Ｓ（ｘ_ｋ）＝ｙ_ｋを用いると、区間［ｘ_ｋ−１，ｘ_ｋ］でＳ（ｘ）は次式のように表される。

【0065】

【数8】

【0066】

よって、後はＭ_ｋ（ｋ＝０，１，…，Ｎ）がわかればＳ（ｘ）が求まる。このＭ_ｋを求めるために１階微分を考えると、上記の結果から、区間［ｘ_ｋ−１，ｘ_ｋ］でのＳ’（ｘ）は次式で表される。

【0067】

【数9】

【0068】

よって、ｘ_ｋ（ｋ＝１，２，…，Ｎ−１）におけるＳ’（ｘ）の左方微分係数と右方微分係数は、次式のように表される。

【0069】

【数10】

【0070】

１回微分の連続性からＳ’（ｘ_ｋ−）＝Ｓ’（ｘ_ｋ＋）なので、次式が導かれる。

【0071】

【数11】

【0072】

未知数がＭ_ｋ（ｋ＝０，１，…，Ｎ）の（Ｎ＋１）個であるのに対し、数１１の方程式は（Ｎ−１）個しかないので、未知数を一意に決めるためにＭ_０＝Ｍ_Ｎ＝０とおく。すると、数１１の方程式から未知数Ｍ_ｋ（ｋ＝１，２，…，Ｎ−１）が求まるので、数８からＳ（ｘ）が求まる。このようにして、信号列Ｆ（ｘ），Ｇ（ｘ）をそれぞれスプライン関数で補間する。

【0073】

なお、本実施形態では一例として３次のスプライン関数を用いるが、次数は何次であってもよい。ただし、Ｓ_Ｆ（ｘ）とＳ_Ｇ（ｘ）は同じ次数にするとよい。また、上記のように標本点とスプライン関数の節点を１対１に対応させる必要はなく、例えばｘ_０からｘ_３までの区間［ｘ_０，ｘ_３］を１つの３次関数で表してもよい。さらに、標本点以外の点を節点としてもよい。ただし、Ｓ_Ｆ（ｘ）とＳ_Ｇ（ｘ）の間では、節点の取り方を一致させるとよい。

【0074】

図５は、スプライン算出部１３の動作例を示したフローチャートである。スプライン算出部１３は、信号取得部１１が取得した信号列Ｆ（ｘ），Ｇ（ｘ）のそれぞれについて、以下の手順でスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）を算出する。

【0075】

まず、標本点ｘ_０，ｘ_１，…，ｘ_Ｎをもとに、区間決定部１２がスプライン関数の節点を定める（ステップ２１）。上記の通り節点を標本点と同じにする必要はないが、ここでは簡単のため、それぞれの標本点を節点として、区間を［ｘ_０，ｘ_１］，…，［ｘ_Ｎ−１，ｘ_Ｎ］と定める。この区間割りは、Ｆ（ｘ）とＧ（ｘ）の間で共通のものとする。

【0076】

次に、スプライン算出部１３は、数１１を解いてＭ_０，…，Ｍ_Ｎを求める（ステップ２２）。その際、数１１においてｈ_ｋ＝ｘ_ｋ−ｘ_ｋ−１とする。また、Ｆ（ｘ）の場合はｙ_ｋ＝Ｆ（ｘ_ｋ）、Ｇ（ｘ）の場合はｙ_ｋ＝Ｇ（ｘ_ｋ）とする。そして、スプライン算出部１３は、１つの区間［ｘ_ｋ−１，ｘ_ｋ］について、数８により３次関数を求める（ステップ２３）。次の区間がある場合はステップ２３に戻り、すべての区間［ｘ_０，ｘ_１］，…，［ｘ_Ｎ−１，ｘ_Ｎ］について３次関数を求め終わったらステップ２５に進む（ステップ２４）。ここまでで、Ｓ_Ｆ（ｘ）は、「［ｘ_０，ｘ_１］において３次関数ｆ_０，…，［ｘ_Ｎ−１，ｘ_Ｎ］において３次関数ｆ_Ｎ−１」と表されている。Ｓ_Ｇ（ｘ）も、「［ｘ_０，ｘ_１］において３次関数ｇ_０，…，［ｘ_Ｎ−１，ｘ_Ｎ］において３次関数ｇ_Ｎ−１」と表されている。そこで、スプライン算出部１３は、このように区間と多項式を対応付けて、Ｓ_Ｆ（ｘ）とＳ_Ｇ（ｘ）をそれぞれ関数格納部１４に格納する（ステップ２５）。以上でスプライン算出部１３の動作は終了する。

【0077】

なお、スプライン補間の計算方法も、これまでに様々なものが提案されている。スプライン算出部１３は、図５の方法に限らず、別の既存の計算方法を利用してもよい。

【0078】

次に、図４のステップ３０で多項式列算出部１５が多項式列を算出する処理について説明する。図６は、第１の実施形態の多項式列算出部１５の動作例を示したフローチャートである。

【0079】

まず、現在の処理対象の区間［ｘ_ｋ−１，ｘ_ｋ］について、多項式列算出部１５は関数格納部１４からスプライン関数Ｓ_Ｆ，Ｓ_Ｇの３次関数ｆ_ｋ，ｇ_ｋを取得し、それぞれをψ_０，ψ_１とする（ステップ３１）。そして、変数ｊをｊ＝１とする（ステップ３２）。多項式列算出部１５はψ_ｊの次数が０でないかどうかを判定し、次数が０でなければステップ３４に進み、０であれば後述するステップ３７に進む（ステップ３３）。

【0080】

ψ_ｊの次数が０でない場合、多項式列算出部１５は、ψ_ｊ−１をψ_ｊで割り（ステップ３４）、その余りを−ψ_ｊ＋１とする（ステップ３５）。そしてｊに１を加算し（ステップ３６）、処理はステップ３３に戻る。一方、ψ_ｊの次数が０である場合は、多項式列算出部１５が多項式列｛ψ_０，…，ψ_ｊ｝を関数格納部１４に格納する（ステップ３７）。以上で多項式列算出部１５の動作は終了する。

【0081】

図６の処理は、ψ_０，ψ_１（すなわち、スプライン関数Ｓ_Ｆ，Ｓ_Ｇの３次関数ｆ_ｋ，ｇ_ｋ）にユークリッドの互除法を適用して多項式列｛ψ_０，…，ψ_ｊ｝を求めるものである。ただし、上記では説明を省略しているが、多項式列｛ψ_０，…，ψ_ｊ｝による数列の符号変化の回数からアンラップ位相を求めるためには、各多項式ψ_ｊがｘのべき乗を因数にもつ（すなわち、０次の項が０になる）場合、そのべき乗を除いた部分をψ_ｊと定義し直して図６の処理を行うとよい。本実施形態でいうユークリッドの互除法は、そのような変形例も含むものとする。

【0082】

一般に、実数値をとる多項式の組（ｆ（ｘ），ｇ（ｘ））が与えられたとき、それらにユークリッドの互除法を適用して得られる多項式剰余列は、スツルム列と呼ばれる条件を満たすことが知られている。（ｆ（ｘ），ｇ（ｘ））から連続なアンラップ位相を求めるには主値にπの適当な整数倍を補正項として加える必要がある。非特許文献１，２では、ｆ（ｘ）が０になるｘの前後の（ｆ（ｘ），ｇ（ｘ））の符号（＋，−）からこの補正項を決定できることが示されている。図６の方法は、スツルム列の性質を利用することにより、ｆ（ｘ）が０になるｘの位置の情報がなくても補正項を正しく求めることを可能にしている。

【0083】

このように、本実施形態の信号処理装置１０は、多項式からアンラップ位相を決定する部分は数学的に厳密な方法を用いている。このため、入力信号列を区分的多項式で上手く近似できれば、精度よく位相アンラップを行うことが可能になる。そして、入力信号列をスプライン関数で近似することにより、標本点間を結ぶ多項式関数を精度よく求めることができ、位相アンラップの精度の向上につながる。

【0084】

＜第２の実施形態＞
第１の実施形態の信号処理装置１０は、アンラップ位相を決定するための多項式列を算出する際に、図６に示したユークリッドの互除法を実行している。しかし、この方法では、多項式の次数が高くなると数値的な不安定性が回避できなくなる。これは、図６の処理に多項式の割り算をするステップ（ステップ３４）があることに起因する（ψ_ｊ＋１（ｘ）は、ψ_ｊ−１（ｘ）をψ_ｊ（ｘ）で割った余りを−１倍して求められる）。計算機で実際に図６の処理を実行しようとすると、この多項式の割り算が正確に行えない場合がある。例えば、１割る１／３の答えは３余り０であるが、計算機では１／３を正確に表すことができないため、実際には余りが０にならない。また、係数を分母と分子に分けてそれぞれを整数として保存した場合でも、ユークリッドの互除法の途中で多項式の係数の桁が大きくなりすぎてしまい、正確に保存できなくなる。つまり、計算機でユークリッドの互除法を素直に実装しようとすると、ある多項式の係数の値が途中で正確でなくなり、その正確でない係数を用いて次の多項式を作るというステップを繰り返してしまうので、途中から多項式の係数が実際のものとずれてしまう。このような数値計算上の誤差が蓄積すると位相アンラップの計算結果に影響を与えるおそれがあるため、ユークリッドの互除法を直接実行せずに多項式列が生成されるようにすることが望ましい。

【0085】

そこで、第２の実施形態の信号処理装置１０は、スプライン算出部１３がスプライン関数を算出した後、図６の処理により得られる多項式剰余列の正の定数倍である多項式列を、多項式列算出部１５が算出する。その多項式列は、以下で説明する部分終結式を計算することにより得られる。ユークリッドの互除法を適用する多項式の組を（ｆ（ｘ），ｇ（ｘ））とすると、部分終結式を計算して得られる多項式の係数は、すべてｆ（ｘ）とｇ（ｘ）の係数の和と差と積を用いて計算される。よって、上記のように正確でない係数の多項式から次の多項式を作り出すという繰り返しは起こらない。第２の実施形態の多項式列算出部１５は、図６のユークリッドの互除法を直接実行しないので、上記のような数値計算上の不安定性が生じることがない。

【0086】

第２の実施形態の信号処理装置１０は、多項式列算出部１５が行う処理内容を除いて、第１の実施形態の信号処理装置１０と同じである。このため、多項式列算出部１５の処理内容以外の部分については説明を省略する。

【0087】

部分終結式の前に、まず終結式について説明する。終結式は、ｍ次多項式ｆ（ｘ）＝ａ_ｍｘ^ｍ＋ａ_ｍ−１ｘ^ｍ−１＋…ａ_０とｎ次多項式ｇ（ｘ）＝ｂ_ｎｘ^ｎ＋ｂ_ｎ−１ｘ^ｎ−１＋…ｂ_０が共通根をもつかどうかを判定するための行列の行列式であり、数１２のような（ｍ＋ｎ）次正方行列の行列式である。なお、行列の中で空欄になっている成分はすべて０である。本実施形態では、数１２の行列を「終結式行列」と呼ぶ。

【0088】

【数12】

【0089】

そして、数１２の終結式行列から、右側２ｊ列と、係数ａ_０，…，ａ_ｍのついての下側ｊ行と、係数ｂ_０，…ｂ_ｎのついての下側ｊ行とを除いて、数１３のような小行列Ｍ_ｊ（ｆ，ｇ）を作る。ここで、ｐ＝ｍｉｎ｛ｍ，ｎ｝として、ｊ＝０，１，…，ｐ−１である。

【0090】

【数13】

【0091】

さらに、数１３の行列Ｍ_ｊ（ｆ，ｇ）について、右端の列の要素を、上から順にｆ（ｘ）ｘ^{ｎ−ｊ−１}，ｆ（ｘ）ｘ^{ｎ−ｊ−２}，…，ｆ（ｘ），ｇ（ｘ）^{ｍ−ｊ−１}，ｇ（ｘ）ｘ^{ｍ−ｊ−２}，…，ｇ（ｘ）で置き換えて、数１４のような行列Ｒ_ｊ（ｆ，ｇ）を作る。この行列の行列式が、ｆ（ｘ）とｇ（ｘ）のｊ次の部分終結式である。Ｒ_ｊ（ｆ，ｇ）は行列の要素に多項式を含むため、その行列式である部分終結式も多項式になる。本実施形態では、数１４の行列を「部分終結式行列」と呼ぶ。

【0092】

【数14】

【0093】

さて、スプライン算出部１３が算出したスプライン関数Ｓ_Ｆ（ｘ），Ｓ_Ｇ（ｘ）の、処理対象の区間における多項式ｆ（ｘ），ｇ（ｘ）が、ともにｍ次であるとする。この多項式ｆ（ｘ），ｇ（ｘ）について、Ψ_０（ｘ）＝ｆ（ｘ），Ψ_１（ｘ）＝ｇ（ｘ）とおく。そして、ｊ＝２，…，ｑについて数１５により多項式列｛Ψ_０（ｘ），Ψ_１（ｘ），…，Ψ_ｑ（ｘ）｝を作る。ただし、ｑはΨ_ｑが０次（定数）となる番号である。

【0094】

【数15】

【0095】

ここで、ｂ_ｍはｇ（ｘ）のｍ次の項の係数である。例えばΨ_０（ｘ）とΨ_１（ｘ）を３次とすると、この多項式列は、例えばΨ_２（ｘ）が２次、Ψ_３（ｘ）が１次、Ψ_４（ｘ）が０次（定数）のように、次数が順に下がる有限個の列となる（この例ではｑ＝４であり、ｑ＝ｍ＋１の関係がある）。

【0096】

すると、こうして得られる多項式Ψ_０（ｘ），Ψ_１（ｘ），…，Ψ_ｑ（ｘ）は、図６の処理により得られる対応する多項式ψ_０（ｘ），…，ψ_ｑ（ｘ）の正の定数倍になる。すなわち、ｃ_０，…，ｃ_ｑを正の定数として、Ψ_０（ｘ）＝ｃ_０ψ_０（ｘ），…，Ψ_ｑ（ｘ）＝ｃ_ｑψ_ｑ（ｘ）の関係がある。さらに、この多項式列｛Ψ_０（ｘ），Ψ_１（ｘ），…，Ψ_ｑ（ｘ）｝を求めるときはユークリッドの互除法を直接実行しないので、数値計算上の不安定性が生じない。

【0097】

本実施形態の信号処理装置１０は、多項式列から作られる数列の符号変化の回数に基づいてアンラップ位相を求めており、正の定数倍の違いはこの符号変化の回数を数える際に影響を与えない。したがって、第２の実施形態では、この多項式列｛Ψ_０（ｘ），Ψ_１（ｘ），…，Ψ_ｑ（ｘ）｝を算出することにより数値計算上の不安定性の問題を回避しつつ、連続的なアンラップ位相を決定する。

【0098】

さらに、第１の実施形態では多項式列を算出してから変数ｘに値を代入して数列を生成していたが、第２の実施形態では、部分終結式を計算する前に変数ｘに値を代入してしまえば、ｄｅｔ（Ｒ_ｊ（ｆ，ｇ））は数値の行列式になりその数列が直接求められる。数値の行列式は容易に計算することができるため、計算量の点でも、第２の実施形態の処理の方が第１の実施形態の処理より簡略化される。

【0099】

図４のステップ３０で第２の実施形態の多項式列算出部１５が多項式列を算出する処理についてまとめておく。図７は、第２の実施形態の多項式列算出部１５の動作例を示したフローチャートである。

【0100】

まず、現在の処理対象の区間［ｘ_ｋ−１，ｘ_ｋ］について、多項式列算出部１５は関数格納部１４からスプライン関数Ｓ_Ｆ，Ｓ_Ｇの３次関数ｆ_ｋ，ｇ_ｋを取得し、それぞれをΨ_０，Ψ_１とする（ステップ４１）。そして、変数ｊをｊ＝１とする（ステップ４２）。多項式列算出部１５はΨ_ｊの次数が０でないかどうかを判定し、次数が０でなければステップ４４に進み、０であれば後述するステップ４７に進む（ステップ４３）。

【0101】

Ψ_ｊの次数が０でない場合、多項式列算出部１５は、数１４により部分終結式行列Ｒ_ｊ＋１（Ψ_０，Ψ_１）を作成する（ステップ４４）。そして、数１５のΨ_ｊ＋１（ｘ）にアンラップ位相を決定するｘ＝ｔの値を代入した上で、数値の行列式Ψ_ｊ＋１（ｔ）を計算する（ステップ４５）。区間［ｘ_ｋ−１，ｘ_ｋ］内の別の点についてもアンラップ位相を求める場合は、ステップ４５の行列式の計算を繰り返してもよい。その後ｊに１を加算し（ステップ４６）、処理はステップ４３に戻る。一方、Ψ_ｊの次数が０である場合は、多項式列算出部１５が数列｛Ψ_０（ｔ），…，Ψ_ｊ（ｔ）｝を出力し、符号カウント部１７に渡す（ステップ４７）。以上で多項式列算出部１５の動作は終了する。

【0102】

なお、第１の実施形態と同様に、各多項式Ψ_ｊ（ｘ）がｘのべき乗を因数にもつ（すなわち、０次の項が０になる）場合、多項式列算出部１５は、そのべき乗を除いた部分をΨ_ｊ（ｘ）と定義し直して多項式列｛Ψ_０（ｘ），Ψ_１（ｘ），…，Ψ_ｑ（ｘ）｝を算出するとよい。

【0103】

このように、本実施形態の信号処理装置１０は、アンラップ位相を決定するための多項式列を算出するときに、部分終結式の方法を用いる。この方法では、ユークリッドの互除法を直接実行しないため、多項式列が数値的に安定に算出される。また、本実施形態の方法は、浮動小数点演算であっても、多倍長整数演算であっても、安定に実現できる。

【0104】

＜数値実験＞
以下では、第１の実施形態の信号処理装置１０を用いた位相アンラップの数値実験について説明する。この数値実験では、ｆ（ｘ）＝ｃｏｓ（Θ（ｘ）），ｇ（ｘ）＝ｓｉｎ（Θ（ｘ））でありΘ（ｘ）が以下の数１６で表される信号を一定の時間間隔にて観測し、観測した信号から元の信号ｆ（ｘ），ｇ（ｘ）の連続位相Θ（ｘ）を推定する。この数値実験では、観測した信号から、信号処理装置１０の位相アンラップ方法と別の位相アンラップ方法とにより、数１６の位相を正しく求められるかどうかを調べる。

【0105】

【数16】

【0106】

その別の方法としては、隣接する標本点で生じるアンラップ位相の変動が±π以内に収まっていることを仮定し、この仮定に基づき、標本点のラップ位相を２πの整数倍ずつ変化させて１つ前の標本点におけるアンラップ位相の値から±π以上離れない値を求め、その値をその標本点のアンラップ位相とする方法を用いる。以下では、この方法のことを「比較法」と呼ぶことにする。

【0107】

図８は、位相アンラップの数値実験に用いた一方の入力信号ｆ（ｘ）のグラフであり、図９は、その数値実験に用いた他方の入力信号ｇ（ｘ）のグラフである。図８（ａ），図８（ｂ），図８（ｃ）は、それぞれ標本間隔を０．２５，０．４，０．５５として信号を観測した場合のグラフであり、これは図９〜図１２でも共通である。それぞれのグラフの中に、標本点を○印で示している。また、それぞれのグラフにおいて、実線が真のｆ（ｘ），ｇ（ｘ）である。鎖線は、各標本点での信号を信号列として信号取得部１１が取得し、その信号列をもとにスプライン算出部１３が算出したスプライン関数を示している。この数値実験では、隣接する標本点間をそれぞれ３次関数で近似している。標本間隔が０．２５と小さい図８（ａ）と図９（ａ）では真のｆ（ｘ），ｇ（ｘ）とスプライン関数との差はほとんどわからないが、標本間隔が０．５５と大きくなる図８（ｃ）と図９（ｃ）では両者の差異がはっきり認められる。

【0108】

図１０は、図８および図９の入力信号のラップ位相のグラフであり、図１１は、図８および図９の入力信号のアンラップ位相のグラフである。図１１のアンラップ位相は数１６から得られる真のアンラップ位相であり、図１０のラップ位相はその値域を（−π，π］に制限して得られるラップ位相である。図１０のラップ位相は値が−πからπまたはπから−πに変化するところで不連続になっているが、図１１のアンラップ位相は連続関数である。なお、図１０（ａ）から図１０（ｃ）のグラフと図１１（ａ）から図１１（ｃ）のグラフは、数１６をプロットしたものであるため、標本間隔の大きさにかかわらずすべて同じである。

【0109】

図１２は、図８および図９の入力信号について位相アンラップを行った数値実験の結果を示したグラフである。実線は比較法により得られたΘ（ｘ）の推定値を表し、鎖線は信号処理装置１０により得られたΘ（ｘ）の推定値を表す。標本間隔が比較的小さい図１２（ａ）と図１２（ｂ）の場合には、信号処理装置１０の方法と比較法のどちらでも、グラフは図１１（ａ）および図１１（ｂ）のものとほとんど同じである。しかし、標本間隔が大きい図１２（ｃ）の場合には、図１１（ｃ）のグラフとの間に差異が見られる。

【0110】

比較法では、隣接する標本点で生じるアンラップ位相の変動が±π以内に収まっていることを仮定しているため、隣接標本点間の真のアンラップ位相が±π以上離れる場合には、正確に位相アンラップを行うことができない。このことは、図１２（ｃ）で実際に確かめられる。図１１（ｃ）を見ると、標本間隔が０．５５の場合は、隣接する２つの標本点であるｘ＝７．７とｘ＝８．２５におけるアンラップ位相Θ（ｘ）の差がπより大きいことがわかる（差は約１．０４πである）。そして、比較法によるアンラップ位相を表した図１２（ｃ）の実線は、その標本点間で図１１（ｃ）の曲線からずれている。このことから、比較法ではｘ＝８．２５における位相アンラップに失敗していることが確かめられる。

【0111】

一方、本実施形態の信号処理装置１０では、図１２（ｃ）の場合であっても、各標本点におけるアンラップ位相の推定値は図１１（ｃ）の真のアンラップ位相に完全に一致しており、全体的に正しく位相アンラップを実行できていることが確かめられる。すなわち、隣接標本点間における真のアンラップ位相の差がπより大きい場合であっても、正しく位相アンラップを行うことができる。このため、信号処理装置１０により出力される位相信号は、比較法による出力結果と異なり、隣接標本点間の位相差がπより大きいような標本点も含まれることがある。

【0112】

以上の数値実験の結果から、本実施形態の信号処理装置１０では、入力信号列の標本間隔を十分小さくとっておけば、入力信号列がスプライン関数によって十分な精度で近似され、位相アンラップが成功することがわかる。

【0113】

ところで、本実施形態の信号処理装置１０は、ＰＣ（Ｐｅｒｓｏｎａｌ Ｃｏｍｐｕｔｅｒ）等の汎用のコンピュータ内で実現するとよい。最後に、このような汎用のコンピュータのハードウェア構成について説明する。

【0114】

図１３は、コンピュータ９０のハードウェア構成を示した図である。図示するように、コンピュータ９０は、プロセッサ９１と、メインメモリ９２と、記憶装置９３と、通信インタフェース９４と、表示機構９５と、入力インタフェース９６とを備える。プロセッサ９１は、記憶装置９３に記憶されるプログラムを実行することにより、上述した各機能を実現する。メインメモリ９２は、プロセッサ９１が実行中のプログラムや、そのプログラムにより一時的に使用されるデータ等を記憶する。記憶装置９３は、プロセッサ９１が実行するプログラムやそのプログラムに関する入出力データ等を記憶する。通信インタフェース９４は、外部装置との間でデータの送受信を行う。表示機構９５は、ビデオメモリやディスプレイ等を含み、ユーザにデータ等を表示する。入力インタフェース９６は、キーボードやマウス等を含み、ユーザからの入力操作を受け付ける。

【0115】

なお、図３を用いて説明した信号処理装置１０の各機能部は、図１３に示すプロセッサ９１がプログラムを記憶装置９３からメインメモリ９２に読み込んで実行することにより実現される。また、関数格納部１４は、例えば図１３に示すメインメモリ９２により実現される。

【0116】

なお、本実施形態を実現するプログラムは、通信手段により提供してもよいし、ＣＤ−ＲＯＭ等の記録媒体に格納して提供してもよい。

【符号の説明】

【0117】

１０ 信号処理装置
１１ 信号取得部
１２ 区間決定部
１３ スプライン算出部
１４ 関数格納部
１５ 多項式列算出部
１６ 数列生成部
１７ 符号カウント部
１８ アンラップ処理部
１９ 信号出力部

【図1】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図2】