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

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

新規登録して質問してみよう
ただいま回答率
86.12%
Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

解決済

Pythonでフーリエ変換をする方法

genuine
genuine

総合スコア20

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

1回答

0グッド

0クリップ

7256閲覧

投稿2018/04/28 04:59

前提・実現したいこと

Pythonでデジタル信号のフーリエ変換に試みています。

発生している問題・エラーメッセージ

https://org-technology.com/posts/smoother.html
↑こちらに示されているサンプルを元に自分のデータで解析しようと考えているのですが、「周波数スケール」や「データ数」を自分のデータの場合、どのように置き換えればよいのかわからないでいます…。
また、フーリエ変換をする際に
yf = fftpack.fft(y)/(n/2)
で、なぜn/2で割るのかがよくわかりません。
y2 = np.real(fftpack.ifft(yf2)*n)
では逆変換では逆にnをかけているのも少し疑問です。

ちなみに一般的な数学や物理の参考書には目を通しており、基本的なフーリエ変換の原理は理解しているつもりです。

該当のソースコード

python

1import os 2import time 3import numpy as np 4from scipy.stats import norm 5from scipy import signal 6import matplotlib.pyplot as plt 7 8os.chdir(r"C:\Users\genuine\Desktop\earth science\20180217~misasa\20180224b") 9files = os.listdir(r"C:\Users\genuine\Desktop\earth science\20180217~misasa\20180224b") 10 11t1 = time.time() 12 13dat_list = [] 14for file in files: 15 if file.count('.dat'): 16 dat_list.append(file) 17 18ampl_list = [] 19for file in dat_list: 20 data = np.loadtxt(file, delimiter=' ') 21 x_value = data[:,0] #1列目を横軸とする 22 y_value = data[:,1] #2列目を縦軸とする 23 yh_value = signal.hilbert(y_value) 24 count = abs(yh_value) 25 26 sum = 0 27 for i in count[12000:13600]: 28 sum = sum + i 29 mean = sum/1600 30 ampl_list.append(mean) 31 32x_axis = np.arange(800, 900.1, 0.2) 33y_axis = np.array(ampl_list) 34plt.plot(x_axis, y_axis, marker="", linewidth=0.2, color='black') 35plt.show() 36 37np.save('ampl_list.npy', y_axis) #.npyで保存する 38 39t2 = time.time() 40proceeding_time = t2 - t1 41print('処理時間'+str(proceeding_time))

イメージ説明

試したこと

ここに問題に対して試したことを記載してください。

補足情報(FW/ツールのバージョンなど)

x_axis = np.arange(800, 900.1, 0.2)
からお分かりかと思いますが、私が解析を試みているのは501個のデータ点からなるグラフのフーリエ変換です。

python 3.6.1
numpy 1.12.1
scipy 0.19.0
matplotlib 2.0.2

以下のような質問にはグッドを送りましょう

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

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

気になる質問をクリップする

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

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

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

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

適切な質問に修正を依頼しましょう。

回答1

3

ベストアンサー

yf = fftpack.fft(y)/(n/2)
で、なぜn/2で割るのかがよくわかりません。

n/2のうちのnはデータ点数で、同じ区間でも解像度が高くなるほど積分値が大きくなってしまうことへの補正です。2は折り返し信号の分をサボったためです。fftshiftを使ったとすれば、正の周波数の負の周波数ではx=0のy軸に対して対称になっていると思います。もう少し細かく見ると、cos項(複素フーリエ係数の実数部分)は正負で同じ値なのに対して、sin項(複素フーリエ係数の虚数部分)は符号が反転していると思います。

少し脱線しましたが、本来ならば同じ周波数の振幅成分が正負に分かれて存在したはずのところを、参考にしている記事では正の部分しか扱っていません。コレ自体はよく行われることです。ただし、本来ならばデータ点数nで割っておけばいいだけのところを、負の周波数成分が無視しているために振幅が合わなくなってしまいます。この点を考慮してnではなくn/2で割っているのです。

y2 = np.real(fftpack.ifft(yf2)*n)
では逆変換では逆にnをかけているのも少し疑問です。

これは元の記事の次の操作が原因です。

python

1yf2[(freq > fs)] = 0 2yf2[(freq < 0)] = 0

折り返し信号の分を残さないと行けないところを、カットオフ周波数以上のところを全て0にしてしまっています。これでは元の振幅成分のうちの半分しか復元されないため、n/2を掛けて補正するところを、(n/2)*2=nを掛けているわけです。


さて、本来やりたかったことですが、カットオフ周波数を決めなければ行けません。そのためにはサンプリング周波数もしくはサンプリング間隔が必要です。それがわかればfftfreqを使えば、その周波数がどこにあるか分かるはずです。

投稿2018/04/28 12:13

tachikoma

総合スコア3597

R.Shigemori, hayataka2049, KSwordOfHaste👍を押しています

良いと思った回答にはグッドを送りましょう。
グッドが多くついた回答ほどページの上位に表示されるので、他の人が素晴らしい回答を見つけやすくなります。

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

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

このような回答には修正を依頼しましょう。

回答へのコメント

genuine

2018/05/01 06:38

ありがとうございます。 とりあえず、いただいたアドバイスを元にやってみようと思います。
tachikoma

2018/05/01 06:41

ほいさ。上手くいかなかったらまた遠慮なく聞いてください。

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

ただいまの回答率
86.12%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問

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

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

NumPy

NumPyはPythonのプログラミング言語の科学的と数学的なコンピューティングに関する拡張モジュールです。

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。