Removal of ECG Artifacts Affects Respiratory Muscle Fatigue Detection-A Simulation Study

Abstract

This work investigates elimination methods for cardiogenic artifacts in respiratory surface electromyographic (sEMG) signals and compares their performance with respect to subsequent fatigue detection with different fatigue algorithms. The analysis is based on artificially constructed test signals featuring a clearly defined expected fatigue level. Test signals are additively constructed with different proportions from sEMG and electrocardiographic (ECG) signals. Cardiogenic artifacts are eliminated by high-pass filtering (HP), template subtraction (TS), a newly introduced two-step approach (TSWD) consisting of template subtraction and a wavelet-based damping step and a pure wavelet-based damping (DSO). Each method is additionally combined with the exclusion of QRS segments (gating). Fatigue is subsequently quantified with mean frequency (MNF), spectral moments ratio of order five (SMR5) and fuzzy approximate entropy (fApEn). Different combinations of artifact elimination methods and fatigue detection algorithms are tested with respect to their ability to deliver invariant results despite increasing ECG contamination. Both DSO and TSWD artifact elimination methods displayed promising results regarding the intermediate, "cleaned" EMG signal. However, only the TSWD method enabled superior results in the subsequent fatigue detection across different levels of artifact contamination and evaluation criteria. SMR5 could be determined as the best fatigue detection algorithm. This study proposes a signal processing chain to determine neuromuscular fatigue despite the presence of cardiogenic artifacts. The results furthermore underline the importance of selecting a combination of algorithms that play well together to remove cardiogenic artifacts and to detect fatigue. This investigation provides guidance for clinical studies to select optimal signal processing to detect fatigue from respiratory sEMG signals.

Keywords: ECG; EMG; MNF; SMR; biomedical signal processing; cardiogenic artifacts; fApEn; fatigue detection; muscle fatigue; neuromuscular fatigue; respiration; respiratory EMG; sEMG.

Figure A1

ECG preprocessing and epoch construction. (A) displays a single QRS complex of baseline-filtered ${\mathrm{ECG}}_{\mathrm{I}}^{\mathrm{BLF}}\left(t\right)$ and ${\mathrm{ECG}}_{\mathrm{I}}^{\mathrm{SG}}\left(t\right)$ pre-processed by the application of a Savitzky–Golay filter. (B) shows the 60% MVC sEMG signal taken from a previous study [16]. (C) shows the artificial test signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}\left(t\right)$. The detected and time-refined position of the QRS is indicated by the orange diamond. Ranges of the signal close to the QRS complex designated to be excluded from the following fatigue analysis are drawn in purple. (D,E) display the epoch construction based on ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}^{\mathrm{TSW}15}\left(t\right)$ obtained by template subtraction. A fatigue value is calculated every 125 $\mathrm{m}$$\mathrm{s}$ (symbolized by dots), yielding a fatigue signal with 8 $\mathrm{Hz}$ sample rate. The epoch that forms the basis of each fatigue value is indicated by the timeline to the left of the dot at the same height. As an example, the fatigue value assigned to time $t=35.2\mathrm{s}$ (as indicated by the red dot) is based on the epoch drawn in blue. This example uses an epoch size of a $0.25$ $\mathrm{s}$ (256 samples). Subsegments as used in the Welch PSD estimation, with $k_w=7$ subsegments shown in green for the exemplary epoch (shown in blue). (D) visualizes the gating strategy GN, which does not exclude former QRS segments. The fatigue evaluated segment is composed of the last 256 samples regardless of remnants of a QRS complex. The gating strategy GO displayed in (E) skips the QRS remnants.

Figure A2

