Wavelet transform-based mode decomposition for EEG signals under general anesthesia

Algorithm for the EWT and the Hilbert transform

In 2013, Gilles (2013) presented the EWT method, in which adaptive wavelet filter banks were built to decompose a given signal into IMFs. Like EMD, the EWT is intended to extract the IMF components from a signal. However, the EWT operates in the frequency domain, whereas EMD is performed in the time domain. Therefore, in the EWT, it is necessary to consider where the construction algorithm used for the filter banks draws boundaries in the frequency content within the spectrum of the signal under analysis (see Fig. 1 and Fig. S1). If we assume that the frequency domain [0, π] is divided into N consecutive segments, we must then extract N −1 boundaries, excluding the 0 and π boundaries. The local maxima in the spectrum are then determined and sorted in descending order to locate the boundaries for the EWT. Each boundary was defined as the average of the consecutive maxima. Let ω_n be the limit between each segments (where ω₀ =0 and ω_N = π) and let each segment be denoted by Λ_n = [ ω_n −1, ω_n]; then, ${\bigcup }_{n=1}^N{\Lambda }_n=$ [0, π]. A transition phase T_n (τ_n = γω_n, where 0 <γ <1) of width 2 τ_n is defined near the center of each Λ_n, as illustrated in (Figs S2A and S2B). An empirical wavelet is defined as a bandpass filter for each Λ_n (Figs S2C and S2S) (Daubechies, 1992). For this purpose, Gilles (2013) used the following steps when constructing the Meyer wavelets.

EWT steps

Step (1)

Apply a fast Fourier transform to the signal f (t)(t = {t_i}, i_=1,2,....,M, where M denotes the number of samples) to obtain the frequency spectrum X(ω). Then, find the set of maxima M ={M_i}_i=1,2,...,N in the Fourier spectrum by applying the convolution integral and estimate their corresponding frequencies ω = {ω_i}_i=1,2....,N, where N denotes both the number of maxima and the number of filter banks introduced from this point.

Step (2)

Perform appropriate segmentation of the Fourier spectrum acquired in Step 1 and define the boundaries Ω_i for each segment as the center points between two consecutive maxima: (1)${\displaystyle {\Omega }_i=\frac{{\omega }_i+{\omega }_{i+1}}{2}}$ where ω_i and ω_i+1 are both frequencies and the set of boundaries is Ω = {Ω_i}_i=1,2,N−1.

Step (3)

