<jupyter_start><jupyter_text>EEG explore2021.6.8<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import librosa as lib<jupyter_output><empty_output><jupyter_text>Overview<jupyter_code>df = pd.read_csv(r'data\Epileptic Seizure Recognition.csv')
df.drop('Unnamed', axis=1, inplace=True)
label = df.loc[:, 'y'].values
label[label > 1] = 0
print('"0":nornal. "1":epilepsy')
label
data = df.drop('y', axis=1)
data
data_nor = df[df['y']!=1].drop('y', axis=1)
data_nor.reset_index(drop=True, inplace = True)
data_epi = df[df['y']==1].drop('y', axis=1)
data_epi.reset_index(drop=True, inplace = True)
data_epi
data_nor
timeline = np.linspace(0, 1, 178)
timeline = np.linspace(0, 1, 178)
count_nor = count_epi =0
for i in np.arange(0, 50, 1):
    if df.iloc[i, 178] != 1:
        count_nor += 1
        plt.plot(timeline, data.iloc[i], color='lightseagreen', alpha=0.1+0.9*(i+1)/50, label='normal', lw=0.3)
    else:
        plt.plot(timeline, data.iloc[i], color='lightcoral', alpha=0.1+0.9*(i+1)/50, label='epilepsy', lw=0.3)
        count_epi += 1
print(count_nor, count_epi)
# plt.savefig(r'output\Summary.png', bbox_inches = 'tight', dpi=300)
fig = plt.figure()
ax = plt.axes(projection = '3d')

x, y = np.meshgrid(timeline, np.arange(0, 10, 1)) 
z = data_epi.iloc[0:10, 0:178].values
print(z)
ax.set_xlabel('t') 
ax.set_ylabel('Sample')
ax.set_zlabel('A') 


ax.plot_surface(x, y, z, rstride = 1, cstride = 1, cmap ='cool', edgecolor ='none') 
ax.set_title('Fig.2 EEG',y=-0.25)<jupyter_output>[[  386   382   356 ...   154   143   129]
 [ -278  -246  -215 ...  -125   -79   -40]
 [ -167  -230  -280 ...   215   165   103]
 ...
 [ -340  -381  -376 ...  -380  -350  -308]
 [  -31   -10     5 ...    38    35    15]
 [  512   351   -90 ...   -99  -652 -1089]]<jupyter_text>Cross zero rate<jupyter_code>x = 0
data_epi_N = data_epi.sub(data_epi.mean(axis=1), axis=0)

data_epi_pos = data_epi_N[data_epi_N > 0]
data_epi_neg = data_epi_N[data_epi_N < 0]

data_epi_pospeek = data_epi_N.sub(1.5*data_epi_pos.mean(axis=1), axis=0)
data_epi_negpeek = data_epi_N.sub(1.5*data_epi_neg.mean(axis=1), axis=0)

pospeek = lib.zero_crossings(np.array(data_epi_pospeek.iloc[x]))
print(pospeek)
print(np.sum(pospeek==1))

negpeek = lib.zero_crossings(np.array(data_epi_negpeek.iloc[x]))
print(negpeek)
print(np.sum(negpeek==1))
print(data_epi_negpeek.iloc[x])
plt.plot(timeline, data_epi_negpeek.iloc[x])
plt.plot(timeline, data_epi.iloc[x])
data_nor_N = data_nor.sub(data_nor.mean(axis=1), axis=0)

data_nor_pos = data_nor_N[data_nor_N > 0]
data_nor_neg = data_nor_N[data_nor_N < 0]

data_nor_pospeek = data_nor_N.sub(1.5*data_nor_pos.mean(axis=1), axis=0)
data_nor_negpeek = data_nor_N.sub(1.5*data_nor_neg.mean(axis=1), axis=0)

pospeek = lib.zero_crossings(np.array(data_nor_pospeek.iloc[x]))
print(pospeek)
print(np.sum(pospeek==1))

