瞳孔特征值提取，blink frequency,fixation frequency,saccad extent, pupil diameter等

进行的分析有，滤波分析，fft，psd

database.py

下面展示一些 内联代码片。

import pandas as pd
import numpy as np

def read_file(raw_path):
    data = pd.DataFrame(pd.read_csv(raw_path, encoding="utf-8", engine='python'))
    LPMM = np.array(data['LPMM'])
    RPMM = np.array(data['RPMM'])
    BKPMIN = np.array(data['BKPMIN'])
    FPOGD = np.array(data['FPOGD'])
    SACCADE_MAG = np.array(data['SACCADE_MAG'])

    NEW_BKPMIN = []
    for i in range(60,len(BKPMIN)):
        new_BKPMIN = BKPMIN[i]/60.0
        NEW_BKPMIN.append(new_BKPMIN)


    return SACCADE_MAG

main.py

import dataset
import iir
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt, hilbert,find_peaks
import numpy as np
import pandas as pd
from scipy import signal
from scipy import fftpack
from scipy.fftpack import fft, fftshift, ifft
import os
def mkdir(path_3):
    folder = os.path.exists(path_3)
    if not folder:  # 判断是否存在文件夹如果不存在则创建为文件夹
        os.makedirs(path_3)  # makedirs 创建文件时如果路径不存在会创建这个路径
        print("---  new folder...  ---")
        print("---  OK  ---")
    else:
        print("---  There is this folder!  ---")
