Python numpyでFFTを実行する.周波数軸はどうしたら
numpyでFFTを実行する場合,周波数軸はnumpy.linspaceで設定するのが便利である.
よーく見るとずれてしまって困った話.
まず周波数解析の例を示す.周波数軸を作る部分は
np.linspace(0, サンプリング周波数, データ長さ)
で大体あっているのだが,微妙にずれる.
このずれの分が気になったのでよーく考えると以下の例の場合には
freq = np.linspace(0, fsmp - 1/(t_len + t[1] - t[0]), num=len(t));
ということになる.
ピったしの周波数軸が作れる.
import numpy as np import matplotlib.pyplot as plt # 例の信号を作成 t = np.linspace(0.001, 4, 4000); z = 0.1 + 0.2 * np.sin(t * 10 * 2 * np.pi) + 0.2 * np.sin(t * 33 * 2 * np.pi); # サンプリング周波数 fsmp = 1 / (t[1] - t[0]); # 解析時間 t_len = t.max() - t.min(); Fz = np.fft.fft(z) / z.shape[0] * 2; # 折り返すのでパワーが2分の1になっている. Fz[0] = Fz[0] / 2; # 平均成分は折り返さない. Fz_abs = np.abs(Fz); freq = np.linspace(0, fsmp, num=len(t)); fname = 'text.csv'; np.savetxt(fname, np.array([freq, np.real(Fz), np.imag(Fz)]).transpose(), fmt='%.18g', delimiter=','); plt.figure(1) plt.subplot(311) plt.plot(t, z) plt.subplot(312) plt.plot(freq, Fz_abs) plt.xlim([0, fsmp/2]); #136/4 plt.subplot(313) plt.plot(range(len(t)), Fz_abs) #plt.xlim([0, 500]); # 後半(高周波側)は前半の折り返し plt.show();
タグ:Python