仕事でFFTを使うのですが、波形の遅れがFFTの位相で理論通り出力されるかを確認したくなりました。正弦波の位相を適当にずらしてFFTするという方法で検証したのですが、位相の出力がうまくいかなったのでその備忘録となります。FFTの使い方の参考にもなるかもしれません。
検証方法
- Pythonで基準となる正弦波と、位相をずらした正弦波を生成する。
- これらのsin波にFFTし、振幅と位相スペクトルを求める。
- 位相スペクトルの変化を比較する。
検証に使用したソースコード
以下のようなソースコードを作成しました。
最初に2つの5Hzの正弦波を作成しています。正弦波のうちの１つは位相がずれていないもの、もう一方は時間軸方向の45度ずらしたものです。それらの正弦波にFFTし振幅と位相スペクトルを求めています。
import matplotlib.pyplot as plt
import numpy as np
import cmath
N = 100 #データ点数
dt = 1/N #Δt(sec)
time = np.arange(0, 1, dt) #時間軸(0-1)
f= 5 #Hz
#sin波の生成
wfm = np.sin(time*2*np.pi*f)
shift_wfm = np.sin(time*2*np.pi*f - 45*np.pi/360) #時間軸方向に45度遅らせる
# シフト前の波形wfmのFFT処理
freq = np.fft.fftfreq(N,dt) #FFT後の周波数軸を求める。
comp_spect = np.fft.fft(wfm) #FFTの結果を保存
comp_spect = comp_spect / (N/2)
amp_spect = np.abs(comp_spect) #振幅スぺクトルを計算
phase_spect = np.angle(comp_spect) * 180 / np.pi #位相スペクトルを計算
#シフト後の波形shift_wfmのFFT処理
shift_comp_spect = np.fft.fft(shift_wfm)
shift_comp_spect = shift_comp_spect / (N/2)
shift_amp_spect = np.abs(shift_comp_spect) #振幅スぺクトルを計算
shift_phase_spect = np.angle(shift_comp_spect) * 180 / np.pi #位相スペクトルを計算
#グラフ化部分
ax1 = fig.add_subplot(3,1,1)
ax2 = fig.add_subplot(3,1,2)
ax3 = fig.add_subplot(3,1,3)
ax1.plot(time,wfm)
ax1.plot(time,shift_wfm)
ax1.set_xlabel("time(s)")
ax1.set_ylabel("Amplitude")
ax2.plot(freq[:N//2], amp_spect[:N//2])
ax2.plot(freq[:N//2], shift_amp_spect[:N//2])
ax2.set_xlabel("Frequency(Hz)")
ax2.set_ylabel("Amplitude")
ax3.plot(freq[:N//2], phase_spect[:N//2])
ax3.plot(freq[:N//2], shift_phase_spect[:N//2])
ax3.set_xlabel("Frequency(Hz)")
ax3.set_ylabel("Phase (deg)")
plt.subplots_adjust(hspace = 0.7)
plt.rcParams["font.size"] = 5
plt.show()
実行結果
上から順に、波形・振幅スペクトル・位相スペクトルのグラフとなります。オレンジの正弦波は、時間的に少し(90度)後ろにずれています。
振幅スペクトルは、狙い通り5Hzで振幅1.0となっており、正しい結果が得られています。一方で、位相スペクトルは全体的にばらついてしまっています。
本来であれば、位相スペクトルも振幅スペクトルと同様に、5Hzの成分のみが値を持ち、それ以外はゼロになってほしいところですがそうなっておりません。
なぜこうなるのか
調べたところ、位相の計算に使用しているnumpy.angle関数は、atan2(Imag, Real)によって計算された値を出力することが分かりました。
単一周波数の正弦波をFFTすると、実部と虚部に微小なノイズのような値が含まれることがあります。そのため、これらの非常に小さな値に基づいてatan2が計算されることで、位相が不安定になり、期待と異なる値が得られてしまうのではないかと考えました。
位相スペクトル計算の修正
位相の計算部分を以下のように修正しました。
虚部と実部の値が閾値以下であれば、0に強制します。その後、atan2(Imag,Real)を使って位相を求めるようにします。ちなみにatan2(0,0)=0となります。
#実部と虚部が閾値以下だったらゼロに矯正する。
for i in range(0,N//2):
#シフト前の波形のスペクトル
if np.abs(comp_spect.real[i]) < 1E-10:
comp_spect.real[i] = 0.0
if np.abs(comp_spect.imag[i]) < 1E-10:
comp_spect.imag[i] = 0.0
#シフト後の波形のスペクトル
if np.abs(shift_comp_spect.real[i]) < 1E-10:
shift_comp_spect.real[i] = 0.0
if np.abs(shift_comp_spect.imag[i]) < 1E-10:
shift_comp_spect.imag[i] = 0.0
phase_spect = np.angle(comp_spect) * 180 / np.pi #位相スペクトルを計算
shift_phase_spect = np.angle(shift_comp_spect) * 180 / np.pi #位相スペクトルを計算
実行結果は下記の通りとなります。期待していた通り、位相スペクトルも5Hzの成分のみが値を持ち、それ以外はゼロになりました。
FFTの位相の出力について
出力された位相スペクトルを見ると、青色に比べオレンジの方が、位相が負に大きくなっています。このことから、時間軸に対して波形が遅くなる（右側にずれる）と位相は負の方向に大きくなることがわかりました。数学的に正弦波\(\sin(\omega t)\)を右にずらすということは、\(\sin(\omega t-\phi)\)の操作になります。\(-\phi\)の値をFFTで出力するため、時間的に右に波形をずらすと、位相が負方向に大きくなるということが改めて確認できました。
また、FFTの位相は余弦波\(\cos\)を基準にした位相である。青色の正弦波の位相が-90度で出力されているのは、\(\sin(\omega t)=\cos(\omega t – 90^\circ)\)の-90度を反映しているからになります。