def main():
    #读取 pupil data
    data = dataset.read_file(raw_path +str(num[v]) + '_all_gaze.csv')
    m= 0; num_fft =64; Mean_psd = [];Std_psd = []

    for i in range(int(len(data)/60)):

        # Bandpass filter
        filter_delta = iir.butter_bandpass_filter_delta(data[i+m:i+m+599], lowcut_delta, highcut_delta, fs, order=5)
        filter_theta = iir.butter_bandpass_filter_theta(data[i+m:i+m+599], highcut_delta, highcut_theta, fs, order=5)
        filter_alpha = iir.butter_bandpass_filter_alpha(data[i+m:i+m+599], highcut_theta, highcut_alpha, fs, order=5)
        filter_beta = iir.butter_bandpass_filter_beta(data[i+m:i+m+599], highcut_alpha, highcut_beta, fs, order=5)
        filter_gamma = iir.butter_bandpass_filter_gamma(data[i+m:i+m+599], highcut_beta, highcut_gamma, fs, order=5)
        filter_gamma_high = iir.butter_bandpass_filter_gamma_high(data[i + m:i + m + 599], highcut_beta_high, highcut_gamma_high, fs,
                                                        order=5)
        filter_total = iir.butter_bandpass_filter_total(data[i+m:i+m+599], lowcut_delta, highcut_gamma, fs, order=5)

        #sliding
        analytical_delta = hilbert(filter_delta)
        analytical_theta = hilbert(filter_theta)
        analytical_alpha = hilbert(filter_alpha)
        analytical_beta = hilbert(filter_beta)
        analytical_gamma = hilbert(filter_gamma)
        analytical_gamma_high = hilbert(filter_gamma_high)
        analytical_total = hilbert(filter_total)
        delta_envelope = np.abs(analytical_delta)
        theta_envelope = np.abs(analytical_theta)
        alpha_envelope = np.abs(analytical_alpha)
        beta_envelope = np.abs(analytical_beta)
        gamma_envelope = np.abs(analytical_gamma)
        gamma_high_envelope = np.abs(analytical_gamma_high)
        total_envelope = np.abs(analytical_total)
        """
                FFT(Fast Fourier Transformation)快速傅里叶变换
       """
        Y_data = fft(data[i+m:i+m+599], num_fft)
        Y_data = np.abs(Y_data)
        psd_data = Y_data** 2
        fft_data = 20 * np.log10(Y_data[:num_fft // 2])
        Psd_data = 20 * np.log10(psd_data[:num_fft // 2])
        Y_delta = fft(filter_delta, num_fft)
        Y_delta = np.abs(Y_delta)
        psd_delta = Y_delta ** 2
        fft_delta = 20 * np.log10(Y_delta[:num_fft // 2])
        Psd_delta = 20 * np.log10(psd_delta[:num_fft // 2])
        Y_theta = fft(filter_theta, num_fft)
        Y_theta = np.abs(Y_theta)
        psd_theta = Y_theta ** 2
        fft_theta = 20 * np.log10(Y_theta[:num_fft // 2])
        Psd_theta = 20 * np.log10(psd_theta[:num_fft // 2])
        Y_alpha = fft(filter_alpha, num_fft)
        Y_alpha = np.abs(Y_alpha)
        psd_alpha = Y_alpha ** 2
        fft_alpha = 20 * np.log10(Y_alpha[:num_fft // 2])
        Psd_alpha = 20 * np.log10(psd_alpha[:num_fft // 2])
        Y_beta = fft(filter_beta, num_fft)
        Y_beta = np.abs(Y_beta)
        psd_beta = Y_beta ** 2
        fft_beta = 20 * np.log10(Y_beta[:num_fft // 2])
        Psd_beta = 20 * np.log10(psd_beta[:num_fft // 2])
        Y_gamma = fft(filter_gamma, num_fft)
        Y_gamma = np.abs(Y_gamma)
        psd_gamma = Y_gamma ** 2
        fft_gamma = 20 * np.log10(Y_gamma[:num_fft // 2])
        Psd_gamma = 20 * np.log10(psd_gamma[:num_fft // 2])
        Y_gamma_high = fft(filter_gamma_high, num_fft)
        Y_gamma_high = np.abs(Y_gamma_high)
        psd_gamma_high = Y_gamma_high ** 2
        fft_gamma_high = 20 * np.log10(Y_gamma_high[:num_fft // 2])
        Psd_gamma_high = 20 * np.log10(psd_gamma_high[:num_fft // 2])
        Y_total = fft(filter_total, num_fft)
        Y_total = np.abs(Y_total)
        psd_total = Y_total ** 2
        fft_total = 20 * np.log10(Y_total[:num_fft // 2])
        Psd_total = 20 * np.log10(psd_total[:num_fft // 2])
        mean_psd = '%.3f' % np.mean(Psd_total)
        std_psd = '%.3f' %np.std(Psd_total)

        Mean_psd.append(mean_psd)
        Std_psd.append(std_psd)
        m = m + 59

    return Mean_psd,Std_psd

num = [13,15,16,17,18,20,22,23]
# num = [1,2,5,6]
M_BLINK = [];S_BLINK = []
path_out = 'E:/ffmpeg-latest-win64-static/2019_11_12_experimental_data/peak_detetion(test)/erd_ers_empathy/blink/'
for v in range(len(num)):
    print(num[v])
    Mean_Psd =[]; Std_Psd =[]

    name = ['오령기'
        , '정호중', '원성재', '최기연', '김예린', '정차윤', '강대원', '박동준', '김관엽', '김은해',
            '신현승', '조민성', '이정인', '박성수', '전인배', '김경빈', '김소의', '성현아', '이연주', '최세영',
            '권세웅', '김초록', '임시현', '하다원', '김예진', '강한결', '이정찬', '주용현', '허수현', '김수연']

    for n in range(len(name)):
        print(name[n])
        raw_path = 'E:/ffmpeg-latest-win64-static/2019_11_12_experimental_data/GAZE_DATA_backup/191202/'+name[n]+'/data_process/Experimenter/all_gaze/'
        path = 'E:/ffmpeg-latest-win64-static/2019_11_12_experimental_data/GAZE_DATA_backup/7_29/pupil/1_pupil/'+str(num[v])+'/'

        lowcut_delta = 0.12
        highcut_delta = 0.48
        highcut_theta = 0.96
        highcut_alpha = 1.56
        highcut_beta = 3.6
        highcut_beta_high = 5
        highcut_gamma = 6
        highcut_gamma_high = 10
        fs = 60

        Mean_psd, Std_psd = main()
        Mean_Psd.append(Mean_psd)
        Std_Psd.append(Std_psd)

    M_Total_ = np.hstack(
        (Mean_Psd[1], Mean_Psd[2], Mean_Psd[3], Mean_Psd[4], Mean_Psd[5],Mean_Psd[6], Mean_Psd[7], Mean_Psd[8], Mean_Psd[9], Mean_Psd[10],
         Mean_Psd[11], Mean_Psd[12], Mean_Psd[13], Mean_Psd[14], Mean_Psd[15],Mean_Psd[16], Mean_Psd[17], Mean_Psd[18], Mean_Psd[19], Mean_Psd[20],
         Mean_Psd[21], Mean_Psd[22], Mean_Psd[23], Mean_Psd[24], Mean_Psd[25],Mean_Psd[26], Mean_Psd[27], Mean_Psd[28], Mean_Psd[29],Mean_Psd[0]))

    S_Total_ = np.hstack(
        (Std_Psd[1], Std_Psd[2], Std_Psd[3], Std_Psd[4], Std_Psd[5], Std_Psd[6], Std_Psd[7], Std_Psd[8],Std_Psd[9], Std_Psd[10],
         Std_Psd[11], Std_Psd[12], Std_Psd[13], Std_Psd[14], Std_Psd[15], Std_Psd[16], Std_Psd[17], Std_Psd[18],
         Std_Psd[19], Std_Psd[20],
         Std_Psd[21], Std_Psd[22], Std_Psd[23], Std_Psd[24], Std_Psd[25], Std_Psd[26], Std_Psd[27], Std_Psd[28],
         Std_Psd[29], Std_Psd[0]))
    M_BLINK.append(M_Total_)
    S_BLINK.append(S_Total_)
    print(len(M_Total_),len(S_Total_))

M_BLINK_ = np.hstack(
    (M_BLINK[0], M_BLINK[1], M_BLINK[2], M_BLINK[3]
     ,M_BLINK[4],
     M_BLINK[5],M_BLINK[6], M_BLINK[7]))
     # M_BLINK[9], M_BLINK[10],

     # M_BLINK[11], M_BLINK[12], M_BLINK[13], M_BLINK[14], M_BLINK[15],M_BLINK[16], M_BLINK[17], M_BLINK[18], M_BLINK[19],M_BLINK[20],
     # M_BLINK[21], M_BLINK[22], M_BLINK[23], M_BLINK[24], M_BLINK[25],M_BLINK[26], M_BLINK[27], M_BLINK[28], M_BLINK[29],M_BLINK[0]))

S_BLINK_ = np.hstack(
    (S_BLINK[0], S_BLINK[1], S_BLINK[2], S_BLINK[3]
     ,M_BLINK[4],
     S_BLINK[5], S_BLINK[6], S_BLINK[7]))
     # S_BLINK[11], S_BLINK[12], S_BLINK[13], S_BLINK[14], S_BLINK[15], S_BLINK[16], S_BLINK[17], S_BLINK[18],
     # S_BLINK[19], S_BLINK[20],
     # S_BLINK[21], S_BLINK[22], S_BLINK[23], S_BLINK[24], S_BLINK[25], S_BLINK[26], S_BLINK[27], S_BLINK[28],
     # S_BLINK[29], S_BLINK[0]))
print(len(M_BLINK_),len(S_BLINK_))
mkdir(path_out)
BLINK_data = {'M_SACCADE': M_BLINK_,'S_SACCADE': S_BLINK_}
BLINK_result = pd.DataFrame(BLINK_data, columns=['M_SACCADE','S_SACCADE'])
BLINK_result.to_csv(path_out + 'SACCADE_NO_empathy.csv')

irr.py

from scipy.signal import butter, lfilter
from scipy import signal

#delta
def butter_bandpass_delta(lowcut_delta, highcut_delta, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut_delta / nyq
    high = highcut_delta / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_delta(data,lowcut_delta, highcut_delta, fs, order=5):
    b, a = butter_bandpass_delta(lowcut_delta, highcut_delta, fs, order=order)
    delta = signal.filtfilt(b, a, data)
    return delta
#theta
def butter_bandpass_theta(highcut_delta, highcut_theta, fs, order=5):
    nyq = 0.5 * fs
    low = highcut_delta / nyq
    high = highcut_theta / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_theta(data,highcut_delta, highcut_theta, fs, order=5):
    b, a = butter_bandpass_theta(highcut_delta, highcut_theta, fs, order=order)
    theta = signal.filtfilt(b, a, data)
    return theta
#alpha
def butter_bandpass_alpha(highcut_theta, highcut_alpha, fs, order=5):
    nyq = 0.5 * fs
    low = highcut_theta / nyq
    high = highcut_alpha / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_alpha(data,highcut_theta, highcut_alpha, fs, order=5):
    b, a = butter_bandpass_alpha(highcut_theta, highcut_alpha, fs, order=order)
    alpha = signal.filtfilt(b, a, data)
    return alpha
#beta
def butter_bandpass_beta(highcut_alpha, highcut_beta, fs, order=5):
    nyq = 0.5 * fs
    low = highcut_alpha / nyq
    high = highcut_beta / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_beta(data, highcut_alpha, highcut_beta, fs, order=5):
    b, a = butter_bandpass_beta(highcut_alpha, highcut_beta, fs, order=order)
    beta = signal.filtfilt(b, a, data)
    return beta
#gamma
def butter_bandpass_gamma(highcut_beta, highcut_gamma, fs, order=5):
    nyq = 0.5 * fs
    low = highcut_beta / nyq
    high = highcut_gamma / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_gamma(data, highcut_beta, highcut_gamma, fs, order=5):
    b, a = butter_bandpass_gamma(highcut_beta, highcut_gamma, fs, order=order)
    gamma = signal.filtfilt(b, a, data)
    return gamma
#high gamma
def butter_bandpass_gamma_high(highcut_beta_high, highcut_gamma_high, fs, order=5):
    nyq = 0.5 * fs
    low = highcut_beta_high / nyq
    high = highcut_gamma_high / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_gamma_high(data, highcut_beta_high, highcut_gamma_high, fs, order=5):
    b, a = butter_bandpass_gamma_high(highcut_beta_high, highcut_gamma_high, fs, order=order)
    gamma_high = signal.filtfilt(b, a, data)
    return gamma_high
#total
def butter_bandpass_total(lowcut_delta, highcut_gamma, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut_delta / nyq
    high = highcut_gamma / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter_total(data, lowcut_delta, highcut_gamma, fs, order=5):
    b, a = butter_bandpass_total(lowcut_delta, highcut_gamma, fs, order=order)
    total = signal.filtfilt(b, a, data)
    return total