A Denoising Algorithm Based on Spatial and Spectral Domain Filtering for Electrically Evoked Compound Action Potential Enhancement

1. Introduction

The electrically evoked compound action potential (ECAP) is the combined response of auditory nerve fibers to electrical stimulation by a cochlear implant (CI) [1,2]. The ECAP is a critical signal for clinicians to evaluate the functionality of a patient's auditory nerve fibers after CI surgery [3,4]. In the absence of feedback from CI users, the ECAP model can be of great value to clinicians assessing the hearing performance of CI users [5,6]. ECAP amplitude can be computerized and estimated using temporal or spatial methods [7,8,9]. The ECAP matrix can be constructed by measuring the ECAPs generated by the different electrode positions of the masker and probe stimuli. Such a method is referred to as the forward-masking technique [10]. The panoramic ECAP (PECAP) method [11,12] in the literature has shown that using a clean ECAP matrix can approximate auditory activity patterns and ECAP magnitudes [11,12], further assisting clinicians in evaluating a patient's speech perception. Nevertheless, the use of noisy ECAP matrices, particularly when the signal-to-noise ratio (SNR) is less than 10 dB, has the potential to introduce inaccuracies in the estimation of auditory activity patterns and ECAP magnitudes [11,12].

The aforementioned reason necessitates the implementation of noise reduction techniques on ECAP matrices prior to estimating the auditory activity patterns and ECAP magnitudes. These techniques can be classified into three primary categories: spatial filtering, temporal filtering, and spectral filtering. In the context of image signal processing, the mean and median filters represent well-established examples of spatial filtering [13,14,15]. The mean filter is a linear estimator that is employed to mitigate the adverse effects of noise in images by eliminating random fluctuations. The principal disadvantage of the mean filter is that it results in image blurring, particularly at low SNRs. The median filter, which demonstrates superior noise removal performance in comparison to the mean filter, is a nonlinear estimator [16]. The median filter's masking size tends to increase with elevated noise levels, which can result in the loss of information from the original images. Thus, an improved median filtering algorithm (I-median) [17] is developed that combines the advantages of mean and median filters with the objective of achieving superior denoising in adverse environments. Specifically, if the current pixel is less than the average of the mask, the pixel is replaced with the median of the mask. Otherwise, the pixel keeps its original value [17].

In the field of time-domain signal analysis, adaptive filtering (AF) can be used to minimize noise. The AF method primarily uses algorithms such as least mean square (LMS) [18] and recursive least square (RLS) [19] to iteratively adjust the weights, thereby reducing the adverse effects of noise on the signal. The advantage of AF methods is its near real-time implementation [20]. However, a critical aspect of this method lies in the reference signal source, as clean signals, such as speech, are almost nonexistent in certain real-world scenarios. In addition, the choice of step size and filter order, which directly affect the rate of error convergence, is an important issue to address [21]. Finding the right balance is essential for noise reduction when AF is used.

Wiener filtering is an effective noise reduction technique in spectral domain signal analysis [22,23]. It operates by estimating the power spectral density (PSD) of the noise to enhance the signal. However, at low SNRs, the accuracy of the estimated noise PSD may decline, directly impacting noise reduction performance. To address this issue, a Wiener-based noise reduction algorithm known as log-spectral amplitude (LSA) Wiener filtering [24,25,26], has been developed to minimize the mean square error of the logarithmic spectrum. This approach helps alleviate the musical noise and improve the speech quality under low SNR conditions.

Recently, deep learning techniques have also been applied in the field of speech enhancement [27,28,29,30]. One such technique is the end-to-end convolutional recurrent neural network (CRN) [31], which enables near real-time implementation using a single microphone. This CRN algorithm can suppress noise under adverse conditions, such as a -5 dB SNR. However, the main drawback of these learning-based approaches is the substantial amount of data acquisition and preprocessing required [27,28,32], which can be both time- and resource- intensive. Additionally, the noise reduction performance may degrade when encountering unseen scenarios.

