# import文
import numpy as np
import matplotlib.pyplot as plt
import time
import csv,os
from scipy import signal
from mpl_toolkits.mplot3d import Axes3D


# 定数宣言
waveMax = 500
waveMin = -500

# 変数宣言
fileName = 'LE029_0514'
position = 0
positionList = []
os.chdir('F:/実験2019/実験A_2019/A_20190514/')


# トリガーのチャンネルを宣言
std = 8
tgt = 9

# 番号入力
answer = input("分析したいチャンネルの番号を入力ください：\n")

print('\n'+"チャンネルデータ読み込み中……")

# 分析するチャンネルを読み込んでからフィルターをかけて行列変換をする
#          トレンド除去                  ファイル名       行頭抜き   分析する列の番号   行列変換をする
chList = signal.detrend(np.loadtxt(fileName+'.m00', skiprows=2, usecols=int(answer)-1).T)        # トレンド除去フィルターで波形を整形

# バンドパス Chebyshev Ⅰ型　リップル1db（第一種チェビシェフフィルタ）
# フィルターパラメータ設定
dt = 0.001     # サンプリング間隔
fn = 1/(2*dt)  # ナイキスト周波数
fp = 2         # 通過域端周波数[Hz]
fs = 3         # 阻止域端周波数[Hz]
gPass = 1      # 通過域最大損失量[dB]
gStop = 40     # 阻止域最小減衰量[dB]
# 正規化
Wp = fp/fn
Ws = fs/fn
N, Wn = signal.cheb1ord(Wp, Ws, gPass, gStop)
b, a = signal.cheby1(N, gPass, Wn, "low")
#chList = signal.filtfilt(b, a, chList)

# トリガー抽出関数
def triggerExtract(trigger):
    # 関数内変数を宣言
    counter = 0
    continuous = False
    position = 0
    tempPosition = 0
    triggerList = []
    tempList = np.loadtxt(fileName+'.m00',skiprows=2, usecols=trigger-1)
    maxValue = max(tempList)
    # 抽出開始
    for i in tempList:

        if i == 0:
            continuous = False
            pass
            
        if i == maxValue and continuous == False:
            if position - tempPosition > 1200:
                counter = counter + 1
                triggerList.append(position)
                continuous = True
                tempPosition = position
                pass
            pass
        # 累進
        position = position + 1
        pass

    print(counter,"個検出された")
    return triggerList

# 外れ値検出関数
def triggerCheck(listName):
    # 局所変数設定
    position = 0
    badList = []
    tempList = []
    tempList1 = []
    trialList = []
    sessionList = []
    # 検出開始

    for i in listName:                     # i はトリガーの時点

        for j in range(i-200,i+1001):      # j はトリガーの時点前200ms~後1000msの時点
            k = chList[j]                  # k はj時点の脳波の値
            if k>waveMax or k<waveMin:
                tempList.append(position)  # k値は過大か過小ならこのトリガーを除外
                pass
            tempList1.append(k)            # tempList1は波形のもの
            pass

        # ベースライン（-200ms間を平均）
        baseLine = np.mean(tempList1[0:200])
        # トライアルリストのデータ（1トリガー前後1200ms間のデータ）とベースラインの差をセッションリストに導入
        for l in tempList1:
            trialList.append(l-baseLine)
            pass
        sessionList.append(trialList)
        # トライアルリストを初期化
        tempList1 = []
        trialList = []
        # ポジションのカウンターを累進
        position = position + 1
        pass

    # 重複要素を削除
    tempList = list(set(tempList))
    # ソーティング
    tempList.sort()
    # 辞書化
    m = {i:element for i, element in enumerate(sessionList)}
    # 外れ値除外
    [m.pop(i) for i in tempList]
    # リストに戻す
    sessionList = list(m.values())
    # 「012345」から「123456」に変換
    for i in tempList:
        i = i + 1
        badList.append(i)
        pass
    print('\n'+"外れるトリガーは",'\n',badList)
    print("計",len(badList),"個")
    return sessionList


# メインプログラム
# Standard
# トリガー抽出
print('\n'+"Standard トリガー抽出開始……")
stdList = triggerExtract(std)
# 外れ値検出 & データ行列化
stdArray = np.array(triggerCheck(stdList))

# Target
# トリガー抽出
print('\n'+"Target トリガー抽出開始……")
tgtList = triggerExtract(tgt)
# 外れ値検出 & データ行列化
tgtArray = np.array(triggerCheck(tgtList))


# 保存
#np.savetxt(str(answer)+'std.csv',stdArray,delimiter=',')
#np.savetxt(str(answer)+'tgt.csv',tgtArray,delimiter=',')

# 加算平均
#stdMean = np.mean(stdArray,axis=0)
#tgtMean = np.mean(tgtArray,axis=0)

# グラフ作成
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.set_xlabel("Time")
ax.set_ylabel("Trials")
ax.set_zlabel("μV")
ax.set_title(fileName+'_Ch'+answer+"_Standard")
ax.view_init(elev=85, azim=-80)
x = np.arange(-200,1001)
channel = 1
plt.tight_layout()
plt.gca(projection='3d').invert_zaxis()
#fig.gca(projection='3d').get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([0.2, 1, 0.2, 1]))

for i in stdArray:
    y = []
    for j in range(1201):
        y.append(channel)
        pass
    ax.plot(x,y,i,alpha=0.5,color='k')
    channel = channel + 1
    pass
#plt.savefig('standardCh'+answer+'.png')
plt.show()

#plt.xlim(-200,1000)
#plt.xticks(np.arange(-200,1001,100))

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.set_xlabel("Time")
ax.set_ylabel("Trials")
ax.set_zlabel("μV")
ax.set_title(fileName+'_Ch'+answer+"_Target")
plt.gca(projection='3d').invert_zaxis()
#fig.gca(projection='3d').get_proj = lambda: np.dot(Axes3D.get_proj(ax), np.diag([0.2, 1, 0.2, 1]))
plt.tight_layout()
ax.view_init(elev=85, azim=-80)
x = np.arange(-200,1001)
channel = 1
for i in tgtArray:
    y = []
    for j in range(1201):
        y.append(channel)
        pass
    ax.plot(x,y,i,alpha=0.5,color='k')
    channel = channel + 1
    pass
#plt.savefig('standardCh'+answer+'.png')
plt.show()