フーリエ変換が上手くいかない

Question

以下のデータをフーリエ変換しています。

date=[0.015625, 0.015625, 0.0625, 0.03515625, 0.09765625, 0.048095703125, 0.10357666015625, 0.07080078125, 0.17388343811035156, 0.13126850128173828, 0.2402327060699463, 0.17284822463989258, 0.26214931160211563, 0.1704627461731434, 0.25162281142547727, 0.15652238787151873, 0.22783483436796814, 0.12990689981961623, 0.1822412261408317, 0.08294962953823415, 0.11317331792179175, 0.0389063490084709, 0.056472000307913106, 0.005247288913015069, 0.009530311218974852, 0.005017450124746858, 0.0015365243368495318, 0.019833988591996476, 0.004779612072076682, 0.03359424834958215, 0.01336109187288527, 0.07296953486201782, 0.07118835430567244, 0.20142962808299264, 0.16327255228856913, 0.2845359127984799, 0.19719090119404797, 0.2510062185645254, 0.13409349144758365, 0.21390037959972297, 0.14605373301204233, 0.22151710012483006, 0.13282687534961274, 0.21408555469675444, 0.11668752550020667, 0.14086572084977306, 0.04890294330442227, 0.0752691630432859, 0.02164490058919605, 0.039182745264478534]

図のように周期的になりますが、途中で切れてしまっているので、窓関数blackmanに通してから、フーリエ変換しましたが、0と50辺りに飛びがでてしまい、また周期的な図と比べても来て欲しい場所の値（26辺り？）にピークがきません。

from scipy import signal
itr = 50
w = signal.blackman(itr)
#plt.figure()
windowdate = date[0:itr]*w[0:itr]
#plt.plot(windowdate)
f = np.fft.fft(windowdate)
f_abs =np.abs(f)
plt.xticks(np.arange(0,itr,10))
plt.bar(time,f_abs)
plt.show()

Answer 1

fft自体の計算や結果の周波数成分はあっていると思いますが、それをプロットするにあたり横軸の次元を間違って扱っておられるようです。

• (A)元の波形の横軸

次元は「時間」です。コードではサンプルリングデータの単位を単純に用いてプロットしておられますが、それをそのまま受け入れるなら波形の一周期は25単位時間ほどであるといえます。

• (B)FFT結果の横軸

次元は「周波数」です。周波数は「単位時間あたりに何周期するか」ということですので、一周期が25単位時間であるならその波の周波数は1/25=0.04ですね。質問者さんは26あたりにスペクトルのピークが出ることを期待しておられますが正しくは周波数1/26=0.038あたりにピークが来ると考えるべきです。

具体的なコードの問題を挙げますと「本来はFFT結果の横軸に周波数を指定すべきところ、元波形の横軸(time)をそのまま設定している」ということになります。

結果のスペクトルの棒グラフをプロットする際に注意すべき点としては

• 横軸の値は基本周波数(※1)の倍数の数列にする
• 基本周波数が1より小さいとき、棒グラフの幅が広すぎて見づらくなるためwidthに適当な値を指定する

くらいかと思います。検証のため、サンプリング周波数=1、サンプル数=50として周期25の正弦波を用いてスペクトルをプロットするコードを書いてみました。なお、本件に直接関係ないですが、matplotlibの関数的インターフェース(plt.xxxを呼び出す)ではなくオブジェクト指向的インターフェース(figureやaxesに対してメソッドを呼び出すスタイル)で記述しています。わかりにくかったらすみません。

import matplotlib.pyplot as plt
import numpy as np

N = 50       # サンプル数
Fs = 1       # サンプリング周波数
Dt = 1 / Fs  # サンプル間の時間差
Fb = Fs / N  # 基本周波数

# 元波形の横軸
ts = np.linspace(0, Dt * N, N, endpoint=False)

# テスト用の波形
Ft = 0.04  # 周波数(周期=25)

data = np.sin([2 * np.pi * Ft * t for t in ts])

# 元波形のプロット
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(ts, data)

# FFT
spectrum = np.fft.fft(data)
abs_spectrum = np.abs(spectrum)

# 周波数成分のプロット
fs = np.linspace(0, Fb * N, N, endpoint=False)
ax2 = fig.add_subplot(212)
ax2.bar(fs, abs_spectrum, width=Fb * 0.8)
# tick間隔をFtにしたほうが結果がよりはっきりわかるのですが、そうしてしまうと
# 横軸ラベルが重なってしまい視認できなかったのでFtの倍の間隔にしてあります
ax2.xaxis.set_ticks(np.arange(0, 1, Ft * 2))

plt.show()

※1:基本周波数
離散フーリエ変換では、サンプリング周波数Fs、サンプル数Nのとき、計算結果のスペクトルは基本周波数(=Fs/N)の倍数の周波数のものが出てきます。

結果

この結果を見ると、質問者さんの棒グラフでいえば横軸の値が2のところ(実際は周波数0.04にあたる場所)に期待どおりのスペクトルが出ていることがわかると思います。

Comments

[Fallout_18]

ありがとうございます、一度こちらを参考にチャレンジしてみます。また他のサイトに書いてあったのですが、0.96のところは無視して良いものですよね。

[Fallout_18]

それと少し、話がずれてしまいますが、この窓関数を通してフーリエ変換してもきれいにいかずに、余計なピークがⅩ軸の端にでてしまうのは何故なのでしょう。。。。？

[KSwordOfHaste]

> 0.96のところの値を無視
DFTの結果はサンプリング周波数の1/2(ナイキスト周波数)を中心に左右対称になります。スペクトルの特徴を観察する目的ではグラフの左半分だけに注目すればよいです。
> 余計なピークがX軸の端に
それは直流成分であり「窓関数」とは無関係で元信号の平均が0付近にないから出てくるのです。正弦波を見ると信号の平均（直流成分）は０になっていますね。だからDFT結果の周波数0.00の成分が０となるのです。

[Fallout_18]

なるほど、それについては理解しました、ありがとうございます。
質問を出し惜しみしてすいません、上の質問のデータ（周期を１００００にしてフーリエ変換した場合、どうやら0.04辺りにピークは来ません、この質問はは新たに設けた方が良いですか？それともこちらでもよろしいでしょうか）、なんどもすいません。

[Fallout_18]

すいません、ちょっともう一回確認してみます。

[Fallout_18]

ただ線が細くて、ピークが見えないだけでした笑
りかいしました、ありがとうございます！