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

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

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

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

NumPy

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

Python 3.x

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

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

1回答

6442閲覧

python パワースペクトルを表示したい

plato

総合スコア44

Matplotlib

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

NumPy

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

Python 3.x

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

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2019/01/23 15:11

前提・実現したいこと

  • 時系列データをパワースペクトル密度として図に表示したい

(リンク先の一番下の図ようなものです)
http://hclab.sakura.ne.jp/stress_novice_LFHF.html

  • パワースペクトル算出に関して

N=512, サンプリング周波数1kHz

  • 時系列データのcsv

csvは下記のような形式です

csv

1time,value 20,561 31,571 42,555 53,547 64,537 75,527 86,530 97,537 108,535 119,531 12(...) 13510,565 14511,566

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

  • 本来ならば、上のリンク先のように、x軸0.2~0.4あたりに山などができるはずなのですが、0付近でy軸が大きい値を示して以降、ほとんど0のままになっています。
  • 該当コードで意図したようにパワースペクトル密度推定ができていない気がしています。

該当のソースコード

python

1import numpy as np 2from scipy import signal 3import matplotlib.pyplot as plt 4import pandas as pd 5 6# サンプル数 7n = 512 8dt = 0.001 9# サンプリング周波数 10fs = 1/dt 11 12ds = pd.read_csv('data.csv') 13ds = ds['value'] 14data = [] 15for i in range(512): 16 data.append(ds[i]) 17y = np.array(data) 18 19freq1, P1 = signal.periodogram(y, fs) 20freq2, P2 = signal.welch(y, fs) 21 22plt.figure() 23plt.plot(freq1*0.001, P1, "b", label="periodogram") 24plt.plot(freq2*0.001, P2, "c", linewidth=2, label="nseg=n/4") 25plt.legend(loc="upper right") 26plt.xlabel("Frequency[kHz]") 27plt.ylabel("PSD") 28plt.show() 29
  • 出力結果

イメージ説明

試したこと

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

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

参考

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

データが一部しかわからないので、実際にサンプルデータを生成してやってみたところコード自体に問題はないように思います。

周波数0.2KHz、振幅20と周波数0.3KHz、振幅10のサンプルデータを生成

python

1import numpy as np 2import pandas as pd 3 4 5fs = 1000 # サンプリング周波数 6n = 512 # サンプル数 7mean = 550 # 平均 8wave_parameters = ((200, 20), (300, 10)) # 生成する波の周波数と振幅 9 10dt = 1 / fs 11ts = np.linspace(0, dt * n, n, endpoint=False) 12values = np.ones(n) * mean 13for freq, amp in wave_parameters: 14 values += amp * np.sin(ts * 2 * np.pi * freq) 15values = values.astype(np.int) 16 17df = pd.DataFrame(values, columns=('value',)) 18print(df) 19df.to_csv('data.csv', index_label='time')

結果のCSVの中身

csv

1time,value 20,550 31,578 42,555 5... 6509,521 7510,549 8511,578

ご質問のコードを実行した結果
Fig.1

つまり質問者さんが用いたデータの周波数成分は実際にご質問のグラフのようになっているのではないでしょうか?

x軸0.2~0.4あたりに山などができるはず

推測ですが元のデータのサンプリング周波数、あるいはデータそのものの周波数を勘違いしているのではないでしょうか?(例えば実際の波の周波数が0.2Hz~0.4Hzだったとか...)

アドバイス:
実際のデータをいきなり分析する前にまずは周波数や振幅などがわかっているテスト用データで確認してみるとよいと思います。

投稿2019/01/24 03:47

KSwordOfHaste

総合スコア18392

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

plato

2019/01/26 10:11

ご丁寧にありがとうございます。 > つまり質問者さんが用いたデータの周波数成分は実際にご質問のグラフのようになっているのではないでしょうか? おっしゃるとおりでした。 また、アドバイスまでいただきありがとうございます。とても助かりました。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問