An evaluation of heartbeat irregularity. (A) visualizes three characteristic properties evaluated for each heart cycle. The evaluation is performed for both ${\mathrm{ECG}}_{\mathrm{I}}\left(t\right)$ and ${\mathrm{ECG}}_{\mathrm{A}}\left(t\right)$. The resulting distributions of these characteristic properties for both ECG channels I and A are shown in (B) for each subject. ECG channel A is artificially mixed and features increased heartbeat irregularity. The violin plot visualizes the variation of the characteristic property within the same channel of one subject. The subject number is stated on the horizontal axis. An alternative way to evaluate the heartbeat irregularity is based on the amount of cardiogenic artifacts that could not be eliminated by the artifact removal step and is still present in ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right)$. The analysis is preferably used in combination with template subtraction TSW15 as the artifact removal method. A rectified segment of ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right)$ is extracted for each heart cycle, as shown in (C). These segments are superimposed by the same QRS times used for template subtraction. A raw version of an average remaining cardiogenic artifact signal ${\mathrm{rARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left({\tau }_c\right)$ is constructed by averaging all segments at the same relative heart cycle time ${\tau }_c$. The average remaining cardiogenic artifact signal ${\mathrm{ARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left({\tau }_c\right)$ is yielded by the application of a 25 tap Blackman smoothing filter to the raw version ${\mathrm{rARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left({\tau }_c\right)$. Subfigures (D,E) visualize the average remaining cardiogenic artifact signals ${\mathrm{ARCA}}_{60,0.02,\mathrm{ch}}^{\mathrm{cr}}\left({\tau }_c\right)$ for different artifact removal methods and ECG channels. The subfigures differ only in the range of the vertical axis. The peak value for each ${\mathrm{ARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left({\tau }_c\right)$ is denoted as ${\mathrm{pARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}$. The ratio of the peak average of remaining cardiogenic artifacts to the peak average of cardiogenic artifacts (based on signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}\left(t\right)$) is used to quantify the heartbeat irregularity. The pARCA ratio is calculated as $\frac{{\mathrm{pARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}}{{\mathrm{pARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{NR}}}$. The artifact removal method NR denotes no treatment of cardiogenic artifacts. pARCA ratios for both ECG channels I and A are shown in (F) for each subject. These ratios are based on template subtraction TSW15 as the artifact removal method and MVC load level $\Lambda =60$. Each violin plot contains pARCA ratios from EMG levels $\eta =0.01,0.02,0.05\mathrm{and}0.1$. For all but one subject, the pARCA ratio indicates a distinctively elevated heartbeat irregularity for EMG channel A.

Figure A3

A block diagram visualizing the test signal generation. See Section 2.1 for the corresponding description.

Figure A4

Table stating the utilized combinations of EMG recordings (differential signals originating from electrodes 1 and 3—subject numbers from [16]) and ECG recordings (identifiers from [20]). Further indicated is the number of heartbeats (# HB) within the analyzed segment of pure ECG. False positive (FP) and false negative (FN) detected QRS complexes from signals ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{A}}\left(t\right)$ are stated for different EMG levels $\eta$ and MVC loads $\Lambda$.

Figure A5

Block diagram displaying different methods to reduce cardiogenic artifacts. See Section 2.2 for a detailed description. For clarity reasons, the TSWD algorithm is only depicted for a single wavelet scale and suggested for two other instead of the utilized 8 scales. Not depicted are the methods TSWD35 and DSO. Both are variants to TSWD15 and differ in the signal that is fed into the stationary wavelet transform (SWT). TSWD35 uses the signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{TSW}35}$ as input. In the case of DSO, the unpreprocessed signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}$ is used as input. The ECG removal methods featuring a high pass with $f_c=35\mathrm{Hz}$ are only relevant for fApEn fatigue detection. PSD-based fatigue detection relies solely on the 15 $\mathrm{Hz}$ version as the PSD lower bound $b_l$ is applied.

Figure A6

(A) displays the template that is subtracted QRS synchronously in the first step. To limit its influence on the current cardiac cycle, values outside the average cardiac cycle length are set to zero, being the additive identity. (B) visualizes the average remaining cardiogenic artifact signal ${\mathrm{wARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{TSW}15\_\mathrm{SWT},i}\left({\tau }_c\right)$ for different, color-coded wavelet levels i. It is constructed by averaging all rectified heart cycle segments of the signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{TSW}15\_\mathrm{SWT},i}\left(t\right)$ at the same relative heart cycle time ${\tau }_c$. The signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{TSW}15\_\mathrm{SWT},i}\left(t\right)$ is obtained by the application of the stationary wavelet transform with eight different levels to the result of the template subtraction from the first step. Dash-dotted lines indicate the threshold ${\theta }_i$. Damping is carried out only at heart cycle times ${\tau }_c$ where the average remaining cardiogenic artifact signal ${\mathrm{wARCA}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{TSW}15\_\mathrm{SWT},i}\left({\tau }_c\right)$ is above the threshold ${\theta }_i$. (C) shows the resulting damping templates $\mathrm{dTemp}\left({\tau }_c,i\right)$ for different, color-coded wavelet levels i. Again, values outside the average cardiac cycle length are set to identity, in this case, 1 due to the multiplication. The remaining subpanels compare the results of different template subtraction methods. (D) displays a segment of ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}$. (E–G) display the corresponding ECG-reduced segments obtained with artifact removal methods DSO, TSW15 and TSWD15, respectively. (H) displays the pure EMG component ${\mathrm{RSM}}_{\Lambda }$ used in the construction of the artificial test signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}$. This signal can be considered as a target value for artifact removal. The example signals illustrate that TSWD can best remove the artifacts. While the pure application of the damping step (DSO) removes too much of the EMG signal in the QRS range, the ordinary template subtraction (TSW) is not able to remove the ECG component completely, and artifacts remain in the QRS range.

Figure A7

Block diagram of applied fatigue detection methods. See Section 2.3 for a detailed description. Fatigue detection is performed based on signals ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr}}$. These signals feature reduced cardiogenic artifacts produced with the cardiogenic artifact elimination method cr. The resulting fatigue index signals are denoted ${\Phi }_{\Lambda ,\eta ,\mathrm{ch}}^{\mathrm{cr},\mathrm{gs},N_e,\mathrm{fd}}$, where gs denotes the gating strategy, $N_e$ the epoch size and fd the fatigue detection method. The prefix se_ in front of MNF and SMR5 fatigue detection methods denotes the applied spectral estimation method (Welch or Burg method).

Figure A8

Detailed description of fatigue detection algorithms and exemplary fatigue detection signals. The mean frequency and spectral moments ratio depend on an estimation of the power spectral density (PSD) of the underlying epoch. The PSD was calculated in a nonparametric way by the Welch method [37]. Epochs are subdivided into subsegments with a quarter, an eighth or a sixteenth of the epoch size. Subsegments had a 50% overlap and were tapered by a Hamming window function. This led to $k_w=7,15\mathrm{or}31$ subsegments for each epoch. In the case of the gating strategy GO, all subsegments affected by gating are excluded from the averaging in the Welch method. Alternatively, a parametric estimate was calculated by the Burg algorithm [37]. The Burg algorithm is performed individually for all fragments separated by gating in the case of the gating strategy GO. A final PSD estimation is obtained as the weighted average of PSDs from all fragments. The weight is chosen according to the length of the fragment in relation to the total length of all fragments. Fragments shorter than a quarter of the intended epoch size are excluded. In practice, the presented cardiogenic artifact removal will only suppress a part of the cardiogenic artifacts. Furthermore, the ECG power spectrum tends to have a larger share in the low-frequency range than the EMG power spectrum. To reflect these two issues, we introduced a PSD lower bound $b_l$. Both PSD-based fatigue detection algorithms only considered PSD bins in between the lower bound $b_l$ and upper bound $b_u$. Fatigue detection was performed with the lower bound $b_l$ at 15, 20, 25, 30, 35, 40, 45 and 50 $\mathrm{Hz}$. In contrast, $b_u$ was fixed at 500 $\mathrm{Hz}$ close to the Nyquist frequency. To facilitate the introduction of $b_l$, all PSD-based fatigue calculations rely on 15 $\mathrm{Hz}$ high-pass filtered signals. The mean frequency (MNF) fatigue identifier [34] is calculated as: $\mathrm{MNF}=\raisebox{1ex}{${\int }_{b_l}^{b_u}f·\mathrm{PSD}\left(f\right)df$}\!\left/ \!\raisebox{-1ex}{${\int }_{b_l}^{b_u}\mathrm{PSD}\left(f\right)df$}\right.$. The spectral moments ratio (SMR) suggested by Dimitrov et al. [35] tries to better incorporate fatigue-related PSD changes at higher frequencies. This is accomplished by multiplying the PSD with a weighting function that amplifies higher PSD bins. The spectral moment M of order p also follows the PSD bounds $b_l$ and $b_u$ and is calculated as $M_p={\int }_{b_l}^{b_u}{f}^p·\mathrm{PSD}\left(f\right)df$. The spectral moments ratio fatigue index of order p is defined as ${\mathrm{SMR}}_p=ln\frac{M_{-1}}{M_p}$. We used an order of $p=5$ and therefore choose the abbreviation SMR5 in this study. The natural logarithm is utilized to prevent a range over several orders of magnitude and ensures comparable values for the following assessment [17]. Fuzzy approximate entropy (fApEn) was introduced by Xie et al. [36]. As an entropy method, it compares one short sequence from the epoch to all other short sequences from the same epoch. The entropy is based on the number of similar sequences of length m compared to the amount of similar sequences with length $m+1$. In contrast to other entropy methods, fApEn altered the decision of whether two sequences are considered similar. The former binary decision was exchanged for a fuzzy similarity $u(d,r)=exp\left(\frac{-d^2}{r}\right)$. d denotes the maximum metric of the difference of both sequences. The method also considers two sequences to be similar if they have an arbitrary constant offset. We used a threshold of $r=0.6$ and a vector length $m=2$. All fragments separated by gating were concatenated in the case of the gating strategy GO for fApEn. In this case, the fApEn method is slightly adapted to take into account the concatenation of fragments. This concerns, on the one hand, the zero-mean: The mean is calculated and compensated separately within each fragment. Furthermore, those short sequences excluded from calculation that are extending over more than one fragment. These short sequences contain the fracture at the concatenation and might disturb the entropy calculation. Subfigure (A) visualizes the signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}^{\mathrm{TSWD}15}\left(t\right)$. Derived are the fatigue index signals with MNF, SMR5 and fApEn fatigue detection methods shown in (B–D), respectively. All fatigue index signals are based on epoch size $N_e=256$. The PSD for MNF was estimated with the Welch method with $k_w=15$, while the Burg method was used for SMR5. SMR5 and fApEn utilized gating and no gating was applied for MNF. All three variants used a PSD lower bound $b_l=35\mathrm{Hz}$ (in the case of fApEn, a corresponding high-pass is applied). (E) displays the corresponding normalized fatigue signals ${\widehat{\Phi }}_{60,0.05,\mathrm{A}}^{\mathrm{TSWD},\mathrm{gs},256,\mathrm{fd}}\left(t\right)$ that are converted to a zero to one range. (F) shows the corresponding fatigue signals that are normalized and filtered to incorporate information from a 1 $\mathrm{s}$ segment of the EMG signal at each fatigue index value. They are denoted ${\widehat{\Phi }}_{60,0.05,\mathrm{A}}^{\mathrm{TSWD},\mathrm{gs},256,\mathrm{fd},\mathrm{MAfilt}}\left(t\right)$.

Figure A9

Temporal and spectral evaluation of ECG “cleaned” signals. The evaluation is based on the different error variants quantifying the deviation of the ECG “cleaned” signal ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}\left(t\right)$ compared to the scaled EMG component ${\mathrm{SMC}}_{\Lambda ,\eta }\left(t\right)$. The properly scaled, pure ECG component is defined as ${\mathrm{SMC}}_{\Lambda ,\eta }\left(t\right)={\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{ch}}\left(t\right)-{\mathrm{ECG}}_{\mathrm{ch}}\left(t\right)=\frac{\eta }{0.2}·{\mathrm{RSM}}_{\Lambda }\left(t\right)$. The squared error of the ECG attenuated signal with removal method “cr” is defined as: $e_{\mathrm{sq}}\left(\eta ,\mathrm{ch},\mathrm{cr}\right)=\frac{\parallel {\mathrm{SMC}}_{60,\eta }\left(t\right)-{\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right){\parallel }_2}{\parallel {\mathrm{SMC}}_{60,\eta }\left(t\right){\parallel }_2}{|}_{0\mathrm{s}<t<60\mathrm{s}}$. This corresponds to the definition of the error for synthetic signals by Petersen et al. [15]. Subpanels (A,B) show the pure EMG component ${\mathrm{SMC}}_{60,0.05}\left(t\right)$ and the ECG “cleaned” signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}^{\mathrm{TSWD}}\left(t\right)$, respectively. The mean removal $e_{\mathrm{sq}}\left(\eta ,\mathrm{A},\mathrm{cr}\right)$ error across all nine subjects is displayed in (C) for different ECG removal logarithms and EMG levels $\eta$. The influence of ECG removal variants on a subsequently calculated envelope signal is quantified by an error $e_{\mathrm{env}}$, as known from Petersen et al. [15]. The mean average value MAV with a centralized sliding window (size 128) is calculated for both the pure EMG component and the cleaned signal. Exemplary envelope signals are shown in (D). The envelope error is defined as: $e_{\mathrm{env}}\left(\eta ,\mathrm{ch},\mathrm{cr}\right)=\frac{\parallel \mathrm{MAV}\left({\mathrm{SMC}}_{60,\eta }\left(t\right)\right)-{\beta }_1^{*}\mathrm{MAV}\left({\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right)\right)-{\beta }_0^{*}{\parallel }_2}{\parallel \mathrm{MAV}\left({\mathrm{SMC}}_{60,\eta }\left(t\right)\right){\parallel }_2}{|}_{0\mathrm{s}<t<60\mathrm{s}}$. With ${\beta }_1^{*}$ and ${\beta }_0^{*}$, a linear factor and an offset are applied to scale and shift the envelope of the cleaned signal to minimize the error. Both values are obtained by a linear regression: ${\beta }_1^{*}{\beta }_0^{*}={argmin}_{{\beta }_0,{\beta }_1}{\parallel \mathrm{MAV}\left({\mathrm{SMC}}_{60,\eta }\left(t\right)\right)-{\beta }_1\mathrm{MAV}\left({\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right)\right)-{\beta }_0\parallel }_2{|}_{0\mathrm{s}<t<60\mathrm{s}}$ The resulting mean error $e_{\mathrm{env}}\left(\eta ,\mathrm{A},\mathrm{cr}\right)$ across all nine subjects is displayed in (E). Spectral properties of the ECG “cleaned” signal are evaluated based on the power spectral density (PSD) (F). The PSD was calculated by the Burg method with a 60 s long segment. Again, a quadratic error was calculated to quantify the deviation of the spectrum: $e_{\mathrm{spec}}\left(\eta ,\mathrm{ch},\mathrm{cr}\right)=\frac{\parallel \mathrm{PSD}\left({\mathrm{SMC}}_{60,\eta }\left(t\right),f\right)-\mathrm{PSD}\left({\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right),f\right){\parallel }_2}{\parallel \mathrm{PSD}\left({\mathrm{SMC}}_{60,\eta }\left(t\right),f\right){\parallel }_2}{|}_{0\mathrm{s}<t<60\mathrm{s},35\mathrm{Hz}<f<512\mathrm{Hz}}$. To exclude effects of the high-pass filter, only frequencies above 35 $\mathrm{Hz}$ are considered. To account for fatigue algorithms that mainly look for spectral changes instead of absolute values of the PSD, an alternative error $e_{\mathrm{SpecShape}}$ is calculated, which allows for a vertical shift and scaling of the PSD of the ECG “cleaned” signal in a similar manner as with $e_{\mathrm{env}}$. It is defined as: $e_{\mathrm{SpecShape}}\left(\eta ,\mathrm{ch},\mathrm{cr}\right)=\frac{\parallel \mathrm{PSD}\left({\mathrm{SMC}}_{60,\eta }\left(t\right),f\right)-{\beta }_3^{*}\mathrm{PSD}\left({\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right),f\right)-{\beta }_2^{*}{\parallel }_2}{\parallel \mathrm{PSD}\left({\mathrm{SMC}}_{60,\eta }\left(t\right),f\right){\parallel }_2}{|}_{0\mathrm{s}<t<60\mathrm{s},35\mathrm{Hz}<f<512\mathrm{Hz}}$. Both values are again obtained by a linear regression: ${\beta }_3^{*}{\beta }_2^{*}={argmin}_{{\beta }_2,{\beta }_3}{\parallel \mathrm{PSD}\left({\mathrm{SMC}}_{60,\eta }\left(t\right),f\right)-{\beta }_3\mathrm{PSD}\left({\mathrm{ATS}}_{60,\eta ,\mathrm{ch}}^{\mathrm{cr}}\left(t\right),f\right)-{\beta }_2\parallel }_2{|}_{0\mathrm{s}<t<60\mathrm{s},35\mathrm{Hz}<f<512\mathrm{Hz}}$. The resulting errors $e_{\mathrm{spec}}\left(\eta ,\mathrm{A},\mathrm{cr}\right)$ and $e_{\mathrm{SpecShape}}\left(\eta ,\mathrm{A},\mathrm{cr}\right)$ are shown in (G,H), respectively. $e_{\mathrm{spec}}\left(\eta ,\mathrm{A},\mathrm{cr}\right)$ for HP15 and no ECG removal at all are so large, that they are not displayed in (H).

Figure A10

A comparison of the best-working ECG removal methods for PSD-based fatigue detection algorithms. Depicted are template subtraction without (TSW15) and with (TSWD15) subsequent wavelet-based damping step. Both methods are shown as variants without (GN) and with gating (GO). The utilized fatigue detection algorithms include MNF with Welch PSD estimation ($k_w=15$) on the left side and SMR5 in combination with Burg spectral estimation on the right side. Both fatigue detection methods are based on an epoch size of $N_e=256$ and incorporate a PSD lower bound $b_l$ of 35 $\mathrm{Hz}$. The white dots indicate the median value.

Figure A11

The left side displays the influence of the epoch size and gating strategy in case of entropy-based fatigue detection with fApEn. The main blocks on the horizontal axis denote different EMG levels $\eta$ again. For each of these main blocks, different epoch sizes are shown on the horizontal axis. To assure a balanced comparison of epoch sizes, results from smaller epochs are filtered the same way as described in Figure A12. The applied gating strategy is denoted by the marker and line style. The three applied evaluation criteria (A, B, C) are again covered by the rows of subfigures. Each data point is obtained by averaging over all subjects. The right side compares the best working ECG removal methods for fApEn fatigue detection. Depicted are template subtraction without (TSW35) and with (TSWD35) subsequent wavelet-based damping step. Both methods are shown as variants without (GN) and with gating (GO). The results are based on an epoch size of $N_e=256$. In contrast to the left side, no filtering of the fatigue index to incorporate one second of EMG was utilized. Because there was no filter, each fatigue index value incorporates only a quarter second of EMG, and the results are therefore not directly comparable to the left side. The white dots indicate the median value.

Figure A12

The influence of epoch size $N_e$ in the case of PSD-based fatigue detection for different combinations of ECG removal method, gating strategy and PSD estimation method. The main blocks on the horizontal axis denote different EMG levels $\eta$. Connected data points denote the same EMG level. For each of these main blocks, different epoch sizes are shown on the horizontal axis. To assure a balanced comparison of epoch sizes, the results from smaller epochs underlying fatigue detection signals $\Phi \left(t\right)$ are filtered. The length of the moving average filter is chosen in a way that each filtered fatigue index data point incorporates EMG from a 1s long window. The results displayed are based on TSWD15 in combination without gating (GN) as well as TSW15 in combination with the gating strategy GO. In the case of pure EMG, a sole high-pass filter (HP15) is applied instead of template subtraction. The applied PSD estimation is color-coded. The Burg method is drawn in red, while Welch methods with $k_w=$ 7, 15 and 31 subsegments are shown in green, blue and purple, respectively. The number of values in each Welch subsegment depends on the epoch size. The three applied evaluation criteria (A, B, C) are again covered by the rows of subfigure, while MNF and SMR5 represent the columns. Each data point is obtained by averaging over all subjects.

Figure A13

The influence of PSD lower bound $b_l$. The applied PSD lower bound $b_l$ is plotted on the horizontal axis. Only PSD bins corresponding to frequencies above $b_l$ are considered in the subsequent calculations. Different EMG levels $\eta$ are color-coded. The results obtained with QRS segments excluded (gating strategy GO) are shown by dashed lines. The three applied evaluation criteria (A, B, C) are again covered by the rows of subfigure, while MNF and SMR5 represent the columns. The PSD for MNF is again estimated by the Welch method with 15 subsegments, while SMR5 relies on the Burg method. The epoch size is 256 values. Each data point is obtained by averaging over all subjects. Due to the limited frequency resolution of the Welch method with 15 subsegments, only two different results are obtained for each EMG level. At higher EMG levels $\eta$, a higher PSD lower bound $b_l$ has a negative effect in terms of the introduced disturbances measured by $R^2$ (criterion C). However, it has hardly any effect in respect to the deviation from the reference fatigue index signal (criterion A) and the separability of load cases (criterion B). Regarding lower EMG levels $\eta$, the deviation to the reference fatigue index signal and the separability intensely improved towards higher $b_l$. In the case of disturbances (criterion C), the results displayed an optimal choice at values around 30 $\mathrm{Hz}$.

Figure A14

A comparison with intra-heartbeat fatigue detection, as suggested by Sinderby et al. [10]. In the case of the intra-heartbeat approach, only a segment extending between 50% and 75% of the RR interval is utilized. Instead of fixed length epochs, this results in a variable epoch size dependent on the heart rate. Prior to fatigue detection, the resulting segment is furthermore detrended. To ensure a fair comparison, the different epoch sizes must be considered. Suitable to the applied epoch size of $N_e=256$, a moving average of the last 7 fatigue index values was applied. This assures that EMG segments of one second are incorporated. This corresponds to the procedure that was used when comparing the length of epochs (Figure A12). The intra-heartbeat approach generates one fatigue value per heart cycle. This value is assigned to the end time of the intra-heartbeat EMG segment. With the help of linear interpolation, these values are mapped to the fixed 8 $\mathrm{Hz}$ fatigue sampling times. Again, the preceding 7 values are moving average filtered to obtain a filtered fatigue index signal for comparison. Due to the interpolation preceding the filtered intra-heartbeat fatigue signal, sometimes it includes EMG from beyond a strict one second interval. Despite the fact that the filtered intra-heartbeat fatigue signal sometimes include EMG from beyond a strict one second interval, its results (purple colored) are generally not better. In most combinations of EMG levels and evaluation criteria, the intra-heartbeat results scored worse compared to our TSW15 with the gating strategy GO and TSWD15 without gating approaches. This applied for both MNF (left column) as well as SMR5 (right column) fatigue detection algorithms. The black dots indicate the mean value, while the white dots symbolize the median.

Figure 1

Artificial signal construction based on electrocardiographic (ECG) and surface electromyographic (sEMG) signals. (A) shows an sEMG signal from a 60% MVC load case (${\mathrm{EMG}}_{60}\left(t\right)$) taken from [16]. Modulation with the weighting signal ${\mathrm{MS}}_{\mathrm{M}}\left(t\right)$ (displayed in (B)) to mimic a breathing-like activity pattern yields ${\mathrm{EMG}}_{60}^{\mathrm{M}}\left(t\right)$ displayed in (C). (D,E) display ECG signals ${\mathrm{ECG}}_{\mathrm{I}}^{\mathrm{SG}}\left(t\right)$ and ${\mathrm{ECG}}_{\mathrm{II}}^{\mathrm{SG}}\left(t\right)$ from channel I and II, respectively. (F) shows the weighting signal ${\mathrm{WS}}_{\mathrm{C}}\left(t\right)$ roughly corresponding to the lung volume utilized to mix both ECG channels. Its application yields ${\mathrm{ECG}}_{\mathrm{A}}\left(t\right)$ with an increased irregularity displayed in (G). The artificial test signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}\left(t\right)$ based on ${\mathrm{ECG}}_{\mathrm{A}}\left(t\right)$ and ${\mathrm{EMG}}_{60}^{\mathrm{M}}\left(t\right)$ is drawn in (H). Note the small amplitude of sEMG components. (I) shows a “cleaned” signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}^{\mathrm{TSW}15}\left(t\right)$ with cardiogenic artifacts removed by template subtraction. Note the change in scale as compared to (H). Some cardiogenic artifacts are still clearly visible in the signal due to the intended ECG irregularity of the artificial signal. (J) shows a signal ${\mathrm{ATS}}_{60,0.05,\mathrm{A}}^{\mathrm{TSWD}15}\left(t\right)$ with further reduced residual cardiogenic artifacts. All physiological signals are displayed in $\mathsf{\mu }$V.

