2.1.

Basic Formulas of the MI-Based NIGM Method

Figure 1 shows a schematic diagram of a near-infrared spectroscopy (NIRS) measurement on a living body. According to the modified Beer-Lambert law (MBLL), by using two different probe wavelengths and solving for the matrix calculation, changes in oxy- and deoxyhemoglobin NIRS signals $\mathrm{\Delta }{N}_{{\mathrm{HbO}}_2}(t)$ and $\mathrm{\Delta }{N}_{\mathrm{Hb}}(t)$ can be expressed as follows:^20–24

Fig. 1

Schematic of NIRS measurement.

Here, both oxy- and deoxyhemoglobin incremental NIRS signals $\mathrm{\Delta }{N}_{{\mathrm{HbO}}_2}(t)$ and $\mathrm{\Delta }{N}_{\mathrm{Hb}}(t)$ can be decomposed into quasi-DC and AC components as follows:

Here, the subscripts AC and DC indicate the AC and quasi-DC components of the corresponding physical quantity, respectively. For reference, a visual explanation of this decomposition operation is shown in Fig. 2.

Fig. 2

Conceptual diagram of NIRS signal decomposition into quasi-DC and AC components.

Here, by applying several assumptions and approximations, the oscillation amplitude of the AC component of deoxyhemoglobin concentration $C_{\mathrm{Hb},\mathrm{AC}}(t)$ is obtained as follows:⁹

Fig. 3

Metabolic index MI behavior versus oxygen saturation.

Here, by applying the assumption about the relationship between $L(t)$ and $L_{\mathrm{AC},\text{amplitude}}(t)$, Eq. (7) can be rewritten as follows:⁹

Here, $L_{\mathrm{AC},\text{amplitude},0}$ and $L_{\mathrm{DC},0}$ are AC and quasi-DC components of the optical path length in the initial condition, respectively. The term $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ denotes an amplitude-corrected metabolic index, and $\alpha (t)$ denotes a dimensionless optical path length correction factor, expressed as follows:⁹

Here, $L_{\mathrm{AC}}(t)$ can be obtained from $\mathrm{\Delta }{N}_{{\mathrm{HbO}}_2,\mathrm{AC},\text{amplitude}}(t)$ and $\mathrm{\Delta }{N}_{\mathrm{Hb},\mathrm{AC},\text{amplitude}}(t)$ by employing the following relationship:

As $C_{\mathrm{Hb},\mathrm{AC}}(t)$ corresponds to the oxygen consumption in each cardiac cycle of the probed region and BGL affects the living body metabolism, $C_{\mathrm{Hb},\mathrm{AC}}(t)$ can be assumed to exhibit a strong correlation with BGL. Furthermore, given $C_0$, $L_{\mathrm{AC},\text{amplitude},0}$, and $L_{\mathrm{DC},0}$ are all constant, BGL can be estimated by monitoring $\mathrm{M}{\mathrm{I}}^{\prime }(t)$.

2.2.

Derivation of the Signal Quality Index for the MI Method

To evaluate the signal quality for the AC oxy- and deoxy hemoglobin NIRS signals $\mathrm{\Delta }{N}_{{\mathrm{HbO}}_2,\mathrm{AC}}(t)$ and $\mathrm{\Delta }{N}_{\mathrm{Hb},\mathrm{AC}}(t)$, the standard deviation concept is introduced. Figure 4(a) shows a typical example of the discretely-sampled oxy- and deoxyhemoglobin NIRS signals, and Fig. 4(b) illustrates the band-pass-filter (BPF)-applied and normalized waveforms of Fig. 4(a), where components around the HR frequency are extracted.

Fig. 4

(a) Typical example of the sampled oxy- and deoxyhemoglobin waveforms. (b) Normalized BPF-applied waveforms of panel (a).

As can be seen in Fig. 4(b), the two waveforms exhibit a high degree of similarity with minimal discernible differences in ideal conditions. In this study, the authors sought to quantify the dissimilarity between the two normalized BPF-applied waveforms by employing the standard deviation $s(t)$, which is defined as

