{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.13","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"},{"sourceId":7611741,"sourceType":"datasetVersion","datasetId":4432380},{"sourceId":7651147,"sourceType":"datasetVersion","datasetId":4460383}],"dockerImageVersionId":30646,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Single numpy file with spectrograms\n\nv1: using only F3 and F4, as two channels.\n\n\nThere is room for tuning when producing the spectrograms: \n  - see influence of window shape and size: https://www.audiolabs-erlangen.de/resources/MIR/FMP/C2/C2_STFT-Window.html  \n  - apply logarithm\n\nTo modify train/val/test sets, go to keras/03_stratified.ipynb\n\n- Band pass applied before the SFFT: Not needed. Comparison applying band pass before with slicing the spectrogram without filtering.\n","metadata":{}},{"cell_type":"code","source":"!pip install --upgrade /kaggle/input/hms-libraries/scipy-1.12.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom scipy.signal import ShortTimeFFT\nfrom scipy.signal.windows import gaussian\nimport random\n\n# base_dir = \"../../kaggle_data/hms\"\n# base_dir = \"../../data/hms\"\nbase_dir = \"/kaggle/input/hms-harmful-brain-activity-classification\"\n\n# data_dir = '../data'\ndata_dir = '/kaggle/input/hms-indices-train-val-test-v1'\n\n# output_dir = '../data/'\noutput_dir = ''\n\nfs = 200  # Sample rate.\n\ndf_traincsv = pd.read_csv(f'{base_dir}/train.csv')\n\nTARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']\n\ndf_traincsv = pd.read_csv(f'{base_dir}/train.csv')\ndf_traincsv.loc[df_traincsv.expert_consensus == 'Seizure', 'target'] = 0\ndf_traincsv.loc[df_traincsv.expert_consensus == 'LPD', 'target'] = 1\ndf_traincsv.loc[df_traincsv.expert_consensus == 'GPD', 'target'] = 2\ndf_traincsv.loc[df_traincsv.expert_consensus == 'LRDA', 'target'] = 3\ndf_traincsv.loc[df_traincsv.expert_consensus == 'GRDA', 'target'] = 4\ndf_traincsv.loc[df_traincsv.expert_consensus == 'Other', 'target'] = 5\n\n# Transform votes into percentages.\ndf_traincsv['sum_votes'] = df_traincsv.seizure_vote + df_traincsv.lpd_vote + df_traincsv.gpd_vote\t+ df_traincsv.lrda_vote + df_traincsv.grda_vote + df_traincsv.other_vote\ndf_traincsv['seizure_vote'] = df_traincsv.seizure_vote/df_traincsv.sum_votes\ndf_traincsv['lpd_vote'] = df_traincsv.lpd_vote/df_traincsv.sum_votes\ndf_traincsv['gpd_vote'] = df_traincsv.gpd_vote/df_traincsv.sum_votes\ndf_traincsv['lrda_vote'] = df_traincsv.lrda_vote/df_traincsv.sum_votes\ndf_traincsv['grda_vote'] = df_traincsv.grda_vote/df_traincsv.sum_votes\ndf_traincsv['other_vote'] = df_traincsv.other_vote/df_traincsv.sum_votes\n\nidxs_train = np.load(f'{data_dir}/03_stratified_v1_idxs_train.npy')\nidxs_val = np.load(f'{data_dir}/03_stratified_v1_idxs_val.npy')\nidxs_test = np.load(f'{data_dir}/03_stratified_v1_idxs_test.npy')\ndf_train = df_traincsv.loc[idxs_train]\ndf_val = df_traincsv.loc[idxs_val]\ndf_test = df_traincsv.loc[idxs_test]\n\nprint(\"Added target column. Transformed into percentages.\")\nprint(\"Train:\", len(df_train))\nprint(\"Val:\", len(df_val))\nprint(\"Test:\", len(df_test))\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check there is no patient_id in more than one set.\nids = np.unique(df_train['patient_id'])\nprint(len(df_val.loc[df_val['patient_id'].isin(ids)]))\nprint(len(df_test.loc[df_test['patient_id'].isin(ids)]))\nids = np.unique(df_val['patient_id'])\nprint(len(df_train.loc[df_train['patient_id'].isin(ids)]))\nprint(len(df_test.loc[df_test['patient_id'].isin(ids)]))\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Version v1\n\nUsing only F3 and F4, in two channels.","metadata":{}},{"cell_type":"code","source":"#\n# Training set\n#\n\nn_channels = 2\nmax_freq = 20  # Only keep freqs below this number.\n# min_freq = 8  # Only keep freqs above this number.\n# time_window = 10  # 10 second event.\n\n#\n# SFT setup: some tuning may be applied.\n#\ng_std = 24  # standard deviation for Gaussian window in samples\nhop = 3\nwin_width = 48  # Pick an odd number.\nmfft = 400\nwin = gaussian(win_width, std=g_std, sym=True)  # symmetric Gaussian wind.\nSFT = ShortTimeFFT(win, hop=hop, fs=fs, mfft=mfft)\n\nn2 = int(max_freq/SFT.delta_f)  # Number of bins below max_freq (30 Hz).\n# Dimensions of Sx, check the code in playing building spectrograms.\ndim1 = 39\ndim2 = 682\n# Each channel is appended to all.\n# all = np.array([]).reshape(0, dim1, dim2, n_channels)\nspecs = np.empty((len(df_train), dim1, dim2, n_channels))\n\n# item: [eeg_id, eeg_sub_id, idx in all (1st index), target,\n#       seizure_vote, lpd_vote, gpd_vote, lrda_vote,\n#       grda_vote, other_vote]\nitems = np.array([], dtype=float).reshape(0,10)\n\nfor i in np.arange(len(df_train)):\n    if i%400 == 0:\n        print(f'{i} files loaded.', end='\\r')\n    item = df_train.iloc[i]\n    eeg = pd.read_parquet(f'{base_dir}/train_eegs/{item.eeg_id}.parquet')\n    eeg = eeg.interpolate(limit_direction='both') # <<<<< Interpolation\n\n    # 10 second eeg sub samples \n    offset = int(item.eeg_label_offset_seconds)\n    start = (offset + 20) * fs\n    end = (offset + 30) * fs\n    eeg_sub_10 = eeg[start:end]\n\n    N = eeg_sub_10.shape[0]\n    t_x = np.arange(N) * 1/fs  # time indexes for signal\n\n    X = np.empty((1, dim1, dim2, n_channels))\n    # X[0,:,:,c] = Sx[n1:n2,:].copy()\n\n    # for c in np.arange(n_channels):\n    x = eeg_sub_10['F3'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,0] = Sx[1:n2,:].copy()\n    x = eeg_sub_10['F4'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,1] = Sx[1:n2,:].copy()\n\n    xitem = np.array([item.eeg_id, item.eeg_sub_id, i, item.target,\n                    item.seizure_vote, item.lpd_vote, item.gpd_vote,\n                    item.lrda_vote, item.grda_vote, item.other_vote],\n                    dtype=float).reshape(1,10)\n    items = np.concatenate([items, xitem])\n\nfilename = '03_single_spectrograms_v1_train'     \nprint(f'Saving to {filename}.npy')\nprint(f'Saving to {filename}_items.npy')\nnp.save(f'{output_dir}{filename}.npy', specs)\nnp.save(f'{output_dir}{filename}_items.npy', items)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#\n# Validation set\n#\n\nn_channels = 2\nmax_freq = 20  # Only keep freqs below this number.\n# min_freq = 8  # Only keep freqs above this number.\n# time_window = 10  # 10 second event.\n\n#\n# SFT setup: some tuning may be applied.\n#\ng_std = 24  # standard deviation for Gaussian window in samples\nhop = 3\nwin_width = 48  # Pick an odd number.\nmfft = 400\nwin = gaussian(win_width, std=g_std, sym=True)  # symmetric Gaussian wind.\nSFT = ShortTimeFFT(win, hop=hop, fs=fs, mfft=mfft)\n\nn2 = int(max_freq/SFT.delta_f)  # Number of bins below max_freq (30 Hz).\n# Dimensions of Sx, check the code in playing building spectrograms.\ndim1 = 39\ndim2 = 682\n# Each channel is appended to all.\n# all = np.array([]).reshape(0, dim1, dim2, n_channels)\nspecs = np.empty((len(df_val), dim1, dim2, n_channels))\n\n# item: [eeg_id, eeg_sub_id, idx in all (1st index), target,\n#       seizure_vote, lpd_vote, gpd_vote, lrda_vote,\n#       grda_vote, other_vote]\nitems = np.array([], dtype=float).reshape(0,10)\n\nfor i in np.arange(len(df_val)):\n    if i%400 == 0:\n        print(f'{i} files loaded.', end='\\r')\n    item = df_val.iloc[i]\n    eeg = pd.read_parquet(f'{base_dir}/train_eegs/{item.eeg_id}.parquet')\n    eeg = eeg.interpolate(limit_direction='both') # <<<<< Interpolation\n\n    # 10 second eeg sub samples \n    offset = int(item.eeg_label_offset_seconds)\n    start = (offset + 20) * fs\n    end = (offset + 30) * fs\n    eeg_sub_10 = eeg[start:end]\n\n    N = eeg_sub_10.shape[0]\n    t_x = np.arange(N) * 1/fs  # time indexes for signal\n\n    X = np.empty((1, dim1, dim2, n_channels))\n    # X[0,:,:,c] = Sx[n1:n2,:].copy()\n\n    # for c in np.arange(n_channels):\n    x = eeg_sub_10['F3'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,0] = Sx[1:n2,:].copy()\n    x = eeg_sub_10['F4'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,1] = Sx[1:n2,:].copy()\n\n    xitem = np.array([item.eeg_id, item.eeg_sub_id, i, item.target,\n                    item.seizure_vote, item.lpd_vote, item.gpd_vote,\n                    item.lrda_vote, item.grda_vote, item.other_vote],\n                    dtype=float).reshape(1,10)\n    items = np.concatenate([items, xitem])\n\nfilename = '03_single_spectrograms_v1_val'     \nprint(f'Saving to {filename}.npy')\nprint(f'Saving to {filename}_items.npy')\nnp.save(f'{output_dir}{filename}.npy', specs)\nnp.save(f'{output_dir}{filename}_items.npy', items)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#\n# Testing set\n#\n\nn_channels = 2\nmax_freq = 20  # Only keep freqs below this number.\n# min_freq = 8  # Only keep freqs above this number.\n# time_window = 10  # 10 second event.\n\n#\n# SFT setup: some tuning may be applied.\n#\ng_std = 24  # standard deviation for Gaussian window in samples\nhop = 3\nwin_width = 48  # Pick an odd number.\nmfft = 400\nwin = gaussian(win_width, std=g_std, sym=True)  # symmetric Gaussian wind.\nSFT = ShortTimeFFT(win, hop=hop, fs=fs, mfft=mfft)\n\nn2 = int(max_freq/SFT.delta_f)  # Number of bins below max_freq (30 Hz).\n# Dimensions of Sx, check the code in playing building spectrograms.\ndim1 = 39\ndim2 = 682\n# Each channel is appended to all.\n# all = np.array([]).reshape(0, dim1, dim2, n_channels)\nspecs = np.empty((len(df_test), dim1, dim2, n_channels))\n\n# item: [eeg_id, eeg_sub_id, idx in all (1st index), target,\n#       seizure_vote, lpd_vote, gpd_vote, lrda_vote,\n#       grda_vote, other_vote]\nitems = np.array([], dtype=float).reshape(0,10)\n\nfor i in np.arange(len(df_test)):\n    if i%400 == 0:\n        print(f'{i} files loaded.', end='\\r')\n    item = df_test.iloc[i]\n    eeg = pd.read_parquet(f'{base_dir}/train_eegs/{item.eeg_id}.parquet')\n    eeg = eeg.interpolate(limit_direction='both') # <<<<< Interpolation\n\n    # 10 second eeg sub samples \n    offset = int(item.eeg_label_offset_seconds)\n    start = (offset + 20) * fs\n    end = (offset + 30) * fs\n    eeg_sub_10 = eeg[start:end]\n\n    N = eeg_sub_10.shape[0]\n    t_x = np.arange(N) * 1/fs  # time indexes for signal\n\n    X = np.empty((1, dim1, dim2, n_channels))\n    # X[0,:,:,c] = Sx[n1:n2,:].copy()\n\n    # for c in np.arange(n_channels):\n    x = eeg_sub_10['F3'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,0] = Sx[1:n2,:].copy()\n    x = eeg_sub_10['F4'].values\n    Sx = SFT.spectrogram(x)  # calculate absolute square of STFT\n    specs[i,:,:,1] = Sx[1:n2,:].copy()\n\n    xitem = np.array([item.eeg_id, item.eeg_sub_id, i, item.target,\n                    item.seizure_vote, item.lpd_vote, item.gpd_vote,\n                    item.lrda_vote, item.grda_vote, item.other_vote],\n                    dtype=float).reshape(1,10)\n    items = np.concatenate([items, xitem])\n\nfilename = '03_single_spectrograms_v1_test'     \nprint(f'Saving to {filename}.npy')\nprint(f'Saving to {filename}_items.npy')\nnp.save(f'{output_dir}{filename}.npy', specs)\nnp.save(f'{output_dir}{filename}_items.npy', items)","metadata":{},"execution_count":null,"outputs":[]}]}