Figure 2

Different EMG levels mixed with ECG. The first main row (A) depicts the reference case, including the sEMG modulation but no ECG contamination. The following main rows (B–D) depict the cases with $\eta =0.01,0.02\mathrm{and}0.05$, respectively. For each situation, the 20% and 60% maximum voluntary contraction (MVC) load cases are shown in the first and second sub row. Column (I) shows the reference signal ${\mathrm{RSM}}_{\Lambda }\left(t\right)$ in the case of row (A) and artificial test signals ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{A}}\left(t\right)$ otherwise. The signals with cardiogenic artifacts removed by template subtraction ${\mathrm{ATS}}_{\Lambda ,\eta ,\mathrm{A}}^{\mathrm{TSW}15}\left(t\right)$ are shown in the column (II). Corresponding fatigue index signals and their evaluation is visualized in Figure 3.

Figure 3

Exemplary fatigue index signals and their evaluation corresponding to exemplary test signals depicted in Figure 2. The first row represents the reference case, while the following rows visualize cases with $\eta =0.01,0.02\mathrm{and}0.05$, respectively. Column (III) shows the resulting filtered fatigue index signals for the 20% and 60% MVC load cases in green and red, respectively. Note that each filtered fatigue index data point incorporates sEMG from a 1s long window. The mean frequency (MNF) is used as fatigue algorithm with an epoch size of $N_e=256$, power spectral density (PSD) estimation with the Welch method ($k_w=15$) and a lower bound $b_l=35\mathrm{Hz}$. Column (IV) shows the distribution of values $\overline{\Phi }$ arising from the last 15 $\mathrm{s}$ of both load cases. The difference to one of the overlapping area of both distributions (light blue) corresponds to the Kolmugorov–Smirnov distance used to quantify the separability of fatiguing and non-fatigued situations by ${\gamma }_B$. The resulting performance evaluation values ${\gamma }_A$, ${\gamma }_B$ and ${\gamma }_C$ for deviation to the reference fatigue index (■), separability of load cases (●) and disturbances measured by $R^2$ (▲), respectively, are shown in column (V). Values arising from the reference case are drawn in light gray for all other EMG levels.

