phase unwrapping

phase unwrappingという処理があります。
これは、振動関数から位相を抽出しよう、ということを実現する処理です。
具体的に問題を書けば、関数
\(\displaystyle
y(x)=\sin(\phi(x))
\)
があり、\(y(x)\)だけが分かっているときに、\(\phi(x)\)を求める、という問題です。
仮定として、\(\phi(x)\)は滑らかな関数であるとします。

逆関数を作用させれば良いということはすぐに思いつきます。つまり、\(\sin\)の逆関数である\(\sin^{-1}\)を作用させればよく、
\(\displaystyle
\begin{align}
\sin^{-1}y(x)&=\sin^{-1}\sin(\phi(x)) \\
&=\phi(x)
\end{align}
\)
を期待します。

しかしながら、ここで問題が生じます。
それは、\(\sin^{-1}(x)\)の値域は\([-\pi/2,\pi/2]\)ですので、この範囲に押し込まれてしまうということです。

もし\(\phi(x)\)が\([-\pi/2,\pi/2]\)の範囲以外にはない、という仮定があるならば話は別ですが、
通常は\(\phi(x)=[-\infty, \infty]\)ですので、\(\sin^{-1}(x)\)を作用させた場合、値域の端で位相に不連続が生じることになります。
つまり、ただ単に\(\sin^{-1}\)を作用させただけでは、\(\phi(x)\)が滑らかである、という条件に矛盾することになります。

しかしながら、\(\sin^{-1}\sin(\phi(x))\)と\(\phi(x)\)は全く無関係なもの、というわけではありません。
\(
\begin{align}
\sin(\theta)&=-\sin(-\theta) \\
\sin(\theta)&=\sin(\theta+2n\pi), ~~~(n=0,\pm 1,\pm 2,\cdots)
\end{align}
\)
の性質があるため、\(\sin^{-1}\cos(\phi(x))\)の角度を負に取ったり、適当に\(2n\pi\)を足したり引いたりすれば、どこかに元の\(\phi(x)\)と一致する組み合わせが見つかるはずです。

このように、関数の自由度を一つに決めて、制限された関数の値域を元に戻すことをunwrappingと言います。
位相の場合は特にphase unwrappingと呼ばれます。

先にどのようなことを考えるか図示してみます。
入力は下の図の黒線で、ここから赤線を復元できるか？という問題です。

どのように\(\phi(x)\)と一致する組み合わせ見つけるかは様々な方法が考えられますが、愚直な方法で実装しましょう。
ここでは一例として、\(\phi(x)\)は滑らかである場合を考えて、補外によってphase unwrappingを行う方法を実装してみます。

方針

滑らかに接続できたことを考えたとしても、絶対的な位相の量は分からないので、適当な\(x_0\)における位相の値\(\sin^{-1}\sin(\phi(x_0))\)を基準として、その位置から滑らかにつながることを想定します。

以下の仮定の下で考えます。

• 離散的な\(y(x_n)=\sin^{-1}\sin(\phi(x_n)), (n=1,2,\cdots,N)\)が等間隔で与えられている。
• 入力\(y(x_n)\)の値域は\([0,\pi]\)または\([-\pi/2,\pi/2]\)
• \(\phi(x)\)は滑らかであること
• \(\phi(x)\)が局所的には2次程度の多項式で書けること

まず、基準点より一つ進んだ点の候補(\(\sin(\theta)=-\sin(-\theta), \sin(\theta)=\sin(\theta+2n\pi)\)の関係から考えられるもの)を挙げ、その中から期待されるであろう補外の値に最も近い点がunwrappingされた点だ、と考えます。

unwrapされる候補の点を挙げる

左側の元の関数\(y(x_n)\)から、右のように候補を生成します。候補は、これまでにunwrapされた関数を値の周辺のみを候補とします。

補外と上手くいかない場合

基本的なアイデアは下図左側で、これまでにunwrapされた点を元に次の点を予測し、それに最も近い候補点がunwrapの値とします。この補外によって、\([-\pi/2,\pi/2]\)ないし\([0,\pi]\)から\([-\infty, \infty]\)の領域に拡張できます。

また、異なる候補と接近する境界付近では間違えて予測する可能性があります。それは右のような状況です。
これは補外が正しく行われないことが理由なので、逆向きにも次の点を予測し、その平均を取ることで精度が上がることを期待します。

コード

Fortran90のコードです。

失敗する場合

離散した点と点の間隔が広い場合

本アルゴリズムの場合、どうしても境界付近で曲がる場合、別の領域のものをトレースしてしまう可能性があります。
厳密に傾きがゼロになり、候補が複数ある場合に正解を選び出すのは不可能です。本アルゴリズムの適用を諦めましょう。
しかしながら、元のデータ点を補間なり何かをしてマシにすることはできますので、データ点を増やして再度実行すると成功するかもしれません。

