ブログ

numpy.fft.fftfreqで作る次のようなコードではfreq2は途中でいきなり負に飛ぶことになり，plt.plot(freq2, Fz_abs);とするとプロットが汚くなる．xlimで範囲制限をしても枠外の点とつながった線が気になる．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 1000);
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; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp - 1.0/(t_len + t[1] - t[0]), num=len(t));
freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);

plt.figure(1)
plt.plot(freq, Fz_abs);

plt.figure(2)
plt.plot(freq2, Fz_abs);
plt.xlim([0, 125])


plt.show();

たとえば，np.where(freq2 >= 0, freq2, np.nan)で選別してもよい．要素がNaNになっているとプロットに失敗するので描画されない（plt.figure(3)に示した）．
もしくは，np.where(freq2 >=0)でインデックス列の状態にしておいてそれでプロット範囲を周波数で制限することも可能でplt.figure(4)に示した．

import numpy as np
import matplotlib.pyplot as plt

# 例の信号を作成
t = np.linspace(0.001, 4, 1000);
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; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);

freq = np.linspace(0, fsmp - 1.0/(t_len + t[1] - t[0]), num=len(t));
freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);

plt.figure(1)
plt.plot(freq, Fz_abs);

plt.figure(2)
plt.plot(freq2, Fz_abs);
plt.xlim([0, 125])

freq3 = np.where(freq2 >= 0, freq2, np.nan)
plt.figure(3)
plt.plot(freq3, Fz_abs);

plt.figure(4)
plotrange = np.where(freq2 >= 0)
plt.plot(freq2[plotrange], Fz_abs[plotrange]);

plt.show();