The Ultra High Speed LMS Algorithm Implemented on Parallel Architecture Suitable for Multidimensional Adaptive Filtering

1. Introduction

Over the past decades a number of new adaptive filter algorithms have been elaborated and applied to meet demands for faster convergence and better tracking properties than earlier techniques could offer. The Filtered LMS algorithm is currently the most popular method for adapting a filter, due to its simplicity and robustness, which have made it widely adopted in many applications. Applications include adaptive channel equalization, adaptive predictive speech coding, Noise Suppression and on-line system identification. Recently, because of the progress of digital signal processors, a variety of selective coefficient update of gradient-based adaptive algorithms could be implemented in practice. Different types of adaptive algorithms have been developed and used in conventional adaptive filters such as, filtered LMS algorithms [1], [2], [3] and [4], filtered X-LMS algorithms [1], [2], [5], and [6], filtered NLMS algorithms and RLS algorithms [1] and [2]. As a result, this chapter surveys sequential filter adaptation techniques and some applications for transversal FIR filter.

In other words filters are devices or systems that processes (or reshapes) the input signal according to some specific rules to generate an output signal Figure 1.

Figure 1.

Filter Representation

Filters could be classified in two categories linear and non-linear filters where in linear filters there is a linear relationship between input and output of the filter that satisfy the following property:

$$\begin{array}{l}{x}_1\to y_1\\ x_2\to y_2\\ \text{then}a{x}_1+b{x}_2=a{y}_1+b{y}_2\end{array}$$E1

meanwhile in non-linear filters; there is a nonlinearity between input and output of the filter that satisfy the following property:

$$\begin{array}{l}{x}_1\to y_1=x_1^2\\ x_2\to y_2=x_2^2\\ \text{then}{x}_1+x_2={\left(x_1+x_2\right)}^2\end{array}$$E2

In this chapter we will be focusing on adaptive filtering which could automatically adjust (or adapt) in the face of changing environments and changing system requirements by training them to perform specific filtering or decision-making tasks and in which they should be some “adaptation algorithm” for adjusting the system’s parameters figure (2) [7].

Figure 2.

General Adaptive Filter Configuration.

Different filter structures could be implemented in the adaptive filter of figure 2 such as:

• Transversal Filter Structure (FIR)

• IIR Filter Structure

• Lattice Filter Structure

• Non-Linear Filter Structure (Volterra Filter)

in which the adaptation approaches are based on:

• Wiener filter theory where the optimum coefficients of a linear filter are obtained by minimization of its mean-square error (MSE)

• Method of least squares where The performance index is the sum of weighted error squares

The main objective of the adaptation algorithm is to set the filter parameters,$$\overline{{\theta }_{\left(n\right)}}$$ in such a way that its output tries to minimize a meaningful objective function F (cost function in some literatures) involving the reference signal (or desired signal) which will be based on

• Mean Square Error

• Least Square Error

• weighted least square

• instantaneous square value

and the minimization algorithms are based on the following methods:

• Newton’s method

• Quasi-Newton method

• Steepest-descent method

2. Signals

Before pursuing the study of adaptive systems, it is important to refresh our memory with some useful definitions from the stochastic process theory. The representation of signals could fall into two categories:

• Deterministic Signals

• Random Signals

2.1. Deterministic signals

A deterministic discrete-time signal is characterized by a defined mathematical function of the time index n, with n= 0, ±1, ±2, ⋯, such as:

$$x_{\left(n\right)}=e^{-\alpha n}\cos \left(wn\right)+u_{\left(n\right)}$$E3

where u_(n) is the unit-step sequence and The response of a linear time-invariant filter to an input x_(n) is given by:

$$y_{\left(n\right)}=x_{\left(n\right)}\ast h_{\left(n\right)}={\displaystyle \sum _{k=-\infty }^{\infty }{x}_{\left(k\right)}{h}_{\left(n-k\right)}}$$E4

where h_(n) is the impulse response of the filter. Knowing that The Z-transform and its inverse of a given sequence x_(n) is defined as:

$$\begin{array}{l}Z\left(x_{\left(n\right)}\right)=X\left(Z\right)={\displaystyle \sum _{n=-\infty }^{\infty }{x}_{\left(n\right)}{Z}^{-n}}\\ x_{\left(n\right)}=\frac{1}{2\pi j}{\displaystyle \underset{c}{\oint }X\left(Z\right)Z^{n-1}dZ}\end{array}$$E5

where C is a counter clockwise closed contour in the region of convergence of X(z) and encircling the origin of the z-plane as a result; by taking the Z-transform of both sides of equation (4), we will obtain

Y(Z) = H(Z) × X(Z)

For finite finite-energy waveform, it is convenient to use the discrete-time Fourier transform defined as:

$$F\left(x_n\right)=X\left(e^{jw}\right)={\displaystyle \sum _{n=-\infty }^{\infty }{x}_{\left(n\right)}{e}^{-jwn}}$$E6

2.2. Random signals

In many real life situations, observations are made over a period of time and they are influenced by random effects, not just at a single instant but throughout the entire interval of time or sequence of times. In a “rough” sense, a random process is a phenomenon that varies to some degree unpredictably as time goes on. If we observed an entire time-sequence of the process on several different occasions, under presumably “identical” conditions, the resulting observation sequences, in general, would be different. A random variable (RV) is a rule (or function) that assigns a real number to every outcome of a random experiment, while a stochastic random process is a rule (or function) that assigns a time function to every outcome of a random experiment. The elements of a stochastic process, {x_(n)}, for different value of time-index n, are in general complex-valued random variable that are characterized by their probability distribution functions. A stochastic random process is called stationary in the strict sense if all of its (single and joint) distribution function is independent of a shift in the time origin.

In this subsection we will be reviewing some useful definitions in stochastic process:

Stochastic Average

$$m_x=E\left[x_{\left(n\right)}\right]$$E7

where E is the expected value of x_(n).

Autocorrelation function for a stochastic process x_(n)

$${\varphi }_{xx}\left(n,m\right)=E\left[x_{\left(n\right)}{x}_m^{\ast }\right]$$E8

where x^* denotes the complex conjugate of x and the symmetry properties of correlation function is:

$${\varphi }_{xx}\left(k\right)={\varphi }_{xx}^{\ast }\left(-k\right)$$E9

Furthermore a stochastic random process is called stationary in the wide sense if $$m_x$$and $${\varphi }_{xx}\left(n,m\right)$$ are independent of a shift in the time origin for any k, m, and n therefore,

$$m_x\left(n\right)=m_x\left(n+k\right)=\text{constantforall}n$$E10

and

$${\varphi }_{xx}\left(n,m\right)={\varphi }_{xx}\left(n+k,m+k\right)$$E11

which means$${\varphi }_{xx}\left(n,m\right)$$depends only on the difference n – m or

$${\varphi }_{xx}\left(k\right)=E\left[x_{\left(n\right)}{x}_{\left(n-k\right)}^{\ast }\right]$$E12

Special case

$${\varphi }_{xx}\left(0\right)=E\left[{\left|x_{\left(n\right)}\right|}^2\right]=\text{meansquareof}{x}_{\left(n\right)}$$E13

The Z-transform of$${\varphi }_{xx}\left(k\right)$$is given by

$${\Phi }_{xx}\left(Z\right)={\displaystyle \sum _{k=-\infty }^{\infty }{\varphi }_{xx}\left(k\right)Z^{-k}}$$E14

where it can be easily shown that

$${\Phi }_{xx}\left(Z\right)={\Phi }_{xx}^{\ast }\left(\frac{1}{Z^{\ast }}\right)$$E15

Equation (16) implies that if$${\Phi }_{xx}\left(Z\right)$$is a rational function of z, then its poles and zeros must occur in complex-conjugate reciprocal pairs. Moreover, the points that belong to the region of convergence of$${\Phi }_{xx}\left(Z\right)$$also occur in complex-conjugate pairs which suggest that the region of convergence of$${\Phi }_{xx}\left(Z\right)$$must be of the form$$\left|a\right|\prec \left|z\right|\prec \frac{1}{\left|a\right|}$$which will cover the unit circle$$\left|z\right|=1$$. By assuming that$${\Phi }_{xx}\left(Z\right)$$is convergent on the unit circle’s contour C and by letting Z = e^jw then,

$${\varphi }_{xx}\left(0\right)=\frac{1}{2\pi }{\displaystyle \underset{-\pi }{\overset{\pi }{\int }}{\Phi }_{xx}\left(e^{jw}\right)dw}$$E16

and if m_x = 0, we will have

$${\sigma }_x^2=E\left[{\left|x_{\left(n\right)}\right|}^2\right]=\frac{1}{2\pi }{\displaystyle \underset{-\pi }{\overset{\pi }{\int }}{\Phi }_{xx}\left(e^{jw}\right)dw}$$E17

Cross-correlation function for two stochastic processes x_(n) and y_(n)

