Homec4science
Diffusion epfl (master) Browse

R2653/writing/sp4comm.multipub/80-stochastic/20-rn-adaptive.tex
b5b12b369fe4master

Last Change
Raw File

20-rn-adaptive.tex
View Options

\section{Adaptive Signal Processing}
\label{sec:rn:adaptive}
The question at the heart of \textit{adaptive signal processing} is whether we can design a digital system that, when processing an unknown signal, can automatically adapt to the input's properties in order to perform a certain task. When the processing system in question is a discrete-time filter, the problem can be simplified to finding a way to adapt the filter coefficients as a function of the difference between the filter's output and a desired ideal output.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t!]
\centering
\begin{dspBlocks}{1.3}{1.2}
$x[n]~~~$ & \BDfilter{$\mathcal{H}$} & \BDsplit & $~~~\hat{d}[n] \approx d[n]$ \\
$d[n]~~~$ & \BDsub & & %
\psset{arrows=->}
\ncline{1,1}{1,2}\ncline{1,2}{1,4}
\ncline{-}{1,3}{2,3}\ncline{2,3}{2,2}
\ncline{2,1}{2,2}
\ncline[arrowscale=1.2,arrowinset=0,ArrowFill=false]{2,2}{1,2}\nbput{$e[n]$}
\end{dspBlocks}
\caption{Standard setup for adaptive signal processing.}\label{fig:rn:adaptive}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
For example, communication systems almost always have to operate over channels with unknown characteristics; if a transmitter sends a signal $\mathbf{d}$, the channel will modify $\mathbf{d}$ in some way and yield a signal $\mathbf x$. The first task for a receiver, therefore, is to compensate as much as possible for the unwanted effects of the channel and to produce a good estimate $\mathbf{\hat{d} = \mathcal{H}\{x\}}$ of the transmitted signal before trying to decode the incoming data; the operation performed by $\mathcal{H}$ is called \textit{channel equalization}. Let's assume that, at least initially, the original signal $\mathbf d$ is known exactly at the receiver; this is usually achieved by designing the transmitter so that it sends a prearranged training (or \textit{bootstrapping}) sequence at the beginning of a communication session. In this case the receiver can compute an exact \textit{error signal}
\[
\mathbf{e = d - \hat{d} = d - \mathcal{H}\{x\}}
\]
and use it to fine-tune the coefficients for $\mathcal{H}$ so that $\mathbf e$ becomes small. This example captures the fundamental setup of most adaptive signal processing scenarios where, as in Figure~\ref{fig:rn:adaptive}, an error signal is used to drive the filter adaptation.
Given an adaptive system and its resulting error signal, in general there is little hope to reduce the error to zero; a more reasonable strategy is trying to \textit{minimize the energy} of the error signal or, in other words, find the system that minimizes the mean squared error
\begin{equation} \label{eq:rn:msefunc}
J = \expt{|e[n]|^2}.
\end{equation}
MSE minimization is certainly satisfying intuitively but, equally importantly, using a quadratic cost function is very convenient because of its very desirable mathematical properties:
\begin{itemize}
\item the minimum exists and it is unique, since the error surface is a paraboloid
\item the expression for error is always differentiable, which leads to simple iterative algorithms to find the minimum
\item as we will see, the signal that minimizes the MSE turns out to be \textit{orthogonal} to the error signal
\item once again, when considering the stochastic nature of the problem, we will only need the first two moments of the processes involved in the estimation.
\end{itemize}
While we have made no assumptions yet on the nature of the system $\mathcal{H}$, in the remainder of this chapter we will restrict ourselves to the case in which $\mathcal{H}$ is a \textit{causal FIR filter} with transfer function $H(z)$; this setup, in fact, covers a large majority of practical applications, particularly in communication system design. With this, our task is to find an algorithm to determine the coefficients of an $N$-tap FIR filter that, given an input signal $\mathbf{x}$ and a desired output $\mathbf{d}$, minimizes the expected squared value of
\begin{equation}\label{eq:rn:errsig}
e[n] = d[n] - \sum_{k=0}^{N-1}h[k]x[n-k].
\end{equation}
We will now analyze in detail two common situations that emerge in adaptive filtering:
\begin{itemize}
\item if all signals can be safely assumed to be stationary then we can solve the problem analytically in one shot and find the optimal filter using classic least squares minimization techniques;
\item if the signals are slowly time-varying, the filter will have to ``track'' the input over time; in this case we will solve~(\ref{eq:rn:msefuncfir}) by an iterative gradient descent method.
\end{itemize}
\subsection{One-Shot Solution: Optimal Least Squares}
\label{sec:rn:leastsquares}
Assume $\mathbf{x}$ and $\mathbf{d}$ are two zero-mean, jointly WSS processes.\footnote{
%
This derivation can also be applied to deterministic signals with expectations replaced by time averages.}
%
Let
\begin{equation}\label{eq:rn:hvec}
\mathbf{h} = \begin{bmatrix}
h[0] & h[1] & h[2] & \ldots & h[N-1]
\end{bmatrix}^T
\end{equation}
be a vector containing the $N$ taps of the FIR filter used to derive an estimate of $\mathbf{d}$ based on the most recent $N$ samples of $\mathbf{x}$. The optimal filter that minimizes the mean squared error can be found by taking partial derivatives of the MSE with respect to each filter coefficient; using the expression for the error signal in~(\ref{eq:rn:errsig}):
\begin{align}
\frac{1}{2}\,\frac{\partial \expt{e^2[n]}}{\partial h[i]} &= \expt{e[n] \,\displaystyle \frac{\partial e[n]}{\partial h[i]}} \label{eq:rn:partialMSE} \\
&= -\expt{e[n] \, x[n-i]} \label{eq:rn:partialMSE2} \\
&= \expt{\displaystyle \sum_{k=0}^{N-1}h[k]x[n-k]x[n-i]} - \expt{d[n]x[n-i]} \nonumber \\
&= \sum_{k=0}^{N-1}h[k]r_x[i-k] - r_{xd}[i]
\end{align}
where we have exploited the fact that $\expt{x[n-k]x[n-i]}$ depends only on the time lag $(n-k) - (n-i) = i-k$ and that, similarly, $\expt{d[n]x[n-i]} = r_{dx}[-i] = r_{xd}[i]$. In order to find the global minimum, we set each partial derivative to zero to obtain the following system of equations that only uses the values of the auto- and cross-correlation:
\[
\sum_{k=0}^{N-1}h[k]r_x[i-k] = r_{xd}[i];
\]
We can write the system in matrix form as
\begin{equation}\label{eq:rn:lsmat}
\mathbf{R\,h = g}
\end{equation}
where $\mathbf{R}$ is the $N\times N$ symmetric Toeplitz matrix
\begin{equation}\label{eq:rn:rmat}
\mathbf{R} = \begin{bmatrix}
r_x[0] & r_x[1] & r_x[2] & \ldots & r_x[N-1] \\
r_x[1] & r_x[0] & r_x[1] & \ldots & r_x[N-2] \\
r_x[2] & r_x[1] & r_x[0] & \ldots & r_x[N-3] \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
r_x[N-1] & r_x[N-2]& \ldots & \ldots & r_x[0]
\end{bmatrix}
\end{equation}
and
\begin{equation}\label{eq:rn:gvec}
\mathbf{g} = \begin{bmatrix}
r_{xd}[0] & r_{xd}[1] & r_{xd}[2] & \ldots & r_{xd}[N-1]
\end{bmatrix}^T.
\end{equation}
From a computational point of view, the particular structure of the autocorrelation matrix allows for a very efficient algorithmic approach to solve the system of equations. While in the general case the solution would require $O(N^3)$ operations, \textit{Levinson's recursion} exploits the Toeplitz nature of $\mathbf{R}$ to build a solution incrementally at a computational cost of just $O(N^2)$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% steepestdescent{a}{b}{x0}{y0}{step}{iterations}{view}{algo}:
% a, b: ellipsoid axes
% x0, y0, step, iterations: gradient descent initial position, learning step and iterations
% view: 3D or 2D (top view)
% algo: LS (least squares) or LMS (stochastic gradient descent)
%
\def\steepestdescent#1#2#3#4#5#6#7#8{%
% cartesian elliptic paraboloid z = f(x, y); expects to find [x y] on the stack
% add an offset to lift the paraboloid off the origin (it's an error surface)
\def\ellparC{#2 div dup mul exch #1 div dup mul add 0.2 add}
% parametric form (x, y, z) = f(theta, rho); expects [theta rho] on the stack
\def\ellparP{/r exch def dup cos #1 r mul mul exch sin #2 r mul mul r r mul 0.2 add}
% x and y updates in steepest descent: x -= stepsize * dz/dx
\ifthenelse{\equal{#8}{LMS}}{\def\rndp{rand 2147483647 div 0.8 sub 0.3 mul}}{\def\rndp{0}}
\def\gradx{ pop dup -2 #5 mul mul #1 dup mul div add \rndp\space add}
\def\grady{exch pop dup -2 #5 mul mul #2 dup mul div add \rndp\space add}
%
% center plot in psbox according to view
\ifthenelse{\equal{#7}{2D}}{\def\la{-2.6}\def\lb{3}}{\def\la{-1}\def\lb{5}}
\begin{pspicture}(-3,\la)(3,\lb)
\ifthenelse{\equal{#7}{2D}}{%
\psset{Beta=90}\psset{Alpha=180}}{
\psset{Beta=30}\psset{Alpha=150}}
\pstThreeDCoor[linecolor=darkgray,%
xMin=-2.5,xMax=2.5,nameX={$h_1$},
yMin=-2.5,yMax=2.5,nameY={$h_0$},
zMin=0,zMax=3,nameZ={}]
% plot the level sets
\parametricplotThreeD[xPlotpoints=180,yPlotpoints=20,linecolor=red,linewidth=0.5pt](0,360)(0, 2){%
t u \ellparP}
\parametricplotThreeD[xPlotpoints=180,yPlotpoints=17,linecolor=red,linewidth=0.5pt](0, 2)(0,360){%
u t \ellparP}
% gradient descent
\pstVerb{/xx #3 def /yy #4 def /xn xx yy \gradx\space def /yn xx yy \grady\space def}
\multido{\nA=0+1}{#6}{%
\pstThreeDLine[linecolor=blue,arrows=->](xx, yy, xx yy \ellparC)(xn, yn, xn yn \ellparC\space)%
\pstVerb{/xx xn def /yy yn def /xn xx yy \gradx\space def /yn xx yy \grady\space def}}
\end{pspicture}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\psset{unit=10mm}
\begin{tabular}{cc}
{\steepestdescent{1}{1}{-1.3}{1.3}{0.1}{0}{3D}{LS}} &
{\steepestdescent{1}{1}{-1.3}{1.3}{0.1}{0}{2D}{LS}}
\end{tabular}
\caption{Perspective and top view of the error surface in~(\ref{eq:rn:paraboloid}) as a function of the filter coefficients for $N=2$ when the input signal is white.}\label{fig:rn:revparaboloid}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{Error surface.} In the general case, the MSE is a high-dimensional quadratic function of the $N$ filter coefficients $\mathbf{h}$. If we restrict ourselves to the case $N=2$ we can actually plot the error surface and gain some intuition on its properties and on the properties of the algorithms that we will use to find its minimum. Indeed, for $N=2$, let $\mathbf{h} = \begin{bmatrix}h_0 & h_1 \end{bmatrix}^T$ be the two-tap optimal filter; we have
\begin{align}
J &= \expt{\left(d[n] - (h_0 x[n] + h_1 x[n-1])\right)^2} \nonumber \\
&= \sigma_d^2 + r_x[0]h_0^2 + r_x[0]h_1^2 + 2r_x[1]h_0h_1 - 2r_{xd}[0]h_0 - 2r_{xd}[1]h_1 \nonumber \\
&= \sigma_d^2 + \begin{bmatrix}h_0 & h_1\end{bmatrix} \begin{bmatrix}r_x[0] & r_x[1] \\ r_x[1] & r_x[0]\end{bmatrix} \begin{bmatrix}h_0 \\ h_1\end{bmatrix} - 2\begin{bmatrix}h_0 & h_1\end{bmatrix}\begin{bmatrix}r_{xd}[0] \\ r_{xd}[1]\end{bmatrix}.
\end{align}
The above expression defines an elliptic paraboloid surface. The last term is nonzero only if input and desired signal are correlated and it simply translates the surface to a different origin; it can be therefore ignored without loss of generality. In compact form we thus have
\begin{equation} \label{eq:rn:paraboloid}
J = \sigma_d^2 + \mathbf{h}^T\mathbf{R\,h}
\end{equation}
where $\mathbf{R}$ is the $2\times 2$ autocorrelation matrix for the input signal. If the latter is uncorrelated (e.g. if $x[n]$ is a white noise sequence) then $\mathbf{R}$ is a multiple of the identity matrix:
\[
r_x[n] = \sigma_x^2\delta[n] \quad\Rightarrow\quad \mathbf{R} = \begin{bmatrix}\sigma_x^2 & 0 \\ 0 & \sigma_x^2\end{bmatrix}
\]
and thus the error surface is a paraboloid of revolution whose cross sections are circles, as in Figure~\ref{fig:rn:revparaboloid}. Otherwise, by appropriately rotating the $h_0, h_1$ coordinate system, we can rewrite~(\ref{eq:rn:paraboloid}) in canonical form as
\begin{equation} \label{eq:rn:paraboloid}
J = \sigma_d^2 + \mathbf{h}^T \begin{bmatrix}\lambda_0 & 0 \\ 0 & \lambda_1 \end{bmatrix} \mathbf{h}
\end{equation}
where $\lambda_{0,1}$ are the eigenvalues of the autocorrelation matrix.\footnote{
%
Since $\mathbf{R}$ is real and symmetric, it can always be diagonalized with real eigenvalues.}.
%
This defines an elliptic paraboloid surface in which the axes of the elliptical cross sections are proportional to $1/\sqrt{\lambda_{0,1}}$, as in Figure~\ref{fig:rn:ellparaboloid}. The consequences of having a correlated input (and therefore a ``squished'' error surface) will be apparent in the next section where we will minimize $J$ via iterative methods.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\psset{unit=10mm}
\begin{tabular}{cc}
{\steepestdescent{1}{0.4}{-1.3}{1.3}{0.1}{0}{3D}{LS}} &
{\steepestdescent{1}{0.4}{-1.3}{1.3}{0.1}{0}{2D}{LS}}
\end{tabular}
\caption{Error surface when the input signal is correlated.}\label{fig:rn:ellparaboloid}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{Orthogonality principle.} The optimal least squares filter corresponds to a point of zero gradient for the MSE; using~(\ref{eq:rn:partialMSE2}) we can thus write
\begin{equation} \label{eq:rn:ortherr}
\expt{e[n] \, x[n-i]}, \qquad i = 0, \ldots, N-1,
\end{equation}
which means that the error signal is orthogonal (in a statistical sense) to the input signal for lags ranging from $0$ to $N-1$. According to the interpretation discussed in Section~\ref{sec:rn:orth}, this means that the error ``shares no information'' with the input samples used to produce the estimate $\mathbf{\hat{d} = h \ast x}$ or, in other words, that all the usable information in the input has been exploited by the optimal filter. As a simple consequence of linearity, note that the error is also orthogonal to the output estimate:
\begin{equation} \label{eq:rn:orthest}
\expt{e[n] \, \hat{d}[n]} = 0.
\end{equation}
\itempar{Segmentation.} The one-shot optimal least squares solution assumes that all signals are stationary so that the auto- and cross-correlations can be computed experimentally over the entire available data set; using all the data obviously produces a more accurate estimate of the correlation and a more precise localization of the minimum of the cost function. Many interesting real-world signals, however, fall in between stationarity and non-stationarity, in the sense that their underlying generative process switches between several distinct modes of operation. These signals are called \textit{piecewise stationary} in the sense that they can be segmented into short pieces so that each piece exhibits enough ``consistency'' as to be considered stationary. Determining the optimal split points can be extremely complicated, however, so the usual strategy is to use a fixed-size window, short enough to ensure consistency and long enough to allow for a reasonably accurate estimation of the correlations. Typical examples of piecewise-stationary signals are speech and music; for those, adaptive algorithms usually operate on segments whose length is between 20 and 30 milliseconds.
\subsection{Iterative Solution: the LMS Filter}
Most function minimization problems can be approached iteratively using a technique called \textit{steepest descent}; the idea is to start with an initial estimate for the location of the minimum and then progressively refine the estimate by computing the local gradient of the function and moving in the opposite direction. Intuitively (at least in three dimensions), gradient descent is akin to letting a ball roll down a terrain with peaks and dips: the ball will follow the slope along a path of maximal steepness until it reaches the bottom of a valley and stop. Formally, given a differentiable multivariate function $f(\mathbf{h})$ and an initial guess $\mathbf{h}_0$, the iterative procedure computes a series of estimates for the function's minimum as
\begin{equation} \label{eq:rn:steepest}
\mathbf{h}_{n+1} = \mathbf{h}_{n} - \alpha_n \, \nabla f(\mathbf{h}_{n})
\end{equation}
where
\[
\nabla f(\mathbf{x}) = \begin{bmatrix}
\displaystyle \frac{\partial f(\mathbf{h})}{\partial h_0} &
\displaystyle \frac{\partial f(\mathbf{h})}{\partial h_1} & \ldots &
\displaystyle \frac{\partial f(\mathbf{h})}{\partial h_{N-1}}
\end{bmatrix}^T
\]
and where $\alpha_n$ is a sequence of small-valued scalar ``learning factors'' that ease the convergence process (in the ball-rolling analogy, think of the learning factor as of some sort of ``brake'' that slows down the ball so that it doesn't overshoot the valley). Under fairly mild assumptions, the sequence of estimates in~(\ref{eq:rn:steepest}) converges reasonably quickly to a \textit{local} minimum of the function.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\psset{unit=10mm}
\begin{tabular}{cc}
{\steepestdescent{1}{1}{-1.3}{1.3}{0.1}{15}{3D}{LS}} &
{\steepestdescent{1}{1}{-1.3}{1.3}{0.1}{15}{2D}{LS}}
\end{tabular}
\caption{Perspective and top view of the steepest descent minimization when the input signal is white.}\label{fig:rn:sdrev}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
The ``local'' qualification is important since a major limitation of the steepest descent method is the fact that, depending on our initial choice for $\mathbf{x}_0$, we may not reach the function's global minimum. The MSE in~(\ref{eq:rn:msefunc}), however, is a quadratic function and therefore, if it possesses a minimum, the minimum will be global. If we apply the steepest descent method, we need to compute the gradient of the MSE as a function of the filter $\mathbf{h}$; we have already done so in~(\ref{eq:rn:partialMSE}) and we found that
\begin{align}
\nabla J(\mathbf{h}) &= \begin{bmatrix}
\displaystyle \frac{\partial \expt{e^2[n]}}{\partial h[0]} &
\displaystyle \frac{\partial \expt{e^2[n]}}{\partial h[1]} & \ldots &
\displaystyle \frac{\partial \expt{e^2[n]}}{\partial h[N-1]}
\end{bmatrix}^T \\
& = 2(\mathbf{R\,h - g}).
\end{align}
where $\mathbf{R}$ and $\mathbf{g}$ are defined in~(\ref{eq:rn:rmat}) and~(\ref{eq:rn:gvec}). We could use this value for the gradient directly in~(\ref{eq:rn:steepest}) and we would arrive, iteratively, at the same values for $\mathbf{h}$ as in the direct least squares solution obtained in the previous section. The procedure is illustrated graphically in Figure~\ref{fig:rn:sdrev} for a white input signal; note that, since the error surface is circularly symmetric, the steepest descent path to the minimum is independent of the initial starting point.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[th!]
\center
\psset{unit=10mm}
\begin{tabular}{cc}
{\steepestdescent{1}{0.4}{-2}{0}{0.1}{12}{2D}{LS}} &
{\steepestdescent{1}{0.4}{-1.1}{0.5}{0.15}{15}{2D}{LS}} \\
(a) & (b)
\end{tabular}
\steepestdescent{1}{0.4}{-1.1}{0.5}{0.17}{10}{2D}{LS} \\
(c)
\caption{Steepest descent minimization when the input signal is correlated: (a) and (b) show the dependence of the minimization path on the initial guess; (c) shows the effects of too large a learning factor, resulting in divergence.}\label{fig:rn:sdell}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
For a correlated input signal, on the other hand, the initial choice plays a much more delicate role, as shown in Figure~\ref{fig:rn:sdell}. Since the trajectory always moves in a direction orthogonal to the slope, a point halfway between the origin and the major semiaxis will determine an iteration step in a direction that is almost parallel to the opposite axis. In this case the convergence will be slower or, if the learning factor is too large, the algorithm may fail to converge.
For these reasons, and since the Levinson-Durbin algorithm is very efficient and numerically stable, there is no real advantage in choosing the iterative route... unless the statistical properties of the signals are slowly time-varying. In this case, we can modify the steepest descent algorithm so that it keeps adapting in real time to the changing properties of the input. For a nonstationary signal, the empirical auto- and cross-correlation are no longer time-invariant; when working with time averages we could try to recompute them based on a sliding data window, but it would be very expensive computationally. Instead, a more pragmatic approach is to do away with expectations and replace the average squared error by its \textit{instantaneous} value:
\begin{equation}
\expt{e^2[n]} \, \leftarrow \, e^2[n];
\end{equation}
this is, for lack of better words, a drastic simplification -- but, remarkably, one that works extremely well in practice. Note that, by removing the expectation operator, the gradient is now a random variable dependent on the index $n$; for this reason, the iterative solution that we will present belongs to the class of \textit{stochastic gradient descent} methods, and it is known with the name of Least Mean Squares, or LMS.
Remember that the error signal, for a FIR adaptive filter, was defined in~(\ref{eq:rn:errsig}) as
\[
e[n] = d[n] - \sum_{k=0}^{N-1}h[k]x[n-k]
\]
so that the $i$-th element of the instantaneous stochastic gradient is simply
\begin{equation}
\frac{\partial e^2[n]}{\partial h[i]} = -2e[n]\,x[n-i].
\end{equation}
We can express the full gradient as
\begin{equation}\label{eq:rn:stochgrad}
\nabla J_n = -2e[n]\, \mathbf{x}_n,
\end{equation}
where we have defined
\begin{equation}
\mathbf{x}_n = \begin{bmatrix}
x[n] & x[n-1] & x[n-2] & \ldots & x[n - N + 1]
\end{bmatrix}^T.
\end{equation}
The above expression can now be used with the steepest descent iteration in~(\ref{eq:rn:steepest}): the resulting LMS update algorithm starts with an initial filter estimate $\mathbf{h}_0$ and, at each step $n$, refine the filter coefficients with the following two-part recursion:\footnote{
%
The factor of 2 in~(\ref{eq:rn:stochgrad}) has been ``absorbed'' in the learning factor $\alpha_n$ for conciseness.}
%
\begin{align}
e[n] &= d[n] - \mathbf{h}_n^T\mathbf{x}_n \label{eq:rn:lmserr} \\
\mathbf{h}_{n+1} &= \mathbf{h}_n + \alpha_n\, e[n] \, \mathbf{x}_n. \label{eq:rn:lmspred}
\end{align}
Note that, contrary to the standard steepest descent minimization in~(\ref{eq:rn:steepest}), where the procedure is deterministic and we expect it to converge to a final estimate for $\mathbf{h}$, the LMS update is a stochastic update and it is meant to run ``forever'' in order to track over time the nonstationary characteristic of the signals. Note also that the algorithm is extremely lightweight computationally speaking, since it only requires on the order of $N$ operations per sample.
Unfortunately, in spite of its algorithmic simplicity, a theoretical analysis of the LMS filter's performance is extremely hard. Convergence is difficult to prove except for the simplest of cases and, even then, the proofs focus on stationary signals and require rather unrealistic hypotheses on their statistical properties. We will therefore skip any attempt at a formal analysis and simply provide some guidelines for practical use cases.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\psset{unit=10mm}
\begin{tabular}{cc}
{\steepestdescent{1}{1}{-1.3}{1.3}{0.1}{25}{2D}{LMS}} &
{\steepestdescent{1}{0.4}{-2}{0}{0.1}{25}{2D}{LMS}} \\
(a) & (b)
\end{tabular}
\caption{LMS error signal for (a) white and (b) correlated inputs.}\label{fig:rn:sdlms}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Intuitively, if we go back to the ball-rolling interpretation of steepest descent, the removal of the expectation operator transforms the smooth, averaged surface of the MSE function into a ragged, spiky terrain; these asperities may momentarily set the descent in the wrong direction and interfere with overall convergence. This is illustrated in Figure~\ref{fig:rn:sdlms}; note how the stochastic nature of the error determines a ``meandering'' trajectory for the gradient descent, which keeps hovering around the minimum. Because of this, the choice of the learning factor becomes crucial; by choosing $\alpha_n$ wisely, we can hope to counteract the effects of using a noisy gradient. One rule of thumb that is often quoted in the literature is to set
\begin{equation}
\alpha_n < 2/\max_i\{\lambda_i\}
\end{equation}
where $\lambda_i$ are the eigenvalues of the autocorrelation matrix. In general this bound is difficult to determine since the autocorrelation matrix is unknown and, for time-varying signals, changing over time; in practice, the learning factor is determined experimentally and kept small to prevent instability at the price of a lower rate of convergence.
\section{Practical Applications}
\subsection{Linear Prediction and Speech Compression}
\label{sec:rn:linpred}
Most physical systems that can be modeled as linear filters are fully characterized by their intrinsic \textit{resonances}, that is, by the frequencies at which the transfer of energy from input to the output is particularly effective. A prototypical example is given by musical instruments, that are designed and built to resonate mechanically over a specific range of pitches; similarly, in communication systems, any transmission medium exhibits a preferred frequency band where the attenuation is minimal. In an idealized model representing the physical system, resonances correspond to the poles of the transfer function and often an all-pole model is sufficient to capture the main characteristics of the system.
Consider a system $\mathcal{H}$ whose unknown transfer function is of the form
\begin{equation}
H(z) = \frac{1}{A(z)} = \frac{1}{1 - a_1 z^{-1} - a_2 z^{-2} - \ldots - a_N z^{-N}};
\end{equation}
such a characteristic, as we have seen in Section~\ref{sec:rn:arma}, is called an \textit{autoregressive} (AR) model. Assume the input to the system is an unknown WSS random excitation $\mathbf{w}$ and that we can observe the resulting output random signal $\mathbf{x = h \ast w}$ but not the excitation signal itself. Can we still determine $H(z)$ from our observations? Let's begin by writing out explicitly the CCDE governing the system, which is
%\begin{equation} \label{eq:rn:arzdoman}
% X(z) = H(z)E(z) = E(z)/A(z)
%\end{equation}
%we can write
\begin{equation}
x[n] = \sum_{k=1}^{N} a_k x[n-k] + w[n]. \label{eq:rn:arccde}
\end{equation}
The equation shows that each output sample is a linear combination of past output samples, plus the unknown excitation term; a good estimate $\hat{a}_k$ of the coefficients in the linear combination would allow us to produce a good \textit{linear prediction} of the current observation $x[n]$. Of course, since the excitation is unknown, the prediction will not necessarily be identical to the observed value but a good strategy to find the set of coefficients is to minimize the prediction error
\begin{equation}
e[n] = x[n] - \sum_{k=1}^{N} \hat{a}_k x[n-k] = w[n]. \label{eq:rn:linpred}
\end{equation}
We see that the error turns out to be equal to the unknown excitation; this makes sense, since the excitation is the part of the output that can not be predicted from past observed samples. With this formulation we can find the set of $\hat{a}_k$'s that minimizes the energy of $\mathbf{e}$ using the Optimal Least Squares approach of Section~\ref{sec:rn:leastsquares}. For this set of optimal coefficients the orthogonality principle in~(\ref{eq:rn:ortherr}) ensures that the excitation $\mathbf{w}$ is orthogonal to the past $N-1$ output values; this means that the linear predictor has extracted all ``structured'' information from the observations and that, indeed, the excitation is uncorrelated from the response of the system.
To determine the optimal coefficients, once again we set to zero the partial derivatives with respect to the $a_k$'s of the expected squared error to obtain
\begin{equation}
\mathbf{R\,\hat{a} = r}
\end{equation}
or, explicitly,
\begin{equation}\label{eq:rn:linpredmat}
\begin{bmatrix}
r_x[0] & r_x[1] & \ldots & r_x[N-1] \\
r_x[1] & r_x[0] & \ldots & r_x[N-2] \\
\vdots & \vdots & \ddots & \vdots \\
r_x[N-1] & r_x[N-2]& \ldots & r_x[0]
\end{bmatrix}
\begin{bmatrix}
\hat{a}_1 \\ \hat{a}_2 \\ \vdots \\ \hat{a}_N
\end{bmatrix} =
\begin{bmatrix}
r_x[1] \\ r_x[2] \\ \vdots \\ r_x[N]
\end{bmatrix}.
\end{equation}
Incidentally, the particular structure of the above system of equations makes it possible to further simplify Levinson's recursion into what is known as \textit{Durbin's algorithm}, making linear prediction remarkably efficient from the computational point of view.
\itempar{LPC.} One of the most successful applications of linear prediction is speech compression for the efficient use of bandwidth in mobile communications. At the heart of the many variants of Linear Prediction Coding (LPC) is a source-filter model for the human vocal apparatus. The source is either a pitched sound, produced by the vibration of the vocal cords when pronouncing vowels, or a noise-like excitation produced by a flow of air from the lungs in the case of consonants; the filter is the time-varying acoustic resonator created by the changing configuration of mouth and lips. Source and filter are independent; as an example, remark that, in order to whisper rather than talk, the pitched excitation normally used for vowels is replaced by the noise-like excitation normally reserved for consonants. A schematic view of the source-filter model is shown in Figure~\ref{fig:rn:lpc}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[tb]
\centering
\includegraphics[height=80mm]{\localpath{figs/lpc.eps}}
\caption{Schematic representation of the source-filter model for the vocal production mechanism. } \label{fig:rn:lpc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The transfer function associated to the voice apparatus can be efficiently described by an all-pole model; the resonances associated to the poles (i.e., the peaks in the magnitude response) are commonly known as the voice \textit{formants} and their shape and relative position determine the acoustic quality of the uttered phonemes. The LPC compression scheme estimates the coefficients of the transfer function via linear prediction; the resulting inverse filter $A(z)$ is used to obtain the prediction error (that is, the excitation) simply by filtering the signal with it. The excitation is then analyzed and it is either classified as noise or further processed to determine its pitch in the case of a voiced sound. Since both the source and the resonator vary quickly over time, LPC compression proceeds by first splitting the voice data into small chunks (about 20ms long) over which the signal can be assumed to be stationary; for each segment, the encoding produces $N$ LPC coefficients $\hat{a}_k$ and a descriptor for the excitation. At the decoder, the source is synthesized from the descriptor and filtered with the transmitted $A(z)$ to produce the reconstructed speech signal. Significant data reduction is achieved by appropriately quantizing the coefficients and the description of the residual, for an average bitrate of about 4Kbit/sec for basic LPC of speech sampled at 8~KHz; at 8~bit/sample, by comparison, the original signal would require 48Kbit/sec. %The perceptual quality of the encoded speech can be improved by modeling the excitation as a linear combination of suitably pre-defined ``snippets'' at the price of a more computationally demanding encoding process; the resulting algorithm is known as Codebook Excited Linear Prediction (CELP) and is used in the GSM standard.
The choice of model order $N$ is crucial to the quality of the reconstructed speech. There should be enough poles to capture the formants of the vocal tract but their number should not be so large as to start modeling the peaks of the excitation in the case of voiced segments. As an example, Figure~\ref{fig:rn:lpcspec} shows the magnitude spectrum of a 20ms voiced speech segment together with the magnitude response of an all-pole filter $1/A(z)$ whose coefficients have been computed using Durbin's algorithm. Clearly, the order-10 filter is not enough to model the overall spectral envelope while the order-30 filter starts to lock on to the first partials of the periodic source. A good compromise is using about 20~LPC coefficients, which is the common choice in commercial voice encoders.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\begin{dspPlot}[xtype=freq,xticks=none,yticks=none]{0, 0.7}{0, 30}
\dspCustomTicks{0 $0$ 0.7 $0.7\pi$}
\dspFuncData[linewidth=0.5pt,linecolor=ColorDF]{0.000 1.032 0.003 3.946 0.006 1.522 0.009 3.315 0.013 9.486 0.016 28.335 0.019 7.367 0.022 5.142 0.025 7.885 0.028 17.662 0.031 24.412 0.034 8.163 0.037 4.155 0.041 3.531 0.044 11.966 0.047 9.724 0.050 6.683 0.053 5.400 0.056 6.014 0.059 24.546 0.062 4.555 0.066 0.697 0.069 4.168 0.072 11.327 0.075 16.578 0.078 4.988 0.081 2.965 0.084 4.597 0.087 6.698 0.091 6.949 0.094 1.626 0.097 1.480 0.100 3.226 0.103 5.024 0.106 2.831 0.109 2.091 0.113 0.326 0.116 2.605 0.119 5.985 0.122 1.631 0.125 1.210 0.128 0.516 0.131 1.910 0.134 3.508 0.138 0.689 0.141 0.984 0.144 1.490 0.147 3.292 0.150 3.718 0.153 1.573 0.156 0.915 0.159 1.669 0.163 2.744 0.166 1.799 0.169 2.157 0.172 1.576 0.175 4.013 0.178 6.629 0.181 1.634 0.184 2.229 0.188 2.203 0.191 4.141 0.194 6.570 0.197 3.931 0.200 3.150 0.203 1.583 0.206 7.722 0.209 12.751 0.212 9.065 0.216 2.051 0.219 1.209 0.222 4.332 0.225 2.593 0.228 5.901 0.231 6.797 0.234 0.700 0.237 2.986 0.241 6.158 0.244 5.257 0.247 2.083 0.250 2.485 0.253 2.002 0.256 1.170 0.259 1.893 0.263 1.858 0.266 3.603 0.269 4.275 0.272 4.068 0.275 1.357 0.278 1.043 0.281 1.866 0.284 1.218 0.287 2.384 0.291 2.011 0.294 1.664 0.297 1.731 0.300 3.189 0.303 4.059 0.306 3.623 0.309 3.590 0.312 2.637 0.316 2.839 0.319 3.392 0.322 3.139 0.325 4.262 0.328 4.209 0.331 6.186 0.334 4.408 0.338 1.870 0.341 3.318 0.344 1.764 0.347 2.931 0.350 1.830 0.353 0.914 0.356 1.243 0.359 0.891 0.362 1.771 0.366 1.737 0.369 1.381 0.372 0.673 0.375 0.820 0.378 1.058 0.381 1.693 0.384 1.248 0.388 0.307 0.391 1.554 0.394 1.075 0.397 1.398 0.400 1.256 0.403 0.054 0.406 1.451 0.409 2.293 0.412 0.960 0.416 1.048 0.419 0.875 0.422 1.148 0.425 1.857 0.428 0.902 0.431 1.045 0.434 0.663 0.438 2.455 0.441 4.253 0.444 2.403 0.447 1.411 0.450 1.962 0.453 2.467 0.456 4.276 0.459 2.663 0.463 0.642 0.466 2.873 0.469 4.175 0.472 4.425 0.475 2.058 0.478 0.394 0.481 3.268 0.484 3.816 0.487 0.949 0.491 1.008 0.494 0.411 0.497 1.317 0.500 2.771 0.503 0.557 0.506 0.668 0.509 1.220 0.512 1.730 0.516 1.706 0.519 0.297 0.522 0.179 0.525 0.317 0.528 0.752 0.531 1.101 0.534 0.575 0.537 0.093 0.541 0.755 0.544 0.973 0.547 0.631 0.550 0.347 0.553 0.474 0.556 0.146 0.559 0.871 0.562 0.410 0.566 0.165 0.569 0.912 0.572 0.760 0.575 1.178 0.578 0.139 0.581 0.637 0.584 0.297 0.588 0.394 0.591 0.289 0.594 0.364 0.597 0.370 0.600 0.270 0.603 0.490 0.606 0.704 0.609 0.340 0.613 0.234 0.616 0.162 0.619 0.317 0.622 0.237 0.625 0.213 0.628 0.158 0.631 0.282 0.634 0.645 0.637 0.064 0.641 0.240 0.644 0.156 0.647 0.187 0.650 0.413 0.653 0.167 0.656 0.249 0.659 0.172 0.662 0.283 0.666 0.004 0.669 0.214 0.672 0.179 0.675 0.156 0.678 0.271 0.681 0.144 0.684 0.185 0.688 0.129 0.691 0.221 0.694 0.096 0.697 0.171 0.700 0.197}
\dspFunc[linecolor=ColorTwo]{x \dspTFM{1}{1 -2.243 2.556 -1.809 0.429 0.545 -0.66 0.183 0.099 -0.113 0.066}}
\dspFunc[linecolor=ColorThree]{x \dspTFM{1}{1 -2.179 2.414 -1.679 0.363 0.562 -0.705 0.196 0.187 -0.288 0.235 -0.058 -0.081 0.095 0.124 -0.136 0.068 -0.062 -0.031 0.043
-0.017}}
\dspFunc[linecolor=ColorFour]{x \dspTFM{1}{1 -2.172 2.4 -1.665 0.36 0.556 -0.692 0.174 0.189 -0.292 0.252 -0.067 -0.074 0.086 0.13 -0.12 0.063 -0.081 0.015 -0.031
0.087 -0.138 0.145 -0.141 0.032 0.067 -0.029 -0.059 0.123 -0.112 0.058}}
\dspLegend(0.5, 25){ColorTwo {LPC 10} ColorThree {LPC 20} ColorFour {LPC 30}}
\end{dspPlot}
\caption{DFT of a 20ms voiced speech segment and LPC spectral envelope $1/A(z)$ for different predictor orders.}\label{fig:rn:lpcspec}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Echo Cancellation}
Hands-free telephony, as provided by VOIP applications such as Skype, relies on robust acoustic echo cancellation techniques to be viable. The problem is illustrated in Figure~\ref{fig:rn:voipecho}: in the absence of echo cancellation, when speaker B talks into the microphone, the resulting voice signal is transmitted over the communication channel and played by a loudspeaker at speaker A's location via an amplifier; however, the transmitted voice is picked up by speaker A's microphone and sent back to the original speaker, thereby generating an audible echo that is extremely annoying. Furthermore, since the signal is amplified at the receiver in order to be played out loud, the echo is picked up by B's microphone again in a positive feedback loop that can quickly send the whole communication system into saturation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b!]
\centering
\def\mike{\raisebox{-0.7em}{\includegraphics[height=2em]{\localpath{figs/micro.eps}}}}
\def\loud#1{\raisebox{-1.1em}{\scalebox{#1}[1]{\includegraphics[height=3em]{\localpath{figs/speaker.eps}}}}}
\def\head#1{\raisebox{-0.8em}{\scalebox{#1}[1]{\includegraphics[height=3em]{\localpath{figs/head.eps}}}}}
\begin{dspBlocks}{1}{0.5}
\head{-1} & \rnode{B}{\mike} & & & & & & \loud{1} & \\
A & & & & & & & & \\
& \rnode{A}{\loud{-1}} & & & & & & \mike & \head{1} \\
& & & & & & & & B\\
\end{dspBlocks}
\ncline{->}{1,2}{1,8}
\ncline{->}{3,8}{3,2}
\psset{linecolor=ColorTwo,linewidth=2pt,linestyle=dashed}
\ncarc[nodesep=3pt,arcangle=90,ncurv=1]{->}{A}{B}
\caption{Acoustic echo in a handsfree telephony system.}\label{fig:rn:voipecho}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The coupling between the loudspeaker and the microphone at speaker A's location can be modeled as a transfer function $H(z)$ that captures the characteristics of microphone and loudspeaker and the reverberation pattern in the room; on top of a direct path between loudspeaker and microphone are the multiple reflections caused by the walls and other sound-reflecting surfaces. An echo canceling system tries to estimate $H(z)$ in order to subtract the echo from the signal transmitted back to B, as in Figure~\ref{fig:rn:echocanc}: if $\mathbf{m}$ is the ``true'' voice signal, the signal captured by the microphone will be
\begin{equation} \label{eq:rn:ac_echo}
\mathbf{m}_e = \mathbf{m + h \ast s}
\end{equation}
where $\mathbf{s}$ is the received signal. A good estimate of $H(z)$ allows A to recover $\mathbf{m}$ before the signal is sent back to B (obviously, an identical mechanism needs to be put in place at speaker B's device, to remove the echo associated to B's location). Since objects and people in the room can move at any time, and since the volume of the speaker can be adjusted by the user, the resulting transfer function is not static and therefore we will need to employ an online adaptive estimation such as the one provided by the LMS filter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t!]
\centering
\def\mike{\raisebox{-0.7em}{\includegraphics[height=2em]{\localpath{figs/micro.eps}}}}
\def\loud#1{\raisebox{-1.1em}{\scalebox{#1}[1]{\includegraphics[height=3em]{\localpath{figs/speaker.eps}}}}}
\def\head#1{\raisebox{-0.8em}{\scalebox{#1}[1]{\includegraphics[height=3em]{\localpath{figs/head.eps}}}}}
%
\begin{dspBlocks}{1}{0.3}
& \rnode{B}{\mike} & \BDsub & & & & & \loud{1} & \\
$H(z)~~~~$ & & \BDfilter{$\hat{H}(z)$} & & & & & & \\
& \rnode{A}{\loud{-1}} & \BDsplit & & & & & \mike & \head{1} \\
\end{dspBlocks}
\ncline{->}{1,2}{1,3}\naput{$m_e[n]$}\ncline{->}{1,3}{1,8}\naput{$\hat{m}[n]$}
\ncline{->}{3,8}{3,2}\nbput{$s[n]$}
\ncline{->}{3,3}{2,3}\ncline{->}{2,3}{1,3}
\psset{linecolor=ColorTwo,linewidth=2pt,linestyle=dashed}
\ncarc[nodesep=3pt,arcangle=90,ncurv=1]{->}{A}{B}
\caption{Adaptive acoustic echo cancellation (shown at one end only).}\label{fig:rn:echocanc}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b!]
\centering
\begin{dspBlocks}{1}{0}
$t[n]~$ & & \\
& \BDfilter{hybrid} & & & [mnode=circle]{network} & $~~~\ldots$ \\
$r[n]~$ & & \\
%
\ncangle[angleA=0,angleB=90,linewidth=\BDwidth]{->}{1,1}{2,2}
\ncangle[angleA=0,angleB=-90,linewidth=\BDwidth]{<-}{3,1}{2,2}
\ncline{<->}{2,2}{2,5}\naput{ $s[n]$}\ncline{<->}{2,5}{2,6}
\end{dspBlocks}
\caption{Full duplex transmission over a single wire.}\label{fig:rn:hybrid}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Another scenario where echo cancellation is of paramount importance is the case of \textit{full-duplex} communication over wired transmission lines, as in the case of voiceband modems and ADSL. The user's device is connected from the home to the telephone network using a single wire which, during full duplex operation, will simultaneously carry both the transmitted and received signals added together. Since the transmitted signal is known, the received signal can be extracted from the single-wire signal via subtraction. The combination and separation of signals from two independent processing lines to a single transmission line is performed by a circuit called the ``hybrid'', as in Figure~\ref{fig:rn:hybrid}. Because of imperfections in the hybrid, as well as impedance mismatches along the analog line, the received signal $\mathbf{r}$ will contain ``traces'' of the transmitted signal $\mathbf{t}$; there will be at least a near echo, due to the non-ideal hybrid, and a far echo due to reflections at the other end of the communication channel that travel back along the line. In other words, we will have to process a received signal
\begin{equation} \label{eq:rn:comm_echo}
\mathbf{r}_e = \mathbf{r + h \ast t}.
\end{equation}
Once again, the goal of the echo canceler is to estimate the transfer function that approximates the effects of near and far echoes in order to recover the true $\mathbf{r}$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t!]
\centering
\begin{dspBlocks}{1.2}{1.2}
$x[n]~$ & \BDsplit & & \BDfilter{$\hat{H}(z)$} & \BDsplit & $~~~~\hat{d}[n]$ \\
& & \BDfilter{$H(z)$} & \BDsub &
\psset{arrows=->}
\ncline{1,1}{1,4}\ncline{1,4}{1,6}
\ncline{-}{1,2}{2,2}\ncline{2,2}{2,3}
\ncline{2,3}{2,4}\taput{$~~~d[n]$}
\ncline{-}{1,5}{2,5}\ncline{2,5}{2,4}
\ncline[arrowscale=1.2,arrowinset=0,ArrowFill=false]{2,4}{1,4}\nbput{$e[n]$}
\end{dspBlocks}
\caption{Echo estimation via adaptive filtering.}\label{fig:rn:echoest}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In general, the echo cancellation problem can be formalized as in Figure~\ref{fig:rn:echoest}. Given a known signal $\mathbf{x}$ and an echo-corrupted version $\mathbf{d[n] = h[n] \ast x}$, the goal is to estimate $\mathbf{h}$ using the LMS algorithm.\footnote{
%
Note that, compared to the channel equalization example illustrated in Section~\ref{sec:rn:adaptive}, in this case the roles of ``input'' and ``desired'' signal are reversed since we want to estimate the echo in order to subtract it.}
%
With respect to the examples we just introduced, we can see that in order to perform the estimation we need a ``clean'' version of the echo signal $\mathbf{d}$. In the case of full duplex data communication this is achieved by having only one transmitter operate at a time during a specific adaptation phase (i.e., while a device adapts its LMS filter, the other device does not transmit, which corresponds to having $\mathbf{r = 0}$ in~(\ref{eq:rn:comm_echo})). For hands-free telephony, the problem is somewhat alleviated by the fact that people rarely speak simultaneously when talking over the phone, which creates periods in which $\mathbf{m = 0}$ in~(\ref{eq:rn:ac_echo}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b]
\centering
\begin{dspBlocks}{1}{1}
% 1 2 3 4 5
$x[n]~~$ & \BDadd & & & \BDsplit & $~~d[n]$ \\
& & \BDfilter{$H(z)$} & \BDdelayN{M} &
\ncline{->}{1,1}{1,2}\ncline{->}{1,2}{1,6}
\ncline{->}{2,2}{1,2}\ncline{2,2}{2,3}\taput{$\alpha$}\ncline{<-}{2,3}{2,4}\ncline{->}{2,5}{2,4}
\ncline{2,5}{1,5}
\end{dspBlocks}
\caption{Echo simulation via a feedback loop and a lowpass filter.}\label{fig:rn:synecho}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
As an example, consider the synthetic echo system shown in Figure~\ref{fig:rn:synecho}; it consists of a feedback loop with a delay of $M$ samples in cascade with a lowpass filter. The delay $M$ models the round-trip delay of the echo, whereas the lowpass filter introduces a progressively increasing attenuation that models the effects of a channel over the successive reflecions of the original signal. Figure Figure~\ref{fig:rn:synechoimpulse} shows the first 500 taps of the impulse response when the lowpass filter is a simple leaky integrator. We can use the LMS algorithm to estimate a portion of the impulse response using the iterative update in~(\ref{eq:rn:lmserr}) and~(\ref{eq:rn:lmspred}). The fastest convergence is achieved by using a white noise sequence as the input $\mathbf{x}$; this corresponds to an autocorrelation matrix proportional to the identity matrix, which makes the least squares problem well-conditioned. Figure~\ref{fig:rn:LMSMSE} shows the corresponding evolution of the mean squared error signal $\expt{e^2[n]}$ as an ensemble average over 200~experiments. Figure~\ref{fig:rn:LMSMSETV} shows a similar error curve for a nonstationary echo signal in which the echo delay $M$ changes abruptly to demonstrate how the continuous LMS adaptation quickly recovers from a change in signal statistics.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b!]
\centering
\begin{dspPlot}[xticks=100,xout=true]{0,500}{-0.4,1.1}
\dspFuncDataAt[linecolor=ColorDT]{0}{%
1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. -0.32 -0.192 -0.115 -0.069 -0.041 -0.025 -0.015 -0.009 -0.005
-0.003 -0.002 -0.001 -0.001 -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. 0.102 0.123 0.111 0.088 0.066 0.048 0.033
0.023 0.015 0.01 0.007 0.004 0.003 0.002 0.001 0.001 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. -0.033 -0.059 -0.071 -0.071 -0.064
-0.054 -0.043 -0.033 -0.025 -0.018 -0.013 -0.009 -0.006 -0.004 -0.003
-0.002 -0.001 -0.001 -0.001 -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
-0. -0. -0. -0. -0. -0. -0. -0. 0.01 0.025
0.038 0.045 0.048 0.046 0.041 0.035 0.029 0.023 0.018 0.014
0.01 0.008 0.006 0.004 0.003 0.002 0.001 0.001 0.001 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.}
\end{dspPlot}
\caption{Impulse response (first 500 taps) of the system in Figure~\ref{fig:rn:synecho} for $M=100$ and $\alpha=-0.8$ and in which $H(z)$ is a leaky integrator with $\lambda=0.6$.}\label{fig:rn:synechoimpulse}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering
\begin{dspPlot}[xticks=custom,xout=true,sidegap=0,xlabel={iterations}]{0,500}{0,1.5}
\dspCustomTicks[]{0 0 250 2500 500 5000}
\dspFuncDataAt[linecolor=ColorDT]{0}{%
1.086 1.414 1.378 1.23 1.195 1.248 0.985 1.318 0.989 1.073
1.057 1.022 0.92 1.017 1.032 1.086 0.968 0.981 1.059 1.121
0.959 0.908 1.021 0.967 0.902 1.02 0.817 1.077 0.891 0.909
0.807 0.862 0.836 0.93 0.973 0.977 0.895 0.698 1.014 0.786
0.753 0.897 0.81 0.821 0.804 0.881 0.755 0.791 0.862 0.669
0.735 0.891 0.732 0.807 0.724 0.638 0.647 0.668 0.741 0.6 0.537
0.691 0.64 0.673 0.601 0.581 0.621 0.673 0.592 0.567 0.511
0.489 0.554 0.468 0.514 0.491 0.489 0.466 0.479 0.564 0.542
0.459 0.529 0.525 0.497 0.458 0.422 0.534 0.457 0.462 0.507
0.411 0.399 0.414 0.367 0.376 0.415 0.372 0.351 0.353 0.411
0.373 0.404 0.364 0.354 0.358 0.34 0.368 0.36 0.298 0.359
0.347 0.27 0.261 0.293 0.284 0.307 0.289 0.297 0.297 0.266
0.238 0.293 0.285 0.3 0.272 0.272 0.233 0.276 0.294 0.231
0.227 0.245 0.226 0.214 0.229 0.223 0.238 0.234 0.249 0.21
0.207 0.198 0.202 0.205 0.239 0.228 0.202 0.231 0.219 0.185
0.203 0.203 0.187 0.17 0.186 0.161 0.174 0.18 0.172 0.162
0.166 0.173 0.176 0.155 0.139 0.17 0.156 0.16 0.17 0.146
0.14 0.15 0.155 0.148 0.119 0.117 0.154 0.146 0.143 0.124
0.147 0.128 0.133 0.112 0.125 0.136 0.13 0.115 0.123 0.123
0.109 0.094 0.119 0.116 0.121 0.121 0.11 0.108 0.11 0.112
0.114 0.11 0.089 0.106 0.125 0.113 0.109 0.109 0.095 0.081
0.084 0.111 0.099 0.082 0.086 0.105 0.089 0.088 0.101 0.088
0.071 0.092 0.074 0.082 0.077 0.078 0.083 0.077 0.089 0.073
0.073 0.091 0.075 0.077 0.061 0.071 0.079 0.063 0.083 0.074
0.065 0.073 0.066 0.064 0.059 0.068 0.069 0.071 0.067 0.076
0.064 0.065 0.076 0.064 0.064 0.066 0.061 0.055 0.065 0.062
0.07 0.065 0.067 0.071 0.059 0.062 0.06 0.061 0.056 0.046
0.061 0.054 0.052 0.059 0.056 0.056 0.058 0.044 0.054 0.056
0.055 0.053 0.047 0.05 0.053 0.05 0.044 0.043 0.051 0.052
0.047 0.054 0.045 0.042 0.049 0.046 0.043 0.047 0.05 0.051
0.043 0.044 0.052 0.042 0.044 0.041 0.041 0.04 0.038 0.043
0.037 0.039 0.04 0.033 0.041 0.039 0.048 0.039 0.042 0.035
0.036 0.047 0.034 0.033 0.036 0.041 0.042 0.033 0.044 0.04
0.035 0.04 0.041 0.032 0.042 0.033 0.032 0.036 0.038 0.036
0.031 0.034 0.036 0.037 0.038 0.04 0.029 0.032 0.038 0.042
0.029 0.036 0.032 0.032 0.038 0.037 0.036 0.031 0.029 0.037
0.031 0.03 0.029 0.034 0.026 0.033 0.026 0.028 0.026 0.038
0.03 0.039 0.029 0.033 0.031 0.028 0.029 0.024 0.026 0.032
0.029 0.033 0.029 0.034 0.031 0.028 0.033 0.028 0.024 0.028
0.027 0.029 0.026 0.036 0.027 0.025 0.031 0.027 0.03 0.026
0.028 0.028 0.025 0.034 0.027 0.029 0.028 0.029 0.029 0.026
0.027 0.029 0.029 0.026 0.026 0.026 0.022 0.026 0.03 0.026
0.028 0.027 0.028 0.024 0.029 0.024 0.023 0.028 0.03 0.027
0.026 0.026 0.03 0.026 0.028 0.027 0.023 0.028 0.027 0.026
0.025 0.029 0.028 0.032 0.027 0.026 0.022 0.024 0.029 0.024
0.023 0.029 0.027 0.027 0.027 0.023 0.025 0.026 0.025 0.025
0.025 0.026 0.026 0.025 0.026 0.022 0.022 0.023 0.022 0.025
0.022 0.025 0.023 0.025 0.025 0.021 0.021 0.022 0.02 0.024
0.023 0.025 0.023 0.024 0.032 0.023 0.02 0.019 0.019 0.023
0.023 0.025 0.027 0.025 0.03 0.024 0.018 0.024 0.025}
\end{dspPlot}
\caption{MSE as a function of the iteration for the LMS estimation of the impulse response in Figure~\ref{fig:rn:synechoimpulse} (ensemble average over 200 experiments). Filter length is 500 taps, $\mathbf{x}$ is a unit-variance Gaussian white noise process and the step size is $\alpha=0.0008$.}\label{fig:rn:LMSMSE}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering
\begin{dspPlot}[xticks=custom,xout=true,sidegap=0,xlabel={iterations}]{0,500}{0,1.5}
\dspCustomTicks[]{0 0 250 2500 500 5000}
\dspFuncDataAt[linecolor=ColorDT]{0}{%
1.259 1.336 1.059 1.116 1.176 1.177 1.414 1.257 1.23 1.095
1.326 1.21 1.119 0.984 1.025 0.96 1.187 0.985 1.132 1.066
0.975 1.175 0.954 1.156 0.949 1.06 0.934 0.877 0.914 0.774
0.832 0.79 0.865 1.043 0.886 0.909 0.826 0.902 0.773 0.742
0.833 0.697 0.718 0.716 0.969 0.839 0.762 0.819 0.742 0.712
0.717 0.717 0.645 0.732 0.742 0.762 0.614 0.764 0.606 0.574
0.682 0.655 0.509 0.587 0.618 0.583 0.703 0.49 0.586 0.497
0.622 0.581 0.531 0.426 0.59 0.425 0.548 0.549 0.483 0.434
0.496 0.468 0.529 0.466 0.385 0.471 0.397 0.471 0.452 0.466
0.422 0.393 0.378 0.385 0.391 0.364 0.316 0.402 0.405 0.351
0.392 0.385 0.386 0.408 0.334 0.338 0.271 0.336 0.335 0.333
0.278 0.33 0.322 0.276 0.292 0.322 0.329 0.329 0.262 0.269
0.264 0.287 0.287 0.272 0.226 0.264 0.244 0.316 0.242 0.264
0.205 0.253 0.269 0.198 0.215 0.269 0.235 0.256 0.244 0.237
0.174 0.23 0.197 0.21 0.237 0.187 0.234 0.178 0.182 0.216
0.24 0.197 0.164 0.182 0.177 0.203 0.184 0.198 0.159 0.177
0.164 0.152 0.152 0.179 0.151 0.189 0.181 0.164 0.187 0.151
0.168 0.147 0.141 0.129 0.162 0.141 0.139 0.149 0.151 0.144
0.154 0.141 0.154 0.132 0.144 0.117 0.116 0.15 0.143 0.144
0.122 0.126 0.139 0.095 0.115 0.113 0.108 0.112 0.113 0.103
0.104 0.129 0.098 0.099 0.121 0.089 0.085 0.114 0.106 0.095
0.093 0.098 0.103 0.099 0.106 0.096 0.084 0.105 0.092 0.095
0.086 0.09 0.086 0.083 0.082 0.064 0.08 0.091 0.1 0.082
0.065 0.095 0.093 0.086 0.065 0.073 0.076 0.094 0.081 0.081
0.062 0.09 0.068 0.081 0.067 0.064 0.072 0.077 0.064 0.062
0.068 0.073 0.071 0.069 0.064 0.067 0.057 0.071 0.06 0.064
0.063 0.068 0.064 0.061 0.059 0.052 0.048 0.052 0.054 0.062
0.048 0.057 0.057 0.047 0.052 0.057 0.049 0.056 0.052 0.051
0.052 0.048 0.047 0.044 0.566 0.466 0.535 0.51 0.445 0.442
0.41 0.461 0.442 0.492 0.396 0.474 0.48 0.391 0.493 0.349
0.415 0.37 0.404 0.397 0.451 0.394 0.443 0.406 0.364 0.366
0.415 0.378 0.414 0.397 0.353 0.368 0.347 0.333 0.354 0.337
0.353 0.353 0.371 0.349 0.32 0.359 0.362 0.338 0.351 0.325
0.325 0.326 0.356 0.354 0.299 0.351 0.35 0.296 0.318 0.316
0.312 0.254 0.258 0.284 0.236 0.306 0.232 0.22 0.251 0.273
0.241 0.274 0.256 0.204 0.229 0.204 0.251 0.225 0.237 0.205
0.239 0.19 0.202 0.198 0.207 0.201 0.215 0.204 0.198 0.241
0.193 0.169 0.183 0.198 0.19 0.18 0.153 0.161 0.152 0.173
0.147 0.152 0.163 0.152 0.19 0.143 0.156 0.144 0.135 0.119
0.148 0.135 0.131 0.157 0.125 0.157 0.124 0.139 0.107 0.128
0.133 0.115 0.138 0.107 0.1 0.128 0.132 0.124 0.124 0.111
0.12 0.111 0.103 0.101 0.085 0.117 0.113 0.116 0.104 0.095
0.105 0.101 0.102 0.091 0.107 0.098 0.085 0.093 0.09 0.097
0.098 0.098 0.08 0.083 0.089 0.09 0.084 0.079 0.086 0.089
0.078 0.091 0.082 0.102 0.083 0.072 0.069 0.071 0.071 0.076
0.073 0.073 0.074 0.061 0.064 0.06 0.068 0.067 0.064 0.062
0.08 0.066 0.06 0.062 0.055 0.061 0.054 0.075 0.06 0.067
0.063 0.06 0.052 0.05 0.056 0.05 0.062 0.056 0.056 0.052
0.054 0.061 0.055 0.051 0.047 0.048 0.049 0.04 0.049 0.048
0.044 0.053 0.044 0.044 0.045 0.043 0.041 0.043 0.044 0.042}
\end{dspPlot}
\caption{MSE for the LMS estimation as in Figure~\ref{fig:rn:LMSMSE} when the echo delay in $\mathbf{d}$ changes from $M=100$ to $M=90$ at $n=3000$.}\label{fig:rn:LMSMSETV}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c4science · Help