{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"}],"dockerImageVersionId":30635,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"This notebook aims to understand and explore data and metadata associated with this competition. For a reminder, the goal of this competition is to detect and classify seizures and other types of harmful brain activity. From electroencephalography (EEG) signals recorded from critically ill hospital patients, we must build a classifier that is able to assign for a given EEG time series the correct class based on the detection and recognition of one of 6 patterns:\n* seizure (SZ)\n* generalized periodic discharges (GPD)\n* lateralized periodic discharges (LPD)\n* lateralized rhythmic delta activity (LRDA)\n* generalized rhythmic delta activity (GRDA)\n* “other” for all signals that do not fit inside above categories\n\nAll patterns have specific characteristics based on their frequency, the existence of time gap between pulse, their temporal extent and their location. For more information about this, you can watch this youtube playlist from Fábio A. Nascimento: [EEG Talk - ACNS Critical Care EEG Terminology 2021](https://www.youtube.com/playlist?list=PL1qAb9U_Ln6EO07t6SuqjZdgkwwRNOECr) (careful for loud sound at beginning of videos ^^). Across those videos, Dr. Lawrence Hirsch will explain each of the pattern of interest and what characterize them. For a more thorough explanation and figures, you can also check this article from ACNS: [American Clinical Neurophysiology Society’s Standardized\nCritical Care EEG Terminology: 2021 Version](https://www.acns.org/UserFiles/file/ACNSStandardizedCriticalCareEEGTerminology_rev2021.pdf) and use your favorite LLM to summarize it.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"
\n We will:\n
    \n
  1. Look at training metadata and understand what information each column gives
    2. \n
  2. Look at an example EEG data, all the signals contained in it, what are the characteristic and plot them on time axis\n
  3. Try to filter the EKG signal of this example record, with the help of scipy signal module and fft module\n
\n
","metadata":{}},{"cell_type":"markdown","source":"
\n \"an\n
","metadata":{}},{"cell_type":"markdown","source":"# 📦 Setup and imports","metadata":{}},{"cell_type":"code","source":"import os\n\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport scipy.fft as fft\nimport scipy.signal as signal\nimport seaborn as sns\n\nmatplotlib.rcParams['font.family'] = 'sans-serif'\nmatplotlib.rcParams['figure.figsize'] = (10, 6)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:16.324843Z","iopub.execute_input":"2024-01-12T00:27:16.325274Z","iopub.status.idle":"2024-01-12T00:27:17.557820Z","shell.execute_reply.started":"2024-01-12T00:27:16.325238Z","shell.execute_reply":"2024-01-12T00:27:17.556645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 📂 Load and explore metadata\n\n## Load metadata","metadata":{}},{"cell_type":"code","source":"dataset_path = \"/kaggle/input/hms-harmful-brain-activity-classification\"\ntraining_metadata_path = os.path.join(dataset_path, \"train.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:17.559456Z","iopub.execute_input":"2024-01-12T00:27:17.559995Z","iopub.status.idle":"2024-01-12T00:27:17.564667Z","shell.execute_reply.started":"2024-01-12T00:27:17.559935Z","shell.execute_reply":"2024-01-12T00:27:17.563659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_meta = pd.read_csv(training_metadata_path)\nprint(f\"Training metadata shape: {training_meta.shape}\")\ntraining_meta.head()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:17.645335Z","iopub.execute_input":"2024-01-12T00:27:17.645749Z","iopub.status.idle":"2024-01-12T00:27:17.961891Z","shell.execute_reply.started":"2024-01-12T00:27:17.645709Z","shell.execute_reply":"2024-01-12T00:27:17.960768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Columns: \")\ntraining_meta.columns.to_list()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:18.614266Z","iopub.execute_input":"2024-01-12T00:27:18.614668Z","iopub.status.idle":"2024-01-12T00:27:18.623354Z","shell.execute_reply.started":"2024-01-12T00:27:18.614636Z","shell.execute_reply":"2024-01-12T00:27:18.622133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Comments","metadata":{}},{"cell_type":"code","source":"print(\"For each eeg and spectrogram, there is a unique patient\")\nprint(training_meta.groupby(\"eeg_id\").patient_id.nunique().value_counts())\nprint(training_meta.groupby(\"spectrogram_id\").patient_id.nunique().value_counts())\n\nprint()\nprint(\"-\" * 80)\nprint(\"But one patient can be recorded several times\")\nprint(training_meta.groupby(\"patient_id\").eeg_id.nunique().sort_values(ascending=False))\nprint(training_meta.groupby(\"patient_id\").spectrogram_id.nunique().sort_values(ascending=False))","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:19.599277Z","iopub.execute_input":"2024-01-12T00:27:19.599678Z","iopub.status.idle":"2024-01-12T00:27:19.674077Z","shell.execute_reply.started":"2024-01-12T00:27:19.599643Z","shell.execute_reply":"2024-01-12T00:27:19.673180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Distribution of patterns is rather balanced\nfig = plt.figure()\nax = sns.countplot(\n x=\"expert_consensus\",\n data=training_meta,\n hue=\"expert_consensus\",\n dodge=False\n)\nax.get_legend().get_frame().set_alpha(0.6)\nplt.title(\"Expert consensus distribution\")\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:20.123797Z","iopub.execute_input":"2024-01-12T00:27:20.124865Z","iopub.status.idle":"2024-01-12T00:27:21.001129Z","shell.execute_reply.started":"2024-01-12T00:27:20.124816Z","shell.execute_reply":"2024-01-12T00:27:21.000007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# only label_id is unique for each row\ntraining_meta.nunique()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:21.003090Z","iopub.execute_input":"2024-01-12T00:27:21.003432Z","iopub.status.idle":"2024-01-12T00:27:21.038622Z","shell.execute_reply.started":"2024-01-12T00:27:21.003397Z","shell.execute_reply":"2024-01-12T00:27:21.037511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# A single spectrogram can regroup several eegs\nprint(training_meta.groupby(\"spectrogram_id\").eeg_id.nunique().sort_values(ascending=False))\nprint()\nprint(\"-\" * 80)\nprint()\ntraining_meta[training_meta.spectrogram_id == 764146759]","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:21.953220Z","iopub.execute_input":"2024-01-12T00:27:21.953626Z","iopub.status.idle":"2024-01-12T00:27:21.992989Z","shell.execute_reply.started":"2024-01-12T00:27:21.953591Z","shell.execute_reply":"2024-01-12T00:27:21.991969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metadata column description","metadata":{}},{"cell_type":"markdown","source":"We will explain meaning of the different columns. Explanations are mainly based on competition data description\navailable at [HMS - Harmful Brain Activity Classification](https://www.kaggle.com/competitions/hms-harmful-brain-activity-classification/data):\n* `eeg_id`: A unique identifier for the entire EEG recording. Those identifies also the parquet files of the eeg in the **train_eegs** folder. Note that as each row denominates a particular subsample of an EEG, there can be several rows with same eeg_id.\n* `eeg_sub_id`: An ID for the specific 50 second long subsample this row's labels apply to. `eeg_sub_id` are assigned in chronological order for each `eeg_id` beginning with 0. The subsample is shifted from the beginning of the `eeg_id` record by `eeg_label_offset_seconds`. Note that even if subsample is 50 seconds long, the annotation is done by looking at the central 10 seconds.\n* `eeg_label_offset_seconds`: The time between the beginning of the consolidated EEG and this subsample. Time shift.\n* `spectrogram_id`: A unique identifier for an entire spectrogram. A spectrogram can regroup several EEGs records. Note that `spectrogram_id` are name of spectrogram files in folder **train_spectrograms**. For more information on what is a [spectrogram](https://en.wikipedia.org/wiki/Spectrogram).\n* `spectrogram_sub_id`: An ID for the specific 10 minute subsample this row's labels apply to. Same concept as for `eeg_sub_id`.\n* `spectogram_label_offset_seconds`: The time between the beginning of the consolidated spectrogram and this subsample. Simple time shift.\n* `label_id`: An ID for this set of labels. Unique for each row.\n* `patient_id`: An ID for the patient who donated the data. An `eeg_id` or `spectrogram_id` is associated with a single `patient_id` but a single patient can be associated with several eegs and spectrograms.\n* `expert_consensus`: The consensus annotator label. Provided for convenience only. Give the final decision for classifying of the 10 seconds central window of a given subsample in one of the 6 available category.\n* `[seizure/lpd/gpd/lrda/grda/other]_vote`: The count of annotator votes for a given brain activity class. The full names of the activity classes are as follows: lpd: lateralized periodic discharges, gpd: generalized periodic discharges, lrd: lateralized rhythmic delta activity, and grda: generalized rhythmic delta activity. Note that total number of annotators can vary for each row, even within the same eeg. Size of cohort of experts range from 1 to 28.","metadata":{}},{"cell_type":"markdown","source":"# 🧠 Let's look at an example eeg","metadata":{}},{"cell_type":"markdown","source":"## Recording parameters","metadata":{}},{"cell_type":"code","source":"fs = 200 # data sampled at 200 Hz\nwindow_length = 10 # 10 seconds central window\nsubsample_length = 50 # 50 seconds subsampling for each row","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:23.884923Z","iopub.execute_input":"2024-01-12T00:27:23.885332Z","iopub.status.idle":"2024-01-12T00:27:23.890262Z","shell.execute_reply.started":"2024-01-12T00:27:23.885299Z","shell.execute_reply":"2024-01-12T00:27:23.889277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Choose particular eeg_id and load eeg","metadata":{}},{"cell_type":"code","source":"eeg_path_template = dataset_path + \"/train_eegs/{eeg_id}.parquet\"","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:24.671870Z","iopub.execute_input":"2024-01-12T00:27:24.672600Z","iopub.status.idle":"2024-01-12T00:27:24.677058Z","shell.execute_reply.started":"2024-01-12T00:27:24.672558Z","shell.execute_reply":"2024-01-12T00:27:24.675942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we chose this id as it is associated with a heavily recognized pattern \neeg_id = 2900632927\ntraining_meta[training_meta.eeg_id == eeg_id]","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:25.762093Z","iopub.execute_input":"2024-01-12T00:27:25.762772Z","iopub.status.idle":"2024-01-12T00:27:25.782527Z","shell.execute_reply.started":"2024-01-12T00:27:25.762736Z","shell.execute_reply":"2024-01-12T00:27:25.781351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_eeg = pd.read_parquet(eeg_path_template.format(eeg_id=eeg_id))\nprint(\"EEG record shape: \", df_eeg.shape)\ndf_eeg.head()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:26.037822Z","iopub.execute_input":"2024-01-12T00:27:26.038202Z","iopub.status.idle":"2024-01-12T00:27:26.239893Z","shell.execute_reply.started":"2024-01-12T00:27:26.038172Z","shell.execute_reply":"2024-01-12T00:27:26.239158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_eeg.columns.to_list())","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:26.515638Z","iopub.execute_input":"2024-01-12T00:27:26.516615Z","iopub.status.idle":"2024-01-12T00:27:26.522354Z","shell.execute_reply.started":"2024-01-12T00:27:26.516578Z","shell.execute_reply":"2024-01-12T00:27:26.521149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have 20 signals, 19 corresponding to the EEG measures, which are voltage fluctuations recorded by electrodes placed on the scalp\nand a last one corresponding to `EKG`, which is an electrocardiogram of the patient.\n\nThe EEG measures give a view of the electrical activity of the brain. The different signals names (`Fp1`, `F3`, ...) corresponds to standardized locations where the electrodes will be attached. The locations are generally associated with a particular area of the brain. The letters design a particular lobe, such as pre-frontal (Fp), frontal (F), temporal (T), parietal (P), occipital (O), and central (C). Even numbers in the electrode names indicate it is placed on the right side of the head. Odd number indicate it is placed on left side. z (for zero) indicate the electrode is placed in the central plane and serves as reference point. For more information, you can see the [10–20 system (EEG)](https://en.wikipedia.org/wiki/10%E2%80%9320_system_(EEG)).\n\nThe unit of measurement is $\\mu V$.","metadata":{}},{"cell_type":"markdown","source":"## Add time","metadata":{}},{"cell_type":"code","source":"df_eeg[\"time\"] = df_eeg.index / fs\ndf_eeg.set_index(\"time\", inplace=True)\ndf_eeg.index","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:27.565938Z","iopub.execute_input":"2024-01-12T00:27:27.566710Z","iopub.status.idle":"2024-01-12T00:27:27.576448Z","shell.execute_reply.started":"2024-01-12T00:27:27.566675Z","shell.execute_reply":"2024-01-12T00:27:27.575336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plot one electrode","metadata":{}},{"cell_type":"code","source":"plt.plot(df_eeg.index, df_eeg[\"Fp1\"])\nplt.xlabel(\"Time (s)\")\nplt.ylabel(\"Amplitude (uV)\")\nplt.grid()\nplt.title(\"Fp1\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:28.325696Z","iopub.execute_input":"2024-01-12T00:27:28.326127Z","iopub.status.idle":"2024-01-12T00:27:28.911275Z","shell.execute_reply.started":"2024-01-12T00:27:28.326092Z","shell.execute_reply":"2024-01-12T00:27:28.910029Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can focus on a particular window of time.","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:11:33.477698Z","iopub.execute_input":"2024-01-12T00:11:33.478099Z","iopub.status.idle":"2024-01-12T00:11:33.486257Z","shell.execute_reply.started":"2024-01-12T00:11:33.478072Z","shell.execute_reply":"2024-01-12T00:11:33.483979Z"}}},{"cell_type":"code","source":"x_start = 25.0\nx_end = 35.0\nplt.plot(df_eeg.loc[x_start:x_end].index, df_eeg.loc[x_start:x_end][\"Fp1\"])\nplt.grid()\nplt.title(\"Fp1\")\nplt.xlabel(\"Time (s)\")\nplt.ylabel(\"Amplitude (uV)\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:28.913565Z","iopub.execute_input":"2024-01-12T00:27:28.914017Z","iopub.status.idle":"2024-01-12T00:27:29.294006Z","shell.execute_reply.started":"2024-01-12T00:27:28.913973Z","shell.execute_reply":"2024-01-12T00:27:29.292809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ⚡ Whole set of electrodes","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=len(df_eeg.columns), ncols=1, sharex=True, figsize=(30, 50))\nfor i, col_name in enumerate(df_eeg.columns):\n ax = sns.lineplot(data=df_eeg, x=df_eeg.index, y=col_name, ax=axes[i])\n ax.set_title(col_name)\n ax.grid()\nplt.subplots_adjust()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:30.001389Z","iopub.execute_input":"2024-01-12T00:27:30.001801Z","iopub.status.idle":"2024-01-12T00:27:40.347244Z","shell.execute_reply.started":"2024-01-12T00:27:30.001766Z","shell.execute_reply":"2024-01-12T00:27:40.345858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not very insightful, signals are very noisy. We can take as extreme example the electrocardiogram `EKG` (or ECG), which is completly unreadable. Yet in the example figures (see one below representing a LPD identified in a portion of an EEG), we have at the bottom right an interpretable EKG in a 10 second window. A possible cause of the noise in our signal is the existence of an interfering signal that add to the true signal.\n\nOne idea would be to smooth the signal with a simple moving average filter. Or filter certain frequencies with a signal processing filter.\nLuckily, the python module [scipy.signal](https://docs.scipy.org/doc/scipy/reference/signal.html) offers us whole range of tools for signal processing.","metadata":{}},{"cell_type":"markdown","source":"
\n \"an\n
","metadata":{}},{"cell_type":"markdown","source":"# 🔎 Analysing the EKG signal","metadata":{}},{"cell_type":"markdown","source":"We will focus on a 20 seconds length centered window of the `EKG` signal.\n\nWe will inspect the frequency components of the EKG signal and use established litterature to assist us in processing this signal into an interpretable and meaningful one.","metadata":{}},{"cell_type":"markdown","source":"## Extract 20 seconds window in the middle of the EKG record","metadata":{}},{"cell_type":"code","source":"# extract EKG signal in a 20 second centered window\nekg_subsample = df_eeg.loc[30.0:50.0][\"EKG\"].values.copy()[:-1]\ntd = np.arange(30., 50., 1 / fs) # associated time domain","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:46.160520Z","iopub.execute_input":"2024-01-12T00:27:46.160933Z","iopub.status.idle":"2024-01-12T00:27:46.167763Z","shell.execute_reply.started":"2024-01-12T00:27:46.160890Z","shell.execute_reply":"2024-01-12T00:27:46.166475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(6, 4))\nplt.plot(td, ekg_subsample)\nplt.xlabel(\"Time (s)\")\nplt.ylabel(\"Amplitude (uV)\")\nplt.grid()\nplt.title(\"EKG\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:46.612309Z","iopub.execute_input":"2024-01-12T00:27:46.613149Z","iopub.status.idle":"2024-01-12T00:27:47.144041Z","shell.execute_reply.started":"2024-01-12T00:27:46.613110Z","shell.execute_reply":"2024-01-12T00:27:47.142734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plot Fourier transform","metadata":{}},{"cell_type":"markdown","source":"NOTE: my memories of my signal processing courses are a little bit old. But I will always try to put resources to give the justification for each step taken 🙃.","metadata":{}},{"cell_type":"markdown","source":"The Fourier transform (FT) is well known in signal processing. It is a transformation that associate to a function defined in the temporal domain a function defined in the frequency domaine. This new function describes the frequency spectrum of the input temporal function. So to make it brief, we can consider our input EKG signal made up of different frequency components (for exemple a 1Hz component for heartbeats rate at 60/min) and the Fourier transform will allow us to decompose the EKG signal into those different frequency components and for each of those components we will have amplitude and phase. For more information about the Fourier transform, you can check out the [wikipedia page](https://en.wikipedia.org/wiki/Fourier_transform).\n\nAs we have discrete signals, we will have to compute a [discrete Fourier transform (DFT)](https://en.wikipedia.org/wiki/Discrete_Fourier_transform), and for that we can make use of [Fast Fourier Transform (FFT)](https://en.wikipedia.org/wiki/Fast_Fourier_transform), which is an efficient algorithm for computing DFT.\n\nScipy offers us a module [fft](https://docs.scipy.org/doc/scipy/reference/fft.html#module-scipy.fft) for Fourier analysis. We will rely upon it.","metadata":{}},{"cell_type":"code","source":"y_ekg = fft.fft(ekg_subsample)\n\n# array of frequencies corresponding to the FT points\n# determined by number of points and sampling interval (inverse of sampling frequency)\nfd = fft.fftfreq(n=len(y_ekg), d=1 / fs)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:48.433908Z","iopub.execute_input":"2024-01-12T00:27:48.434294Z","iopub.status.idle":"2024-01-12T00:27:48.445266Z","shell.execute_reply.started":"2024-01-12T00:27:48.434262Z","shell.execute_reply":"2024-01-12T00:27:48.444019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to plot the amplitude of each frequency, we take the absolute value of the complex number\n# also as the FT is symmetric between (-fs/2, fs/2), we only plot the first half\nfig = plt.figure(figsize=(10, 6))\nplt.plot(fd[:len(fd) // 2], np.abs(y_ekg)[:len(fd) // 2])\nplt.grid()\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Amplitude\")\nplt.title(\"EKG spectrum\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:48.825290Z","iopub.execute_input":"2024-01-12T00:27:48.825706Z","iopub.status.idle":"2024-01-12T00:27:49.130117Z","shell.execute_reply.started":"2024-01-12T00:27:48.825674Z","shell.execute_reply":"2024-01-12T00:27:49.129006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ok so we have some good information here. So we can see that we have the amplitude of the frequency components from 0 Hz to 100 Hz. This value of 100 Hz is quite important, it is related to the [Nyquist–Shannon sampling theorem](https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem) which states that a continuous signal $x(t)$ which does not contain frequency component higher than $B\\; Hz$ can be fully reconstructed by a discrete sampling $\\{x(n \\times \\frac{1}{fs})\\}_n$ if the sampling frequency $f_s$ is at least two times the higher $x(t)$ frequency: $f_s > 2B$.\n\nOur original `EKG` signal is sampled at 200 Hz. So the highest frequency we can have reliable information on is: $200/2 = 100 \\;Hz$\n\nWe observe two major amplitude peaks:\n* one at 0 Hz: it is the \"non-oscillating\" component, the mean component. As we are measuring potential in volt (V), the measure is equivalent up to a constant. So we can subtract mean from the signal for convenience.\n* one huge at approximately 60 Hz: it is likely to be due to powerline interference (50 Hz or 60 Hz). We will surely want to filter this component.\n\nWe can also zoom on low frequencies (< 2 Hz)","metadata":{}},{"cell_type":"code","source":"delta_f = 1 / len(y_ekg) * fs # frequency resolution (step between each frequency point)\nprint(f\"Frequency resolution: {delta_f} Hz\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:49.928655Z","iopub.execute_input":"2024-01-12T00:27:49.929630Z","iopub.status.idle":"2024-01-12T00:27:49.935487Z","shell.execute_reply.started":"2024-01-12T00:27:49.929591Z","shell.execute_reply":"2024-01-12T00:27:49.934404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(10, 6))\nplt.plot(fd[:40], np.abs(y_ekg)[:40])\nplt.grid()\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Amplitude\")\nplt.title(\"EKG spectrum\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:51.148381Z","iopub.execute_input":"2024-01-12T00:27:51.148760Z","iopub.status.idle":"2024-01-12T00:27:51.523738Z","shell.execute_reply.started":"2024-01-12T00:27:51.148728Z","shell.execute_reply":"2024-01-12T00:27:51.522619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Apart from the mean component (0 Hz), there are very low frequencies (0.05-0.5 Hz) components\nHowever heart rate typically ranges between 60 and 100 bpm (1-1.6 Hz). So those very low frequencies are \nlikely not to be related to the heart rate but maybe to the respiration rate or slow movements of the patient.\n\nWe can filter them with a highpass filter.","metadata":{}},{"cell_type":"markdown","source":"# 🔉 Filtering !","metadata":{}},{"cell_type":"markdown","source":"For this part, I used this resource from GE Healthcare: [A Guide to ECG Signal Filtering](https://www.gehealthcare.com/insights/article/a-guide-to-ecg-signal-filtering), which gives insight on which frequencies to filter.","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:16:55.680211Z","iopub.execute_input":"2024-01-12T00:16:55.680576Z","iopub.status.idle":"2024-01-12T00:16:55.688568Z","shell.execute_reply.started":"2024-01-12T00:16:55.680547Z","shell.execute_reply":"2024-01-12T00:16:55.686964Z"}}},{"cell_type":"markdown","source":"We first remove the mean component.","metadata":{}},{"cell_type":"code","source":"ekg_subsample_centered = ekg_subsample - ekg_subsample.mean()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:53.486469Z","iopub.execute_input":"2024-01-12T00:27:53.487046Z","iopub.status.idle":"2024-01-12T00:27:53.493080Z","shell.execute_reply.started":"2024-01-12T00:27:53.486984Z","shell.execute_reply":"2024-01-12T00:27:53.491924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then a highpass filter to remove frequency components below 0.5 Hz.","metadata":{}},{"cell_type":"code","source":"f_cutoff = 0.5 # cut-off frequency in Hz\n\n# generate butterworth filter (by generating the coefficients)\nb, a = signal.butter(N=4, Wn=f_cutoff, btype='highpass', fs=fs, analog=False, output='ba')","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:54.424198Z","iopub.execute_input":"2024-01-12T00:27:54.425017Z","iopub.status.idle":"2024-01-12T00:27:54.434809Z","shell.execute_reply.started":"2024-01-12T00:27:54.424968Z","shell.execute_reply":"2024-01-12T00:27:54.433453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ekg_subsample_high = signal.filtfilt(b, a, ekg_subsample_centered) ","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:54.446342Z","iopub.execute_input":"2024-01-12T00:27:54.446745Z","iopub.status.idle":"2024-01-12T00:27:54.454406Z","shell.execute_reply.started":"2024-01-12T00:27:54.446710Z","shell.execute_reply":"2024-01-12T00:27:54.453292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# notch filter to remove power line interference\nf0 = 60.0 # Frequency to be removed from signal (Hz)\nQ = 30.0 # Quality factor\n\nb_notch, a_notch = signal.iirnotch(w0=f0, Q=Q, fs=fs)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:54.484253Z","iopub.execute_input":"2024-01-12T00:27:54.484681Z","iopub.status.idle":"2024-01-12T00:27:54.490508Z","shell.execute_reply.started":"2024-01-12T00:27:54.484646Z","shell.execute_reply":"2024-01-12T00:27:54.489201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ekg_subsample_filtered = signal.filtfilt(b_notch, a_notch, ekg_subsample_high)","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:54.510078Z","iopub.execute_input":"2024-01-12T00:27:54.510448Z","iopub.status.idle":"2024-01-12T00:27:54.516384Z","shell.execute_reply.started":"2024-01-12T00:27:54.510419Z","shell.execute_reply":"2024-01-12T00:27:54.515295Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# lets plot the Fourier transform of the filtered signal\ny_filtered = fft.fft(ekg_subsample_filtered)\n\nfig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8))\n\naxes[0].plot(fd[:len(fd) // 2], np.abs(y_filtered)[:len(fd) // 2])\naxes[0].set_title(\"EKG spectrum after filtering\")\naxes[0].grid()\naxes[0].set_xlabel(\"Frequency (Hz)\")\naxes[0].set_ylabel(\"Amplitude\")\n\naxes[1].plot(fd[:40], np.abs(y_filtered)[:40])\naxes[1].set_title(\"EKG spectrum after filtering (zoomed on low frequencies)\")\naxes[1].grid()\naxes[1].set_xlabel(\"Frequency (Hz)\")\naxes[1].set_ylabel(\"Amplitude\")\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:54.543987Z","iopub.execute_input":"2024-01-12T00:27:54.544368Z","iopub.status.idle":"2024-01-12T00:27:55.327333Z","shell.execute_reply.started":"2024-01-12T00:27:54.544328Z","shell.execute_reply":"2024-01-12T00:27:55.326238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have successfully filtered our original signal. We can now plot the temporal filtered EKG signal and compare it to the unfiltered one.","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:18:07.223114Z","iopub.execute_input":"2024-01-12T00:18:07.223630Z","iopub.status.idle":"2024-01-12T00:18:07.230151Z","shell.execute_reply.started":"2024-01-12T00:18:07.223597Z","shell.execute_reply":"2024-01-12T00:18:07.228785Z"}}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8), sharex=True)\naxes[0].plot(td, ekg_subsample)\naxes[0].set_title(\"EKG signal\")\naxes[0].grid()\naxes[0].set_ylabel(\"Amplitude\")\n\naxes[1].plot(td, ekg_subsample_filtered)\naxes[1].set_title(\"EKG signal after filtering\")\naxes[1].grid()\naxes[1].set_xlabel(\"Time (s)\")\naxes[1].set_ylabel(\"Amplitude\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:55.329315Z","iopub.execute_input":"2024-01-12T00:27:55.329631Z","iopub.status.idle":"2024-01-12T00:27:56.255047Z","shell.execute_reply.started":"2024-01-12T00:27:55.329602Z","shell.execute_reply":"2024-01-12T00:27:56.253899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are some boundary effect. For the moment we will ignore those extreme points.","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:18:25.151501Z","iopub.execute_input":"2024-01-12T00:18:25.153050Z","iopub.status.idle":"2024-01-12T00:18:25.159648Z","shell.execute_reply.started":"2024-01-12T00:18:25.152995Z","shell.execute_reply":"2024-01-12T00:18:25.158258Z"}}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8), sharex=True)\naxes[0].plot(td[100:-100], ekg_subsample[100:-100])\naxes[0].set_title(\"EKG signal\")\naxes[0].grid()\naxes[0].set_ylabel(\"Amplitude\")\n\naxes[1].plot(td[100:-100], ekg_subsample_filtered[100:-100])\naxes[1].set_title(\"EKG signal after filtering\")\naxes[1].grid()\naxes[1].set_xlabel(\"Time (s)\")\naxes[1].set_ylabel(\"Amplitude\")","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:27:56.256835Z","iopub.execute_input":"2024-01-12T00:27:56.257214Z","iopub.status.idle":"2024-01-12T00:27:57.135518Z","shell.execute_reply.started":"2024-01-12T00:27:56.257181Z","shell.execute_reply":"2024-01-12T00:27:57.134131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Great 🙌 !","metadata":{"execution":{"iopub.status.busy":"2024-01-12T00:20:20.419870Z","iopub.execute_input":"2024-01-12T00:20:20.420256Z","iopub.status.idle":"2024-01-12T00:20:20.426503Z","shell.execute_reply.started":"2024-01-12T00:20:20.420224Z","shell.execute_reply":"2024-01-12T00:20:20.424914Z"}}},{"cell_type":"markdown","source":"# 🌛 TO BE CONTINUED (it's 1 AM in France)","metadata":{}}]}