# -*- coding: utf-8 -*-
"""
Created on Fri Dec 27 19:38:40 2019

@author:Dr. Raphael Vallat (using with permission) edited by Kira Bukharina
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import signal


sf=250.
eeg = pd.read_csv("s013.csv")
#data = np.loadtxt('s01.txt')
#sns.set(font_scale=1.6)


#reading T-electrodes:
eeg_sampleT3=np.squeeze(np.asarray(eeg.iloc[1:75000,[7]]))
eeg_sampleT4=np.squeeze(np.asarray(eeg.iloc[1:75000,[2]]))

eeg_sampleT5=np.squeeze(np.asarray(eeg.iloc[1:75000,[8]]))
eeg_sampleT6=np.squeeze(np.asarray(eeg.iloc[1:75000,[3]]))


#eeg_sampleT4=eeg.iloc[1:211200,[0]]
t=np.linspace(0,eeg_sampleT6.size/250,eeg_sampleT6.size)

fig, ax=plt.subplots()
ax.plot(t,eeg_sampleT3,'-b',label='eeg_sampleT4',linewidth=3)
#ax.plot(t,eeg_sampleT4,'--r',label='eeg_sampleT4')
axes=plt.gca()
axes.set_xlim([0,600])

plt.xlabel('Time(s)')
plt.ylabel('Volts')
plt.title('T3')



#Spectral Analysis
win = 2 * sf
freqs, psd = signal.welch(eeg_sampleT6, sf, nperseg=win)


sns.set(font_scale=1.2, style='white')
plt.figure(figsize=(8, 4))

plt.plot(freqs, psd, color='k', lw=2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (V^2 / Hz)')
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
plt.xlim([0, freqs.max()])
plt.xticks(np.arange(0,100,2))
plt.grid(True)

# Define delta lower and upper limits
low, high = 1, 4

#Define normalization range
low_1, high_1 = 1, 30

# Find intersecting values in frequency vector
idx_delta = np.logical_and(freqs >= low, freqs <= high)
norm_total = np.logical_and(freqs >= low_1, freqs <= high_1)




# Plot the power spectral density and fill the delta area
plt.figure(figsize=(7, 4))
plt.plot(freqs, psd, lw=2, color='k')
plt.fill_between(freqs, psd, where=idx_delta, color='skyblue')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (uV^2 / Hz)')
plt.xlim([0, 100])
plt.ylim([0, psd.max() * 1.1])
plt.title("Welch's periodogram")
sns.despine()

from scipy.integrate import simps

# Frequency resolution
freq_res = freqs[1] - freqs[0]  # = 1 / 4 = 0.25

# Compute the absolute power by approximating the area under the curve
delta_power = simps(psd[idx_delta], dx=freq_res)
print('Absolute delta power: %.15f uV^2' % delta_power)

#total for normalization
total_power_1_30 = simps(psd[norm_total], dx=freq_res)

# Relative delta power (expressed as a percentage of total power)
total_power = simps(psd, dx=freq_res)
delta_rel_power = delta_power / total_power
delta_rel_power_1_30 = delta_power / total_power_1_30
print('Relative delta power: %.3f' % delta_rel_power)
print('Relative delta power, normilized with 1-30Hz range: %.3f' % delta_rel_power_1_30)