Given the compromised estimation of activation patterns and ECAP magnitudes using PECAP under low SNR conditions, this work proposes a method that combines PECAP with a two-stage preprocessing denoising algorithm (TSPD), referred to as PECAP-TSPD. Specifically, the ECAP signals are first denoised by TSPD and then used to estimate the neural parameters by PECAP. In the first stage of TSPD, the I-median algorithm is used to reduce the random noise of an ECAP matrix, since the ECAP matrix can be considered as an image. In the second step of TSPD, the denoised ECAP matrix is expanded and rearranged as a $1\times NN$ vector, where $N$ is the total number of electrodes. The rearranged vector can then be used as the input of the one-dimensional signal to LSA Wiener filtering for residual noise reduction. The PECAP-TSPD algorithm improves the accuracy of the estimates for key parameter, e.g., neural health and current spread. As a result, the refined approach leads to more reliable estimates of neural activation patterns and ECAP magnitudes, improving the effectiveness of neural assessments by clinicians. The normalized root mean square error (RMSE) serves as an objective measure in this work to evaluate the performance of the unprocessed data, PECAP-processed data, and PECAP-TSPD-processed data under various SNR conditions, densities, and scenarios. It also serves to highlight the accuracy of estimating ECAP magnitudes and activation patterns.

2. Panoramic ECAP Method

The panoramic ECAP (PECAP) method consists of two procedures: the forward-masking method and neural parameter estimation using sequential quadratic programming (SQP) [10,11,12]. The forward-masking method uses a probe and a masker to stimulate auditory neurons, generating the ECAP signal and reducing artifact. This process constructs a matrix in which each element represents the measured ECAP signal from each pair of probe and masker positions, as illustrated in Figure 1. In Figure 1, the x-axis represents the position of the masker, and the y-axis represents the position of the probe. Figure 1 shows the simulated ECAP signal, illustrating the inverse relationship between the amplitude of the ECAP and the distance from the probe to the masker. When the positions of the probe and the masker are closer, the amplitude of the ECAP signal is larger. Conversely, when the probe and masker are farther apart, the amplitude of the ECAP signal is smaller. This phenomenon can be explained in Figure 2 [6]. In Figure 2, the neural activity pattern is assumed to follow a Gaussian distribution. The ECAP signal can be regarded as the overlapping area between two Gaussian distributions, corresponding to the neural activity patterns stimulated by the probe and masker, respectively. The overlap is maximized when the positions of the probe and the masker coincide. In contrast, the overlap becomes minimal when the positions of the probe and the masker are distant from each other.

The auditory neural activity pattern is assumed to be a Gaussian distribution as follows:

$$a_i\left(k\right)={\alpha }_i{\eta }_i\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\displaystyle \frac{{(k-{\mathsf{\mu }}_i)}^2}{2{\sigma }_i^2}}\right\}$$

(1)

where $i=\mathrm{1,2},\dots ,N$ denotes the $i$th electrode. $N$ denotes the total number of electrodes. $k$ is the position along the cochlea. Ideally, $k$ is a range from $-\infty$ to $+\infty$. In this work, the range of $k$ is from $1$ to $K$ in which $K$ is set to $N$, the total number of electrode. ${\alpha }_i$ and ${\eta }_i$ are the auditory neuron amplitude and neural health of the $i$th electrode, respectively. $\mathrm{e}\mathrm{x}\mathrm{p}\left\{·\right\}$ is an exponential function. $u_i$ and ${\sigma }_i$ represent the mean and standard deviation (also referred to as current spread¹⁰ in physiological signals) of the Gaussian pattern for the $i$th electrode, respectively. In this work, $u_i=i$ and ${\alpha }_i$ is a prior information. Under an assumption that the ECAP signal is the overlap between the two auditory neuron responses of the probe and masker stimuli, the ECAP signal can be formulated as

$$M\left(p,m\right)=\sum _{k=1}^K{a}_p\left(k\right)a_m\left(k\right),$$

(2)

where $p=\mathrm{1,2},\dots ,N$ is the index indicating that the $p$th electrode is stimulated by the probe. $m=\mathrm{1,2},\dots ,N$ is the index indicating that the $m$th electrode is stimulated by the masker. Equation (2) can also be expressed as a matrix form [12].