Here, consider the scenario where noise is added to the NIRS waveforms. Figure 5(a) shows an example of the discretely-sampled oxy- and deoxyhemoglobin signals with noise, and Fig. 5(b) illustrates the normalized BPF-applied waveforms of Fig. 5(a), where components around the HR frequency are extracted. In this example, a small spike is superimposed on the deoxyhemoglobin waveform around data number 200. Figure 5(b) illustrates that a minor fluctuation in the NIRS waveform can have a significant impact on the BPF-applied waveform, resulting in an increase in the standard deviation $s(t)$.

Fig. 5

(a) Example of the sampled oxy- and deoxyhemoglobin waveforms with noise. (b) Normalized BPF-applied waveforms of (a).

Nevertheless, even under optimal conditions, the standard deviation $s(t)$ is constrained by the inherent delay between the oxy- and deoxyhemoglobin NIRS waveforms, which is on the order of several tens of milliradians. Given a specific value for the phase delay $\mathrm{\Delta }\theta (t)$, the reachable standard deviation limit $\sigma (t)$ can be expressed as follows:

Here, the derivation process of Eq. (16) is described in detail in Sec. S1 of the Supplementary Material. Subsequently, by applying a deformation to Eq. (16), the following is derived:

In this study, each calculated $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ is initially screened by ${\epsilon }_{\sigma }(t)$ to be less than a predetermined value.

2.3.

Derivation of Theoretical Phase Estimation Errors Defined by the Background Noise Level

Although ${\epsilon }_{\sigma }(t)$ in Eq. (18) is effective in identifying and rejecting visually recognizable errors in NIRS signals, it lacks the capacity to assess potential errors that are not visually apparent. In a strict sense, the FFT-based phase delay, $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$, is inherently affected by perturbation from the background noise, which is typically defined according to the noise floor level. This results in an error in phase detection. In this subsection, the authors attempt to formulate the phase estimation error defined by SNR.

Figure 6 illustrates the basic idea of the phase estimation error caused by the random background noise. Here, $\overrightarrow{a}$ represents the normalized main peak vector of the FFT-applied NIRS signal, and $\overrightarrow{b}$ represents the background noise vector with amplitude and argument of $1/\mathrm{SNR}(t)$ and $\varphi$, respectively. In this representation, the phase estimation error, which is denoted by ${\psi }_{\mathrm{SNR},\varphi }$, can be replaced by the argument of the combined vector $\overrightarrow{a}+\overrightarrow{b}$ and can be expressed as follows:

Fig. 6

Basic concept of the phase estimation error caused by the background noise.

Then, the typical phase estimation error due to the background noise, denoted by ${\epsilon }_{\mathrm{SNR}}(t)$, can be calculated as the square root of the mean value of ${\psi }_{\mathrm{SNR},\varphi }^2$ over $-\pi \le \varphi \le \pi$, which is expressed as follows:

Here, Eq. (20) cannot be solved analytically. However, by assuming that $\nu$, the inverse value of the SNR, is sufficiently less than 1, Eq. (20) can be approximated as

The derivation process of Eq. (21) is described in detail in Sec. S2 of the Supplementary Material. Here, $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ is practically calculated by comparing the FFT phase of the oxy- and deoxyhemoglobin signals at the HR frequency. Therefore, the total phase estimation error, denoted by ${\epsilon }_{\mathrm{SNR}(t),\text{total}}$, should be expressed as the resultant error as follows:

Typically, the signal amplitude of the deoxyhemoglobin NIRS signal is approximately one-tenth of that of the oxyhemoglobin NIRS signal. Assuming that the background noise amplitude is almost identical between oxy- and deoxyhemoglobin NIRS signals, ${\mathrm{SNR}}_{\mathrm{Hb}}(t)$ is $\sim 10$ times worse than ${\mathrm{SNR}}_{{\mathrm{HbO}}_2}(t)$. Therefore, ${\epsilon }_{\mathrm{SNR},\mathrm{Hb}}(t)$ is dominant in Eq. (22), and Eq. (23) can be approximated as follows:

As indicated by Eq. (24), ${\epsilon }_{\mathrm{SNR},\text{total}}(t)$ becomes as high as 70 mrad when ${\mathrm{SNR}}_{\mathrm{Hb}}(t)$ is 10, or 20 dB in decibels. In the previous research, $40\mathrm{mg}/\mathrm{dL}$ of BGL change induced $\sim 50\mathrm{mrad}$ of phase delay in $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$.⁹ Consequently, the suboptimal quality of ${\mathrm{SNR}}_{\mathrm{Hb}}(t)$ may readily compromise the precision of the BGL estimation.

