Enhancing PDC Functional Connectivity Analysis for Subjects with Dyslexia Using Artifact Cancellation Techniques

BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER Enhancing PDC Functional Connectivity Analysis for Subjects with Dyslexia Using Artifact Cancellation Techniques Taha Mahmoud AL-NAIMI 1, Shanthini Chandra Sekaran NAIDU 2, Ahmad Zuri SHA’AMERI 1, Norlaili Mat SAFRI 1, Narina Abu SAMAH 2 1Department of Electronic and Computer Engineering, School of Electrical Engineering, Faculty of Engineering, Universiti Teknologi Malaysia, Skudai, 81310 Johor Bahru, Johor, Malaysia 2School of Education, Faculty of Social Sciences and Humanities, Universiti Teknologi Malaysia, Skudai, 81310 Johor Bahru, Johor, Malaysia [email protected]y, cshan[email protected]y, [email protected]y, [email protected]y, [email protected]y DOI: 10.15598/aeee.v20i4.4525 Article history: Received Apr 11, 2022; Revised Aug 02, 2022; Accepted Aug 30, 2022; Published Dec 31, 2022. This is an open access article under the BY-CC license. Abstract. The neurobiological origin of dyslexia allows the study of this disorder by examining functional connectivity between regions of the brain. During rest-state or at task completion, Electroencephalograms (EEG) are used to observe brain signals. By using Partial Directed Coherence (PDC) analysis, the correct analysis of functional connectivity was assessed. In spite of that, the estimation of functional connectivity can be inaccurate due to the presence of artifacts. Several methods have been employed by researchers to remove artifacts, including Moving Average Filters (MAF), Wiener Filters (WF), Wavelet Transforms (WT), and hybrid filters. Despite this, no research has been conducted on the effects of artifact removal methods on functional connectivity. Consequently, Artifact Cancellation (AC) algorithms are developed to reduce the effects of eye blinks, eye movements, and muscle movements on functional connectivity estimation. In this work, the denoising filters discussed earlier are utilized as part of the AC algorithm. Additionally, a comparison was conducted to determine the effectiveness of the filters. According to the results, AC-MAF removed all artifacts with the least computational complexity after improving the MAF. In order to test its efficacy in real-world conditions, it was applied to the real signals recorded while children with dyslexia were participating in rapid automatized naming activities. Utilizing the PDC approach, the developed algorithm accurately assessed functional connectivity. Keywords Artifacts Cancellation (AC), computational complexity, dyslexia, functional connectivity, Partial Directed Coherence (PDC), Rapid Automatized Naming (RAN). 1. Introduction Dyslexia is commonly conceptualised as a specific learning deficit that affects the ability to identify lettersound correspondence (phonemic awareness), which in turn affects decoding ability. It is considered a neurobiological disorder due to the differences in the network of neurons and how one part of the brain communicates with other regions [1] and [2]. Hence, brain imagery techniques, such as Electroencephalography (EEG), Magnetoencephalography (MEG), Positron Emission Tomography (PET), and Functional Magnetic Resonance Imaging (fMRI) are required to examine the complex cortical brain connectivity [3]. EEG was selected due to its excellent temporal resolution that displays the connectivity within the brain regions during any given time intervals [4]. Further, EEG can be useful in identifying brain functional connectivity in children with dyslexia when they perform Rapid Automatized Naming (RAN) tasks. RAN has been used widely as it is one of the predictors for early reading development [5]. It is based on how rapid and accurate ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 592 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER a person names a string of digits, objects, colours, and letters [6]. The brain functional connectivity can be analyzed using methods such as Granger Causality (GC), Direct Transfer Function (DTF), and Partial Directed Coherence (PDC) were developed to detect the functional connectivity between two or more channels [7]. Moreover, the GC method depends on the predicted error for the model prediction, whereas the DTF and PDC methods vary depending on the predicted coefficients [8]. However, since EEG signals can be contaminated with a non-neural signal known as artifacts, they must be removed during the pre-processing step before the clean signal can be analysed to see the functional connectivity. The artifacts always corrupt EEG signals due to eye movement and blinking, as well as muscular movements such as the jaw, neck movements and swallowing, causing additional interaction between the channels [9] and [10]. Also, it can be caused by the interference of power lines (50/60 Hz), substandard technology, and electronic devices that create electromagnetic field interferences [9]. Despite that, there are various methods used to remove artifacts from EEG signals and one of the commonly used methods based on Blind-Source Separation (BSS) is known as the Independent Component Analysis (ICA) method [11], [12] and [13]. However, this method still requires a reference signal to select the right source of artifacts, or it can be done by asking the subjects to produce artifacts [14]. Therefore, one-channel methods using denoising filters such as: Discrete Wavelet-based Denoising (DWD), Stationary Wavelet-based Denoising (SWD), Wiener Filter (WF), Median Filter (MF), and Moving Average Filter (MAF), were developed by researchers in [15], [16], [17], [18] and [19] to remove the artifacts. Despite the successful results achieved by these techniques in removing either ocular or muscular artifacts, each technique has its own limitations. More recently, a combination of methods (known as a hybrid method) was used to enhance the performance of the artifacts removal technique, and to evaluate a single EEG channel [20] and [21]. One of these methods is the combination of DWD and SWD. Nonetheless, the computational complexity of using a hybrid technique can be complex and take longer process time compared to the standard denoising techniques [15]. Also, there are no studies conducted to investigate the effectiveness of removing EEG artifacts on EEG functional connectivity. As a consequence, it is crucial to improve artifact removal methods and ensure proper functional connectivity estimation. In this paper, the Artifact Cancellation (AC) algorithm is developed and implemented with various filters, those are: MF, MAF, DWD, SWD, WF, and hybrid of DWD/SWD with Linear Prediction Filter (LPC), as well, with WF. Additionally, this study compares the computational complexity of the modified MAF to the other filters in the proposed AC algorithm. Furthermore, simulated and real EEG signals were analyzed from a PDC perspective before and after the artifact removal process. This is to ensure estimating the correct functional connectivity between the brain areas produced by the brain (clean EEG) signal. The remaining of this paper is organized as follows: Sec. 2. covers the brain anatomy, including the areas involved in the learning of the reading process. Section 3. describes the signal mode and includes EEG signal characteristics, characteristics of artifacts, and EEG signal measurement procedure. Section 4. presents artifacts removal methods, including the AC algorithm. Section 5. presents the results for simulated and real signals, including the effects of the artifacts on the estimation of functional connectivity using PDC, the consequences of using the proposed method, and the computational complexity for each method. Conclusions are stated in Sec. 6. 2. Brain and Dyslexia The reading process is a recent cultural invention and the human brain adapted to this activity by making connections between several regions [22]. Although the entire brain works during any given task, multiple studies have revealed that there are essential areas in the left hemisphere that play significant roles in acquiring reading skills as well as reaching reading mastery [23], [24], [25], [26], [27] and [28]. Dyslexia occurs across different levels of intelligence, but the prominent identification factor is the tendency for an individual to struggle in mastering the reading process or to break the ‘reading code’ despite mainstream educational experiences [29]. At the early stage of learning to read, children master how letters (grapheme) relate to their corresponding sounds (phoneme) [30]. As children with dyslexia face challenges in grasping these basic reading skills, they are unable to gain sufficient knowledge since the education system relies heavily on text mastery. Therefore, early detection is vital to provide a systematic intervention to help children with dyslexia to learn to read [31]. A pathway is also formed in the left hemisphere that is dominant in the reading process. In this pathway, as one views the visual stimuli, the occipital lobe becomes activated and sends signals to the temporal lobe and Wernicke’s area. The temporal area is important to access known words from memory, whereas the Wernicke’s area is vital to make connections between letter and sound based on the visual stimuli (word or string of words). Finally, the signal is sent to the Broca’s ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 593 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER area (phonological processor) before proceeding to the motor cortex, which controls the muscles in the lips, tongue, and face [32]. The brain regions crucial during the reading process are shown in Fig. 1. Thus, any disconnection in the pathway may lead to dyslexia or reading difficulties, depriving their ability to decode words [24]. Temporo-parietal cortex Visual cortex Occipito-temporal cortex Inerior frontal cortex Fig. 1: Reading circuit [6]. 3. Signal Model and Analysis In this section, EEG signals were modelled and described as a Multivariate Autoregressive (MVAR) process, followed by the description of artifacts, their characteristics and causes. Also, in the last sub-section, the effect of artifacts on EEG signals is presented. 3.1. Multivariate Autoregressive Model The EEG signals that are measured and collected by attaching electrodes to the subjects’ heads are nonstationary signals – time-varying signals that will arise during any activation in the brain regions. It can be simulated by adding signals epochs with different frequencies chosen randomly in the range between 0.04–120 Hz that includes the bands: Delta, Theta, Alpha, Beta, and Gamma [33]. Nevertheless, there are no correlations between those signals in terms of time and frequency. Therefore, the MVAR model was chosen for this matter due to its similarity to the EEG signals [34] and [35]. It is identical to an Autoregressive (AR) signal, but usefull for more than one-time series model (Mchannels) and all the channels interact with each other by a set of coefficients. These coefficients will reflect the relationship between the channels and the strength of functional connectivity between two or more channels [36]. The MVAR model can be written in the following form:      x1(n) x2(n) . . . xM(n)      = p X r=1 Ar      x1(n−1) x2(n−1) . . . xM(n−1)      +     w1(n) w2(n) . . . wM(n)      ,(1) where wi(n)is an uncorrelated white noise source, xi(n)is the stationary processes for each channel, pis the order. The coefficients Arin Eq. (1) is defined as: Ar=          a11 (r)a12 (r). . . a1M(r) . . .. . .. . .. . . . . .. . .aij (r). . . . . .. . .. . .. . . aM1(r). . . . . . aMM (r)          ,(2) where the coefficients aij (r)are the interaction effect due to xj(n−r)on xi(n)for each delay point r. 3.2. Artifacts Model The MVAR models simulate the functional connectivity between two or more channels, where the parameters (the optimum order of MVAR model and the coefficients) are not constant, and they will change for each second because those channels are time-variant [37]. Apart from that, artifacts will result in new interaction between channels leading to the wrong estimation of MVAR parameters and hence the wrong functional connectivity [38]. Consequently, the artifacts are modelled according to their type and added to the MVAR model’s channels. The MVAR model for five channels as explained in [39] is used in this paper, with artifacts added as described in the following equation: x1(n) = 0.95√2x1(n−1) + 0.9025 x1(n−2)+ +s(n) + w1(n), x2(n) = 0.5x1(n−2) + w2(n), x3(n) = −0.4x1(n−3) + w3(n), x4(n) = −0.5x1(n−2) + 0.25√2x4(n−1)+ + 0.25√2x5(n−1) + w4(n), x5(n) = −0.25√2x4(n−1)+ + 0.25√2x5(n−1) + s(n) + w5(n), (3) where xi(n)represents the simulated EEG signals for 5 channels with sampling frequency of 250 Hz that follows the Nyquist–Shannon sampling theorem for sampling the continuous-time signals with a finite bandwidth. The artifacts (eye blinking/movement and muscular artifacts) are modelled by s(n)that is defined as: ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 594 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER s(n) = b(n) + o(n) + m(n),(4) where b(n),o(n), and m(n)represent the blinking, eye-movement, and muscular artifacts, respectively. The eye blinking artifacts occur in a period of 15–20 times per minute [40], when the duration of each blink is 150–400 msec with high magnitude (800 µV and above) [41]. On the other hand, eye-ball movement has a tremendous amplitude compared to blinking at a frequency of below 5 Hz [16]. The characteristics of blinking signals are described in [42]. Thus, this signal is a convolution of a unit sample response h(n)with pulse train x(n)that can be expressed as: b(n) = M−1 X m=0 h(m)x(n−m),(5) where h(n)is a Hamming function with size (M= 80) samples, while the pulse train x(n)is designated as: x(n) =   10 0 . . . 0 | {z } ksegments 1. . .  ,(6) where the space (segments) is randomly selected for each period (k∈[38 : 50]) to follow the characteristics of eye blinking artifacts. Consequently, each peak’s duration is 100 samples, while the period between two blinking peaks is [750 : 1000] samples. As a result, simulated eye blinking is shown in Fig. 2. Fig. 2: Simulated eye-blinking signal [42]. The eye-movement signal is simulated using half a period of a sine wave with immense magnitude, along with a frequency lower than 5 Hz [16]. The signal can be defined as: o(n) =    Asin(2πf1n) 0 <x<(N−1) 2, 0otherwise, (7) where Ais amplitude, Nis the number of samples within an observed interval, and f1is frequency. For simulation purposes, the amplitude and frequency were selected as 20 and 4 Hz, respectively [16]. The frequency range of muscular movements ranges from 20 to 50 Hz, and it is difficult to remove them from the EEG signals because they overlap one another [43], [43] and [45]. Muscular signals were simulated according to the attributes of Electromyography (EMG) signals [46] and [47] where the muscle fibre membranes create an action potential during depolarization and repolarization of the muscle. Figure 3 shows an example of a Motor Unit Action Potential (MUAP) and its firing frequency as one important factor influencing both magnitude and density. The mathematical model for the EMG signal can be expressed as: m(n) = N−1 X r=0 MU(r)w(n−r),(8) where m(n)is the simulated EMG signal, MU(n)is the MUAP, w(n)is point processed (firing impulse), and Nis the number of motor unit firing. The MUAP can further be defined as: MU(n) = M−1 X l=0 A(l) sin(2πfln),(9) where Mis the number of muscle-fiber action potentials, and A(l)represents the l-th amplitude for each action potential that would be randomly selected – positive or negative for each (l). The firing frequency is defined in (fl), which is randomly chosen in the interval (fl∈[20 : 50]). Superposed signal of the whole motor unit Muscle Fiber Detection Site Motor Endplate Action Potentials α motonuron 1 2 3 n + + + = (a) Generation of MUAP. Voltage (mV) Motor Unit Firing MU1 (3 Hz) MU1 (3 Hz) MU2 (4 Hz) MU3 (6 Hz) MU4 (8 Hz) Superposed Surface Signal = + + + Motor Unit Recruitment (b) Requirement and firing frequency. Fig. 3: Requirements and firing frequencies of MUAP [44]. Figure 4 shows the simulated signal with effects from ocular and muscular artifacts on time and frequency domains. The simulated blinking effect appears in the same time domain as sharp peaks, while the simulated eye movements appear as the typical shape of a half sine wave. On the other hand, the simulated muscular artifacts overlap with the simulated EEG signal but have slightly higher amplitude compared to the simulated EEG signal, as shown in the time intervals ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 595 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER (10–20), (40–60), and (90–110). In the frequency domain, major effects of ocular artifacts are seen within the (0–10) Hz range, while for muscular artifacts, the effect overlaps with the simulated EEG within (20–50) Hz range. (a) Time representation. (b) Power spectrum. Fig. 4: Simulated signal for channel x5(n)with artifacts. 3.3. Functional Connectivity In brain neural networks, nodes (channels) are connected to each other and transmit signals from one node to another for example, 1 →2, 2 →3, and 3 →4 with arrows showing the direction of information flow. This flow can be detected by using statistical relationship methods such as the cross-correlation method, by which leading and lagging signals are detected from output delays. However, it cannot differentiate direct and indirect interactions between the channels, such as the indirect interaction between 1 →4. Moreover, correlation or coherence methods can be used only between two channels [48], while in brain neural networks, each channel can be affected by more than one channel. Furthermore, as stated in Granger Causality (GC) [49] and [50] , if a process in channel Yis caused by another process in channel X, then the prediction of the earlier one is enhanced by using information from the latter one. Therefore, methods such as GC, Direct Transfer Function (DTF), and Partial Directed Coherence (PDC) were developed to detect the functional connectivity between two or more channels [7]. Moreover, the GC method depends on the predicted error for the model prediction, whereas the DTF and PDC methods vary depending on the predicted coefficients [8]. Both DTF and PDC can estimate the functional connectivity since both are applicable to EEG signals. However, the latter was selected due to its ability to detect only direct interactions for more than two channels [34]. Additionally, applying the DTF method to more than two channels cannot detect the correct direction of information [39]. The PDC is computed for the estimated MVAR coefficients in Eq. (2), as defined in [39] using the following equation: πij(f) = |aij(f)| qaH j(f)aj(f) ,(10) where aij(f)represents the difference between the Fourier transform of the coefficient series Arin Eq. (2), and an identity matrix of the same size as Ar. While πij is the PDC result - between 0 and 1 where the higher value represents stronger functional connectivity between i←jat a frequency f. The average functional connectivity between each of the two channels among the frequencies was calculated using Eq. (11) from [51], as described: πij(f) = Z∞ 0 f πij(f)df, (11) where πij(f)is the average of the functional connectivity between i←j. The result from this equation was then compared to a threshold γto determine if functional connectivity was present. The threshold range of (0.05–0.2) was described in [52]; thus, the selected threshold in this work is (γ= 0.2). The number of coefficients for each channel depends on the order of the AR signal [53] and the parameters (the optimum order of MVAR model and the coefficients), which can be estimated using the ARFIT algorithm [36]. The algorithm uses different methods, such as Akaike Information Criterion (AIC) and Schwarz’s Bayesian Criterion (SBC), for the MVAR model estimation. Nevertheless, the latter is preferable for time series analysis [54]. Since the characteristics of the EEG signals changes in time, the signals are segmented into epochs to satisfy local stationary conditions and compute the PDC accordingly. Meanwhile, the number of samples Nwas discussed in [37] and [49] and described as: N > Mp, (12) ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 596 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER where Mis the number of channels and pis the order as defined in Eq. (1). Therefore, the selected size for each epoch is 250 samples (1 sec for a sampling frequency of 250 Hz), and the number of parameters for each delay (r) should be as low as possible (p≤5) [55]. The PDC analysis results for Eq. (3) can be plotted as shown in Fig. 5, where the vertical and horizontal axis represents the magnitude and frequency in samples, respectively. Fig. 5: PDC analysis of the simulated signal from Eq. (3) [39]. From Fig. 5, the functional connectivity was considered when the average was equal to or higher than the selected threshold (γ= 0.2). Consequently, it can be observed that the diagonal plots are higher in magnitude as compared to the other plots due to the representation of how much a channel is connected to itself. The off-diagonal plots represents the functional connectivity of the four channels, namely: 1 to 2, 1 to 3, 1 to 4, and 4 and 5 from/to each (the desired results of the PDC). Moreover, the functional connectivity caused by artifacts’ effect such as 4 to 1, 5 to 1, and 1 to 5, were demonstrated. Therefore, artifacts need to be removed to ensure the correct estimation of the parameters in Eq. (1). 4. Artifacts Removal Method Artifacts Cancellation (AC) algorithm was developed based on a previous approach [42] to remove the major artifacts such as ocular and muscular artifacts. Also, most of the artifacts cancellation methods, as described in [16], [17], [33] and [40] follows the structure shown in Fig. 6. The denoising filter used in the AC structure can be any filtering or denoising method. In Fig. 6, x(n)is the input signal to the denoising filter – by which, it is the raw EEG recording received from one of the channels or the simulated signals as defined in Eq. (3). The output of the denoising filter is y(n), while the clean EEG recording z(n)is the difference between system output y(n)and the input signal x(n). The various AC algorithm described in this section can be implemented using Median Filter (MF), Denoising Filter ∑ ()yn AC x(n)_ +z(n) Fig. 6: AC algorithm structure for single EEG channel. Moving Average Filter (MAF), Discrete Wavelet-based Denoising (DWD), Stationary Wavelet-based Denoising (SWD), Wiener Filter (WF), and hybrid of wavelet denoising methods with linear prediction filter or with Wiener filter. 4.1. Wiener Filter A Wiener filter is an optimum filter that derives coefficients from the desired and noisy observed signals by minimizing the mean square error (accuracy depends on the signal-to-noise ratio). It is used in [17] to remove eye-blinking artifacts by estimating the artifacts signal and subtracting the estimated signal from the raw EEG data, as shown in Fig. 7. In other words, simulated EEG data was considered as a noisy signal, while blinking signal is approximately deterministic since it would be modelled as a mathematical function as described in Eq. (5). y(n) ()hn Reference signal [s(n)] ∑ _ + x(n)= s(n)+v(n)e(n) Fig. 7: Block diagram for Wiener filter. In Fig. 7, s(n)represents the signal – in this case the blinking signal, while v(n)is the white noise signal and x(n)is the combination of both. The clean artifacts represented by y(n)is the output of the optimum linear filter h(n), and e(n)is the error that is used to optimize the impulse response filter h(n)by minimizing the mean-square error. Moreover, for the Wiener filter to function effectively, the reference signal must be the same as the signal of interest. 4.2. Wavelet Denoising Non-stationary signals such as EEG can be analyzed using Wavelet Transform (WT) [15]. For EEG signals, the WT can be applied using either Discrete Wavelet Transform (DWT) or Stationary Wavelet Transform (SWT) [16]. The DWT decomposes the signals to L ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 597 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER level through high-pass filters and low-pass filters using bases functions known as wavelets [16]. Next, downsampling is applied to reduce computational complexity and obtain the detail coefficient (DL) as well as approximate coefficient (AL). Next, a thresholding algorithm is staged at each (DL) to ensure the correct selection of artifacts’ coefficients without selecting the EEG signal coefficients. After that, upsampling was applied before reconstructing the decomposed coefficients into a single signal. This process is known as Discrete Wavelet-based Denoising (DWD). The block diagram of DWD is described in Fig. 8. DWT Decomposition ()xn ()yn Threshold DWT Reconstruction DWD Fig. 8: Block diagram for discrete wavelet denoising. EEG signals have a band frequency range of 0.5 to 120 Hz, while the band frequency for eye movement, blinking, and muscular activity are 0 to 7 Hz, 8 to 13 Hz and 20 to 50 Hz, respectively. Moreover, the DWD decomposition should have five levels (as shown in Fig. 9), and the artifacts caused by eye movement, blinking, and muscular contraction should be at decomposition levels 5, 4, and 3, respectively. After that, the statistical thresholding is used for denoising at each decomposition level with artifacts. In the current case, the selected threshold was applied on the detailed coefficients to allow the artifacts signals’ coefficients to pass each decomposition level and cancel the EEG signals. Moreover, the wavelet used in this work was symlet (sym3) due to its effectiveness, as stated in [12], [15], [16] and [21]. D1D2D3D4D5A5 125 Hz 62.5 Hz 31.2 Hz15.6 Hz 7.8 Hz 3.9 Hz 0 Eye movement Blinking Muscular movement Fig. 9: Wavelet decomposition levels and the presence of artifacts [42]. In contrast to DWD, the SWD decomposes signals without a downsampling algorithm to maintain all the coefficients within the decomposition levels. This step is helpful for essential wave shapes (especially with white noises) [18] and pulses – in this case, blink-signal shape and occurring time for localization calculations. However, the computing time is longer due to the exclusion of the downsampling algorithm, which retains the redundancy in the calculations. Similar to the DWD, the wavelet used is (sym3) and the decomposition level for the SWD is five, as shown in Fig. 9. 4.3. Moving Average Filter The MAF employs a window function to smooth the signal from noise. It is frequently used in the time domain due to the weak frequency response (unable to separate one band of frequency from another) caused by the weak stopband attenuation [56]. Yet, it provides the faster step response for a given noise reduction. After using this filter in a three-stage algorithm called gradient artifact correction in [57], the clean EEG signals were restored, where (N=155) was the estimated filter size. However, the size of the filter in this case does not follow the characteristics of the MAF to remove the artifacts. The filter length is inversely proportional to the frequency; therefore, increasing the length of the filter will maintain low frequencies and vice versa. For example, if the selected filter size is (N=10) and the sampling frequency is 250 Hz, then the cut-off frequency is 25 Hz, as described in [58] using the following equation: N=fs f,(13) where Nis filter size, fsis the sampling frequency, and the artifacts’ frequency is f. Figure 10 shows the frequency representation for various filter lengths following Eq. (13). Ocular artifacts Muscular artifacts Fig. 10: Frequency response of moving average filter for various filter lengths. The length for the MAF can be determined using the presence of the artifacts in the time and frequency domains (as described earlier in Sec. 3.2. ), where the magnitude of the artifacts is higher than the magnitude of the clean EEG signals in the time domain. Also, since the muscular artifacts are within the range of 20–50 Hz and the sampling frequency is 250 Hz, the selected size of the filter is (N= 5), which covers up both ocular and muscular artifacts. ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 598 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER 4.4. Median Filter Median Filter (MF) is a non-linear filter that relies on a sorting algorithm to select the middle-ranked value within a set of data. Since MF is inherently an imaging technique [59], it is usually used for two-dimensional signals to maintain white noise-like signals and to remove impulse artifacts. Nevertheless, the filter can be used in the process of one-dimension signals [60]. For example, if the data set is x={6,5,4,9,3}, the data will be sorted following the merge sort algorithm (shown in Fig. 11) to become {3,4,5,6,9}, and the middle value is (y= 5), as described in [60] according to the equation: y(n) = median (X(n)) , X(n) = x(n−N 2), . . . , x(n), . . . x(n+N 2), (14) where y(n)is the output of the median filter, x(n)is the input data, and Nis the window size of the filter. { 6 , 5 , 4 , 9 , 3 } { 6 , 5 , 4 } { 9 , 3 } { 3 , 9 } { 9 } { 3 } { 6 , 5 } { 5 , 6 } { 6 } { 5 } { 4 } { 4 , 5 , 6 } { 3 , 4 , 5 , 6 , 9 } Unsorted Array Sorted Array Fig. 11: Visualization of merge sort algorithm [64]. Although the MF filter is similar to the MAF, it is deemed to be better than the latter one due to the fact that it does not average the values specified in the window [62], [63] and [64]. The MF was used previously in [64], [65] and [66] to remove the ocular artifacts from the EEG signals where the window size used to remove the artifacts was three (N= 3). Therefore, this work will be performed with the same window size. 4.5. Hybrid Methods In recent years, researchers have focused on using a hybrid strategy to remove the artifacts by merging two approaches [19], [20] and [21]. These techniques remove artifacts more accurately from the EEG signals compared to a single method. Besides, hybrid methods can estimate the reference signals and eliminate the need for these signals. In this work, the hybrid methods are Linear Predictive Coding (LPC) with DWD and LPC with SWD, as shown in Fig. 12. Also, the hybrid WF with DWD was used in the proposed study, as shown in Fig. 13. ()xn LPC coefficients Whitening Dewhitening ()yn DWD/SWD Fig. 12: Block diagram for (LPC-WD) hybrid denoising filter. Based on Fig. 12, the whitening filter (whitening transform) uses the inverse of an all-pole system [58] – in this case, derived from the LPC predictor – to whiten the signals. This step increases the level of the EEG signal to ensure the removal of the EEG signal from the raw signals at the wavelet denoising stage. Subsequently, de-whitening is applied to retrieve the clean artifacts’ signals to be subtracted from the raw EEG signals in the AC algorithm. In Fig. 13, the output of the DWD filter was used as a reference signal for Wiener filter, to ensure all the EEG signals were removed and retrieve clean artifacts in the output. Subsequently, the clean artifacts’ signals will be subtracted from the raw EEG signals in the AC algorithm. ()yn ()en ()hn DWD (reference) ∑ _ + x(n)= s(n)+v(n) ˆ()sn Fig. 13: Block diagram for (WF-DWD) hybrid denoising filter. 5. Results In this section, functional connectivity using PDC was used to compare the effect of artifacts on simulated and real EEG signals. Also, the performance measures used were based on the Probability Density Function (PDF) and Cumulative Density Function (CDF) for the signals before and after removing the artifacts. ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 599 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER 5.1. Performance for Simulated Signal From the characteristics of the EEG signals modelled as MVAR, the distribution of the signal’s amplitude measured at each channel should be symmetric about the mean and approximately follows a Gaussian distribution. (a) PDF before using AC. (b) PDF after using AC DWD. (c) PDF after using AC SWD. (d) PDF after using AC MAF. (e) PDF after using AC LPC-DWD. (f) PDF after using AC LPC-SWD. (g) PDF after using AC WF-DWD. (h) PDF after using AC MF. Fig. 14: PDF comparison for various AC algorithms. Figure 14 shows the PDF applied before and after removing the artifacts using MF, MAF, DWD, SWD, and hybrid methods, such as LPC-DWD/SWD and WF-WT. The artifacts caused a non-symmetric distribution for the amplitude at the lower range from −20 to −8and upper range from 8 to 40, as shown in Fig. 14(a). Nevertheless, the distribution of amplitude is symmetric at around zero mean after the artifacts are removed, as shown in Fig. 14(b), Fig. 14(c), Fig. 14(d), Fig. 14(e), Fig. 14(f), Fig. 14(g) and Fig. 14(h) In order to verify conformity to Gaussian PDF, a comparison was made between the distribution of the measured amplitude and a standard normal Gaussian distribution, as shown in Fig. 15. Similarly, this is also true for the CDF after artifact removal, resulting in a CDF closer to the standard normal CDF. The AC MAF and AC DWD results were also closest to standard normal CDF compared to other methods. (a) CDF before using AC. (b) CDF after using AC DWD. (c) CDF after using AC SWD. (d) CDF after using AC MAF. (e) CDF after using AC LPC-DWD. (f) CDF after using AC LPC-SWD. (g) CDF after using AC WF-DWD. (h) CDF after using AC MF. Fig. 15: CDF comparison for various AC algorithms. Based on Fig. 15 the removal of artifacts does not guarantee the restoration of functional connectivity between the various channels. Therefore, a comparison of functional connectivity between the various channels for the simulated signals before and after removing the artifacts is necessary. Figure 16 shows the functional connectivity at various threshold ranges from 0.05 to ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 600 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER B. GUSCHLBAUER, C. H. LALCKING, R. DAHLHAUS and J. TIMMER. Testing for directed influences among neural signals using partial directed coherence. Journal of Neuroscience Methods. 2006, vol. 152, iss. 1, pp. 210–219. ISSN 0165-0270. DOI: 10.1016/j.jneumeth.2005.09.001. [39] BACCALA, L. A. and K. SAMESHIMA. Partial directed coherence: A new concept in neural structure determination. Biological Cybernetics. 2001, vol. 84, iss. 6, pp. 463–474. ISSN 14320770. DOI: 10.1007/PL00007990. [40] SREEJA, S. R., R. R. SAHAY, D. SAMANTA and P. MITRA. Removal of Eye Blink Artifacts from EEG Signals Using Sparsity. IEEE Journal of Biomedical and Health Informatics. IEEE, 2018, vol. 22, iss. 5, pp. 1362–1372. ISSN 2168-2208. DOI: 10.1109/JBHI.2017.2771783. [41] YANG, B., K. DUAN, C. FAN, C. HU and J. WANG. Automatic ocular artifacts removal in EEG using deep learning. Biomedical Signal Processing and Control. 2018, vol. 43, iss. 1. pp. 148–158. ISSN 1746-8094. DOI: 10.1016/j.bspc.2018.02.021. [42] AL-NAIMI, T. M., A. Z. SHA’AMERI and N. M. SAFRI. Effect of Artifacts on the Interpretation of EEG-Based Functional Connectivity Estimation Using Partial Directed Coherence. In: Advances in Intelligent Systems and Computing. Singapore: Springer, 2022, pp. 253–268. ISBN 978-981-16-2123-9. DOI: 10.1007/978-98116-2123-9_19. [43] DIESSEN, E. van., T. NUMAN, E. van DELLEN, A. W. van der KOOI, M. BOERSMA, D. HOFMAN, R. van LUTTERVELD, B. W. van DIJK, E. C. W. van STRAATEN, A. HILLEBRAND and C. J. STAM. Opportunities and methodological challenges in EEG and MEG resting state functional brain network research. Clinical Neurophysiology. 2015, vol. 126,iss. 8, pp. 1468–1481. ISSN 1872-8952. DOI: 10.1016/j.clinph.2014.11.018. [44] CHEN, X., A. LIU, J. CHIANG, Z. J. WANG, M. J. MCKEOWN and R. K. WARD. Removing Muscle Artifacts from EEG Data: Multichannel or Single-Channel Techniques? IEEE Sensors Journal. IEEE, 2016, vol. 16, iss. 7, pp.1986–1997. ISSN 1558-1748. DOI: 10.1109/JSEN.2015.2506982. [45] FROLICH, L. and I. DOWDING. Removal of muscular artifacts in EEG signals: a comparison of linear decomposition methods. Brain Informatics. 2018, vol. 6, pp. 13–22. ISSN 2198-4026. DOI: 10.1007/s40708-017-0074-6. [46] KONRAD, P. The ABC of EMG: A Practical Introduction to Kinesiological Electromyography. 1st ed. Scottsdale: Jones &Bartlett Learning, 2010. ISBN 0-977162-21-4. Available at: https: //www.noraxon.com/wp-content/ uploads/2014/12/ABC-EMG-ISBN.pdf. [47] CRISWELL, E. Cram’s Introduction to Surface Electromyography. 2nd ed. Massachusetts: Jones &Bartlett Learning, 2010. ISBN 978-0-76-3732745. [48] PEEBLES, P. Z. Probability, Random Variables, and Random Signal Principles. 4th ed. Singapore: McGraw-Hill, 2000. ISBN 978-0-07-118181-5. [49] BLINOWSKA, K. J. Review of the methods of determination of directed connectivity from multichannel data. Medical &Biological Engineering &Computing. 2011, vol. 49, iss. 5, pp. 521–529. ISSN 1741-0444. DOI: 10.1007/s11517-011-0739x. [50] BACCALA, L. A., D. Y. TAKAHASHI and K. SAMESHIMA. Directed Transfer Function: Unified Asymptotic Theory and Some of Its Implications. IEEE Transactions on Biomedical Engineering. 2016, vol. 63, iss. 12, pp. 2450–2460. ISSN 1558-2531. DOI: 10.1109/TBME.2016.2550199. [51] NAVIDI, W. C. Statistics for Engineers and Scientists. 3rd ed. New York: McGraw-Hill, 2011. ISBN 978-0-07-337633-2. [52] DJORDJEVIC, Z., A. JOVANOVIC and A. PEROVIC. Brain connectivity measure — the direct transfer function — advantages and weak points. 2012 IEEE 10th Jubilee International Symposium on Intelligent Systems and Informatics. Subotica: IEEE, 2012, pp. 93–97. ISBN 978-1-4673-4749-5. DOI: 10.1109/SISY.2012.6339493. [53] ANDERSON, C. W., E. A. STOLZ and S. SHAMSUNDER. Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks. IEEE Transactions on Biomedical Engineering. 1998, vol. 45, iss. 3, pp. 277–286. ISSN 1558-2531. DOI: 10.1109/10.661153. [54] SCHNEIDER, T. and A. NEUMAIER. Algorithm 808: ARfit—a matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Transactions on Mathematical Software. 2001, vol. 27, iss. 1, pp. 58–65. ISSN 0098-3500. DOI: 10.1145/382043.382316. [55] CHIANG, J., Z. J. WANG and M. J. MCKEOWN. Sparse multivariate autoregressive ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 607 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER (MAR)-based partial directed coherence (PDC) for electroencephalogram (EEG) analysis. In: 2009 IEEE International Conference on Acoustics, Speech and Signal Processing. Taipei: IEEE, 2009, pp. 457–460. ISSN 2379-190X. DOI: 10.1109/icassp.2009.4959619. [56] SMITH, S. W. The Scientist and Engineer’s Guide to Digital Signal Processing. 2nd ed. San Diego: California Technical Publishing, 1999. ISBN 09660176-7-6. [57] FERREIRA, J. L., R. M. AARTS and P. J. M. CLUITMANS. Optimized movingaverage filtering for gradient artefact correction during simultaneous EEG-fMRI. In: 5th ISSNIP-IEEE Biosignals and Biorobotics Conference (2014): Biosignals and Robotics for Better and Safer Living (BRC). Salvador: IEEE, 2014, pp. 1–6. ISBN 978-1-4799-5689-0. DOI: 10.1109/BRC.2014.6880955. [58] PROAKIS, J. G. and D. K. MANOLAKIS. Digital Signal Processing: Principles, Algorithms, and Applications. 3rd ed. New Jersey: Prentice Hall, 2007. ISBN 978-0-131-87374-2. [59] KAWALA-STERNIUK, A., M. PODPORA, M. PELC, M. BLASZCZYSZYN, E. J. GORZELANCZYK, R. MARTINEK and S. OZANA. Comparison of Smoothing Filters in Analysis of EEG Data for the Medical Diagnostics Purposes. Sensors. 2020, vol. 20, iss. 3, pp. 1–19. ISSN 1424-8220. DOI: 10.3390/s20030807. [60] PASQUINI, C., G. BOATO, N. ALAJLAN and F. G. B. DE NATALE. A Deterministic Approach to Detect Median Filtering in 1D Data. IEEE Transactions on Information Forensics and Security. IEEE, 2016, vol. 11, iss. 7, pp. 1425–1437. ISSN 1556-6021. DOI: 10.1109/TIFS.2016.2530636. [61] CORMEN, T. H., C. E. LEISERSON, R. L. RIVEST and C. STEIN. Introduction to Algorithms. 4th ed. Cambridge: MIT Press, 2022. ISBN 978-0-262-04630-5. [62] YADAV, S., J. PAWAR, G. PATIL and S. AGARWAL. EEG Signal Enhancement Using OWA Filter. International Conference on Automation, Computing and Communication (ICACC). 2021, vol. 40, iss. 1. pp.1–7. ISSN 2271-2097. DOI: 10.1051/itmconf/20214001010. [63] BAE, T. W., S. H. LEE and K. K. KWON. An Adaptive Median Filter Based on Sampling Rate for R-Peak Detection and Major-Arrhythmia Analysis. Sensors. 2020, vol. 20, iss. 21, pp. 1–21. ISSN 1424-8220. DOI: 10.3390/s20216144. [64] LIN, I. C., W. T. CHEN, Y. C. CHOU and P. Y. CHEN. A Novel Comparison-Free 1-D Median Filter. IEEE Transactions on Circuits and Systems II: Express Briefs. IEEE, 2020, vol. 67, iss. 7, pp. 1329–1333. ISSN 1558-3791. DOI: 10.1109/TCSII.2019.2936230. [65] LIN, I. C., W. T. CHEN, Y. C. CHOU and P. Y. CHEN. An Advanced EEG Motion Artifacts Eradication Algorithm. The Computer Journal. 2021, vol. 1, iss. 1. ISSN 0010-4620. DOI: 10.1093/comjnl/bxab170. [66] SUNDARAM, K., MARICHAMY and PRADEEPA. FPGA based filters for EEG pre-processing. 2016 Second International Conference on Science Technology Engineering and Management (ICONSTEM). Chennai: IEEE, 2016, pp. 572–576. ISBN 978-1-5090-1706-5. DOI: 10.1109/ICONSTEM.2016.7560958. [67] MARION, A. In: Introduction to Image Processing. 1st ed. New York: Springer, 1991. ISBN 9781-4899-3186-3. DOI: 10.1007/978-1-4899-3186-3. [68] SENDALL, J. L. and W. P. Du PLESSIS. Optimised Wiener Filtering in Overdetermined Systems. In: 2018 International Conference on Radar (RADAR). Brisbane: IEEE, 2018, pp. 1–6. ISBN 978-1-5386-7217-4. DOI: 10.1109/RADAR.2018.8557253. [69] SHAUGHNESSY, D. O. Linear predictive coding. IEEE Potentials. 1988, vol. 7, iss. 1, pp. 29–32. ISSN 1558-1772. DOI: 10.1109/45.1890. [70] MARIN-HURTADO, J. I. and D. V. ANDERSON. Reduced-bandwidth and low-complexity multichannel wiener filter for binaural hearing aids. In: IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA). New Paltz: IEEE, 2011, pp. 85–88. ISBN 978-14577-0693-6. DOI: 10.1109/ASPAA.2011.6082342. [71] SHUKLA, K. K. and A. K. TIWARI. Efficient algorithms for discrete wavelet transform: With applications to denoising and fuzzy inference systems. 1st ed. London: Springer, 2013. ISBN 978-1-4471-4941-5. [72] DALY, I., F. PICHIORRI, J. FALLER, V. KAISER, A. KREILINGER, R. SCHERER and G.MULLER-PUTZ. What does clean EEG look like? In: Annual International Conference of the IEEE Engineering in Medicine and Biology Society. San Diego: IEEE, 2012, pp. 3963–3966. ISBN 978-1-4577-17871. DOI: 10.1109/EMBC.2012.6346834. ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 608 BIOMEDICAL ENGINEERING VOLUME: 20 |NUMBER: 4 |2022 |DECEMBER About Authors Taha Mahmoud AL-NAIMI received his B.Sc. degree in Computer Communication Engineering from Al-Mansour University College in Baghdad, Iraq in 2015. He received his M.Eng. in Electronics and Telecommunication from Universiti Teknologi Malaysia (UTM) in 2018 and he is a Ph.D. student in UTM now. His research intersts include signal processing, signal analysis, biomedical signals, radar application, and microcontroller programming. Shanthini Chandra Sekaran NAIDU received her B.Ed. in Social Science from the teachers’ training college in Malaysia. She obtained her Master’s degree in Early Childhood Education from University Malaya. At present, she is pursuing her doctoral degree in Educational Psychology focusing on the reading acquisition among children with dyslexia and its relation to the brain functional connectivity. Ahmad Zuri SHA’AMERI obtained his B.Sc. in Electrical Engineering from the University of Missouri, Columbia, USA in 1984, and M.Eng. Electrical Engineering and Ph.D. both from UTM in 1991 and 2000 respectively. At present, he is an associate professor, Coordinator for the Digital Signal and Image Processing (DSIP) Research Group and Academic Coordinator for the DSP Lab, Electronic and Computer Engineering Department, Faculty of Electrical Engineering, UTM. His research interest includes signal theory, signal processing for radar and communication, signal analysis and classification, and information security. He has also conducted short courses for both government and private sectors. At present, he has published 165 papers in his areas of interest at both national and international levels in conferences and journals. Norlaili Mat SAFRI received the M.Eng. degree in electrical engineering from the Universiti Teknologi Malaysia, Malaysia, in 1999 and a Ph.D. degree in Systems and Information from the Kumamoto University, Japan, in 2006. She currently works as an associate professor at School of Electrical Engineering, Faculty of Engineering, Universiti Teknologi Malaysia. She has published over 85 refereed articles in international journals and conferences in her areas of interest. Her current research interests include the electrophysiological signal processing (EEG, EMG, ECG), human brain-muscle communication and control, cortical networks, brain-computer interface and cardiac arrhythmias detection. Narina Abu SAMAH is an associate professor at the School of Education, Faculty of Social Sciences and Humanities, Universiti Teknologi Malaysia (UTM) Johor Bahru. She obtained her Bachelor of Human Sciences in Psychology (Hons) from the International Islamic University Malaysia (IIUM) in 1997, and her Master of Human Sciences in Psychology from the same university in 2001. Dr. Narina received her Doctor of Philosophy (Ph.D.) from the University of Bristol, United Kingdom in 2011 and was later awarded a Graduate Certificate in Tertiary Education Management, from the University of Melbourne, Australia in 2015. Currently, she teaches educational foundation courses for undergraduate programmes, and cognitive psychology and qualitative research methodology for post-graduate programmes. Her research projects focused on the cognitive aspects of learning and teaching, which include critical reflective thinking, metacognition, and complex problem solving in learning, as well as brain functional connectivity during the learning process. ©2022 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING 609