Homec4science
Diffusion epfl (master) Browse

R2653/writing/sp4comm.multipub/50-filters/00-lf-intro.tex
e6025b0b48e7master

Last Change
Raw File

00-lf-intro.tex
View Options

%5a
\chapter{Discrete-Time Filters}
\label{ch:linearfilters}\label{ch:lf}
In the previous chapters we tackled discrete-time signals from an analytical point of view: we put in place a taxonomy and we described the tools used to unlock a frequency domain representation. We will now move on to the {\em processing} part of our journey and will start by introducing the most fundamental device for modifying the content of a discrete-time signal, namely, the linear time-invariant system. As we will see momentarily, LTI systems (or filters, for short) admit a full, explicit mathematical characterization of their properties -- a remarkable exception in the universe of all possible signal processing devices. Because of this, they occupy a fundamental position in the arsenal of the signal processing engineer.
In this chapter we will begin by illustrating the fundamental properties of LTI systems and by introducing a pair of simple intuitive examples that we will use repeatedly in the rest of this book both for theoretical and for applied purposes. In the next two chapters, we will continue our exploration by looking at the algorithmic nature of practical discrete-time filters and at some of the established techniques that are used in filter design.
%The previous chapters gave us a thorough overview on both
%the nature of discrete-time signals and on the tools used in
%analyzing their properties. In the next few Chapters, we will
%study the fundamental building block of any digital signal
%processing system, that is, the linear filter. In the
%discrete-time world, filters are nothing but procedures
%which store and manipulate mathematically the numerical
%samples appearing at their input and their output; in other
%words, any discrete-time filter can be described
%procedurally in the form of an algorithm. In the special
%case of linear and time-invariant filters, such an algorithm
%can be concisely described mathematically by a
%constant-coefficient difference equation.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Linear Time-Invariant Systems}\index{system!LTI|(}
In its most general form, a {\em discrete-time system} is a processing unit that accepts one or more discrete-time sequences as inputs and produces one or more discrete-time sequences as outputs. The inner workings of such a device can be arbitrarily complex but a first basic requirement here is that the system acts synchronously, that is, at each increment of the ``time'' index $n$ for the input sequences, a new set of results appears at the outputs. This is in contrast with {\em multirate}\/ systems, which we will study in Chapter~\ref{ch:mr}, where the data rates at the input and at the output can differ.
A special subset of discrete-time system, and the one we will focus on, is the class of single-input, single-output {\em linear and time-invariant} (LTI) systems, more commonly known as ``filters'', and represented schematically in Figure~\ref{fig:lf:system}. A discrete-time filter is an operator that transforms an input sequence into an output sequence:
\begin{equation}
\mathbf{y} = \mathcal{H}\mathbf{x};
\end{equation}
linearity implies that for any signal $\mathbf{x}$ and $\mathbf{y}$ and any scalar $\alpha$ we have
\begin{align}
\mathcal{H}(\alpha\mathbf{x}) &= \alpha(\mathcal{H}\mathbf{x}) \\
\mathcal{H}(\mathbf{x + w}) &= \mathcal{H}\mathbf{x} + \mathcal{H}\mathbf{w}
\end{align}
whereas time-invariance implies that
\begin{equation}
\mathcal{H}(\mathcal{D}^k\mathbf{x}) = \mathcal{D}^k(\mathcal{H}\mathbf{x})
\end{equation}
where $\mathcal{D}^k$ is the delay-by-$k$ operator\footnote{
See Section~\ref{dt:notation} for a recap on the notation used for signal operators.
}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center\small
\begin{dspBlocks}{1.2}{0}
$\mathbf{x}$~~ & \rnode{F}{\BDfilter{$\mathcal{H}$}} & ~~$\mathbf{y}$
\ncline{->}{1,1}{1,2}\ncline{->}{1,2}{1,3}
\end{dspBlocks}
\caption{Representation of a discrete-time filter.}\label{fig:lf:system}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In intuitive terms, linearity expresses the natural requirement that if we filter a source that is the combination of independent signals, the result is the superposition of what we obtain by filtering each signal independently; a physical analogy is that of, say, an amplifier connected to an electric guitar, where the sound obtained by plucking multiple strings at once is the sum of the sound of each string plucked independently. Similarly, time-invariance describes the natural requirement that, regardless of when the filter is ``switched on,'' its effects will remain the same, only shifted in time. Again, this is the idealized behavior of most physical systems such as an amplifier, once we neglect the slight changes due to the aging of the internal electronic components.
As we mentioned earlier, the linearity and time invariance conditions are not always fulfilled or even desired in a processing system: a distorted electric guitar, for instance, which is the defining sound of over half a century of rock music, is achieved by connecting the instrument to a highly nonlinear device. Non-LTI building blocks are also used ubiquitously in communication system devices; even a simple radio transmitters involves at least a time-variant operation as we will see later in the chapter. Yet, the versatility and mathematical attractiveness of LTI systems make them an extremely fruitful case study among signal processing tools.
\subsection{Two Examples}\label{sec:lf:examples}
To better understand the properties of LTI system that we will study in the rest of the chapter, let's consider a practical example:
the data shown in Figure \ref{fig:lf:noisyUnfiltered} can be interpreted as the sum of a clean, smooth, and slow-varying signal plus a fast-varying, noisy perturbation that appears as little wiggles in the plot; we would like to process the signal to recover its slow-varying component and eliminate the additive noise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center\small
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-1.txt}}
\end{dspPlot}
\caption{Portion of a discrete-time noisy signal.}\label{fig:lf:noisyUnfiltered}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\itempar{The Moving Average.}
\label{sec:lf:MovingAverage}
\index{moving average|(}
An intuitive and basic approach to removing noise from a data set is to replace each sample by a
\emph{local average} of neighboring samples, which will work if the underlying trend is smooth and if the noise is varies independently of and much more rapidly than the clean data. For instance, we could build a denoised sequence $\mathbf{y}$ by computing each output sample as
\begin{equation} \label{eq:lf:ma}
y[n] = \frac{1}{M}\sum_{k = 0}^{M-1} x[n-k],
\end{equation}
that is, by taking the average of $M$ consecutive input samples. It is easy to verify numerically that the approach does work, and that the smoothing power of the system is dependent on $M$, as shown in Figure~\ref{fig:lf:noisyFIRFiltered}. Note that the algorithm in~(\ref{eq:lf:ma}) uses only the current and \textit{past} samples, so that it can work in ``on the fly'' by simply keeping track of the past $M-1$ input values; note also that, since the output depends only on a finite number of input samples, the algorithm exhibits so-called finite-memory effects, that is, input values farther in the past than $M$ samples do not influence the output anymore.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\center\small
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-2.1.txt}}
\rput[tr](500,7){$M=2$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-2.2.txt}}
\rput[tr](500,7){$M=4$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-2.3.txt}}
\rput[tr](500,7){$M=12$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-2.4.txt}}
\rput[tr](500,7){$M=100$}
\end{dspPlot}
\caption{Effects of the moving average filter for different values of $M$.}\label{fig:lf:noisyFIRFiltered}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Equation~(\ref{eq:lf:ma}) defines a linear, time-invariant transformation $\mathcal{H}$ called the \textit{moving average} (MA) filter; the LTI nature of the operation can be easily confirmed by writing the input-output relationship as
\begin{equation}
\mathcal{H}\mathbf{x} = \frac{1}{M}\sum_{k = 0}^{M-1} \mathcal{D}^k\mathbf{x}
\end{equation}
and noticing that it involves only sums, scalar multiplications (that is, linear operations) and delays (that is, time-invariant operations). Alternatively, we can verify the LTI properties explicitly like so:
\[
(\mathcal{H}(\alpha \mathbf{x}))[n] = \frac{1}{M}\sum_{k = 0}^{M-1} \alpha x[n-k] = \alpha ((\mathcal{H}\mathbf{x})[n])
\]
\[
(\mathcal{H}(\mathbf{x + w}))[n] = \frac{1}{M}\sum_{k = 0}^{M-1} (x[n-k] + w[n-k]) = (\mathcal{H}\mathbf{x})[n] + (\mathcal{H}\mathbf{w})[n]
\]
\[
(\mathcal{H}(\mathcal{D}^m \mathbf{x}))[n] = \frac{1}{M}\sum_{k = 0}^{M-1} x[(n - m) - k] = (\mathcal{H}\mathbf{x})[n-m]
\]
\itempar{The Leaky Integrator.}
\label{sec:lf:leakyIntegrator}%
\index{leaky integrator|(}\index{integrator!leaky}
The moving average filter, while effective, has a computational cost that is linear in $M$, both in terms of number of operations per output sample and in terms of memory storage; since effective noise removal seems to require large values of $M$, we should be mindful of these costs. An alternative, and ``cheaper'', smoothing algorithm can be derived by writing out the formula for a moving average filter of length $M$
\[
y_M[n] = \frac{1}{M}\sum_{k=0}^{M-1} x[n-k]
\]
and by observing that
\begin{align*}
y_{M}[n] & = \frac{M-1}{M} \, y_{M-1}[n-1] + \frac{1}{M}\, x[n] \\
& = \lambda y_{M-1}[n-1] + (1-\lambda ) x[n]
\end{align*}
where we have defined $\lambda = (M-1)/M$. Now, as $M$ grows large, we can safely assume that the averages over $M-1$ and over $M$ points are going to be very similar; in other words, for $M$ large, we can say that $y_{M-1}[n] \approx y_{M}[n]$. This leads to the following \emph{recursive} algorithm
\begin{equation}\label{eq:lf:leaky}
y[n] = \lambda y[n-1] + (1-\lambda) x[n]
\end{equation}
in which the output is computed as a linear combination of the current input sample and the past \textit{output} sample and whose computational cost is fixed and very small (two multiplications and one addition). Before analyzing the structure and the implications of this new approach, we can easily verify numerically that it works; Figure~\ref{fig:lf:noisyIIRFiltered} shows the effects of the algorithm on the input sequence in Figure~\ref{fig:lf:noisyUnfiltered} for different values of $\lambda$ and we can notice that the smoothing power increases as $\lambda$ approaches unity.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[p]
\center\small
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-3.1.txt}}
\rput[tr](500,7){$\lambda=0.2$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-3.2.txt}}
\rput[tr](500,7){$\lambda=0.5$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-3.3.txt}}
\rput[tr](500,7){$\lambda=0.8$}
\end{dspPlot}
\begin{dspPlot}[xticks=100,yticks=2,sidegap=0,xout=true]{-50,550}{-8, 8}
\dspFuncFile[linecolor=ColorDT]{\localpath{data/filtering-example-3.4.txt}}
\rput[tr](500,7){$\lambda=0.98$}
\end{dspPlot}
\caption{Effects of the leaky integrator for different values of $\lambda$.}\label{fig:lf:noisyIIRFiltered}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Equation~(\ref{eq:lf:leaky}) defines a linear, time-invariant transformation $\mathcal{H}$ called the \textit{leaky integrator} but the LTI nature of the algorithm is a bit trickier to show in this case because the equation is recursive. Let's begin by asking ourselves how we are going to compute the output sequence if every output sample depends on its predecessor: isn't this a classic chicken and egg situation? The way out of this impasse is to impose so-called \textit{zero initial conditions:} we assume there is a starting time $n_0$ at which we ``switch on'' the filter and begin to compute the output and that, for all $n < n_0$, both $y[n]$ and $x[n]$ are identically zero. With these initial conditions in place, we can first proceed to show that if $\mathbf{y} = \mathcal{H}\mathbf{x}$, then $\mathbf{y}$ is the \textit{unique} response of the system to $\mathbf{x}$. Indeed, we first observe that the equation
\[
y[n] - \lambda y[n-1] - (1-\lambda)x[n] = 0
\]
is satisfied for all values of $n$. Assume now there were another sequence, distinct from $\mathbf{y}$, such that $\mathbf{w} = \mathcal{H}\mathbf{x}$; this implies that also
\[
w[n] - \lambda w[n-1] - (1-\lambda)x[n] = 0
\]
is valid for all $n$ so that it must be
\[
(y[n] - w[n]) = \lambda(y[n-1] - w[n-1]) \qquad \forall n.
\]
Because of zero initial conditions, there exists a time $n_0$ so that $y[n] = w[n] = 0$ for $n < n_0$ and therefore $y[n] = w[n]$ for all $n \ge n_0$ as well.
With this in place, and keeping $\mathbf{y} = \mathcal{H}\mathbf{x}$, consider now the response to a scaled input, $\mathbf{y}_1 = \mathcal{H}(\alpha \mathbf{x})$. We have
\[
y_1[n] - \lambda y_1[n-1] - \alpha(1-\lambda)x[n] = 0 \qquad \forall n
\]
and, if we replace $y_1[n] = \alpha y[n]$ in the above expression, we can see that the equation is satisfied so that $\mathbf{y}_1 = \alpha \mathbf{y}$ is indeed the unique solution. With the same technique we can prove linearity with respect to the sum of two sequences and demonstrate time invariance.
Finally, a word on the origin of the name ``Leaky Integrator''. We have seen in Section~\ref{sec:dt:operators} the discrete-time integration operator, defined as
\[
\left( \mathcal{S}\mathbf{x} \right)[n] = \sum_{m = -\infty}^{n} x[m];
\]
clearly the sequence $\mathbf{y} = \mathcal{S}\mathbf{x}$ can be computed recursively as
\[
y[n] = y[n-1] + x[n]
\]
where $y[n]$ represent the running total accumulated up to $n$. In the Leaky Integrator the factor $\lambda$ (which is usually close to one) acts as a ``leakage'' or ``forgetting'' factor for the running total, which is also incremented at each step not by the full value of the input sample but by a fraction $1-\lambda$ of it.
\begin{comment}
In digital filters, time-invariance implies that the input samples used to determine the output at time $n$ are at a fixed distance from $n$; visually, imagine to line up the input and output sequences and draw a pattern that connects each output sample to all the input samples used in its computation; in a time-invariant system, as you move along the time axis, the pattern is shifted but it never changes shape. For instance, here is the input-output pattern in a simple moving average filter whose output can be described by the equation $y[n] = (x[n-1] + x[n])/2$, and it is obvious that the shape remains the same for all values of $n$:
\begin{center}
\vspace{1ex}
\footnotesize
\begin{dspBlocks}{0.4}{1}
$\ldots$ & $x[-4]$ & $x[-3]$ &$x[-2]$ & $x[-1]$ & $x[0]$ & $x[1]$ & $x[2]$ & $x[3]$ & $x[4]$ & $x[5]$ & $\ldots$ \\
$\ldots$ & $y[-4]$ & $y[-3]$ &$y[-2]$ & $y[-1]$ & $y[0]$ & $y[1]$ & $y[2]$ & $y[3]$ & $y[4]$ & $y[5]$ & $\ldots$
\psset{linecolor=ColorOne,linewidth=2pt}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,3}{2,3}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,2}{2,3}
\psset{linecolor=ColorTwo}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,6}{2,6}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,5}{2,6}
\psset{linecolor=ColorThree}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,9}{2,9}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,8}{2,9}
\end{dspBlocks}
\vspace{1ex}
\end{center}
By contrast, if we consider the transformation $y[n] = x[n] + x[2n]$, then the shape of the pattern changes at each step since the underlying operator is not time-invariant:
\begin{center}
\vspace{1ex}
\footnotesize
\begin{dspBlocks}{0.4}{1}
$\ldots$ & $x[-4]$ & $x[-3]$ &$x[-2]$ & $x[-1]$ & $x[0]$ & $x[1]$ & $x[2]$ & $x[3]$ & $x[4]$ & $x[5]$ & $\ldots$ \\
$\ldots$ & $y[-4]$ & $y[-3]$ &$y[-2]$ & $y[-1]$ & $y[0]$ & $y[1]$ & $y[2]$ & $y[3]$ & $y[4]$ & $y[5]$ & $\ldots$
\psset{linecolor=ColorOne,linewidth=2pt}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,4}{2,4}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,2}{2,4}
\psset{linecolor=ColorTwo}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,6}{2,6}
%q\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,5}{2,6}
\psset{linecolor=ColorThree}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,7}{2,7}
\ncangle[angleA=90,angleB=90,nodesepA=-1.2em,nodesepB=.2em,armB=1em]{->}{1,8}{2,7}
\end{dspBlocks}
\vspace{1em}
\end{center}
\end{comment}
c4science · Help