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

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

ただいまの
回答率

90.43%

  • Python

    9763questions

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

  • Python 3.x

    7901questions

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

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

受付中

回答 0

投稿

  • 評価
  • クリップ 0
  • VIEW 82

kamikura

score 6

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

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

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import pycwt as wavelet
from pycwt.helpers import find
import pycwt as mothers
from scipy import signal
import pandas as pd

#分析データ読み込み
url = open("kindeni_data_daitai.txt")
dat = np.genfromtxt(url)

#メタデータ設定
title = "EMG Time chart"
label = "EMG[V]"
units = "[V]"

t0 = 0
dt = 0.1

msec = "msec"
gyaku = np.reciprocal(dt)
timetable = ("msec*{0}".format(gyaku))

#時間データ生成
N = dat.size
t = np.arange(0,N) * dt + t0

#print("{}".format(t))

#分析対象データの整形
p = np.polyfit(t - t0, dat, 1)
dat_notrend = dat - np.polyval(p, t - t0)
std = dat_notrend.std()#標準偏差
var = std ** 2#分散
dat_norm = dat_notrend / std#正規化

#ウェーブレット変換のパラメタ設定
mother = mothers.Morlet(6)
s0 = 2 * dt
dj = 1 / 12
J = 7 / dj
alpha, _, _ = wavelet.ar1(dat)

#ウェーブレット変換と逆ウェーブレット変換
wave, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(dat_norm, dt, dj, s0, J, mother)
iwave = wavelet.icwt(wave, scales, dt, dj, mother) * std

#ウェーブレットとフーリエの各スペクトル算出
power = (np.abs(wave)) ** 2
fft_power = np.abs(fft) ** 2
period = 1 / freqs#適当な間隔

period = period * 12

#規格化
power /= scales[:, None]

#パワースペクトル95%信頼区間での有意性検証
signif, fft_theor = wavelet.significance(1.0, dt, scales, 0, alpha, significance_level=0.95, wavelet=mother)
sig95 = np.ones([1, N]) * signif[:, None]
sig95 = power / sig95

#グローバルウェーブレットスペクトルとその有意性の算出→出力潰れる→いる?
glbl_power = power.mean(axis=1)
dof = N - scales
glbl_signif, tmp = wavelet.significance(var, dt, scales, 1, alpha, significance_level=0.95, dof=dof, wavelet=mother)

#スケールの平均値とその有意性の算出
sel = find((period >=0) & (period < 350))#いつからいつまでの有意性を求めるか
Cdelta = mother.cdelta
scale_avg = (scales * np.ones((N, 1))).transpose()
scale_avg = power / scale_avg
scale_avg = var * dj * dt / Cdelta * scale_avg[sel, :].sum(axis=0)
scale_avg_signif, tmp = wavelet.significance(var, dt, scales, 2, alpha, significance_level=0.95, dof=[scales[sel[0]],scales[sel[-1]]], wavelet=mother)
#print("{}".format(scale_avg))
summ = np.array(scale_avg).T.tolist()
#print(type(summ))
#print(scale_avg_signif)

data = summ * t
#print(data)
np.savetxt("writing.txt",data)

#ウェーブレット解析結果の可視化
plt.close("all")#現在表示しているプロット結果(グラフ)を全て閉じる
plt.ioff()
figprops = dict(figsize=(11, 8), dpi=72)
fig = plt.figure(**figprops)

#変換などをかけない生の波形
ax = plt.axes([0.1, 0.75, 0.55, 0.2])
ax.plot(t, iwave, "-", linewidth=1, color=[0.5, 0.5, 0.5])
ax.plot(t, dat, "k", linewidth=1.5)
ax.set_title("a) {}".format(title))
ax.set_ylabel(r"{}".format(label))