The empirical scaling function and the empirical wavelets with ∀_n >0 are defined in Eqs. (2) and (3), respectively, as shown in Fig. S2C, where ∀ denotes universal quantification and is read as “for all”. (2)${\displaystyle {\hat{\varphi }}_n\left(\omega \right)=\left\{\begin{array}{cc}\left|\omega \right|\le {\omega }_n-{\tau }_n\hfill & 1\hfill \\ {\omega }_n-{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n\hfill & cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)\right]\hfill \\ \text{otherwise}\hfill & 0\hfill \end{array}\right.}$ (3)${\displaystyle {\hat{\psi }}_n\left(\omega \right)=\left\{\begin{array}{cc}{\omega }_n+{\tau }_n\le \left|\omega \right|\le {\omega }_n-{\tau }_n\hfill & 1\hfill \\ {\omega }_{n+1}-{\tau }_{n+1}\le \left|\omega \right|\le {\omega }_{n+1}+{\tau }_{n+1}\hfill & cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_{n+1}}\left(\left|\omega \right|-{\omega }_{n+1}+{\tau }_{n+1}\right)\right)\right]\hfill \\ {\omega }_n-{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n\hfill & sin\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)\right]\hfill \\ \text{otherwise}\hfill & 0\hfill \end{array}\right.}$ where β(x) is any function C_k ([0,1]), such that: (4)${\displaystyle \beta \left(x\right)=\left\{\begin{array}{cc}x\le 0\hfill & 0\hfill \\ x\ge 1\hfill & 1\hfill \\ x\in \left(0,1\right)\hfill & \left(x\right)+\beta \left(1-x\right)=1\hfill \end{array}\right.}$ where (5)${\displaystyle \beta \left(\mathrm{x}\right)={\mathrm{x}}^4\left(35-84\mathrm{x}+70{\mathrm{x}}^2+20{\mathrm{x}}^3\right)}$ Regarding the selection of τ_n, τ_n is proportional to ω_n:τ_n = γω_n, where 0 <γ < 1. Consequently, when ∀_n >0, Eqs. (6) and (7) can be redefined as follows: (6)${\displaystyle {\hat{\varphi }}_n\left(\omega \right)=\left\{\begin{array}{cc}⌈\omega ⌉\le \left(1-\gamma \right){\omega }_n\hfill & 1\hfill \\ \left(1-\gamma \right){\omega }_n\le \left|\omega \right|\le \left(1+\gamma \right){\omega }_n\hfill & cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2\gamma {\omega }_n}\left(\left|\omega \right|-\left(1-\gamma \right){\omega }_n\right)\right)\right]\hfill \\ \text{otherwise}\hfill & 0\hfill \end{array}\right.}$ (7)${\displaystyle {\hat{\psi }}_{\mathbf{n}}\left(\omega \right)=\left\{\begin{array}{cc}\left(1+\gamma \right){\omega }_n\le \left|\omega \right|\le \left(1-\gamma \right){\omega }_{n+1}\hfill & 1\hfill \\ \left(1-\gamma \right){\omega }_{n+1}\le \left|\omega \right|\le \left(1+\gamma \right){\omega }_{n+1}\hfill & cos\left[\frac{\pi }{2}\beta \left(\frac{1}{2\gamma {\omega }_{n+1}}\left(\left|\omega \right|-\left(1-\gamma \right){\omega }_{n+1}\right)\right)\right]\hfill \end{array}\right.}$

Proposition 1: if $\gamma <mi{n}_n\left(\frac{{\omega }_{n+1}-{\omega }_n}{{\omega }_{n+1}+{\omega }_n}\right)$.

The concept for construction of the Meyer wavelet is based on the set $\left\{{\varphi }_1\left(t\right),{\left\{{\psi }_n\left(t\right)\right\}}_{n=1}^N\right\}$, which is a tight frame when (8)${\displaystyle \sum _{k=-\infty }^{+\infty }\left({\left|{\hat{\varphi }}_1\left(\omega +2k\pi \right)\right|}^2+\sum _{n=1}^N{\left|{\hat{\psi }}_n\left(\omega +2k\pi \right)\right|}^2\right)=1}$ as shown in Fig. S2D.

Step (4)

Empirical wavelet transform

The EWT $W_f^{\varepsilon }\in \left(n,t\right)$ of the signal f (t) is defined as the inner product of the signal to be analyzed f with the scaling function and the empirical wavelets. (9)${\displaystyle W_f^{\varepsilon }\left(0,t\right)=〈f,{\varphi }_1〉=\int f\left(\tau \right)\overline{{\varphi }_1\left(\tau -t\right)}d\tau }$ (10)${\displaystyle W_f^{\varepsilon }\left(n,t\right)=〈f,{\psi }_n〉=\int f\left(\tau \right)\overline{{\psi }_n\left(\tau -t\right)}d\tau }$ where $W_f^{\varepsilon }$(i, t) denotes the detailed coefficient of the i- th filter bank at a time t.

The total mode f (t) can then be reconstructed as follows: (12)${\displaystyle f\left(t\right)=W_f^{ℰ}\left(0,t\right)\ast ,{\varphi }_1\left(t\right)+\sum _{n=1}^N{W}_f^{ℰ}\left(n,t\right)\ast {\psi }_n\left(t\right)={\left(\widehat{W_f^{ℰ}}\left(0,\omega \right){\hat{\varphi }}_1\left(\omega \right)+\sum _{n=1}^N\widehat{W_f^{ℰ}}\left(n,\omega \right)\ast {\hat{\psi }}_n\left(\omega \right)\right)}^{\vee }}$ Here, ∧ denotes a Fourier transform and ∨ denotes an inverse Fourier transform.

Finally, the empirical mode f_k is defined as follows: (13)${\displaystyle f_0\left(t\right)=W_f^{ℰ}\left(0,t\right)\ast ,{\varphi }_1\left(t\right)}$ (14)${\displaystyle f_k\left(t\right)=W_f^{ℰ}\left(k,t\right)\ast ,{\psi }_k\left(t\right)}$

Algorithm for WMD

In the EWT, the boundary frequencies that separate the IMFs are determined as variable values for each analysis epoch by searching for local polar regions of power via application of an empirical method in the Fourier domain. (15)${\displaystyle EEG=\mathrm{IMF}{1}_{0.125\mathrm{Hz}\sim {\Omega }_1}+\mathrm{IMF}{2}_{{\Omega }_1\sim {\Omega }_2}+\mathrm{IMF}{3}_{{\Omega }_2\sim {\Omega }_3}+\mathrm{IMF}{4}_{{\Omega }_3\sim {\Omega }_4}+\mathrm{IMF}{5}_{{\Omega }_4\sim {\Omega }_5}+\mathrm{IMF}{6}_{{\Omega }_5\sim 64\mathrm{Hz}}}$ (where Ω_i denotes the boundary between the (i −1)th segment and the i-th segment)

However, if the boundary frequencies are determined in advance and the transformation performed by the Meyer wavelet is used as a bandpass filter, the EEG can be then decomposed into multiple IMF modes that fall within the range of the target frequencies. This method is a modified EWT method; however, we call it the WMD method because it discards the empirical approach and decomposes the process into modes with specific target frequencies via the Meyer wavelet transform (see Fig. 1), as indicated by the following equation: (16)${\displaystyle EEG=\mathrm{IMF}{1}_{0.125\sim 4\mathrm{Hz}}+\mathrm{IMF}{2}_{4\sim 8\mathrm{Hz}}+\mathrm{IMF}{3}_{8\sim 14\mathrm{Hz}}+\mathrm{IMF}{4}_{14\sim 20\mathrm{Hz}}+\mathrm{IMF}{5}_{20\sim 30\mathrm{Hz}}+\mathrm{IMF}{6}_{30\sim 64\mathrm{Hz}}.}$ To set the preset boundaries Ω₁₋₅ required in Step 2 of the EWT, the EWT algorithm was modified to decompose six IMFs in the WMD. Ω₁, Ω₂, Ω₃, Ω₄, and Ω₅ were set at 4 Hz, 8 Hz, 14 Hz, 20 Hz, and 30 Hz, respectively.

Anesthesia management and data acquisition

As reported previously (Yamada et al., 2023), the original EEG data were acquired according to the principles of the Declaration of Helsinki. The Institutional Review Board at the Kyoto Prefectural University of Medicine (IRB of KPUM) approved the current study (No. ERB-C-1074) for human experiments. The requirement for informed patient consent was waived by the IRB of the KPUM for this noninterventional and noninvasive retrospective observational study (all patients were provided with an opt-out option, of which they were notified in the preoperative anesthesia clinic). We analyzed the intraoperative EEG data (the sampling frequency is 128 Hz and the data length is 8 s, therefore, the number of data points per epoch is 1,024) acquired from 10 patients that were maintained under sevoflurane anesthesia. Data acquisition was performed via using EEG Analyzer software (ver. 54_GP; https://anesth-kpum.org/blog_ts/?p=3169) (Hayase et al., 2021) linked to a VISTA A-3000 BIS monitor with a BIS Quatro sensor. For further details about patient information, anesthesia management, and the data acquisition procedure, including raw EEG data, the BIS index, the 95% spectral edge frequency (SEF₉₅), the total power (TP), and the absolute power derived from electromyography in the 70–110 Hz range (EMG _low), please see our previous report (Yamada et al., 2023). This analysis used the sevoflurane GA dataset (Supplemental File S3) that we reported previously (Yamada et al., 2023).

EEG Mode Decomposition for EMD, VMD, EWT, and WMD

The EMD analysis software EEG Mode Decompositor (ver. 2.7N, https://anesth-kpum.org/blog_ts/?p=4114) was initially developed to perform EMD (Obata et al., 2023) and was subsequently updated by adding a VMD mode (Yamada et al., 2023) and, in this case, EWT and WMD modes (Gilles, 2013). Supplemental File S4 (Python programming code for EWT) shows the Python code that was described in the report by Carvalho et al. (2020) and used to perform EWT analysis. Supplemental File S5 (Processing programming code for EWT and WMD) shows the EWT and WMD classes of the Processing and Java programming codes, respectively. We intended to observe significant changes in specific narrow frequency bands within the 0–64 Hz range, such as δ (0.0625<f ≤4 Hz), θ (4<f ≤8 Hz), α (8<f ≤14 Hz), low- β (14<f ≤20 Hz), high- β (20<f ≤30 Hz), and γ waves (30<f ≤64 Hz).

Data processing and statistics

The changes in the various EEG parameters between two time points were compared statistically by performing Kruskal–Wallis tests with SPSS Statistics (ver. 28.0.1.0, IBM Corp., Armonk, NY, USA), and p values <0.05 were considered statistically significant. Multiple linear regression (MLR) analysis involving between the BIS index and the parameters of the IMFs was performed via the scikit-learn and StatsModel libraries in Python (ver 3.8), as we reported recently (Sawa, Yamada & Obata, 2022).

Case analysis of the EWT application

First, we present the VMD, EWT, and WMD analyses that were conducted using the EEG data acquired from Patient #6 (a 27-year-old female with an ovarian tumor undergoing laparoscopic oophorectomy) (Yamada et al., 2023) when anesthetized with sevoflurane (video files for a 30 min period until emergence: Supplemental Files S6–S9: Sev6_EMD.mp4, Sev6_VMD.mp4, Sev6_EWT.mp4, and Sev6_WMD.mp4, respectively). For a comparative analysis with the previously reported VMD adaptation (see Figs. 3, and 4 in Yamada et al. 2023), we analyzed 4-s-long epochs of EEGs that were recorded at 128 Hz during three distinct phases of the GA up to the point of emergence: the maintenance (30 min before emergence), transition (1 min before emergence), and emergence phases. Figure 2 presents the raw EEGs (4 s), the IMFs obtained over 4 s via the three decomposition methods, and the density spectral arrays (DSAs) (64 s, starting with a 4-s-long segment of raw EEG data) across these GA phases. Given the challenges involved in visually discerning the differences among the characteristics of the IMFs from such waveform diagrams, a Hilbert transform was applied to derive both the instantaneous frequency (IF) and instantaneous amplitude (IA). The results are shown in the Hilbert spectrum in Fig. 3. During the maintenance phase of anesthesia, IMFs 1 to 4 were analyzed using the VMD method and were predominantly in the frequency region around 20 Hz or lower. In contrast, the six IMFs spanned a broader frequency range from 0 Hz up to 50–60 Hz in both the EWT and WMD methods. During the transition period, IMF-6 notably shifted noticeably to the 30-Hz frequency during the VMD analysis, whereas the peaks of IMF-5 and IMF-6 in the EWT analysis moved toward higher frequencies above 40 Hz. Upon emergence, the peaks of IMF-5 and IMF-6 in the EWT increased to frequencies above 50 Hz, whereas they remained below 40 Hz in the VMD results. The WMD method displays frequency histograms consistently for each IMF within a narrow band corresponding to predefined target frequencies regardless of the GA phase considered. This distinction arises because the EWT and VMD do not impose frequency band constraints on the decomposition of IMFs; in contrast, WMD requires decomposition into IMFs within specified frequency bands.

Hilbert spectrograms of the IMFs in the EWT and WMD methods

To gain further understanding of the frequency characteristics of each IMF produced by VMD, the EWT, and WMD across the 30 min phase before emergence, a Hilbert spectral analysis of the IMFs was performed. This analysis assessed the variations in the IF and IA values derived from the Hilbert transforms of the IMFs during each phase for all ten patients. Fig. S3 shows the Hilbert spectrograms for each IMF obtained via the VMD, EWT, and WMD methods, the summed IMF (IMF-all, which corresponds to the original EEG), and the DSAs for 30 min before emergence for all patients. In the EWT method, IMF-1, IMF-2, and IMF-3 effectively capture the δ, θ, and α frequency bands, respectively. After emergence, however, a notable reduction in the power of IMF-3, which primarily extracts the α-band, was observed. The central frequencies of IMF-4, IMF-5, and IMF-6 span a broad spectrum within the β- and γ-bands. Conversely, in the WMD approach, IMF-1, IMF-2, and IMF-3 effectively extract the δ, θ, and α bands, respectively, and this is consistent with their preset boundaries. In contrast to the EWT results, IMF-4, IMF-5, and IMF-6 are characterized by narrow bandwidths within specific frequency ranges, i.e., low β, high β, and γ, respectively. Notably, the WMD method demonstrated more stable frequency band retention across the separated IMFs than either of the EWT and VMD methods.

Multiple linear regression models to predict the BIS

Next, we analyzed the processed EEG data, which were sampled every 8 s, and were acquired from 10 patients who had been administered a sevoflurane GA 30 min before emergence. We explored the relationship between the BIS index (a clinical gauge of the depth of anesthesia) and the IMF parameters, including the center frequency (i.e., the average of the Hilbert spectrum) and TP, which were converted from the power in P (µV²) to decibels (dB), using dB = 10 ⋅ log₁₀P). MLR models were developed using the IMF parameters to predict the BIS index based on the Data_EME30s dataset (Supplemental File S11). These models incorporated the center frequencies and TP values of the IMFs as variables and were tested using three modal decomposition methods: VMD, EWT, and WMD. The comprehensive MLR formula included individual terms for the center frequencies and the TPs of the six IMFs, thus leading to a model with 36 explanatory variables. Figure 4 shows the time-trend graphs for the BIS index, the SEF₉₅, the TP, EMG_low, and the corresponding IMFs and TPs across the decomposition methods. The analysis yielded robust correlations, with the WMD-based MLR model using all 36 parameters (18 center frequencies and 18 TPs) and demonstrating the strongest correlation (R ² = 0.899) and the lowest errors (mean absolute error: MAE = 0.244; root mean squared error: RMSE = 0.318) among the four models examined; the results are detailed in Table 1. A simplified WMD model that used only 12 parameters ranked second. For a more targeted approach, we evaluated an MLR model using only seven IMF parameters (IMF-1_freq, IMF-3_freq, IMF-5_freq, IMF-6_freq, IMF-2_TP, IMF-4_TP, and IMF-6_TP) that presented significant p values (less than 0.05) but still presented excellent fit metrics (R ² =0.897, MAE = 0.244, RMSE = 0.321), as shown in Fig. 5 and Table S1. Further analyses were conducted by using separate models for the central frequencies and TPs, and the results confirmed that combining the two parameters offered superior predictive accuracy, as shown in Fig. S4 and Table S2. Notably, in the WMD approach, the center frequencies of IMF-1, IMF-3, IMF-5, and IMF-6, and the TPs of IMF-2, IMF-4, and IMF-6, were identified as significant predictors. This finding indicates that WMD effectively captures the narrow band changes in the IMF parameters correlated with the shifts in the BIS index that occur during emergence from GA, thus indicating the superior sensitivity of WMD for monitoring changes in depth of anesthesia.