## 2 Methods
### Real-Time Signal Processing  
Figure 2 shows a flow diagram of our present method (rt-ShortPCA GLM based on the method of off-line ShortPCA GLM). The key aspect of ShortPCA GLM is the design matrix, which contains a global temporal pattern of SBF explained by the first principal component of PCA. However, due to real-time restrictions, we cannot define this element before measuring samples. For this reason, we applied SWA in each signal step individually.  
When a new raw sample was measured, at first, only \(L\)-sample-length-windowed long- and short-channel samples had a 2 nd-order Butterworth bandpass filtering with cut-off frequency of 0.01-0.3 Hz applied (these values are similar to those used in previous studies; for review, see [1]), then the newest sample was stored as filtered data (Figure 2I-II). Before storing, z-scoring (mean \(\mu\) and standard deviation \(\sigma\) during all filtered samples were 0 and 1, respectively) was applied only to the newest samples of short-channels (Figure 2 II), in order to prevent extracted principal components for the next PCA step from being biased to a specific channel signal because of motion artifacts, signal-to-noise ratios, or different path lengths [14; 18]. Furthermore, to estimate a global temporal pattern of SBF (global scalp-hemodynamics model; GSHM), we applied PCA to the windowed and preprocessed samples and extracted principal components. The newest sample of the first principal component was then stored to the design matrix of GLM (Figure 2 III-IV).  
For the final step, we employed GLM with SWA to estimate the newest sample of CBF by removing the weighted GSHM. The GLM expresses response variable \(\mathbf{y}_{\text{Long}}\) (time sample \(N_{T}\times\) number of channels \(N_{\text{CH}}\)) in terms of a linear combination of explanatory variables (design matrix) \(\mathbf{X}\) (\(N_{T}\times\) number of elements \(N\)):  
\[\mathbf{y}_{\text{Long}}=\mathbf{X}\mathbf{\beta}+\mathbf{\varepsilon}, \tag{1}\]  
where \(\mathbf{\beta}\) is an unknown weight matrix (\(N\times N_{\text{CH}}\)) as relative contributions of the assumed elements in long-distance channels. \(\mathbf{\varepsilon}\) is a normally distributed model error term with a mean of zero and variance of \(\sigma^{2}\) (i.e., \(\mathbf{\varepsilon}\sim N(0,\sigma^{2})\)). In this study, we constructed a design matrix with five elements similar to the off-line ShortPCA GLM: cerebral-hemodynamic model (CHM) with its temporal and dispersion derivatives (CHM\({}^{(1)}\),CHM\({}^{(2)}\)), constants (Const.) and a global temporal pattern of SBF (GSHM) estimated in real-time PCA. Here, the CHM was generated by convolving the mathematical hemodynamic response function [19] and experimental task function (box-car function). All elements without the GSHM were normalized by scaling the maximum peak amplitudes to one. Note that all elements without the GSHM were able to prepare before real-time processing, thus we only had to update the GSHM for each measuring time. The optimized estimation of weight vector at time \(k\), \(\mathbf{\hat{\beta}}_{k}\), was calculated by ordinary least squares (Figure 2 IV) as follows:  
\[\mathbf{\hat{\beta}}_{k} = \mathbf{X}_{k}^{\dagger}\mathbf{y}_{\text{Long},k}, \tag{2}\] \[\mathbf{X}_{k}^{\dagger} = (\mathbf{X}_{k}^{T}\mathbf{X}_{k})^{-1}\mathbf{X}_{k}^{T}. \tag{3}\]  
Here, \(N_{T}\) was restricted to \(L\)-samples by SWA so that each element of the matrices in Equations (2) and (3) had the samples ranging from \(k\)-\(L\)+1 to \(k\) (see Figure 1 again). Owing to the changing data samples in the sliding-window and the transition of \(k\), we could estimate weight \(\mathbf{\beta}\) dynamics corresponding to the amplitudes at time \(k\) of the assumed elements in the design matrix. After the weight estimation, we can remove global scalp-hemodynamics at time \(k\) from each long channel by subtracting the real-time estimated GSHM multiplied by the \(\mathbf{\hat{\beta}}_{\text{GSHM},k}\) for the corresponding channel.