$$\mathbf{M}=\sqrt{\mathbf{A}{\mathbf{A}}^T},$$

(3)

in which $\mathbf{A}={\left[\begin{array}{ccc}{\mathbf{a}}_1& {\mathbf{a}}_2& \begin{array}{cc}\cdots & {\mathbf{a}}_N\end{array}\end{array}\right]}^T\in {\mathfrak{R}}^{N\times K}$ is an auditory neuron pattern matrix for each auditory neuron vector ${\mathbf{a}}_i={\left[\begin{array}{ccc}{a}_i\left(1\right)& a_i\left(2\right)& \begin{array}{cc}\cdots & a_i\left(K\right)\end{array}\end{array}\right]}^T\in {\mathfrak{R}}^K$. Because $\mathbf{A}$ is a symmetric matrix, the matrix of the ECAP signal in Equation (3) is also a symmetric matrix. The part of random noise can be alleviated using the following equation:

$${\mathbf{M}}_o={\displaystyle \frac{\mathbf{M}+{\mathbf{M}}^T}{2}},$$

(4)

where ${\mathbf{M}}_o$ represents the denoised matrix of the ECAP signal. To further estimate ${\eta }_i$ and ${\sigma }_i$, the SQP, a constrained optimization algorithm, is utilized. Computerization of the parameter estimation of the ECAP matrix is illustrated in Figure 3.

In Figure 3, ${\mathsf{\epsilon }}_M$, which denotes the normalized root mean square error between the measured and estimated ECAP matrices, is formulated as

$${\mathsf{\epsilon }}_M={\displaystyle \frac{\sqrt{{\frac{1}{N^2}\left|M_s\left(p,m\right)-{\widehat{M}}_s\left(p,m\right)\right|}^2}}{{\overline{M}}_s}},$$

(5)

where ${\overline{M}}_s$ represents the maximum absolute value of ${\mathbf{M}}_s.$ Different SNR and density conditions are also introduced to assess the parameter estimation capability using SQP in this work.

Figure 4 depicts a structure designed to test the noise resistance performance using the SQP algorithm. The ${\mathbf{M}}_{\mathit{s}}$ in Figure 3 and Figure 4 can be replaced by ${\mathbf{M}}_{\mathit{o}}$, as described in Equation (4) for random noise reduction. In addition, the normalized root mean square error between the ground truth $\mathbf{A}$ and the estimated $\widehat{\mathbf{A}}$ can be calculated, as shown in Equation (6).

$${\mathsf{\epsilon }}_A=\sqrt{{{\displaystyle \frac{1}{N^2}}\left|A\left(p,m\right)-\widehat{A}\left(p,m\right)\right|}^2}/\overline{A}$$

(6)

where $\overline{\mathit{A}}$ denotes the maximum absolute value of $\mathbf{A}$. After estimating the parameters of ${\mathit{\eta }}_{\mathit{i}}$ and ${\mathit{\sigma }}_{\mathit{i}}$, the denoised ECAP matrix can be reconstructed using Equation (3). The above procedures can eliminate most part of the noise. However, in low SNR scenarios, the distortion increases due to the error estimation of ${\mathit{\eta }}_{\mathit{i}}$ and ${\mathit{\sigma }}_{\mathit{i}}$. Therefore, the preprocessing algorithm is developed below.

3. Proposed Method

3.1. First Stage Noise Reduction Processing

In light of the detrimental effect of noise on the estimation of neural parameters, this work proposes a two-stage preprocessing denoising (TSPD) algorithm for the ECAP matrix, prior to the PECAP algorithm. First, the ECAP matrix is treated as an image. In the first stage of TSPD, the improved median filtering (I-median) algorithm is applied to reduce noise, as shown in Equation (7).