$${\varphi }_{xy}\left(n,m\right)=E\left[x_{\left(n\right)}{y}_m^{\ast }\right]$$E18

and the symmetry properties of the correlation function is:

$${\varphi }_{xy}\left(k\right)={\varphi }_{yx}^{\ast }\left(-k\right)$$E19

The Z-transform of$${\varphi }_{xy}\left(k\right)$$is given by

$${\Phi }_{xy}\left(Z\right)={\displaystyle \sum _{k=-\infty }^{\infty }{\varphi }_{xy}\left(k\right)Z^{-k}}$$E20

where it can be easily shown that

$${\Phi }_{xy}\left(Z\right)={\Phi }_{xy}^{\ast }\left(\frac{1}{Z^{\ast }}\right)$$E21

Auto-covariance function for stationary process is defined as:

$${\gamma }_{xx}\left(k\right)=E\left[\left(x_{\left(n\right)}-m_x\right){\left(x_{\left(n-k\right)}-m_x\right)}^{\ast }\right]={\varphi }_{xx}\left(k\right)-{\left|m_x\right|}^2$$E22

and the symmetry properties of auto-covariance function

$${\gamma }_{xx}\left(k\right)={\gamma }_{xx}^{\ast }\left(-k\right)$$E23

Special case

$${\gamma }_{xx}\left(0\right)={\sigma }_x^2=\mathrm{var}iance\text{of}{x}_{\left(n\right)}$$E24

Cross-covariance function is defined as:

$${\gamma }_{xy}\left(k\right)=E\left[\left(x_{\left(n\right)}-m_x\right){\left(y_{\left(n-k\right)}-m_y\right)}^{\ast }\right]={\varphi }_{xy}\left(k\right)-m_x{m}_y^{\ast }$$E25

and the symmetry properties of the cross-covariance function is:

$${\gamma }_{xy}\left(k\right)={\gamma }_{yx}^{\ast }\left(-k\right)$$E26

2.2.1. Power spectral density

Consider the wide-sense stationary random process {x_(n)} from which we will consider a window of 2N+1 elements of x_(n) such as

