位相アンラッピング(Phase unwrapping)

1D Phase Unwappingをサクッと行う

1.元データの作成
まず元データを作成。適当なデータを作る。

In [104]:

from numpy import *
N = 500

In [105]:

x = linspace(0, 2*pi, N)

In [106]:

import matplotlib.pyplot as plt
%matplotlib inline
y = 6*sin(x)*x
plt.plot(x, y)
plt.show()

2.関数をWrapする　
ラッピングされたデータを作る

In [107]:

yw = arctan2(sin(y), cos(y))
plt.plot(x, yw)

Out[107]:

[<matplotlib.lines.Line2D at 0x10bd8cc50>]

[-pi, pi]の範囲で折り返された(wrapされた)データが出来た。

3.Unwrap
ここまででデータ作成終了。連続だったデータもラップ(wrap)されて位相のジャンプが現れる。 これをUnwrapする。

In [108]:

def mywrap(y):
    for i in range(1, len(y)):
        div = y[i] - y[i-1]
        if div < -pi:
            y[i:] += 2*pi
        elif div > pi:
            y[i:] -= 2*pi
        else:
            pass
    
    return y
mywrap(yw)
plt.plot(x, yw)
plt.show()

元のデータが再現された。 でもこれはサンプリングが十分でかつノイズがない場合に限られる。
実はnumpyにはunwrapが準備されている。wrapメソッドはないようだ。

In [109]:

yw = unwrap(yw)
plt.plot(x, yw)
plt.show()

と同じ結果になる。

5.ノイズののったデータの場合

In [112]:

y = 6*sin(x)*x + random.normal(0, 1, len(x))
plt.plot(x,y)
plt.show()

In [113]:

yw = arctan2(sin(y), cos(y))
yw = unwrap(yw)
plt.plot(x, yw)
plt.show()

というふうに関数の値自体に対しては小さいノイズにも関わらず位相ジャンプが発生してしまう。
次にサンプリングを少なくしてみる。

In [127]:

N=40
x = linspace(0, 2*pi, N)
y = 6*sin(x)*x 
plt.plot(x, y, "o-")
plt.show()

In [132]:

yw = arctan2(sin(y), cos(y))
plt.plot(x, yw, "-o")

Out[132]:

[<matplotlib.lines.Line2D at 0x10bdbd190>]

In [133]:

yw = unwrap(yw)
plt.plot(x, yw, "o-")
plt.show()

序盤は再現できているが変化率の大きいところではまったく反対方向に増加してしまっている。 このようにノイズやサンプリング数によってunwrapは影響を受けてしまう。
以上、ちょっとした実験でした。