ブログ

numpy.fft.fftfreqによって生成した周波数軸とnumpy.whereを使うと周波数でフィルタリングが行える．
例えば，
Fzf = np.fft.fft(z);
Fz2f = np.where(np.abs(freq2) > 12, 0, Fzf)
などとすれば，周波数成分を取り除くことができる（０で埋めることができる）．
np.fft.ifftで時間領域の信号に戻せば0で埋めた成分は消えている．

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();

Fzf = np.fft.fft(z); # Fzfはfftした結果
Fz = Fzf / z.shape[0] * 2; # 折り返すのでパワーが２分の１になっている．
Fz[0] = Fz[0] / 2; # 平均成分は折り返さない．
Fz_abs = np.abs(Fz);


freq2 = np.fft.fftfreq(len(t), d=1.0/fsmp);# 周波数軸の生成

# abs(freq2)が12以上の成分をカット
Fz2f = np.where(np.abs(freq2) > 12, 0, Fzf) # filter
z_filtered = np.fft.ifft(Fz2f);
Fz2 = Fz2f / z.shape[0] * 2;
Fz2[0] = Fz2[0] / 2;
Fz2_abs = np.abs(Fz2);


plt.figure(1)
# 比較FFTプロット
plt.subplot(211);
plotrange = np.where(freq2 >= 0)
plt.plot(freq2[plotrange], Fz_abs[plotrange]);
plt.plot(freq2[plotrange], Fz2_abs[plotrange]);

# 比較波形プロット
plt.subplot(212);
plt.plot(t, z);
plt.plot(t, z_filtered);

plt.savefig("valuestar.work.177.png")

plt.show();