{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"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":"gpu","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"I noticed that in clinical setting, the Bispectrogram is used. I wondered if this might create a more visual interpretation of the situation. In early AI approaches various parameters from the bispectrum are used to train systems. But deep learning systems are superior to these approaches. I wondered if using a visual deep learning system on the image output from bispectrum might be useful. So far the results are inconclusive. But it may be that other aspects of the bispectrum are useful. \n\nThere are many articles on the use of bispectrum in both monitoring and analysis of epilepsy. There appear to be no standard libraries for calculating bispectrum. \n\nMy plan is to create an alternative to spectrograms and see if this gives superior results. Note that the images here are a different size to the competition spectrograms, so some small changes will be needed.\n\nI have of course used the excellent Chris notebooks here. ","metadata":{}},{"cell_type":"markdown","source":"Update 6/3/24 : I tried to classify the images produced here, but the results were quite poor. So I have abandoned this direction. ","metadata":{}},{"cell_type":"markdown","source":"**Important note**  I have sampled a small time slice in the center of the eeg waveform, in order to get results in a reasonable time. If you have access to a powerful machine you could enlarge this window. From my visual inspection I cannot see any difference in the images created from the whole waveform. ","metadata":{}},{"cell_type":"markdown","source":"bispectrum2D is from https://gist.github.com/jkmackie  - I have modified the code to make use of GPU. I found the documentation on how to do this appalling. How can this company be worth trillions of dollars? ","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom matplotlib import pyplot as plt\n\nimport librosa\nimport os,gc\n\n\n","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:29.980965Z","iopub.execute_input":"2024-03-04T05:58:29.981850Z","iopub.status.idle":"2024-03-04T05:58:29.986384Z","shell.execute_reply.started":"2024-03-04T05:58:29.981806Z","shell.execute_reply":"2024-03-04T05:58:29.985463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"USING_GPU = True\n\nif USING_GPU:\n    import cupy as cp","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:29.988171Z","iopub.execute_input":"2024-03-04T05:58:29.989117Z","iopub.status.idle":"2024-03-04T05:58:29.998402Z","shell.execute_reply.started":"2024-03-04T05:58:29.989078Z","shell.execute_reply":"2024-03-04T05:58:29.997231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from dataclasses import dataclass, field\nimport numpy as np\nfrom scipy.fft import fftshift, fft2, ifftshift, fft\nfrom scipy.linalg import toeplitz\nimport matplotlib.pyplot as plt\n\nif USING_GPU:\n    from cupy.fft import fftshift, fft2, ifftshift, fft\n","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:29.999750Z","iopub.execute_input":"2024-03-04T05:58:30.000519Z","iopub.status.idle":"2024-03-04T05:58:30.010269Z","shell.execute_reply.started":"2024-03-04T05:58:30.000493Z","shell.execute_reply":"2024-03-04T05:58:30.009465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n@dataclass(repr=False)\nclass Bispectrum2D:\n    '''\n    Make Bispectrum dataclass from 1D signal with frequency in Hertz.\n    \n    Parameters:\n    -----------\n        signal: The 1-Dimensional signal in (n,) numpy.ndarray\n        freqsample:  The sample frequency in Hertz.  \n        window: 'None', 'hanning', 'triangular'\n    \n    References\n    ----------\n    Matteo Bachetti, et al., stingray v1.0 code, DOI: https://zenodo.org/record/6394742, 2022.\n\n    Plot bispectrum magnitude as follows: bs2D.plot_bispec_magnitude()\n    \n    '''\n    signal: np.ndarray\n    freqsample: float\n    window_name: str\n    dt: float = field(init=False)\n    n: int = field(init=False)  # number of data points in signal\n    maxlag: int = field(init=False)\n    lagindex: np.ndarray = field(init=False)\n    cum3_dim: int = field(init=False)\n\n    def __post_init__(self):\n        self.dt = 1 / self.freqsample\n        self.n = self.signal.shape[0]\n        self.maxlag = int(self.n / 2)\n        self.lagindex = np.arange(-self.maxlag, self.maxlag + 1)\n        self.cum3_dim = 2 * self.maxlag + 1\n        self._calc_bispectrum()\n\n    def _calc_bispectrum(self):\n        self._cumulant3()\n        self._window()\n        self._bispectrum()\n\n    def _cumulant3(self):\n        '''Biased cumulant estimate.'''\n        self.cum3 = np.zeros((self.cum3_dim, self.cum3_dim))  # include zeros matrix to reset calc\n        ind = np.arange((self.n - self.maxlag) - 1, self.n)  # consecutive idx from (n-maxlag-1) to n\n        ind_t = np.arange(self.maxlag, self.n)\n        zero_maxlag = np.zeros((1, self.maxlag))\n        zero_maxlag_t = zero_maxlag.T\n        sig = np.reshape(self.signal, (1, len(self.signal)))  # Reshape original self.sig\n        sig = sig - np.mean(sig)  # sig is 1xn row vector of counts.\n        rev_sig = np.array([sig[0][::-1]])\n        col = np.concatenate((sig.T[ind], zero_maxlag_t), axis=0)\n        row = np.concatenate((rev_sig[0][ind_t], zero_maxlag[0]), axis=0)\n        toep = toeplitz(col, row)\n        rev_sig_repeat = np.repeat(rev_sig, [2 * self.maxlag + 1], axis=0)  # n repeats\n        # toep is n x (n-1).  It must be square to be a circulant.\n        self.cum3 = (self.cum3 + np.matmul(np.multiply(toep, rev_sig_repeat), toep.T)) / self.n\n\n    def _window(self):\n        n = np.arange(self.cum3_dim)  # Total wind data points matches cum3_dim\n        self.window = np.zeros(self.cum3_dim)\n        if self.window_name == 'None':\n            return\n        if self.window_name == 'hanning':\n            hanning = 0.5 * (1 - np.cos(2 * np.pi * n / (self.cum3_dim - 1)))\n            wind2D = np.tile(hanning, (self.cum3_dim, 1))  # Make 2D wind by repeating rows N times.\n            self.window[:self.maxlag + 1] = hanning[self.maxlag:]\n        if self.window_name == 'triangular':\n            N_div_2 = int((np.floor((self.cum3_dim - 1) / 2)))\n            triangular = 1 - np.abs((n - (N_div_2)) / self.cum3_dim)\n            wind2D = np.tile(triangular, (self.cum3_dim, 1))  # Make 2D wind by repeating rows N times.\n            self.window[:self.maxlag + 1] = triangular[self.maxlag:]\n        self.window[self.maxlag:] = 0\n        # Put wind in toeplitz.  Each row of final wind is sliding hanning.\n        row = np.concatenate(([self.window[0]], np.zeros(2 * self.maxlag)))\n        toep_matrix = toeplitz(self.window, row)\n        toep_matrix += np.tril(toep_matrix, -1).transpose()\n        self.window = toep_matrix[..., ::-1] * wind2D * wind2D.T\n\n    def _bispectrum(self):\n        if USING_GPU:\n            xg = cp.array(self.cum3)\n            xwg = cp.array(self.cum3*self.window)\n            if self.window_name == 'None':\n                bispecg = cp.fft.fftshift(cp.fft.fft2(cp.fft.ifftshift(xg)))\n                self.bispec = cp.asnumpy(bispecg)\n            else:\n                bispecg = cp.fftshift(cp.fft2(cp.ifftshift(xwg)))\n                self.bispec = cp.asnumpy(bispecg)\n        else:\n            if self.window_name == 'None':\n                self.bispec = fftshift(fft2(ifftshift(self.cum3)))\n            else:\n                self.bispec = fftshift(fft2(ifftshift(self.cum3*self.window)))    \n\n        self.freqvals = 0.5 * self.freqsample * self.lagindex / self.maxlag\n        self.bispec_mag = np.abs(self.bispec)\n        self.bispec_phase = np.angle(self.bispec)\n\n    def plot_cum3(self):\n        lags = self.lagindex * self.dt\n        fig, ax1 = plt.subplots(1, 1, figsize=(6, 6))  # gist has (11,11)\n        contplot1 = ax1.contourf(lags, lags, self.cum3, levels=100, cmap=plt.cm.Spectral_r)\n        ax1.set_title('Third Order Cumulant');\n        ax1.set_xlabel('lag 1 values');\n        ax1.set_ylabel('lags 2 values')\n\n    def plot_bispec_magnitude(self):\n        fig, ax1 = plt.subplots(1, 1, figsize=(6, 6))\n        contplot1 = ax1.contourf(self.freqvals, self.freqvals, self.bispec_mag, levels=100, cmap=plt.cm.Spectral_r)\n        ax1.set_title('Signal Bispectrum Magnitude');\n        ax1.set_xlabel('freq 1 Hz');\n        ax1.set_ylabel('freq 2 Hz')\n        plt.colorbar(contplot1)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:30.011200Z","iopub.execute_input":"2024-03-04T05:58:30.011440Z","iopub.status.idle":"2024-03-04T05:58:30.039189Z","shell.execute_reply.started":"2024-03-04T05:58:30.011419Z","shell.execute_reply":"2024-03-04T05:58:30.038210Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:30.041076Z","iopub.execute_input":"2024-03-04T05:58:30.041329Z","iopub.status.idle":"2024-03-04T05:58:30.203846Z","shell.execute_reply.started":"2024-03-04T05:58:30.041308Z","shell.execute_reply":"2024-03-04T05:58:30.202969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NAMES = ['LL','LP','RP','RR']\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\n","metadata":{"collapsed":false,"ExecuteTime":{"end_time":"2024-02-26T01:08:08.780218Z","start_time":"2024-02-26T01:08:08.773590Z"},"jupyter":{"outputs_hidden":false},"execution":{"iopub.status.busy":"2024-03-04T05:58:30.205310Z","iopub.execute_input":"2024-03-04T05:58:30.205613Z","iopub.status.idle":"2024-03-04T05:58:30.210859Z","shell.execute_reply.started":"2024-03-04T05:58:30.205588Z","shell.execute_reply":"2024-03-04T05:58:30.209795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is where I slice the waveform and only look at the centre. I did this because it was taking too long to calculate. If you have access to a powerful\nmachine you could enlarge this window. \n","metadata":{}},{"cell_type":"code","source":"def mid_slice(x,size):\n    mid = x.shape[0]//2\n    low = mid-size//2\n    high = mid+size//2\n    return x[low:high]\n    ","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:30.211981Z","iopub.execute_input":"2024-03-04T05:58:30.212308Z","iopub.status.idle":"2024-03-04T05:58:30.220854Z","shell.execute_reply.started":"2024-03-04T05:58:30.212277Z","shell.execute_reply":"2024-03-04T05:58:30.220069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def spectrogram_from_eeg(parquet_path, chunk_size, 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:middle+10_000]\n    \n    # VARIABLE TO HOLD SPECTROGRAM\n    img = np.zeros((64,64,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        \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            signals.append(x)\n            \n            eeg_length = x.shape[0]\n            x = mid_slice(x,chunk_size)\n            eeg_length = x.shape[0]\n    \n            # RAW BI-SPECTROGRAM    \n            bi_spec = Bispectrum2D(signal=x, freqsample=100, window_name='None')\n            \n            low_index = eeg_length//2 - 32\n            high_index = eeg_length//2 + 32\n            \n            freq = bi_spec.freqvals[low_index:high_index]         \n            bi_spec_mag = bi_spec.bispec_mag[low_index:high_index,low_index:high_index]\n            \n            img[:,:,k] += bi_spec_mag\n                \n        # AVERAGE THE 4 MONTAGE DIFFERENCES\n        img[:,:,k] /= 4.0\n        \n        if display:\n            plt.subplot(2,2,k+1)\n            plt.imshow(img[:,:,k],aspect='auto',origin='lower')\n            plt.title(f'EEG {eeg_id} - Bi-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(10_000),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":{"execution":{"iopub.status.busy":"2024-03-04T05:58:30.236670Z","iopub.execute_input":"2024-03-04T05:58:30.237280Z","iopub.status.idle":"2024-03-04T05:58:30.250432Z","shell.execute_reply.started":"2024-03-04T05:58:30.237257Z","shell.execute_reply":"2024-03-04T05:58:30.249532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nPATH = '/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/'\nDISPLAY = 4\nEEG_IDS = train.eeg_id.unique()\nall_eegs = {}\nchunk = 500 \n\nimage_save = '/kaggle/working/bi_spectrograms'\nif not os.path.exists(image_save):\n    os.makedirs(image_save)\n\nsample_eeg = EEG_IDS\n\nfor i,eeg_id in enumerate(sample_eeg):\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', chunk,  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        \n    filename = image_save + '/' + str(eeg_id)\n    np.save(filename,img)\n    all_eegs[eeg_id] = img\n    \n    gc.collect()\n   \n","metadata":{"execution":{"iopub.status.busy":"2024-03-04T05:58:30.252244Z","iopub.execute_input":"2024-03-04T05:58:30.252836Z","iopub.status.idle":"2024-03-04T05:58:39.027445Z","shell.execute_reply.started":"2024-03-04T05:58:30.252805Z","shell.execute_reply":"2024-03-04T05:58:39.026528Z"},"trusted":true},"execution_count":null,"outputs":[]}]}