ローレンツプロット - yana_lab

ローレンツプロット

ローレンツプロットは時系列信号のn秒のデータをx軸，n+1軸をy軸にとり，信号をプロットする方法です．プロットの面積を求めることで，信号のバラツキ（カオス性）を求めることができます．例えば，心拍数でローレンツプロットの面積を求めると，ストレス状態のほうが面積が小さくなることが知られています．

栁澤研究室では発達障害の１つであるASDに注目し，ASD者の安静時のNIRS信号のローレンツプロットの面積を求めた結果，定型発達者（健常者）の面積より大きい傾向にあることを確認しています．

プログラム関連

ローレンツプロットの面積を求めるプログラム

import math
import numpy as np
 
def lorenz_plot_area(restdata):
    nn = len(restdata) #データの長さ
    x_data = restdata[0:nn-1] #x軸となるデータ
    y_data = restdata[1:nn] #y軸になるデータ（x軸のデータを1点ずらしたデータ）
    rotation_x_data = [0]; #回転後ｘ軸データ
    rotation_y_data = [0]; #回転後y軸データ

    #すべてのデータを-45度回転させる．###########################################
    sita = -1*(math.pi/4)
    for n in range(0,nn-1): #データ数が1点減るので、-1してる
        xx = x_data[n] * math.cos(sita) -1 * y_data[n] * math.sin(sita)
        rotation_x_data = np.hstack([rotation_x_data, xx])
        yy = x_data[n] * math.sin(sita) + y_data[n] * math.cos(sita)
        rotation_y_data = np.hstack([rotation_y_data, yy])
    ########################################################################
        
    #配列の初期値を0にしていたので、そこを捨てる
    rotation_x_data = rotation_x_data[1:nn+1]
    rotation_y_data = rotation_y_data[1:nn+1]
 
    std_x = np.std(rotation_x_data) #X軸の標準偏差
    std_y = np.std(rotation_y_data) #X軸の標準偏差

    area = (std_x * std_y * math.pi)/4 #楕円として面積を求める

    return area
 

ローレンツプロットのグラフ作成プログラム

import glob
import lorenz_plot
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
data=glob.glob("*.csv")
 
for i in data:
    print(i)
 
    plt.rcParams['font.family'] = 'Times New Roman'#フォントをTimes New Romanにする
    
    data = np.loadtxt(i, delimiter=',')#列を飛ばす
    data = data[:360,3]
    name=i.split("\\")[-1].split(".")[0]+"lorenz.png"#名前付け

    xdata = data[0:359]#nデータ
    ydata = data[1:360]#n+1データ

    area = lorenz_plot.lorenz_plot_area(data)#lorenz_plotのプログラムから面積を求める
    print(area)
 
    fig=plt.figure(figsize=(9,9))#グラフの大きさ
    ax1 = fig.add_subplot(1,1,1)
 
    ax1.scatter(
        xdata,#x軸
        ydata,#y軸
        color="r",
        label = "oxy-Hb",
        linewidth = 5,
        marker = '.')
 
##    plt.title("Rest of lorenz (Left outer portion)")#グラフのタイトル

    ax1.set_xlabel("n")#x軸ラベル
    ax1.set_ylabel("n+1")#y軸ラベル

    plt.show()#グラフ確認コマンド

##    plt.savefig(name)#保存コマンド
##    plt.close("all")

print("完了")