(python) 位相アンラップしてみた

ちょっと便利すぎやしないかい？python。
10年くらい前にC++で書いたオレオレ信号処理ライブラリが完全に陳腐化したわ。

sin関数の位相にガウス関数を加えたもの。

これから

を抽出するのが目標。

FFTして負の周波数成分をカット＝ヒルベルト変換

IFFTして、元の関数の虚数部を復元（解析信号の取得）

解析信号の位相。±πにラップされている。

アンラップした位相。ほぼ傾きの直線だが、
少しもっこりしている。

傾きを引く。
ピークがx=5のガウス関数が得られている。

以下ソース。

import numpy as np
from scipy.fftpack import fft
from scipy.fftpack import ifft
from scipy.fftpack import fftfreq
from scipy.fftpack import fftshift
import matplotlib.pyplot as plt
import pandas as pd
import math

#データ数
N = 1000

#サンプリング周期(sec)
dt = 0.01

#周波数  
f = 1.0

#時間軸       
t = np.linspace(1,N,N) * dt - dt

#時系列データ
#y = np.sin(2 * np.pi * f * t)

#f(x)を自分で定義
#ガウス関数
fx = lambda x: (math.exp(-(x-5)**2/2)) /  math.sqrt(2*math.pi) * 10

y = []
for i in range(len(t)):
    y.append(np.sin(2 * np.pi * f * t[i] + fx(t[i])))
    #y.append(fx(t[i]))

#周波数軸(MATLAB等で言うfftshiftがかかっている)
freq = fftfreq(N,dt)

#FFT
yf = fft(y)

#解析信号を得る
yf[int(N/2):N] = 0.0 + 0.0j
af = ifft(yf)

#unwrap
unwrap = np.unwrap(np.angle(af))

phase = []
for i in range(len(t)):
    phase.append(unwrap[i] - 2*np.pi * f * t[i]) #オフセットを引く

plt.figure(1)
plt.xlabel("time(sec)")
plt.ylabel("value(a.u.)")
plt.plot(t,y)

plt.figure(2)
plt.xlabel("freq(Hz)")
plt.ylabel("amplitude(a.u.)")
plt.plot(freq,np.abs(yf))

plt.figure(3)
plt.xlabel("time(sec)")
plt.ylabel("value(a.u.)")
plt.plot(t,np.real(af))
plt.plot(t,np.imag(af))

plt.figure(4)
plt.xlabel("time(sec)")
plt.ylabel("phase(rad)")
plt.plot(t,np.angle(af))

plt.figure(5)
plt.xlabel("time(sec)")
plt.ylabel("phase(rad)")
plt.plot(t,unwrap)

plt.figure(6)
plt.xlabel("time(sec)")
plt.ylabel("phase(rad)")
plt.plot(t,phase)

plt.show()

python(scipy?)にunwrap関数なんてあるのね。
今回はFFTして負の周波数成分をカットして・・・という方法で実装しているけど、
hilbert関数というのがあるので、そちらを使った方が早いかも知れない。
やってることは同じはずだけども。