$$M_I\left(p,m\right)=\left\{\begin{array}{cc}o\left(p,m\right)& M\left(p,m\right)<\mu \left(p,m\right)\\ M\left(p,m\right)& otherwise\end{array}\right.,$$

(7)

where $o(p,m)$ and $\mu (p,m)$ denote the processed values using the median and mean filters at the positions of $(p,m)$, respectively. $M_I(p,m)$ denotes the processed result using I-median. The masker sizes of median and mean filtering are set to $3\times 3$ in this work. Equation (7) describes that if the ECAP value at the positions $(p,m)$ is less than the mean filter processing value, the ECAP value is considered noise and can be replaced by the median filter processing value. Conversely, if the ECAP value is greater, it is retained. The processed ECAP matrix is expressed as

$${\mathbf{M}}_I=\left[\begin{array}{ccc}{M}_I\left(\mathrm{1,1}\right)& \cdots & M_I(1,N)\\ ⋮& \ddots & ⋮\\ M_I(N,1)& \cdots & M_I(N,N)\end{array}\right],$$

(8)

The I-Median algorithm is suitable for removing low to medium density noise from an image. In some cases, however, the noise is distributed across all pixels of an image. To further deal with this type of noise, and assuming that the ECAP signals are the overlaps between the probe and masker stimuli. This work expands ${\mathbf{M}}_{\mathit{I}}$ into a $1\times \mathit{N}\mathit{N}$ vector in accordance with the following rule:

1.

Calculate the absolute value of index $\mathit{p}$ minus index $\mathit{m}.$

$$\mathbf{B}=\left[\begin{array}{ccc}\left|1-1\right|& \cdots & \left|1-N\right|\\ ⋮& \ddots & ⋮\\ \left|N-1\right|& \cdots & \left|N-N\right|\end{array}\right]=\left[\begin{array}{ccc}0& \cdots & \left|1-N\right|\\ ⋮& \ddots & ⋮\\ \left|N-1\right|& \cdots & 0\end{array}\right],$$

(9)

where $\mathbf{B}$ is an index matrix that records the absolute value of position difference between $\mathit{p}$ and $\mathit{m}.$

2.

record $|\mathit{p}-\mathit{m}|$ for each element in $\mathbf{B}$ and concatenate each row into a long vector.

$$\mathbf{b}=\left[b\left(\mathrm{1,1}\right),\dots ,b\left(1,N\right),b\left(\mathrm{2,1}\right),\dots ,b\left(2,N\right),\dots ,b\left(N,1\right),\dots ,b(N,N)\right]\in {\mathfrak{R}}^{1\times NN},$$

(10)

where $\mathit{b}(\mathit{p},\mathit{m})=\left|\mathit{p}-\mathit{m}\right|.$

3.

Conduct a descending order of $\mathbf{b}$ and record the descending order index.

$$\left({\mathbf{i}}_d,\overline{\mathbf{b}}\right)=descend\left(\mathbf{b}\right),$$

(11)

where $\mathit{d}\mathit{e}\mathit{s}\mathit{c}\mathit{e}\mathit{n}\mathit{d}\left(·\right)$ represents the descending order operation. $\overline{\mathbf{b}}$ is vector based on the descending order results of $\mathbf{b}$. ${\mathbf{i}}_{\mathit{d}}$ is the index vector corresponding to $\overline{\mathbf{b}}$.

4.

The desired vector then can be obtained using the following equation:

$${\mathit{m}}_I=\left[M_I\left({\mathit{i}}_d\left(1\right)\right),M_I\left({\mathit{i}}_d\left(2\right)\right),\dots ,M_I\left({\mathit{i}}_d\left(NN\right)\right)\right]\in {\mathfrak{R}}^{1\times NN},$$

(12)

By examining the ECAP signal in Figure 1 and Figure 2, the maximum can be reached when the probe and masker positions are close, and the minimum can be reached when the probe and masker positions are far apart. Therefore, in the previous steps, the noise (where ECAP signal is weak) would be placed in the first part of the vector. The target signal (where strong ECAP signal) can then be placed in the last part of the vector. The reordering processing is part of the second stage of TSPD, which employs log-spectral amplitude (LSA) Wiener filtering for improved noise reduction.

3.2. Second Stage Noise Reduction Processing

In the second stage of TSPD, LSA Wiener filtering is utilized to address the residual noise in the vector ${\mathbf{m}}_I$, which is treated as a discrete-time series input. The first step is to transform ${\mathbf{m}}_I$ into a short-time Fourier transform (STFT) domain, as shown in Equation (13).

$$Y_m\left(l,k\right)=\sum _{n=0}^{N_{FFT}-1}{M}_I(n+l{N}_{hop})w\left(n\right)W_{N_{FFT}}^{kn},$$

(13)

where $Y_m\left(l,k\right)$ is the STFT signal of ${\mathbf{m}}_I$. $l$ and $k$ denote the time frame and frequency indices, respectively. $N_{FFT}$ is the size of fast Fourier transform (FFT), $N_{hop}$ is the hop size, and $w\left(n\right)$ is the window for short-time signal analysis. $W_{N_{FFT}}^k=\mathrm{e}\mathrm{x}\mathrm{p}\left\{-j2\pi k/N_{FFT}\right\}$ is the $N_{FFT}$th root of unity with $k=\mathrm{0,1},\dots ,N_{FFT}-1$. The LSA Wiener filter aims at minimizing the log-spectral amplitude

$$J=\mathrm{a}\mathrm{r}\mathrm{g}\underset{H(l,k)}{\min }E\left\{{\left|{\log }_e\left|X\left(l,k\right)\right|-{\log }_e\left|\widehat{X}\left(l,k\right)\right|\right|}^2\right\},$$

(14)

where $J$ is the cost function to minimize the mean square error of log-spectral amplitudes ${\log }_e\left|X\left(l,k\right)\right|$ and ${\log }_e\left|\widehat{X}\left(l,k\right)\right|$. $X\left(l,k\right)$ and $X\left(l,k\right)$ are the clean signal and estimated clean signal spectrums, respectively. The optimal solution is

$$H_{LSA}\left(l,k\right)={\displaystyle \frac{\xi (l,k)}{1+\xi (l,k)}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{{\displaystyle \frac{1}{2}}{\int }_{\nu (l,k)}^{\infty }{\displaystyle \frac{e^{-t}}{t}}dt\right\},$$

(15)

where $\xi \left(l,k\right)=E\left\{{\left|X(l,k)\right|}^2\right\}/P_u\left(l,k\right)$ is the prior SNR with $P_x\left(l,k\right)=E\left\{{\left|X(l,k)\right|}^2\right\}$ being a clean signal PSD and $P_u\left(l,k\right)=E\left\{{\left|U_r(l,k)\right|}^2\right\}$ being a residual noise PSD. $\upsilon \left(l,k\right)=\left(\xi \left(l,k\right)/1+\xi \left(l,k\right)\right)\gamma (l,k)$ with $\gamma \left(l,k\right)={\left|Y_m(l,k)\right|}^2/P_u\left(l,k\right)$ representing the posterior SNR. $\xi \left(l,k\right)$ can be estimated using the decision-directed approach as described below:

$$\widehat{\xi }\left(l,k\right)=\alpha {\displaystyle \frac{{\left|X(l-1,k)\right|}^2}{P_u\left(l-1,k\right)}}+(1-\alpha )\mathrm{m}\mathrm{a}\mathrm{x}\left\{\gamma \left(l,k\right)-1,{\epsilon }_0\right\},$$

(16)

where $\alpha$ is a forgetting factor. $P_u\left(l,k\right)$ can be estimated and updated using the following log-likelihood ratio criterion:

$$\mathsf{\Lambda }\left(l,k\right)=\mathrm{l}\mathrm{n}{\displaystyle \frac{f\left(Y_m\right(l,k\left)\right|H_1)}{f\left(Y_m\right(l,k\left)\right|H_0)}}=-\ln \left(1+\xi \left(l,k\right)\right)+\upsilon (l,k),$$

(17)

where

$$f\left(Y_m\left(l,k\right)|H_1\right)={\displaystyle \frac{1}{\pi \left({P_x\left(l,k\right)+P}_u\left(l,k\right)\right)}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\displaystyle \frac{{\left|Y_m(l,k)\right|}^2}{P_x\left(l,k\right)+P_u\left(l,k\right)}}\right\}$$

(18)

is the conditional probability density function (PDF) of $Y_m\left(l,k\right)$ given that the event ($H_1$) of the ECAP signal occurs.

$$f\left(Y_m\left(l,k\right)|H_0\right)={\displaystyle \frac{1}{\pi P_u\left(l,k\right)}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\displaystyle \frac{{\left|Y_m(l,k)\right|}^2}{P_u\left(l,k\right)}}\right\}$$