2.4.

Derivation of Theoretical Phase Estimation Errors Defined by the Data Sampling Frequency

It is also crucial to consider the sampling frequency of the NIRS signals to estimate the phase estimation error. When the sampling frequency and the HR frequency are defined as $f_s$ and $f_{\mathrm{HR}}(t)$, respectively, the mathematical phase division step size of the discretely sampled NIRS signals, denoted by ${\theta }_{\mathrm{div}}(t)$, can be expressed as follows:

Here, Figs. 7(a) and 7(b) illustrate examples of the potential range of temporal uncertainty resulting from discrete sampling, where Fig. 7(b) is an enlarged view of Fig. 7(a) at $t=0$. For the sake of simplicity, the values of $f_{\mathrm{s}}$ and $f_{\mathrm{HR}}(t)$ are set to 100 and 1.0 Hz in Figs. 7(a) and 7(b). In addition, the term ${\theta }_{\mathrm{OS}}$ represents the phase offset from the original waveform, the dashed vertical lines in Fig. 7(b) represent the horizontal sampling step size determined by $f_{\mathrm{s}}$ and $f_{\mathrm{HR}}(t)$, and the dash-dot lines in Fig. 7(b) represent a waveform reconstructed from the discretely sampled data.

Fig. 7

(a) Basic concept of the phase estimation error resulting from discrete sampling. (b) An enlarged view of panel (a) at $t=0$.

Given that the temporal phase uncertainty resulting from discrete sampling is constrained by the boundaries of $\pm {\theta }_{\mathrm{div}}(t)/2$, its probability can be approximated by the continuous uniform distribution. Subsequently, by applying the formula for the standard deviation of the continuous uniform distribution, the typical phase uncertainty due to discrete sampling ${\epsilon }_{\mathrm{s}}(t)$ can be expressed as follows:

As with the case of ${\epsilon }_{\mathrm{SNR}}$, both the oxy- and deoxyhemoglobin NIRS signals have their own sampling uncertainty. Consequently, the total sampling uncertainty, denoted by ${\epsilon }_{\mathrm{s},\text{total}}(t)$, can be given as follows:

2.5.

Total Estimation Error of the Metabolic Index Resulting from Phase Estimation Errors

From Eqs. (8), (10), (11), (24), and (27), the combined value of the latent estimation errors of the metabolic index $\mathrm{M}{\mathrm{I}}^{\prime }(t)$, denoted by $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$, can be approximated as follows:

In this study, the value of $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ screened by ${\epsilon }_{\sigma }(t)$ in Eq. (18) is subjected to further examination based on $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$, to ensure the estimation quality.

3.1.

Smartwatch-Based Prototype

In this study, the Samsung Galaxy Watch 4 44 mm (SM-R870) was utilized as the experimental unit in a manner consistent with that employed in the previous research.⁹ Figure 8(a) shows the schematic rear view of the SM-R870. Red and infrared (IR) LEDs are located in the center, and their center wavelengths are $\sim 650$ and 930 nm, respectively. Eight photodetectors are arranged radially around the LEDs. The sampling frequency of the red and IR LED signals is 100 Hz. Here, no hardware modifications were made to the prototype, and the BGL estimation function was implemented in software. Figure 8(b) shows a schematic of the application style of the prototype device. Here, the device was wrapped around a finger pad by taking account of the capillary density of the measurement portion, which correlates strongly with signal quality.

Fig. 8

(a) Schematic rear view of the smartwatch experimental unit. (b) A schematic diagram of the application style of the smartwatch prototype on a fingertip.

3.1.1.

Experimental protocol

Figure 9(a) shows the schematic of the clinical test setup. A protocol was established while referencing factors that can affect reading values in the case of a pulse oximeter.²⁵ To ensure consistency in the acquired data, a single, healthy, and non-diabetic male was selected as the sole test subject. The test subject was asked to sit still during the experiment with the smartwatch resting on their finger pad to reduce fluctuations in peripheral blood flow caused by minor changes in posture and local blood pressure. The subject was also asked to rest his elbows on armrests. In addition, before starting the measurement, the subject held a hand warmer for a few minutes to ensure adequate peripheral blood flow.