$$x_{N_{\left(n\right)}}=\{\begin{array}{l}{x}_{\left(n\right)},-N\le n\le N\\ \mathrm{0,}otherwise\end{array}$$E27

and the discrete-time Fourier transform of$$x_{N_{\left(n\right)}}$$is computed as:

$$X_N\left(e^{jw}\right)={\displaystyle \sum _{n=-\infty }^{\infty }{x}_{N_{\left(n\right)}}{e}^{-jwn}}={\displaystyle \sum _{n=-N}^N{x}_{\left(n\right)}{e}^{-jwn}}$$E28

By conjugating both sides of equation 29 with:

$$X_N^{\ast }\left(e^{jw}\right)={\displaystyle \sum _{n=-N}^N{x}_{\left(m\right)}^{\ast }{e}^{jwm}}$$E29

we will obtain

$${\left|X_M\left(e^{jw}\right)\right|}^2={\displaystyle \sum _{n=-N}^N{\displaystyle \sum _{m=-N}^N{x}_{\left(n\right)}{x}_{\left(m\right)}^{\ast }{e}^{-jw\left(n-m\right)}}}$$E30

and by taking the expected value of equation 31, we will have:

$$E\left[{\left|X_M\left(e^{jw}\right)\right|}^2\right]={\displaystyle \sum _{n=-N}^N{\displaystyle \sum _{m=-N}^{N}E\left[x_{\left(n\right)}{x}_{\left(m\right)}^{\ast }\right]e^{-jw\left(n-m\right)}}}$$E31

which after simplification by imposing k = n – m will yield:

$$\frac{1}{2N+1}E\left[{\left|X_N\left(e^{jw}\right)\right|}^2\right]={\displaystyle \sum _{k=-2N}^{2N}\left(1-\frac{\left|k\right|}{2N+1}\right){\varphi }_{xx}\left(k\right)e^{-jwk}}$$E32

By assuming that the summation on the right hand side of equation 33 will converge for large k then

$$\underset{N\to \infty }{{\displaystyle \lim }}\left(\frac{1}{2N+1}E\left[{\left|X_N\left(e^{jw}\right)\right|}^2\right]\right)={\displaystyle \sum _{k=-2N}^{2N}{\varphi }_{xx}\left(k\right)e^{-jwk}}={\Phi }_{xx}\left(e^{jw}\right)$$E33

The function$${\Phi }_{xx}\left(e^{jw}\right)$$is called the power spectral density of the stochastic, wide-sense stationary process {x_(n)} which is always real and non-negative.

2.2.2. Response of linear system to stochastic processes

Given the Linear Time-Invariant System (LTI) illustrated in figure 3 where we will assume that$${\varphi }_{xx}\left(k\right)$$is known therefore, the Z-transform of$${\varphi }_{xy}\left(k\right)$$could be elaborated as:

$${\Phi }_{xy}\left(Z\right)={\displaystyle \sum _{k=-\infty }^{\infty }{\varphi }_{xy}\left(k\right)Z^{-k}}={\displaystyle \sum _{k=-\infty }^{\infty }E\left[x_{\left(n\right)}{y}_{\left(n-k\right)}^{\ast }\right]Z^{-k}}$$E34

and

$$y_{\left(n\right)}={\displaystyle \sum _{l=-\infty }^{\infty }{h}_{\left(l\right)}{x}_{\left(n-l\right)}}$$E35

Figure 3.

Structure of LTI system

To be noted that the following expressions could be easily derived:

$${\Phi }_{yy}\left(Z\right)=H\left(Z\right)H^{\ast }\left(\frac{1}{Z}\right){\Phi }_{xx}\left(Z\right)$$E36

and

$$\begin{array}{l}{\Phi }_{xy}\left(e^{jw}\right)=H^{\ast }\left(e^{jw}\right){\Phi }_{xx}\left(e^{jw}\right)\\ {\Phi }_{yx}\left(e^{jw}\right)=H\left(e^{jw}\right){\Phi }_{xx}\left(e^{jw}\right)\\ {\Phi }_{yy}\left(e^{jw}\right)=\left|H\left(e^{jw}\right)\right|{\Phi }_{xx}\left(e^{jw}\right)\end{array}$$E37

3. Regression

Data fitting is one of the oldest adaptive systems which are a powerful tool, allowing predictions of the present or future events to be made based on information about the past or present events. The two basic types of regression are linear regression and multiple regressions.

3.1. Linear regression

The goal of linear regression is to adjust the values of slope and intercept to find the line that best predicts d from x (Figure 4) or in other words linear regression estimates how much d changes when x changes by one unit.

The slope quantifies the steepness of the line which is equals to the change in d for each unit change in x. If the slope is positive, d increases as d increases. If the slope is negative, d decreases as d increases. The d intercept is the d value of the line when x equals zero. It defines the elevation of the line.

Figure 4.

Slope and Intercept.

The deviation from the straight line which represents the best linear fits of a set of data as shown in figure 5 is expressed as:

d ≈ w × x + b

or more specifically,

d_(n) = w × x_(n) +b + e_(n) = y_(n) + e_(n)

where e_(n) is the instantaneous error that is added to y_(n) (the linearly fitted value), w is the slope and b is the y intersect (or bias). More precisely, the goal of regression is to minimize the sum of the squares of the vertical distances of the points from the line. The problem can be solved by a linear system with only two free parameters, the slope w and the bias b.

Figure 5.

Example of linear regression with one independent variable

Figure 6.

Linear Regression Processing Element.

Figure 6 is called the Linear Regression Processing Element (LRPE) which is built from two multipliers and one adder. The multiplier w scales the input, and the multiplier b is a simple bias, which can also be thought of as an extra input connected to the value +1.The parameters (b, w) have different functions in the solution.

3.1.1. Least squares for linear model

Least squares solves the problem by finding the best fitted line to a set of data; for which the sum of the square deviations (or residuals) in the d direction are minimized (figure 7).

Figure 7.

Regression line showing the deviations.

The goal is to find a systematic procedure to find the constants b and w which minimizes the error between the true value d_(n) and the estimated value y_(n) which is called the linear regression.

$$d_{\left(n\right)}–\left(b+w{x}_{(n)}\right)=d_{(n)}– y_{(n)}=e_{(n)}$$E38

The best fitted line to the data is obtained by minimizing the error e_(n) which is computed by the mean square error (MSE) that is a widely utilized as performance criterion

$$\xi =\frac{1}{2N}{\displaystyle \sum _{n=1}^N{e}_{\left(n\right)}^2}$$E39

where N in the number of observations and ξ is the mean square error.

Our goal is to minimize ξ analytically, which can be achieved by taking the partial derivative of this quantity with respect to the unknowns and equate the resulting equations to zero, i.e.

$$\{\begin{array}{c}\frac{\delta \xi }{\delta b}=0\\ \frac{\delta \xi }{\delta w}=0\end{array}$$E40

which yields after some manipulation to:

$$b=\frac{{\displaystyle \sum _n{x}_{\left(n\right)}^2{\displaystyle \sum _n{d}_{\left(n\right)}-{\displaystyle \sum _n{x}_{\left(n\right)}{\displaystyle \sum _n{x}_{\left(n\right)}{d}_{\left(n\right)}}}}}}{N\left[{\displaystyle \sum _n{\left(x_{\left(n\right)}-\overline{x}\right)}^2}\right]}w=\frac{{\displaystyle \sum _n\left(x_{\left(n\right)}-\overline{x}\right)\left(d_{\left(n\right)}-\overline{d}\right)}}{{\displaystyle \sum _n{\left(x_{\left(n\right)}-\overline{x}\right)}^2}}$$E41

where the bar over the variable means the mean value and the procedure to determine the coefficients of the line is called the least square method.

3.1.2. Search procedure

The purpose of least squares is to find parameters (b, w₁, w₂, …, w_p) that minimize the difference between the system output y_(n) and the desired response d_(n). So, regression is effectively computing the optimal parameters of an interpolating system which predicts the value of d from the value of x.

Figure 8 shows graphically the operation of adapting the parameters of the linear system in which the system output y is always a linear combination of the input x with a certain bias b according to the equation y= wx + b. Changing b modifies the y intersect, while changing w modifies the slope. The goal of linear regression is to adjust the position of the line such that the average square difference between the y values (on the line) and the cloud of points d_(n) i.e. the error e_(n), is minimized.

Figure 8.

Regression as a linear system design problem.

The key point is to recognize the information transmitted by the error which can be used to optimally place the line and this could be achieved by including a subsystem that accepts the error and modifies the parameters of the system. Thus, the error e_(n) is fed-back to the system and indirectly affects the output through a change in the parameters (b,w). With the incorporation of the mechanism that automatically modifies the system parameters, a very powerful linear system can be built that will constantly seek optimal parameters. Such systems are called Adaptive systems, and are the focus of this chapter.

3.2. Multiple regression

With the same reasoning as above, the regression for multiple variables could be derived as:

$$e_{\left(i\right)}=d_{\left(i\right)}-\left(b+{\displaystyle \sum _{k=0}^p{w}_{\left(k\right)}{x}_{\left(i,k\right)}}\right)=d_{\left(i\right)}-{\displaystyle \sum _{k=0}^p{w}_{\left(k\right)}{x}_{\left(i,k\right)}}$$E42

where for this case is to find the coefficient vector W that minimizes the MSE of e_(i) over the i samples (Figure 9).

Figure 9.

Regression system for multiple inputs.

The mean square error (MSE) becomes for this case:

$$\xi =\frac{1}{2N}{\displaystyle \sum _n{\left(d_{\left(n\right)}-{\displaystyle \sum _{k=0}^p{w}_{\left(k\right)}{x}_{\left(n,k\right)}}\right)}^2}$$E43

where the solution of this equation can be found by taking the derivatives of ξ with respect to the unknowns (w_(k)), and equating the result to zero which will yield to the famous normal matrix equation expressed as:

$${\displaystyle \sum _n{x}_{\left(n,j\right)}{d}_{\left(n\right)}}={\displaystyle \sum _{k=0}^p{w}_{\left(k\right)}}{\displaystyle \sum _n{x}_{\left(n,k\right)}{x}_{\left(n,j\right)}}$$E44

for j = 0, 1, …, p.

By defining:

$$R_{\left(n,j\right)}=\frac{1}{N}{\displaystyle \sum _n{x}_{\left(n,k\right)}{x}_{\left(n,j\right)}}$$E45

as the autocorrelation function,

$$P_{\left(j\right)}=\frac{1}{N}{\displaystyle \sum _n{x}_{\left(n,j\right)}{d}_{\left(n\right)}}$$E46

as the cross-correlation of the input x for index j and the desired response y and substituting these definitions into Eq. (47), the set of normal equations can be written simply as:

$$P=R{W}^{\ast }\text{or}{W}^{\ast }=R^{-1}P$$E47

where W is a vector with the p+1 weights w_i in which W* represents the value of the vector for the optimum (minimum) solution. The solution of the multiple regression problems can be computed analytically as the product of the inverse of the autocorrelation of the input samples multiplied by the cross-correlation vector of the input and the desired response.

All the concepts previously mentioned for linear regression can be extended to the multiple regression case where J in matrix notation is illustrated as:

$$\xi =\left[W^{T}RW-2P^{T}W+{\displaystyle \sum _n\frac{d_{\left(n\right)}^2}{N}}\right]$$E48

where T means the transpose and the values of the coefficients that minimize the solution are:

$$\nabla \xi =0=RW-P\text{or}{W}^{\ast }=R^{-1}P$$E49

4. Wiener filters

The Wiener filter is a filter proposed by Norbert Wiener during the 1940s and published in 1949 [8]. Its purpose was to reduce the amount of noise present in a signal in comparison with an estimation of the desired noiseless signal. This filter is an MSE-optimal stationary linear filter which was mainly used for images degraded by additive noise and blurring. The optimization of the filter is achieved by minimizing mean square error defined as the difference between the output filter and the desired response (Figure 10) which is known as the cost function expressed as:

$$\xi =E\left[{\left|e_n\right|}^2\right]$$E50

Figure 10.

block schematic of a linear discrete-time filter W(z) for estimating a desired signal d_(n) based on an excitation x_(n) where d_(n) and x_(n) are random processes.

In signal processing, a causal filter is a linear and time-invariant causal system. The word causal indicates that the filter output depends only on past and present inputs. A filter whose output also depends on future inputs is non-causal. As a result two cases should be considered for the optimization of the cost function (equation 53).

• The filter W(Z) is causal and or FIR (Finite Impulse Response).

• The filter W(Z) is non-causal and or IIR (Infinite Impulse Response Filter).

4.1. Wiener filter – the transversal filter

Let$$\overline{W}$$be the transversal filter’s tap weights illustrated in figure 11 which is defined as:

$$\overline{W}={\left[\begin{array}{cccc}{w}_0& w_1& \cdots & w_{N-1}\end{array}\right]}^T$$E51

where T denotes the vector transpose and

$$\overline{X}={\left[\begin{array}{cccc}{x}_0& x_1& \cdots & x_{N-1}\end{array}\right]}^T$$E52

be the input vector signal where two cases should be treated separately depending on the required application:

• Real valued input signal

• Complex valued input signal

Figure 11.

Block Diagram of the Transversal Filter.

4.1.1. Wiener filter – the transversal filter - real-valued input signal

The output filter could be expressed as:

$$y_{\left(n\right)}={\displaystyle \sum _{i=0}^{N-1}{w}_i{x}_{\left(n-i\right)}}={\overline{W}}^T\overline{X_{\left(n\right)}}={\overline{X_{\left(n\right)}}}^T\overline{W}$$E53

Where

$$e_{\left(n\right)}=d_{\left(n\right)}-y_{\left(n\right)}=d_{\left(n\right)}-{\overline{X_{\left(n\right)}}}^T\overline{W}$$E54

The performance or cost function expressed in equation 53 could be elaborated as:

$$\begin{array}{l}\xi =E\left[{\left|e_n\right|}^2\right]=E\left[\left(d_{\left(n\right)}-{\overline{W}}^T\overline{X_{\left(n\right)}}\right)\left(d_{\left(n\right)}-{\overline{X_{\left(n\right)}}}^T\overline{W}\right)\right]\\ =E\left[d_{\left(n\right)}^2\right]-{\overline{W}}^{T}E\left[\overline{X_{\left(n\right)}}{d}_{\left(n\right)}\right]-E\left[{\overline{X_{\left(n\right)}}}^T{d}_{\left(n\right)}\right]\overline{W}+{\overline{W}}^{T}E\left[\overline{X_{\left(n\right)}}{\overline{X_{\left(n\right)}}}^T\right]\overline{W}\end{array}$$E55

By examining equation 58 we could easily notice that$$E\left[\overline{X_{\left(n\right)}}{d}_{\left(n\right)}\right]$$is the cross correlation function which will be defined as:

$$\overline{P}\equiv E\left[\overline{X_{\left(n\right)}}{d}_{\left(n\right)}\right]={\left[\begin{array}{cccc}{p}_0& p_1& \cdots & p_{N-1}\end{array}\right]}^T$$E56

and with the same reasoning as above the autocorrelation function$$E\left[\overline{X_{\left(n\right)}}{\overline{X_{\left(n\right)}}}^T\right]$$could be expressed as:

$$R\equiv E\left[\overline{X_{\left(n\right)}}{\overline{X_{\left(n\right)}}}^T\right]=\left[\begin{array}{cccc}{r}_{\mathrm{0,0}}& r_{\mathrm{0,1}}& \cdots & r_{\mathrm{0,}N-1}\\ r_{\mathrm{1,0}}& r_{\mathrm{1,1}}& \cdots & r_{\mathrm{1,}N-1}\\ ⋮& ⋮& ⋮& ⋮\\ r_{N-\mathrm{1,0}}& r_{N-\mathrm{1,1}}& \cdots & r_{N-\mathrm{1,}N-1}\end{array}\right]$$E57

Knowing that$$E\left[d_{\left(n\right)}{\overline{X_{\left(n\right)}}}^T\right]={\overline{P}}^T$$and$${\overline{W}}^T\overline{P}={\overline{P}}^T\overline{W}$$therefore, the cost function would be:

$$\xi =E\left[d_{\left(n\right)}^2\right]-2{\overline{W}}^T\overline{P}+{\overline{W}}^{T}R\overline{W}$$E58

where$$\overline{W}$$in the quadratic function expressed in equation 47 will have a global minimum if and only if R is a positive definite Matrix.

The gradient method is the most commonly used method to compute the tap weights that minimizes the cost function therefore,

$$\nabla \xi =\frac{\partial \xi }{\partial w_i}=\left[\begin{array}{cccc}\frac{\partial \xi }{\partial w_0}& \frac{\partial \xi }{\partial w_1}& \cdots & \frac{\partial \xi }{\partial w_{N-1}}\end{array}\right]=0$$E59

which yield after expanding equation 61 to:

$$\begin{array}{l}\frac{\partial \xi }{\partial w_i}=0\\ 0=-2p_i+{\displaystyle \sum _{l=0}^{N-1}{w}_l\left(r_{li}+r_{il}\right)}\\ {\displaystyle \sum _{l=0}^{N-1}{w}_l\left(r_{li}+r_{il}\right)}=2p_i\end{array}$$E60

where we have assumed in the expanded equation$$\xi =E\left[d_{\left(n\right)}^2\right]-2{\displaystyle \sum _{l=0}^{N-1}{p}_l{w}_l}+{\displaystyle \sum _{l=0}^{N-1}{\displaystyle \sum _{m=0}^{N-1}{w}_l{w}_l{r}_{lm}}}$$that:

$${\displaystyle \sum _{l=0}^{N-1}{\displaystyle \sum _{m=0}^{N-1}{w}_l{w}_l{r}_{lm}}}={\displaystyle \sum _{\begin{array}{l}l=0\\ l\ne i\end{array}}^{N-1}{\displaystyle \sum _{\begin{array}{l}m=0\\ m\ne i\end{array}}^{N-1}{w}_l{w}_l{r}_{lm}}}+w_i{\displaystyle \sum _{\begin{array}{l}m\ne 0\\ l\ne i\end{array}}^{N-1}{w}_m{r}_{im}+w_i^2{r}_{ii}}$$E61

Knowing that r_li = r_il due to the symmetry property of the autocorrelation function of real-valued signal equation 49 could be expressed as:

$${\displaystyle \sum _{l=0}^{N-1}{w}_l{r}_{il}}=p_i$$E62

for i = 0, 1, …, N – 1 which could be expressed in matrix notation as:

$$R{\overline{W}}_{op}=\overline{P}$$E63

known as Weiner-Hopf equation, whose solution for the optimal tap-weight vector$${\overline{W}}_{op}$$and by assuming that R has an inverse matrix, is:

$${\overline{W}}_{op}=R^{-1}\overline{P}$$E64

Finally, the minimum value of the cost function can be expressed as:

$${\xi }_{\min }=E\left[d_{\left(n\right)}^2\right]-{\overline{W}}_{op}^T\overline{P}=E\left[d_{\left(n\right)}^2\right]-{\overline{W}}_{op}^{T}R{\overline{W}}_{op}=E\left[d_{\left(n\right)}^2\right]-{\overline{P}}^T{R}^{-1}\overline{P}$$E65

4.1.2. Wiener filter – the transversal filter - complex-valued input signal

The data transmission of the basebands QPSK and the QAM is complex valued random signals where the filter’s tap weight vector is also assumed to be complex therefore, the cost function for this case could be expressed as:

$$\xi =E\left[{\left|e_{\left(n\right)}\right|}^2\right]=E\left[e_{\left(n\right)}{e}_{\left(n\right)}^{\ast }\right]$$E66

and the gradient of the cost function would be:

$${\nabla }_{w_i}^c=E\left[e_{\left(n\right)}{\nabla }_{w_i}^c{e}_{\left(n\right)}^{\ast }+{\nabla }_{w_i}^c{e}_{\left(n\right)}^{\ast }{e}_{\left(n\right)}\right]$$E67

By examining figure 11 we could easily notice that equation 57 could be formulated as:

$$e_{\left(n\right)}=d_{\left(n\right)}-{\displaystyle \sum _{k=0}^{N-1}{w}_{\left(k\right)}{x}_{\left(n-k\right)}}$$E68

and since d is independent of w and since

$$\begin{array}{l}{\nabla }_{w_i}^c{e}_{\left(n\right)}=-x_{\left(n-i\right)}{\nabla }_{w_i}^c{w}_i\\ {\nabla }_{w_i}^c{e}_{\left(n\right)}^{\ast }=-x_{\left(n-i\right)}^{\ast }{\nabla }_{w_i}^c{w}_i^{\ast }\\ {\nabla }_{w_i}^c{w}_i=\frac{\partial w_i}{\partial w_{i,R}}+j\frac{\partial w_i}{\partial w_{i,I}}=1+j\left(j\right)=0\\ {\nabla }_{w_i}^c{w}_i^{\ast }=\frac{\partial w_i^{\ast }}{\partial w_{i,R}^{\ast }}+j\frac{\partial w_i^{\ast }}{\partial w_{i,I}^{\ast }}=1+j\left(-j\right)=2\end{array}$$E69

Where the sub-indices R and I refers to the real part and imaginary part of the complex number therefore, the gradient of the cost function (Equation 70) would be:

$${\nabla }_{w_i}^c\xi =-2E\left[e_{\left(n\right)}{x}_{\left(n-i\right)}^{\ast }\right]$$E70

The optimum filter’s tap-weights e_0(n) are obtained by setting the complex gradient of equation 73 to zero yielding:

$$E\left[e_{o\left(n\right)}{x}_{\left(n-i\right)}^{\ast }\right]=0$$E71

By defining

$$\begin{array}{l}{\overline{x}}_{\left(n\right)}\equiv {\left[\begin{array}{cccc}{x}_{\left(n\right)}& x_{\left(n-1\right)}& \cdots & x_{\left(n-N+1\right)}\end{array}\right]}^T\equiv {\left[\begin{array}{cccc}{x}_{\left(n\right)}^{*}& x_{\left(n-1\right)}^{*}& \cdots & x_{\left(n-N+1\right)}^{*}\end{array}\right]}^H\\ {\overline{w}}_{\left(n\right)}\equiv {\left[\begin{array}{cccc}{w}_0^{*}& w_1^{*}& \cdots & w_{N-1}^{*}\end{array}\right]}^T\equiv {\left[\begin{array}{cccc}{w}_{\left(0\right)}& w_1& \cdots & x_{N-1}\end{array}\right]}^H\end{array}$$E72

where H denotes the complex-conjugate transpose or Hermitian and by re-writing equation 71 as:

$$e_{\left(n\right)}=d_{\left(n\right)}-{\overline{w}}_o^H{x}_{\left(n\right)}$$E73

whose solution is:

$$R{\overline{w}}_o=\overline{p}$$E74

where

$$\begin{array}{l}R=E\left[{\overline{x}}_{\left(n\right)}{\overline{x}}_{\left(n\right)}^H\right]\\ \overline{p}=E\left[{\overline{x}}_{\left(n\right)}{d}_{\left(n\right)}^{\ast }\right]\end{array}$$E75

Equation 77 is known as Weiner-Hopf equation for the complex valued signals and the minimum of the cost function will be:

$${\xi }_{\min }=E\left[d_{\left(n\right)}^2\right]-{\overline{w}}_o^{H}R{\overline{w}}_o$$E76

5. Least Mean Square algorithm (LMS algorithm)

The purpose of least squares is to is to find the optimal filter’s tap weights that that minimize the difference between the system output y_(n) and the desired response d_(n) or in other words to minimize the cost function. Instead of solving the Weiner-Hopf equation in order to obtain the optimal filter’s tap weights as seen previously, their exists other iterative methods which employ an iterative search method that starts with an arbitrary initial weight$${\overline{W}}_{\left(0\right)}$$, then a recursive search method that may require many iterations in order to converge to the optimal filter’s tap weights$${\overline{W}}_o$$. The most important iterative methods are the gradient based iterative methods which are listed below:

• Steepest Descent Method

• Newton’s Method

5.1. Steepest descent method

The steepest descent method (also known as the gradient method) is the simplest example of a gradient based method that minimizes a function of several variables. This process is employing an iterative search method that starts with an arbitrary initial weight$${\overline{W}}_{\left(0\right)}$$, and then at the k^th iteration$${\overline{W}}_{\left(k\right)}$$is updated according to the following equation:

$${\overline{W}}_{\left(k+1\right)}={\overline{W}}_{\left(k\right)}-\mu {\nabla }_k\xi$$E77

where µ is positive known as the step size and$${\nabla }_k\xi$$ denotes the gradient vector $$\nabla \xi =2R\overline{W}-2\overline{p}$$evaluated at the point$$\overline{W}={\overline{W}}_{\left(k\right)}$$. Therefore, equation 80 could be formulated as:

$${\overline{W}}_{\left(k+1\right)}={\overline{W}}_{\left(k\right)}-2\mu \left(R{\overline{W}}_{\left(k\right)}-{\overline{p}}_{\left(k\right)}\right)=\left(1-2\mu R\right)\left({\overline{W}}_{\left(k\right)}-{\overline{W}}_0\right)$$E78

Knowing that the auto-correlation matrix R may be diagonalized by using the unitary similarity decomposition:

$$R=Q\Lambda Q^T$$E79

where$$\Lambda$$is a diagonal matrix consisting of the auto-correlation matrix R eigenvalues $${\lambda }_0{\lambda }_1\cdots {\lambda }_{N-1}$$ and the columns of Q contain the corresponding orthonormal eigenvectors and by defining the vector$${\overline{v}}_{\left(k\right)}$$as:

$${\overline{v}}_{\left(k\right)}={\overline{W}}_{\left(k\right)}-{\overline{W}}_0$$E80

As a result equation 81 could be expressed as:

$${\overline{v}}_{\left(k+1\right)}=\left(I-2\mu Q\Lambda Q^T\right){\overline{v}}_{\left(k\right)}$$E81

which after mathematical manipulation and by using the vector transformation $${\overline{v}}_{\left(k\right)}^{\text{'}}=Q^T{\overline{v}}_{\left(k\right)}$$could be written as:

$$v_{i_{\left(k+1\right)}}^{\text{'}}=\left(1-2\mu {\lambda }_i\right)v_{i_{\left(k\right)}}^{\text{'}}$$E82

For i = 0, 1, 2,..., N – 1 and$${\overline{v}}_{\left(k\right)}^{\text{'}}=\left[\begin{array}{cccc}{v}_{0_{\left(k\right)}}^{\text{'}}& v_{1_{\left(k\right)}}^{\text{'}}& \cdots & v_{N-1_{\left(k\right)}}^{\text{'}}\end{array}\right]$$that yields:

$$v_{i_{\left(k\right)}}^{\text{'}}={\left(1-2\mu {\lambda }_i\right)}^k{v}_{i_{\left(0\right)}}^{\text{'}}$$E83

Equation 86 converges to zero if and only if the step-size parameter μ is selected so that:

$$\left|1-2\mu {\lambda }_i\right|\prec 1$$E84

which means

$$0\prec \mu \prec {\lambda }_i$$E85

for all i or equivalently

$$0\prec \mu \prec {\lambda }_{\max }$$E86

where λ_max is the maximum of the eigenvalues$${\lambda }_0{\lambda }_1\cdots {\lambda }_{N-1}$$.

5.2. Newton’s method

By replacing the scalar step-size μ with a matrix step-size given by μR^-1in the steepest descent algorithm (equation 81) and by using$$\overline{p}=R{\overline{W}}_0$$, the resulting algorithm is:

$${\overline{W}}_{\left(k+1\right)}={\overline{W}}_{\left(k\right)}-\mu R^{-1}{\nabla }_k\xi$$E87

Substituting$$\nabla \xi =2R\overline{W}-2\overline{p}$$in equation 90 we will get:

$$\begin{array}{l}{\overline{W}}_{\left(k+1\right)}={\overline{W}}_{\left(k\right)}-2\mu R^{-1}\left(R{\overline{W}}_{\left(k\right)}-\overline{p}\right)\\ =\left(1-2\mu \right){\overline{W}}_{\left(k\right)}+2\mu R^{-1}\overline{p}\\ =\left(1-2\mu \right){\overline{W}}_{\left(k\right)}+2\mu {\overline{W}}_0\end{array}$$E88

and by subtracting$${\overline{W}}_0$$ from both sides of equation 91 will yield:

$${\overline{W}}_{\left(k+1\right)}-{\overline{W}}_0=\left(1-2\mu \right)\left({\overline{W}}_{\left(k\right)}-{\overline{W}}_0\right)$$E89

where in actual implementation of adaptive filters, the exact values of$${\nabla }_k\xi$$ and R ^{– 1} are not available and have to be estimated.

5.3. Least Mean Square (LMS) algorithm

The LMS algorithm is a practical scheme for realizing Wiener filters, without explicitly solving the Wiener-Hopf equation and this was achieved in the late 1960’s by Widrow [2] who proposed an extremely elegant algorithm to estimate the gradient that revolutionized the application of gradient descent procedures by using the instantaneous value of the gradient as the estimator for the true quantity which means replacing the cost function$$\xi =E\left[e_{\left(n\right)}^2\right]$$by its instantaneous coarse estimate$$\widehat{\xi }=e_{\left(n\right)}^2$$into steepest descent algorithm. By doing so, we will obtain

$${\overline{W}}_{\left(n+1\right)}={\overline{W}}_{\left(n\right)}-\mu \nabla e_{\left(n\right)}^2$$E90

where$${\overline{W}}_{\left(n\right)}={\left[\begin{array}{cccc}{w}_{0_{\left(n\right)}}& w_{1_{\left(n\right)}}& \cdots & w_{N-1_{\left(n\right)}}\end{array}\right]}^T$$and

$$\nabla ={\left[\begin{array}{cccc}\frac{\partial }{\partial w_0}& \frac{\partial }{\partial w_1}& \cdots & \frac{\partial }{\partial w_{N-1}}\end{array}\right]}^T$$E91

.

The i^th element of the gradient vector$$\nabla e_{\left(n\right)}^2$$is:

$$\frac{\partial e_{\left(n\right)}^2}{\partial w_i}=2e_{\left(n\right)}\frac{\partial e_{\left(n\right)}}{\partial w_i}=-2e_{\left(n\right)}\frac{\partial y_{\left(n\right)}}{\partial w_i}=-2e_{\left(n\right)}{x}_{\left(n-i\right)}$$E92

which leads to:

$$\nabla e_{\left(n\right)}^2=-2e_{\left(n\right)}{\overline{x}}_{\left(n\right)}$$E93

that will be replaced in equation 93 to obtain:

$${\overline{W}}_{\left(n+1\right)}={\overline{W}}_{\left(n\right)}+2\mu e_{\left(n\right)}{\overline{x}}_{\left(n\right)}$$E94

where

$${\overline{x}}_{\left(n\right)}={\left[\begin{array}{cccc}{x}_{\left(n\right)}& x_{\left(n-1\right)}& \cdots & x_{\left(n-N+1\right)}\end{array}\right]}^T$$E95

.

Figure 12.

LMS Filter Structure

Equation 96 is known as the LMS recursion and the summary of the LMS algorithm illustrated on figure 12 will be as follow:

1. Input

2. tap-weight vector,$${\overline{W}}_{\left(n\right)}$$

3. input vector, $${\overline{x}}_{\left(n\right)}$$

4. desired output, d_(n)

5. Output

6. Filter output, y_(n)

7. Tap-weight vector update, $${\overline{W}}_{\left(n+1\right)}$$

8. Filtering:$$y_{\left(n\right)}={\overline{W}}_{\left(n\right)}^T{\overline{x}}_{\left(n\right)}$$

9. Error estimation: e_(n) = d_(n) - y_(n)

10. Tap-weight vector adaptation:$${\overline{W}}_{\left(n+1\right)}={\overline{W}}_{\left(n\right)}+2\mu e_{\left(n\right)}{\overline{x}}_{\left(n\right)}$$

6. Classifying adaptive filtering applications

Various applications of adaptive filtering differ in the manner in which the desired response is extracted. In this context, we may distinguish four basic classes of adaptive filtering applications (depicted in Figures 13 to 16, which follow):

• Identification

• Inverse Modeling

• Prediction

• Interference Cancelling

Adaptive Filtering Class	Application
Identification	System Identification; Layered Earth Modeling
Inverse Modeling	Predictive Convolution; Adaptive Equalization
Prediction	Linear Prediction Coding; Adaptive Differential PCM; Auto Regressive Spectrum Analysis; Signal Detection
Interference Cancelling	Adaptive Noise Cancelling; Echo Cancellation; Radar Polarimetry; Adaptive Beam-forming

Table 1.

Adaptive Filtering Applications

The following notations are used in Figures 11-15:

u = input applied to the adaptive filter

y = output of the adaptive filter

d = desired response

e = d – y = estimation error

The functions of the four basic classes of adaptive filtering applications are summarized as follow:

Identification (Figure 13).

The notion of a mathematical model is fundamental to sciences and engineering. In the class of applications dealing with identification, an adaptive filter is used to provide a linear model that represents the best fit to an unknown plant. The plant and the adaptive filter are driven by the same input. The plant output supplies the desired responses for the adaptive filter. If the plant is dynamic in nature, the model will be time varying.

Figure 13.

Identification

Inverse Modeling (Figure 14).

In this second class of applications, the adaptive filter provides an inverse model representing the best fit to an unknown noisy plant. Ideally, the inverse model has a transfer function equal to the reciprocal of the plant’s transfer function. A delayed version of the plant input constitutes the desired response for the adaptive filter. In some applications, the plant input is used without delay as the desired response.

Figure 14.

Inverse Modeling

Prediction (Figure 15).

In this example, the adaptive filter provides the best prediction of the present value of a random signal. The present value of the signal serves the purpose of a desired response for the adaptive filter. Past values of the signal supply the input applied to the adaptive filter. Depending on the application of interest, the adaptive filter output or the estimation error may service as the system output. In the first case, the system operates as a predictor; in the latter case, it operates as a prediction error filter.

Figure 15.

Prediction

Interference Cancelling (Figure 16).

In this final class of applications, the adaptive filter is used to cancel unknown interference contained in a primary signal, with the cancellation being optimized in some sense. The primary signal serves as the desired response for the adaptive filter. A reference signal is employed as the input to the adaptive filter. The reference signal is derived from a sensor or set of sensors located in relation to the sensor(s) supplying the primary signal in such a way that the information-bearing signal component is weak or essentially undetectable.

Figure 16.

Interference Cancelling

7. Active noise suppressor ANC

In this subsection we propose a new approach of noise control named ANS, where we propose stability-guaranteed algorithm for an adaptive filters, which can be derived on a basis of the strictly positive real property of the error model treated in adaptive system theory. It is important to assure the stability of the adaptive system especially in presence of unknown disturbances and mismatch in the order of the adaptive filter. Experimental results, performed on real mining noise, validate the effectiveness of the proposed stable algorithm [9].

7.1. Sounds and hearing

When sound waves from a point source strike a plane wall, they produce reflected circular wave fronts as if there were an "image" of the sound source at the same distance on the other side of the wall. If something obstructs the direct sound from the source from reaching your ear, then it may sound as if the entire sound is coming from the position of the "image" behind the wall. This kind of sound imaging follows the same laws of reflection as your image in a plane mirror (figure 17). The reflection of sound follows the law "angle of incidence equals angle of reflection" as the light does and other waves or the bounce of a billiard ball off the bank of a table figure (18).

Figure 17.

Point source of sound reflecting from a plane surface.

Figure 18.

Wave reflection

The main item of note about sound reflections off of hard surfaces is the fact that they undergo a 180-degree phase change upon reflection. This can lead to resonance such as standing waves in rooms. It also means that the sound intensity near a hard surface is enhanced because the reflected wave adds to the incident wave, giving pressure amplitude that is twice as great in a thin "pressure zone" near the surface. This is used in pressure zone microphones to increase sensitivity. The doubling of pressure gives a 6-decibel increase in the signal picked up by the microphone figure (19). Since the reflected wave and the incident wave add to each other while moving in opposite directions, the appearance of propagation is lost and the resulting vibration is called a standing wave. The modes of vibration associated with resonance in extended objects like strings and air columns have characteristic patterns called standing waves. These standing wave modes arise from the combination of reflection and interference such that the reflected waves interfere constructively with the incident waves. An important part of the condition for this constructive interference is the fact that the waves change phase upon reflection from a fixed end. Under these conditions, the medium appears to vibrate in segments or regions and the fact that these vibrations are made up of traveling waves is not apparent - hence the term "standing wave" figure (19).

Figure 19.

Resultant wave in pressure zone.

Two traveling waves, which exist in the same medium, will interfere with each other figure (20). If their amplitudes add, the interference is said to be constructive interference, and destructive interference if they are "out of phase" and subtract figure (21). Patterns of destructive and constructive interference may lead to "dead spots" and "live spots" in auditorium acoustics.

Interference of incident and reflected waves is essential to the production of resonant standing waves figure (22).

Figure 20.

Interference of Sound

Figure 21.

Interference and Phase

Figure 22.

The fundamental and second harmonic standing waves for a stretched string.

The sound intensity from a point source of sound will obey the inverse square law if there are no reflections or reverberation figure (23). Any point source, which spreads its influence equally in all directions without a limit to its range, will obey the inverse square law. This comes from strictly geometrical considerations. The intensity of the influence at any given radius r is the source strength divided by the area of the sphere. Being strictly geometric in its origin, the inverse square law applies to diverse phenomena. Point sources of gravitational force, electric field, light, sound or radiation obey the inverse square law.

Figure 23.

Inverse Square Law

A plot of this intensity drop shows that it drops off rapidly figure (24). A plot of the drop of sound intensity according to the inverse square law emphasizes the rapid loss associated with the inverse square law. In an auditorium, such a rapid loss is unacceptable. The reverberation in a good auditorium mitigates it. This plot shows the points connected by straight lines but the actual drop is a smooth curve between the points.

Figure 24.

Inverse Square Law Plot

Reverberation is the collection of reflected sounds from the surfaces in an enclosure like an auditorium. It is a desirable property of auditoriums to the extent that it helps to overcome the inverse square law drop-off of sound intensity in the enclosure. However, if it is excessive, it makes the sounds run together with loss of articulation - the sound becomes muddy, garbled figure (25).

Figure 25.

Reverberant sound is the collection of all the reflected sounds in an auditorium.

7.2. ANC

In order to cancel unwanted noise, it is necessary to obtain an accurate estimate of the noise to be cancelled. In an open environment, where the noise source can be approximated as a point source, microphones can be spaced far apart as necessary and each will still receive a substantially similar estimate of the background noise. However in a confined environment containing reverberation noise caused by multiple sound reflections, the sound field is very complex and each point in the environment has a very different background noise signal. The further apart the microphones are, the more dissimilar the sound field. As a result, it is difficult to obtain an accurate estimate of the noise to be cancelled in a confined environment by using widely spaced microphones.

The proposed model is embodied in a dual microphone noise suppression system in which the echo between the two microphones is substantially cancelled or suppressed. Reverberations from one microphone to the other are cancelled by the use of first and second line echo cancellers. Each line echo canceller models the delay and transmission characteristics of the acoustic path between the first and second microphones Figure (26).

Figure 26.

A pictorial representation of the sound field reaching an ear set in accordance with the proposed model intended to be worn in the human ear.

If the two microphones are moved closer together, the second microphone should provide a better estimate of the noise to be cancelled in the first microphone. However, if the two microphones are placed very close together, each microphone will cause an additional echo to strike the other microphone. That is, the first microphone will act like a speaker (a sound source) transmitting an echo of the sound field striking the second microphone. Similarly, the second microphone will act like a speaker (a sound source) transmitting an echo of the sound field striking the first microphone. Therefore, the signal from the first microphone contains the sum of the background noise plus a reflection of the background noise, which results in a poorer estimate of the background noise to be cancelled figures (27) and (28).

Figure 27.

Simplified model for one direction.

Figure 28.

Simplified model for the opposite direction.

In a first embodiment, a noise suppression system in accordance with the proposed model acts as an ear protector, cancelling substantially all or most of the noise striking the dual microphones of the ear set figure (29). In a second embodiment, a noise suppression system in accordance with the present invention acts a noise suppression communication system, suppressing background noise while allowing communication signals to be heard by the wearer figures (30 and 31).

Figure 29.

A noise suppression System in accordance with a first embodiment of the proposed model.

Figure 30.

A noise suppression communications system in accordance with a second embodiment of the proposed model.

Figure 31.

An alternate scheme for a noise suppression communications system in accordance with a second embodiment of the proposed model.

Figure 32.

Simulink block diagram of the proposed noise suppressor in figure 29.

The conceptual key to this proposed model is that the signals received at two closely spaced microphones in a multi-path acoustic environment are each made up of a sum of echoes of the signal received at the other one. This leads to the conclusion that the difference between the two microphone signals is a sum of echoes of the acoustic source in the environment. In the absence of a speech source, the ANS scheme proposed by Jaber first attempts to isolate the difference signal at each of the microphones by subtracting from it an adaptively predicted version of the other microphone signal. It then attempts to adaptively cancel the two difference signals. When speech is present (as detected by some type of VAD-based strategy), the adaptive cancellation stage has its adaptivity turned off (i.e. the impulse responses of the two FIR filters, one for each microphone, are unchanged for the duration of the speech). The effect here is that the adaptive canceller does not end up cancelling the speech signal contained in the difference between the two microphone signals.

The Simulink implementation of the noise suppression system illustrated in figure 29 is displayed in figure 32 where figures 33 and 34 display the noise captured by Mic.1 and Mic. 2 respectively, meanwhile figure 35 shows the filter output.

Figure 33.

Noise Captured by Mic. 1.

Figure 34.

Noise Captured By Mic. 2.

The Simulink implementation of the Active Noise Suppressor illustrated in Figure 30 with no Voice Activity Detection (VAD) attached is sketched in 36 and the taking off plane’s noise captured by Mic.1 and Mic. 2 are presented in figures 37 and 38 respectively meanwhile figure 39 shows the system output.

Figure 35.

The Filtered Noise output.

With an accurate front-point and end-point detection VAD algorithm [10] implemented in our proposed models illustrated in figures 30-31, an active noise suppressor is obtained.

The simulation results have been obtained from 25 seconds of noisy voice conditions as shown in Fig. 40 at Mic. 1 and the clean speech of our output filter is shown in figure 41 mean while figure 42 will illustrate the clean speech on a reduced scale [11].

Figure 36.

Simulink block diagram of the proposed noise suppressor in figure 30

Figure 37.

Taking off Airplane’s Noise Captured by Mic. 1.

Figure 38.

Taking off Airplane’s Noise Captured by Mic. 2.

Figure 39.

The system’s output

Figure 40.

Noisy speech captured by Mic 1.

Figure 41.

The clean speech obtained at the output of our proposed ANC (Fig. 30) by using a VAD.

Figure 42.

The clean speech obtained at the output of our proposed ANC (Fig. 30) by reducing the time scale

8. The ultra high speed LMS algorithm implemented on parallel architecture

There are many problems that require enormous computational capacity to solve, and therefore the success of computational science to accurately describe and model the real world has helped to fuel the ever increasing demand for cheap computing power. Scientists are eager to find ways to test the limits of theories, using high performance computing to allow them to simulate more realistic systems in greater detail. Parallel computing offers a way to address these problems in a cost effective manner. Parallel Computing deals with the development of programs where multiple concurrent processes cooperate in the fulfilment of a common task. Finally, in this section we will develop the theory of the parallel computation of the widely used algorithms named the least-mean-square (LMS) algorithm

M Jaber “Method and apparatus for enhancing processing speed for performing a least mean square operation by parallel processing” US patent No. 7,533,140, 2009

by its originators, Widrow and Hoff (1960) [2].

8.1. The spatial radix-r factorization

This section will be devoted in proving that discrete signals could be decomposed into r partial signals and whose statistical properties remain invariant therefore, given a discrete signal x_n of size N

$$x_n=\left[x_0{x}_1{x}_2\cdots x_{\left(N-1\right)}\right]$$E96

and the identity matrix I_r of size r

$$I_r=\left[\begin{array}{cccc}1& 0& \mathrm{...}& 0\\ 0& 1& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ 0& 0& \mathrm{...}& 1\end{array}\right]=\left[I_{\left(l,c\right)}\right]=\{\begin{array}{l}1forl=c\\ 0elsewhere\end{array}$$E97

for l = c = 0, 1,..., r – 1.

Based on what was proposed in [2]-[9]; we can conclude that for any given discrete signal x_(n) we can write:

$$x_{\left(l,n\right)}=\left[I_{{}_{\left(l,c\right)}}\right]\times \left[x_{\left(rn\right)}{x}_{\left(rn+1\right)}\cdots x_{\left(rn+p\right)}\right]$$E98

is the product of the identity matrix of size r by r sets of vectors of size N/r (n = 0,1,.., N/r -1) where the l^th element of the n^th product is stored into the memory address location given by

l = (rn + p)

for p = 0, 1, …, r – 1.

The mean (or Expected value) of x_n is given as:

$$E\left[x\right]=\frac{{\displaystyle \sum _{n=0}^{N-1}{x}_n}}{N}$$E99

which could be factorizes as:

$$\begin{array}{l}E\left[x\right]=\frac{{\displaystyle \sum _{n=0}^{\left(N/r\right)-1}{x}_{\left(rn\right)}}+\cdots +{\displaystyle \sum _{n=0}^{\left(N/r\right)-1}{x}_{\left(rn+\left(r-1\right)\right)}}}{N}={\mu }_{\left(x\right)}\\ \text{=}\frac{\frac{1}{r}\times {\displaystyle \sum _{n=0}^{\left(N/r\right)-1}{x}_{\left(rn+p\right)}}}{\frac{N}{r}}\text{=}\frac{{\mu }_{\left(x_{\left(rn+p\right)}\right)}}{r}\end{array}$$E100

therefore, the mean of the signal x_n is equal to sum of the means of its r partial signals divided by r for p = 0, 1, …, r – 1.

Similarly to the mean, the variance of the signal x_n equal to sum of the variances of its r partial signals according to:

$$\begin{array}{l}Var\left(x\right)={\sigma }_x^2=E\left[{\left(x-\mu \right)}^2\right]={\displaystyle \sum _{n=0}^{N-1}{\left(x_{\left(n\right)}-\mu \right)}^2}\\ \text{=}{\displaystyle \sum _{n=0}^{\left(N/r\right)-1}{\left(x_{\left(rn\right)}-{\mu }_{\left(rn\right)}\right)}^2}+\cdots +{\displaystyle \sum _{n=0}^{\left(N/r\right)-1}{\left(x_{\left(rn+\left(r-1\right)\right)}-{\mu }_{\left(rn+\left(r-1\right)\right)}\right)}^2}\\ \text{=}{\displaystyle \sum _{p=0}^{r-1}{\sigma }_{x_{\left(rn+p\right)}}^2}\end{array}$$E101

8.2. The parallel implementation of the least squares method

The method of least squares assumes that the best-fit curve of a given type is the curve that has the minimal sum of the deviations squared (least square error) from a given set of data. Suppose that the N data points are (x₀, y₀), (x₁, y₁)… (x_{(n – 1)}, y_{(n – 1)}), where x is the independent variable and y is the dependent variable. The fitting curve d has the deviation (error) σ from each data point, i.e., σ₀ = d₀ – y₀, σ₁ = d₁ – y₁... σ_{(n – 1)} = d_{(n – 1)} – d_{(n – 1)} which could be re-ordered as:

$$\begin{array}{l}{\sigma }_0=d_0-y_0,{\sigma }_r=d_r-y_r,{\sigma }_{2r}=d_{2r}-y_{2r}\mathrm{,...,}\\ {\sigma }_{rn}=d_{rn}-y_{rn},{\sigma }_1=d_1-y_1\mathrm{,...,}{\sigma }_{\left(rn+1\right)}=d_{\left(rn+1\right)}-y_{\left(rn+1\right)}\mathrm{,...,}{\sigma }_{\left(rn+\left(r-1\right)\right)}=d_{\left(rn+\left(r-1\right)\right)}-y_{\left(rn+\left(r-1\right)\right)}\end{array}$$E102

for n = 0, 1, …, (N/r) – 1.

According to the method of least squares, the best fitting curve has the property that:

$$\begin{array}{l}J={\sigma }_0^2+\mathrm{...}+{\sigma }_{rn}^2+{\sigma }_1^2+\mathrm{...}+{\sigma }_{\left(rn+1\right)}^2+\mathrm{...}+{\sigma }_{\left(rn+\left(r-1\right)\right)}^2\\ ={\displaystyle \sum _{j_0=0}^{r-1}{\displaystyle \sum _{n=0}^{\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}{\left[d_{\left(rn+j_0\right)}-y_{\left(rn+j_0\right)}\right]}^2}}\\ ={\displaystyle \sum _{j_0=0}^{r-1}{\displaystyle \sum _{n=0}^{\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}{\sigma }_{\left(rn+j_0\right)}^2}}\text{=aminimum}\end{array}$$E103

The parallel implementation of the least squares for the linear case could be expressed as:

$$\begin{array}{l}{e}_{\left(n\right)}=d_{\left(n\right)}-\left(b+w{x}_{\left(n\right)}\right)=d_{\left(n\right)}-y_{\left(n\right)}\\ \text{=}\left(I_r\times \left[d_{\left(rn+j_0\right)}-y_{\left(rn+j_0\right)}\right]\right)\\ \text{=}\left(I_r\times \left[e_{\left(rn+j_0\right)}\right]\right)\end{array}$$E104

for j₀ = 1, …, r – 1 and in order to pick the line which best fits the data, we need a criterion to determine which linear estimator is the “best”. The sum of square errors (also called the mean square error (MSE)) is a widely utilized performance criterion.

$$\begin{array}{l}J=\frac{1}{2N}{\displaystyle \sum _{n=0}^{N-1}{e}_{\left(n\right)}^2}\\ =\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)$}\right.\right){\displaystyle \sum _{n=0}^{\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}{\left(\left(I_r\times \left[e_{\left(rn+j_0\right)}\right]\right)\right)}^2}}\end{array}$$E105