Figure 4

Fatigue detection results depend on the magnitude and characteristics of the cardiogenic artifacts. This figure juxtaposes the results originating from ECG channel I (${\Phi }_{\Lambda ,\eta ,\mathrm{I}}^{\mathrm{TSW}15,\mathrm{GN},256,\mathrm{fd}}\left(t\right)$) and the artificially mixed ECG channel A (${\Phi }_{\Lambda ,\eta ,\mathrm{A}}^{\mathrm{TSW}15,\mathrm{GN},256,\mathrm{fd}}\left(t\right)$). Each subfigure includes different EMG levels $\eta$ plotted on the horizontal axis. Additionally, two cases without ECG contamination are included for reference. The first uncontaminated case contains breathing phase modulation ${\Xi }_{\Lambda }^{\mathrm{GN},256,\mathrm{fd}}\left(t\right)$ (based on ${\mathrm{RSM}}_{\Lambda }\left(t\right)$), while the later ${\mathsf{\Psi }}_{\Lambda }^{\mathrm{GN},256,\mathrm{fd}}\left(t\right)$ (based on ${\mathrm{RS}}_{\Lambda }\left(t\right)$) does not. Each row resembles one of the applied evaluation criteria. The discrepancy of fatigue indexes originating from the contaminated signals in relation to the reference fatigue index with modulated EMG only (quantified by ${\gamma }_A$) is shown in the first row (A,D). The smaller the deviation, the better the fatigue algorithm’s performance. The second row (B,E) visualizes the separability of the 20% and 60% MVC load cases (quantified by ${\gamma }_B$, the higher the better). The level of disturbances measured by the coefficient of determination $R^2$ (quantified by ${\gamma }_C$) is depicted in the last row (C,F). The higher $R^2$, the better the fatigue detection method. The template subtraction method TSW15 was applied to compensate cardiogenic artifacts. The subfigures on the left are based on MNF with Welch PSD estimation ($k_w=15$). The subfigures to the right are based on spectral moments ratio of order five (SMR5) with Burg PSD estimation. Both fatigue detection methods are based on an epoch size $N_e=256$ and utilized a PSD lower bound $b_l=15\mathrm{Hz}$. The black dots indicate the mean value, while the white/orange dots symbolize the median.