Fig. 9

(a) Schematic diagram of the clinical trial setup. (b) Clinical trial procedure.

Figure 9(b) shows the typical course of a clinical test. A CGM device (Abbot, FreeStyle Libre^®), which records BGL values at 15-min intervals, was used as a reference for the BGLs. To use CGM sensors at their best performance, CGM sensors were installed at least two days before the experiments for aging.²⁶ Furthermore, the sensors were not utilized for clinical testing during the final two days of their specified expiration period. In addition, to correct for the systematic offset and delay specific to each CGM sensor, SMBG strips (Abbot, FreeStyle Precision^®, Chicago, Illinois, United States) were used as needed before and after data recording by the smartwatch. To prevent possible distortion of the PPG signal due to body movements resulting from SMBG usage, no SMBG strips were used during the recording period. Data recording then started after at least 2 h of fasting, oral challenges were given $\sim 20$ to 30 min after the start of recording, and data recording continued $\sim 60\min$ after the oral challenges until the subject’s BGL had mostly returned to the initial level. As a result, $\sim 90\min$ of PPG data were collected in each data recording. The clinical tests were conducted over a period of four months, with a total of 30 repeatability tests performed.

In this research, sugar-containing carbonated beverages and glucose-containing jelly beverages were used for oral challenges. For reference, Table S1 in Sec. S3 of the Supplementary Material shows the main nutritional values of the oral challenges. When administering oral challenges, care was taken to ensure that all oral challenges were not excessively cold, thus avoiding peripheral blood flow reduction due to the lowering of the body temperature.

All clinical trials described in this paper were conducted in accordance with the Clinical Trials Act of the Ministry of Health, Labor and Welfare of Japan, published on the basis of the Declaration of Helsinki, and were approved by the Ethical Committee of Hamamatsu Photonics K. K. Informed consent was obtained from the subject before measurements were performed. All clinical tests were performed under the supervision of the co-author having a medical license.

3.1.2.

Data processing

Figure 10 depicts the flowchart for the calculation and data screening process of the modified metabolic index $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ based on the data quality metrics mentioned in Eqs. (18) and (28). Here, ${\epsilon }_{\sigma ,\lim }$ and $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }$ represent arbitrary-defined thresholds for ${\epsilon }_{\sigma }(t)$ and $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$, respectively.

Fig. 10

Process flowchart of the modified metabolic index $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ calculation and screening.

Initially, the retrieved data underwent preprocessing. Raw PPG data from red and IR LEDs were retrieved and accumulated to a certain data length. Then, the raw PPG data were converted to the incremental NIRS signals $\mathrm{\Delta }{N}_{{\mathrm{HbO}}_2}(t)$ and $\mathrm{\Delta }{N}_{\mathrm{Hb}}(t)$. Furthermore, BPF was then applied to each NIRS signal to remove unnecessary high-frequency components and low-frequency components due to respiratory cycles. Here, two different BPF parameter sets were used for different purposes. One parameter set was for $\mathrm{\Delta }{\theta }_{\sigma }(t)$ calculation, which extracts around the HR components of NIRS signals. The other parameter set was for FFT calculation, which extracts 0.8 to 10 Hz components. After applying the BPFs, the excess end portions in each waveform were trimmed so that the first and last points of the data corresponded to the beginning and end of the pulse wave. Furthermore, FFT was applied to the filtered-and-trimmed NIRS signals. Here, each signal for FFT calculation was resampled to make the data length the smallest power of two greater than or equal to the original length. Subsequently, the BPF-applied signals and FFT spectra were employed to calculate $\mathrm{\Delta }{\theta }_{\sigma }(t)$, $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$, ${\mathrm{SNR}}_{\mathrm{Hb}}(t)$, $\alpha (t)$, and ${\mathrm{SaO}}_2(t)$, which comprise ${\epsilon }_{\sigma }(t)$, $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$, and $\mathrm{M}{\mathrm{I}}^{\prime }(t)$. Finally, only those $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ elements that satisfy the specified screening conditions were retained as valid data.

Once the preprocessing was complete, the recorded dataset of $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ was subjected to further postprocessing. Figure 11 shows a flowchart and its visual explanation of the postprocessing.

