import sys import os import re import pandas as pd import numpy as np from datetime import datetime, timedelta from tsfresh import extract_features from tsfresh.feature_extraction import MinimalFCParameters from sklearn.preprocessing import MinMaxScaler from sklearn.linear_model import LogisticRegression from sklearn.feature_selection import SelectFromModel import scipy.signal from scipy.stats import entropy #--------------------------------------- # Preprocesamiento de datos #--------------------------------------- #Pasamos como argumento el numero de milisegundos de pre-decisión milisegundos = sys.argv[1] #Se elimina RAW, ya que el resto de canales derivan de RAW canales = ["Delta", "Theta", "Alpha", "Beta", "Gamma"] #En muchos casos, los sensores TP9 y TP10 fallan demasiado, por los que los descartaremos no_valid_columns = ['Delta_TP9', 'Delta_TP10', 'Theta_TP9', 'Theta_TP10', 'Alpha_TP9', 'Alpha_TP10', 'Beta_TP9', 'Beta_TP10', 'Gamma_TP9', 'Gamma_TP10', 'RAW_TP9', 'RAW_AF7', 'RAW_AF8', 'RAW_TP10', 'AUX_RIGHT', 'Accelerometer_X', 'Accelerometer_Y', 'Accelerometer_Z', 'Gyro_X', 'Gyro_Y', 'Gyro_Z', 'HeadBandOn', 'HSI_TP9', 'HSI_AF7', 'HSI_AF8', 'HSI_TP10', 'Battery', 'Mellow', 'Concentration'] X_data = [] #este array contiene todos los datos de cada intento/sujeto (eso son las filas); y en las columnas tendrá las distintas variables que deseemos (por ejemplo, datos en bruto de cada canal; o la media/min/max/etc. de cada canal; o fórmulas más complejas como el cálculo de la PSD/entropía,FFT,etc. por cada canal) all_labels = [] #necesitamos un array con los labels 0/1 (p/q) que identifica a cada fila según si se pulsó p/q para ese dato subjects = [] #necesitamos un array con IDs de los sujetos para cada fila (este array engloba a todos los datos que tengamos, para identificar cada dato a un sujeto, solo incluye enteros para decir si una fila pertenece al sujeto 0, 1, 2, 3...N) nan_columns0 = [] nan_columns1 = [] directorio_local = 'data/' # Directorio donde se encuentran los archivos csv locales directorio_muse = 'GoogleDrive/museData/' directorio_data_procesed = 'data-preprocesed/' archivos = os.listdir(directorio_local) # lista de todos los nombres de los archivos en el directorio regex = re.compile(r'results(\d+)\.csv') #extraer el user_id de los nombres de los archivos n_intentos = 0 def rescale_data(dataframe): scaler = MinMaxScaler() columns = dataframe.columns # Obtener las columnas del DataFrame original # Escalar los datos del DataFrame dataframe_scaled = scaler.fit_transform(dataframe) # Crear un nuevo DataFrame con los datos escalados y las mismas columnas dataframe_scaled = pd.DataFrame(dataframe_scaled, columns=columns) return dataframe_scaled def _phi(m, r, time_series): """ Calculate the ratio of the number of matches for m-length sequences. Args: m: Integer, length of sequences to be compared. r: Float, tolerance for accepting matches. time_series: List of integers or floats representing the time series. Returns: Float representing the ratio. """ x = np.array([time_series[i: i + m] for i in range(len(time_series) - m + 1)]) B = np.sum([np.sum(np.abs(x[:, None] - x[i, :]) <= r, axis=0) - 1 for i in range(x.shape[0])]) return B / ((len(time_series) - m + 1) * (len(time_series) - m)) def sample_entropy(time_series, m=2, r=None): """ Calculate the Sample Entropy of a time series. Args: time_series: List of integers or floats representing the time series. m: Integer, length of sequences to be compared (default: 2). r: Float, tolerance for accepting matches (default: 0.2 * standard deviation of the time series). Returns: Float representing the Sample Entropy. """ if r is None: r = 0.2 * np.std(time_series) return -np.log(_phi(m+1, r, time_series) / _phi(m, r, time_series)) def calculate_differential_entropy(time_series): """ Calculate the Differential Entropy of a Gaussian-distributed random variable. Args: time_series: List of integers or floats representing the time series. Returns: Float representing the Differential Entropy. """ sigma = np.std(time_series) if sigma == 0: sigma = 1**-10 return 0.5 * np.log(2 * np.pi * np.e * sigma**2) def calculate_psd(time_series, nperseg, fs=1.0, noverlap=0): """ Calculate the Power Spectral Density (PSD) of a time series. Args: time_series: List of integers or floats representing the time series. fs: Float, the sample rate of the time series (default: 1.0). nperseg: Integer, the length of each segment for the STFT (default: 256). noverlap: Integer, the number of points to overlap between segments (default: 0). Returns: freqs: Array of floats, the frequencies at which the PSD was calculated. psd: Array of floats, the calculated PSD. """ _, _, Zxx = scipy.signal.stft(time_series, fs=fs, nperseg=nperseg, noverlap=noverlap, window='hann') psd = np.abs(Zxx)**2 mean_psd = np.mean(psd) freqs = np.fft.fftfreq(nperseg, 1/fs) return freqs, mean_psd def hjorth_params(time_series): """ Calculate Hjorth parameters (Activity, Mobility, and Complexity) of a time series. Args: time_series: List of integers or floats representing the time series. Returns: activity: Float, the signal power. mobility: Float, the standard deviation of the power spectrum. complexity: Float, indicating the similarity of the signal and a sinusoidal signal. """ # First derivative of the time series first_derivative = np.diff(time_series) # Second derivative of the time series second_derivative = np.diff(time_series, 2) # Activity activity = np.var(time_series) # Mobility mobility = np.sqrt(np.var(first_derivative) / activity) # Complexity complexity = np.sqrt(np.var(second_derivative) / np.var(first_derivative)) / mobility return activity, mobility, complexity for archivo in archivos: if os.path.isfile(os.path.join(directorio_local, archivo)): match = regex.match(archivo) if match and int(match.group(1)) == 0 : # Si el nombre del archivo coincide con nuestro patrón y el numero es positivo (results-1.csv es mi archivo de pruebas) user_id = match.group(1) # El user_id es el primer (y único) grupo capturado por la expresión regular print(f"*** Procesando usuario: {user_id} ***") dataframe_local = pd.read_csv(os.path.join(directorio_local, archivo)) #'data/results0.csv') dataframe_local["user_id"] = np.repeat(user_id, dataframe_local.shape[0]) #añadimos el id de usuario como una columna más, tantas veces como filas tengamos de datos dataframe_muse = pd.read_csv(os.path.join(directorio_muse, f'museData{user_id}.csv')) dataframe_muse = dataframe_muse.drop(columns=no_valid_columns) dataframe_muse = dataframe_muse.drop(dataframe_muse[(dataframe_muse['Elements'].apply(lambda x: isinstance(x, str)))].index) #eliminamos aquellas filas en las que solo hay un evento (como pestañear) dataframe_muse = dataframe_muse.drop(columns=['Elements']) X_data_single_subject = [] y_single_subject = [] subjects_single_subject = [] for i, row in dataframe_local.iloc[0:].iterrows(): if row['Letra observada'] != '#': # Si el tiempo de reacción no es demasiado alto tiempo_impulso = datetime.strptime(row['Tiempo de aparición de la letra observada'], '%Y-%m-%d %H:%M:%S.%f') tiempo_pre_decision = tiempo_impulso - timedelta(milliseconds=int(milisegundos)) tiempo_impulso = tiempo_impulso.strftime('%Y-%m-%d %H:%M:%S.%f') tiempo_pre_decision = tiempo_pre_decision.strftime('%Y-%m-%d %H:%M:%S.%f') #Nos quedamos con los valores del intento dataframe_trial = dataframe_muse.drop(dataframe_muse[(dataframe_muse['TimeStamp'] < tiempo_pre_decision) | (dataframe_muse['TimeStamp'] > tiempo_impulso)].index) if len(dataframe_trial) > 0: #Si el valor de las ondas es 'Infinity' o 0.0, ha habido un error en los sensores filtro = (dataframe_trial == 0.0) | (dataframe_trial == float('-inf')) | pd.isna(dataframe_trial) lineas_filtradas = dataframe_trial.loc[filtro.any(axis=1)] #Si no hay celdas con valores invalidos, se procede a la extracción de medidas del intento if len(lineas_filtradas) == 0: n_intentos += 1 user_id = row['ID del participante'] dataframe_normalizado = rescale_data(dataframe_trial.drop("TimeStamp", axis=1)) dataframe_normalizado['id'] = user_id # Se necesita un id para la extracción de características, todos los puntos de datos en el dataframe_trial tienen el mismo id features = extract_features(dataframe_normalizado, column_id='id') _, psd = calculate_psd(dataframe_normalizado.values, nperseg=len(dataframe_normalizado.values)) features['psd'] = psd.flatten() features['differential_entropy'] = calculate_differential_entropy(dataframe_normalizado.values).flatten() features['sample_entropy'] = sample_entropy(dataframe_normalizado.values).flatten() _, _, hjorth_complexity = hjorth_params(dataframe_normalizado.values) features['hjorth'] = hjorth_complexity.flatten() X_data.append(features) X_data_single_subject.append(features) label = 0 if row['Tecla elegida'] == 'p': label = 1 all_labels.append(label) y_single_subject.append(label) subjects.append(user_id) subjects_single_subject.append(user_id) X_data_reshaped_single_subject = np.array(X_data_single_subject).reshape(np.array(X_data_single_subject).shape[0], np.array(X_data_single_subject).shape[2]) df_single_subject = pd.DataFrame(X_data_reshaped_single_subject) df_single_subject_clean = df_single_subject.dropna(axis=1) X_data_single_subject_clean = df_single_subject_clean.to_numpy() ''' sfm_selector = SelectFromModel(estimator=LogisticRegression(max_iter=4000)) sfm_selector.fit(df_single_subject_clean, y_single_subject) columnas_seleccionadas = df_single_subject_clean.columns[sfm_selector.get_support()] df_columnas_seleccionadas_single_subject = df_single_subject_clean[columnas_seleccionadas].copy() X_data_single_subject_selected = df_columnas_seleccionadas_single_subject.to_numpy() ''' os.makedirs(f"{directorio_data_procesed}{user_id}/") np.save(f'{directorio_data_procesed}{user_id}/x.npy', np.array(X_data_single_subject_clean)) #guardamos parcialmente los datos de cada sujeto np.save(f'{directorio_data_procesed}{user_id}/y.npy', np.array(y_single_subject)) #guardamos parcialmente los datos de cada sujeto np.save(f'{directorio_data_procesed}{user_id}/subject_ids.npy', np.array(subjects_single_subject)) #guardamos parcialmente los sujetos de cada sujeto X_data_reshaped = np.array(X_data).reshape(np.array(X_data).shape[0], np.array(X_data).shape[2]) df = pd.DataFrame(X_data_reshaped) df_clean = df.dropna(axis=1) X_data_clean = df_clean.to_numpy() ''' sfm_selector = SelectFromModel(estimator=LogisticRegression(max_iter=50000)) sfm_selector.fit(df_clean, all_labels) columnas_seleccionadas = df_clean.columns[sfm_selector.get_support()] df_columnas_seleccionadas = df_clean[columnas_seleccionadas].copy() X_data_selected = df_columnas_seleccionadas.to_numpy() ''' np.save(f'{directorio_data_procesed}x.npy', np.array(X_data_clean)) np.save(f'{directorio_data_procesed}y.npy', np.array(all_labels)) np.save(f'{directorio_data_procesed}subject_ids.npy', np.array(subjects)) print(f"*** Dimensiones de x_data: {np.array(X_data_clean).shape} ***") print(f"*** Dimensiones de y_data: {np.array(all_labels).shape} ***") print(f"*** Dimensiones de subject_ids: {np.array(subjects).shape} ***") print(f"*** Dimensiones de x_data SINGLE: {np.array(X_data_single_subject_clean).shape} ***") print(f"*** Dimensiones de y_data SINGLE: {np.array(y_single_subject).shape} ***") print(f"*** Numero de intentos: {n_intentos} ***")