Pythonでフーリエ変換をする方法
解決済
回答 1
投稿
- 評価
- クリップ 0
- VIEW 4,356
前提・実現したいこと
Pythonでデジタル信号のフーリエ変換に試みています。
発生している問題・エラーメッセージ
https://org-technology.com/posts/smoother.html
↑こちらに示されているサンプルを元に自分のデータで解析しようと考えているのですが、「周波数スケール」や「データ数」を自分のデータの場合、どのように置き換えればよいのかわからないでいます…。
また、フーリエ変換をする際に
yf = fftpack.fft(y)/(n/2)
で、なぜn/2で割るのかがよくわかりません。
y2 = np.real(fftpack.ifft(yf2)*n)
では逆変換では逆にnをかけているのも少し疑問です。
ちなみに一般的な数学や物理の参考書には目を通しており、基本的なフーリエ変換の原理は理解しているつもりです。
該当のソースコード
import os
import time
import numpy as np
from scipy.stats import norm
from scipy import signal
import matplotlib.pyplot as plt
os.chdir(r"C:\Users\genuine\Desktop\earth science\20180217~misasa\20180224b")
files = os.listdir(r"C:\Users\genuine\Desktop\earth science\20180217~misasa\20180224b")
t1 = time.time()
dat_list = []
for file in files:
if file.count('.dat'):
dat_list.append(file)
ampl_list = []
for file in dat_list:
data = np.loadtxt(file, delimiter=' ')
x_value = data[:,0] #1列目を横軸とする
y_value = data[:,1] #2列目を縦軸とする
yh_value = signal.hilbert(y_value)
count = abs(yh_value)
sum = 0
for i in count[12000:13600]:
sum = sum + i
mean = sum/1600
ampl_list.append(mean)
x_axis = np.arange(800, 900.1, 0.2)
y_axis = np.array(ampl_list)
plt.plot(x_axis, y_axis, marker="", linewidth=0.2, color='black')
plt.show()
np.save('ampl_list.npy', y_axis) #.npyで保存する
t2 = time.time()
proceeding_time = t2 - t1
print('処理時間'+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日5回まで評価を下げられます
- 1日に1ユーザに対して2回まで評価を下げられます
質問の評価を下げる
teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。
- プログラミングに関係のない質問
- やってほしいことだけを記載した丸投げの質問
- 問題・課題が含まれていない質問
- 意図的に内容が抹消された質問
- 過去に投稿した質問と同じ内容の質問
- 広告と受け取られるような投稿
評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。
質問の評価を下げたことを取り消します
この機能は開放されていません
評価を下げる条件を満たしてません
質問の評価を下げる機能の利用条件
この機能を利用するためには、以下の事項を行う必要があります。
- 質問回答など一定の行動
-
メールアドレスの認証
メールアドレスの認証
-
質問評価に関するヘルプページの閲覧
質問評価に関するヘルプページの閲覧
checkベストアンサー
+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をかけているのも少し疑問です。
これは元の記事の次の操作が原因です。
yf2[(freq > fs)] = 0
yf2[(freq < 0)] = 0折り返し信号の分を残さないと行けないところを、カットオフ周波数以上のところを全て0にしてしまっています。これでは元の振幅成分のうちの半分しか復元されないため、n/2を掛けて補正するところを、(n/2)*2=nを掛けているわけです。
さて、本来やりたかったことですが、カットオフ周波数を決めなければ行けません。そのためにはサンプリング周波数もしくはサンプリング間隔が必要です。それがわかればfftfreqを使えば、その周波数がどこにあるか分かるはずです。
投稿
-
回答の評価を上げる
以下のような回答は評価を上げましょう
- 正しい回答
- わかりやすい回答
- ためになる回答
評価が高い回答ほどページの上位に表示されます。
-
回答の評価を下げる
下記のような回答は推奨されていません。
- 間違っている回答
- 質問の回答になっていない投稿
- スパムや攻撃的な表現を用いた投稿
評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。
15分調べてもわからないことは、teratailで質問しよう！
- ただいまの回答率 88.33%
- 質問をまとめることで、思考を整理して素早く解決
- テンプレート機能で、簡単に質問をまとめられる
2018/05/01 15:38
とりあえず、いただいたアドバイスを元にやってみようと思います。
2018/05/01 15:41