{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"}],"dockerImageVersionId":30635,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"papermill":{"default_parameters":{},"duration":102.485133,"end_time":"2024-01-16T12:14:21.10455","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-01-16T12:12:38.619417","version":"2.4.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd, numpy as np, os\nimport matplotlib.pyplot as plt, gc\n\ntrain = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')\nprint('Train shape', train.shape )\ndisplay( train.head() )","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.693827,"end_time":"2024-01-16T12:12:42.606147","exception":false,"start_time":"2024-01-16T12:12:41.91232","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:15:31.849494Z","iopub.execute_input":"2024-03-04T08:15:31.849962Z","iopub.status.idle":"2024-03-04T08:15:32.427644Z","shell.execute_reply.started":"2024-03-04T08:15:31.849921Z","shell.execute_reply":"2024-03-04T08:15:32.426641Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NAMES = ['LL','LP','RP','RR']\n\n\nFEATS = [['Fp1','F7','T3','T5','O1'],\n         ['Fp1','F3','C3','P3','O1'],\n         ['Fp2','F8','T4','T6','O2'],\n         ['Fp2','F4','C4','P4','O2'],]\n\ndirectory_path = '/new/EEG_Spectrograms/'\nif not os.path.exists(directory_path):\n    os.makedirs(directory_path) ","metadata":{"papermill":{"duration":0.017131,"end_time":"2024-01-16T12:12:42.644893","exception":false,"start_time":"2024-01-16T12:12:42.627762","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:15:32.429678Z","iopub.execute_input":"2024-03-04T08:15:32.430095Z","iopub.status.idle":"2024-03-04T08:15:32.436056Z","shell.execute_reply.started":"2024-03-04T08:15:32.430059Z","shell.execute_reply":"2024-03-04T08:15:32.435197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Optional Signal Denoising with Wavelet transform\nWe can optionally denoise the signal before creating the spectrogram. I'm not sure yet if this creates better or worse spectrograms. We can experiment with this. This code comes from Yusaku5738 notebook [here][1] and was suggested by SeshuRajuP in the comments. We have many parent functions to use for denoising. Yusaku5738 suggests using `wavelet = db8`.\n\n[1]: https://www.kaggle.com/code/yusaku5739/eeg-signal-denosing-using-wavelet-transform","metadata":{"papermill":{"duration":0.006935,"end_time":"2024-01-16T12:12:42.659252","exception":false,"start_time":"2024-01-16T12:12:42.652317","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import pywt\nprint(\"The wavelet functions we can use:\")\nprint(pywt.wavelist())\n\nUSE_WAVELET = None","metadata":{"papermill":{"duration":0.540601,"end_time":"2024-01-16T12:12:43.206759","exception":false,"start_time":"2024-01-16T12:12:42.666158","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:15:32.437312Z","iopub.execute_input":"2024-03-04T08:15:32.437845Z","iopub.status.idle":"2024-03-04T08:15:32.844011Z","shell.execute_reply.started":"2024-03-04T08:15:32.437813Z","shell.execute_reply":"2024-03-04T08:15:32.843109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DENOISE FUNCTION\ndef maddest(d, axis=None):\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef denoise(x, wavelet='haar', level=1):    \n    coeff = pywt.wavedec(x, wavelet, mode=\"per\")\n    sigma = (1/0.6745) * maddest(coeff[-level])\n\n    uthresh = sigma * np.sqrt(2*np.log(len(x)))\n    coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])\n\n    ret=pywt.waverec(coeff, wavelet, mode='per')\n    \n    return ret","metadata":{"papermill":{"duration":0.018392,"end_time":"2024-01-16T12:12:43.23335","exception":false,"start_time":"2024-01-16T12:12:43.214958","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:15:32.846181Z","iopub.execute_input":"2024-03-04T08:15:32.846669Z","iopub.status.idle":"2024-03-04T08:15:32.854724Z","shell.execute_reply.started":"2024-03-04T08:15:32.846637Z","shell.execute_reply":"2024-03-04T08:15:32.853656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Spectrograms with Librosa\nWe can use library librosa to create spectrograms. We will save them to disk. For each `eeg_id` we will make 1 spectrogram from the middle 50 seconds. We don't want to use more information than 50 seconds at a time because during test inference, we only have access to 50 seconds of EEG for each test `eeg_id`. We will create spectrograms of `size = 128x256 (freq x time)`.\n\nThe main function is \n\n    mel_spec = librosa.feature.melspectrogram(y=x, sr=200, hop_length=len(x)//256, \n              n_fft=1024, n_mels=128, fmin=0, fmax=20, win_length=128)\n              \nLet's explain these variables.\n* `y` is the input time series signal\n* `sr` is the sampling frequency. In this competition EEG is sample 200 times per sec\n* `hop_length` produces image with `width = len(x)/hop_length`\n* `n_fft` controls vertical resolution and quality of spectrogram\n* `n_mels` produces image with `height = n_mels`\n* `fmin` is smallest frequency in our spectrogram\n* `fmax` is largest frequency in our spectrogram\n* `win_length` controls hortizonal resolution and quality of spectrogram","metadata":{"papermill":{"duration":0.007002,"end_time":"2024-01-16T12:12:43.247431","exception":false,"start_time":"2024-01-16T12:12:43.240429","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import librosa\n\ndef spectrogram_from_eeg(parquet_path, display=False):\n    \n    # LOAD MIDDLE 50 SECONDS OF EEG SERIES\n    eeg = pd.read_parquet(parquet_path)\n    middle = (len(eeg)-10_000)//2\n    eeg = eeg.iloc[middle+3000:middle+7000]\n    \n    # VARIABLE TO HOLD SPECTROGRAM\n    img = np.zeros((128,256,4),dtype='float32')\n    \n    if display: plt.figure(figsize=(10,7))\n    signals = []\n    for k in range(4):\n        COLS = FEATS[k]\n        if k < 4:\n            for kk in range(4):\n        \n                # COMPUTE PAIR DIFFERENCES\n                x = eeg[COLS[kk]].values - eeg[COLS[kk+1]].values\n\n                # FILL NANS\n                m = np.nanmean(x)\n                if np.isnan(x).mean()<1: x = np.nan_to_num(x,nan=m)\n                else: x[:] = 0\n\n                # DENOISE\n                if USE_WAVELET:\n                    x = denoise(x, wavelet=USE_WAVELET)\n                signals.append(x)\n\n                # RAW SPECTROGRAM\n                mel_spec = librosa.feature.melspectrogram(y=x, sr=200, hop_length=len(x)//256, \n                     n_fft=1024, n_mels=128, fmin=0, fmax=20, win_length=128)\n\n                # LOG TRANSFORM\n                width = (mel_spec.shape[1]//32)*32\n                mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max).astype(np.float32)[:,:width]\n\n                # STANDARDIZE TO -1 TO 1\n                mel_spec_db = (mel_spec_db+40)/40 \n                img[:,:,k] += mel_spec_db\n                \n        # AVERAGE THE 4 MONTAGE DIFFERENCES\n            img[:,:,k] /= 4.0\n            \n            \n        \n        if display:\n            plt.subplot(3,2,k+1)\n            plt.imshow(img[:,:,k],aspect='auto',origin='lower')\n            plt.title(f'EEG {eeg_id} - Spectrogram {NAMES[k]}')\n            \n    if display: \n        plt.show()\n        plt.figure(figsize=(10,5))\n        offset = 0\n        for k in range(4):\n            if k>0: offset -= signals[3-k].min()\n            plt.plot(range(4000),signals[k]+offset,label=NAMES[3-k])\n            offset += signals[3-k].max()\n        plt.legend()\n        plt.title(f'EEG {eeg_id} Signals')\n        plt.show()\n        print(); print('#'*25); print()\n        \n    return img","metadata":{"papermill":{"duration":0.032361,"end_time":"2024-01-16T12:12:43.287163","exception":false,"start_time":"2024-01-16T12:12:43.254802","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:16:20.975135Z","iopub.execute_input":"2024-03-04T08:16:20.975889Z","iopub.status.idle":"2024-03-04T08:16:20.988596Z","shell.execute_reply.started":"2024-03-04T08:16:20.975855Z","shell.execute_reply":"2024-03-04T08:16:20.987878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nPATH = '/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/'\nDISPLAY = 4\nEEG_IDS = train.eeg_id.unique()\nall_eegs = {}\n\nfor i,eeg_id in enumerate(EEG_IDS):\n    if (i%100==0)&(i!=0): print(i,', ',end='')\n        \n    # CREATE SPECTROGRAM FROM EEG PARQUET\n    img = spectrogram_from_eeg(f'{PATH}{eeg_id}.parquet', i<DISPLAY)\n    \n    # SAVE TO DISK\n    if i==DISPLAY:\n        print(f'Creating and writing {len(EEG_IDS)} spectrograms to disk... ',end='')\n    np.save(f'{directory_path}{eeg_id}',img)\n    all_eegs[eeg_id] = img\n   \n# SAVE EEG SPECTROGRAM DICTIONARY\nnp.save('eeg_specs',all_eegs)","metadata":{"papermill":{"duration":66.745724,"end_time":"2024-01-16T12:13:50.040166","exception":false,"start_time":"2024-01-16T12:12:43.294442","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-03-04T08:16:21.255527Z","iopub.execute_input":"2024-03-04T08:16:21.256412Z"},"trusted":true},"execution_count":null,"outputs":[]}]}