(19)

is the conditional PDF of $Y_m\left(l,k\right)$ assuming only noise occurs at event $H_0$. $\upsilon \left(l,k\right)=\gamma \left(l,k\right)\xi \left(l,k\right)/(1+\xi \left(l,k\right))$. If $\sum _{l=1}^{N_t}\mathsf{\Lambda }\left(l,k\right)$ is less than a small value ${\epsilon }_1$, the noise PSD can be updated using the following recursive averaging.

$${\widehat{P}}_u\left(l,k\right)=\alpha {\widehat{P}}_u\left(l-1,k\right)+\left(1-\alpha \right){\left|Y_m\left(l,k\right)\right|}^2.$$

(20)

The estimated clean ECAP signal can be obtained using the following equation:

$$\widehat{X}\left(l,k\right)=H_{LSA}\left(l,k\right)Y_m(l,k).$$

(21)

The complete LAS Wiener filtering procedures used in the second noise reduction stage are listed in Table 1.

In Table 1, $N_u$ is the number of time frames used to represent the presence of noise. Next, the processed row vector $\widehat{X}\left(l,k\right)$ is conversed into using the inverse STFT the time domain signal $\widehat{x}\left(n\right)$ which can be reconstructed into the matrix format, as described below:

$${\tilde{\mathbf{M}}}_I=\left[\begin{array}{ccc}\widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=1)& \cdots & \widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=N)\\ \begin{array}{c}\widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=N+1)\\ ⋮\end{array}& \ddots & \begin{array}{c}\widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=2N)\\ ⋮\end{array}\\ \widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=\left(N-1\right)N+1)& \cdots & \widehat{x}\left(n^{\text{'}}\right|{\mathbf{i}}_d\left(n^{\text{'}}\right)=NN)\end{array}\right].$$

