Enhanced Improved Proportionate Set-Membership Fast NLMS Using the l_0 Norm

The aim of this research work is to address the challenging problem of AEC for different sparsity levels of acoustic impulse responses (AIRs) [1]. Sparse system identification is a prevalent issue in various practical applications [1-3]. For instance, in AEC the AIR often exhibits sparsity, where only a small percentage of components have significant magnitudes while the rest are close to zero. The commonly used algorithms, NLMS and Recursive Least Squares (RLS) do not consider the sparse nature of such systems. This led to the development of "proportionate" approaches, resulting in the PNLMS algorithms [2-4]. The idea behind these algorithms is to individually update each coefficient by adjusting the adaptation gain proportionally to the estimated coefficient's magnitude. This prioritizes coefficients with high magnitudes, leading to faster initial convergence and an overall improved convergence rate. Unfortunately, the enhanced performance of the PNLMS algorithms is not sustained throughout the adaptation process. This is because the PNLMS algorithm accelerates the convergence of high-magnitude coefficients at the expense of slowing down the convergence of low-magnitude coefficients. Furthermore, when the Impulse Response (IR) is dispersive, PNLMS performs worse than the NLMS [3-5]. Multiple variations of the PNLMS algorithm have been introduced to enhance the resilience of proportionate algorithms across diverse system types. Among these alternatives, the IPNLMS algorithm [3-5] has emerged as a successful option. It utilizes the $\ell_1$ norm to leverage the sparsity of the IR along with the combination of the NLMS gain and a proportionality rule in the PNLMS. The IPNLMS-L0 algorithm, proposed in studies [6-8], replaces the $\ell_1$ norm with the $\ell_0$ norm in the IPNLMS algorithm. This substitution is based on the recognition of the $\ell_0$ norm as an effective and intuitive mathematical measure of sparseness. The IPNLMS-L0 algorithm introduces a new parameter β for calculating the $\ell_0$ norm and incorporates an exponential term in the gain control matrix equation. It effectively utilizes the sparsity of the IR and performs well, particularly when dealing with sparse IR. However, it has practical limitations [4, 7]. Firstly, choosing an appropriate value for the parameter β is crucial and depends on the sparseness of the IR. Additionally, the exponential term in the gain control matrix equation may pose challenges during real-world implementations.

In addition, recent studies have revealed that algorithms with proportionate characteristics, such as the IPNLMS algorithm, exhibit inadequate performance in numerous sparse systems or when confronted with non-stationary and compressible systems. Specifically, these systems refer to those that do not possess sufficient sparsity [9].

On a different note, while the RLS algorithm and its variants demonstrate a considerably faster convergence rate than the NLMS algorithm, they suffer from high computational complexity. To address this issue, researchers have conducted several experiments to achieve the fast convergence of the RLS algorithm while reducing its computational complexity. The result of these efforts is the development of the FNLMS algorithm [10-13]. When compared to the NLMS algorithm, the FNLMS algorithm exhibits a faster convergence rate and enhanced tracking capability.

The FNLMS algorithm has recently integrated the Set-Membership (SM) approach [14-16]. This integration involves updating the step size based on the estimated output error, resulting in the development of an ISM-FNLMS algorithm.

These enhancements lead to reduced computational update costs and improved tracking capability. However, the algorithm still requires additional consideration for sparseness and an upgrade to incorporate it. In a more recent development, a hybrid sparseness-aware algorithm called Improved Proportionate FNLMS (IPFNLMS) has been created [17]. This algorithm incorporates the gain control matrix of the IPNLMS algorithm into the Kalman-based adaptation gain of the FNLMS algorithm. In comparison to the IPNLMS algorithm, the IPFNLMS algorithm exhibits similar computational costs but demonstrates superior behavior.