具体的に例を挙げてみます。

関数
\(\displaystyle
\begin{align}
y(x)&=acos(cos(\phi(x))) \\
& phi(x)=3.1cos(2*x)**2+exp(-0.2x)-0.8
\end{align}
\)
を考えます（適当に試していた時に失敗した関数なので、これより単純な関数でも細かく離散化されていなければ失敗するはずです）。

離散化する間隔\(\Delta x\)を\(\Delta x=0.1\)でとり、unwrapすると失敗します。

上図の横軸は離散点のインデックス、縦軸は候補やその他の情報を示しています。\(cos(x)\)の任意性は無視しています。
数値的な解は緑線であり、9,10番目の離散点（\(x=0.9, 1.0\)）で失敗しており、本来の解（赤線）とは異なる点を結んでしまっています。

しかしながら、離散化する間隔\(\Delta x\)を\(\Delta x=0.05\)でとり、unwrapすると成功します。

成功させるためには、できるだけ細かい間隔で離散化しましょう。

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
２種類（1．初期値近傍の極値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法）
があるようですが、ここでは1. 初期値近傍の極値を求める方のニュートン法についてのお話です。

まとめ

ベクトル\(\mathbf{x}=(x_1, x_2, \cdots, x_N)\)で指定される実数値関数\(f(\mathbf{x})\)の極値を求める事を考えます。
ここで、初期値近傍の点\(\mathbf{x}=\mathbf{x}^{(0)}\)が分かっている時、その近傍にある\(f(\mathbf{x})\)の極小値を与える\(\mathbf{x}=\mathbf{a}\)は適当な反復計算を行い、
\begin{equation}
\mathbf{a}=\lim_{k\to\infty}\mathbf{x}^{(k)}
\end{equation}
で与えられます。極値を求めるニュートン法を用いる場合、\(\mathbf{x}^{(k)}\)は漸化式として与えられ、
\begin{align}
\mathbf{x}^{(k+1)}&=\mathbf{x}^{(k)}+\mathbf{d}^{(k)} \\
\mathbf{d}^{(k)}&=-\left[H f\bigl(\mathbf{x}^{(k)}\bigr)\right]^{-1}\cdot \left[\nabla f\bigl(\mathbf{x}^{(k)}\bigr)\right]
\end{align}
として逐次的に与えられます。ここで、\(H f(\mathbf{x}), \nabla f\bigl(\mathbf{x}\bigr)\)は関数\(f(\mathbf{x})\)のヘッセ行列、傾き(Gradient)を表し、
\begin{equation}
\hat{H}f=
\begin{pmatrix}
\frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_1 \partial x_n} \\
\frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} &\cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\
\vdots & \vdots &\ddots & \vdots \\
\frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} &\cdots &\frac{\partial^2 f}{\partial x_n^2} \\
\end{pmatrix}
,~~\nabla f=
\begin{pmatrix}
\frac{\partial f}{\partial x_1} \\
\frac{\partial f}{\partial x_2} \\
\vdots \\
\frac{\partial f}{\partial x_n}
\end{pmatrix}
\end{equation}
で与えられます。

特に\(f\)が\(x, y\)の２変数関数\(f(x,y)\)であるならば、
\begin{align}
\begin{pmatrix}
x^{(k+1)} \\
y^{(k+1)}
\end{pmatrix}
&=
\begin{pmatrix}
x^{(k)} \\
y^{(k)}
\end{pmatrix}
-\frac{1}{f_{xx}f_{yy}-f_{xy}f_{yx}}
\begin{pmatrix}
f_{yy} & -f_{xy} \\
-f_{yx} & f_{xx}
\end{pmatrix}
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{align}
ここで、
\begin{equation}
\hat{H}f
=
\begin{pmatrix}
f_{xx} & f_{xy} \\
f_{yx} & f_{yy}
\end{pmatrix}
,~~
\nabla f=
\begin{pmatrix}
f_x \\
f_y
\end{pmatrix}
\end{equation}
です。ここで、\(f_p\equiv\frac{\partial f}{\partial p}, ~f_{pq}\equiv\frac{\partial^2 f}{\partial p\partial q}, (p, q=x,y)\)と書きました。
微分を差分で近似するとします。\(x, y\)方向の刻み幅を\(h=h_x=h_y\)とすればそれぞれの偏微分は点\((x,y)\)近傍の7点を用いて
\begin{align}
f_x &= \frac{f(x+h, y) – f(x-h, y)}{2h} +O(h^2)\\
f_y &= \frac{f(x, y+h) – f(x, y-h)}{2h} +O(h^2)\\
f_{xx} &= \frac{f(x+h, y) – 2f(x, y)+ f(x-h, y)}{h^2} +O(h^2)\\
f_{xy} &=f_{yx}= \frac{1}{2h^2}\bigl[f(x+h, y) + f(x-h, y) + f(x, y+h) \\
&\hspace{1em}+ f(x, y-h)-2f(x,y)-f(x+h,y+h)-f(x-h,y-h)\bigr]+O(h^2) \\
f_{yy} &= \frac{f(x, y+h) – 2f(x, y)+ f(x, y-h)}{h^2} +O(h^2)
\end{align}
と書けます[2]。

