2017-01-31 5 views
3

私は時系列から周波数成分を削除する方法は?

enter image description here

以下の時系列を持っている私は(時系列の通常の前処理の一部として)それらを除去するために、サイクルをチェックしたいので、私はFFTを適用します。

# Number of samplepoints 
N = len(y) 
# sample spacing 
T = 1.0 # 1 day 
x = np.linspace(0.0, N*T, N) 
yf = scipy.fftpack.fft(y) 
xf = np.linspace(0.0, 1.0/(2.0*T), N/2) 
components = 2.0/N * np.abs(yf[:N//2]) 

fig, ax = plt.subplots(1, 1, figsize=(10, 5)) 
ax.plot(xf, components) 

この結果、以下のプロットが得られます。

enter image description here

私は4つの、最大のコンポーネントを削除します。これを行うために、私は以下の式を実装しています。

enter image description here

max_components = sorted(components, reverse=True)[:4] 
idx_max_comp = [] 

for comp in max_components: 
    for i in range(len(components)): 
     if components[i] == comp: 
      idx_max_comp.append(i) 
      break 

cycle_signal = np.zeros(len(y)) 
for idx in idxs: 
    a, b = (2.0/N) * np.real(yf[idx]), (2.0/N) * np.imag(yf[idx]) 
    fi = xf[idx] 
    cycle_signal += (a * np.cos(2 * np.pi * fi * x)) + (b * np.sin(2 * np.pi * fi * x)) 

y = y - cycle_signal 

しかし、私は再びFFTを適用するとき、それは動作しませんでしたが簡単にわかります。

enter image description here

なぜ?

+0

statsmodelモジュールを使用しようとしましたか?それはここに示されているようにhttp://stackoverflow.com/questions/20672236/time-series-decomposition-function-in-python? –

+0

別の投稿で見つかったこの素晴らしい例を確認してください: http://stackoverflow.com/questions/36968418/python-designing-a-time-series-filter-after-fourier-analysis –

+0

FFT配列 'yfピークに対応する値をゼロに設定する。これらのピークは単一点ではないので、値の小さな範囲をゼロに設定する必要があります。その後、逆FFTを実行するだけで、目的の結果が得られます。 –

答えて

0

私は問題は次のようだと思う。

T = 1.0 # 1 day

サンプリング周波数を使用すると、1つのサンプル日を持っている場合は、あなたのサンプリング周波数をf =(1/24である1秒あたりのサンプル数として定義されます* 60 * 60)は約11.57407 uHz(マイクロヘルツ)で、ナイキスト周波数は5.787035 uHzになります。約2日間です。つまり、2日に1回より頻繁にサイクルの発生を確認することはできません。

+0

あなたの答えをより簡単に記述できますか? –

関連する問題