#パワースペクトル可視化
bx = plt.axes([0.1, 0.37, 0.55, 0.28], sharex=ax)
levels = [0.03125, 0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32]
bx.contourf(t, np.log2(period), np.log2(power), np.log(levels), extend="both", cmap=plt.cm.viridis)#log2倍したperiodが縦軸になってる、横軸は時間で間違いない
extent = [t.min(), t.max(), 0, max(np.log2(period))]
bx.contour(t, np.log2(period), sig95, [-99, 1], colors="k", linewidths=2, extent=extent)
bx.fill(np.concatenate([t, t[-1:] + dt, t[-1:] + dt,
                        t[:1] - dt, t[:1] - dt]), 
        np.concatenate([np.log2(coi), [1e-9], np.log2(period[-1:]),
                       np.log2(period[-1:]), [1e-9]]),
        "k", alpha=0.3, hatch="x")
bx.set_title("b){} Wavelet Power Spectrum ({})".format(label, mother.name))
bx.set_ylabel("Period(syuki)")#period = 周期

Yticks = 2 ** np.arange(np.ceil(np.log2(period.min())), np.ceil(period.max()))
bx.set_yticks(np.log2(Yticks))
bx.set_yticklabels(Yticks)

#グローバルウェーブレットスペクトル可視化
cx = plt.axes([0.72, 0.37, 0.20, 0.28], sharey=bx)
cx.plot(glbl_signif, np.log2(period), "k--")
cx.plot(var * fft_theor, np.log2(period), "--", color="#cccccc")
cx.plot(var * fft_power, np.log2(1./fftfreqs), "-", color="#cccccc", linewidth=1.)
cx.plot(var * glbl_power, np.log2(period), "k-", linewidth=1.5)
cx.set_title("c) Global Wavelet Spectrum")
cx.set_xlabel(r"Power {}^2".format(units))
cx.set_xlim([0, glbl_power.max() + var])
cx.set_ylim(np.log2([period.min(), period.max()]))
cx.set_yticks(np.log2(Yticks))
cx.set_yticklabels(Yticks)
plt.setp(cx.get_yticklabels(), visible=False)

#振幅スペクトル可視化
dx = plt.axes([0.1, 0.07, 0.40, 0.2], sharex=ax)#グラフの配置指定、
dx.axhline(scale_avg_signif, color="k", linestyle="--", linewidth=1.)#xに対する平行な線
dx.plot(t, scale_avg, "k-", linewidth=1.5)

dx.set_title("d) {}--{} EMG[V] power spectrum".format(t0, max(t)))
dx.set_xlabel("Time [{}]".format(timetable))
dx.set_ylabel(r"Scale Avg")
ax.set_xlim([t.min(), t.max()])

#ピーク検出
freq = np.linspace(0, 16.0/dt, N) #周波数軸

maximal_idx = signal.argrelmax(data, order=1)[0]#ピーク検出、オーダー値によって両側いくつを対象にするか
peak_cut = 1.5#ピーク閾値、検出感度調整用
maximal_idx = maximal_idx[(data[maximal_idx] > peak_cut) & (maximal_idx <= N/2)]

#グラフ表示位置と縦横のラベル
ex = plt.axes([0.575, 0.07, 0.40, 0.20])#グラフ位置
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

#ピーク検出用グラフについて
plt.axis([t0, max(np.log2(period)), t0, max(scale_avg)])#座標軸の設定
plt.plot(freq, scale_avg)#縦軸=振幅、ここの横軸は周波数だけど…
plt.plot(freq[maximal_idx], scale_avg[maximal_idx],'ro')#ピーク点表示


print('peak', freq[maximal_idx])#ピーク出現位置の表示
print("5秒までのピーク個数", sum(1 for i in freq[maximal_idx]), "個")#検出ピーク個数の表示→全体での個数は?

plt.show()

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

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

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

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

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

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

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

    クリップを取り消します

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

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

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

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

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

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

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

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

    質問の評価を下げる

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

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

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

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

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

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

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

    詳細な説明はこちら

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

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

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

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

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

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

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

  • Python

    9763questions

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

  • Python 3.x

    7901questions

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