which yields after simplification to:

$$\begin{array}{l}J=\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\left(I_r\times \left[\left(\frac{1}{2\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)}\right){\displaystyle \sum _{n=0}^{\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}{e}_{\left(rn+j_0\right)}^2}\right]\right)}\\ \text{=}\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\left(I_r\times \left[J_{j_0}\right]\right)}\text{=}\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}{J}_{j_0}}\end{array}$$E106

where $J_{j_0}$ is the partial MSE applied on the subdivided data.

Our goal is to minimize J analytically, which according to Gauss can be done by taking its partial derivative with respect to the unknowns and equating the resulting equations to zero:

$$\{\begin{array}{c}\frac{\partial J}{\partial b}=0\\ \frac{\partial J}{\partial w}=0\end{array}$$E107

(109)

which yields:

$$\{\begin{array}{c}\frac{\partial J}{\partial b}=\frac{\partial \left(\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\left(I_r\times \left[J_{j_0}\right]\right)}\right)}{\partial b}=\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\frac{\partial \left(J_{j_0}\right)}{\partial b}=0}\\ \frac{\partial J}{\partial w}=\frac{\partial \left(\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\left(I_r\times \left[J_{j_0}\right]\right)}\right)}{\partial w}=\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right){\displaystyle \sum _{j_0=0}^{r-1}\frac{\partial \left(J_{j_0}\right)}{\partial w}=0}\end{array}$$E108