(22)

where ${\mathbf{i}}_d\left(n^{\text{'}}\right)$ is defined in Equation (11) with $n^{\text{'}}$ being a range of $1$ to $NN.$ In this work, the mean filter is applied in ${\tilde{\mathbf{M}}}_I$ if the SNR value is below 4 dB, in order to leverage the advantage of random noise removal at low SNRs [15].

4. Settings and Results

Two types of noise – random noise and impulse noise – are used to evaluate the performance of the proposed PECAP-TSPD algorithm. Twelve SNR levels, -5 dB, -2 dB, 1 dB, 4 dB, 7 dB, 10 dB, 13 dB, 16 dB, 19 dB, 22 dB, 25 dB, and 100 dB (representing the clean ECAP signal situation for random noise case) are used in the random noise. Four densities, 10%, 20%, 30%, and 40%, are used in the impulse noise. Normalized RMSE is used as an objective quality measure to calculate the error between the ground truth and the estimated results (neural health and ECAP amplitude). Seven different combinations of neural health and current spread are listed in Table 2. The results of the ECAP matrices before and after processing at -5 dB of SNR are shown in Figure 5.

In Figure 5(b), the ECAP matrix is filled with random noise, thereby, making it difficult to observe the pristine measured ECAP data, as depicted in Figure 5(a). The PECAP and PECAP-TSPD methods can mitigate the detrimental effects of extremely noisy environments to recover the clean ECAP matrix shown in Figure 5(c) and (d). The images of Figure. 5(a) and (d) are almost identical, further presenting the satisfactory performance of PECAP-TSPD. The results of the normalized RMSE of the ECAP magnitude (${\mathsf{\epsilon }}_M$) and neural activity pattern (${\mathsf{\epsilon }}_A$) are depicted in Figure 6.

Figure 6(a) shows the normalized RMSE of the unprocessed ECAP signals, ECAP signals processed by PECAP and PECAP-TSPD algorithms for different SNRs in Scenario 1. The normalized RMSE of the magnitude of the unprocessed ECAP signals at -5 dB SNR increases to 83.39%, which is comparably higher than those of the processed ECAP signals by PECAP (16.17%) and PECAP-TSPD (5.23%), indicating the need for ECAP signal processing. When comparing the ${\mathsf{\epsilon }}_A$ between PECAP and PECAP-TSPD processing ECAP signals, the values of ${\mathsf{\epsilon }}_A$ by PECAP-TSPD are all smaller than those by PECAP except the 100 dB SNR case, where ${\mathsf{\epsilon }}_A$ are $0.0043\%$ and $0.1885\%$ for PECAP and PECAP-TSPD, respectively. Figure 6(b) presents the curve of the average normalized RMSE, which is formulated as follows for each scenario:

$${\overline{\mathsf{\epsilon }}}_M=\sum _{k=1}^{N_{SNR}}{\displaystyle \frac{\sqrt{{\frac{1}{N^2}\left|M_{s,k}\left(p,m\right)-{\widehat{M}}_{s,k}\left(p,m\right)\right|}^2}}{{\overline{M}}_{s,k}}},$$

(23)

where $N_{SNR}=12$ represents the total number of SNR conditions in this work. $M_{s,k}\left(p,m\right)$ and ${\widehat{M}}_{s,k}\left(p,m\right)$ represent each index ($p,m$) value of each SNR condition ($k$) in the pristine and estimated ECAP matrices, respectively. ${\overline{M}}_{s,k}$ is the maximum absolute value of ${\mathbf{M}}_{s,k}$ which denotes the pristine ECAP matrix of each SNR condition. Similar to ${\overline{\mathsf{\epsilon }}}_M$, ${\overline{\mathsf{\epsilon }}}_A$ is described as

$${\overline{\mathsf{\epsilon }}}_A=\sum _{k=1}^{N_{SNR}}\sqrt{{{\displaystyle \frac{1}{N^2}}\left|A_{s,k}\left(p,m\right)-{\widehat{A}}_{s,k}\left(p,m\right)\right|}^2}/{\overline{A}}_{s,k},$$

(24)

in which $A_{s,k}\left(p,m\right)$ and ${\widehat{A}}_{s,k}\left(p,m\right)$ denote each index value of each SNR condition in the pristine and estimated neural activity patterns, respectively. ${\overline{A}}_{s,k}$ is the maximum absolute value of ${\mathbf{A}}_{s,k}$ which denotes the pristine neural activity pattern of each SNR condition. Compared to the unprocessed and processed ECAP signals, the maximum values of ${\overline{\mathsf{\epsilon }}}_M$ are 6.17% and 5.48% for PECAP and PECAP-TSPD, respectively. In contrast, the maximum value of the unprocessed ECAP signals is 28.97% for ${\overline{\mathsf{\epsilon }}}_M$, showing that the noise resistance capability of the PECAP and PECAP-TSPD algorithms. The performance of the PECAP-TSPD algorithm is superior to that of the PECAP algorithm, as the values of ${\overline{\mathsf{\epsilon }}}_M$ and ${\overline{\mathsf{\epsilon }}}_A$ obtained with PECAP-TSPD are lower than those obtained with PECAP. The average of ${\overline{\mathsf{\epsilon }}}_M$ from Scenarios 1 to 7 can be ranked as: PECAP-TSPD (3.83%), PECAP (5.14%), and unprocessed (23.07%). The averages of ${\overline{\mathsf{\epsilon }}}_A$ from Scenarios 1 to 7 for PECAP-TSPD and PECAP are 3.01% and 4.64%, respectively.

Next, the impulse noise is added into the EACP matrix to evaluate the performance of PECAP-TSPD under four different densities. Figure 7 illustrates the results of the ECAP matrices before and after using the PECAP and PECAP-TSPD algorithms at the 40% density of the impulse noise.

The impulse noise with 40% density heavily contaminates the original ECAP matrix (Figure 7(a)), as shown in Figure 7(b), emphasizing the importance of signal processing. The processed ECAP matrices using the PECAP and PECAP-TSPD algorithms are depicted in Figure 7(c) and (d), where the impulse noise is most reduced. When comparing Figure 7(c) to Figure 7(d), the restored ECAP matrix in Figure 7(d) shows more resemblance to Figure 7(a) than to that of Figure 7(c), exhibiting the satisfactory performance of PECAP-TSPD under adverse noisy environments. The normalized RMSE results of the ECAP magnitude and neural activity pattern of the impulse noise case are depicted in Figure 8.

Figure 8(a) presents the ${\mathsf{\epsilon }}_M$ and ${\mathsf{\epsilon }}_A$ curves at four distinct densities in Scenario 2 for the unprocessed, PECAP, and PECAP-TSPD approaches, all of which increase as the impulse noise density increases. In the case of 40% impulse noise density, the ${\mathsf{\epsilon }}_M$ values for unprocessed, PECAP, and PECAP-TSPD are 34.07%, 8.44%, and 2.96%, respectively, displaying the effectiveness of PECAP and PECAP-TSPD in highly noisy environments. The average normalized RMSE results are depicted in Figure 6(b), in which ${\overline{\mathsf{\epsilon }}}_M$ and ${\overline{\mathsf{\epsilon }}}_A$ are similar to Eqs. (23) and (24), except that $N_{SNR}$ is replaced by $N_{density}=4$, denoting the total number of impulse noise densities in this work. The maximum value of ${\overline{\mathsf{\epsilon }}}_M$ for the PECAP and PECAP-TSPD algorithms is 6.77% and 4.22%, respectively. For the unprocessed ECAP matrices, the maximum value of ${\overline{\mathsf{\epsilon }}}_M$ is 27.59%, which suggests that PECAP and PECAP-TSPD are robust against noise. The PECAP-TSPD algorithm performs better than the PECAP algorithm because the values of ${\overline{\mathsf{\epsilon }}}_M$ and ${\overline{\mathsf{\epsilon }}}_A$ calculated by PECAP-TSPD are lower than those calculated by PECAP. The mean values of ${\overline{\mathsf{\epsilon }}}_M$ from Scenarios 1 through 7 can be arranged in ascending order as follows: PECAP-TSPD (3.13%), PECAP (5.57%), and unprocessed (21.97%). The mean values of ${\overline{\mathsf{\epsilon }}}_A$ from Scenarios 1 to 7 for PECAP-TSPD and PECAP are 2.39% and 5.23%, respectively.

5. Conclusions

The PECAP-TSPD algorithm, which integrates the improved median filter, the log-spectral Wiener filter, and PECAP, has been developed for noise reduction in ECAP data to accurately estimate neural health and current spread from severely noisy ECAP matrices. The normalized RMSE for ECAP magnitude (${\mathsf{\epsilon }}_M$) and neural activity pattern (${\mathsf{\epsilon }}_A$) demonstrates that both the PECAP and PECAP-TSPD algorithms effectively reduce ${\mathsf{\epsilon }}_M$ and ${\mathsf{\epsilon }}_A$ compared to unprocessed data under various SNRs, noise densities, and experimental scenarios. Notably, PECAP-TSPD outperforms PECAP in terms of both ${\mathsf{\epsilon }}_M$ and ${\mathsf{\epsilon }}_A$. For ECAP matrices corrupted by random noise, the average ${\mathsf{\epsilon }}_M$ across seven different scenarios and twelve SNR levels is as follows: PECAP-TSPD (3.83%), PECAP (5.14%), and unprocessed (23.07%). In the case of impulse noise, the ranking is: PECAP-TSPD (3.13%), PECAP (5.57%), and unprocessed (21.97%). Similarly, the average ${\mathsf{\epsilon }}_A$ for PECAP-TSPD and PECAP under random noise conditions are 3.01% and 4.64%, respectively. Under impulse noise, the values are 2.39% for PECAP-TSPD and 5.23% for PECAP.