Figure 5

The effect of cardiogenic artifact removal on the subsequent fatigue detection. Different methods to reduce the cardiogenic artifacts are compared next to omitting any artifact treatment (gray). Included is the sole application of a high-pass filter (HP15 with $f_c=15\mathrm{Hz}$) and two different variants of template subtraction (TS15 and TSW15). The latter template subtraction utilizes the additional wavelet denoising of the template. Additionally shown is heartbeat synchronous damping in different wavelet scales applied solely (DSO) and subsequent to template subtraction (TSWD15). The results considering gating are shown for each method separately (dashed line-gating strategy GO). In addition to the different EMG levels $\eta$, each figure also includes two cases without any ECG contamination. The first case (gray filled) evaluates ${\Xi }_{\Lambda }^{\mathrm{GN},256,\mathrm{fd}}\left(t\right)$, while the second one evaluates ${\Xi }_{\Lambda }^{\mathrm{GO},256,\mathrm{fd}}\left(t\right)$. The second case (gray shaded) introduces a gating based on time information taken over from pure ECG signal. It shows the impact of causeless gating (despite cardiogenic artifacts lacking). This allows an assessment of how much the fatigue detection quality drops by the sole interruptions due to gating. The three applied evaluation criteria are again covered by the rows of subfigures ((A,D), (B,E) and (C,F), respectively), while MNF and SMR5 represent the columns. The Welch method ($k_w=15$) is utilized for MNF, while SMR5 relies on the Burg method. $N_e=256$ and $b_l=35\mathrm{Hz}$ are applied for both methods. Each data point is obtained by averaging over all subjects. A more detailed comparison for a selection of ECG removal algorithms in combination with SMR5 and no gating (corresponding to the solid lines of the right column) can be found in Figure 6. Further comparisons for the algorithms TSW15 and TSWD15 can be found in the Appendix A Figure A10 and Figure A11).