With the same reasoning as above the MSE could be obtained for multiple variables by:

$$\begin{array}{l}J=\frac{1}{2N}{{\displaystyle \sum _n\left(d_{\left(n\right)}-{\displaystyle \sum _{k=0}^p{w}_{\left(k\right)}{x}_{\left(n,k\right)}}\right)}}^2\\ \text{=}\frac{1}{r}\left({\displaystyle \sum _{j_0=0}^{r-1}\frac{1}{2\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)}{\displaystyle \sum _{n=0}^{\left(\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$r$}\right.\right)-1}{\left[d_{\left(rn+j_0\right)}-{\displaystyle \sum _{k=0}^p{w}_{\left(\left(rn+j_0\right),k\right)}{x}_{\left(\left(rn+j_0\right),k\right)}}\right]}^2}}\right)\\ =\frac{1}{r}{\displaystyle \sum _{j_0=0}^{r-1}{J}_{j_0}}\end{array}$$E109

for j₀ = 1, …, r – 1 and where $J_{j_0}$ is the partial MSE applied on the subdivided data.

The solution to the extreme (minimum) of this equation can be found in exactly the same way as before, that is, by taking the derivatives of $J_{j_0}$with respect to the unknowns (w_k), and equating the result to zero.

Instead of solving equations 110 and 111 analytically, a gradient adaptive system can be used which is done by estimating the derivative using the difference operator. This estimation is given by:

$${\nabla }_{w}J=\frac{\delta J}{\delta w}$$E110

where in this case the bias b is set to zero.

8.3. Search of the performance surface with steepest descent

The method of steepest descent (also known as the gradient method) is the simplest example of a gradient based method for minimizing a function of several variables [12]. In this section we will be elaborating the linear case.

Since the performance surface for the linear case implemented in parallel, are r paraboloids each of which has a single minimum, an alternate procedure to find the best value of the coefficient $$w_{j_{0\left(k\right)}}$$ is to search in parallel the performance surface instead of computing the best coefficient analytically by Eq. 110. The search for the minimum of a function can be done efficiently using a broad class of methods that use gradient information. The gradient has two main advantages for search.

• The gradient can be computed locally.

• The gradient always points in the direction of maximum change.

If the goal is to reach the minimum in each parallel segment, the search must be in the direction opposite to the gradient. So, the overall method of search can be stated in the following way:

Start the search with an arbitrary initial weight$$w_{j_{0\left(0\right)}}$$, where the iteration is denoted by the index in parenthesis (Fig 43). Then compute the gradient of the performance surface at$$w_{j_{0\left(0\right)}}$$, and modify the initial weight proportionally to the negative of the gradient at$$w_{j_{0\left(0\right)}}$$. This changes the operating point to$$w_{j_{0\left(1\right)}}$$. Then compute the gradient at the new position$$w_{j_{0\left(1\right)}}$$, and apply the same procedure again, i.e.

