質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.57%

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

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 1
  • VIEW 2,157

Fallout_18

score 106

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

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()
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 1

checkベストアンサー

+3

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にあたる場所)に期待どおりのスペクトルが出ていることがわかると思います。

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/10/21 15:25

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

    キャンセル

  • 2018/10/21 16:56

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

    キャンセル

  • 2018/10/21 17:51

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.57%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る