diff --git a/.idea/Loop System.iml b/.idea/Loop System.iml new file mode 100644 index 00000000..811de46d --- /dev/null +++ b/.idea/Loop System.iml @@ -0,0 +1,12 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 00000000..dd4c951e --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,7 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 00000000..4502a191 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,15 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 00000000..94a25f7f --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/loop_system/Baseline.py b/loop_system/Baseline.py index 67549921..8aba133a 100644 --- a/loop_system/Baseline.py +++ b/loop_system/Baseline.py @@ -7,12 +7,14 @@ warnings.filterwarnings("ignore") - def select_file(title, filetypes): root = tk.Tk() root.withdraw() # hide the main Tkinter window - file_path = filedialog.askopenfilename(title=title, filetypes=filetypes) + file_path = filedialog.askopenfilename( + title=title, + filetypes=filetypes + ) if not file_path: print("No file selected. Exiting...") @@ -20,12 +22,11 @@ def select_file(title, filetypes): return file_path - def main(): start_time = time.time() # Prompt user to select the XDF file - file_path = select_file("Select the Baseline data file", [("XDF files", "*.xdf")]) + file_path = select_file("Select the Baseline data file",[("XDF files", "*.xdf")]) # Extract folder and condition name from the file path folder = os.path.dirname(file_path) @@ -70,6 +71,5 @@ def main(): except Exception as e: print(f"An error occurred saving the Baseline Dataframe: {e}") - if __name__ == "__main__": main() diff --git a/loop_system/Manager.py b/loop_system/Manager.py index af4b2419..fea7aa70 100644 --- a/loop_system/Manager.py +++ b/loop_system/Manager.py @@ -58,9 +58,7 @@ def run(self): modelTrainer = ModelTrainer() """Load and Select Baseline Data file""" - baseline_file = select_file( - "Select the Baseline data file.", [("CSV files", "*.csv")] - ) + baseline_file = select_file("Select the Baseline data file.", [("CSV files", "*.csv")]) if not baseline_file: print("No baseline file selected. Exiting...") @@ -78,9 +76,7 @@ def run(self): """Load and Select Training Data file""" - training_data = select_file( - "Select the training data file.", [("CSV files", "*.csv")] - ) + training_data = select_file("Select the training data file.", [("CSV files", "*.csv")]) if not training_data: print("No training data file selected. Exiting...") @@ -118,9 +114,7 @@ def run(self): os.makedirs(log_folder, exist_ok=True) # Ensure directory exists log_file_path = os.path.join(log_folder, "output.txt") - print( - "Logging output to {log_file_path}...".format(log_file_path=log_file_path) - ) + print("Logging output to {log_file_path}...".format(log_file_path=log_file_path)) log_file = open(log_file_path, "w") sys.stdout = TeeLogger(log_file) @@ -245,8 +239,11 @@ def run(self): pass finally: try: - model_path = os.path.join(base_dir, f"model_v{self.model_version}.pkl") - joblib.dump(self.model, model_path) + model_path = os.path.join(base_dir, f"model_v{self.model_version}.pkl" ) + joblib.dump( + self.model, + model_path + ) print(f"Model saved successfully to {model_path}") except Exception as e: print(f"Error saving the model: {e}") @@ -265,7 +262,6 @@ def run(self): data_sender.join() print("Closed all processes.") - if __name__ == "__main__": multiprocessing.freeze_support() manager = Manager() diff --git a/loop_system/Training_Model.py b/loop_system/Training_Model.py index d7084558..22a088d7 100644 --- a/loop_system/Training_Model.py +++ b/loop_system/Training_Model.py @@ -9,7 +9,6 @@ warnings.filterwarnings("ignore") - def main(): start_time = time.time() @@ -20,9 +19,7 @@ def main(): """Load Baseline Data file""" - baseline_file = select_file( - "Select the Baseline data file", [("CSV files", "*.csv")] - ) + baseline_file = select_file("Select the Baseline data file", [("CSV files", "*.csv")]) """Baseline Dataframe""" try: @@ -61,6 +58,7 @@ def main(): print(f"An error occurred loading the N-Back data: {e}") sys.exit(1) + """Signals Processing""" signals = getSignals( data, "OpenSignals", "PsychoPy Markers", "PsychoPy Ratings", sensors=sensors @@ -106,6 +104,7 @@ def main(): print("Performing GridSearchCV to find the best models...") best_models = gridSearchCV(X, Y) + print("GridSearch completed.") # Sort the models by their best score in descending order sorted_models = sorted( @@ -138,6 +137,5 @@ def main(): print(f"Elapsed time = {(time.time()-start_time)/60:.2f} minutes") - if __name__ == "__main__": main() diff --git a/physiological_data/Epochs.py b/physiological_data/Epochs.py index 8195059b..a7f0e6a8 100644 --- a/physiological_data/Epochs.py +++ b/physiological_data/Epochs.py @@ -7,7 +7,7 @@ def getMarkers(marker, timestamps): category, onset, offset = list(), list(), list() for timestamp, markers in zip(timestamps, marker): - if markers[1] == "1": + if markers[1] == '1': offset.append(timestamp) elif markers[0] == "end": offset.append(timestamp) diff --git a/physiological_data/Load.py b/physiological_data/Load.py index 6601ea93..b5533922 100644 --- a/physiological_data/Load.py +++ b/physiological_data/Load.py @@ -43,12 +43,8 @@ def Load_Ratings(data, stream_name: str): else: arousal.append("High") - ratings = { - "Valence": valence, - "Arousal": arousal, - "Valence Timestamps": valence_timestamps, - "Arousal Timestamps": arousal_timestamps, - } + ratings = {"Valence": valence, "Arousal": arousal, "Valence Timestamps": valence_timestamps, + "Arousal Timestamps": arousal_timestamps} return ratings diff --git a/physiological_data/PsychoPy.py b/physiological_data/PsychoPy.py index d9cebf62..cbb4cc8c 100644 --- a/physiological_data/PsychoPy.py +++ b/physiological_data/PsychoPy.py @@ -1,184 +1,150 @@ from csv_load import * - def Horror(data): horror = {} for users in sorted(data.keys()): x = pd.DataFrame(columns=["videos"]) for index, row in data[users].iterrows(): - if "Horror" in row["videos"]: - d = { - "videos": row["videos"], - "valence": row["slider_valence.response"], - "arousal": row["slider_arousal.response"], - "dominance": row["slider_dominance.response"], - "anticipatory": row["anticipatory_slider.response"], - "perceived_control": row["perceived_obstacles_slider.response"], - "novelty": row["novelty.response"], - "unexpectedness": row["unexpectedness.response"], - "intrisic_goal": row["intrisic_goal.response"], - "perceived_obstacles": row["perceived_obstacles.response"], - "control_stress": row["control_stress.response"], - } + if ('Horror' in row['videos']): + d = {"videos": row['videos'], 'valence': row['slider_valence.response'], + 'arousal': row['slider_arousal.response'], + 'dominance': row['slider_dominance.response'], 'anticipatory': row['anticipatory_slider.response'], + 'perceived_control': row['perceived_obstacles_slider.response'], + 'novelty': row['novelty.response'], + 'unexpectedness': row['unexpectedness.response'], 'intrisic_goal': row['intrisic_goal.response'], + 'perceived_obstacles': row['perceived_obstacles.response'], + 'control_stress': row['control_stress.response']} x = x.append(d, ignore_index=True) horror[users] = x for users in horror.keys(): - horror[users] = horror[users].sort_values("videos") + horror[users] = horror[users].sort_values('videos') return horror - def Erotic(data): erotic = {} for users in sorted(data.keys()): x = pd.DataFrame(columns=["videos"]) for index, row in data[users].iterrows(): - if "Erotic" in row["videos"]: - d = { - "videos": row["videos"], - "valence": row["slider_valence.response"], - "arousal": row["slider_arousal.response"], - "dominance": row["slider_dominance.response"], - "anticipatory": row["anticipatory_slider.response"], - "perceived_control": row["perceived_obstacles_slider.response"], - "novelty": row["novelty.response"], - "unexpectedness": row["unexpectedness.response"], - "intrisic_goal": row["intrisic_goal.response"], - "perceived_obstacles": row["perceived_obstacles.response"], - "control_stress": row["control_stress.response"], - } + if ('Erotic' in row['videos']): + d = {"videos": row['videos'], 'valence': row['slider_valence.response'], + 'arousal': row['slider_arousal.response'], + 'dominance': row['slider_dominance.response'], 'anticipatory': row['anticipatory_slider.response'], + 'perceived_control': row['perceived_obstacles_slider.response'], + 'novelty': row['novelty.response'], + 'unexpectedness': row['unexpectedness.response'], 'intrisic_goal': row['intrisic_goal.response'], + 'perceived_obstacles': row['perceived_obstacles.response'], + 'control_stress': row['control_stress.response']} x = x.append(d, ignore_index=True) erotic[users] = x for users in erotic.keys(): - erotic[users] = erotic[users].sort_values("videos") + erotic[users] = erotic[users].sort_values('videos') return erotic - def Scenery(data): scenery = {} for users in sorted(data.keys()): x = pd.DataFrame(columns=["videos"]) for index, row in data[users].iterrows(): - if "Scenery" in row["videos"]: - d = { - "videos": row["videos"], - "valence": row["slider_valence.response"], - "arousal": row["slider_arousal.response"], - "dominance": row["slider_dominance.response"], - "anticipatory": row["anticipatory_slider.response"], - "perceived_control": row["perceived_obstacles_slider.response"], - "novelty": row["novelty.response"], - "unexpectedness": row["unexpectedness.response"], - "intrisic_goal": row["intrisic_goal.response"], - "perceived_obstacles": row["perceived_obstacles.response"], - "control_stress": row["control_stress.response"], - } + if ('Scenery' in row['videos']): + d = {"videos": row['videos'], 'valence': row['slider_valence.response'], + 'arousal': row['slider_arousal.response'], + 'dominance': row['slider_dominance.response'], 'anticipatory': row['anticipatory_slider.response'], + 'perceived_control': row['perceived_obstacles_slider.response'], + 'novelty': row['novelty.response'], + 'unexpectedness': row['unexpectedness.response'], 'intrisic_goal': row['intrisic_goal.response'], + 'perceived_obstacles': row['perceived_obstacles.response'], + 'control_stress': row['control_stress.response']} x = x.append(d, ignore_index=True) scenery[users] = x for users in scenery.keys(): - scenery[users] = scenery[users].sort_values("videos") + scenery[users] = scenery[users].sort_values('videos') return scenery - def Social_positive(data): social_positive = {} for users in sorted(data.keys()): x = pd.DataFrame(columns=["videos"]) for index, row in data[users].iterrows(): - if "Social Positive" in row["videos"]: - d = { - "videos": row["videos"], - "valence": row["slider_valence.response"], - "arousal": row["slider_arousal.response"], - "dominance": row["slider_dominance.response"], - "anticipatory": row["anticipatory_slider.response"], - "perceived_control": row["perceived_obstacles_slider.response"], - "novelty": row["novelty.response"], - "unexpectedness": row["unexpectedness.response"], - "intrisic_goal": row["intrisic_goal.response"], - "perceived_obstacles": row["perceived_obstacles.response"], - "control_stress": row["control_stress.response"], - } + if ('Social Positive' in row['videos']): + d = {"videos": row['videos'], 'valence': row['slider_valence.response'], + 'arousal': row['slider_arousal.response'], + 'dominance': row['slider_dominance.response'], 'anticipatory': row['anticipatory_slider.response'], + 'perceived_control': row['perceived_obstacles_slider.response'], + 'novelty': row['novelty.response'], + 'unexpectedness': row['unexpectedness.response'], 'intrisic_goal': row['intrisic_goal.response'], + 'perceived_obstacles': row['perceived_obstacles.response'], + 'control_stress': row['control_stress.response']} x = x.append(d, ignore_index=True) social_positive[users] = x for users in social_positive.keys(): - social_positive[users] = social_positive[users].sort_values("videos") + social_positive[users] = social_positive[users].sort_values('videos') return social_positive - def Social_negative(data): social_negative = {} for users in sorted(data.keys()): x = pd.DataFrame(columns=["videos"]) for index, row in data[users].iterrows(): - if "Social Negative" in row["videos"]: - d = { - "videos": row["videos"], - "valence": row["slider_valence.response"], - "arousal": row["slider_arousal.response"], - "dominance": row["slider_dominance.response"], - "anticipatory": row["anticipatory_slider.response"], - "perceived_control": row["perceived_obstacles_slider.response"], - "novelty": row["novelty.response"], - "unexpectedness": row["unexpectedness.response"], - "intrisic_goal": row["intrisic_goal.response"], - "perceived_obstacles": row["perceived_obstacles.response"], - "control_stress": row["control_stress.response"], - } + if ('Social Negative' in row['videos']): + d = {"videos": row['videos'], 'valence': row['slider_valence.response'], + 'arousal': row['slider_arousal.response'], + 'dominance': row['slider_dominance.response'], 'anticipatory': row['anticipatory_slider.response'], + 'perceived_control': row['perceived_obstacles_slider.response'], + 'novelty': row['novelty.response'], + 'unexpectedness': row['unexpectedness.response'], 'intrisic_goal': row['intrisic_goal.response'], + 'perceived_obstacles': row['perceived_obstacles.response'], + 'control_stress': row['control_stress.response']} x = x.append(d, ignore_index=True) social_negative[users] = x for users in social_negative.keys(): - social_negative[users] = social_negative[users].sort_values("videos") + social_negative[users] = social_negative[users].sort_values('videos') return social_negative - -def Category_Dataframe(category, dimension): +def Category_Dataframe(category,dimension): df1 = pd.DataFrame(columns=[category["P10_S1"]["videos"]]) df2 = pd.DataFrame(columns=[category["P10_S2"]["videos"]]) users_list = list() for users in category.keys(): - if "EMDB/A" in str(category[users]["videos"]): + if ("EMDB/A" in str(category[users]["videos"])): users_list.append(users.split("_")[0]) df1.loc[len(df1)] = category[users][dimension].values - if "EMDB/B" in str(category[users]["videos"]): + if ("EMDB/B" in str(category[users]["videos"])): df2.loc[len(df2)] = category[users][dimension].values df1.insert(0, "Users", users_list) df = (df1.join(df2)).set_index("Users") return df - -def Full_Dataframe(data, dimension): +def Full_Dataframe(data,dimension): horror = Horror(data) erotic = Erotic(data) scenery = Scenery(data) social_positive = Social_positive(data) social_negative = Social_negative(data) - df_erotic = Category_Dataframe(erotic, dimension) - df_horror = Category_Dataframe(horror, dimension) - df_scenery = Category_Dataframe(scenery, dimension) - df_positive = Category_Dataframe(social_positive, dimension) - df_negative = Category_Dataframe(social_negative, dimension) + df_erotic = Category_Dataframe(erotic,dimension) + df_horror = Category_Dataframe(horror,dimension) + df_scenery = Category_Dataframe(scenery,dimension) + df_positive = Category_Dataframe(social_positive,dimension) + df_negative = Category_Dataframe(social_negative,dimension) - df = (((df_erotic.join(df_horror)).join(df_negative)).join(df_scenery)).join( - df_positive - ) + df = (((df_erotic.join(df_horror)).join(df_negative)).join(df_scenery)).join(df_positive) return df - # path = 'G:\\O meu disco\\PhD\\1st Study\\PsychoPy Data\\' # data = CSV_Load(path) # horror = Horror(data) + diff --git a/physiological_data/csv_load.py b/physiological_data/csv_load.py index 8b291cd5..403ccc34 100644 --- a/physiological_data/csv_load.py +++ b/physiological_data/csv_load.py @@ -10,4 +10,4 @@ def CSV_Load(path): for files in glob.glob("*.csv"): data[files.split(".")[0]] = pd.read_csv(files) - return data + return data \ No newline at end of file diff --git a/physiological_data/lib/biosignals.py b/physiological_data/lib/biosignals.py index 5cd9f270..a8658f4f 100644 --- a/physiological_data/lib/biosignals.py +++ b/physiological_data/lib/biosignals.py @@ -65,20 +65,22 @@ def readHeader(path): def convertAndfilterSignal(data, sensor, fs, resolution): if sensor == "ECG": low_f, high_f = 5, 15 - data = ((data / (2**resolution)) - 1 / 2) * 3000 / 1019 # millivolts + data = ((data / (2 ** resolution)) - 1 / 2) * 3000 / 1019 # millivolts elif sensor == "hSpO2": low_f, high_f = 0.05, 1 - data = (0.15 * data) / (2**resolution) # microAmpere + data = (0.15 * data) / (2 ** resolution) # microAmpere elif sensor == "EDA": low_f, high_f = 0.016, 35 - data = 1e6 * (data / (2**resolution)) * 3 / 0.12 # microSiemens + data = 1e6 * (data / (2 ** resolution)) * 3 / 0.12 # microSiemens elif sensor == "EEG": # CONFIRM low_f, high_f = 1, 30 - data = 1e6 * (data / ((2**resolution) - 1) - 0.5) * 3 / 40000 # microVolts + data = ( + 1e6 * (data / ((2 ** resolution) - 1) - 0.5) * 3 / 40000 + ) # microVolts elif sensor == "RESPIRATION": low_f, high_f = 0.01, 1 - data = ((data / (2**resolution)) - 1 / 2) * 100 # % of displacement + data = ((data / (2 ** resolution)) - 1 / 2) * 100 # % of displacement elif sensor == "XYZ": # CONFIRM - Movements happening at more than 1 per second is very unlikely... low_f, high_f = 0.01, 10 diff --git a/physiological_data/lib/biosignals_classifier.py b/physiological_data/lib/biosignals_classifier.py index 2578084f..9e68965c 100644 --- a/physiological_data/lib/biosignals_classifier.py +++ b/physiological_data/lib/biosignals_classifier.py @@ -19,18 +19,12 @@ def __init__(self, **kwargs): super().__init__() self._features = None - self.encoder_hidden_layer = nn.Linear( - in_features=kwargs["input_shape"], out_features=128 - ) + self.encoder_hidden_layer = nn.Linear(in_features=kwargs["input_shape"], out_features=128) self.encoder_output_layer = nn.Linear(in_features=128, out_features=128) self.decoder_hidden_layer = nn.Linear(in_features=128, out_features=128) - self.decoder_output_layer = nn.Linear( - in_features=128, out_features=kwargs["input_shape"] - ) + self.decoder_output_layer = nn.Linear(in_features=128, out_features=kwargs["input_shape"]) - self.encoder_output_layer.register_forward_hook( - self.save_outputs_hook(self.encoder_output_layer) - ) + self.encoder_output_layer.register_forward_hook(self.save_outputs_hook(self.encoder_output_layer)) def forward(self, features): activation = self.encoder_hidden_layer(features) @@ -46,7 +40,6 @@ def forward(self, features): def save_outputs_hook(self, layer_id): def fn(_, __, output): self._features = output.detach() - return fn @@ -112,15 +105,14 @@ def train_classifier(epochs, train_loader): return model, features, labels -if __name__ == "__main__": - groups = [ - "05_11_2020_1", - "05_11_2020_2", - "10_11_2020_2", - "27_11_2020", - "11_12_2020_1", - "11_12_2020_2", - ] +if __name__ == '__main__': + groups = ['05_11_2020_1', + '05_11_2020_2', + '10_11_2020_2', + '27_11_2020', + '11_12_2020_1', + '11_12_2020_2' + ] data = SensorsDataset(sensors=[1], groups=groups[:-1]) diff --git a/physiological_data/lib/classification.py b/physiological_data/lib/classification.py index 0bda224c..955316c0 100644 --- a/physiological_data/lib/classification.py +++ b/physiological_data/lib/classification.py @@ -3,28 +3,22 @@ except (ImportError, ModuleNotFoundError): from lib.acquisition import * -folders = [ - "05_11_2020_1", - "05_11_2020_2", - "10_11_2020_2", - "27_11_2020", - "11_12_2020_1", - "11_12_2020_2", -] +folders = ['05_11_2020_1', + '05_11_2020_2', + '10_11_2020_2', + '27_11_2020', + '11_12_2020_1', + '11_12_2020_2' + ] features, labels, participants = [], [], [] for folder in folders: - path = os.path.join("..\..\acquisitions\Acquisitions", folder) - acquisition = Acquisition(r"..\acquisitions\Acquisitions\11_12_2020_1") - timestamps, initialLabels = ( - acquisition.getPsychoTimestamps(), - acquisition.getPsychoActivity(), - ) + path = os.path.join('..\..\acquisitions\Acquisitions', folder) + acquisition = Acquisition(r'..\acquisitions\Acquisitions\11_12_2020_1') + timestamps, initialLabels = acquisition.getPsychoTimestamps(), acquisition.getPsychoActivity() acquisition.getBiosignalsSensors() - acquisition.segmentWindowingBiosignals( - acquisition.signal, timestamps, initialLabels - ) + acquisition.segmentWindowingBiosignals(acquisition.signal, timestamps, initialLabels) acquisition.extractFeatures(acquisition.segmentedBiosignals) if len(features) == 0: features, labels, participants = acquisition.getDataset() @@ -39,14 +33,8 @@ from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import accuracy_score, f1_score - ################# Classification ################# -from sklearn.model_selection import ( - train_test_split, - ShuffleSplit, - KFold, - LeaveOneGroupOut, -) +from sklearn.model_selection import train_test_split, ShuffleSplit, KFold, LeaveOneGroupOut # model = KFold(n_splits=10, shuffle=True, random_state=42) model = LeaveOneGroupOut() @@ -58,9 +46,7 @@ labels_train, labels_test = labels[train_index], labels[test_index] # Build the random forest clasdsifier. - random_forest = RandomForestClassifier( - n_estimators=100, criterion="entropy", n_jobs=-1, random_state=42 - ) + random_forest = RandomForestClassifier(n_estimators=100, criterion='entropy', n_jobs=-1, random_state=42) # Train the classifier on the training set. random_forest = random_forest.fit(samples_train, labels_train) @@ -70,14 +56,8 @@ # This step is not necessary for the classification procedure, but is important to store the values # of accuracy to calculate the mean and standard deviation values and evaluate the performance of the classifier. - acc.append(accuracy_score(labels_test, results) * 100) - f1.append( - f1_score( - [0 if x == "task" else 1 for x in labels_test], - [0 if x == "task" else 1 for x in results], - ) - * 100 - ) + acc.append(accuracy_score(labels_test, results)*100) + f1.append(f1_score([0 if x =='task' else 1 for x in labels_test], [0 if x =='task' else 1 for x in results])*100) print(f"Accuracy: {np.mean(acc):.2f} +- {np.std(acc):.2f}%") print(f"F1-Score: {np.mean(f1):.2f} +- {np.std(f1):.2f}%") diff --git a/physiological_data/lib/dataLoader.py b/physiological_data/lib/dataLoader.py index 7fc58e28..1b79afe3 100644 --- a/physiological_data/lib/dataLoader.py +++ b/physiological_data/lib/dataLoader.py @@ -6,27 +6,20 @@ class SensorsDataset(Dataset): - def __init__( - self, - numpy_file="..\\datasets\\dataset_1_024s.npy", - sensors=[1, 2, 3], - groups=[1], - normalize=True, - quantize_data=True, - ): + def __init__(self, numpy_file='..\\datasets\\dataset_1_024s.npy', sensors=[1,2,3], groups=[1], normalize=True, quantize_data=True): # Segments, label, group if type(sensors) == int: sensors = list(sensors) if type(numpy_file) == str: - print(f"Reading {numpy_file} file...", "\r") + print(f"Reading {numpy_file} file...", '\r') data = np.load(numpy_file) - print(f"Finished loading {numpy_file} file!", "\r") + print(f"Finished loading {numpy_file} file!", '\r') else: - num = str(numpy_file).replace(".", "_") - print(f"Reading .\\datasets\\dataset_{num}s.npy file...", end="\r") - data = np.load(f".\\datasets\\dataset_{num}s.npy") - print(f"Finished loading .\\datasets\\dataset_{num}s.npy file!", end="\r") + num = str(numpy_file).replace('.', '_') + print(f"Reading .\\datasets\\dataset_{num}s.npy file...", end='\r') + data = np.load(f'.\\datasets\\dataset_{num}s.npy') + print(f"Finished loading .\\datasets\\dataset_{num}s.npy file!", end='\r') if type(groups) == str: groups = [groups] @@ -41,7 +34,7 @@ def __init__( all_data = [] num_cols = data.shape[1] // 12 for s in sensors: - aux = data[samples, (s - 1) * num_cols : s * num_cols] + aux = data[samples, (s-1) * num_cols:s * num_cols] if len(all_data) == 0: all_data = aux else: @@ -49,20 +42,16 @@ def __init__( all_data = np.array(all_data, dtype=np.float) if quantize_data: - self.normalized_data = np.array( - normalizeFeatures(all_data)[0], dtype=np.float - ) - self.quantized_data = np.array( - quantize(self.normalized_data), dtype=np.float - ) + self.normalized_data = np.array(normalizeFeatures(all_data)[0], dtype=np.float) + self.quantized_data = np.array(quantize(self.normalized_data), dtype=np.float) if normalize and not quantize_data: self.normalized_data, _, _ = normalizeFeatures(all_data) - + self.labels = data[samples, -2] self.groups = data[samples, -1] - + def __len__(self): return len(self.labels) - + def __getitem__(self, idx): return self.quantized_data[idx], self.labels[idx], self.groups[idx] diff --git a/physiological_data/lib/latent.py b/physiological_data/lib/latent.py index 1f64fcc0..7113975a 100644 --- a/physiological_data/lib/latent.py +++ b/physiological_data/lib/latent.py @@ -9,11 +9,11 @@ from scipy.io import wavfile from scipy.signal import spectrogram -AUDIO = "audio.wav" -KEYBOARD = "keyboard.csv" -MOUSE = "mouse.csv" -SCREENSHOT = "screenshot.csv" -SNAPSHOT = "snapshot.csv" +AUDIO = 'audio.wav' +KEYBOARD = 'keyboard.csv' +MOUSE = 'mouse.csv' +SCREENSHOT = 'screenshot.csv' +SNAPSHOT = 'snapshot.csv' files = [AUDIO, KEYBOARD, MOUSE, SCREENSHOT, SNAPSHOT] @@ -25,15 +25,15 @@ def __init__(self, pathFolder, fs): def getPath(self, pathFolder): for folder in listdir(pathFolder): - if "ip" in folder: + if 'ip' in folder: return join(pathFolder, folder) def readLatentFile(self, file): end = file[-4:] - if ".csv" == end: + if '.csv' == end: data = read_csv(join(self.path, file)) return data - elif ".wav" == end: + elif '.wav' == end: sampling_ratem, data = wavfile.read(join(self.path, AUDIO)) return data, sampling_ratem return [] @@ -41,16 +41,16 @@ def readLatentFile(self, file): def readData(self): for file in files: if AUDIO != file: - self.data[file.split(".")[0]] = self.readLatentFile(file) + self.data[file.split('.')[0]] = self.readLatentFile(file) else: - self.data[file.split(".")[0]], self.audioFS = self.readLatentFile(file) + self.data[file.split('.')[0]], self.audioFS = self.readLatentFile(file) @staticmethod def makeSpectrogram(signal, fs): f, t, Sxx = spectrogram(signal, fs) - plt.pcolormesh(t, f, Sxx, shading="gouraud") - plt.ylabel("Frequency [Hz]") - plt.xlabel("Time [sec]") + plt.pcolormesh(t, f, Sxx, shading='gouraud') + plt.ylabel('Frequency [Hz]') + plt.xlabel('Time [sec]') plt.show() return f, t, Sxx @@ -61,57 +61,32 @@ def interpolateInterfaces(self, interface): columns, keys = self.identifyTypeInterp(interface) # time = np.arange(self.data[interface].time[0], self.data[interface].time[1]+1, T) time_DF = np.array(self.data[interface].time) - time = np.linspace( - time_DF[0], time_DF[-1], int(np.round((time_DF[-1] - time_DF[0]) * self.fs)) - ) + time = np.linspace(time_DF[0], time_DF[-1], int(np.round((time_DF[-1] - time_DF[0]) * self.fs))) time[0] = self.data[interface].time[0] data = np.empty((len(time), len(self.data[interface].keys())), dtype=np.object) data[:] = np.nan data[:, 0] = time - data = self.duplicateValues( - self.data[interface], - interface, - data, - time_DF, - keys["duplicate"], - columns["duplicate"], - ) - data = self.fillValues( - self.data[interface], - interface, - data, - time_DF, - keys["fill"] + keys["interpolate"], - columns["fill"] + columns["interpolate"], - ) - data = self.interpolateValues( - self.data[interface], - interface, - data, - time_DF, - keys["interpolate"], - columns["interpolate"], - ) + data = self.duplicateValues(self.data[interface], interface, data, time_DF, keys['duplicate'], columns['duplicate']) + data = self.fillValues(self.data[interface], interface, data, time_DF, keys['fill']+keys['interpolate'], columns['fill'] + columns['interpolate']) + data = self.interpolateValues(self.data[interface], interface, data, time_DF, keys['interpolate'], columns['interpolate']) return data def duplicateValues(self, duplicate_data, interface, new_data, data, keys, columns): time = new_data[:, 0] duplicate_index = 0 for i in range(len(time)): - print(f"Duplicating {interface} Data: {i*100/len(time):.2f}%", end="\r") + print(f"Duplicating {interface} Data: {i*100/len(time):.2f}%", end='\r') try: - while np.abs(data[duplicate_index + 1] - time[i]) < np.abs( - data[duplicate_index] - time[i] - ): + while np.abs(data[duplicate_index + 1] - time[i]) < np.abs(data[duplicate_index] - time[i]): duplicate_index += 1 except IndexError: duplicate_index = len(data) - 1 # duplicate_index += closestIndex(data[duplicate_index:], time[i]) for j, key in enumerate(keys): new_data[i, columns[j]] = duplicate_data[key][duplicate_index] - print(" ") + print(' ') return new_data def fillValues(self, fill_data, interface, new_data, data, keys, columns): @@ -119,11 +94,9 @@ def fillValues(self, fill_data, interface, new_data, data, keys, columns): time = new_data[:, 0] fill_index = 0 for i, value in enumerate(data): - print(f"Filling {interface} Data: {i*100/len(data):.2f}%", end="\r") + print(f"Filling {interface} Data: {i*100/len(data):.2f}%", end='\r') try: - while np.abs(time[fill_index + 1] - value) < np.abs( - time[fill_index] - value - ): + while np.abs(time[fill_index + 1] - value) < np.abs(time[fill_index] - value): fill_index += 1 except IndexError: fill_index = len(time) - 1 @@ -131,108 +104,64 @@ def fillValues(self, fill_data, interface, new_data, data, keys, columns): for j, key in enumerate(keys): new_data[fill_index, columns[j]] = fill_data[key][i] prev_fill_index = fill_index - print(" ") + print(' ') return new_data - - def interpolateValues( - self, interpolate_data, interface, new_data, data, keys, columns - ): + + def interpolateValues(self, interpolate_data, interface, new_data, data, keys, columns): prev_index, index = -1, 0 time = new_data[:, 0] not_nan_indexes = self.not_nan(new_data[:, columns[1]]) - T = 1 / self.fs + T = 1/self.fs m, b, x0, x1, y0, y1 = 0, 0, 0, 0, 0, 0 for j, key in enumerate(columns[1:]): for i in range(1, len(new_data)): - print( - f"Interpolating {interface} Data: {i*100/len(new_data):.2f}%", - end="\r", - ) - + print(f"Interpolating {interface} Data: {i*100/len(new_data):.2f}%", end='\r') + if np.isnan(new_data[i, key]): new_value = m * time[i] + b new_data[i, key] = new_value else: if index + 1 < len(not_nan_indexes): x0 = time[not_nan_indexes[index]] - x1 = time[not_nan_indexes[index + 1]] + x1 = time[not_nan_indexes[index+1]] y0 = new_data[not_nan_indexes[index], key] - y1 = new_data[not_nan_indexes[index + 1], key] + y1 = new_data[not_nan_indexes[index+1], key] m = (y1 - y0) / (x1 - x0) b = y0 - m * x0 - + try: prev_index = index if i >= not_nan_indexes[index]: index += 1 except IndexError: index = len(not_nan_indexes) - 1 - print(" ") + print(' ') return new_data - + @staticmethod def not_nan(data): index = [] for i, value in enumerate(data): if not np.isnan(value): index.append(i) - return index + return index def identifyTypeInterp(self, interface): - if interface.lower() == "mouse": - return { - "interpolate": [0, 7, 8, 9, 10, 14], - "duplicate": [1, 2, 13], - "fill": [3, 4, 5, 6, 11, 12], - }, { - "interpolate": [ - "time", - "screen_x", - "screen_y", - "page_x", - "page_y", - "event_timestamp", - ], - "duplicate": ["type", "tab_id", "xpath"], - "fill": [ - "button", - "shift_key", - "alt_key", - "ctrl_key", - "delta_y", - "scrolltop", - ], - } - elif interface.lower() == "keyboard": - return { - "interpolate": [0, 8], - "duplicate": [1, 2], - "fill": [3, 4, 5, 6, 7], - }, { - "interpolate": ["time", "event_timestamp"], - "duplicate": ["type", "tab_id"], - "fill": [ - "key_code", - "shift_key", - "alt_key", - "ctrl_key", - "capslock_key", - ], - } - - -if __name__ == "__main__": - latent = Latent( - r"D:\Google Drive\Faculdade\Doutoramento\Acquisitions\new_biosignalsnotebooks\acquisitions\Acquisitions\03_11_2020", - 100, - ) + if interface.lower() == 'mouse': + return {'interpolate': [0, 7, 8, 9, 10, 14], 'duplicate': [1, 2, 13], 'fill': [3, 4, 5, 6, 11, 12]}, {'interpolate': ['time', 'screen_x', 'screen_y', 'page_x', 'page_y', 'event_timestamp'], 'duplicate': ['type', 'tab_id', 'xpath'], 'fill': ['button', 'shift_key', 'alt_key', 'ctrl_key', 'delta_y', 'scrolltop']} + elif interface.lower() == 'keyboard': + return {'interpolate': [0, 8], 'duplicate': [1, 2], 'fill': [3, 4, 5, 6, 7]}, {'interpolate': ['time', 'event_timestamp'], 'duplicate': ['type', 'tab_id'], 'fill': ['key_code', 'shift_key', 'alt_key', 'ctrl_key', 'capslock_key']} + +if __name__ == '__main__': + latent = Latent(r"D:\Google Drive\Faculdade\Doutoramento\Acquisitions\new_biosignalsnotebooks\acquisitions\Acquisitions\03_11_2020", 100) latent.readData() - data = latent.interpolateInterfaces("keyboard") + data = latent.interpolateInterfaces('keyboard') print(data[:, 7]) plt.figure() - plt.scatter(latent.data["keyboard"].time, latent.data["keyboard"].key_code) - plt.scatter(data[:, 0], data[:, 3], color="r") + plt.scatter(latent.data['keyboard'].time, latent.data['keyboard'].key_code) + plt.scatter(data[:, 0], data[:, 3], color='r') plt.show() # latent.interpolateInterfaces(1000, 'keyboard') # data, fs = latent.readLatentFile(AUDIO) # latent.makeSpectrogram(data[::1000, 0], fs/1000) + diff --git a/physiological_data/lib/main.py b/physiological_data/lib/main.py index afd2a24f..3c302ce4 100644 --- a/physiological_data/lib/main.py +++ b/physiological_data/lib/main.py @@ -9,8 +9,8 @@ import matplotlib.pylab as plt -psycho = Psycho(r"..\..\acquisitions\Acquisitions\11_12_2020_1\results_3.csv") -device = Devices(r"..\..\acquisitions\Acquisitions\11_12_2020_1") +psycho = Psycho(r'..\..\acquisitions\Acquisitions\11_12_2020_1\results_3.csv') +device = Devices(r'..\..\acquisitions\Acquisitions\11_12_2020_1') timestamps = psycho.getTimestamps() labels = psycho.getActivity() diff --git a/physiological_data/lib/novainstrumentation/ARparameters.py b/physiological_data/lib/novainstrumentation/ARparameters.py index c173b223..7c772f5e 100644 --- a/physiological_data/lib/novainstrumentation/ARparameters.py +++ b/physiological_data/lib/novainstrumentation/ARparameters.py @@ -1,44 +1,46 @@ from scipy import signal - def lpcar2cc(ar): sh = shape(ar) - - if len(sh) == 1: - p1 = len(ar) - nf = 1 - b = zeros(len(ar)) + + if len(sh) == 1 : + p1 = len(ar) + nf = 1 + b = zeros(len(ar)) else: - p1 = len(ar[0]) - nf = len(ar) - b = zeros(len(ar[0])) - p = p1 - 1 + p1 = len(ar[0]) + nf = len(ar) + b = zeros(len(ar[0])) + p = p1-1 - cm = arange(1, p + 1) ** (-1.0) + cm = arange(1,p+1)**(-1.) - xm = -arange(1, p + 1) + xm = -arange(1,p+1) + b[0] = 1 cc = [] - + for k in xrange(nf): if nf > 1: - cc += [signal.lfilter(b, ar[k], ar[k, 1:p1] * xm) * cm] + cc += [ signal.lfilter(b,ar[k],ar[k,1:p1]*xm)*cm ] else: - cc = signal.lfilter(b, ar, ar[1:p1] * xm) * cm + cc = signal.lfilter(b,ar,ar[1:p1]*xm)*cm return cc - - -def lpc_coef(s, coefNumber): - + +def lpc_coef(s,coefNumber): + from scikits.talkbox import * - - # you have to install scikits.talkbox library before you can use it. - + + #you have to install scikits.talkbox library before you can use it. + order = coefNumber - 1 - - est, e, k = lpc(s, order) - + + est,e,k = lpc(s,order) + return est + + + \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/__init__.py b/physiological_data/lib/novainstrumentation/__init__.py index bffcf2dc..6d28563a 100644 --- a/physiological_data/lib/novainstrumentation/__init__.py +++ b/physiological_data/lib/novainstrumentation/__init__.py @@ -14,4 +14,4 @@ def peakdelta(first_derivative, param): - return None + return None \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/figtools.py b/physiological_data/lib/novainstrumentation/figtools.py index 59be6f62..1650be29 100644 --- a/physiological_data/lib/novainstrumentation/figtools.py +++ b/physiological_data/lib/novainstrumentation/figtools.py @@ -4,77 +4,61 @@ import tkinter.filedialog as tkFileDialog import pandas - # added import - David Belo from matplotlib import pyplot as plot from numpy import arange -from pylab import ( - rc, - close, - figure, - axes, - subplot, - plot, - axis, - show, - grid, - savefig, - text, -) +from pylab import rc, close, figure, axes, subplot, plot, axis, show, grid, savefig, text def load_data_dialog(path): - # root = Tkinter.Tk() - # root.withdraw() - + #root = Tkinter.Tk() + #root.withdraw() + # Make it almost invisible - no decorations, 0 size, top left corner. - # root.overrideredirect(True) - # root.geometry('0x0+0+0') - + #root.overrideredirect(True) + #root.geometry('0x0+0+0') + # Show window again and lift it to top so it can get focus, # otherwise dialogs will end up behind the terminal. - # root.deiconify() - # root.lift() - # root.focus_force() - - # filename = tkFileDialog.askopenfile(parent=root) # Or some other dialog - filename = tkFileDialog.askopenfile() # Or some other dialog - + #root.deiconify() + #root.lift() + #root.focus_force() + + #filename = tkFileDialog.askopenfile(parent=root) # Or some other dialog + filename = tkFileDialog.askopenfile() # Or some other dialog + # Get rid of the top-level instance once to make it actually invisible. - # root.destroy() - return pandas.read_csv(filename, sep=" ", header=None) - + #root.destroy() + return pandas.read_csv(filename, sep=' ', header=None) + ########## -# Initial configurations +#Initial configurations def pylabconfig(): - rc("lines", linewidth=2, color="k") - # rc('lines', linewidth=1, color='k') + rc('lines', linewidth=2, color='k') + #rc('lines', linewidth=1, color='k') - rc("font", **{"family": "serif", "serif": ["Palatino"]}) - rc("font", style="italic", size=10) + rc('font', **{'family': 'serif', 'serif': ['Palatino']}) + rc('font', style='italic', size=10) - rc("text", color="grey") + rc('text', color='grey') # rc('text', usetex=True) - rc("text", usetex=False) + rc('text', usetex=False) + + rc('figure', figsize=(8, 5), dpi=80) + rc('axes', grid=True, edgecolor='grey', labelsize=10,) + rc('grid', color='grey') + rc('xtick', color='grey', labelsize=10) + rc('ytick', color='grey', labelsize=10) - rc("figure", figsize=(8, 5), dpi=80) - rc( - "axes", - grid=True, - edgecolor="grey", - labelsize=10, - ) - rc("grid", color="grey") - rc("xtick", color="grey", labelsize=10) - rc("ytick", color="grey", labelsize=10) + close('all') - close("all") + def plotwithhist(t, s, bins=50): from matplotlib.ticker import NullFormatter @@ -86,7 +70,8 @@ def plotwithhist(t, s, bins=50): ax1.plot(t, s) ax1.set_xticks(ax1.get_xticks()[:-1]) - ax2.hist(s, bins, normed=True, facecolor="white", orientation="horizontal", lw=2) + ax2.hist(s, bins, normed=True, facecolor='white', + orientation='horizontal', lw=2) ax2.axis([0, 1, ax1.axis()[2], ax1.axis()[3]]) @@ -95,11 +80,11 @@ def plotwithhist(t, s, bins=50): return ax1, ax2 - ########### + ########### def plotwithstats(t, s): - + from matplotlib.ticker import NullFormatter nullfmt = NullFormatter() @@ -116,9 +101,9 @@ def plotwithstats(t, s): mx = s.max() sd = s.std() - ax2.bar(-0.5, mx - mi, 1, mi, lw=2, color="#f0f0f0") - ax2.bar(-0.5, sd * 2, 1, meanv - sd, lw=2, color="#c0c0c0") - ax2.bar(-0.5, 0.2, 1, meanv - 0.1, lw=2, color="#b0b0b0") + ax2.bar(-0.5, mx - mi, 1, mi, lw=2, color='#f0f0f0') + ax2.bar(-0.5, sd * 2, 1, meanv - sd, lw=2, color='#c0c0c0') + ax2.bar(-0.5, 0.2, 1, meanv - 0.1, lw=2, color='#b0b0b0') ax2.axis([-1, 1, ax1.axis()[2], ax1.axis()[3]]) ax2.yaxis.set_major_formatter(nullfmt) @@ -127,42 +112,43 @@ def plotwithstats(t, s): return ax1, ax2 -def multilineplot(signal, linesize=250, events=None, title="", dir="", step=1): +def multilineplot(signal, linesize=250, events=None, title = '', dir = '', step=1): from pylab import rc + rc('axes', labelcolor='#a1a1a1', edgecolor='#a1a1a1', labelsize='xx-small') + rc('xtick',color='#a1a1a1', labelsize='xx-small') - rc("axes", labelcolor="#a1a1a1", edgecolor="#a1a1a1", labelsize="xx-small") - rc("xtick", color="#a1a1a1", labelsize="xx-small") - - grid("off") - nplots = len(signal) // linesize + 1 + + grid('off') + nplots=len(signal)//linesize + 1 ma_x = max(signal) mi_x = min(signal) - f = figure(figsize=(20, 1.5 * nplots), dpi=80) + f=figure(figsize=(20, 1.5*nplots), dpi=80) for i in range(nplots): - ax = subplot(nplots, 1, i + 1) - - start = i * len(signal) / nplots - end = (i + 1) * len(signal) / nplots - plot(arange(start, end, step), signal[start:end:step], "k") - axis((start, end, mi_x, ma_x)) + ax=subplot(nplots,1,i+1) + + start = i*len(signal)/nplots + end = (i+1)*len(signal)/nplots + plot(arange(start,end,step),signal[start:end:step],'k') + axis((start,end,mi_x,ma_x)) ax.set_yticks([]) ax.set_xticks(ax.get_xticks()[1:-1]) - ax.spines["right"].set_visible(False) - ax.spines["top"].set_visible(False) - ax.spines["left"].set_visible(False) - - # ax.spines['bottom'].set_visible(False) - # ax.spines['left'].set_bounds(mi_x,ma_x) - ax.xaxis.set_ticks_position("bottom") + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.spines['left'].set_visible(False) + + #ax.spines['bottom'].set_visible(False) + #ax.spines['left'].set_bounds(mi_x,ma_x) + ax.xaxis.set_ticks_position('bottom') if events != None: e = events[(events >= start) & (events < end)] - if len(e) > 0: - plot.vlines(e, mi_x, ma_x - (ma_x - mi_x) / 4.0 * 3.0, lw=2) + if len(e)>0: + plot.vlines(e,mi_x,ma_x-(ma_x-mi_x)/4.*3., lw=2) if title != None: - text(start, ma_x, title) - grid("off") + text(start,ma_x,title) + grid('off') f.tight_layout() - savefig(dir + title + ".pdf") + savefig(dir+title+'.pdf') close() + diff --git a/physiological_data/lib/novainstrumentation/filter.py b/physiological_data/lib/novainstrumentation/filter.py index eddda81f..1e717aa6 100644 --- a/physiological_data/lib/novainstrumentation/filter.py +++ b/physiological_data/lib/novainstrumentation/filter.py @@ -5,7 +5,7 @@ def lowpass(s, f, order=2, fs=1000.0, use_filtfilt=False): - """ + ''' @brief: for a given signal s rejects (attenuates) the frequencies higher then the cuttof frequency f and passes the frequencies lower than that value by applying a Butterworth digital filter @@ -29,8 +29,8 @@ def lowpass(s, f, order=2, fs=1000.0, use_filtfilt=False): signal: array-like filtered signal - """ - b, a = signal.butter(order, f / (fs / 2)) + ''' + b, a = signal.butter(order, f / (fs/2)) if use_filtfilt: return filtfilt(b, a, s) @@ -39,7 +39,7 @@ def lowpass(s, f, order=2, fs=1000.0, use_filtfilt=False): def highpass(s, f, order=2, fs=1000.0, use_filtfilt=False): - """ + ''' @brief: for a given signal s rejects (attenuates) the frequencies lower then the cuttof frequency f and passes the frequencies higher than that value by applying a Butterworth digital filter @@ -63,9 +63,9 @@ def highpass(s, f, order=2, fs=1000.0, use_filtfilt=False): signal: array-like filtered signal - """ + ''' - b, a = signal.butter(order, f * 2 / (fs / 2), btype="highpass") + b, a = signal.butter(order, f * 2 / (fs/2), btype='highpass') if use_filtfilt: return filtfilt(b, a, s) @@ -73,7 +73,7 @@ def highpass(s, f, order=2, fs=1000.0, use_filtfilt=False): def bandstop(s, f1, f2, order=2, fs=1000.0, use_filtfilt=False): - """ + ''' @brief: for a given signal s rejects (attenuates) the frequencies within a certain range (between f1 and f2) and passes the frequencies outside that range by applying a Butterworth digital filter @@ -100,15 +100,15 @@ def bandstop(s, f1, f2, order=2, fs=1000.0, use_filtfilt=False): signal: array-like filtered signal - """ - b, a = signal.butter(order, [f1 * 2 / fs, f2 * 2 / fs], btype="bandstop") + ''' + b, a = signal.butter(order, [f1 * 2 / fs, f2 * 2 / fs], btype='bandstop') if use_filtfilt: return filtfilt(b, a, s) return signal.lfilter(b, a, s) def bandpass(s, f1, f2, order=2, fs=1000.0, use_filtfilt=False): - """ + ''' @brief: for a given signal s passes the frequencies within a certain range (between f1 and f2) and rejects (attenuates) the frequencies outside that range by applying a Butterworth digital filter @@ -135,8 +135,8 @@ def bandpass(s, f1, f2, order=2, fs=1000.0, use_filtfilt=False): signal: array-like filtered signal - """ - b, a = signal.butter(order, [f1 * 2 / fs, f2 * 2 / fs], btype="bandpass") + ''' + b, a = signal.butter(order, [f1 * 2 / fs, f2 * 2 / fs], btype='bandpass') if use_filtfilt: return filtfilt(b, a, s) diff --git a/physiological_data/lib/novainstrumentation/freq_analysis.py b/physiological_data/lib/novainstrumentation/freq_analysis.py index 0b499350..1a316584 100644 --- a/physiological_data/lib/novainstrumentation/freq_analysis.py +++ b/physiological_data/lib/novainstrumentation/freq_analysis.py @@ -7,9 +7,9 @@ from novainstrumentation.tools import plotfft -def fundamental_frequency(s, FS): - # TODO: review fundamental frequency to guarantee that f0 exists - # suggestion peak level should be bigger +def fundamental_frequency(s,FS): + # TODO: review fundamental frequency to guarantee that f0 exists + # suggestion peak level should be bigger # TODO: explain code """Compute fundamental frequency along the specified axes. @@ -18,37 +18,37 @@ def fundamental_frequency(s, FS): s: ndarray input from which fundamental frequency is computed. FS: int - sampling frequency + sampling frequency Returns ------- f0: int its integer multiple best explain the content of the signal spectrum. """ - + s = s - mean(s) f, fs = plotfft(s, FS, doplot=False) - - # fs = smooth(fs, 50.0) - - fs = fs[1 : len(fs) / 2] - f = f[1 : len(f) / 2] - + + #fs = smooth(fs, 50.0) + + fs = fs[1:len(fs) / 2] + f = f[1:len(f) / 2] + cond = mlab.find(f > 0.5)[0] - + bp = bigPeaks(fs[cond:], 0) - - if bp == []: - f0 = 0 + + if bp==[]: + f0=0 else: - + bp = bp + cond - + f0 = f[min(bp)] - - return f0 + + return f0 -def max_frequency(sig, FS): +def max_frequency (sig,FS): """Compute max frequency along the specified axes. Parameters @@ -56,22 +56,22 @@ def max_frequency(sig, FS): sig: ndarray input from which max frequency is computed. FS: int - sampling frequency + sampling frequency Returns ------- f_max: int 0.95 of max_frequency using cumsum. """ - - f, fs = plotfft(sig, FS, doplot=False) + + f, fs = plotfft(sig, FS, doplot=False) t = cumsum(fs) - - ind_mag = mlab.find(t > t[-1] * 0.95)[0] - f_max = f[ind_mag] + + ind_mag = mlab.find (t>t[-1]*0.95)[0] + f_max=f[ind_mag] return f_max -def median_frequency(sig, FS): +def median_frequency(sig,FS): """Compute median frequency along the specified axes. Parameters @@ -79,16 +79,17 @@ def median_frequency(sig, FS): sig: ndarray input from which median frequency is computed. FS: int - sampling frequency + sampling frequency Returns ------- f_max: int 0.50 of max_frequency using cumsum. """ - - f, fs = plotfft(sig, FS, doplot=False) + + f, fs = plotfft(sig, FS, doplot=False) t = cumsum(fs) - - ind_mag = mlab.find(t > t[-1] * 0.50)[0] - f_median = f[ind_mag] + + ind_mag = mlab.find (t>t[-1]*0.50)[0] + f_median=f[ind_mag] return f_median + diff --git a/physiological_data/lib/novainstrumentation/niplot.py b/physiological_data/lib/novainstrumentation/niplot.py index e81d16db..dc47a90b 100644 --- a/physiological_data/lib/novainstrumentation/niplot.py +++ b/physiological_data/lib/novainstrumentation/niplot.py @@ -9,15 +9,15 @@ def zoom(event): # edit the scale if needed base_scale = 1.1 - xdata = event.xdata # get event x location - ydata = event.ydata # get event y location + xdata = event.xdata # get event x location + ydata = event.ydata # get event y location # performs a prior check in order to not exceed figure margins if xdata != None and ydata != None: - if event.button == "up": + if event.button == 'up': # deal with zoom in scale_factor = 1 / base_scale - elif event.button == "down": + elif event.button == 'down': # deal with zoom out scale_factor = base_scale else: @@ -28,11 +28,11 @@ def zoom(event): new_width = (cur_xlim[1] - cur_xlim[0]) * scale_factor new_height = (cur_ylim[1] - cur_ylim[0]) * scale_factor - relx = (cur_xlim[1] - xdata) / (cur_xlim[1] - cur_xlim[0]) - rely = (cur_ylim[1] - ydata) / (cur_ylim[1] - cur_ylim[0]) + relx = (cur_xlim[1] - xdata)/(cur_xlim[1] - cur_xlim[0]) + rely = (cur_ylim[1] - ydata)/(cur_ylim[1] - cur_ylim[0]) - ax.set_xlim([xdata - new_width * (1 - relx), xdata + new_width * (relx)]) - ax.set_ylim([ydata - new_height * (1 - rely), ydata + new_height * (rely)]) + ax.set_xlim([xdata - new_width * (1-relx), xdata + new_width * (relx)]) + ax.set_ylim([ydata - new_height * (1-rely), ydata + new_height * (rely)]) ax.figure.canvas.draw() return zoom @@ -41,49 +41,49 @@ def zoom(event): def on_key_press(event): # keyboard zoom-in - if event.key == "+": + if event.key == '+': a = axis() w = a[1] - a[0] - axis([a[0] + w * 0.2, a[1] - w * 0.2, a[2], a[3]]) + axis([a[0] + w * .2, a[1] - w * .2, a[2], a[3]]) draw() # keyboard zoom-out - if event.key in ["-", "'"]: + if event.key in ['-', '\'']: a = axis() w = a[1] - a[0] axis([a[0] - w / 3.0, a[1] + w / 3.0, a[2], a[3]]) draw() # right displacement - if event.key in [".", "right"]: + if event.key in ['.', 'right']: a = axis() w = a[1] - a[0] - axis([a[0] + w * 0.2, a[1] + w * 0.2, a[2], a[3]]) + axis([a[0] + w * .2, a[1] + w * .2, a[2], a[3]]) draw() # left displacement - if event.key in [",", "left"]: + if event.key in [',', 'left']: a = axis() w = a[1] - a[0] - axis([a[0] - w * 0.2, a[1] - w * 0.2, a[2], a[3]]) + axis([a[0] - w * .2, a[1] - w * .2, a[2], a[3]]) draw() # up displacement - if event.key == "up": + if event.key == 'up': a = axis() w = a[3] - a[2] - axis([a[0], a[1], a[2] + w * 0.2, a[3] + w * 0.2]) + axis([a[0], a[1], a[2] + w * .2, a[3] + w * .2]) draw() # down displacement - if event.key == "down": + if event.key == 'down': a = axis() w = a[3] - a[2] - axis([a[0], a[1], a[2] - w * 0.2, a[3] - w * 0.2]) + axis([a[0], a[1], a[2] - w * .2, a[3] - w * .2]) draw() # close figure - if event.key == "q": + if event.key == 'q': close() # NOTE: We should make the disconnect (mpl_disconect(cid) # But since the figure is destroyed we may keep this format @@ -110,6 +110,9 @@ def niplot(): """ fig = gcf() - cid = fig.canvas.mpl_connect("key_press_event", on_key_press) # @UnusedVariable - cid = fig.canvas.mpl_connect("key_release_event", on_key_release) # @UnusedVariable - cid = fig.canvas.mpl_connect("scroll_event", zoom) + cid = fig.canvas.mpl_connect('key_press_event', # @UnusedVariable + on_key_press) + cid = fig.canvas.mpl_connect('key_release_event', # @UnusedVariable + on_key_release) + cid = fig.canvas.mpl_connect('scroll_event', zoom) + diff --git a/physiological_data/lib/novainstrumentation/panthomkins/detect_panthomkins_peaks.py b/physiological_data/lib/novainstrumentation/panthomkins/detect_panthomkins_peaks.py index 57b555bf..65a35add 100644 --- a/physiological_data/lib/novainstrumentation/panthomkins/detect_panthomkins_peaks.py +++ b/physiological_data/lib/novainstrumentation/panthomkins/detect_panthomkins_peaks.py @@ -9,17 +9,9 @@ __license__ = "MIT" -def detect_panthomkins_peaks( - x, - mph=None, - mpd=1, - threshold=0, - edge="rising", - kpsh=False, - valley=False, - show=False, - ax=None, -): +def detect_panthomkins_peaks(x, mph=None, mpd=1, threshold=0, edge='rising', + kpsh=False, valley=False, show=False, ax=None): + """Detect peaks in data based on their amplitude and other features. Parameters @@ -62,10 +54,10 @@ def detect_panthomkins_peaks( References ---------- - .. [1] http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/DetectPeaks.ipynb - """ + .. [1] http://nbviewer.ipython.org/github/demotu/BMC/blob/master/notebooks/DetectPeaks.ipynb""" + - x = np.atleast_1d(x).astype("float64") + x = np.atleast_1d(x).astype('float64') if x.size < 3: return np.array([], dtype=int) if valley: @@ -81,30 +73,26 @@ def detect_panthomkins_peaks( if not edge: ine = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) > 0))[0] else: - if edge.lower() in ["rising", "both"]: + if edge.lower() in ['rising', 'both']: ire = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0] - if edge.lower() in ["falling", "both"]: + if edge.lower() in ['falling', 'both']: ife = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) >= 0))[0] ind = np.unique(np.hstack((ine, ire, ife))) # handle NaN's if ind.size and indnan.size: # NaN's and values close to NaN's cannot be peaks - ind = ind[ - np.in1d( - ind, np.unique(np.hstack((indnan, indnan - 1, indnan + 1))), invert=True - ) - ] + ind = ind[np.in1d(ind, np.unique(np.hstack((indnan, indnan-1, indnan+1))), invert=True)] # first and last values of x cannot be peaks if ind.size and ind[0] == 0: ind = ind[1:] - if ind.size and ind[-1] == x.size - 1: + if ind.size and ind[-1] == x.size-1: ind = ind[:-1] # remove peaks < minimum peak height if ind.size and mph is not None: ind = ind[x[ind] >= mph] # remove peaks - neighbors < threshold if ind.size and threshold > 0: - dx = np.min(np.vstack([x[ind] - x[ind - 1], x[ind] - x[ind + 1]]), axis=0) + dx = np.min(np.vstack([x[ind]-x[ind-1], x[ind]-x[ind+1]]), axis=0) ind = np.delete(ind, np.where(dx < threshold)[0]) # detect small peaks closer than minimum peak distance if ind.size and mpd > 1: @@ -113,9 +101,8 @@ def detect_panthomkins_peaks( for i in range(ind.size): if not idel[i]: # keep peaks with the same height if kpsh is True - idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) & ( - x[ind[i]] > x[ind] if kpsh else True - ) + idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) \ + & (x[ind[i]] > x[ind] if kpsh else True) idel[i] = 0 # Keep current peak # remove the small peaks and sort back the indices by their occurrence ind = np.sort(ind[~idel]) @@ -135,34 +122,24 @@ def _plot(x, mph, mpd, threshold, edge, valley, ax, ind): try: import matplotlib.pyplot as plt except ImportError: - print("matplotlib is not available.") + print('matplotlib is not available.') else: if ax is None: _, ax = plt.subplots(1, 1, figsize=(8, 4)) - ax.plot(x, "b", lw=1) + ax.plot(x, 'b', lw=1) if ind.size: - label = "valley" if valley else "peak" - label = label + "s" if ind.size > 1 else label - ax.plot( - ind, - x[ind], - "+", - mfc=None, - mec="r", - mew=2, - ms=8, - label="%d %s" % (ind.size, label), - ) - ax.legend(loc="best", framealpha=0.5, numpoints=1) - ax.set_xlim(-0.02 * x.size, x.size * 1.02 - 1) + label = 'valley' if valley else 'peak' + label = label + 's' if ind.size > 1 else label + ax.plot(ind, x[ind], '+', mfc=None, mec='r', mew=2, ms=8, + label='%d %s' % (ind.size, label)) + ax.legend(loc='best', framealpha=.5, numpoints=1) + ax.set_xlim(-.02*x.size, x.size*1.02-1) ymin, ymax = x[np.isfinite(x)].min(), x[np.isfinite(x)].max() yrange = ymax - ymin if ymax > ymin else 1 - ax.set_ylim(ymin - 0.1 * yrange, ymax + 0.1 * yrange) - ax.set_xlabel("Data #", fontsize=14) - ax.set_ylabel("Amplitude", fontsize=14) - mode = "Valley detection" if valley else "Peak detection" - ax.set_title( - "%s (mph=%s, mpd=%d, threshold=%s, edge='%s')" - % (mode, str(mph), mpd, str(threshold), edge) - ) + ax.set_ylim(ymin - 0.1*yrange, ymax + 0.1*yrange) + ax.set_xlabel('Data #', fontsize=14) + ax.set_ylabel('Amplitude', fontsize=14) + mode = 'Valley detection' if valley else 'Peak detection' + ax.set_title("%s (mph=%s, mpd=%d, threshold=%s, edge='%s')" + % (mode, str(mph), mpd, str(threshold), edge)) diff --git a/physiological_data/lib/novainstrumentation/panthomkins/panthomkins.py b/physiological_data/lib/novainstrumentation/panthomkins/panthomkins.py index 9f195ec3..35d7dfc5 100644 --- a/physiological_data/lib/novainstrumentation/panthomkins/panthomkins.py +++ b/physiological_data/lib/novainstrumentation/panthomkins/panthomkins.py @@ -3,9 +3,7 @@ import numpy as np from novainstrumentation.panthomkins.butterworth_filters import butter_bandpass_filter -from novainstrumentation.panthomkins.detect_panthomkins_peaks import ( - detect_panthomkins_peaks, -) +from novainstrumentation.panthomkins.detect_panthomkins_peaks import detect_panthomkins_peaks from novainstrumentation.panthomkins.rr_update import rr_1_update, rr_2_update, sync @@ -19,7 +17,7 @@ def panthomkins(ecg_signal, fs, butterlow=8, butterhigh=30): ecg_filter = butter_bandpass_filter(ecg, butterlow, butterhigh, fs) # Squaring - ecg_filter = 50.0 * ecg_filter**2.0 + ecg_filter = 50.0 * ecg_filter ** 2.0 # Find Peaks pksInd = detect_panthomkins_peaks(ecg_filter, mpd=35) @@ -91,18 +89,17 @@ def panthomkins(ecg_signal, fs, butterlow=8, butterhigh=30): if NFound_Old != NFound - 1: rr_1, rr_average_1 = rr_1_update(rr_1, NFound - 1, Found) rr_2, rr_average_2, flag, rr_low_limit, rr_high_limit = rr_2_update( - rr_2, NFound - 1, Found, rr_low_limit, rr_high_limit - ) + rr_2, NFound - 1, Found, rr_low_limit, rr_high_limit) NFound_Old = NFound - 1 - # if np.mod(NFound, 8) == 0: - # print(['Average of the 8 most recent HR is ', - # str(rr_average_1 / fs * 60.0), ' (BPM)']) - # print('') + #if np.mod(NFound, 8) == 0: + #print(['Average of the 8 most recent HR is ', + #str(rr_average_1 / fs * 60.0), ' (BPM)']) + #print('') if flag: - print("Gap Found") + print('Gap Found') flag = 0 back = ii diff --git a/physiological_data/lib/novainstrumentation/panthomkins/rr_update.py b/physiological_data/lib/novainstrumentation/panthomkins/rr_update.py index 0620dcfb..c05aa0e6 100644 --- a/physiological_data/lib/novainstrumentation/panthomkins/rr_update.py +++ b/physiological_data/lib/novainstrumentation/panthomkins/rr_update.py @@ -3,9 +3,9 @@ def rr_1_update(rr_1, NFound, Found): if np.logical_and(NFound <= 7, NFound > 0): - rr_1[0 : NFound - 1] = np.ediff1d(Found[0:NFound, 0]) + rr_1[0:NFound - 1] = np.ediff1d(Found[0:NFound, 0]) elif NFound > 7: - rr_1 = np.ediff1d((Found[NFound - 7 : NFound - 1, 0])) + rr_1 = np.ediff1d((Found[NFound - 7:NFound - 1, 0])) rr_average_1 = np.mean(rr_1) @@ -14,10 +14,10 @@ def rr_1_update(rr_1, NFound, Found): def rr_2_update(rr_2, NFound, Found, rr_low_limit, rr_high_limit): if NFound > 0: - delta = np.ediff1d(Found[NFound - 1 : NFound, 0]) + delta = np.ediff1d(Found[NFound - 1:NFound, 0]) if np.logical_and(delta >= rr_low_limit, delta <= rr_high_limit): pos = np.mod(NFound, 7) - if pos == 0: + if (pos == 0): rr_2[7] = delta else: rr_2[pos] = delta @@ -42,12 +42,12 @@ def sync(Found, NFound, ecg, N): for ii in range(0, NFound): xtemp = Found[ii, 0] - if xtemp - 60 > 0: + if (xtemp - 60 > 0): indInf = xtemp - 60 else: indInf = 0 - if xtemp + 60 < N: + if (xtemp + 60 < N): indSup = xtemp + 60 else: indSup = N diff --git a/physiological_data/lib/novainstrumentation/peakdelta.py b/physiological_data/lib/novainstrumentation/peakdelta.py index b24c20a2..60249d2e 100644 --- a/physiological_data/lib/novainstrumentation/peakdelta.py +++ b/physiological_data/lib/novainstrumentation/peakdelta.py @@ -14,7 +14,6 @@ ########################### Peaks Detection ################################## ############################################################################## - def peakdelta(v, delta, x=None): """ Returns two arrays @@ -25,7 +24,7 @@ def peakdelta(v, delta, x=None): % maxima and minima ("peaks") in the vector V. % MAXTAB and MINTAB consists of two columns. Column 1 % contains indices in V, and column 2 the found values. - % + % % With [MAXTAB, MINTAB] = peakdelta(V, DELTA, X) the indices % in MAXTAB and MINTAB are replaced with the corresponding % X-values. @@ -38,21 +37,21 @@ def peakdelta(v, delta, x=None): % This function is released to the public domain; Any use is allowed. """ maxtab = [] - mintab = [] - + mintab = [] + if x is None: x = arange(len(v)) - + v = asarray(v) - + if len(v) != len(x): - sys.exit("Input vectors v and x must have same length") + sys.exit('Input vectors v and x must have same length') if not isscalar(delta): - sys.exit("Input argument delta must be a scalar") + sys.exit('Input argument delta must be a scalar') if delta <= 0: - sys.exit("Input argument delta must be positive") + sys.exit('Input argument delta must be positive') mn, mx = Inf, -Inf mnpos, mxpos = NaN, NaN @@ -82,3 +81,4 @@ def peakdelta(v, delta, x=None): lookformax = True return array(maxtab), array(mintab) + diff --git a/physiological_data/lib/novainstrumentation/peaks.py b/physiological_data/lib/novainstrumentation/peaks.py index c4d16615..c9a3546f 100644 --- a/physiological_data/lib/novainstrumentation/peaks.py +++ b/physiological_data/lib/novainstrumentation/peaks.py @@ -9,13 +9,15 @@ from numpy import array, clip, argsort, sort -# from scipy.signal import argrelmax +#from scipy.signal import argrelmax #### FROM SCIPY #### Given that the new version failed to compile in some linux versions -def _boolrelextrema(data, comparator, axis=0, order=1, mode="clip"): + +def _boolrelextrema(data, comparator, + axis=0, order=1, mode='clip'): """ Calculate the relative extrema of `data`. @@ -54,8 +56,8 @@ def _boolrelextrema(data, comparator, axis=0, order=1, mode="clip"): array([False, False, True, False, False], dtype=bool) """ - if (int(order) != order) or (order < 1): - raise ValueError("Order must be an int >= 1") + if((int(order) != order) or (order < 1)): + raise ValueError('Order must be an int >= 1') datalen = data.shape[axis] locs = np.arange(0, datalen) @@ -67,12 +69,12 @@ def _boolrelextrema(data, comparator, axis=0, order=1, mode="clip"): minus = data.take(locs - shift, axis=axis, mode=mode) results &= comparator(main, plus) results &= comparator(main, minus) - if ~results.any(): + if(~results.any()): return results return results -def argrelmin(data, axis=0, order=1, mode="clip"): +def argrelmin(data, axis=0, order=1, mode='clip'): """ Calculate the relative minima of `data`. @@ -110,7 +112,7 @@ def argrelmin(data, axis=0, order=1, mode="clip"): return argrelextrema(data, np.less, axis, order, mode) -def argrelmax(data, axis=0, order=1, mode="clip"): +def argrelmax(data, axis=0, order=1, mode='clip'): """ Calculate the relative maxima of `data`. @@ -148,7 +150,7 @@ def argrelmax(data, axis=0, order=1, mode="clip"): return argrelextrema(data, np.greater, axis, order, mode) -def argrelextrema(data, comparator, axis=0, order=1, mode="clip"): +def argrelextrema(data, comparator, axis=0, order=1, mode='clip'): """ Calculate the relative extrema of `data`. @@ -182,18 +184,20 @@ def argrelextrema(data, comparator, axis=0, order=1, mode="clip"): argrelmin, argrelmax """ - results = _boolrelextrema(data, comparator, axis, order, mode) + results = _boolrelextrema(data, comparator, + axis, order, mode) if ~results.any(): return (np.array([]),) * 2 else: return np.where(results) -#### FROM SCIPY +#### FROM SCIPY + def peaks(signal, tol=None): - """This function detects all the peaks of a signal and returns those time + """ This function detects all the peaks of a signal and returns those time positions. To reduce the amount of peaks detected, a threshold is introduced so only the peaks above that value are considered. @@ -221,24 +225,24 @@ def peaks(signal, tol=None): array([], dtype=int32) """ - if tol is None: + if (tol is None): tol = min(signal) pks = argrelmax(clip(signal, tol, signal.max())) return pks[0] def post_peak(v, v_post): - """Detects the next peak""" + """ Detects the next peak """ return array([v_post[find(v_post > i)[0]] for i in v]) def prior_peak(v, v_prior): - """Detects the previous peak""" + """ Detects the previous peak """ return array([v_prior[find(v_prior < i)[-1]] for i in v]) def clean_near_peaks(signal, peaks_, min_distance): - """Given an array with all the peaks of the signal ('peaks') and a + """ Given an array with all the peaks of the signal ('peaks') and a distance value ('min_distance') and the signal, by argument, this function erases all the unnecessary peaks and returns an array with only the maximum peak for each period of the signal (the period is given by the @@ -261,15 +265,15 @@ def clean_near_peaks(signal, peaks_, min_distance): See also: clean_near_events() """ - # order all peaks - + #order all peaks + ars = argsort(signal[peaks_]) - + pp = peaks(ars) fp = [] - # clean near peaks + #clean near peaks while len(pp) > 0: fp += [pp[-1]] pp = pp[abs(pp - pp[-1]) > min_distance] @@ -278,7 +282,7 @@ def clean_near_peaks(signal, peaks_, min_distance): def clean_near_events(points, min_distance): - """Given an array with some specific points of the signal and a distance + """ Given an array with some specific points of the signal and a distance value, this function erases all the surplus points and returns an array with only one point (the first point) per distance samples values @@ -308,26 +312,24 @@ def clean_near_events(points, min_distance): points = points[abs(points - points[0]) > min_distance] return array(fp) - - + from numpy import ceil - -def bigPeaks(s, th, min_peak_distance=5, peak_return_percentage=0.1): - p = peaks(s, th) - pp = [] - if len(p) == 0: - pp = [] +def bigPeaks(s,th,min_peak_distance=5,peak_return_percentage=0.1): + p=peaks(s,th) + pp=[] + if len(p)==0: + pp=[] else: - p = clean_near_peaks(s, p, min_peak_distance) - - if len(p) != 0: - ars = argsort(s[p]) - pp = p[ars] - - num_peaks_to_return = ceil(len(p) * peak_return_percentage) - - pp = pp[-num_peaks_to_return:] + p=clean_near_peaks(s,p,min_peak_distance) + + if len(p)!=0: + ars=argsort(s[p]) + pp=p[ars] + + num_peaks_to_return=ceil(len(p)*peak_return_percentage) + + pp=pp[-num_peaks_to_return:] else: - pp == [] + pp==[] return pp diff --git a/physiological_data/lib/novainstrumentation/sandbox/TestClustering.py b/physiological_data/lib/novainstrumentation/sandbox/TestClustering.py index bf145eb7..3dcc0ab5 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/TestClustering.py +++ b/physiological_data/lib/novainstrumentation/sandbox/TestClustering.py @@ -1,260 +1,201 @@ import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D -from sklearn.cluster import ( - MiniBatchKMeans, - KMeans, - SpectralClustering, - DBSCAN, - AgglomerativeClustering, - AffinityPropagation, -) +from sklearn.cluster import MiniBatchKMeans, KMeans, SpectralClustering, DBSCAN, AgglomerativeClustering, \ + AffinityPropagation from sklearn.neighbors import kneighbors_graph from sklearn.preprocessing import StandardScaler -def MultiDimensionalClusteringKmeans( - Xmatrix, time, xdata, n_clusters=2, ax=None, show=False -): - - seed = np.random.seed(0) - colors = np.array([x for x in "bgrcmykbgrcmykbgrcmykbgrcmyk"]) - colors = np.hstack([colors] * 20) +def MultiDimensionalClusteringKmeans(Xmatrix, time, xdata, n_clusters=2, ax = None, show=False): - # normalize dataset for easier parameter selection - X = StandardScaler().fit_transform(Xmatrix) + seed = np.random.seed(0) + colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk']) + colors = np.hstack([colors] * 20) - # algorithm KMeans - kmeans = KMeans(n_clusters=n_clusters, random_state=seed) - - # Apply algorithm - kmeans.fit(X) + #normalize dataset for easier parameter selection + X = StandardScaler().fit_transform(Xmatrix) - y_pred = kmeans.labels_.astype(np.int) - centers = kmeans.cluster_centers_ - center_colors = colors[: len(centers)] + #algorithm KMeans + kmeans = KMeans(n_clusters=n_clusters, random_state=seed) - # Representation + #Apply algorithm + kmeans.fit(X) - if np.logical_and(show, ax != None): + y_pred = kmeans.labels_.astype(np.int) + centers = kmeans.cluster_centers_ + center_colors = colors[:len(centers)] - ax.set_title( - "Clustering Tech: " - + "KMEANS; " - + "Number of Clusters = " - + str(n_clusters), - fontsize=15, - ) + #Representation - ax.plot(time, xdata, color="lightgray", alpha=0.4) - ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - ax.set_xlabel("time (ms)") - ax.set_ylabel("Amplitude") - - return X, y_pred, centers, center_colors - - # fig = plt.figure(5, figsize=(4, 3)) - # plt.clf() - # ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134) - # plt.cla() - # - # ax.scatter(X[:, 2], X[:, 0], X[:, 1], color=colors[y_pred].tolist()) - - elif np.logical_and(show, ax == None): - fig, axis = plt.subplots(2, 1) - fig.tight_layout() - axis[0].set_title( - "Clustering Tech: " - + "KMEANS; " - + "Number of Clusters = " - + str(n_clusters), - fontsize=15, - ) - axis[0].plot(time, xdata, color="lightgray", alpha=0.4) - axis[0].scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - axis[0].set_xlabel("time (ms)") - axis[0].set_ylabel("Amplitude") - axis[1].plot(y_pred) - - return X, y_pred, centers, center_colors - - else: - - return X, y_pred, centers, center_colors - - -def MultiDimensionalClusteringSPCL( - Xmatrix, time, xdata, eigen_solver="arpack", n_clusters=2, ax=None, show=False -): - - seed = np.random.seed(0) - colors = np.array([x for x in "bgrcmykbgrcmykbgrcmykbgrcmyk"]) - colors = np.hstack([colors] * 20) - - # normalize dataset for easier parameter selection - X = StandardScaler().fit_transform(Xmatrix) - # algorithm SpectralClustering - SC = SpectralClustering( - n_clusters=n_clusters, eigen_solver=eigen_solver, affinity="nearest_neighbors" - ) - - # Apply algorithm - fit = SC.fit(X) - - y_pred = fit.labels_.astype(np.int) - - # Representation - if np.logical_and(show, ax == None): - - ax.set_title( - "Clustering Tech: " - + "SpectralClustering; " - + "Number of Clusters = " - + str(n_clusters), - fontsize=15, - ) - ax.plot(time, xdata, color="lightgray", alpha=0.4) - ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - ax.set_xlabel("time (ms)") - ax.set_ylabel("Amplitude") - - return X, y_pred - - elif np.logical_and(show, ax == None): - - fig, axis = plt.subplots(1, 1) - fig.tight_layout() - axis.plot(time, xdata, color="lightgray", alpha=0.4) - axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - axis.set_xlabel("time (ms)") - axis.set_ylabel("Amplitude") - - return X, y_pred - - else: - return X, y_pred - - -def MultiDimensionalClusteringAGG( - Xmatrix, - time, - xdata, - n_clusters=2, - Linkage="ward", - Affinity="euclidean", - ax=None, - show=False, -): - - # normalize dataset for easier parameter selection - X = StandardScaler().fit_transform(Xmatrix) - - # connectivity matrix for structured Ward - connectivity = kneighbors_graph(X, n_neighbors=10, include_self=False) - # make connectivity symmetric - connectivity = 0.5 * (connectivity + connectivity.T) - - seed = np.random.seed(0) - colors = np.array([x for x in "bgrcmykbgrcmykbgrcmykbgrcmyk"]) - colors = np.hstack([colors] * 20) - - if Linkage is "ward": - ward = AgglomerativeClustering( - n_clusters=n_clusters, linkage="ward", connectivity=connectivity - ) - # Apply algorithm - fit = ward.fit(X) - else: - Agg = AgglomerativeClustering( - n_clusters=n_clusters, - linkage=Linkage, - connectivity=connectivity, - affinity=Affinity, - ) - # Apply algorithm - fit = Agg.fit(X) - - y_pred = fit.labels_.astype(np.int) - - # Representation - if np.logical_and(show, ax != None): - - ax.set_title( - "Clustering Tech: " - + "Agglomerative Clustering " - + "Number of Clusters = " - + str(n_clusters), - fontsize=15, - ) - ax.plot(time, xdata, color="lightgray", alpha=0.4) - ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - ax.set_xlabel("time (ms)") - ax.set_ylabel("Amplitude") - - return X, y_pred - - elif np.logical_and(show, ax == None): - - fig, axis = plt.subplots(1, 1) - fig.tight_layout() - axis.set_title( - "Clustering Tech: " - + "Agglomerative Clustering " - + "Number of Clusters = " - + str(n_clusters), - fontsize=15, - ) - axis.plot(time, xdata, color="lightgray", alpha=0.4) - axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - axis.set_xlabel("time (ms)") - axis.set_ylabel("Amplitude") - - else: - - return X, y_pred + if np.logical_and(show, ax != None): + ax.set_title('Clustering Tech: ' + "KMEANS; " + 'Number of Clusters = ' + str(n_clusters), fontsize=15) + + ax.plot(time, xdata, color='lightgray', alpha=0.4) + ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + ax.set_xlabel("time (ms)") + ax.set_ylabel("Amplitude") + + return X, y_pred, centers, center_colors + + # fig = plt.figure(5, figsize=(4, 3)) + # plt.clf() + # ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134) + # plt.cla() + # + # ax.scatter(X[:, 2], X[:, 0], X[:, 1], color=colors[y_pred].tolist()) + + elif np.logical_and(show, ax == None): + fig, axis = plt.subplots(2,1) + fig.tight_layout() + axis[0].set_title('Clustering Tech: ' + "KMEANS; " + 'Number of Clusters = ' + str(n_clusters), fontsize=15) + axis[0].plot(time, xdata, color='lightgray', alpha=0.4) + axis[0].scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + axis[0].set_xlabel("time (ms)") + axis[0].set_ylabel("Amplitude") + axis[1].plot(y_pred) + + return X, y_pred, centers, center_colors + + else: + + return X, y_pred, centers, center_colors + + +def MultiDimensionalClusteringSPCL(Xmatrix, time, xdata, eigen_solver = 'arpack', n_clusters=2, ax = None, show=False): + + seed = np.random.seed(0) + colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk']) + colors = np.hstack([colors] * 20) + + # normalize dataset for easier parameter selection + X = StandardScaler().fit_transform(Xmatrix) + # algorithm SpectralClustering + SC = SpectralClustering(n_clusters=n_clusters, eigen_solver=eigen_solver, affinity="nearest_neighbors") + + # Apply algorithm + fit = SC.fit(X) + + y_pred = fit.labels_.astype(np.int) + + # Representation + if np.logical_and(show, ax == None): + + ax.set_title('Clustering Tech: ' + "SpectralClustering; " + 'Number of Clusters = ' + str(n_clusters), fontsize=15) + ax.plot(time, xdata, color='lightgray', alpha=0.4) + ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + ax.set_xlabel("time (ms)") + ax.set_ylabel("Amplitude") + + return X, y_pred + + elif np.logical_and(show, ax == None): + + fig, axis = plt.subplots(1, 1) + fig.tight_layout() + axis.plot(time, xdata, color='lightgray', alpha=0.4) + axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + axis.set_xlabel("time (ms)") + axis.set_ylabel("Amplitude") + + return X, y_pred + + else: + return X, y_pred + +def MultiDimensionalClusteringAGG(Xmatrix, time, xdata, n_clusters=2, Linkage = 'ward', Affinity = 'euclidean', ax = None, show=False): + + # normalize dataset for easier parameter selection + X = StandardScaler().fit_transform(Xmatrix) + + # connectivity matrix for structured Ward + connectivity = kneighbors_graph(X, n_neighbors=10, include_self=False) + # make connectivity symmetric + connectivity = 0.5 * (connectivity + connectivity.T) + + seed = np.random.seed(0) + colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk']) + colors = np.hstack([colors] * 20) + + + if(Linkage is 'ward'): + ward = AgglomerativeClustering(n_clusters=n_clusters, linkage='ward', connectivity=connectivity) + # Apply algorithm + fit = ward.fit(X) + else: + Agg = AgglomerativeClustering(n_clusters=n_clusters, linkage=Linkage, connectivity=connectivity, affinity=Affinity) + # Apply algorithm + fit = Agg.fit(X) + + y_pred = fit.labels_.astype(np.int) + + # Representation + if np.logical_and(show, ax != None): + + ax.set_title('Clustering Tech: ' + "Agglomerative Clustering " + 'Number of Clusters = ' + str(n_clusters), fontsize=15) + ax.plot(time, xdata, color='lightgray', alpha=0.4) + ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + ax.set_xlabel("time (ms)") + ax.set_ylabel("Amplitude") + + return X, y_pred + + elif np.logical_and(show, ax == None): + + fig, axis = plt.subplots(1, 1) + fig.tight_layout() + axis.set_title('Clustering Tech: ' + "Agglomerative Clustering " + 'Number of Clusters = ' + str(n_clusters), fontsize=15) + axis.plot(time, xdata, color='lightgray', alpha=0.4) + axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + axis.set_xlabel("time (ms)") + axis.set_ylabel("Amplitude") + + else: + + return X, y_pred def MultiDimensionalClusteringDBSCAN(Xmatrix, time, xdata, eps, ax=None, show=False): - # normalize dataset for easier parameter selection - X = StandardScaler().fit_transform(Xmatrix) + # normalize dataset for easier parameter selection + X = StandardScaler().fit_transform(Xmatrix) - colors = np.array([x for x in "bgrcmykbgrcmykbgrcmykbgrcmyk"]) - colors = np.hstack([colors] * 20) + colors = np.array([x for x in 'bgrcmykbgrcmykbgrcmykbgrcmyk']) + colors = np.hstack([colors] * 20) - # normalize dataset for easier parameter selection - dbscan = DBSCAN(eps=eps) + # normalize dataset for easier parameter selection + dbscan = DBSCAN(eps=eps) - # Apply algorithm - fit = dbscan.fit(X) + # Apply algorithm + fit = dbscan.fit(X) - y_pred = fit.labels_.astype(np.int) + y_pred = fit.labels_.astype(np.int) - # Representation - if np.logical_and(show, ax != None): + # Representation + if np.logical_and(show, ax != None): - ax.set_title("Clustering Tech: " + "DBSCAN", fontsize=15) + ax.set_title('Clustering Tech: ' + "DBSCAN", fontsize=15) - ax.plot(time, xdata, color="lightgray", alpha=0.4) - ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - ax.set_xlabel("time (ms)") - ax.set_ylabel("Amplitude") + ax.plot(time, xdata, color='lightgray', alpha=0.4) + ax.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + ax.set_xlabel("time (ms)") + ax.set_ylabel("Amplitude") - return X, y_pred + return X, y_pred - elif np.logical_and(show, ax == None): - fig, axis = plt.subplots(1, 1) - fig.tight_layout() - axis.set_title("Clustering Tech: " + "DBSCAN", fontsize=15) + elif np.logical_and(show, ax == None): + fig, axis = plt.subplots(1, 1) + fig.tight_layout() + axis.set_title('Clustering Tech: ' + "DBSCAN", fontsize=15) - axis.plot(time, xdata, color="lightgray", alpha=0.4) - axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) - axis.set_xlabel("time (ms)") - axis.set_ylabel("Amplitude") + axis.plot(time, xdata, color='lightgray', alpha=0.4) + axis.scatter(time, xdata, color=colors[y_pred].tolist(), s=10) + axis.set_xlabel("time (ms)") + axis.set_ylabel("Amplitude") - else: + else: - return X, y_pred + return X, y_pred # def MultidimensionalClustering(FeaturesMatrix, Xarray, Yarray, n_clusters = 2, ClusteringMethod = 'kmeans', ax= None, show = False): @@ -282,4 +223,4 @@ def MultiDimensionalClusteringDBSCAN(Xmatrix, time, xdata, eps, ax=None, show=Fa # """ # # # normalize dataset for easier parameter selection -# X = StandardScaler().fit_transform(FeaturesMatrix) +# X = StandardScaler().fit_transform(FeaturesMatrix) \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/sandbox/fastperlin.py b/physiological_data/lib/novainstrumentation/sandbox/fastperlin.py index 72524e49..35ef6421 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/fastperlin.py +++ b/physiological_data/lib/novainstrumentation/sandbox/fastperlin.py @@ -1,22 +1,7 @@ -# @PydevCodeAnalysisIgnore - -from numpy import ( - array, - abs, - arange, - dot, - int8, - int32, - floor, - fromfunction, - hypot, - ones, - prod, - random, - indices, - newaxis, - poly1d, -) +#@PydevCodeAnalysisIgnore + +from numpy import array, abs, arange, dot, int8, int32, floor, fromfunction,\ + hypot, ones, prod, random, indices, newaxis, poly1d class PerlinNoise(object): @@ -27,61 +12,59 @@ def noise(self, coords): uvw = coords - ijk - indexes = self.P[ijk[:, :, self.order - 1]] + indexes = self.P[ijk[:,:, self.order - 1]] for i in range(self.order - 1): - indexes = self.P[(ijk[:, :, i] + indexes) % len(self.P)] + indexes = self.P[(ijk[:,:, i] + indexes) % len(self.P)] gradiens = self.G[indexes % len(self.G)] - # gradiens = self.G[(ijk[:,:, 0] + indexes) % len(self.G)] - - res = ( - self.drop(abs(uvw)).prod(axis=2) * prod([gradiens, uvw], axis=0).sum(axis=2) - ).sum(axis=1) +# gradiens = self.G[(ijk[:,:, 0] + indexes) % len(self.G)] + + res = (self.drop(abs(uvw)).prod(axis=2)*prod([gradiens, uvw], axis=0).sum(axis=2)).sum(axis=1) res[res > 1.0] = 1.0 res[res < -1.0] = -1.0 - return ((res + 1) * 128).astype(float) + return ((res + 1)*128).astype(float) def getData(self, scale=32.0): return self.noise(indices(self.size).reshape(self.order, 1, -1).T / scale) + def __init__(self, size=None, n=None): - n = n if n else 256 + n = n if n else 256 self.size = size if size else (256, 256) self.order = len(self.size) - + # Generate WAY more numbers than we need # because we are throwing out all the numbers not inside a unit # sphere. Something of a hack but statistically speaking # it should work fine... or crash. - G = (random.uniform(size=2 * self.order * n) * 2 - 1).reshape(-1, self.order) + G = (random.uniform(size=2*self.order*n)*2 - 1).reshape(-1, self.order) # GAH! How do I generalize this?! - # length = hypot(G[:,i] for i in range(self.order)) + #length = hypot(G[:,i] for i in range(self.order)) if self.order == 1: - length = G[:, 0] + length = G[:,0] elif self.order == 2: - length = hypot(G[:, 0], G[:, 1]) + length = hypot(G[:,0], G[:,1]) elif self.order == 3: - length = hypot(G[:, 0], G[:, 1], G[:, 2]) - - self.G = (G[length < 1] / (length[length < 1])[:, newaxis])[:n,] + length = hypot(G[:,0], G[:,1], G[:,2]) + + self.G = (G[length < 1] / (length[length < 1])[:,newaxis])[:n,] self.P = arange(n, dtype=int32) - + random.shuffle(self.P) - - self.idx_ar = ( - indices(2 * ones(self.order), dtype=int8).reshape(self.order, -1).T - ) + + self.idx_ar = indices(2*ones(self.order), dtype=int8).reshape(self.order, -1).T self.drop = poly1d((-6, 15, -10, 0, 0, 1.0)) - def perlin(n, scale=32.0): pn = PerlinNoise(size=(1, n)) return pn.getData(scale) + + diff --git a/physiological_data/lib/novainstrumentation/sandbox/randomemg.py b/physiological_data/lib/novainstrumentation/sandbox/randomemg.py index ec170a13..49a5bbdd 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/randomemg.py +++ b/physiological_data/lib/novainstrumentation/sandbox/randomemg.py @@ -2,7 +2,7 @@ from pylab import * -def smooth(x, window_len=10, window="hanning"): +def smooth(x,window_len=10,window='hanning'): """smooth the data using a window with requested size. This method is based on the convolution of a scaled window with the signal. @@ -34,43 +34,41 @@ def smooth(x, window_len=10, window="hanning"): """ if x.ndim != 1: - raise (ValueError, "smooth only accepts 1 dimension arrays.") + raise(ValueError, "smooth only accepts 1 dimension arrays.") if x.size < window_len: - raise (ValueError, "Input vector needs to be bigger than window size.") + raise(ValueError, "Input vector needs to be bigger than window size.") - if window_len < 3: + + if window_len<3: return x - if not window in ["flat", "hanning", "hamming", "bartlett", "blackman"]: - raise ( - ValueError, - "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'", - ) - s = numpy.r_[2 * x[0] - x[window_len:1:-1], x, 2 * x[-1] - x[-1:-window_len:-1]] - # print(len(s)) - if window == "flat": # moving average - w = numpy.ones(window_len, "d") - else: - w = eval("numpy." + window + "(window_len)") + if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: + raise(ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'") + - y = numpy.convolve(w / w.sum(), s, mode="same") - return y[window_len - 1 : -window_len + 1] + s=numpy.r_[2*x[0]-x[window_len:1:-1],x,2*x[-1]-x[-1:-window_len:-1]] + #print(len(s)) + if window == 'flat': #moving average + w=numpy.ones(window_len,'d') + else: + w=eval('numpy.'+window+'(window_len)') + y=numpy.convolve(w/w.sum(),s,mode='same') + return y[window_len-1:-window_len+1] def randomEMG(identifier): - + seed(identifier) t = numpy.arange(1000) - on_i = numpy.random.randint(100, 300) - off_i = numpy.random.randint(500, 600) - gain = numpy.random.randint(10, 20) - + on_i = numpy.random.randint(100,300) + off_i = numpy.random.randint(500,600) + gain = numpy.random.randint (10,20) + sinal = numpy.random.randn(1000) - envelope = numpy.concatenate( - (numpy.ones(on_i), gain * numpy.ones(off_i - on_i), numpy.ones(1000 - off_i)) - ) - e = smooth(envelope, 50) + envelope = numpy.concatenate((numpy.ones(on_i),gain*numpy.ones(off_i-on_i),numpy.ones(1000-off_i))) + e = smooth(envelope,50) + + return e*sinal - return e * sinal diff --git a/physiological_data/lib/novainstrumentation/sandbox/sync.py b/physiological_data/lib/novainstrumentation/sandbox/sync.py index 4cb89a8b..b845c64f 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/sync.py +++ b/physiological_data/lib/novainstrumentation/sandbox/sync.py @@ -1,7 +1,6 @@ from pylab import * - -def sync(e1, s1, e2, s2, type="timelinear"): +def sync(e1, s1, e2, s2, type='timelinear'): # type : 'timelinear' - First event vs last # 'locally linear' linear inside a sequential set of events # spline @@ -9,6 +8,7 @@ def sync(e1, s1, e2, s2, type="timelinear"): # Sync s2 to s1 # return the new s2 + ti1 = find(diff(e1) == 1)[0] tf1 = find(diff(e1) == -1)[-1] @@ -16,11 +16,11 @@ def sync(e1, s1, e2, s2, type="timelinear"): tf2 = find(diff(e2) == -1)[-1] -def isync(s1, s2): +def isync(s1,s2): - # considering f1 = 1.0 and f2 = f1 * len(s1)/len(s2) - t1 = arange(0, len(s1)) - t2 = arange(0, len(s1), len(s1) / float(len(s2))) + # considering f1 = 1.0 and f2 = f1 * len(s1)/len(s2) + t1 = arange( 0, len(s1) ) + t2 = arange( 0, len(s1), len(s1)/float(len(s2)) ) # TODO test with scipy spline return interp(t1, t2, s2) @@ -28,44 +28,44 @@ def isync(s1, s2): # signal frequencies f1 = 1000.0 -f2 = 1050.0 # +f2 = 1050.0 # # signal duration -duration = 100 # 60 seconds +duration = 100 # 60 seconds # temporal basis -t1 = arange(0, duration, 1.0 / f1) -t2 = arange(0, duration, 1.0 / f2) +t1 = arange(0, duration, 1.0/f1) +t2 = arange(0, duration, 1.0/f2) # Signals out of phase -s1 = sin(2 * pi * 5 * t1) -s2 = sin(2 * pi * 5 * t2) +s1 = sin(2*pi*5*t1) +s2 = sin(2*pi*5*t2) # Events -e1 = ((t1 % 2) > 1) * 1 -e2 = ((t2 % 2) > 1) * 1 +e1 = ((t1 % 2) > 1 ) * 1 +e2 = ((t2 % 2) > 1 ) * 1 # sync signals -is1 = isync(s1, s2) -is2 = isync(s2, s1) +is1 = isync(s1,s2) +is2 = isync(s2,s1) print(len(is1)) -print(sum(abs(s1 - is1)) / len(s1)) +print(sum(abs(s1-is1))/len(s1)) close("all") # show data f = figure() -plot(s1, "g", label="s1") -plot(s2, "r", label="s2") +plot(s1,'g', label = "s1") +plot(s2,'r', label = "s2") title("s1 and s2") legend() f.show() f = figure() -plot(s2, label="s2", color="black", lw=2) -plot(is2, label="is2", color="g", alpha=0.5, lw=5) +plot(s2, label="s2", color='black' ,lw = 2) +plot(is2, label="is2", color= 'g', alpha = 0.5, lw = 5) title("s2 and is2") legend() f.show() @@ -73,4 +73,4 @@ def isync(s1, s2): f = figure() plot(s2 - is2) title("s2 - is2") -f.show() +f.show() \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/sandbox/testwfdb.py b/physiological_data/lib/novainstrumentation/sandbox/testwfdb.py index 1a45673f..60cc91d5 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/testwfdb.py +++ b/physiological_data/lib/novainstrumentation/sandbox/testwfdb.py @@ -1,19 +1,18 @@ from novainstrumentation.code.sandbox.wfdbtools import rdsamp, rdann, plot_data +#from pprint import pprint + + # Record is a format 212 record from physiobank. + # Note that name of record does not include extension. +record = '100' -# from pprint import pprint - -# Record is a format 212 record from physiobank. -# Note that name of record does not include extension. -record = "100" - -# Read in the data from 0 to 10 seconds + # Read in the data from 0 to 10 seconds data, info = rdsamp(record, 0, 10) - -# returned data is an array. The columns are time(samples), -# time(seconds), signal1, signal2 + + # returned data is an array. The columns are time(samples), + # time(seconds), signal1, signal2 print(data.shape) -# info is a dictionary containing header information + # info is a dictionary containing header information print(info) # {'first_values': [995, 1011], # 'gains': [200, 200], @@ -21,15 +20,15 @@ # 'samp_freq': 360, # 'signal_names': ['MLII', 'V5'], # 'zero_values': [1024, 1024]} + +ann = rdann(record, 'atr', 0, 10) -ann = rdann(record, "atr", 0, 10) - -print(ann[:4, :]) +print(ann[:4,:]) # array([[ 18. , 0.05 , 28. ], # [ 77. , 0.214, 1. ], # [ 370. , 1.028, 1. ], # [ 662. , 1.839, 1. ]]) -# - -# Plot the data and the mark the annotations +# + + # Plot the data and the mark the annotations plot_data(data, info, ann) diff --git a/physiological_data/lib/novainstrumentation/sandbox/testwfdbsyntrh.py b/physiological_data/lib/novainstrumentation/sandbox/testwfdbsyntrh.py index 05f27169..e9555fd4 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/testwfdbsyntrh.py +++ b/physiological_data/lib/novainstrumentation/sandbox/testwfdbsyntrh.py @@ -4,20 +4,20 @@ import novainstrumentation as ni -# from pprint import pprint +#from pprint import pprint + + # Record is a format 212 record from physiobank. + # Note that name of record does not include extension. +record = '100' -# Record is a format 212 record from physiobank. -# Note that name of record does not include extension. -record = "100" - -# Read in the data from 0 to 10 seconds + # Read in the data from 0 to 10 seconds data, info = rdsamp(record, 0, 10) - -# returned data is an array. The columns are time(samples), -# time(seconds), signal1, signal2 + + # returned data is an array. The columns are time(samples), + # time(seconds), signal1, signal2 print(data.shape) -# info is a dictionary containing header information + # info is a dictionary containing header information print(info) # {'first_values': [995, 1011], # 'gains': [200, 200], @@ -25,38 +25,41 @@ # 'samp_freq': 360, # 'signal_names': ['MLII', 'V5'], # 'zero_values': [1024, 1024]} + +ann = rdann(record, 'atr', 0, 10) -ann = rdann(record, "atr", 0, 10) - -print(ann[:4, :]) +print(ann[:4,:]) # array([[ 18. , 0.05 , 28. ], # [ 77. , 0.214, 1. ], # [ 370. , 1.028, 1. ], # [ 662. , 1.839, 1. ]]) -# +# + + # Plot the data and the mark the annotations +#plot_data(data, info, ann) -# Plot the data and the mark the annotations -# plot_data(data, info, ann) +s = data[936+100:1222+100,2] -s = data[936 + 100 : 1222 + 100, 2] +duration = 30 # in seconds -duration = 30 # in seconds +signal, beats = ni.synthbeats(duration, meanhr = 60, stdhr = 3, samplingfreq = 500) -signal, beats = ni.synthbeats(duration, meanhr=60, stdhr=3, samplingfreq=500) +time = np.arange(len(signal))/250. -time = np.arange(len(signal)) / 250.0 +ns=np.convolve(signal,s+0.4) -ns = np.convolve(signal, s + 0.4) +s=ns+pl.randn(len(ns))*0.01 -s = ns + pl.randn(len(ns)) * 0.01 +s = s+sin(arange(len(ns))/1000.+pl.randn(len(ns))/5.)*.1 -s = s + sin(arange(len(ns)) / 1000.0 + pl.randn(len(ns)) / 5.0) * 0.1 pl.plot(s) + pl.draw() -# pl.show() +#pl.show() + diff --git a/physiological_data/lib/novainstrumentation/sandbox/wfdbtools.py b/physiological_data/lib/novainstrumentation/sandbox/wfdbtools.py index a1ceffb5..f8efd016 100644 --- a/physiological_data/lib/novainstrumentation/sandbox/wfdbtools.py +++ b/physiological_data/lib/novainstrumentation/sandbox/wfdbtools.py @@ -1,4 +1,4 @@ -# @PydevCodeAnalysisIgnore +#@PydevCodeAnalysisIgnore #!/usr/bin/env python @@ -23,7 +23,7 @@ """ -A pure python module for accessing and using the waveform data in `Physiobank `_. +A pure python module for accessing and using the waveform data in `Physiobank `_. Provides `rdsamp` and `rdann` which are the python equivalents of the wfdb applications of similar names. A deliberate attempt is made to try to keep names and usage similar to the @@ -37,14 +37,14 @@ >> from wfdbtools import rdsamp, rdann, plot_data >> from pprint import pprint - + # Record is a format 212 record from physiobank. # Note that name of record does not include extension. >> record = 'samples/format212/100' # Read in the data from 0 to 10 seconds >> data, info = rdsamp(record, 0, 10) - + # returned data is an array. The columns are time(samples), # time(seconds), signal1, signal2 >> print data.shape @@ -58,7 +58,7 @@ 'samp_freq': 360, 'signal_names': ['MLII', 'V5'], 'zero_values': [1024, 1024]} - + # And now read the annotation >> ann = rdann(record, 'atr', 0, 10) @@ -68,8 +68,8 @@ [ 77. , 0.214, 1. ], [ 370. , 1.028, 1. ], [ 662. , 1.839, 1. ]]) - - + + # Plot the data and the mark the annotations >> plot_data(data, info, ann) @@ -85,57 +85,57 @@ import numpy -# import pylab -# from pprint import pprint +#import pylab +#from pprint import pprint -__author__ = "Raja Selvaraj " -__version__ = "0.2" +__author__ = 'Raja Selvaraj ' +__version__ = '0.2' ## Annotation codes CODEDICT = { - 0: "NOTQRS", # not-QRS (not a getann/putann codedict) */ - 1: "NORMAL", # normal beat */ - 2: "LBBB", # left bundle branch block beat */ - 3: "RBBB", # right bundle branch block beat */ - 4: "ABERR", # aberrated atrial premature beat */ - 5: "PVC", # premature ventricular contraction */ - 6: "FUSION", # fusion of ventricular and normal beat */ - 7: "NPC", # nodal (junctional) premature beat */ - 8: "APC", # atrial premature contraction */ - 9: "SVPB", # premature or ectopic supraventricular beat */ - 10: "VESC", # ventricular escape beat */ - 11: "NESC", # nodal (junctional) escape beat */ - 12: "PACE", # paced beat */ - 13: "UNKNOWN", # unclassifiable beat */ - 14: "NOISE", # signal quality change */ - 16: "ARFCT", # isolated QRS-like artifact */ - 18: "STCH", # ST change */ - 19: "TCH", # T-wave change */ - 20: "SYSTOLE", # systole */ - 21: "DIASTOLE", # diastole */ - 22: "NOTE", # comment annotation */ - 23: "MEASURE", # measurement annotation */ - 24: "PWAVE", # P-wave peak */ - 25: "BBB", # left or right bundle branch block */ - 26: "PACESP", # non-conducted pacer spike */ - 27: "TWAVE", # T-wave peak */ - 28: "RHYTHM", # rhythm change */ - 29: "UWAVE", # U-wave peak */ - 30: "LEARN", # learning */ - 31: "FLWAV", # ventricular flutter wave */ - 32: "VFON", # start of ventricular flutter/fibrillation */ - 33: "VFOFF", # end of ventricular flutter/fibrillation */ - 34: "AESC", # atrial escape beat */ - 35: "SVESC", # supraventricular escape beat */ - 36: "LINK", # link to external data (aux contains URL) */ - 37: "NAPC", # non-conducted P-wave (blocked APB) */ - 38: "PFUS", # fusion of paced and normal beat */ - 39: "WFON", # waveform onset */ - # WFON : 'PQ', # PQ junction (beginning of QRS) */ - 40: "WFOFF", # waveform end */ - # WFOFF : 'JPT', # J point (end of QRS) */ - 41: "RONT", # R-on-T premature ventricular contraction */ -} + 0 : 'NOTQRS', # not-QRS (not a getann/putann codedict) */ + 1 : 'NORMAL', # normal beat */ + 2 : 'LBBB', # left bundle branch block beat */ + 3 : 'RBBB', # right bundle branch block beat */ + 4 : 'ABERR', # aberrated atrial premature beat */ + 5 : 'PVC', # premature ventricular contraction */ + 6 : 'FUSION', # fusion of ventricular and normal beat */ + 7 : 'NPC', # nodal (junctional) premature beat */ + 8 : 'APC', # atrial premature contraction */ + 9 : 'SVPB', # premature or ectopic supraventricular beat */ + 10 : 'VESC', # ventricular escape beat */ + 11 : 'NESC', # nodal (junctional) escape beat */ + 12 : 'PACE', # paced beat */ + 13 : 'UNKNOWN', # unclassifiable beat */ + 14 : 'NOISE', # signal quality change */ + 16 : 'ARFCT', # isolated QRS-like artifact */ + 18 : 'STCH', # ST change */ + 19 : 'TCH', # T-wave change */ + 20 : 'SYSTOLE', # systole */ + 21 : 'DIASTOLE', # diastole */ + 22 : 'NOTE', # comment annotation */ + 23 : 'MEASURE', # measurement annotation */ + 24 : 'PWAVE', # P-wave peak */ + 25 : 'BBB', # left or right bundle branch block */ + 26 : 'PACESP', # non-conducted pacer spike */ + 27 : 'TWAVE', # T-wave peak */ + 28 : 'RHYTHM', # rhythm change */ + 29 : 'UWAVE', # U-wave peak */ + 30 : 'LEARN', # learning */ + 31 : 'FLWAV', # ventricular flutter wave */ + 32 : 'VFON', # start of ventricular flutter/fibrillation */ + 33 : 'VFOFF', # end of ventricular flutter/fibrillation */ + 34 : 'AESC', # atrial escape beat */ + 35 : 'SVESC', # supraventricular escape beat */ + 36 : 'LINK', # link to external data (aux contains URL) */ + 37 : 'NAPC', # non-conducted P-wave (blocked APB) */ + 38 : 'PFUS', # fusion of paced and normal beat */ + 39 : 'WFON', # waveform onset */ + #WFON : 'PQ', # PQ junction (beginning of QRS) */ + 40 : 'WFOFF', # waveform end */ + #WFOFF : 'JPT', # J point (end of QRS) */ + 41 : 'RONT' # R-on-T premature ventricular contraction */ + } def rdsamp(record, start=0, end=-1, interval=-1): @@ -145,7 +145,7 @@ def rdsamp(record, start=0, end=-1, interval=-1): Only 2 channel records in format 212 are supported. This is the most common record in the Physionet database(http://www.physionet.org/physiobank/). - + Parameters ---------- record : str @@ -165,7 +165,7 @@ def rdsamp(record, start=0, end=-1, interval=-1): col 1 - Elapsed time in samples col 2 - Elapsed time in milliseconds col 3,4 - The two signals - Signal amplitude is in physical units (mV) + Signal amplitude is in physical units (mV) info : dict Dictionary containing header information keys : @@ -175,17 +175,16 @@ def rdsamp(record, start=0, end=-1, interval=-1): 'first_values' - First value of each signal 'gains' - Gain for each signal 'zero_values' - Zero value for each signal - + """ # read the header file - output is a dict info = rdhdr(record) # establish start and end in samples start, end = _get_read_limits(start, end, interval, info) # read the data - data = _read_data(record, start, end, info) + data = _read_data(record, start, end, info) return data, info - def rdann(record, annotator, start=0, end=-1, types=[]): """ Reads annotations for given record by specified annotator. @@ -206,7 +205,7 @@ def rdann(record, annotator, start=0, end=-1, types=[]): Types are input as annotation code (integer from 0 to 49) Annotation types not in list will be ignored. Default is empty list, which results in all types being read. - + Returns ------- data : (N, 3) ndarray @@ -218,10 +217,10 @@ def rdann(record, annotator, start=0, end=-1, types=[]): """ # get header data info = rdhdr(record) - - annfile = "".join((record, ".", annotator)) - with open(annfile, "rb") as f: - arr = numpy.fromstring(f.read(), dtype=numpy.uint8).reshape((-1, 2)) + + annfile = ''.join((record, '.', annotator)) + with open(annfile, 'rb') as f: + arr = numpy.fromstring(f.read(), dtype = numpy.uint8).reshape((-1, 2)) rows = arr.shape[0] annot = [] @@ -231,14 +230,9 @@ def rdann(record, annotator, start=0, end=-1, types=[]): while i < rows: anntype = arr[i, 1] >> 2 if anntype == 59: - annot.append(arr[i + 3, 1] >> 2) - annot_time.append( - arr[i + 2, 0] - + (arr[i + 2, 1] << 8) - + (arr[i + 1, 0] << 16) - + arr[i + 1, 1] - << 24 - ) + annot.append(arr[i+3, 1] >> 2) + annot_time.append(arr[i+2, 0] + (arr[i+2, 1] << 8) + + (arr[i+1, 0] << 16) + arr[i+1, 1] << 24) i += 3 elif anntype in [60, 61, 62]: pass @@ -254,27 +248,26 @@ def rdann(record, annotator, start=0, end=-1, types=[]): # last values are EOF indicator annot_time = annot_time[:-1] annot = annot[:-1] - + # annot_time should be total elapsed samples annot_time = numpy.cumsum(annot_time) - annot_time_ms = annot_time / info["samp_freq"] # in seconds - + annot_time_ms = annot_time / info['samp_freq'] # in seconds + # limit to requested interval start, end = _get_read_limits(start, end, -1, info) ann = numpy.array([annot_time, annot_time_ms, annot]).transpose() - + # filter by annot_time in interval - ann = ann[start <= ann[:, 0]] + ann = ann[start <= ann[:, 0]] ann = ann[ann[:, 0] <= end] # filter by type if types != []: ann = ann[numpy.logical_or.reduce([ann[:, 2] == x for x in types])] - # ann = ann[numpy.array([ann[x, 2] in types for x in range(len(ann))])] + #ann = ann[numpy.array([ann[x, 2] in types for x in range(len(ann))])] return ann - - + def plot_data(data, info, ann=None): """ Plot the signal with annotations if available. @@ -299,27 +292,27 @@ def plot_data(data, info, ann=None): except ImportError: warnings.warn("Could not import pylab. Abandoning plotting") return - - time = data[:, 1] # in seconds. use data[:, 0] to use sample no. + + time = data[:, 1] #in seconds. use data[:, 0] to use sample no. sig1 = data[:, 2] sig2 = data[:, 3] - + pylab.subplot(211) - pylab.plot(time, sig1, "k") + pylab.plot(time, sig1, 'k') pylab.xticks([]) - pylab.ylabel("%s (mV)" % (info["signal_names"][0])) - + pylab.ylabel('%s (mV)' %(info['signal_names'][0])) + pylab.subplot(212) - pylab.plot(time, sig2, "k") - pylab.ylabel("%s (mV)" % (info["signal_names"][1])) - pylab.xlabel("Time (seconds)") + pylab.plot(time, sig2, 'k') + pylab.ylabel('%s (mV)' %(info['signal_names'][1])) + pylab.xlabel('Time (seconds)') if ann != None: # annotation time in samples from start - ann_x = (ann[:, 0] - data[0, 0]).astype("int") - pylab.plot(ann[:, 1], data[ann_x, 3], "xr") + ann_x = (ann[:, 0] - data[0, 0]).astype('int') + pylab.plot(ann[:, 1], data[ann_x, 3], 'xr') pylab.subplot(211) - pylab.plot(ann[:, 1], data[ann_x, 2], "xr") + pylab.plot(ann[:, 1], data[ann_x, 2], 'xr') pylab.show() @@ -348,122 +341,88 @@ def rdhdr(record): 'gains' - Gain for each signal 'zero_values' - Zero value for each signal 'signal_names' - Name/Descr for each signal - + """ - info = { - "signal_names": [], - "gains": [], - "units": [], - "first_values": [], - "zero_values": [], - } - - RECORD_REGEX = re.compile( - r"".join( - [ - "(?P\d+)\/*(?P\d*)\s", - "(?P\d+)\s*", - "(?P\d*)\/?(?P\d*)\(?(?P\d*)\)?\s*", - "(?P\d*)\s*", - "(?P\d{,2}:*\d{,2}:*\d{,2})\s*", - "(?P\d{,2}\/*\d{,2}\/*\d{,4})", - ] - ) - ) - - SIGNAL_REGEX = re.compile( - r"".join( - [ - "(?P[0-9a-zA-Z\._/-]+)\s+", - "(?P\d+)x{,1}(?P\d*):*", - "(?P\d*)\+*(?P\d*)\s*", - "(?P\d*)\(?(?P\d*)\)?\/?", - "(?P\w*)\s*(?P\d*)\s*", - "(?P\d*)\s*(?P[\d-]*)\s*", - "(?P[0-9-]*)\s*(?P\d*)\s*", - "(?P[a-zA-Z0-9\s]*)", - ] - ) - ) + info = {'signal_names':[], 'gains':[], 'units':[], + 'first_values':[], 'zero_values':[]} + + RECORD_REGEX = re.compile(r''.join([ + "(?P\d+)\/*(?P\d*)\s", + "(?P\d+)\s*", + "(?P\d*)\/?(?P\d*)\(?(?P\d*)\)?\s*", + "(?P\d*)\s*", + "(?P\d{,2}:*\d{,2}:*\d{,2})\s*", + "(?P\d{,2}\/*\d{,2}\/*\d{,4})"])) + + SIGNAL_REGEX = re.compile(r''.join([ + "(?P[0-9a-zA-Z\._/-]+)\s+", + "(?P\d+)x{,1}(?P\d*):*", + "(?P\d*)\+*(?P\d*)\s*", + "(?P\d*)\(?(?P\d*)\)?\/?", + "(?P\w*)\s*(?P\d*)\s*", + "(?P\d*)\s*(?P[\d-]*)\s*", + "(?P[0-9-]*)\s*(?P\d*)\s*", + "(?P[a-zA-Z0-9\s]*)"])) header_lines, comment_lines = _getheaderlines(record) - ( - record_name, - seg_count, - signal_count, - samp_freq, - counter_freq, - base_counter, - samp_count, - base_time, - base_date, - ) = RECORD_REGEX.findall(header_lines[0])[0] + (record_name, seg_count, signal_count, samp_freq, + counter_freq, base_counter, samp_count, + base_time, base_date) = RECORD_REGEX.findall(header_lines[0])[0] # use 250 if missing - if samp_freq == "": + if samp_freq == '': samp_freq = 250 - if samp_count == "": + if samp_count == '': samp_count = 0 - - info["samp_freq"] = float(samp_freq) - info["samp_count"] = int(samp_count) - + + + info['samp_freq'] = float(samp_freq) + info['samp_count'] = int(samp_count) + for sig in range(2): - ( - file_name, - file_format, - samp_per_frame, - skew, - byte_offset, - gain, - baseline, - units, - resolution, - zero_value, - first_value, - checksum, - blocksize, - signal_name, - ) = SIGNAL_REGEX.findall(header_lines[sig + 1])[0] + (file_name, file_format, samp_per_frame, skew, + byte_offset, gain, baseline, units, + resolution, zero_value, first_value, + checksum, blocksize, signal_name) = SIGNAL_REGEX.findall( + header_lines[sig+1])[0] # replace with defaults for missing values - if gain == "" or gain == 0: + if gain == '' or gain == 0: gain = 200 - if units == "": - units = "mV" - if zero_value == "": + if units == '': + units = 'mV' + if zero_value == '': zero_value = 0 - if first_value == "": - first_value = 0 # do not use to check - - info["gains"].append(float(gain)) - info["units"].append(units) - info["zero_values"].append(float(zero_value)) - info["first_values"].append(float(first_value)) - info["signal_names"].append(signal_name) + if first_value == '': + first_value = 0 # do not use to check + + info['gains'].append(float(gain)) + info['units'].append(units) + info['zero_values'].append(float(zero_value)) + info['first_values'].append(float(first_value)) + info['signal_names'].append(signal_name) return info - - + def _getheaderlines(record): """Read the header file and separate comment lines and header lines""" - hfile = record + ".hea" - all_lines = open(hfile, "r").readlines() + hfile = record + '.hea' + all_lines = open(hfile, 'r').readlines() comment_lines = [] header_lines = [] # strip newlines - all_lines = [l.rstrip("\n").rstrip("\r") for l in all_lines] + all_lines = [l.rstrip('\n').rstrip('\r') for l in all_lines] # comments for l in all_lines: - if l.startswith("#"): + if l.startswith('#'): comment_lines.append(l) - elif l.strip() != "": + elif l.strip() != '': header_lines.append(l) - + return header_lines, comment_lines - + def _get_read_limits(start, end, interval, info): """ Given start time, end time and interval @@ -485,76 +444,73 @@ def _get_read_limits(start, end, interval, info): >>> _get_read_limits(-1, -1, -1, {'samp_count':100, 'samp_freq':10}) (0, 100) """ - start *= info["samp_freq"] - end *= info["samp_freq"] - - if start < 0: # If start is negative, start at 0 + start *= info['samp_freq'] + end *= info['samp_freq'] + + if start < 0: # If start is negative, start at 0 start = 0 - if end < 0: # if end is negative, use end of record - end = info["samp_count"] - if end < start: # if end is before start, swap them + if end < 0: # if end is negative, use end of record + end = info['samp_count'] + if end < start: # if end is before start, swap them start, end = end, start - interval_end = start + interval * info["samp_freq"] # end det by interval + interval_end = start + interval * info['samp_freq'] # end det by interval if interval_end < start: - interval_end = info["samp_count"] - end = min(end, interval_end, info["samp_count"]) # use earlier end + interval_end = info['samp_count'] + end = min(end, interval_end, info['samp_count']) # use earlier end return int(start), int(end) - - + def _read_data(record, start, end, info): """Read the binary data for each signal""" - datfile = record + ".dat" + datfile = record + '.dat' samp_to_read = end - start # verify against first value in header - with open(datfile, "rb") as f: - data = _arr_to_data( - numpy.fromstring(f.read(3), dtype=numpy.uint8).reshape(1, 3) - ) - - if [data[0, 2], data[0, 3]] != info["first_values"]: - warnings.warn("First value from dat file does not match value in header") - + with open(datfile, 'rb') as f: + data = _arr_to_data(numpy.fromstring(f.read(3), + dtype=numpy.uint8).reshape(1,3)) + + if [data[0, 2], data[0, 3]] != info['first_values']: + warnings.warn( + 'First value from dat file does not match value in header') + # read into an array with 3 bytes in each row - with open(datfile, "rb") as f: - f.seek(start * 3) - arr = numpy.fromstring(f.read(3 * samp_to_read), dtype=numpy.uint8).reshape( - (samp_to_read, 3) - ) + with open(datfile, 'rb') as f: + f.seek(start*3) + arr = numpy.fromstring(f.read(3*samp_to_read), + dtype=numpy.uint8).reshape((samp_to_read, 3)) data = _arr_to_data(arr) # adjust zerovalue and gain - data[:, 2] = (data[:, 2] - info["zero_values"][0]) / info["gains"][0] - data[:, 3] = (data[:, 3] - info["zero_values"][1]) / info["gains"][1] + data[:, 2] = (data[:, 2] - info['zero_values'][0]) / info['gains'][0] + data[:, 3] = (data[:, 3] - info['zero_values'][1]) / info['gains'][1] # time columns data[:, 0] = numpy.arange(start, end) # elapsed time in samples - data[:, 1] = (numpy.arange(samp_to_read) + start) / info["samp_freq"] # in sec + data[:, 1] = (numpy.arange(samp_to_read) + start + ) / info['samp_freq'] # in sec return data - def _arr_to_data(arr): """From the numpy array read from the dat file using bit level operations, extract the 12-bit data""" - second_col = arr[:, 1].astype("int") - bytes1 = second_col & 15 # bytes belonging to first sample - bytes2 = second_col >> 4 # belongs to second sample - sign1 = (second_col & 8) << 9 # sign bit for first sample - sign2 = (second_col & 128) << 5 # sign bit for second sample + second_col = arr[:, 1].astype('int') + bytes1 = second_col & 15 # bytes belonging to first sample + bytes2 = second_col >> 4 # belongs to second sample + sign1 = (second_col & 8) << 9 # sign bit for first sample + sign2 = (second_col & 128) << 5 # sign bit for second sample # data has columns - samples, time(ms), signal1 and signal2 - data = numpy.zeros((arr.shape[0], 4), dtype="float") + data = numpy.zeros((arr.shape[0], 4), dtype='float') data[:, 2] = (bytes1 << 8) + arr[:, 0] - sign1 data[:, 3] = (bytes2 << 8) + arr[:, 2] - sign2 return data - def get_annotation_code(code=None): """Returns the symbolic definition for the wfdb annotation code. See http://www.physionet.org/physiotools/wpg/wpg_31.htm for details. Based on ecgcodes.h from wfdb. - + Parameters ---------- code : int @@ -568,16 +524,14 @@ def get_annotation_code(code=None): """ return CODEDICT[code] - def main(): """Run tests when called directly""" import nose - print("-----------------") print("Running the tests") print("-----------------") nose.main() - -if __name__ == "__main__": +if __name__ == '__main__': main() + diff --git a/physiological_data/lib/novainstrumentation/smooth.py b/physiological_data/lib/novainstrumentation/smooth.py index 7c927d74..b06e397b 100644 --- a/physiological_data/lib/novainstrumentation/smooth.py +++ b/physiological_data/lib/novainstrumentation/smooth.py @@ -1,7 +1,7 @@ import numpy -def smooth(input_signal, window_len=10, window="hanning"): +def smooth(input_signal, window_len=10, window='hanning'): """ @brief: Smooth the data using a window with requested size. @@ -48,23 +48,19 @@ def smooth(input_signal, window_len=10, window="hanning"): if window_len < 3: return input_signal - if window not in ["flat", "hanning", "hamming", "bartlett", "blackman"]: - raise ValueError( - """Window is on of 'flat', 'hanning', 'hamming', -'bartlett', 'blackman'""" - ) + if window not in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: + raise ValueError("""Window is on of 'flat', 'hanning', 'hamming', +'bartlett', 'blackman'""") - sig = numpy.r_[ - 2 * input_signal[0] - input_signal[window_len:0:-1], - input_signal, - 2 * input_signal[-1] - input_signal[-2 : -window_len - 2 : -1], - ] + sig = numpy.r_[2 * input_signal[0] - input_signal[window_len:0:-1], + input_signal, + 2 * input_signal[-1] - input_signal[-2:-window_len-2:-1]] - if window == "flat": # moving average - win = numpy.ones(window_len, "d") + if window == 'flat': # moving average + win = numpy.ones(window_len, 'd') else: - win = eval("numpy." + window + "(window_len)") + win = eval('numpy.' + window + '(window_len)') - sig_conv = numpy.convolve(win / win.sum(), sig, mode="same") + sig_conv = numpy.convolve(win / win.sum(), sig, mode='same') - return sig_conv[window_len:-window_len] + return sig_conv[window_len: -window_len] diff --git a/physiological_data/lib/novainstrumentation/stat_analysis.py b/physiological_data/lib/novainstrumentation/stat_analysis.py index 57713701..863a7509 100644 --- a/physiological_data/lib/novainstrumentation/stat_analysis.py +++ b/physiological_data/lib/novainstrumentation/stat_analysis.py @@ -5,15 +5,15 @@ def scott_density(signal): - # Based on seaborn 0.4 + #Based on seaborn 0.4 s_k = np.asarray(signal, np.float) - bw, widths, gridsize, cut = "scott", 0.8, 100, 3 + bw, widths, gridsize, cut = 'scott', .8, 100, 3 kde = gaussian_kde(s_k, bw) if isinstance(bw, str): - bw = "scotts" if bw == "scott" else bw - bw = getattr(kde, "%s_factor" % bw)() + bw = 'scotts' if bw == 'scott' else bw + bw = getattr(kde, '%s_factor' % bw)() support_min = max(min(s_k) - bw * cut, -np.inf) - support_max = min(max(s_k) + bw * cut, np.inf) + support_max = min(max(s_k) + bw * cut, np.inf) y = np.linspace(support_min, support_max, gridsize) dens = kde.evaluate(y) scl = 1 / (dens.max() / (widths / 2)) @@ -22,13 +22,13 @@ def scott_density(signal): return dens, y -def matrix_recurrence(signal): +def matrix_recurrence(signal): N = len(signal) - S = np.ones((N, N)) + S = np.ones((N,N)) for j in np.arange(N): for i in np.arange(N): - ij = abs(signal[j] - signal[i]) - if ij <= 0.1 * np.std(signal): - S[i, j] = 0 - - return S + ij = abs(signal[j]-signal[i]) + if ij <= 0.1*np.std(signal): + S[i,j] = 0 + + return S \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/tests/__init__.py b/physiological_data/lib/novainstrumentation/tests/__init__.py index d918fc63..c9ac9ee7 100644 --- a/physiological_data/lib/novainstrumentation/tests/__init__.py +++ b/physiological_data/lib/novainstrumentation/tests/__init__.py @@ -2,4 +2,4 @@ from novainstrumentation.test_peaks import * from novainstrumentation.test_smooth import * from novainstrumentation.test_tools import * -from novainstrumentation.test_waves import * +from novainstrumentation.test_waves import * \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/tests/test_filter.py b/physiological_data/lib/novainstrumentation/tests/test_filter.py index 7df48893..7b9c067d 100644 --- a/physiological_data/lib/novainstrumentation/tests/test_filter.py +++ b/physiological_data/lib/novainstrumentation/tests/test_filter.py @@ -2,7 +2,8 @@ import numpy from scipy import signal -from sklearn.utils.testing import assert_array_almost_equal, assert_true, assert_less +from sklearn.utils.testing import assert_array_almost_equal, \ + assert_true, assert_less import novainstrumentation import novainstrumentation as ni @@ -12,30 +13,29 @@ base_dir = os.path.dirname(novainstrumentation.__file__) - def test_emg_filter_params(): - fname = base_dir + "/code/data/emg.txt" + fname = base_dir + '/code/data/emg.txt' t, emg = numpy.loadtxt(fname) - env_ni = ni.lowpass(abs(emg), order=2, f=10, fs=1000) + env_ni = ni.lowpass(abs(emg), order=2, f=10,fs=1000) - b, a = signal.butter(2, Wn=10 / (1000.0 / 2.0)) - env_ref = signal.lfilter(b, a, abs(emg)) + b,a = signal.butter(2,Wn=10/(1000.0/2.0)) + env_ref = signal.lfilter(b,a,abs(emg)) - assert_array_almost_equal(env_ni, env_ref) + assert_array_almost_equal(env_ni,env_ref) def test_emg_filter_values(): - f_name = base_dir + "/code/data/emg.txt" + f_name = base_dir + '/code/data/emg.txt' t, emg = numpy.loadtxt(f_name) - env_ni = ni.lowpass(abs(emg), order=2, f=10, fs=1000.0) + env_ni = ni.lowpass(abs(emg),order=2,f=10,fs=1000.) - b, a = signal.butter(2, Wn=10 / (1000.0 / 2.0)) - env_ref = signal.lfilter(b, a, abs(emg)) + b,a = signal.butter(2,Wn=10/(1000.0/2.0)) + env_ref = signal.lfilter(b,a,abs(emg)) - assert_array_almost_equal(env_ni, env_ref) + assert_array_almost_equal(env_ni,env_ref) -# def test_fail(): +#def test_fail(): # assert False diff --git a/physiological_data/lib/novainstrumentation/tests/test_peaks.py b/physiological_data/lib/novainstrumentation/tests/test_peaks.py index 23e878db..1ce71b5c 100644 --- a/physiological_data/lib/novainstrumentation/tests/test_peaks.py +++ b/physiological_data/lib/novainstrumentation/tests/test_peaks.py @@ -5,64 +5,31 @@ from novainstrumentation import peaks, peakdelta -path = ".." + os.path.sep + "data" + os.path.sep +path = '..' + os.path.sep + 'data' + os.path.sep class TestPeaks(TestCase): def test_peaks_cleanecg(self): - t, s = loadtxt(path + "cleanecg.txt") - pr = array( - [20, 87, 156, 225, 291, 355, 418, 482, 550, 624, 694, 764, 834, 905, 978] - ) + t, s = loadtxt(path + 'cleanecg.txt') + pr = array([20, 87, 156, 225, 291, 355, 418, 482, 550, 624, 694, 764, 834, 905, 978]) assert_array_equal(pr, peaks(s, 75)) def test_peaks_eda(self): - t, s = loadtxt(path + "eda.txt") + t, s = loadtxt(path + 'eda.txt') pr = array([95, 146, 346]) assert_array_equal(pr, peaks(s, 600)) class TestPeakdelta(TestCase): def test_peakdelta_cleanecg(self): - t, s = loadtxt(path + "cleanecg.txt") - pr = array( - [ - 20, - 40, - 87, - 107, - 156, - 175, - 225, - 244, - 291, - 311, - 355, - 375, - 418, - 438, - 482, - 501, - 550, - 569, - 624, - 644, - 694, - 713, - 764, - 784, - 834, - 854, - 905, - 925, - 978, - ] - ) + t, s = loadtxt(path + 'cleanecg.txt') + pr = array([20, 40, 87, 107, 156, 175, 225, 244, 291, 311, 355, 375, 418, 438, 482, + 501, 550, 569, 624, 644, 694, 713, 764, 784, 834, 854, 905, 925, 978]) assert_array_equal(pr, peakdelta(s, 50)[0][:, 0].astype(int)) def test_peakdelta_eda(self): - t, s = loadtxt(path + "eda.txt") + t, s = loadtxt(path + 'eda.txt') pks = peakdelta(s, 50) prmax = array([0, 95, 278, 346]) @@ -71,6 +38,5 @@ def test_peakdelta_eda(self): assert_array_equal(prmax, pks[0][:, 0].astype(int)) assert_array_equal(prmin, pks[1][:, 0].astype(int)) - if __name__ == "__main__": run_module_suite() diff --git a/physiological_data/lib/novainstrumentation/tests/test_smooth.py b/physiological_data/lib/novainstrumentation/tests/test_smooth.py index c6be891c..f9292ee4 100644 --- a/physiological_data/lib/novainstrumentation/tests/test_smooth.py +++ b/physiological_data/lib/novainstrumentation/tests/test_smooth.py @@ -5,10 +5,9 @@ def test_symmetric_window(): - x = arange(1.0, 11.0) - sx = smooth(x, window_len=5, window="flat") + x = arange(1., 11.) + sx = smooth(x, window_len=5, window='flat') assert_allclose(x, sx) - if __name__ == "__main__": run_module_suite() diff --git a/physiological_data/lib/novainstrumentation/tests/test_tools.py b/physiological_data/lib/novainstrumentation/tests/test_tools.py index a9c9775c..db092ebb 100644 --- a/physiological_data/lib/novainstrumentation/tests/test_tools.py +++ b/physiological_data/lib/novainstrumentation/tests/test_tools.py @@ -1,7 +1,8 @@ import os import numpy -from sklearn.utils.testing import assert_array_almost_equal, assert_true, assert_less +from sklearn.utils.testing import assert_array_almost_equal, \ + assert_true, assert_less import novainstrumentation from novainstrumentation import load_with_cache @@ -13,20 +14,20 @@ def test_load_with_cache(): - fname = base_dir + "/code/data/bvp.txt" - fname_npy = fname + ".npy" + fname = base_dir + '/code/data/bvp.txt' + fname_npy = fname + '.npy' direct_d = numpy.loadtxt(fname) - cache_d = load_with_cache(fname, data_type="float") + cache_d = load_with_cache(fname, data_type='float') - # Check that the data was correctly read + #Check that the data was correctly read assert_array_almost_equal(direct_d, cache_d, decimal=5) - # Guarantee that the file was created + #Guarantee that the file was created assert_true(os.path.exists(fname_npy)) - # The file should be smaller that the original + #The file should be smaller that the original assert_less(os.path.getsize(fname_npy), os.path.getsize(fname)) os.remove(fname_npy) diff --git a/physiological_data/lib/novainstrumentation/tests/test_waves.py b/physiological_data/lib/novainstrumentation/tests/test_waves.py index 2635dbf3..7b3c3e80 100644 --- a/physiological_data/lib/novainstrumentation/tests/test_waves.py +++ b/physiological_data/lib/novainstrumentation/tests/test_waves.py @@ -2,7 +2,8 @@ import matplotlib.pyplot as plt import numpy -from sklearn.utils.testing import assert_array_almost_equal, assert_true, assert_less +from sklearn.utils.testing import assert_array_almost_equal, \ + assert_true, assert_less import novainstrumentation from novainstrumentation.waves import stdwave @@ -10,65 +11,27 @@ base_dir = os.path.dirname(novainstrumentation.__file__) - -def test_stdwaves(show_graph=0): - fname = base_dir + "/code/data/xyzcal.txt" +def test_stdwaves(show_graph = 0): + fname = base_dir + '/code/data/xyzcal.txt' t, xcal, ycal, zcal = numpy.loadtxt(fname) if show_graph == 1: - plt.plot(t, xcal, "r-", t, ycal, "b-", t, zcal, "k-") + plt.plot(t, xcal, 'r-', t, ycal, 'b-', t, zcal, 'k-') plt.show() test_matrix = numpy.vstack([t, xcal, ycal, zcal]).transpose() - reference_values = numpy.std(test_matrix, 0) + reference_values = numpy.std(test_matrix,0) test_values = stdwave(test_matrix) - assert_array_almost_equal( - test_values, - reference_values, - 6, - "\n ERROR: waves.stdwaves function not working properly", - ) - + assert_array_almost_equal(test_values, reference_values, 6 ,"\n ERROR: waves.stdwaves function not working properly") def test_waves(): - t, s = numpy.loadtxt(base_dir + "/code/data/cleanecg.txt") - test_matrix = numpy.array( - [ - 20, - 40, - 87, - 107, - 156, - 175, - 225, - 244, - 291, - 311, - 355, - 375, - 418, - 438, - 482, - 501, - 550, - 569, - 624, - 644, - 694, - 713, - 764, - 784, - 834, - 854, - 905, - 925, - 978, - ] - ) - x = waves(s, test_matrix, -10, 100) - plt.plot(x[0, :]) + t, s = numpy.loadtxt(base_dir + '/code/data/cleanecg.txt') + test_matrix = numpy.array([20, 40, 87, 107, 156, 175, 225, 244, 291, 311, 355, 375, 418, 438, 482, + 501, 550, 569, 624, 644, 694, 713, 764, 784, 834, 854, 905, 925, 978]) + x = waves(s, test_matrix, -10, 100) + plt.plot(x[0,:]) plt.show() - test_waves() + diff --git a/physiological_data/lib/novainstrumentation/tictac.py b/physiological_data/lib/novainstrumentation/tictac.py index ee550080..c0e5ca3f 100644 --- a/physiological_data/lib/novainstrumentation/tictac.py +++ b/physiological_data/lib/novainstrumentation/tictac.py @@ -10,7 +10,7 @@ def tic(): _tic = time.time() -def tac(label=""): +def tac(label=''): """@brief This function prints in the screen the difference between the time saved with function tic.py and current time. @@ -23,7 +23,7 @@ def tac(label=""): delta_t = time.time() - _tic - if label != "": - print("%s - %3.4f s" % (label, delta_t)) + if label != '': + print('%s - %3.4f s' % (label, delta_t)) else: - print("%3.4f s" % delta_t) + print('%3.4f s' % delta_t) diff --git a/physiological_data/lib/novainstrumentation/tools.py b/physiological_data/lib/novainstrumentation/tools.py index 8b0c1ec2..37f46a16 100644 --- a/physiological_data/lib/novainstrumentation/tools.py +++ b/physiological_data/lib/novainstrumentation/tools.py @@ -7,7 +7,7 @@ def plotfft(s, fmax, doplot=False): - """This functions computes the fft of a signal, returning the frequency + """ This functions computes the fft of a signal, returning the frequency and their magnitude values. Parameters @@ -30,15 +30,15 @@ def plotfft(s, fmax, doplot=False): fs = abs(np.fft.fft(s)) f = linspace(0, fmax / 2, len(s) / 2) if doplot: - pl.plot(f[1 : len(s) / 2], fs[1 : len(s) / 2]) - return (f[1 : len(s) / 2].copy(), fs[1 : len(s) / 2].copy()) + pl.plot(f[1:len(s) / 2], fs[1:len(s) / 2]) + return (f[1:len(s) / 2].copy(), fs[1:len(s) / 2].copy()) def synthbeats2(duration, meanhr=60, stdhr=1, samplingfreq=250): - # Minimaly based on the parameters from: - # http://physionet.cps.unizar.es/physiotools/ecgsyn/Matlab/ecgsyn.m - # Inputs: duration in seconds - # Returns: signal, peaks + #Minimaly based on the parameters from: + #http://physionet.cps.unizar.es/physiotools/ecgsyn/Matlab/ecgsyn.m + #Inputs: duration in seconds + #Returns: signal, peaks ibi = 60 / float(meanhr) * samplingfreq @@ -50,7 +50,7 @@ def synthbeats2(duration, meanhr=60, stdhr=1, samplingfreq=250): if peaks[-1] >= duration * samplingfreq: peaks = peaks[:-1] - peaks = peaks.astype("int") + peaks = peaks.astype('int') signal = np.zeros(duration * samplingfreq) signal[peaks] = 1.0 @@ -58,11 +58,11 @@ def synthbeats2(duration, meanhr=60, stdhr=1, samplingfreq=250): def synthbeats(duration, meanhr=60, stdhr=1, samplingfreq=250, sinfreq=None): - # Minimaly based on the parameters from: - # http://physionet.cps.unizar.es/physiotools/ecgsyn/Matlab/ecgsyn.m - # If freq exist it will be used to generate a sin instead of using rand - # Inputs: duration in seconds - # Returns: signal, peaks + #Minimaly based on the parameters from: + #http://physionet.cps.unizar.es/physiotools/ecgsyn/Matlab/ecgsyn.m + #If freq exist it will be used to generate a sin instead of using rand + #Inputs: duration in seconds + #Returns: signal, peaks t = np.arange(duration * samplingfreq) / float(samplingfreq) signal = np.zeros(len(t)) @@ -75,26 +75,25 @@ def synthbeats(duration, meanhr=60, stdhr=1, samplingfreq=250, sinfreq=None): npeaks = 1.2 * (duration * meanhr / 60) # add 20% more beats for some cummulative error hr = pl.randn(npeaks) * stdhr + meanhr - peaks = pl.cumsum(60.0 / hr) * samplingfreq - peaks = peaks.astype("int") + peaks = pl.cumsum(60. / hr) * samplingfreq + peaks = peaks.astype('int') peaks = peaks[peaks < t[-1] * samplingfreq] else: hr = meanhr + sin(2 * pi * t * sinfreq) * float(stdhr) - index = int(60.0 / hr[0] * samplingfreq) + index = int(60. / hr[0] * samplingfreq) peaks = [] while index < len(t): peaks += [index] - index += int(60.0 / hr[index] * samplingfreq) + index += int(60. / hr[index] * samplingfreq) signal[peaks] = 1.0 return t, signal, peaks -def load_with_cache( - file_, recache=False, sampling=1, columns=None, temp_dir=".", data_type="int16" -): +def load_with_cache(file_, recache=False, sampling=1, + columns=None, temp_dir='.', data_type='int16'): """@brief This function loads a file from the current directory and saves the cached file to later executions. It's also possible to make a recache or a subsampling of the signal and choose only a few columns of the signal, @@ -114,7 +113,9 @@ def load_with_cache( TODO: receive a file handle """ - cfile = "%s.npy" % file_ + + + cfile = '%s.npy' % file_ if (not path.exists(cfile)) or recache: if columns == None: @@ -132,5 +133,7 @@ def load_data(filename): """ :rtype : numpy matrix """ - data = pandas.read_csv(filename, header=None, delimiter="\t", skiprows=9) + data = pandas.read_csv(filename, header=None, delimiter='\t', skiprows=9) return data.as_matrix() + + diff --git a/physiological_data/lib/novainstrumentation/waves/alignWaves.py b/physiological_data/lib/novainstrumentation/waves/alignWaves.py index 877d45e3..65f6d490 100644 --- a/physiological_data/lib/novainstrumentation/waves/alignWaves.py +++ b/physiological_data/lib/novainstrumentation/waves/alignWaves.py @@ -3,32 +3,32 @@ from ..peaks import peaks -def align(w, signals, mode): - """Align various waves in a specific point. - - Given a mode of alignment, this function computes the specific time point of +def align(w,signals,mode): + """ Align various waves in a specific point. + + Given a mode of alignment, this function computes the specific time point of a wave where all the waves would be aligned. With the difference between the - time point of the reference wave and the time points of all the other waves, - we have the amount of samples the waves will move to align in the specific - point computed. - + time point of the reference wave and the time points of all the other waves, + we have the amount of samples the waves will move to align in the specific + point computed. + Parameters ---------- w: array-like the input signal to use as a reference of alignment (all the other signals will be aligned with this one). signals: array-like or matrix-like - the input signals to align. + the input signals to align. mode: string - the mode used in the alignment, from 'max', 'min', 'peak', 'peak_neg', - 'infMaxAlign' and 'infMinAlign'. + the mode used in the alignment, from 'max', 'min', 'peak', 'peak_neg', + 'infMaxAlign' and 'infMinAlign'. Returns ------- nw: a masked array a new set of aligned signals in a masked array (some cells have NAN values due to the alignment). - + Example ------- >>> align([6,3,4,5,2,2],[10,30,28,26,13,20],'max') @@ -49,45 +49,45 @@ def align(w, signals, mode): mask = [ True False False False False False], fill_value = 1e+20) ] - """ - + """ + nw = [] - - if len(shape(signals)) == 1: + + if len(shape(signals))==1: signals = [signals] - + for i in range(len(signals)): - - if mode == "max": - al = maxAlign(w, signals[i]) - elif mode == "min": - al = minAlign(w, signals[i]) - elif mode == "peak": - al = peakAlign(w, signals[i]) - elif mode == "peak_neg": - al = peakNegAlign(w, signals[i]) - elif mode == "infMaxAlign": - al = infMaxAlign(w, signals[i]) - elif mode == "infMinAlign": - al = infMinAlign(w, signals[i]) - - nw += [moveWave(signals[i], al)] - + + if (mode == 'max'): + al = maxAlign(w,signals[i]) + elif (mode == 'min'): + al = minAlign(w,signals[i]) + elif (mode == 'peak'): + al = peakAlign(w,signals[i]) + elif (mode == 'peak_neg'): + al = peakNegAlign(w,signals[i]) + elif (mode == 'infMaxAlign'): + al = infMaxAlign(w,signals[i]) + elif (mode == 'infMinAlign'): + al = infMinAlign(w,signals[i]) + + nw += [ moveWave(signals[i],al) ] + return nw -def moveWave(w, move): - """Move a signal in time. - +def moveWave(w,move): + """ Move a signal in time. + This function returns a signal created by a shifting in time on the original - signal. - + signal. + Parameters ---------- w: array-like the input signal to move. move: int - the ammount of samples to shift the signal (if <0 the signal moves back, + the ammount of samples to shift the signal (if <0 the signal moves back, if >0 the signal moves forward). Returns @@ -95,34 +95,34 @@ def moveWave(w, move): nw: a masked array a new aligned signal in a masked array (some cells have NAN values due to the alignment). - - + + See also: align() - """ - + """ + nw = ma.masked_all(len(w)) - - if move == 0: - # don't move - nw[0:] = w[0:] - elif move > 0: - # forward - nw[move : len(w)] = w[: len(w) - move] + + if (move == 0): + #don't move + nw[0:] = w[0:] + elif (move>0): + #forward + nw[move:len(w)] = w[:len(w)-move] else: - # backwards - nw[: len(w) - abs(move)] = w[abs(move) : len(w)] - + #backwards + nw[:len(w)-abs(move)] = w[abs(move):len(w)] + return nw -def maxAlign(refw, w): - """Difference between the maximums positions of the signals. - - Given the position in time of each maximum value (argmax) of the signals, +def maxAlign(refw,w): + """ Difference between the maximums positions of the signals. + + Given the position in time of each maximum value (argmax) of the signals, this function returns the difference, in samples, between those two events. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -134,26 +134,26 @@ def maxAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> maxAlign([5,7,3,20,13],[0,5,0,4,7]) -1 - + See also: minAlign(), peakAlign(), peakNegAlign(), infMaxAlign(), infMinAlign() - """ - - return argmax(refw) - argmax(w) + """ + + return argmax(refw)-argmax(w) -def minAlign(refw, w): - """Difference between the minimums positions of the signals. - - Given the position in time of each minimum value (argmin) of the signals, +def minAlign(refw,w): + """ Difference between the minimums positions of the signals. + + Given the position in time of each minimum value (argmin) of the signals, this function returns the difference, in samples, between those two events. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -165,29 +165,29 @@ def minAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> minAlign([5,7,3,20,13],[0,5,-4,4,7]) 0 - + See also: maxAlign(), peakAlign(), peakNegAlign(), infMaxAlign(), infMinAlign() - """ - - return argmin(refw) - argmin(w) - + """ + + return argmin(refw)-argmin(w) -def peakAlign(refw, w): - """Difference between the maximum peak positions of the signals. - This function returns the difference, in samples, between the peaks position - of the signals. If the reference signal has various peaks, the one - chosen is the peak which is closer to the middle of the signal, and if the +def peakAlign(refw,w): + """ Difference between the maximum peak positions of the signals. + + This function returns the difference, in samples, between the peaks position + of the signals. If the reference signal has various peaks, the one + chosen is the peak which is closer to the middle of the signal, and if the other signal has more than one peak also, the chosen is the one closer to - the reference peak signal. + the reference peak signal. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -199,45 +199,43 @@ def peakAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> peakAlign([5,7,3,20,13,5,7],[5,1,8,4,3,10,3]) 1 - + See also: maxAlign(), minAlign(), peakNegAlign(), infMaxAlign(), infMinAlign() - """ - - p_mw = array(peaks(array(refw), min(refw))) - p_w = array(peaks(array(w), min(w))) - - if len(p_mw) > 1: - min_al = argmin( - abs((len(refw) / 2) - p_mw) - ) # to choose the peak closer to the middle of the signal - p_mw = p_mw[min_al] - if list(p_w) == []: + """ + + p_mw = array ( peaks(array(refw),min(refw)) ) + p_w = array ( peaks(array(w),min(w)) ) + + + if (len(p_mw)>1): + min_al = argmin(abs( (len(refw)/2) - p_mw)) #to choose the peak closer to the middle of the signal + p_mw=p_mw[min_al] + if (list(p_w) == [] ): p_w = p_mw - elif len(p_w) > 1: - min_al = argmin( - abs(p_w - p_mw) - ) # to choose the peak closer to the peak of the reference signal - p_w = p_w[min_al] - - return int(array(p_mw - p_w)) - - -def peakNegAlign(refw, w): - """Difference between the minimum peak positions of the signals. - - This function returns the difference, in samples, between the minimum peaks - position of the signals. If the reference signal has various peaks, the one - chosen is the peak which is closer to the middle of the signal, and if the + elif (len(p_w)>1): + min_al = argmin(abs(p_w - p_mw)) #to choose the peak closer to the peak of the reference signal + p_w=p_w[min_al] + + + return int(array(p_mw-p_w)) + + +def peakNegAlign(refw,w): + """ Difference between the minimum peak positions of the signals. + + This function returns the difference, in samples, between the minimum peaks + position of the signals. If the reference signal has various peaks, the one + chosen is the peak which is closer to the middle of the signal, and if the other signal has more than one peak also, the chosen is the one closer to - the reference peak signal. + the reference peak signal. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -249,43 +247,39 @@ def peakNegAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> peakNegAlign([5,7,3,20,13,5,7],[5,1,8,4,3,10,3]) 1 - + See also: maxAlign(), minAlign(), peakAlign(), infMaxAlign(), infMinAlign() - """ - - p_mw = array(peaks(-array(refw), min(-array(refw)))) - p_w = array(peaks(-array(w), min(-array(w)))) - - if len(p_mw) > 1: - min_al = argmin( - abs((len(refw) / 2) - p_mw) - ) # to choose the peak closer to the middle of the signal - p_mw = p_mw[min_al] - - if list(p_w) == []: + """ + + p_mw = array ( peaks(-array(refw),min(-array(refw)) ) ) + p_w = array ( peaks(-array(w),min(-array(w)) ) ) + + if (len(p_mw)>1): + min_al = argmin(abs( (len(refw)/2) - p_mw)) #to choose the peak closer to the middle of the signal + p_mw=p_mw[min_al] + + if (list(p_w) == [] ): p_w = p_mw - elif len(p_w) > 1: - min_al = argmin( - abs(p_w - p_mw) - ) # to choose the peak closer to the peak of the reference signal - p_w = p_w[min_al] - - return int(array(p_mw - p_w)) - + elif (len(p_w)>1): + min_al = argmin(abs(p_w - p_mw)) #to choose the peak closer to the peak of the reference signal + p_w=p_w[min_al] + + return int(array(p_mw-p_w)) -def infMaxAlign(refw, w): - """Difference between the maximum peak positions of the signal's derivative - This function returns the difference, in samples, between the maximum peaks +def infMaxAlign(refw,w): + """ Difference between the maximum peak positions of the signal's derivative + + This function returns the difference, in samples, between the maximum peaks position of the derivative signal. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -297,27 +291,27 @@ def infMaxAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> infMaxAlign([97,87,45,65,34],[57,10,84,93,32]) 1 - - See also: maxAlign(), minAlign(), peakAlign(), peakNegAlign(), + + See also: maxAlign(), minAlign(), peakAlign(), peakNegAlign(), infMinAlign() - """ - - return argmax(diff(refw)) - argmax(diff(w)) + """ + + return argmax(diff(refw))-argmax(diff(w)) -def infMinAlign(refw, w): - """Difference between the minimum peak positions of the signal's derivative - - This function returns the difference, in samples, between the minimum peaks +def infMinAlign(refw,w): + """ Difference between the minimum peak positions of the signal's derivative + + This function returns the difference, in samples, between the minimum peaks position of the derivative signal. The first signal introduced is the reference signal. - - + + Parameters ---------- refw: array-like @@ -329,14 +323,15 @@ def infMinAlign(refw, w): ------- al: int the difference between the two events position - + Example ------- >>> infMinAlign([67,87,45,65,34],[57,63,84,12,32]) -1 - - See also: maxAlign(), minAlign(), peakAlign(), peakNegAlign(), + + See also: maxAlign(), minAlign(), peakAlign(), peakNegAlign(), infMaxAlign() - """ + """ + + return argmin(diff(refw))-argmin(diff(w)) - return argmin(diff(refw)) - argmin(diff(w)) diff --git a/physiological_data/lib/novainstrumentation/waves/computemeanwave.py b/physiological_data/lib/novainstrumentation/waves/computemeanwave.py index c1da1bf5..3c217ca8 100644 --- a/physiological_data/lib/novainstrumentation/waves/computemeanwave.py +++ b/physiological_data/lib/novainstrumentation/waves/computemeanwave.py @@ -1,13 +1,13 @@ -def computemeanwave(signal, events, fdist, lmin=0, lmax=0): - # Output: Array with mean wave, standard deviation and distance to the mean wave - # acordding to fdist +def computemeanwave(signal, events, fdist, lmin=0,lmax=0): +# Output: Array with mean wave, standard deviation and distance to the mean wave +# acordding to fdist - if (lmin == 0) & (lmax == 0): - lmax = mean(diff(events)) / 2 - lmin = -lmax - w = waves(signal, events, lmin, lmax) - w_ = meanwave(w) - d = wavedistance(w_, w, fdist) - ws = stdwave(w) - - return (w_, ws, d) + if (lmin==0) & (lmax==0): + lmax=mean(diff(events))/2 + lmin=-lmax + w=waves(signal, events, lmin, lmax) + w_=meanwave(w) + d=wavedistance(w_,w,fdist) + ws=stdwave(w) + + return (w_, ws, d) \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/waves/frequencyAlignDistance.py b/physiological_data/lib/novainstrumentation/waves/frequencyAlignDistance.py index 0647ce9f..89e3f41e 100644 --- a/physiological_data/lib/novainstrumentation/waves/frequencyAlignDistance.py +++ b/physiological_data/lib/novainstrumentation/waves/frequencyAlignDistance.py @@ -3,40 +3,39 @@ from scipy import interpolate -def frequencyAlignDistance(f1, w1, f2, w2): - - w1 = smooth(w1) - w2 = smooth(w2) - - # return msedistance(w2[:100],w1[:100]) - - m1 = argmax(w1) - m2 = argmax(w2) - - f1 = f1 / float(f1[m1]) - f2 = f2 / float(f2[m2]) - - w1 = w1 / w1[m1] - w2 = w2 / w2[m2] - - # plot(f1,w1) - # plot(f2,w2) - - nf = f1[3:103] - - f_interp_w2 = interpolate.interp1d(f2, w2) - - iw2 = f_interp_w2(nf) - - plot(nf, iw2) - plot(nf, w1[3:103]) - - s2 = iw2 - s1 = w1[3:103] - - r = np.sqrt(np.sum(((s1 - s2) ** 2) / s1)) - +def frequencyAlignDistance(f1,w1,f2,w2): + + w1=smooth(w1) + w2=smooth(w2) + +# return msedistance(w2[:100],w1[:100]) + + m1=argmax(w1) + m2=argmax(w2) + + f1=f1/float(f1[m1]) + f2=f2/float(f2[m2]) + + w1=w1/w1[m1] + w2=w2/w2[m2] + +# plot(f1,w1) +# plot(f2,w2) + + nf=f1[3:103] + + f_interp_w2=interpolate.interp1d(f2,w2) + + iw2=f_interp_w2(nf) + + plot(nf,iw2) + plot(nf,w1[3:103]) + + s2=iw2 + s1=w1[3:103] + + r=np.sqrt(np.sum(((s1-s2)**2)/s1 )) + return r - - -# return msedistance(iw2,w1[3:103]) + +# return msedistance(iw2,w1[3:103]) \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/waves/getarray.py b/physiological_data/lib/novainstrumentation/waves/getarray.py index da304b07..f8708d61 100644 --- a/physiological_data/lib/novainstrumentation/waves/getarray.py +++ b/physiological_data/lib/novainstrumentation/waves/getarray.py @@ -1,8 +1,8 @@ def getarray(d): - # select the array variable from a matlab structure - - p = ["_" not in k for k in d.keys()] - +#select the array variable from a matlab structure + + p=['_' not in k for k in d.keys()] + for i in range(len(p)): if p[i]: - return d[d.keys()[i]] + return d[d.keys()[i]] \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/waves/meanwave.py b/physiological_data/lib/novainstrumentation/waves/meanwave.py index 4cfa3912..d405745f 100644 --- a/physiological_data/lib/novainstrumentation/waves/meanwave.py +++ b/physiological_data/lib/novainstrumentation/waves/meanwave.py @@ -1,13 +1,12 @@ from numpy import mean - def meanwave(signals): - """This function computes the meanwave of various signals. - - Given a set of signals, with the same number of samples, this function + """ This function computes the meanwave of various signals. + + Given a set of signals, with the same number of samples, this function returns an array representative of the meanwave of those signals - which is - a wave computed with the mean values of each signal's samples. - + a wave computed with the mean values of each signal's samples. + Parameters ---------- signals: matrix-like @@ -16,7 +15,7 @@ def meanwave(signals): Returns ------- mw: array-like - the resulted meanwave + the resulted meanwave """ - - return mean(signals, 0) + + return mean(signals,0) diff --git a/physiological_data/lib/novainstrumentation/waves/plotheatmap.py b/physiological_data/lib/novainstrumentation/waves/plotheatmap.py index f8607f84..b5f95531 100644 --- a/physiological_data/lib/novainstrumentation/waves/plotheatmap.py +++ b/physiological_data/lib/novainstrumentation/waves/plotheatmap.py @@ -4,16 +4,16 @@ from novainstrumentation.waves.waves import waves -def plotheatmap(signal, events, lmin=0, lmax=0, dt=0.01, color="r"): - - w = waves(signal, events, lmin, lmax) - w_ = meanwave(w) - # ws_=stdwave(w) - t_ = (arange(len(w_)) + lmin) * dt +def plotheatmap( signal, events,lmin=0,lmax=0,dt=0.01,color='r'): + w=waves(signal, events, lmin, lmax) + w_=meanwave(w) + # ws_=stdwave(w) + t_=(arange(len(w_))+lmin)*dt + for iw in w: - plot(t_, iw, lw=3, alpha=0.15, color=color) + plot(t_,iw,lw=3,alpha=.15,color=color) - plot(t_, w_, lw=2, color="k") - # plot(t_,w_+ws_,lw=0.5,color='k',ls='--') - # plot(t_,w_-ws_,lw=0.5,color='k',ls='--') + plot(t_,w_,lw=2,color='k') + #plot(t_,w_+ws_,lw=0.5,color='k',ls='--') + #plot(t_,w_-ws_,lw=0.5,color='k',ls='--') diff --git a/physiological_data/lib/novainstrumentation/waves/plotheatmap2.py b/physiological_data/lib/novainstrumentation/waves/plotheatmap2.py index 9cced7e2..117f65d5 100644 --- a/physiological_data/lib/novainstrumentation/waves/plotheatmap2.py +++ b/physiological_data/lib/novainstrumentation/waves/plotheatmap2.py @@ -1,19 +1,21 @@ -def plotheatmap2(s, valmin, valmax, points=100): +def plotheatmap2(s,valmin,valmax,points=100): + import numpy as np import matplotlib.cm as cm import matplotlib.mlab as mlab import matplotlib.pylab as plt - x = range(1, len(s)) + x = range(1,len(s)) y = np.linspace(valmin, valmax, points) X, Y = np.meshgrid(x, y) - Z = mlab.bivariate_normal(X, Y) * 0.0 + Z = mlab.bivariate_normal(X, Y)*0.0 for i in range(len(s)): Z1 = mlab.bivariate_normal(X, Y, 3.0, 3.0, i, s[i]) - Z = Z1 + Z + Z=Z1+Z - im = plt.imshow(Z, interpolation="bilinear", cmap=cm.gray, origin="lower") - # extent=[-3,3,-3,3]) - return Z + im = plt.imshow(Z, interpolation='bilinear', cmap=cm.gray, + origin='lower') +#extent=[-3,3,-3,3]) + return Z \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/waves/stdwave.py b/physiological_data/lib/novainstrumentation/waves/stdwave.py index c920d886..2714d974 100644 --- a/physiological_data/lib/novainstrumentation/waves/stdwave.py +++ b/physiological_data/lib/novainstrumentation/waves/stdwave.py @@ -1,12 +1,11 @@ from numpy import std - def stdwave(signal): - """This function computes the standard deviation error wave of various + """ This function computes the standard deviation error wave of various signals. - - Given a set of signals, with the same number of samples, this function - returns an array representative of the error wave of those signals - which + + Given a set of signals, with the same number of samples, this function + returns an array representative of the error wave of those signals - which is a wave computed with the standard deviation error values of each signal's samples. @@ -20,6 +19,6 @@ def stdwave(signal): Returns ------- stdw: array-like - the resulting error wave. - """ - return std(signal, 0) + the resulting error wave. + """ + return std(signal,0) diff --git a/physiological_data/lib/novainstrumentation/waves/sumvolve.py b/physiological_data/lib/novainstrumentation/waves/sumvolve.py index 225f98f9..73d8effc 100644 --- a/physiological_data/lib/novainstrumentation/waves/sumvolve.py +++ b/physiological_data/lib/novainstrumentation/waves/sumvolve.py @@ -1,7 +1,7 @@ -def sumvolve(x, window): - lw = len(window) - res = zeros(len(x) - lw, "d") - for i in range(len(x) - lw): - res[i] = sum(abs(x[i : i + lw] - window)) / float(lw) - # res[i]=sum(abs(x[i:i+lw]*window)) +def sumvolve(x,window): + lw=len(window) + res=zeros(len(x)-lw,'d') + for i in range(len(x)-lw): + res[i]=sum(abs(x[i:i+lw]-window))/float(lw) + #res[i]=sum(abs(x[i:i+lw]*window)) return res diff --git a/physiological_data/lib/novainstrumentation/waves/wavedistance.py b/physiological_data/lib/novainstrumentation/waves/wavedistance.py index 17433c3c..374ae886 100644 --- a/physiological_data/lib/novainstrumentation/waves/wavedistance.py +++ b/physiological_data/lib/novainstrumentation/waves/wavedistance.py @@ -1,6 +1,5 @@ import numpy as np - -def wavedistance(meanwave, waves, fdistance): - - return np.array([fdistance(wave, meanwave) for wave in waves]) +def wavedistance(meanwave,waves,fdistance): + + return np.array([fdistance(wave,meanwave) for wave in waves]) \ No newline at end of file diff --git a/physiological_data/lib/novainstrumentation/waves/waves.py b/physiological_data/lib/novainstrumentation/waves/waves.py index 97edbe43..50420736 100644 --- a/physiological_data/lib/novainstrumentation/waves/waves.py +++ b/physiological_data/lib/novainstrumentation/waves/waves.py @@ -1,19 +1,16 @@ import numpy as np -# testar -# escrever comentários: -# dizer o que a função faz -# referir valores negativos faz com que a função segmente de valores anteriores ao evento até valores superiores +#testar +#escrever comentários: + # dizer o que a função faz + # referir valores negativos faz com que a função segmente de valores anteriores ao evento até valores superiores # sugestões de melhoria: -# testar se lowerBound < upperBound -# aviso de segmentos ignorados quando [(upperBound - lowerBound) == len(i)] == False -# bloquear o inicio e o fim da funcao - + # testar se lowerBound < upperBound + # aviso de segmentos ignorados quando [(upperBound - lowerBound) == len(i)] == False + # bloquear o inicio e o fim da funcao def waves(signal, events, lowerBound, upperBound): - signal = [ - signal[(center + lowerBound) : (center + upperBound)] for center in events - ] + signal = [signal[(center + lowerBound):(center + upperBound)] for center in events] x = np.array(list(filter(lambda i: (upperBound - lowerBound) == len(i), signal))) - return x + return x \ No newline at end of file diff --git a/physiological_data/lib/peak_detector.py b/physiological_data/lib/peak_detector.py index cde7bc58..921cb9f5 100644 --- a/physiological_data/lib/peak_detector.py +++ b/physiological_data/lib/peak_detector.py @@ -2,24 +2,21 @@ from scipy.signal import decimate try: - from physiological_data.lib.signal_processing import ( - peakDetection, - pan_tompkins_filter, - ) + from physiological_data.lib.signal_processing import peakDetection, pan_tompkins_filter except (ImportError, ModuleNotFoundError): from physiological_data.lib.signal_processing import peakDetection def peak_detector_RT(signal, fs, step=1, win_size=1.5, x_thr=0.3): if fs > 50: - signal = decimate(signal, int(fs / 50)) - fs /= int(fs / 50) + signal = decimate(signal, int(fs/50)) + fs /= int(fs/50) win_size *= fs x_thr *= fs win = [] - pro_th = np.inf # Update the value after the first 3 seconds + pro_th = np.inf # Update the value after the first 3 seconds - has_valley = False # A valley can only be found after a peak + has_valley = False # A valley can only be found after a peak valleys, peaks = [[], []], [[], []] peak_index, valley_index = 0, 0 @@ -45,10 +42,10 @@ def peak_detector_RT(signal, fs, step=1, win_size=1.5, x_thr=0.3): win_aux = pan_tompkins_filter(win, fs) pro = 0.3 * np.ptp(win_aux) - sup_thr = min(win_aux) + 0.5 * np.ptp(win_aux) - inf_thr = min(win_aux) + 0.1 * np.ptp(win_aux) + sup_thr = min(win_aux) + .5 * np.ptp(win_aux) + inf_thr = min(win_aux) + .1 * np.ptp(win_aux) - value = win_aux[len(win) // 2] + value = win_aux[len(win)//2] ###################################################################################### ################################## VALLEY DETECTION ################################## @@ -56,35 +53,24 @@ def peak_detector_RT(signal, fs, step=1, win_size=1.5, x_thr=0.3): valley_condition = True if len(valleys[1]) > 0: valley_condition = index > valleys[1][-1] + x_thr - + # print(valley_value, value, inf_thr, valley_condition, np.abs(peak_value - value)) - if ( - value < valley_value - and value < inf_thr - and valley_condition - and np.abs(peak_value - value) > pro - ): + if value < valley_value and value < inf_thr and valley_condition and np.abs(peak_value - value) > pro: valley_value = value - valley_index = int(index + len(win) // 2) + valley_index = int(index + len(win)//2) has_valley = True - + if value > sup_thr and old_has_valley: has_valley = False - + # Peak Detection between two valleys - if ( - not has_valley - and value > peak_value - and np.abs(value - valley_value) > pro - ): + if not has_valley and value > peak_value and np.abs(value - valley_value) > pro: peak_value = value - peak_index = int(index + len(win) // 2) - + peak_index = int(index + len(win)//2) + valley_condition = True if len(peaks[1]) > 0 and len(valleys[1]) > 0: - valley_condition = ( - peaks[1][-1] < valley_index and peaks[1][-1] > valleys[1][-1] - ) + valley_condition = peaks[1][-1] < valley_index and peaks[1][-1] > valleys[1][-1] if old_has_valley != has_valley and not has_valley and valley_condition: valleys[0].append(valley_value) valleys[1].append(valley_index) @@ -103,59 +89,54 @@ def peak_detector_RT(signal, fs, step=1, win_size=1.5, x_thr=0.3): return peaks, valleys -def peak_detector(signal, fs, window=3, overlap=0.3, x_thr=0.3, pro=0.5): +def peak_detector(signal, fs, window=3, overlap=.3, x_thr=.3, pro=0.5): peaks, peaksPos, valleys, valleysPos = [], [], [], [] window *= fs - overlap = int(fs * overlap) - + overlap = int(fs*overlap) + for i in range(overlap, len(signal), window): - segment = signal[i - overlap : i + window] - - peaks, peakPos_det, _, valleyPos_det = peakDetection( - segment, x_threshold=int(overlap * fs), proeminence=pro * np.ptp(segment) - ) + segment = signal[i - overlap:i+window] + + peaks, peakPos_det, _, valleyPos_det = peakDetection(segment, x_threshold=int(overlap * fs), proeminence=pro*np.ptp(segment)) if len(peakPos_det) > 1: peakPos_det, valleyPos_det = peakPos_det[:-1], valleyPos_det[:-1] - + for p in range(len(peakPos_det)): peaksPos += [peakPos_det[p] + i - overlap] valleysPos += [valleyPos_det[p] + i - overlap] overlap = len(segment) - peakPos_det[-1] + 1 - + return (peaks, peaksPos), (valleys, valleysPos) -if __name__ == "__main__": + +if __name__ == '__main__': from biosignals import Devices import matplotlib.pyplot as plt - device = Devices(r"..\..\acquisitions\Acquisitions\03_11_2020") - data_ecg = device.getSensorsData(["ECG"]) - ecg_signal = device.convertAndfilterSignal( - data_ecg["data"][:, 1], "ECG", device.fs, device.resolution - ) + device = Devices(r'..\..\acquisitions\Acquisitions\03_11_2020') + data_ecg = device.getSensorsData(['ECG']) + ecg_signal = device.convertAndfilterSignal(data_ecg['data'][:, 1], 'ECG', device.fs, device.resolution) - peaks, valleys = peak_detector_RT(ecg_signal, device.fs, x_thr=0.1) + peaks, valleys = peak_detector_RT(ecg_signal, device.fs, x_thr=.1) plt.figure() - plt.plot(decimate(ecg_signal, int(device.fs / 50))) + plt.plot(decimate(ecg_signal, int(device.fs/50))) # plt.plot(ecg_signal) # plt.scatter(peaks[1], peaks[0], c='r') - plt.vlines(peaks[1], -0.5, 0.75, "r") + plt.vlines(peaks[1], -.5, .75, 'r') plt.show() - data_resp = device.getSensorsData(["RESPIRATION"]) - resp_signal = device.convertAndfilterSignal( - data_resp["data"][:, 1], "RESPIRATION", device.fs, device.resolution - ) + data_resp = device.getSensorsData(['RESPIRATION']) + resp_signal = device.convertAndfilterSignal(data_resp['data'][:, 1], 'RESPIRATION', device.fs, device.resolution) peaks, valleys = peak_detector_RT(resp_signal, device.fs, win_size=20) plt.figure() - plt.plot(decimate(resp_signal, int(device.fs / 50))) + plt.plot(decimate(resp_signal, int(device.fs/50))) # plt.plot(ecg_signal) # plt.scatter(peaks[1], peaks[0], c='r') - plt.vlines(peaks[1], -0.5, 0.75, "r") + plt.vlines(peaks[1], -.5, .75, 'r') plt.show() diff --git a/physiological_data/lib/psychopy.py b/physiological_data/lib/psychopy.py index 81377975..ce953939 100644 --- a/physiological_data/lib/psychopy.py +++ b/physiological_data/lib/psychopy.py @@ -5,11 +5,11 @@ class Psycho: def __init__(self, file_path, task=None): self.file_path = file_path - if task is "N-Back": + if task is 'N-Back': self.data = self.getNBack() - elif task is "Subtraction": + elif task is 'Subtraction': self.data = self.getSubtraction() - elif task is "Baseline": + elif task is 'Baseline': self.data = self.getBaseline() else: self.data = self.getPsychoData() @@ -19,54 +19,46 @@ def getPsychoData(self): def getNBack(self): data = self.getPsychoData() - indexes = np.where(np.array(data["Activity"]) == "N-back")[0] + indexes = np.where(np.array(data['Activity']) == 'N-back')[0] return data.loc[np.append(indexes, indexes[-1] + 1)] def getSubtraction(self): data = self.getPsychoData() - indexes = np.sort( - np.concatenate( - [ - np.where(np.array(data["Activity"]) == "o")[0], - np.where(np.array(data["Activity"]) == "sub_baseline")[0], - ] - ) - ) - return data.loc[np.append(indexes, indexes[-1] + 1)] + indexes = np.sort(np.concatenate([np.where(np.array(data['Activity']) == 'o')[0], + np.where(np.array(data['Activity']) == 'sub_baseline')[0]])) + return data.loc[np.append(indexes, indexes[-1] +1)] def getBaseline(self): data = self.getPsychoData() - indexes = np.sort( - np.concatenate( - [ - np.where(np.array(data["Activity"]) == "Begin of Baseline")[0], - np.where(np.array(data["Activity"]) == "End of Baseline")[0], - ] - ) - ) + indexes = np.sort(np.concatenate([np.where(np.array(data['Activity']) == 'Begin of Baseline')[0], + np.where(np.array(data['Activity']) == 'End of Baseline')[0]])) return data.loc[np.append(indexes, indexes[-1])] def getTimestamps(self): - return np.array(self.data["Timestamp"]) + return np.array(self.data['Timestamp']) def getActivity(self): - return np.array(self.data["Activity"]) + return np.array(self.data['Activity']) def getResult(self): - return np.array(self.data["Result"]) + return np.array(self.data['Result']) def getDifficulty(self): - return np.array(self.data["Difficulty"]) + return np.array(self.data['Difficulty']) def getTimeAnswer(self): - return np.array(self.data["Time to Answer"]) + return np.array(self.data['Time to Answer']) def getKeyAnswer(self): - return np.array(self.data["Key Answer"]) + return np.array(self.data['Key Answer']) + + -if __name__ == "__main__": - psycho = Psycho(r"..\..\acquisitions\Acquisitions\03_11_2020\results_3.csv") +if __name__ == '__main__': + psycho = Psycho(r'..\..\acquisitions\Acquisitions\03_11_2020\results_3.csv') print(psycho.getNBack()) print(psycho.getSubtraction()) print(psycho.getBaseline()) + + diff --git a/physiological_data/lib/respRT.py b/physiological_data/lib/respRT.py index 7d74ca4f..ca9c5a22 100644 --- a/physiological_data/lib/respRT.py +++ b/physiological_data/lib/respRT.py @@ -1,10 +1,9 @@ + import numpy as np from scipy.signal import decimate -def peak_detector_Resp( - signal, fs, window=20, x_threshold=0.3, prominence_threshold=0.1 -): +def peak_detector_Resp(signal, fs, window=20, x_threshold=.3, prominence_threshold=.1): if fs > 50: signal = decimate(signal, int(fs / 50)) fs /= int(fs / 50) @@ -25,7 +24,7 @@ def peak_detector_Resp( index = 0 for i in range(0, len(signal)): - print(f"{i * 100 / len(signal):.2f}%", end="\r") + print(f"{i * 100 / len(signal):.2f}%", end='\r') old_has_peak = has_peak if len(win) >= win_size: @@ -55,11 +54,7 @@ def peak_detector_Resp( valley_value = value if len(peaks[1]) > 0: - if ( - old_has_peak != has_peak - and not has_peak - and peak_value - peaks[1][-1] > prominence - ): + if old_has_peak != has_peak and not has_peak and peak_value - peaks[1][-1] > prominence: peaks[0].append(peak_index) peak_value = -np.inf @@ -72,15 +67,13 @@ def peak_detector_Resp( return peaks, signal[peaks[0]] -if __name__ == "__main__": +if __name__ == '__main__': from biosignals import Devices import matplotlib.pyplot as plt - device = Devices(r"..\..\acquisitions\Acquisitions\03_11_2020") - data = device.getSensorsData(["RESPIRATION"]) - signal = device.convertAndfilterSignal( - data["data"][:, 1], "RESPIRATION", device.fs, device.resolution - ) + device = Devices(r'..\..\acquisitions\Acquisitions\03_11_2020') + data = device.getSensorsData(['RESPIRATION']) + signal = device.convertAndfilterSignal(data['data'][:, 1], 'RESPIRATION', device.fs, device.resolution) peaks, _ = peak_detector_Resp(signal, device.fs) @@ -88,5 +81,5 @@ def peak_detector_Resp( plt.figure() plt.plot(signal, zorder=-1) - plt.scatter(peaks[0], signal[peaks[0]], c="r", marker="x", zorder=1) + plt.scatter(peaks[0], signal[peaks[0]], c='r', marker='x', zorder=1) plt.show() diff --git a/physiological_data/lib/save_data.py b/physiological_data/lib/save_data.py index aa7a03a7..dc9c7451 100644 --- a/physiological_data/lib/save_data.py +++ b/physiological_data/lib/save_data.py @@ -1,6 +1,5 @@ from acquisition import * - def concatenateSegments(segments): out = np.empty((segments.shape[0], segments.shape[1] * segments.shape[2])) for i in range(segments.shape[0]): @@ -8,35 +7,29 @@ def concatenateSegments(segments): return out -folders = [ - "05_11_2020_1", - "05_11_2020_2", - "10_11_2020_2", - "27_11_2020", - "11_12_2020_1", - "11_12_2020_2", -] +folders = ['05_11_2020_1', + '05_11_2020_2', + '10_11_2020_2', + '27_11_2020', + '11_12_2020_1', + '11_12_2020_2' + ] # folders = ['05_11_2020_1',] segments, labels, group = [], [], [] for i, folder in enumerate(folders): - path = os.path.join(r"..\acquisitions\Acquisitions", folder) - + path = os.path.join(r'..\acquisitions\Acquisitions', folder) + acquisition = Acquisition(path) - - timestamps, intialLabels = ( - acquisition.getPsychoTimestamps(), - acquisition.getPsychoActivity(), - ) - + + timestamps, intialLabels = acquisition.getPsychoTimestamps(), acquisition.getPsychoActivity() + acquisition.deleteArtifacts() acquisition.getBiosignalsSensors(ALL, rightPos=True) acquisition.convertSensors() - acquisition.segmentWindowingBiosignals( - acquisition.signal, timestamps, intialLabels, timeWindow=0.128 - ) + acquisition.segmentWindowingBiosignals(acquisition.signal, timestamps, intialLabels, timeWindow=.128) if len(segments) == 0: segments = concatenateSegments(acquisition.segmentedBiosignals) @@ -44,13 +37,9 @@ def concatenateSegments(segments): labels = np.array(acquisition.labels) group = np.array([folder] * segments.shape[0]) else: - segments = np.vstack( - [segments, concatenateSegments(acquisition.segmentedBiosignals)] - ) + segments = np.vstack([segments, concatenateSegments(acquisition.segmentedBiosignals)]) labels = np.concatenate([labels, acquisition.labels]) - group = np.concatenate( - [group, [folder] * acquisition.segmentedBiosignals.shape[0]] - ) + group = np.concatenate([group, [folder] * acquisition.segmentedBiosignals.shape[0]]) dataset = np.hstack([segments, labels.reshape(-1, 1), group.reshape(-1, 1)]) -np.save("dataset_0_128s", dataset) +np.save('dataset_0_128s', dataset) diff --git a/physiological_data/lib/signal_processing.py b/physiological_data/lib/signal_processing.py index b94515c0..cd5a8255 100644 --- a/physiological_data/lib/signal_processing.py +++ b/physiological_data/lib/signal_processing.py @@ -1,38 +1,33 @@ import numpy as np from biosignalsnotebooks import bsnb - def almostEqual(v1, v2, delta): if isinstance(v1, int): return abs(v1 - v2) < delta else: return all(abs(v1 - v2) < delta) - def pan_tompkins_filter(data, fs): # Pan & Tompkins signal conditioning: https://www.robots.ox.ac.uk/~gari/teaching/cdt/A3/readings/ECG/Pan+Tompkins.pdf # (used the normal differentiation, as the 4-point differentiation gave equal results) - red_diff = differentation(data) # differentiate data - red_sqr = square(red_diff) # square it - red_sqr_int_DO = integration_DO(red_sqr, fs) # and integrate it + red_diff = differentation(data) # differentiate data + red_sqr = square(red_diff) # square it + red_sqr_int_DO = integration_DO(red_sqr, fs) # and integrate it return red_sqr_int_DO - def getValuesArray(array, index_s, index_e): index_s, index_e = int(index_s), int(index_e) - new_array = np.zeros((index_e - index_s,)) + new_array = np.zeros((index_e - index_s, )) for i in range(index_s, index_e): new_array[i - index_s] = array[i] return new_array - def differentation(signal): diff_signal = np.empty(len(signal) - 1) for i in range(1, len(signal)): - diff_signal[i - 1] = signal[i] - signal[i - 1] + diff_signal[i-1] = signal[i] - signal[i-1] return diff_signal - def mean(data): mean_value = 0 for i in range(len(data)): @@ -40,7 +35,6 @@ def mean(data): mean_value = mean_value / len(data) return mean_value - def std(data): mean_value = mean(data) std_value = 0 @@ -56,83 +50,73 @@ def std(data): def square(data): square_data = np.zeros(len(data)) for i in range(len(data)): - square_data[i] = data[i] * data[i] + square_data[i] = data[i]*data[i] return square_data - def smooth_(signal, n): smoothen = np.zeros(len(signal)) - aux_signal = np.zeros(2 * n + len(signal)) + aux_signal = np.zeros(2*n + len(signal)) for i in range(n): # mirror the start of signal - aux_signal[i] = signal[n - i] + aux_signal[i] = signal[n-i] # mirror the end of signal - aux_signal[len(signal) + n + i] = signal[len(signal) - i - 2] + aux_signal[len(signal)+n+i] = signal[len(signal)-i-2] for i in range(len(signal)): - aux_signal[i + n] = signal[i] - for i in range(n, len(aux_signal) - n): - smoothen[i - n] = mean(getValuesArray(aux_signal, i - n // 2, i + n // 2)) + aux_signal[i+n] = signal[i] + for i in range(n, len(aux_signal)-n): + smoothen[i-n] = mean(getValuesArray(aux_signal, i-n//2, i+n//2)) return smoothen - def smooth(signal, n): smoothen = np.zeros(len(signal)) - aux_signal = np.zeros(2 * n + len(signal)) + aux_signal = np.zeros(2*n + len(signal)) for i in range(len(signal)): if i < n: # mirror the start of signal - aux_signal[i] = signal[n - i] + aux_signal[i] = signal[n-i] # mirror the end of signal - aux_signal[len(signal) + n + i] = signal[len(signal) - i - 2] - aux_signal[i + n] = signal[i] - for i in range(n, len(aux_signal) - n): - smoothen[i - n] = mean(getValuesArray(aux_signal, i - n // 2, i + n // 2)) + aux_signal[len(signal)+n+i] = signal[len(signal)-i-2] + aux_signal[i+n] = signal[i] + for i in range(n, len(aux_signal)-n): + smoothen[i-n] = mean(getValuesArray(aux_signal, i-n//2, i+n//2)) return smoothen +def bsnbSmooth(input_signal, window_len=10, window='hanning'): + sig = np.r_[input_signal[window_len:0:-1], input_signal, input_signal[-2:-window_len-2:-1]] + win = np.ones(window_len, 'd') + sig_conv = np.convolve(win / win.sum(), sig, mode='same') -def bsnbSmooth(input_signal, window_len=10, window="hanning"): - sig = np.r_[ - input_signal[window_len:0:-1], - input_signal, - input_signal[-2 : -window_len - 2 : -1], - ] - win = np.ones(window_len, "d") - sig_conv = np.convolve(win / win.sum(), sig, mode="same") - - return sig_conv[window_len:-window_len] + return sig_conv[window_len: -window_len] # return sig - def cumsum(data): value = data for i in range(1, len(data)): - value[i] = value[i - 1] + data[i] + value[i] = value[i-1] + data[i] return value - def integrate(data, rate): - wind_size = int(0.15 * rate) - int_ecg = np.zeros(len(data) - wind_size) + wind_size = int(0.15*rate) + int_ecg = np.zeros(len(data)-wind_size) cum_sum = cumsum(data) for i in range(len(int_ecg)): if i >= wind_size: - int_ecg[i] = (cum_sum[i] - cum_sum[i - wind_size]) / wind_size + int_ecg[i] = (cum_sum[i] - cum_sum[i-wind_size]) / wind_size else: - int_ecg[i] = cum_sum[i] / (i + 1) + int_ecg[i] = cum_sum[i] / (i+1) return int_ecg - def integration_DO(data, fs): # Moving window integration. N is the number of samples in the width of the integration window - out = [] - points = int(fs * 0.15) - for i in range(0, len(data) - points): - temp = 0 + out = [ ] + points = int(fs*0.15) + for i in range(0, len(data)-points): + temp = 0; for j in range(0, points): - temp += data[i + j] - - out.append(temp / points) - return out + temp += data[i+j]; + + out.append(temp/points); + return out; def peakDetection_(data, x_threshold=0, y_threshold=-np.inf, delta=0, proeminence=0): @@ -145,20 +129,12 @@ def peakDetection_(data, x_threshold=0, y_threshold=-np.inf, delta=0, proeminenc if prev_value > 0 and diff_data[i] < 0 and data[i] > y_threshold or i == 0: if x_threshold > 0: if i + x_threshold < len(diff_data): - pos, peak = higherPeak( - getValuesArray(data, i, i + x_threshold), - getValuesArray(diff_data, i, i + x_threshold), - y_threshold, - ) + pos, peak = higherPeak(getValuesArray(data, i, i+x_threshold), getValuesArray(diff_data, i, i+x_threshold), y_threshold) else: - pos, peak = higherPeak( - getValuesArray(data, i, len(diff_data) - 1), - getValuesArray(diff_data, i, len(diff_data) - 1), - y_threshold, - ) + pos, peak = higherPeak(getValuesArray(data, i, len(diff_data)-1), getValuesArray(diff_data, i, len(diff_data)-1), y_threshold) peaks.append(peak) - peaks_pos.append(i + pos) + peaks_pos.append(i+pos) i += pos + x_threshold else: peaks.append(data[i]) @@ -191,11 +167,11 @@ def peakDetection(data, x_threshold=0, y_threshold=-np.inf, delta=0, proeminence def median(array): array = np.sort(array) length = int(len(array)) - if length % 2 == 0: - sumOfMiddleElements = array[int(length / 2)] + array[int(length / 2 - 1)] - medians = (sumOfMiddleElements) / 2 + if(length%2 == 0): + sumOfMiddleElements = array[int(length / 2)] + array[int(length / 2 - 1)] + medians = ( sumOfMiddleElements) / 2; else: - medians = array[int(length / 2)] + medians = array[int(length/2)] return medians @@ -218,8 +194,8 @@ def deletePos(data, deleteIndexes): def detectAllPeaks(data): peaks, peaksPos = [], [] prev_value = data[0] - for i in range(1, len(data) - 1): - if data[i - 1] < data[i] and data[i + 1] < data[i]: + for i in range(1, len(data)-1): + if data[i-1] < data[i] and data[i+1] < data[i]: peaks.append(data[i]) peaksPos.append(i) return peaksPos, peaks @@ -228,7 +204,7 @@ def detectAllPeaks(data): def higherPeak(data, diff_data, y_threshold): peak, pos = data[0], 0 prev_value = -1 - for i in range(len(diff_data) - 1): + for i in range(len(diff_data)-1): if prev_value > 0 and diff_data[i] < 0 and data[i] > y_threshold: if peak != -1000 and data[i] > peak: peak = data[i] @@ -236,15 +212,11 @@ def higherPeak(data, diff_data, y_threshold): prev_value = diff_data[i] return pos, peak - def limitProeminences(data, peaksPos, threshold): peaksToDelete = [] aux = 0 - for i in range(1, len(peaksPos) - 1): - proeminence = calculateProeminencePeaks( - getValuesArray(data, peaksPos[aux], peaksPos[i + 1]), - peaksPos[i] - peaksPos[aux], - ) + for i in range(1, len(peaksPos)-1): + proeminence = calculateProeminencePeaks(getValuesArray(data, peaksPos[aux], peaksPos[i+1]), peaksPos[i] - peaksPos[aux]) print(proeminence, threshold) if threshold > 0 and proeminence < threshold: # print(" i: ", i, peaksPos[i] - peaksPos[aux], x_threshold) @@ -253,7 +225,6 @@ def limitProeminences(data, peaksPos, threshold): aux = i return peaksToDelete - def limitDistancePeaks(peaksPos, threshold): peaksToDelete = [] aux = 0 @@ -265,7 +236,6 @@ def limitDistancePeaks(peaksPos, threshold): aux = i return peaksToDelete - def correctPeaks(data, peaksPos): if len(peaksPos) == 0: return [], [], [], [] @@ -275,7 +245,7 @@ def correctPeaks(data, peaksPos): if len(peaksPos) == 0: return [], [], [], [] if peaksPos[len(peaksPos) - 1] == len(peaksPos) - 1: - peaksPos = getValuesArray(peaksPos, 0, len(peaksPos) - 1) + peaksPos = getValuesArray(peaksPos, 0, len(peaksPos)-1) if len(peaksPos) == 0: return [], [], [], [] elif len(peaksPos) == 1: @@ -290,34 +260,25 @@ def correctPeaks(data, peaksPos): for i in range(len(peaksPos)): if i == 0: if len(peaksPos) > 1: - new_data = getValuesArray(data, 0, peaksPos[i + 1]) + new_data = getValuesArray(data, 0, peaksPos[i+1]) else: new_data = data peakPos, peak, valleyPos, valley = getCorrectPeak(new_data, peaksPos[i]) elif i == len(peaksPos) - 1: - new_data = getValuesArray(data, peaksPos[i - 1], len(data)) - peakPos, peak, valleyPos, valley = getCorrectPeak( - new_data, peaksPos[i] - peaksPos[i - 1] - ) - peakPos += peaksPos[i - 1] - valleyPos += peaksPos[i - 1] + new_data = getValuesArray(data, peaksPos[i-1], len(data)) + peakPos, peak, valleyPos, valley = getCorrectPeak(new_data, peaksPos[i] - peaksPos[i-1]) + peakPos += peaksPos[i-1] + valleyPos += peaksPos[i-1] else: - new_data = getValuesArray(data, peaksPos[i - 1], peaksPos[i + 1]) - peakPos, peak, valleyPos, valley = getCorrectPeak( - new_data, peaksPos[i] - peaksPos[i - 1] - ) - peakPos += peaksPos[i - 1] - valleyPos += peaksPos[i - 1] + new_data = getValuesArray(data, peaksPos[i-1], peaksPos[i+1]) + peakPos, peak, valleyPos, valley = getCorrectPeak(new_data, peaksPos[i] - peaksPos[i-1]) + peakPos += peaksPos[i-1] + valleyPos += peaksPos[i-1] new_peaksPos[i] = peakPos new_peaks[i] = peak valleysPos[i] = valleyPos valleys[i] = valley - return ( - [int(v) for v in new_peaksPos], - new_peaks, - [int(v) for v in valleysPos], - valleys, - ) + return [int(v) for v in new_peaksPos], new_peaks, [int(v) for v in valleysPos], valleys def getCorrectPeak(data, peakPos): @@ -335,15 +296,10 @@ def getCorrectPeak(data, peakPos): def calculateProeminencePeaks(data, peakPos): - first_proeminence = ( - data[int(peakPos)] - minimumValue(getValuesArray(data, 0, peakPos))[1] - ) - second_proeminence = ( - data[int(peakPos)] - minimumValue(getValuesArray(data, peakPos, len(data)))[1] - ) + first_proeminence = data[int(peakPos)] - minimumValue(getValuesArray(data, 0, peakPos))[1] + second_proeminence = data[int(peakPos)] - minimumValue(getValuesArray(data, peakPos, len(data)))[1] return minimumValue([first_proeminence, second_proeminence])[1] - def calculateProeminencePeaks(data, threshold): peaksPos, peaks = [], [] min_abs_value = minimumValue(data)[1] @@ -355,11 +311,7 @@ def calculateProeminencePeaks(data, threshold): min_value = data[i] if data[i] > max_value: max_value = data[i] - if ( - data[i - 1] < data[i] - and data[i] > data[i + 1] - and data[i] - min_value > threshold - ): + if data[i-1] < data[i] and data[i] > data[i+1] and data[i] - min_value > threshold: peaksPos.append(i) peaks.append(data[i]) max_value = min_abs_value @@ -370,14 +322,13 @@ def calculateProeminencePeaks(data, threshold): # plt.show() return peaksPos, peaks - + def invertSignal(data): mean = np.mean(data) data = (-(data - mean)) + mean return data - def minimumValue(data): minimum = data[0] minPos = 0 @@ -387,7 +338,6 @@ def minimumValue(data): minPos = i return minPos, minimum - def maximumValue(data): maximum = data[0] maxPos = 0 @@ -395,29 +345,26 @@ def maximumValue(data): if data[i] > maximum: maximum = data[i] maxPos = i - + return maxPos, maximum - def detect_SpO2Peaks(red, infrared, fs): # filt_red = bsnb.lowpass(red, 5, 3, fs, True) # filt_infrared = bsnb.lowpass(infrared, 5, 3, fs, True) - filt_red = bsnb.bandpass(red, 0.5, 5, 2, fs, True) - filt_infrared = bsnb.bandpass(infrared, 0.5, 5, 2, fs, True) + filt_red = bsnb.bandpass(red, .5, 5, 2, fs, True) + filt_infrared = bsnb.bandpass(infrared, .5, 5, 2, fs, True) aux = filt_red + filt_infrared # Filter Signal if len(aux) > 1 * fs: - baseline = bsnbSmooth(aux, fs * 1) + baseline = bsnbSmooth(aux, fs*1) else: baseline = mean(aux) signal = aux - baseline - aux = bsnbSmooth(signal, int(0.3 * fs)) + aux = bsnbSmooth(signal, int(.3*fs)) # aux = bsnb.lowpass(signal, 3, 2, fs, True) - peak, peakPos, valley, valleyPos = peakDetection( - aux, proeminence=0.9 * np.std(aux), x_threshold=0.3 * fs - ) + peak, peakPos, valley, valleyPos = peakDetection(aux, proeminence=.9*np.std(aux), x_threshold=.3*fs) # plt.figure() # plt.plot(aux) @@ -436,12 +383,12 @@ def calculate_SpO2(Rpeak, Rvalley, IRpeak, IRvalley): # ir_ratio = (IRpeak - IRvalley) / (IRvalley + (IRpeak - IRvalley)/2) # ir_ratio = (IRpeak - IRvalley) / (IRpeak + IRvalley) # ir_ratio = (IRpeak - IRvalley) / (IRvalley) - ir_ratio = np.log10(IRpeak - IRvalley) - + ir_ratio = np.log10(IRpeak - IRvalley) + R = r_ratio / ir_ratio # R = r_ratio - ir_ratio print(R) - percentage = 110 - 25 * R + percentage = 110 - 25*R # percentage = 103.32 - 6.2033 * R # percentage = 109.06 + 6.3985 * (R**2) - 19.018 * R # r_o = 319.6 @@ -451,31 +398,27 @@ def calculate_SpO2(Rpeak, Rvalley, IRpeak, IRvalley): # percentage = np.abs((r_d - R * ir_d) * 100 / (R * (ir_o - ir_d) + (r_o - r_d))) return percentage - def calculate_SpO2_(Rpeak, Rvalley, IRpeak, IRvalley): red_AC = Rpeak - Rvalley red_DC = ((Rpeak - Rvalley) / 2) + Rvalley infrared_AC = IRpeak - IRvalley infrared_DC = ((IRpeak - IRvalley) / 2) + IRvalley R = (red_AC / red_DC) / (infrared_AC / infrared_DC) - return max(0, min(100, 110 - 25 * R)) - + return max(0, min(100, 110 - 25*R)) def calculate_SpO2_____(Rpeak, Rvalley, IRpeak, IRvalley): - r_ratio = np.log(np.abs(Rpeak / Rvalley)) + r_ratio = np.log(np.abs(Rpeak/Rvalley)) ir_ratio = np.log(np.abs(IRpeak / IRvalley)) R = r_ratio / ir_ratio # percentage = 110 - 25*R - percentage = 113.8 - 24.87 * R + percentage = 113.8 - 24.87*R # print("Percentage: ", percentage) return percentage - -if __name__ == "__main__": +if __name__ == '__main__': from scipy.misc import electrocardiogram - t = np.linspace(-50, 50, 1000) - signal = np.sin(t * 2 * np.pi / 1) * 0.3 + signal = np.sin(t*2*np.pi/1) * .3 plt.figure() plt.plot(t, signal) @@ -483,6 +426,6 @@ def calculate_SpO2_____(Rpeak, Rvalley, IRpeak, IRvalley): plt.show() plt.figure() - plt.plot(electrocardiogram()[:3600]) - plt.plot(invertSignal(electrocardiogram()[:3600])) - plt.show() + plt.plot(electrocardiogram()[: 3600]) + plt.plot(invertSignal(electrocardiogram()[: 3600])) + plt.show() \ No newline at end of file diff --git a/physiological_data/lib/visualization.py b/physiological_data/lib/visualization.py index 463cd1ea..341e649e 100644 --- a/physiological_data/lib/visualization.py +++ b/physiological_data/lib/visualization.py @@ -8,18 +8,14 @@ import acquisition as acq import tools -parser = argparse.ArgumentParser( - description="Identify the path to the folder that contains all signals, photos and PsychoPy files." -) +parser = argparse.ArgumentParser(description="Identify the path to the folder that contains all signals, photos and PsychoPy files.") parser.add_argument("-P", "--path") args = parser.parse_args() if args.path: - folder = rf"D:\Google Drive\Faculdade\Doutoramento\Acquisitions\new_biosignalsnotebooks\acquisitions\Acquisitions\{args.path}" + folder = fr"D:\Google Drive\Faculdade\Doutoramento\Acquisitions\new_biosignalsnotebooks\acquisitions\Acquisitions\{args.path}" else: - raise Exception( - "Please specify the path to the folder containing all signals, photos and PsychoPy files." - ) + raise Exception("Please specify the path to the folder containing all signals, photos and PsychoPy files.") acquisition = acq.Acquisition(folder) acquisition.getBiosignalsSensors(acq.ALL, True) @@ -32,12 +28,12 @@ pictures_filenames = [] timestamps = [] for name in os.listdir(folder): - if "ip_" in name and "-acquisition" in name: - folder_snapshots = os.path.join(os.path.join(folder, name), "snapshot") + if 'ip_' in name and '-acquisition' in name: + folder_snapshots = os.path.join(os.path.join(folder, name), 'snapshot') for filename in os.listdir(folder_snapshots): - if filename.endswith(".jpeg"): + if filename.endswith('.jpeg'): pictures_filenames.append(filename) - timestamps.append(float(filename.split(".jpeg")[0])) + timestamps.append(float(filename.split('.jpeg')[0])) index = tools.closestIndex(timestamps, time_answer[0][0]) timestamps = np.array(timestamps)[np.where(timestamps > time_answer[0][0])[0]] @@ -48,31 +44,31 @@ gs = fig.add_gridspec(6, 12) ax1 = fig.add_subplot(gs[0, 0:3]) -ax1.set_title("fNIRS 1") +ax1.set_title('fNIRS 1') acquisition.getBiosignalsSensors([acq.FNIRS], True) ax1.plot(time_axis, acquisition.signal[:, 1]) ax2 = fig.add_subplot(gs[0, 3:6], sharex=ax1) -ax2.set_title("fNIRS 2") +ax2.set_title('fNIRS 2') ax2.plot(time_axis, acquisition.signal[:, 2]) ax3 = fig.add_subplot(gs[0:3, 6:12]) img = mpimg.imread(os.path.join(folder_snapshots, pictures_filenames[index])) picture = ax3.imshow(img) -ax3.axis("off") +ax3.axis('off') ax4 = fig.add_subplot(gs[1, :3], sharex=ax1) -ax4.set_title("fNIRS 3") +ax4.set_title('fNIRS 3') acquisition.getBiosignalsSensors([acq.FNIRS], True) ax4.plot(time_axis, acquisition.signal[:, 3]) ax5 = fig.add_subplot(gs[1, 3:6], sharex=ax1) -ax5.set_title("fNIRS 4") +ax5.set_title('fNIRS 4') ax5.plot(time_axis, acquisition.signal[:, 4]) ax6 = fig.add_subplot(gs[2, :3], sharex=ax1) -ax6.set_title("EEG 1") +ax6.set_title('EEG 1') acquisition.getBiosignalsSensors([acq.EEG], True) ax6.plot(time_axis, acquisition.signal[:, 1]) @@ -80,21 +76,21 @@ ax7.plot(time_axis, acquisition.signal[:, -1]) ax8 = fig.add_subplot(gs[3, :6], sharex=ax1) -ax8.set_title("EDA") +ax8.set_title('EDA') acquisition.getBiosignalsSensors([acq.EDA], True) # eda = np.where(acquisition.sensors == acq.EDA)[0] ax8.plot(time_axis, acquisition.signal[:, -1]) ax9 = fig.add_subplot(gs[3, 6:], sharex=ax1) -ax9.set_title("ECG") +ax9.set_title('ECG') # ecg = np.where(acquisition.sensors == acq.ECG)[0] # ax2.plot(time_axis, biosignals[:, ecg]) acquisition.getBiosignalsSensors([acq.ECG], True) ax9.plot(time_axis, acquisition.signal[:, -1]) ax10 = fig.add_subplot(gs[4, :6], sharex=ax1) -ax10.set_title("RESP") +ax10.set_title('RESP') # resp = np.where(acquisition.sensors == acq.RESP)[0] # ax4.plot(time_axis, biosignals[:, resp]) acquisition.getBiosignalsSensors([acq.RESP], True) @@ -102,23 +98,23 @@ ax11 = fig.add_subplot(gs[4, 6:], sharex=ax1) -ax11.set_title("Time to Answer") +ax11.set_title('Time to Answer') ax11.scatter(time_answer[0], time_answer[1]) ax12 = fig.add_subplot(gs[5, :4], sharex=ax1) -ax12.set_title("X") +ax12.set_title('X') # acc = np.where(acquisition.sensors == acq.ACC)[0] # ax10.plot(time_axis, biosignals[:, acc[0]]) acquisition.getBiosignalsSensors([acq.ACC], True) ax12.plot(time_axis, acquisition.signal[:, 1]) ax13 = fig.add_subplot(gs[5, 4:8], sharex=ax1) -ax13.set_title("Y") +ax13.set_title('Y') # ax11.plot(time_axis, biosignals[:, acc[1]]) ax13.plot(time_axis, acquisition.signal[:, 2]) ax14 = fig.add_subplot(gs[5, 8:], sharex=ax1) -ax14.set_title("Z") +ax14.set_title('Z') # ax12.plot(time_axis, biosignals[:, acc[2]]) ax14.plot(time_axis, acquisition.signal[:, 3]) @@ -127,35 +123,35 @@ # Comment these lines if you intend to continue the labelling process -with open(os.path.join(folder, "custom_labels.txt"), "w") as f: - f.write("Timestamp,Label\n") +with open(os.path.join(folder, 'custom_labels.txt'), 'w') as f: + f.write('Timestamp,Label\n') -label = "rest" +label = 'rest' i = 0 while i < len(timestamps): timestamp = timestamps[i] - ax1.set_xlim(timestamp - 10, timestamp + 10) - - img = mpimg.imread(os.path.join(folder_snapshots, f"{timestamp}.jpeg")) + ax1.set_xlim(timestamp-10, timestamp+10) + + img=mpimg.imread(os.path.join(folder_snapshots, f"{timestamp}.jpeg")) # ax3.imshow(img) picture.set_data(img) print(timestamp) plt.draw() - plt.show(block=False) + plt.show(block = False) new_label = input("Type a label if something changed, else press Enter: ") - if new_label == "previous": - with open(os.path.join(folder, "custom_labels.txt"), "r") as f: + if new_label == 'previous': + with open(os.path.join(folder, 'custom_labels.txt'), 'r') as f: lines = f.readlines() print(lines[:-1]) - with open(os.path.join(folder, "custom_labels.txt"), "w") as f: + with open(os.path.join(folder, 'custom_labels.txt'), 'w') as f: f.writelines(lines[:-1]) i -= 1 continue - if new_label is not "": + if new_label is not '': label = new_label - - with open(os.path.join(folder, "custom_labels.txt"), "a") as f: + + with open(os.path.join(folder, 'custom_labels.txt'), 'a') as f: f.write(f"{timestamp},{label}\n") - i += 1 + i+=1