Fig. 11

Flowchart of the postprocessing for the stored $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ data.

First, the $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ values were divided into 1-min chunks. Second, the representative value within each chunk was calculated. Third, linear interpolation was employed to address the absence of representative values resulting from the low-quality data rejection in the preprocessing stage. Finally, the moving average was applied to the representative values. In this study, the median value was employed as the representative value for each 1-min interval, and the window length for the moving average was set at 30 points, which is equivalent to 30 min for a 1-min interval dataset. This two-stage smoothing approach is beneficial for identifying and eliminating residual outlier values, which have been shown to have a favorable correlation with the reference CGM values from previous research.⁹

3.2.

First-Stage Screening Results for the Obtained Data

Figure 12(a) depicts an example scatter plot of $\mathrm{\Delta }{\theta }_{\sigma }(t)$ versus $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ generated from a typical oral challenge test. Here, for illustrative purposes, the area of the $\pm 30\mathrm{mrad}$ error region is depicted in Fig. 12(a), and the color scale in the plot indicates the spatial density of nearby points. Although the majority of data points were distributed around the diagonal line, notable amounts of outliers were confirmed. Figure 12(b) shows a time-series plot of $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ generated from the same data as Fig. 12(a). Here, the color scale in each data point indicates ${\epsilon }_{\sigma }(t)$, and the dashed red curve shows the smoothed values. Figure 12(b) shows that data points close to the smoothed curve have lower ${\epsilon }_{\sigma }(t)$ and vice versa. This result suggests that outlier data can be excluded by examining ${\epsilon }_{\sigma }(t)$. Then, an appropriate value for the acceptable ${\epsilon }_{\sigma }(t)$ limit ${\epsilon }_{\sigma ,\lim }$ needs to be determined. Figure 12(c) shows the typical transition of the pass rate and the $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$ fluctuation level according to different ${\epsilon }_{\sigma ,\lim }$ values. Here, the typical $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$ fluctuation level was derived by calculating the standard deviations of the difference between adjacent time-domain $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$ data after applying the ${\epsilon }_{\sigma }$-screening. Although it is possible to establish strict data screening by applying a tighter ${\epsilon }_{\sigma ,\lim }$, eliminating the excessive amount of data may result in the generation of inaccurate representative values in the post-processing. To achieve a balance between the two opposing factors of stricter screening and higher pass rate, 10 mrad of ${\epsilon }_{\sigma ,\lim }$, which yields $\sim 50\%$ of the pass rate and a quasi-minimum value in the $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$ fluctuation level, has been adopted in this study. Finally, Fig. 12(d) shows the binarized result of Fig. 12(b) for the passed and rejected data when ${\epsilon }_{\sigma ,\lim }=10\mathrm{mrad}$. With reference to Fig. 12(d), it was revealed that a sufficient amount of outlier values could be eliminated by ${\epsilon }_{\sigma }$-screening while maintaining an adequate quantity of data for the post-processing. Here, it is more probable that the remaining outliers in Fig. 12(d) are the result of a random correlation caused by fluctuations in both $\mathrm{\Delta }{\theta }_{\sigma }(t)$ and $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$. Therefore, it was not possible to eliminate them by applying a stricter threshold. For reference, Figs. S2(a)–S2(d) in the Supplementary Material show binarized results of Fig. 12(b) for various ${\epsilon }_{\sigma ,\lim }$ values.

Fig. 12

(a) Example scatter plot of $\mathrm{\Delta }{\theta }_{\sigma }(t)$ versus $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$. (b) Time-series plot of $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}(t)$ generated from the data of panel (a). (c) Typical transition of the pass rate and the $\mathrm{\Delta }{\theta }_{\mathrm{FFT}}$ fluctuation level according to ${\epsilon }_{\sigma ,\lim }$. (d) Binarized result of panel (b) for the passed and rejected data.

3.3.

Post-Processed Results for the Obtained Data