negpeek = lib.zero_crossings(np.array(data_nor_negpeek.iloc[x]))
print(negpeek)
print(np.sum(negpeek==1))
print(data_nor_pospeek.iloc[x])
plt.plot(timeline, data_nor_pospeek.iloc[x])
plt.plot(timeline, data_nor.iloc[x])
x = 300
plt.plot(timeline, data_nor.iloc[1310])
plt.plot(timeline, data_epi.iloc[1310])
plt.grid()

# 先绝对值，然后0.9max值，然后过零，测量峰宽
zero_array = lib.zero_crossings(np.array(data_epi.iloc[28]))
count_zero = 0
for i in zero_array:
    if i == 1:
        count_zero += 1
print(count_zero)
print(zero_array)
print(np.mean(zero_array))
zero_array = lib.zero_crossings(np.array(data_nor.iloc[x]))
count_zero = 0
for i in zero_array:
    if i == 1:
        count_zero += 1

print(np.mean(zero_array))
print(count_zero)
print(zero_array)   
for i in np.arange(50, 60, 1):
    plt.plot(timeline, data_epi.iloc[i], color='lightseagreen', label='normal', lw=0.3)
from scipy.fftpack import fft, fftfreq
def FFT_contrast(y):
    fs = 179
    # signal
    t = timeline
    N = len(t)
    # fft
    yf = fft(y)
    xf = fftfreq(N, 1/fs)[:N//2]
    yf_nor = 2.0/N * np.abs(yf[0:N//2]) # normalization
    # return np.mean(yf_nor[np.argsort(yf_nor)[::-1][:2]])/np.mean(yf_nor) # contrast 
    return np.sum(yf_nor[0:40]>2*np.mean(yf_nor[0:40]))

data_nor_N = data_nor.sub(data_nor.mean(axis=1), axis=0)

data_nor_pos = data_nor_N[data_nor_N > 0]
data_nor_neg = data_nor_N[data_nor_N < 0]

data_nor_pospeek = data_nor_N.sub(data_nor_pos.mean(axis=1), axis=0)
data_nor_negpeek = data_nor_N.sub(data_nor_neg.mean(axis=1), axis=0)

pospeek = lib.zero_crossings(np.array(data_nor_pospeek.iloc[x]))
print(pospeek)
print(np.sum(pospeek==1))

negpeek = lib.zero_crossings(np.array(data_nor_negpeek.iloc[x]))
print(negpeek)
print(np.sum(negpeek==1))
def nor_SUM_peek(i):
    pospeek = lib.zero_crossings(np.array(data_nor_pospeek.iloc[i]))
    negpeek = lib.zero_crossings(np.array(data_nor_negpeek.iloc[i]))
    return max(np.sum(pospeek==1), np.sum(negpeek==1))
def nor_SUM_peek(i):
    pospeek = lib.zero_crossings(np.array(data_nor_pospeek.iloc[i]))
    negpeek = lib.zero_crossings(np.array(data_nor_negpeek.iloc[i]))
    return max(np.sum(pospeek==1), np.sum(negpeek==1))
nor_SUM_peek(0)<jupyter_output><empty_output><jupyter_text>Feature<jupyter_code>feature = pd.DataFrame(columns=['peek', 'RMS', 'label'])
for i in np.arange(0, len(data_epi), 1):
    peek = SUM_peek(np.array(data_epi.iloc[i]).astype(float))
    rms = lib.feature.rms(y=np.array(data_epi.iloc[i]).astype(float))[0, 0]
    feature = feature.append(pd.DataFrame({'peek':peek, 'RMS':rms, 'label':1}, index=[i]))
for i in np.arange(0, len(data_nor), 1):
    apen = ApEn(np.array(data_nor.iloc[i]).astype(float))
    rms = lib.feature.rms(y=np.array(data_nor.iloc[i]).astype(float))[0, 0]
    feature = feature.append(pd.DataFrame({'ApEn':apen, 'RMS':rms, 'label':0}, index=[i+2300]))
feature.to_csv(r'data\feature.csv')
feature.to_csv(r'data\feature.csv')<jupyter_output><empty_output><jupyter_text>Root mean squre<jupyter_code>lib.feature.rms(y=np.array(data_epi.iloc[0]).astype(float))
lib.feature.rms(y=np.array(data_nor.iloc[500]).astype(float))
f_rms = np.array([])
for i in np.arange(0, len(data_epi), 1):
    rms = lib.feature.rms(y=np.array(data_nor.iloc[i]).astype(float))[0, 0]
    np.append(f_rms, rms))
f_rms<jupyter_output><empty_output><jupyter_text>Apen<jupyter_code>def ApEn(U, m=30):
    r = np.std(U)
    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(U)

    return abs(_phi(m+1) - _phi(m))

# Usage example
U = np.array([85, 80, 89] * 17)
print (ApEn(U))
# 1.0996541105257052e-05
ApEn(np.array(data_nor.iloc[1359]).astype(float))
ApEn(np.array(data_epi.iloc[1359]).astype(float))
# SampEn calculation

import numpy as np

def SampEn(U, m=30):
    r = np.std(U)
    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        B = [(len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) - 1.0) / (N - m) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(B)

    N = len(U)

    return -np.log(_phi(m+1) / _phi(m))

# Usage example
U = np.array([85, 80, 89] *17)
print(SampEn(U))

print(SampEn(np.array(data_nor.iloc[13]).astype(float)))
print(SampEn(np.array(data_epi.iloc[13]).astype(float)))<jupyter_output>0.03069977797529917
0.05399986210963958<jupyter_text>fft<jupyter_code>from scipy.fftpack import fft, fftfreq
def FFT(y):
    fs = 179
    # signal
    t = timeline
    N = len(t)
    # fft
    yf = fft(y)
    xf = fftfreq(N, 1/fs)[:N//2]
    yf_nor = 2.0/N * np.abs(yf[0:N//2]) # normalization
    return xf, yf_nor
from scipy.fftpack import fft, fftfreq
def mean_freq(y):
    fs = 178
    # signal
    t = np.linspace(0, 1, 178)
    N = len(t)
    # fft
    yf = fft(y)
    xf = fftfreq(N, 1/fs)[:N//2]
    yf_nor = 2.0/N * np.abs(yf[0:N//2]) # normalization
    xf_half = xf[:len(yf_nor)]
    return np.mean(xf_half[yf_nor.argsort()[-2:]])
y = np.array(data_nor.iloc[1342]).astype(float)
fs = 178
# signal
t = np.linspace(0, 1, 178)
N = len(t)
# fft
yf = fft(y)
xf = fftfreq(N, 1/fs)[:N//2]
yf_nor = 2.0/N * np.abs(yf[0:N//2]) # normalization
xf_half = xf[:len(yf_nor)]
print(xf)
print(len(yf_nor))
print(len(xf_half))
print(xf_half[yf_nor.argsort()[-2:]])
print(np.mean(xf_half[np.argsort(yf_nor)][-2:]))
mean_freq(np.array(data_nor.iloc[1342]).astype(float))
fft_norx, fft_nory = FFT(np.array(data_nor.iloc[1342]).astype(float))
plt.plot(fft_norx[0:40], fft_nory[0:40])
print(FFT_contrast(np.array(data_nor.iloc[1342]).astype(float)))
print(fft_nory)
fft_epix, fft_epiy = FFT(np.array(data_epi.iloc[1254]).astype(float))
plt.plot(fft_epix[0:40], fft_epiy[0:40])
print(FFT_contrast(np.array(data_epi.iloc[1254]).astype(float)))
print(fft_epiy)<jupyter_output>6
[34.65168539  8.34874383 47.62397666 40.02758655 41.53623786 77.87816412
 85.44306086 49.71271296 47.70148595 28.62239794 37.87876372 48.13446872
 37.00970892 24.110315   10.62052923 39.42931347 13.82459076 17.63452522
  7.993031   16.80284824 12.85368511  7.4575432  15.20446353  3.2365041
  3.29955808  9.21621398 10.63017161 12.99331991  2.78586478  4.4133429
  5.31359074  3.48503118  5.02166932  2.20604125  5.92915275  2.70454021
  3.17061878  2.61435491  1.6813102   4.24074176  2.10983996  3.72117623
  2.92120203  2.82684456  2.66885045  2.61211919  2.39021015  2.07327282
  2.56921399  2.81457424  2.30704833  1.77427012  2.13951193  2.50974593
  2.4518736   2.29173541  2.06199433  2.54783293  2.48184537  2.93430016
  2.26771263  1.70005912  1.55586434  2.35781406  2.22205326  2.75340152
  1.93877559  1.95011065  1.50794653  2.11329198  1.86362451  1.56262266
  2.26709317  1.91756717  2.26467481  2.80446891  2.5796651   1.7215242
  2.10329399  2.40633444  2.25756609  1.61394842  1.[...]<jupyter_text>other<jupyter_code>lib.feature.spectral_centroid(y=np.array(data_epi.iloc[245]).astype(float), sr=1/179)
lib.feature.spectral_centroid(y=np.array(data_nor.iloc[245]).astype(float), sr=1/179)
lib.feature.spectral_flatness(y=np.array(data_epi.iloc[1300]).astype(float))
lib.feature.spectral_flatness(y=np.array(data_nor.iloc[1300]).astype(float))
def stEnergyEntropy(frame, n_short_blocks=10):
    """Computes entropy of energy"""
    Eol = np.sum(frame ** 2)    # total frame energy
    L = len(frame)
    sub_win_len = int(np.floor(L / n_short_blocks))
    if L != sub_win_len * n_short_blocks:
            frame = frame[0:sub_win_len * n_short_blocks]
    # sub_wins is of size [n_short_blocks x L]
    sub_wins = frame.reshape(sub_win_len, n_short_blocks, order='F').copy()

    # Compute normalized sub-frame energies:
    s = np.sum(sub_wins ** 2, axis=0) / (Eol)

    # Compute entropy of the normalized sub-frame energies:
    Entropy = -np.sum(s * np.log2(s))
    return Entropy
stEnergyEntropy(np.array(data_nor.iloc[1300]).astype(float))
stEnergyEntropy(np.array(data_epi.iloc[1300]).astype(float))
x = np.array([3, 4, 1, 6, 2, 8])
print(x.argsort()[-2:])
df1 = pd.read_csv(r'data\feature.csv', index_col=0)
df2 = pd.read_csv(r'data\feature2.csv', index_col=0)
df2 = df2.drop('label', axis=1)
df2
df1
new = pd.concat([df1, df2], axis=1)
new = new[['ApEn', 'peek', 'RMS', 'mean frequency', 'label']]
new.to_csv(r'data\feature.csv')
x = lib.util.frame(np.pad(np.array(data_nor.iloc[80]).astype(float), int(2048 // 2), mode="reflect"), frame_length=2048, hop_length=512)
plt.plot(np.linspace(0, 1, 2048), x.flatten())
print(x.shape)
plt.plot(np.linspace(0, 1, 178), np.array(data_epi.iloc[600]))
plt.grid()
plt.plot(np.linspace(0, 1, 178), np.array(data_epi.iloc[1310]))
plt.grid()
plt.plot(np.linspace(0, 1, 89)[0:30], FFT(np.array(data_epi.iloc[130]))[0:30])
plt.grid()
plt.plot(np.linspace(0, 1, 178), np.array(data_nor.iloc[3]))
plt.grid()
plt.plot(np.linspace(0, 1, 89)[0:30], FFT(np.array(data_nor.iloc[3]))[0:30])
plt.grid()
plt.plot(np.linspace(0, 1, 178), data_nor_pospeek.iloc[3])
plt.grid()
plt.plot(np.linspace(0, 1, 178), data_nor_negpeek.iloc[3])
plt.grid()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17966436327739388, 'fraction_numerical': 0.11154985192497532, 'mean_word_length': 4.543547651618787, 'pattern_counts': {'":': 2, '<': 35, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '460993', 'license': 'MIT', 'n_tokens_mistral': 6455, 'n_tokens_neox': 5646, 'n_words': 872, 'path': 'EXPLORE/EEG_explore.ipynb', 'repo_name': 'H-cyan/EEG-Epilepsy-Detection', 'size': '12156'}