$$w_{j_{0_{\left(k+1\right)}}}=w_{j_{0_{\left(k\right)}}}-\eta \nabla J_{j_{0_{\left(k\right)}}}$$E111

where η is a small constant and $$\nabla J_{j_{0_{\left(k\right)}}}$$ denotes the gradient of the performance surface at the k^th iteration of j₀ parallel segment. η is used to maintain stability in the search by ensuring that the operating point does not move too far along the performance surface. This search procedure is called the steepest descent method (fig 43)

Figure 43.

The search using the gradient information [13].

If one traces the path of the weights from iteration to iteration, intuitively we see that if the constant η is small, eventually the best value for the coefficient w* will be found. Whenever w>w*, we decrease w, and whenever w<w*, we increase w.

8.4. The radix- r parallel LMS algorithm

Based on what was proposed in [2] by using the instantaneous value of the gradient as the estimator for the true quantity which means by dropping the summation in equation 108 and then taking the derivative with respect to w yields:

$$\begin{array}{l}\nabla J_{\left(k\right)}=\nabla \left(\frac{1}{r}\times \left(I_r\times \left[J_{j_0}{}_{\left(k\right)}\right]\right)\right)\\ =\frac{\partial }{\partial w}{J}_{\left(k\right)}=\frac{\partial }{\partial w}\left(\frac{1}{r}\times \left(I_r\times \left[J_{j_0}{}_{\left(k\right)}\right]\right)\right)\\ \approx \frac{\partial }{\partial w}\frac{1}{2Nr}{\displaystyle \sum e^2}=\frac{\partial }{\partial w}\frac{1}{2Nr}{\displaystyle \sum _{j_0=0}^{r-1}{\displaystyle \sum _{n=0}^{\left(\frac{N}{r}\right)-1}{\left(I_r\times \left[e_{\left(rn+j_0\right)}\right]\right)}^2}}\\ \approx \frac{1}{r}\frac{\partial }{\partial w}\left(I_r\times \left[e_{{\left(rn+j_0\right)}_{\left(k\right)}}{x}_{{\left(rn+j_0\right)}_{\left(k\right)}}\right]\right)\end{array}$$E112