Figures 13(a)–13(d) show typical post-processed results of oral challenge tests. For comparison, the $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ evaluation result without $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening was indicated in each plot. The CGM delays were adjusted by using the sensor-specific constant value of each sensor, which was determined through comparison with SMBG readings during the pre- or post-experimental period, and the error bars calculated using the typical values of $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$ are indicated in each $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ curve. For the reference CGM values, the 10% vertical error bar range was applied, taking into account the typical accuracy of the sensor.⁴ In addition, the ranges of the left and right vertical axes are adjusted based on the relationship between $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ and BGL as established in the previous research.⁹ In this study, $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }=10\times {10}^{-3}$ of threshold and $n=0.5$ of power exponent in $\alpha$-correction have been applied, and the optimization process will be explained later on. Looking at each of Figs. 13(a)–13(d), it can be seen that the error bar width and the overestimation of $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ have been effectively mitigated through the implementation of $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening.

Fig. 13

(a)–(d) Typical post-processed results of oral challenge tests.

To illustrate the impact of the $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening on $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$ and $\mathrm{M}{\mathrm{I}}^{\prime }(t)$, Figs. 14(a) and 14(b) present histograms of $\delta \mathrm{M}{\mathrm{I}}^{\prime }$ and $\mathrm{M}{\mathrm{I}}^{\prime }$ generated from $\sim 2800$ data points across 30 repeatability tests, with and without the $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening, respectively. Note that the y-axes of Figs. 14(a) and 14(b) are presented on a logarithmic scale to enhance the visibility of outlier values and that $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }=10\times {10}^{-3}$ of threshold and $n=0.5$ of power exponent in $\alpha$-correction have been applied, as is the case in Fig. 13. Moreover, the standard deviation of each histogram, indicated by STD, is presented in the corresponding plot. Looking at each of Figs. 14(a) and 14(b), it was demonstrated that $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening can effectively remove outlier values in $\delta \mathrm{M}{\mathrm{I}}^{\prime }$, thereby reducing the standard deviation of $\delta \mathrm{M}{\mathrm{I}}^{\prime }$, and resulting in eliminating inappropriately large $\mathrm{M}{\mathrm{I}}^{\prime }$ values.

Fig. 14

(a) Comparison of the $\delta \mathrm{M}{\mathrm{I}}^{\prime }$ histograms between the case with and without $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening. (b) Comparison of the $\mathrm{M}{\mathrm{I}}^{\prime }$ histograms between the case with and without $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening.

3.4.

Repeatability Test Results

Figure 15(a) shows a scatter plot of the acquired $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ values versus reference CGM values, generated from the entire repeatability test results, without applying both ${\epsilon }_{\sigma }$ and $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening processes. Here, $n=0.5$ of power exponent was applied in the $\alpha$-correction. For illustrative purposes, the linear approximation and the correlation coefficient $r$ calculated by the linear least squares (LLS) fitting are plotted in Fig. 15(a). Furthermore, the highlighted region around the linear approximation indicates a range of $\pm 5\times {10}^{-3}$. Figure 15(b) shows a Parkes error grid for type 1 diabetes²⁷ that was generated from the results and the conversion coefficients obtained by the LLS fitting presented in Fig. 15(a). Here, each numerical data for $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ corresponding to each CGM data point has been calculated by interpolation, resulting in the generation of 183 data point pairs of $\mathrm{M}{\mathrm{I}}^{\prime }(t)$ and CGM values. The color scale in the plot indicates the spatial density of nearby points. In this study, CGM values were utilized as the reference values for error grid analysis, in lieu of blood glucose values obtained via venous blood glucose testing. Because of this substitution, note that error grids presented in this paper have combined errors from $\mathrm{M}{\mathrm{I}}^{\prime }$ and the CGM. In addition, the mean absolute relative difference (MARD) and the root-mean-square error (RMSE) are also plotted in Fig. 15(b) as BGL-estimation accuracy metrics. Here, in fact, both Figs. 15(a) and 15(b) contain values significantly outside the vertical plot ranges. For reference, Fig. S3 in the Supplementary Material shows zoom-out views of Figs. 15(a) and 15(b).

Fig. 15

(a) Scatter plot of the obtained $\mathrm{M}{\mathrm{I}}^{\prime }$ versus reference CGM values without $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening. (b) Parkes error grid for type 1 diabetes generated from the results presented in panel (a). (c) Scatter plot of the obtained $\mathrm{M}{\mathrm{I}}^{\prime }$ versus reference CGM values with $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening. (d) Parkes error grid for type 1 diabetes generated from the results presented in panel (c).