Figure 6

A detailed comparison of selected cardiogenic artifact removal methods regarding their impact on the subsequent fatigue detection by SMR5 method and no gating. The shown violin plots correspond to solid lines in the same color from the right column of Figure 5. The best performing ECG removal method was determined for each group of four violin plots sharing the same EMG level and evaluation category. The existence of significant differences to the other three algorithms were evaluated for each group. If significant differences exist, the corresponding p-values are stated. SMR5 is calculated with $N_e=256$ and $b_l=35\mathrm{Hz}$ and utilized the Burg method for spectral estimation.

Figure 7

A comparison of fatigue detection algorithms. The results of MNF, SMR5 and fuzzy approximate entropy (fApEn) methods are shown in combination with ECG removal by TSWD and no gating (solid lines) as well as TSW and gating (dashed dotted lines). MNF is utilized with Welch ($k_w=15$) and SMR5 with Burg PSD. Both PSD-based algorithms include a PSD lower bound $b_l$ of 35 $\mathrm{Hz}$, while fApEn includes a 35 $\mathrm{Hz}$ high-pass filter within the ECG removal step. All fatigue detection methods are based on epoch size $N_e=256$. The results are shown as the mean value obtained by averaging over all subjects. Note that the distributions underlying the shown mean values are visualized in Figure 8.

Figure 8

A detailed comparison of fatigue detection algorithms. The left column displays the results based on TSWD artifact removal method (including the wavelet-based damping step) without QRS gating (corresponding to solid lines in Figure 7). The right visualizes the results based on template subtraction TSW15 with QRS gating (GO), corresponding to dashed dotted lines in Figure 7. Fatigue algorithms include MNF, SMR5 and fApEn with the same parameters as in Figure 7. The best-performing fatigue detection method was determined for each group of three violin plots sharing the same EMG level and evaluation category. The existence of significant differences to the other two methods were evaluated for each group. If significant differences were found, the corresponding p-values are stated. The white dots symbolize the median value. The horizontal axis denotes EMG levels $\eta$. The three applied evaluation criteria (A, B, C) are covered by the rows.