What this equation tells us is that an instantaneous estimate of the gradient is simply the product of the input to the weight times the error at iteration k. This means that with one multiplication per weight the gradient can be estimated. This is the gradient estimate that leads to the famous Least Means Square (LMS) algorithm (Fig 44).

If the estimator of Eq.114 is substituted in Eq.113, the steepest descent equation becomes

$$I_r\times \left[w_{j_{0_{\left(k+1\right)}}}\right]=I_r\times \left(w_{j_{0_{\left(k\right)}}}+\left(\left(\frac{1}{r}\right)\left({\eta }_{j_0}{e}_{{\left(rn+j_0\right)}_{\left(k\right)}}{x}_{{\left(rn+j_0\right)}_{\left(k\right)}}\right)\right)\right)$$E113

for j₀ = 0,1, …, r -1.

This equation is the r Parallel LMS algorithm, which is used as predictive filter, is illustrated in Figure 45. The small constant η is called the step size or the learning rate.

Figure 44.

LMS Filter

Figure 45.

r Parallel LMS Algorithm Used in Predictive Filter.

8.5. Simulation results

The notion of a mathematical model is fundamental to sciences and engineering. In class of applications dealing with identification, an adaptive filter is used is used to provide a linear model that represents the best fit (in some sense) to an unknown signal. The LMS Algorithm which is widely used is an extremely simple and elegant algorithm that is able to minimize the external cost function by using local information available to the system parameters. Due to its computational burden and in order to speed up the process, this paper has presented an efficient way to compute the LMS algorithm in parallel where it follows from the simulation results that the stability of our models relies on the stability of our r parallel adaptive filters. It follows from figures 47 and 48 that the stability of r parallel LMS filters (in this case r = 2) has been achieved and the convergence performance of the overall model is illustrated in figure 49. The complexity of the proposed method will be reduced by a factor of r in comparison to the direct method illustrated in figure 46. Furthermore, the simulation result of the channel equalization is illustrated in figure 50 in which the blue curves represents our parallel implementation (2 LMS implemented in parallel) compared to the conventional method where the curve is in a red colour.

Figure 46.

Simulation Result of the original signal

Figure 47.

Simulation Result of the first partial LMS Algorithm.

Figure 48.

Simulation Result of the second partial LMS Algorithm.

Figure 49.

Simulation Result of the Overall System.

Figure 50.

Simulation Result of the channel equalization (blue curve two LMS implemented in Parallel red curve one LMS).