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

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

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

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

Python

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

Q&A

0回答

4111閲覧

pythonで周波数の強さをグラフにプロットする

kamikura

総合スコア16

Python 3.x

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

Python

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

0グッド

0クリップ

投稿2019/01/11 03:50

こんにちは。
次元の変換が全く分からなくなってしまったため、ご質問です。

先にコードを載せておきます。少々長くなりますが…

python

1from __future__ import division 2import numpy as np 3import matplotlib.pyplot as plt 4import pycwt as wavelet 5from pycwt.helpers import find 6import pycwt as mothers 7from scipy import signal 8import pandas as pd 9 10#分析データ読み込み 11url = open("kindeni_data_daitai.txt") 12dat = np.genfromtxt(url) 13 14#メタデータ設定 15title = "EMG Time chart" 16label = "EMG[V]" 17units = "[V]" 18 19t0 = 0 20dt = 0.1 21 22msec = "msec" 23gyaku = np.reciprocal(dt) 24timetable = ("msec*{0}".format(gyaku)) 25 26#時間データ生成 27N = dat.size 28t = np.arange(0,N) * dt + t0 29 30#print("{}".format(t)) 31 32#分析対象データの整形 33p = np.polyfit(t - t0, dat, 1) 34dat_notrend = dat - np.polyval(p, t - t0) 35std = dat_notrend.std()#標準偏差 36var = std ** 2#分散 37dat_norm = dat_notrend / std#正規化 38 39#ウェーブレット変換のパラメタ設定 40mother = mothers.Morlet(6) 41s0 = 2 * dt 42dj = 1 / 12 43J = 7 / dj 44alpha, _, _ = wavelet.ar1(dat) 45 46#ウェーブレット変換と逆ウェーブレット変換 47wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J, mother) 48iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std 49 50#ウェーブレットとフーリエの各スペクトル算出 51power = (np.abs(wave)) ** 2 52fft_power = np.abs(fft) ** 2 53period = 1 / freqs#適当な間隔 54 55period = period * 12 56 57#規格化 58power /= scales[:, None] 59 60#パワースペクトル95%信頼区間での有意性検証 61signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha, significance_level=0.95, wavelet=mother) 62sig95 = np.ones([1, N]) * signif[:, None] 63sig95 = power / sig95 64 65#グローバルウェーブレットスペクトルとその有意性の算出→出力潰れる→いる? 66glbl_power = power.mean(axis=1) 67dof = N - scales 68glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha, significance_level=0.95, dof=dof, wavelet=mother) 69 70#スケールの平均値とその有意性の算出 71sel = find((period >=0) & (period < 350))#いつからいつまでの有意性を求めるか 72Cdelta = mother.cdelta 73scale_avg = (scales * np.ones((N, 1))).transpose() 74scale_avg = power / scale_avg 75scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0) 76scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha, significance_level=0.95, dof=[scales[sel[0]],scales[sel[-1]]], wavelet=mother) 77#print("{}".format(scale_avg)) 78summ = np.array(scale_avg).T.tolist() 79#print(type(summ)) 80#print(scale_avg_signif) 81 82data = summ * t 83#print(data) 84np.savetxt("writing.txt",data) 85 86#ウェーブレット解析結果の可視化 87plt.close("all")#現在表示しているプロット結果(グラフ)を全て閉じる 88plt.ioff() 89figprops = dict(figsize=(11, 8), dpi=72) 90fig = plt.figure(**figprops) 91 92#変換などをかけない生の波形 93ax = plt.axes([0.1, 0.75, 0.55, 0.2]) 94ax.plot(t, iwave, "-", linewidth=1, color=[0.5, 0.5, 0.5]) 95ax.plot(t, dat, "k", linewidth=1.5) 96ax.set_title("a) {}".format(title)) 97ax.set_ylabel(r"{}".format(label)) 98 99#パワースペクトル可視化 100bx = plt.axes([0.1, 0.37, 0.55, 0.28], sharex=ax) 101levels = [0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32] 102bx.contourf(t, np.log2(period), np.log2(power), np.log(levels), extend="both", cmap=plt.cm.viridis)#log2倍したperiodが縦軸になってる、横軸は時間で間違いない 103extent = [t.min(), t.max(), 0, max(np.log2(period))] 104bx.contour(t, np.log2(period), sig95, [-99, 1], colors="k", linewidths=2, extent=extent) 105bx.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt, 106 t[:1] - dt, t[:1] - dt]), 107 np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]), 108 np.log2(period[-1:]), [1e-9]]), 109 "k", alpha=0.3, hatch="x") 110bx.set_title("b){} Wavelet Power Spectrum ({})".format(label, mother.name)) 111bx.set_ylabel("Period(syuki)")#period = 周期 112 113Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())), np.ceil(period.max())) 114bx.set_yticks(np.log2(Yticks)) 115bx.set_yticklabels(Yticks) 116 117#グローバルウェーブレットスペクトル可視化 118cx = plt.axes([0.72, 0.37, 0.20, 0.28], sharey=bx) 119cx.plot(glbl_signif, np.log2(period), "k--") 120cx.plot(var * fft_theor, np.log2(period), "--", color="#cccccc") 121cx.plot(var * fft_power, np.log2(1./fftfreqs), "-", color="#cccccc", linewidth=1.) 122cx.plot(var * glbl_power, np.log2(period), "k-", linewidth=1.5) 123cx.set_title("c) Global Wavelet Spectrum") 124cx.set_xlabel(r"Power {}^2".format(units)) 125cx.set_xlim([0, glbl_power.max() + var]) 126cx.set_ylim(np.log2([period.min(), period.max()])) 127cx.set_yticks(np.log2(Yticks)) 128cx.set_yticklabels(Yticks) 129plt.setp(cx.get_yticklabels(), visible=False) 130 131#振幅スペクトル可視化 132dx = plt.axes([0.1, 0.07, 0.40, 0.2], sharex=ax)#グラフの配置指定、 133dx.axhline(scale_avg_signif, color="k", linestyle="--", linewidth=1.)#xに対する平行な線 134dx.plot(t, scale_avg, "k-", linewidth=1.5) 135 136dx.set_title("d) {}--{} EMG[V] power spectrum".format(t0, max(t))) 137dx.set_xlabel("Time [{}]".format(timetable)) 138dx.set_ylabel(r"Scale Avg") 139ax.set_xlim([t.min(), t.max()]) 140 141#ピーク検出 142freq = np.linspace(0, 16.0/dt, N) #周波数軸 143 144maximal_idx = signal.argrelmax(data, order=1)[0]#ピーク検出、オーダー値によって両側いくつを対象にするか 145peak_cut = 1.5#ピーク閾値、検出感度調整用 146maximal_idx = maximal_idx[(data[maximal_idx] > peak_cut) & (maximal_idx <= N/2)] 147 148#グラフ表示位置と縦横のラベル 149ex = plt.axes([0.575, 0.07, 0.40, 0.20])#グラフ位置 150plt.xlabel('Frequency(Hz)') 151plt.ylabel('Amplitude') 152 153#ピーク検出用グラフについて 154plt.axis([t0, max(np.log2(period)), t0, max(scale_avg)])#座標軸の設定 155plt.plot(freq, scale_avg)#縦軸=振幅、ここの横軸は周波数だけど… 156plt.plot(freq[maximal_idx], scale_avg[maximal_idx],'ro')#ピーク点表示 157 158 159print('peak', freq[maximal_idx])#ピーク出現位置の表示 160print("5秒までのピーク個数", sum(1 for i in freq[maximal_idx]), "個")#検出ピーク個数の表示→全体での個数は? 161 162plt.show()

上記のコードの中の最後の方、#ピーク検出用グラフについて、という項目についての質問になります。
ピークを求めたいのですが、その元となるグラフがうまく出力出来ない状態です。
縦軸にデータの振幅が来るのは分かったのですが、横軸にどうしても周波数を落とし込めません。
b)のグラフの周波数情報を用いれればと思ったのですが、変換もうまくいかず詰まっております。

今回使用しているデータをdrop_box様に載せております↓。
https://www.dropbox.com/h

一応、freqsやscale_avgの中身も見ては見ましたが…なぜうまくいかないのか、としか…。

質問内容が非常に不明瞭とは思いますが、どうかご回答いただければ、もしくはご教授いただければと思います。
また質問内容に不備や疑問点がありましたらこちらにご質問ください。
回答出来る限りでお答えします。

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

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

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

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

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

guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問