By applying both ${\epsilon }_{\sigma }$- and $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening to the repeatability test data presented in Fig. 15(a), the scatter plot and the Parkes error grid can be obtained as Figs. 15(c) and 15(d), respectively. Here, the same parameters for $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening and $\alpha$-correction applied in Fig. 13, along with the plotting methods employed in Figs. 15(a) and 15(b), were applied to Figs. 15(c) and 15(d). A comparison of Figs. 15(a) and 15(c) reveals that the number of data points within the $\pm 5\times {10}^{-2}$ region and the correlation coefficient $r$ are greater when the two-stage screening is applied. In the same way, a comparison of Figs. 15(b) and 15(d) reveals that the number of data points within Zone A and the accuracy metrics are improved by applying the screening process. For reference, zoomed-out views of Figs. 15(c) and 15(d) are shown in Fig. S4 in the Supplementary Material. In addition, Fig. S5 in the Supplementary Material presents a transition movie between without and with the two-stage screening process (Video 1, Mp4, 191 KB [URL: https://doi.org/10.1117/1.JBO.29.10.107001.s1]). Moreover, Figs. S6–S9 in the Supplementary Material illustrate the individual results of the entire oral challenge test.

Finally, Table 1 shows the comparison of key performance metrics with and without the two-stage screening process. For purposes of comparison, the performance metrics observed in the previous research using a smartphone camera-based prototype⁹ are also listed in the table. As demonstrated in Table 1, the smartwatch prototype exhibited comparable or enhanced BGL estimation performance by adopting the proposed two-stage screening process. Moreover, although the BGL range and the number of test subjects are limited and cannot be compared directly, the performance metrics obtained through the proposed screening process are comparable to those obtained through other NIGM methods, including Raman spectroscopy²⁸ and radio frequency (RF) spectroscopy.²⁹

Table 1

Comparison of key performance indicators of BGL estimation with and without screening process.

3.5.

Parameter Optimization

Figure 16(a) shows a transition of the correlation coefficient $r$ computed by the LLS fitting according to different $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening threshold $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }$, with the power exponent $n$ for the $\alpha$-correction fixed at 0.5. In this plot, the region below a typical least value of $\delta \mathrm{M}{\mathrm{I}}^{\prime }(t)$ is filled in gray. Figure 16(a) demonstrates that the correlation coefficient $r$ reaches its maximum value at approximately $x=10\times {10}^{-3}$. Figure 16(b) shows a transition of MARD and RMSE according to different $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }$. Similar to Fig. 16(a), both MARD and RMSE exhibited their optimal values around $x=10\times {10}^{-3}$. In general, a more strict threshold leads to more effective control of data variations. However, in this case, an excessively strict threshold leads to a thorough rejection of the obtained data. This results in insufficient data points in the post-processing, which in turn leads to a deterioration in accuracy. Therefore, it can be inferred that the modest screening threshold exhibited the best performance.

Fig. 16

(a) Transition of the correlation coefficient $r$ computed by the LLS fitting according to different $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening threshold $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }$. (b) Transition of MARD and RMSE according to different $\delta {\mathrm{M}{\mathrm{I}}^{\prime }}_{\lim }$.

Next, Fig. 17(a) shows a transition of the correlation coefficient $r$ computed by the LLS fitting according to different power exponent $n$ for $\alpha$-correction, with the threshold for the $\delta \mathrm{M}{\mathrm{I}}^{\prime }$-screening fixed at $10\times {10}^{-3}$. Here, the area where $\alpha$-correction is non-applicable due to significant deterioration in accuracy is filled in gray. In this case, the correlation coefficient $r$ attained its optimal value at $n=0.5$, which implies that the AC amplitude of the optical path length, $L_{\mathrm{AC}}(t)$, is approximately proportional to the square root of the total optical path length, $L(t)$. Similarly, Fig. 17(b) shows a transition of MARD and RMSE according to different power exponents $n$. The optimal values for MARD and RMSE were observed at $n=0.5$ in this case as well.

Fig. 17

(a) Transition of the correlation coefficient $r$ computed by the LLS fitting according to different power exponent $n$ for the $\alpha$-correction. (b) Transition of MARD and RMSE according to different power exponents $n$.