導出

詳しくは述べませんが、1次元の問題における極小値を求めるニュートン法くらいは導いておきましょう。
多次元の場合は拡張すれば良く、考え方は変わりません。多次元では他の文献、例えば[1]を参考にして下さい。

関数\(f(x)\)を\(x=x^{(k)}\)周りでテイラー展開を行えば、
\begin{equation}
f(x)=f\bigl(x^{(k)}\bigr)+f’\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+\frac{1}{2}f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)^2+O(\Delta^3)
\end{equation}
となります。ここで、\(\Delta=|x-x^{(k)}|\)とおきました。
今、欲しい解\(a\)は関数\(f(x)\)の極値なので、
\(\frac{df}{dx}=0\)を満たす\(x\)を探せばよいのです。よって、両辺を\(x\)で微分して、
\begin{align}
\frac{df(x)}{dx}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x-x^{(k)}\bigr)+O(\Delta^2)
\end{align}
を得ます。これが\(0\)になる\(x\)が解\(a\)なので、
\begin{equation}
\left. \frac{df(x)}{dx}\right|_{x=a}=f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(a-x^{(k)}\bigr)+O(\Delta^2) = 0
\end{equation}
となります。今、\(O(\Delta^2)\)を無視すれば得られる解\(a\)は近似解\(a\approx x^{(k+1)}\)となるので、
\begin{equation}
f’\bigl(x^{(k)}\bigr)+f^{\prime\prime}\bigl(x^{(k)}\bigr)\bigl(x^{(k+1)}-x^{(k)}\bigr) = 0 \\
\end{equation}
を解いて
\begin{equation}
x^{(k+1)} = x^{(k)}-\frac{f’\bigl(x^{(k)}\bigr)}{f^{\prime\prime}\bigl(x^{(k)}\bigr)}
\end{equation}
となります。…ここまで来てしまうとニュートン・ラフソン法と同じですね。関数値がゼロになる点を探すのか、導関数値がゼロになる点を探すのかの違いなだけです。
多次元だと行列になってしまいシンプルではない形になり難しく見えますが、１次元では非常に単純です。

問題は1階(傾き)、2階導関数(ヘッセ行列)をどのように求めるかです。厳密な値が無い場合は近似的にしか得るしかありません。基本的にこれらの計算量は膨大になります。
最も簡単な方法は微分を差分に置き換えてしまうことですが、差分は余り良い方法とは言えません。

厳密なヘッセ行列は、例えば自動微分の考えで得ることが出来ます。例えば、Hyper-dual numbersによる二階偏微分の計算を参考にして得ることが出来ます。

また、引数が複素数で計算結果が実数で返ってくる関数\(f(z)\)の最小値を求めたいのであれば、\(z=x+iy\)として、\(f(x,y)\)の二変数関数としてその極致を探せば良いです。

2次元のプログラム

ここでは2次元の問題に限り、傾き、ヘッセ行列を差分で近似した場合のプログラムを記述します。

下のプログラムは関数
\begin{equation}
f(x,y)=x e^{-(x^2+y^2)/2}
\end{equation}
の、点\((x_0, y_0)=(-1.2, -0.3)\)近傍にある極値を求めるプログラムです。
グラフで示せば、以下の通りになります。
この関数の極値は黒点、初期値は赤点で示しています。

以下に実行例を示します。

1行目1,2,3列目はそれぞれ\(x,y\)の初期値, 関数\(f(x,y)\)の値を示しています。
また、2行目1,2,3,4列目は求まった\(x,y\)の値, 関数\(f(x,y)\)の値、そして収束に至るまでに行った繰り返しの回数を示しています。

参考文献

[1]塩浦昭義, 数理計画法 第13回
[2]Abramowitz and Stegun, “Handbook of Mathematical Functions”, http://people.math.sfu.ca/~cbm/aands/page_883.htm, http://people.math.sfu.ca/~cbm/aands/page_884.htm