The practical limitations of the IPNLMS-L0 algorithm are well-known [4, 7]: the determination of the parameter β and the exponential term evaluation. To address these limitations, the proposed EIPSMFNLMS-L0 algorithm provides a solution by introducing a method that incorporates sparseness measure based on the $\ell_1, \ell_2$, and $\ell_{\infty}$ norms which accurately captures both sparse and dispersive IRs to determine the parameter β and employing an efficient approximation of the exponential term through a Maclaurin expansion. This approach yields to an effective gain control matrix that not only mitigates the shortcomings of the IPNLMS-L0 algorithm but also works effectively and robustly to the variations in the sparseness of the IR. To further reduce the computational complexity of our proposed algorithm and enhance the algorithm's ability to track time-varying changes in the system. we take a step further by integrating our new gain control matrix obtained into the concept of ISMFNLMS. This integration involves incorporating the gain control matrix into the Kalman-based adaptation gain of the FNLMS algorithm and subsequently integrating the SM concept. The resulting EIPSMFNLMS-L0 algorithm achieves a desirable balance between convergence speed, steady-state performance, lower MSE level and robustness to sparse and dispersive AIRs with low-cost complexity, surpassing NLMS, IPNLMS and IPNLMS-L0 and significantly surpasses the recent ISMFNLMS and IPFNLMS algorithms. This study offers comprehensive derivations and explanations of the proposed algorithm, along with the results of numerical simulations. The focus is on comparing the performance of the proposed algorithm with existing algorithms in terms of convergence speed, tracking capability, final mean square error (MSE), and computational complexity. The evaluation covers scenarios involving both sparse and dispersive IRs, as well as stationary and non-stationary signals and systems.

3.1 The standard FNLMS algorithm

The FNLMS algorithm, a novel algorithm obtained by simplifying the Fast RLS algorithm, is introduced. This algorithm demonstrates rapid convergence and offers low computational complexity. The adaptation gain for the FNLMS algorithm is given in studies [12, 13] as:

$\boldsymbol{g}(n)=\gamma(n) \tilde{\boldsymbol{c}}(n)$               (14)

where, the likelihood variable is $\gamma(n)$, and the dual Kalman gain is $\tilde{\boldsymbol{c}}(n)$. By completely disregarding the backward and forward predictors and solely utilizing the forward prediction error of the input signal, the evaluation of the dual Kalman gain is conducted. This evaluation process is described in studies [11-13] as follows:

$\left[\begin{array}{c}\tilde{\boldsymbol{c}}(n) \\ c(n)\end{array}\right]=\left[\begin{array}{c}-\frac{\varepsilon(n)}{\lambda \alpha(n-1)+C_0} \\ \tilde{\boldsymbol{c}}(n-1)\end{array}\right]$            (15)

A small positive constant, $C_0$, is employed to prevent division by extremely small values or zero, especially in cases where the input signal remains absent for a prolonged period. The forgetting factor, λ, falls within the range of 0 to 1, while a(n) represents the variance of the forward prediction error.

$\alpha(n)=\lambda \alpha(n-1)+\varepsilon^2(n)$            (16)

The prediction error $\varepsilon(n)$ in Eq. (15) and Eq. (16) is computed by a first-order prediction model.

$\varepsilon(n)=x(n)-a(n) x(n-1)$                (17)

In scenarios where the input statistics are often unknown or subject to variations over time, the prediction parameter a(n) needs to be determined based on the input signal. To estimate a(n), the following equation is employed:

$a(n)=\frac{r_1(n)}{r_0(n)+C_a}$             (18)

To prevent division by values that are very close to zero, $C_a$ is introduced as a small positive constant. The estimate of the first lag correlation function of $x(n)$ is denoted as $r_1(n)$, while $r_0(n)$ represents the estimate of the power of the input signal. Recursive formulas are used to estimate these two functions, and they are given by:

$r_1(n)=\lambda_a r_1(n-1)+x(n) x(n-1)$             (19)

$r_0(n)=\lambda_a r_0(n-1)+x^2(n)$               (20)

where, $\lambda_a$ is an exponential forgetting factor. The likelihood variable $\gamma(n)$ is recursively defined as [13].

$\gamma(n)=\frac{\gamma(n-1)}{1+\gamma(n-1) o(n)}$             (21)

where,

$o(n)=c(n) x(n-L)+\frac{x(n) \varepsilon(n)}{\lambda \cdot \alpha(n-1)+C_0}$               (22)

In order to have better control over the adaptation gain of the filtering component, a fixed step-size μ has been incorporated. The update equation for this algorithm is defined as:

$\widehat{\boldsymbol{h}}(n)=\widehat{\boldsymbol{h}}(n-1)-\mu \boldsymbol{g}(n) e(n)$               (23)

3.2 The ISMFNLMS algorithm

The SM adaptive filtering presents a new algorithm known as SM-NLMS. Its update equation is similar to that of the NLMS algorithm, but with an "optimal" adaptive step size. This means that the step size assumes a value of zero when the optimization problem described below has been verified, indicating that there is no need for further updates [14-16]. This SM concept was expanded to include the FNLMS algorithm. This extension involved replacing the a priori error with its estimates, which resulted in the development of the ISM-FNLMS algorithm, as described in the study [16]. The objective of this algorithm is to minimize the optimization problem represented by $\|\widehat{\boldsymbol{h}}(n)-\widehat{\boldsymbol{h}}(n-1)\|^2$, while adhering to the constraint that the estimated impulse response vector $\widehat{\boldsymbol{h}}(n)$ belongs to the set $\mathcal{H}(n)$ which is a set comprising all vectors $\widehat{\boldsymbol{h}}(n)$ for which the error is within by a certain upper bound $\omega$. In simpler terms, $\boldsymbol{\mathcal { H }}(n)$ encompasses the vectors whose estimation error is limited by $\omega$. This set can be defined as [14-16]:

$\mathcal{H}(n)=\left\{\widehat{\boldsymbol{h}} \in \mathbb{R}^{L+1}:|d(n)-\widehat{\boldsymbol{h}}(n-1) \boldsymbol{x}(n)| \leq \omega\right\}$                (24)

The achievement of the constraint involves orthogonal projection of the previous estimate $\widehat{\boldsymbol{h}}(n-1)$ onto the closest boundary of the set $\boldsymbol{\mathcal { H }}(\mathrm{n})$. As a result, the ISM-FNLMS algorithm has an updated equation given as follows [14]:

$\widehat{\boldsymbol{h}}(n)=\widehat{\boldsymbol{h}}(n-1)-\mu(n) \boldsymbol{g}(n) e(n)$           (25)

The step-size $\mu(n)$ in Eq. (25) is defined as follows:

$\mu(n)=\left\{\begin{array}{cc}1-\frac{\omega}{\sigma_e(n)} & \text { if }|e(n)|>\omega \\ 0 & \text { otherwise }\end{array}\right.$           (26)

Here, $\omega$ represents the given bound, and $\sigma_e(n)$ is an estimation of the modulus of the error signal. Initially, $\sigma_e(n)$ is set to a value close to the input signal power $\sigma_x$. The parameter $\zeta$ is a forgetting factor, and Eq. (27) defines the update rule for $\sigma_e(n)$ [16]:

$\sigma_e(n)=\zeta \sigma_e(n-1)+(1-\zeta)|e(n)|$            (27)

Our study involved a comprehensive performance evaluation of the proposed EIPFSMNLMS-L0 algorithm in comparison to several other algorithms, including ISMFNLMS, IPNLMS, IPNLMS-L0, IPFNLMS, and classical NLMS algorithms.

5.1 Experimental setup

The schematic diagram in Figure 1 was employed to investigate an AEC application. We conducted simulations under various operating conditions. These conditions encompassed the use of two distinct correlated stationary input signals, namely USA Standards Institute (USASI) depicted in Figure 2 (a) which is a zero-mean correlated noise, featuring a spectral dynamic range of 29dB, and auto-regressive (AR20) which is a simulated by AR stationary process of order 20 with a spectral dynamic range of 42dB shown in Figure 2 (b), as well as a speech input signal, displayed in Figure 2 (c), is a real non-stationary signal characterized by a spectral dynamic range of 46dB. All signals used in this simulation are sampled at a rate of 16 kHz.

Additionally, we considered different impulse responses with different levels of sparseness which were synthetically generated using the method described in previous studies [5, 17] to generate output signal given in Figures 3-5. These synthetic IRs are subsequently named SIRs.

5.1.1 Time-varying SIRs generation

The time-varying scenarios were generated by multiplying the SIRs by an artificial gain depicted in Figure 6. This gain variation was applied between 60000 and 80000 samples for USASI noise and in the interval 90000 to 110000 samples for the AR20 input noise. It is important to note that the artificial gain was only applied to the SIRs corresponding to stationary input signals.

5.1.2 Generation of desired output signals

The time-varying SIRs were convolved with the three input signals. We used three sets of two pairs and one triad of SIRs, each set of equal length consisting of sparse and dispersive SIRs. The first pair, shown in Figures 3 (a) and 3 (b), was used in conjunction with the USASI input signal and length L=256. The second pair, shown in Figures 4 (c) and 4 (d), with a length of L=512, was used with the AR (20) input signal. Finally, the third triplet, shown in Figures 5 (e), 5 (f) and 5 (g), with a length of L=1024, was used with the speech input signal. To include the noise component to the desired signal $d(n)$, we added a white Gaussian noise $b(n)$ with an SNR of 50 dB or 30 dB to the echo signal $y(n)$.

2.png

Figure 2. The used input signals

3.png

Figure 3. Two types of SIRs with L=256: (a) $\hat{\xi}_{12}(n)=0.8296$, (b) $\hat{\xi}_{12}(n)=0.3028$

5.2 Performance metric

The Mean-Square Error value (MSE) is the performance metric employed in all our simulations. The MSE is expressed in decibels (dB) and is calculated as follows:

$M S E(d B)=10 \log _{10}\left(<e^2(n)>\right)$             (46)

where, < .> indicates time averages over blocks of 1024 samples for USASI noise and SPEECH signals, and 2048 samples for AR (20).

4.png

Figure 4. Two types of SIRs with L=512: (c) $\hat{\xi}_{12}(n)=0.8637,(d) \hat{\xi}_{12}(n)=0.4080$

5.png

Figure 5. Three types of SIRs with L=1024: (e) $\hat{\xi}_{12}(n)=0.8789$, (f) $\hat{\xi}_{12}(n)=0.6506,(\mathrm{~g}) \hat{\xi}_{12}(n)=0.4552$

6.png

Figure 6. Artificial gain for tracking capability tests

5.3 Selection of simulation parameters

Regarding the algorithm parameters, we selected a common step-size of $\mu=0.5$ for all algorithms to achieve a roughly consistent steady-state MSE value. The remaining parameters were chosen to provide good performance with a trade-off between convergence, tracking ability, and the minimum steady-state MSE value.

For the ISMFNLMS, IPFNLMS, and EIPSMFNLMS-L0 algorithms, two forgetting factors were employed. The first one $\lambda_a$ is used to track changes in the statistics of the input signal [9-13] and it is evaluated over a rectangular window of 35 ms (equivalent to a 560 -sample window) $\lambda_a \approx \zeta \approx 1-$ $1 / 560=0.9982$. The second forgetting factor $\lambda$ captures variations in the unknown system and is determined by $\lambda=$ $1-1 / L$. We considered three different lengths of impulse response (L=256, 512 and 1024).

The regularization constants $C_0$ and $C_a$ were initially set to values comparable to the input signal power and were subsequently modified to achieve equivalent steady-state MSE levels for ISMFNLMS, IPFNLMS, and EIPSMFNLMS-L0 algorithms. The constant $E_0$ is utilized during initialization and set to $E_0=1$. The initial value of $\sigma_e(0)$, is set to $\sigma_x / 100$ at an SNR of 50 dB or $\sigma_x / 10$ when the SNR is 30 dB.

The calculation of regularization parameters differs among the algorithms. The IPFNLMS and EIPSMFNLMS-L0 algorithms utilize Eq. (40), whereas the NLMS algorithm uses $\delta_{N L M S} \approx \sigma_x^2$. In the case of the IPNLMS-L0 algorithm, the value of $\beta$ is adjusted based on the type of input signal 0.5 for stationary inputs, 50 for sparse IRs and 500 for dispersive IRs with a speech input signal. For IPNLMS-L0, a specific value of $\rho=0.5$ is assigned. On the other hand, for IPNLMS, IPFNLMS, and EIPNLMS-L0, $\rho$ is consistently set to -0.5.

Table 3 summarizes the length of the SIRs utilized for each input signal and the associated forgetting factors for each length.

Table 3. Signals, synthetic impulse response lengths, and forgetting factors

5.4 Results and discussions

We began the performance comparisons using the USASI input with the first pair of time-varying SIRs (for both sparse and dispersive cases) of length L=256, as illustrated in Figure 7 for sparse case and Figure 8 for dispersive case. The results reveal that IPNLMS-L0 outperforms IPNLMS in handling sparse IRs but suffers from slower convergence and reduced tracking as sparseness decreases. While ISMFNLMS shows improvement for dispersive IRs, all algorithms fall behind our proposed algorithm EIPSMFNLMS-L0, which excels in faster convergence, superior tracking and lower error for sparse and dispersive cases with a lower computational complexity, surpassing even the recently introduced IPFNLMS algorithm.

In the case of the AR (20) input signal convoluted with the time-varying second pair of SIRs of L=512, the proposed algorithm demonstrates superior overall performance in situations with sparse IRs, as shown in Figure 9. It achieves the fastest convergence speed and the strongest tracking capability compared to the NLMS, IPNLMS, IPNLMS-L0, ISMFNLMS and IPFNLMS algorithms. Concerning dispersive SIRs, shown in Figure 10, the EIPSMFNLMS-L0 algorithm exhibits a fastest convergence speed approximately equivalent to that of the ISMFNLMS algorithm. However, it surpasses the other algorithms in terms of convergence speed, tracking ability, and final MSE.

7.png

Figure 7. Performance comparison under time-varying- in sparse $\hat{\xi}_{12}(n)=0.8296$ SIRs of length $L=256$. Using USASI noise Input signal. $C_0=C_a=1.2, \omega=0.0025$ and SNR=50dB

8.png

Figure 8. Performance comparison under time varying-time in dispersive $\hat{\xi}_{12}(n)=0.3028$ SIRs of length $L=256$. Using USASI noise Input signal. $C_0=C_a=1.2, \omega=0.0025$ and SNR=50dB

9.png

Figure 9. Performance comparison under time varying- in sparse $\hat{\xi}_{12}(n)=0.8637$ SIRs of length $L=512$. Using AR20 noise input signal. $C_0=C_a=1.8, \omega=0.0025$ and SNR=50dB

10.png

Figure 10. Performance comparison under time varying- in dispersive$\hat{\xi}_{12}(n)=0.4080$ SIRs of length $L=512$. Using AR (20) noise Input signal. $C_0=C_a=1.8, \omega=0.0025$ and SNR=50dB

11.png

Figure 11. Performance comparison in sparse$\hat{\xi}_{12}(n)=0.6539$ SIRs of length $L=1024$. Using SPEECH Input signal. $C_0=C_a=1.2, \omega=0.0025$ and SNR=50dB

12.png

Figure 12. Performance comparison in dispersive $\hat{\xi}_{12}(n)=0.4552$ synthetic AIRs of length $L=1024$. Using SPEECH Input signal. $C_0=C_a=1.2, \omega=0.0025$ and SNR=50dB

The comparison with speech input signals and the third pair of SIRs of length L = 1024, with signal-to-noise ratio SNR = 50 dB. Figures 11 and 12 show the results obtained for these cases, as it can be seen that these tests demonstrate the superior performance of the proposed algorithm in both sparse and dispersive systems, reinforcing its effectiveness compared to other algorithms.

Figure 13 shows an additional experiment in a rather difficult situation with non-stationary speech input convolved with a large SIR with L=1024 points and very sparse $\hat{\xi}_{12}(n)=0.8789$ with an SNR of 30 dB. This last experiment confirms the superiority of the proposed algorithm compared to other algorithms.

Based on these simulations, it is clear that our proposed algorithm outperforms all algorithms used in this study in terms of MSE, convergence speed and tracking capability, at a reduced computational cost.

13.png

Figure 13. Performance comparison in very sparse $\hat{\xi}_{12}(n)=0.8789$ SIRs of length $L=1024$. Using SPEECH Input signal. $C_0=C_a=1.2, \omega=0.015$ and SNR=30dB