Lectures Lecture Topics

MODULE: 1 Multirate Digital Signal Processing

1. Introduction
2. Decimation by a factor D, Interpolation by a factor I
3. Sampling rate conversion by rational factor I/D
4. Filter Design and Implementation for sampling-rate
5. Multistage implementation of sampling rate conversion
6. Sampling rate conversion of Band pass signal
7. Application of multi rate signal processing: design of phase shifters
8. Implementation of narrowband lowpass filters
9. Implementation of Digital filter banks.
10. Filter Bank and Sub-band Filters and its applications.

MODULE: 2 Linear prediction and Optimum Linear Filters

11. Innovations Representation of a stationary random process
12. Forward Linear Prediction
13. Backward Linear Prediction
14. Solution of the normal equations
15. Solution of the normal equations
16. Properties of the linear prediction-error filters
17. AR lattice- ladder filters
18. ARMA lattice- ladder filters
19. Wiener filter for filtering and Prediction: FIR Wiener Filter
20. Orthogonality Principle in linear mean square estimation.
MODULE: 3 Power Spectrum Estimation
21. Estimation of spectra from finite-duration observation of signals
22. Non parametric method for power spectrum estimation
23. Bartlett method
24. Blackman method
25. Turkey method
26. parametric method for power estimation
27. Yuke-Walker method
 28. Burg method

29. MA model and ARMA model

MODULE: 4 Adaptive Signal Processing

30. Introduction to adaptive signal processing

31. Application of adaptive signal processing: System identification, Channel

equalization, adaptive noise cancellation, adaptive line enhancer
32. Application of adaptive signal processing: System identification, Channel
equalization, adaptive noise cancellation, adaptive line enhancer
33. Adaptive linear combiner

34. Basics of Wiener filtering

35. Widrow-Hopf Equation

36. Least mean square algorithm

37. Recursive least square

algorithm
38. variants of LMS algorithm: FX-LMS, Fast LMS, N-LMS, PN-LMS

39. Continued... (FX-LMS, Fast LMS, N-LMS, PN-LMS)

40. Design of Adaptive FIR & IIR filters

Text Book
1. Digital Signal Processing, Third Edition, J.G. Proakis and D.G. Manolakis, Prentice Hall.
2. Adaptive Signal Processing, B. Widrow and Stern,

Reference Book
1. Digital Signal Processing, by Sanjit K Mitra, new edition, TMH.
2. Digital Signal Processing, by Salivahanan, new edition, TMH.
 Module 1.
Multirate Digital Signal
Processing
INTRODUCTION
The process of converting a signal from a given rate to a different rate is called sampling rate
conversion. Systems that employ multiple sampling rates in the processing of digital signals are
called multirate digital signal processing systems. The two basic operations in a multirate system
are decreasing (decimation) and increasing (interpolation) the sampling-rate of a signal.

DECIMATION
It is also called as down sampling. In this,the sampling rate of a discrete-time signal x(n) with
sampling frequency Fs is reduced to a discrete-time signal y(n) of sampling frequency Fs/D, D is
down sampling factor .The simplest way of doing so is to discard (D-1) samples for every D
samples in original sequence.

Mathematically:

If x(n) is the original discrete-time signal

Then y(n) = x(nD)

So, in z domain,

Y(z) = x(nD) z-n

Let nD = p => n= p/D

/
So Y(z) = x(p)(z )

=>Y(z) = X(z 1/D)

On unit circle, Y( ) = X( /D)

It seems hence that spectrum of original signal gets stretched as an effect of down sampling.
Hence, an anti-aliasing digital filter is used before down-sampling to prevent aliasing.

The input sequence x(n) is passed through a lowpass filter, characterized by the impulse
response h(n) and a frequency response HD which ideally satisfies the condition

1, | | /D
HD
0, otherwise
 x(n) v(n) y(m)
h(n) Downsampler

Fx = 1/T x Fy = Fx/D

Decimation by a factor D

Time domain expression:

v(n) = x(n) * h(n)

= ( ) ( )

y(m) = v(mD)

= ( ) ( )

( ) | ( )|

- 0 0

Spectra of signals during decimation by factor

 | ( )| ( )

0 - 0

Spectra of signals during decimation by factor D

INTERPOLATION
It is also called as up sampling. In this,the sampling rate of a discrete-time signal x(n) with
sampling frequency Fs is increased to a discrete-time signal y(n) of sampling frequency IFs, I is
up sampling factor .The simplest way of doing so is to insert (I - 1) samples for every I samples
in original sequence.

Mathematically:

If x(n) is the original discrete-time signal

Then y(n) = x(n/I)

So, in z domain,

Y(z) = x(n/I) z-n

Let n/I= p => n= pI

So Y(z) = x(p)(z )

=>Y(z) = X(z I)

So interpolation causes the original spectrum to get compressed by I-folds. It yields undesirable
replicas in the signal’s frequency spectrum. Hence, it is necessary to remove these replicas from
the frequency spectrum. So the expansion process is followed by a unique digital low-pass filter
called an anti-imaging filter.
 C ,0 | | /I
HI
0 , otherwise

CX( I) ,0 | | /I
Y( ) =
0 , otherwise

x(n) Upsampler v(m) y(m)

h(m)

Fx = 1/Tx Fy = IFx

Interpolation of signal by a factor I

Time domain expression:

, = 0, ± , ±2 , … .
v(m) =
0 ,

y(m) = v(m) * h(m)

= ( ) ( )

Since v(k)=0 except at multiples of I, where v(kI)=x(k) , therefore

( )= ( ) ( )

| ( )|

Spectra of x(n)
 ( )

0 =

Spectra of v(n)

Sampling Rate Conversion by a Rational Factor I/D

We can achieve sampling rate conversion by a rational factor I/D by first performing
interpolation by a factor I and then decimating the output of the interpolator by the factor D. The
importance of performing interpolation first and decimation second , is to preserve the desired
spectral characteristics of x(n). Since the two filters are operated at the same sampling rate IF x ,
the two filters can be replaced by a single lowpass filter with impulse response h(l) .The
v) of the combined filter is given as

,0 | | min ( , )
v) =
0 ,

x(n) Upsampler y(m)

Filter Filter Downsampler

hu(l) hd(l)

Rate F x Rate = =Fy

Rate = IFx

Method of sampling rate conversion by factor I/D

 Upsampler Low pass Downsampler
x(n) v(k) w(l) y(m)
filter

h(l)

Rate = Fx Rate = =Fy

Rate = IFx =FV

Method of sampling rate conversion by factor I/D

v= 2 F/Fv= 2 F/IFx x/ I

, = 0, ± , ±2 , …
v(l) =
0 ,

w(l) = ( ) ( )

= ( ) ( )

y(m) = w(mD)

= ( ) ( )

Filter Design and Implementation for Sampling-rate Conversion

It includes various structures like:

Direct form FIR filter

Polyphase filter
Time-variant filter

Direct Form FIR Filter Structure

This is the simplest realization with system function

H(z) = ( )

where h(k) is the unit sample response of the FIR filter.

Downsampler

h(0)
y(m)
z-1
x(n) Upsampler

h(1)

z-1

h(2)

h(M-2)

z-1
h(M-1)

Although the direct form FIR realization is simple ,it is also very inefficient. The inefficiency results
from the fact that the upsampling introduces I-1 samples between successive points of the input signal. If
I is large, most of the signal components in the FIR filter are zero. Consequently, most of the
multiplications and additions result in zeros. Furthermore, the downsampling process at the output of the
filter implies that only one out of every D output samples is required at the output of the filter.

Polyphase Filter Structure

The polyphase structure is based on the fact that any system function can be split as

H(z) = .......+ h(0) + h(M) z-M +…………

 ……+ h(1)z-1 + h(M+1) z -(M+1) +……

……+ h(M-1) z-(M-1) + h(2M-1) z-(2M-1) +………

If we next factor out the term z-(i-1) at the ith row, we obtain

H(z) = [……+ h(0) + h(M) z -M + ……]

+ z-1[…..+ h(1) + h(M+1) z-M +…….]

+ z-(M-1) […..+ h(M-1) + h(2M-1) z-M +………..]

H(z) = ( )
This implies

where Pi(z) = ( + )

P 0(z3)
x(n) y(n)

z-1

P1(z3)

z-1

P2(z3)

Block diagram of polyphase filter structure for M = 3

The output sequence for M = 3 will be:

Y(z) = H(z)X(z)

=P0(z3) X(z) + z -1 P1(z 3) X(z) + z-2 P2(z 3) X(z)

Implementation of a decimation system using a polyphase structure

P 0(z3)
x(n) y(m)
z-1

P 1(z3)

z-1

P2(z3)

Implementation of a interpolation system using a polyphase structure

P0(z3)

x(n) y(m)

z-1

P1(z3)

z-1

P2(z3)
Multistage Implementation of Sampling-rate Conversion
Let us consider interpolation by a factor I>>1 and let us assume that I can be factored into a
product of positive integers as

I=

x(n) I1Fx I 1l 2F x y(m)

1 h 1(n) 2 h2(n) L hL(n)
Fx Fy = IFx

Stage 1 Stage 2 Stage L

Multistage implementation of interpolation by a factor I

Interpolation by a factor I can be accomplished by cascading L stages of interpolation and

filtering. The filter in each of the interpolators eliminates the images introduced by the
upsampling process in the corresponding interpolator.

In a similar manner, decimation by a factor D , where D may be factored into a product of

positive integers as

D=

can be implemented as a cascade of J stages of filtering and decimation.

x(n) h1(n) 1 h2(n) 2 hJ(n) J

y(m)

Fx Stage 1 Stage 2 Stage J

Multistage implementation of decimation by a factor D

Sampling Rate Conversion of Bandpass Signals

Any band pass signal has an equivalent low pass representation , obtained by a simple frequency
translation of the band pass signal.

An analog band pass signal can be represented as

x(t) = A(t)cos[2 + ( )]

= A(t) cos ( ) cos(2 ) A(t) sin ( ) sin (2 )

= uc(t) cos(2 ) u s(t) sin(2 )

= Re[xl(t) ]

where uc(t) = A(t) cos ( )

u s(t) = A(t) sin ( )

xl(t) = u c(t) + j us(t)

A(t) is called the amplitude or envelope of the signal , ( ) is the phase , uc(t) and us(t) are called
the quadrature components of the signal.

Physically the translation of x(t) to low pass involves multiplying x(t) by the quadrature carriers
cos(2 ) and sin (2 ) and then low pass filtering the two products to eliminate the
frequency components generated around the frequency 2 .

Low pass
Filter uc(n)

cos (2 )

x(n)
Oscillator
Band pass signal

sin (2 )

us(n)
Low pass
Filter

(Conversion of a band pass signal to low pass)

The mathematical equivalence between the band pass signal and its low pass representation
provides one method for altering the sampling rate of the signal .

x(n) Frequency uc(n) Filter

Translation
Band pass signal us(n)

Sampling rate conversion of a band pass signal

The filter in the above diagram has the frequency response characteristics :

,0 | | min ( , )
0 ,

where is the bandwidth of the discrete – time band pass signal ( ).

Applications of Multirate Signal Processing

1. Design of Phase Shifters:
A delay of (k/I)Tx can be achieved by sample rate conversion method without introducing
any significant distortion in the signal. The sampling rate is increased by a factor I using
standard interpolator. The lowpass filter eliminates the images in the spectrum of the
interpolated signal , and its output is delayed by k samples at the sampling rate IF x . The
delayed signal is decimated by a factor D = I. Thus we have achieved the desired delay of
(k/I)Tx.

x(n) Delay by k y(n)

Lowpass
filter samples
Fx IF x IFx IFx Fx

x(n) upsampled= x(n/I)

x(n/I) delayed by k samples= x((n-k)/I)
x((n-k)/I) downsampled = x(n-(k/I))

If [ ( )] = ( ) ,Then [ ( )] = ( )
 /
So [ ( / )] = ( )
=

2. Interfacing of Digital Systems with Different Sampling Rates :

Let us consider interfacing the two systems with independent clocks. The output of the
system A at rate Fx is fed to an interpolator which increases the sampling rate by I. The
output of the interpolator is fed at the rate lF x to a digital sample-and-hold which serves
as the interface to system B at the high sampling rate IFx . Signals from the digital
sample-and-hold are read out into system B at the clock rate DF y of system B . Thus the
output rate from the sample-and-hold is not synchronized with the input rate.

x(n) IF x Digital IFy x(m)

System
System sample-
A Interpolation and-hold Decimator B

Fx DFy Fy

Clock Clock
IF x
A B

3. Implementation of Narrowband Lowpass Filters:

A lowpass , linear-phase FIR filter may be more efficiently implemented in a multistage
decimator-interpolator configuration. To be more specific , we can employ a multistage
implementation of a decimator of size D , followed by a multistage implementation of an
interpolator of size I , where I = D.

4. Implementation of Digital Filter Banks:

Filter banks are generally categorized as two types , analysis filter banks and synthesis
filter banks . An analysis filter bank consists of a set of filters , with system function
{Hk(z)} , arranged in parallel. The frequency response characteristics of this filter bank
split the signal into a corresponding number of subbands. On the other hand, a synthesis
filter consists of a set of filters with system function {G k(z)}, with corresponding inputs
{yk(n)}. The outputs of filters are summed to form the synthesized signal {x(n)}.
 y1(n)
H0(z)

H1(z) y2(n)
x(n)
.
.
.
yN-1(n)
HN-1(z)

(Analysis filter bank)

y1(n) G0(z)

y2(n) G1(z)
x(n)
.
.
.

yN-1(n) GN-1(z)

(Synthesis filter bank)

 H0 1 2 N-1

… ......
( )
0 2

(Frequency response characteristics of N filters)

5. Subband Filters :
Subband coding is a method where the signal is subdivided into several frequency bands
and each band is digitally encoded separately.
Let us assume a signal with sampling rate Fs . The first frequency subdivision splits
the signal spectrum into two equal-width segments , a lowpass signal (0 F F s/4) and a
highpass signal (Fs/4 Fs/2). . The second frequency subdivision splits the
lowpass signal from the first stage into two equal bands , a lowpass signal (0 F Fs/8)
and a highpass signal (Fs/8 Fs/4). Finally, the third subdivision splits the
lowpass signal from the second stage into two equal-bandwidth signals. Thus the signal is
subdivided into four frequency bands , covering three octaves.
 Module 2.
Linear Prediction and
Optimum Linear Filters
 Innovations Representation of a stationary random process:
A wide-sense stationary random process can be represented as the output of a causal and casually
invertible linear system excited by a white noise process. Let the wide-sense stationary random
xx xx (f),

f < . Then it can be shown that :

( )= ( ) ( )

Or on the unit circle,

( )= ( )

H(f) hence represents a filter which when excited by white noise w(n) of power spectral density
2 2 2
w , gives an output {x(n)} with power spectral density w . The random process x(n) is
generated by passing white noise input sequence w(n) through a linear causal filter H(z) :

( )
H(z) = = | |> <1
( )

where the polynomials B(z) and A(z) have roots that fall inside the unit circle in the z-plane.
{bk} and {ak} are filter coefficients that determine the location of the zeros and poles of H(z) ,
respectively. Thus the output x(n) is related to the input w(n) by the difference equation

x(n) + ( )= ( )

The representation of the stationary random process be {x(n)}, as the output of an IIR filter with
transfer function H(z) and excited by a white noise process {w(n)} is called as Wold
representation. The stationary random process {x(n)} can be transformed into white noise
process by passing {x(n)} through a linear filter with system function 1/H(z) . This filter is called
a noise whitening filter. Its output, denoted as {w(n)} is called the innovations process
associated with the stationary random process {x(n)}.

Linear
w(n) causal filter x(n) = ( ) ( )
H(z)
White noise

Filter for generating the random process x(n) from white noise
 Linear
x(n) causal filter w(n)
1/H(z)
white noise

Noise Whitening Filter

Rational Power Spectra

When the power spectral density of a stationary random process is a rational function given by:

B(z)B(z )
( )=
A(z)A(z )

and B(z) and A(z) have roots that fall inside the unit circle, then the filter H(z) for generation of
{x(n)} is also rational expressed as :

( )
H(z) = = | |> <1
( )

The H(z) is causal, stable and minimum phase linear filter. The noise whitening filter 1/H(z) is
also causal, stable and minimum phase linear filter. We can have three special cases on the basis
of coesfficients ak’s and bk’s.

Autoregressive (AR) process:

b0 = 1 , bk = 0 , k > 0. In this case , the linear filter H(z) = 1/A(z) is an all-pole filter and the
difference equation for the input-output relationship is

x(n) + ( ) = w(n)

The noise-whitening filter for generating the innovations process is an all-zero filter.

Moving average (MA) process:

ak = 0 , k . In this case , the linear filter H(z) = B(z) is an all-zero filter and the difference
equation for the input-output relationship is

x(n) = ( )
The noise-whitening filter for the MA process is an all-pole filter.

Autoregressive, moving average (ARMA) process:

In this case , the linear filter H(z) = B(z)/A(z) has both finite poles and zeros in the z-plane and
the corresponding difference equation for the input-output relationship is

x(n) + ( )= ( )

The inverse system for generating the innovations process from x(n) is a pole-zero system of the
form 1/H(z) = A(z)/B(z).

Forward and Backward Linear Prediction :

Forward Linear Prediction:

In forward linear prediction , a future value of stationary random process is predicted from
observation of past values of the process . The one-step forward linear predictor forms the
prediction of the value x(n) by weighted linear combination of the past values x(n-1) , x(n-2) , . .
. . , x(n-p) . Hence the linearly predicted value of x(n) is

( )= ( ) ( )

where the {- ( )} represent the weights in the linear combination . These weights are called
the prediction coefficients of the one-step forward linear predictor of order p .

The difference between the value x(n) and the predicted value x(n) is called the forward
prediction error , denoted as fp(n):

fp(n) = ( ) ( )

= ( )+ ( ) ( )

The direct form FIR filter realization to find the forward prediction error can hence be written as:

( )= ( )

where by definition a p(0)=1. The equivalent direct for FIR filter realization structure is drawn
below :
 z-1 z-1 z-1 z-1

fp(n)

Direct Form Structure Of Prediction Error Filter

fp(n)
x(n)

x(n-1) Forward ( )
z -1 linear
predictor

Forward linear prediction

Backward Linear Prediction

Suppose we have the data sequence x(n) , x(n-1) , . . . . , x(n - p + 1) from a stationary random
process and we wish to predict the value x(n-p) of the process . In this case a one-step backward
linear predictor of order p is employed. Hence

( )= ( ) ( )

The difference between the value x(n-p) and the estimate ( ) is called the backward
prediction error , denoted as gp(n) :

gp(n) = x(n - p) + ( ) ( )

= ( ) ( ), b p(p) = 1
 First Stage Second Stage Pth Stage fp(n)
x(n)
gp(n)

fm-1(n) + fm(n)

Km *

Km
gm-1(n) +
z-1 gm(n)
(n)
p-stage Lattice filter for forward and backward prediction

Normal Equations :

The forward prediction error is given as

fp(n) = ( ) ( )

The corresponding z-transform relationship is

Fp(z) = Ap(z) X(z)

( ) ( )
=> Ap(z) = =
( ) ( )

The mean square value of forward linear prediction error fp(n) is

= [ ( ) ]

= (0) + 2 [ ( ) (k)] + () ( ) (l k)
 is a quadratic function of the predictor coefficients and its minimum leads to the set of linear
equations

(l) = ( ) (l k) , l=1, 2, . . . , p

These are called the normal equations for the coefficients of the linear predictor.

The minimum mean-square prediction error is simply

min[ ] = (0) + ( ) (l k)

Solutions of the Normal Equations:

The normal equations may be expressed in the compact form

( ) (l k) = 0 l = 1, 2, … . , p

(0) = 1

If we augment the minimum MSE expression with the above equation , we get

, =0
( ) (l k) =
0 , = 1, 2, … . , p

The Levinson-Durbin Algorithm:

This algorithm exploits the symmetry in the autocorrelation matrix

(0) (1) … ( 1)
(1) (0) … ( 2)
=
…
(p 1) (p 2) … (0)

The solution to the first-order predictor is

( )
(1) = ( )
The next step is to solve for the coefficients { (1) , (2)} of the second order predictor. The
two equations obtained from the normal equation are:

(1) (0) + (2) (1) = (1)

(1) (1) + (2) (0) = (2)

By using above equation and the expression for (1) , we obtain the solution

(2) + (1) (1)

(2) =
(0)[1 | (1)| ]

Properties of the linear prediction-error filters :

1. Minimum-phase property of the forward prediction-error filter:

If a random process consists of a mixture of a continuous power spectral density and
a discrete spectrum , the prediction-error filter must have all its roots inside the unit
circle.

2. Maximum-phase property of the backward prediction-error filter:

The system function for the backward prediction-error filter of order p is
Bp(z) = z-p (z-1)
where Ap(z) is the system function for the forward prediction filter.
Consequently, the roots of Bp(z) are the reciprocals of the roots of the forward
prediction-error filter with system function Ap(z). Hence if Ap(z) is minimum phase ,
then Bp(z) is maximum phase.

3. Orthogonality of the backward prediction errors:

The backward prediction errors {gm(k)} from different stages in the FIR lattice filter
are orthogonal. That is ,
0 ,0 1
[ ( ) ( )] =
, =
4. Whitening property:
The prediction-error filter attempts to remove the correlation among the signal
samples of the input process. The response of the prediction-error filter is a white
noise sequence . Thus it whitens the input random process and is called a whitening
filter.
AR Lattice Structure:
The difference equation for AR process is given by:

y(n) = x(n) + ( )

If fp(n) is the forward prediction error , gp(n) is the backward prediction error and {Km} are the
reflection coefficients then we have:

x(n) = fp(n)
fm-1(n) = fm(n) – Km g m-1(n-1) , m = p, p-1, . . . , 1
gm(n) = ( ) + g m-1(n-1)
y(n) = f0(n) = g0(n)

x(n)=fp(n) fp-1(n) f2(n) f1(n) f0(n)=y(n)

Input output

-Kp -K2 -K1

z-1 z-1 z-1

gp(n) g2(n) g1(n) g0(n)

Lattice structure of an all pole system

ARMA Lattice-Ladder Filters:

We consider an all-pole lattice filter with coefficients Km , 1
by taking as the output a weighted linear combination of {gm(n)}.The result is a pole-zero filter
that has the lattice-ladder structure. Its output is

y(n) = ( )

where { } are the parameters that determine the zeros of the system.
x(n) = fp(n) fp-1(n) fp-2(n) f1(n) f0(n)
Stage Stage Stage
gp(n) p gp-1(n) p-1 gp-2(n) g1(n) 1 g0(n)

output

Pole-zero system

fm(n) fm-1(n)

- Km

gm(n) gm-1(n)
z-1

mth stage lattice

Wiener filter for filtering and prediction:

It deals with the problem of estimating a desired signal {s(n)} in the presence of an undesired
additive noise disturbance {w(n)}. The estimator is constrained to be a linear filter with impulse
response {h(n)}, designed so that the output approximates some specified desired signal
sequence {d(n)}.
 d(n)

s(n) x(n) y(n) e(n)

Optimum
linear filter

Noise w(n)

Weiner Filter

The criterion selected for optimizing the filter impulse response {h(n)} is the minimization of the
mean-square error. The optimum linear filter, in the sense of minimum mean-square error
(MMSE), is called a Wiener filter.

Orthogonality Principle in Linear Mean-square Estimation

The mean- M is a minimum if the filter coefficients {h(k)} are selected such that the
error is orthogonal to each of the data points in estimate,

E[e(n) x*(n-l)] = 0, l = 0, 1, . . . , M-1

where
e(n) = d(n) – ( ) ( )

Since the MSE is minimized by selecting the filter coefficients to satisfy the orthogonality
principle, the residual minimum MSE is simply

MMSE M = E[e(n) d*(n)]

 Module 3.
Power spectrum estimation
 Power Spectrum Estimation
The power spectrum estimation deals with the estimation of the spectral characteristics of signals
characterized as random processes. Many of the phenomena that occur in nature are best
characterized statistically in terms of averages. For example, meteorological phenomena such as
the fluctuations in air temperature and pressure are best characterized statistically as random
processes.

Due to random fluctuations in such signals, we must adopt a statistical view point, which deals
with the average characteristics of random signals. In particular, the autocorrelation function of
random process is the appropriate statistical average that we will use for characterizing random
signals in the time domain, and the Fourier transform of the autocorrelation function, which
yields the power density spectrum, provides the transform from the time domain to frequency
domain.

Estimation of Spectra from Finite-Duration Observation of Signals:

Computation of the Energy Density Spectrum:
Let ( ) be a finite-duration sequence obtained by sampling a continuous-time signal ( ) at
some uniform sampling rate .

If x(t) is a finite-energy signal, that is

= | ( )| <

then its Fourier transform exists and is given as

( )= ( )

From Parseval’s theorem we have

= | ( )| = | ( )|

The quantity | ( )| represents the distribution of signal energy as a function of frequency,

and hence it is called the energy density spectrum of the signal, that is,

( )=| ( )|
Thus the total energy in the signal is simply the integral of ( ) over all [i.e., the total area
under ( )].

( ) can be viewed as the Fourier transform of another function, ( ) , called the

autocorrelation function of the finite-energy signal ( ) , defined as

( )= ( ) ( + )

It follows that

( ) = ( )=| ( )|

Hence ( ) and ( ) are a Fourier transform pair.

Now the Fourier transform (voltage spectrum) of ( ):

( )= ( )

or, equivalently,

( )= ( )

, = ( )
( )

where = is the normalized frequency variable.

In absence of aliasing, within the fundamental range | | , we have

= ( ) ,| |
( ) 2

Hence the voltage spectrum of the sampled signal is identical to the voltage spectrum of the
analog signal. As a consequence, the energy density spectrum of the sampled signal is
 = = | ( )|
( ) ( )

The autocorrelation of the sampled signal ( ) is defined as:

( )= ( ) ( + )

its Fourier transform (Wiener-Khintchine theorem):

( )= ( )

Hence the energy density spectrum can be obtained by the Fourier transform of the
autocorrelation of the sequence { ( )} , that is,

( ) = | ( )|

= ( )

Estimation of the Autocorrelation and Power Spectrum of Random Signals:

The Periodogram:
The finite-energy signals possess a Fourier transform and are characterized in the spectral
domain by their energy density spectrum. On the other hand, the important class of signals
characterized as stationary random processes do not have finite energy and hence do not posses a
Fourier transform. Such signals have finite average power and hence are characterized by a
power density spectrum.

If ( ) is a stationary random process, its autocorrelation function is

( ) = [ ( ) ( + )]

where [. ] denotes the statistical average. Then by Wiener-Khintchine theorem, the power
density spectrum of the stationary random process is the Fourier transform of the autocorrelation
function:
 ( )= ( )

But we do not know the true autocorrelation function ( ) and as a consequence, we cannot
compute the Fourier transform in (1.13) to obtain ( ). On the other hand, from a single
realization of the random process we can compute the time-average autocorrelation function:

1
( )= () ( + )
2

where 2 is the observation interval.

The Fourier transform of ( ) provides an estimate ( ) of the power density spectrum, that
is,

( )= ( )

1
= ( ) ( + )
2

1
= ()
2

The actual power density spectrum is the expected value of ( ) in the limit as ,

( ) = lim [ ( )]

1
= lim ( )
2

The estimate ( ) can also be expressed as

1 1
( )= ( ) = | ( )|

where ( ) is the Fourier transform of the finite duration sequence ( ) , 0 1. This

form of the power density spectrum estimate is called the periodogram.
Nonparametric Methods for Power Spectrum Estimation:
The nonparametric methods make no assumption about how the data were generated.

The Bartlett Method: Averaging Periodograms:

It reduces the variance in the periodogram. The -point sequence is subdivided into
nonoverlapping segments, where each segment has length . This results in the data segments

= 0, 1, … , 1
( )= ( + ),
= 0, 1, … , 1

For each segment, we compute the periodogram

( ) 1
( )= ( ) , = 0, 1, … , 1

The Bartlett power spectrum estimate obtained by averaging the periodograms for the
segments is
1 ()
( )= ( )

The mean value:

1 ()
[ ( )] = [ ( )]

()
= ( )
The expected value of single periodogram:

() | |
( ) = 1 ( )
( )

1 sin ( )
= ( )
sin ( )
where
1 sin
( )=
sin

is the frequency characteristic of the Bartlett window

 | |
( )= 1 , | | 1
0,

The variance of the Bartlett estimate:

1 ()
[ ( )] = [ ( )]

1 ()
= ( )

But
sin 2
[ ( )] = ( ) 1+
sin 2

Putting the above value in (2.8), we get

1 sin 2
[ ( )] = ( ) 1+
sin 2

Therefore, the variance of the Bartlett power spectrum estimate has been reduced by the factor
.

The Blackman and Tukey Method: Smoothing the Periodogram:

In this method the sample autocorrelation sequence is windowed first and then Fourier
transformed to yield the estimate of the power spectrum.

The Blackman-Tukey estimate is

( )= ( ) ( )
( )

where the window function ( ) has length 2 1 and is zero for | | .

The frequency domain equivalent expression:

 /
( )= ( ) ( )
/

where ( ) is the periodogram.

The expected value of the Blackman-Tukey power spectrum estimate:

/
[ ( )] = [ ( )] ( )
/

where
/
[ ( )] = ( ) ( )
/

and ( ) is the Fourier transform of the Bartlett window. Substitution of (2.14) into (2.13)
yields

/ /
[ ( )] = ( ) ( ) ( )
/ /

In time domain we have:

[ ( )] = [ ( )] ( )
( )

= ( ) ( ) ( )
( )

where the Bartlett window is,

| |
( )= 1 , | |<
0 ,

The variance of the Blackman-Tukey power spectrum estimate is

[ ( )] = {[ ( )] } { [ ( )]}
Assuming ( ) is narrow compared to the true power spectrum ( )

/
1
[ ( )] ( ) ( )
/

1
( ) ( )
( )

Parametric Methods for Power Spectrum Estimation:

Parametric methods avoid the problem of spectral leakage and provide better frequency
resolution than do the nonparametric methods. Parametric methods also eliminate the need for
window functions.

The Yuke-Walker Method for the AR Model Parameters:

This method is used to estimate the autocorrelation from the data and use the estimates to solve
for the AR model parameters.
The autocorrelation estimate is given by

1
( )= ( ) ( + ), 0

The corresponding power spectrum estimate is

( )=
1+ ( )

where ( ) are estimates of the AR parameters and

= = (0) [1 | ( )| ]

is the estimated minimum mean-square value for the pth-order predictor.

The Burg method for the AR Model Parameters:
Suppose that we are given the data ( ), = 0,1, … , 1. The forward and backward linear
prediction estimates of order are given as:

( )= ( ) ( )

( )= ( ) ( + )

and the corresponding forward and backward errors ( ) and ( ) defined as

( )= ( ) ( ) and ( )= ( ) ( ) where ( ), 0 1,
= 1, 2, … , , are prediction coefficients.

The least-squares error is

= [| ( )| + | ( )| ]

The power spectrum estimate is

( )=
1+ ( )

where is the total least-squares error.

Advantages of the Burg method:

results in high frequency resolution
yields a stable AR model
computationally efficient

Disadvantages of the Burg method:

exhibits spectral line splitting at high signal-to-noise ratios
exhibits a sensitivity to the initial phase of sinusoid for sinusoidal signals in noise
resulting in a phase-dependent frequency bias.
The MA Model for Power Spectrum Estimation :

The parameters of MA model are related to the statistical autocorrelation ( ) by

, 0
( )=
0, >
( ), <0
However,

( ) ( )= ( )=

where coefficients { } are related to the MA parameters by the expression

| |

= , | |

Clearly, then,
, | |
( )=
0, | |>

and the power spectrum for MA process is

( )= ( )

It is apparent from these expressions that we do not have to solve for the MA parameters { }
to estimate the power spectrum. The estimates of the autocorrelation ( ) for | | suffice.
From such estimates we compute the estimated MA power spectrum , given as

P ( )= ( )
The ARMA Model for Power Spectrum Estimation :

An ARMA model provides us with an opportunity to improve on the AR spectrum estimate by

using few model parameters. The ARMA model is particularly appropriate when the signal has
been corrupted by noise.

The sequence ( ) can be filtered by an FIR filter to yield the sequence

( )= ( )+ ( ), = 0, 1, … , 1

The filtered sequence ( ) for 1 is used to form the estimated correlation

sequences ( ) , from which we obtain the MA spectrum

P ( )= ( )

Finally, the estimated ARMA power spectrum is

P ( )
( )=
1+
 Module 4.
Adaptive Signal Processing
ADAPTIVE NOISE CANCELLATION
In speech communication from a noisy acoustic environment such as a moving car or train,
or over a noisy telephone channel, the speech signal is observed in an additive random noise.

In signal measurement systems the information-bearing signal is often contaminated by noise

from its surrounding environment. The noisy observation, y(m), can be modelled as
y(m) = x(m)+n(m)
where x(m) and n(m) are the signal and the noise, and m is the discrete-time index. In
some situations, for example when using a mobile telephone in a moving car, or when using
a radio communication device in an aircraft cockpit, it may be possible to measure and
estimate the instantaneous amplitude of the ambient noise using a directional microphone.
The signal, x(m), may then be recovered by subtraction of an estimate of the noise from the
noisy signal.

(Configuration of a two-microphone adaptive noise canceller.)

Figure shows a two-input adaptive noise cancellation system for enhancement of noisy speech.
In this system a directional microphone takes as input the noisy signal x(m)+n(m), and a second

-simplified model of the effects of

propagation of the noise to different positions in the space where the microphones are placed.
The noise from the second microphone is processed by an adaptive digital filter to make it equal
to the noise contaminating the speech signal, and then subtracted from the noisy signal to cancel
out the noise. The adaptive noise canceller is more effective in cancelling out the low-frequency
part of the noise, but generally suffers from the non-stationary character of the signals, and from
the over-simplified assumption that a linear filter can model the diffusion and propagation of the
noise sound in the space.
ADAPTIVE NOISE REDUCTION
In many applications, for example at the receiver of a telecommunication system, there is no
access to the instantaneous value of the contaminating noise, and only the noisy signal is
available. In such cases the noise cannot be cancelled out, but it may be reduced, in an average
sense, using the statistics of the signal and the noise process.

A frequency-domain Wiener filter for reducing additive noise

Figure shows a bank of Wiener filters for reducing additive noise when only the noisy signal is
available. The filter bank coefficients attenuate each noisy signal frequency in inverse proportion
to
the signal-to-noise ratio at that frequency. The Wiener filter bank coefficients, are calculated
from estimates of the power spectra of the signal and the noise processes.
BLIND CHANNEL EQUALISATION
Channel equalisation is the recovery of a signal distorted in transmission through a
communication channel with a nonflat magnitude or a nonlinear phase response. When the
channel response is unknown, the process of signal recovery is called ‘blind equalisation’. Blind
equalisation has a wide range of applications, for example in digital telecommunications for
removal of inter-symbol interference due to nonideal channel and multipath propagation, in
speech recognition for removal of the effects of the microphones and communication channels,
in correction of distorted images, in analysis of seismic data and in de-reverberation of acoustic
gramophone recordings. In practice, blind equalisation is feasible only if some useful statistics of
the channel input are available. The success of a blind equalisation method depends on how
much is known about the characteristics of the input signal and how useful this knowledge can
be in the channel identification and equalisation process. Figure 1.6 illustrates the configuration
of a decision-directed equaliser. This blind channel equaliser is composed of two distinct
sections: an adaptive equaliser that removes a large part of the channel distortion, followed by a
nonlinear decision device for an improved estimate of the channel input. The output of the
decision device is the final estimate of the channel input, and it is used as the desired signal to
direct the equaliser adaptation process.

Configuration of a decision-directed blind channel equaliser.

Adaptive Linear Combiner

Adaptive linear combiner showing the combiner and the adaption process. k = sample number,

ower box=adaption algorithm.

Adaptive linear combiner, compact representation. k = sample number, n=input variable

The adaptive linear combiner (ALC) resembles the adaptive tapped delay line FIR filter except
that there is no assumed relationship between the X values. If the X values were from the outputs
of a tapped delay line, then the combination of tapped delay line and ALC would comprise an
adaptive filter. However, the X values could be the values of an array of pixels. Or they could be
the outputs of multiple tapped delay lines. The ALC finds use as an adaptive beam former for
arrays of hydrophones or antennas.
 where refers to the 'th weight at k'th time.

LMS algorithm

Main article: Least mean squares filter

If the variable filter has a tapped delay line FIR structure, then the LMS update algorithm is
especially simple. Typically, after each sample, the coefficients of the FIR filter are adjusted as
follows:[5](Widrow)

for

convergence factor.

The LMS algorithm does not require that the X values have any particular relationship; therefor
it can be used to adapt a linear combiner as well as an FIR filter. In this case the update formula
is written as:

The effect of the LMS algorithm is at each time, k, to make a small change in each weight. The
direction of the change is such that it would decrease the error if it had been applied at time k.

at time k. The weights making the largest contribution to the output, , are changed the most. If
the error is zero, then there should be no change in the weights. If the associated value of X is
zero, then changing the weight makes no difference, so it is not changed.

Convergence

oo large to prevent
convergence, the algorithm reaches steady state rapidly but continuously overshoots the optimum

minimize overshoot.

Widrow and Stearns state in 1985 that they have no knowledge of a proof that the LMS
algorithm will converge in all cases.
However under certain assumptions about stationarity and independence it can be shown that the
algorithm will converge if

where

= sum of all input power

is the RMS value of the 'th input

In the case of the tapped delay line filter, each input has the same RMS value because they are
simply the same values delayed. In this case the total power is

where
is the RMS value of , the input stream.

This leads to a normalized LMS algorithm:

in which case the convergence criteria becomes:

.

Least mean squares filter

Least mean squares (LMS) algorithms are a class of adaptive filter used to mimic a desired
filter by finding the filter coefficients that relate to producing the least mean squares of the error
signal (difference between the desired and the actual signal). It is a stochastic gradient descent
method in that the filter is only adapted based on the error at the current time. It was invented in
1960 by Stanford University professor Bernard Widrow and his first Ph.D. student, Ted Hoff.
Problem formulation

Relationship to the least mean squares filter

The realization of the causal Wiener filter looks a lot like the solution to the least squares
estimate, except in the signal processing domain. The least squares solution, for input matrix
and output vector is

The FIR least mean squares filter is related to the Wiener filter, but minimizing the error
criterion of the former does not rely on cross-correlations or auto-correlations. Its solution
converges to the Wiener filter solution. Most linear adaptive filtering problems can be
formulated using the block diagram above. That is, an unknown system is to be identified
and the adaptive filter attempts to adapt the filter to make it as close as possible to ,
while using only observable signals , and ; but , and are not
directly observable. Its solution is closely related to the Wiener filter.

Definition of symbols

is the number of the current input sample

is the number of filter taps

(Hermitian transpose or conjugate transpose)

 estimated filter; interpret as the estimation of the filter coefficients after n samples

Idea

The basic idea behind LMS filter is to approach the optimum filter weights , by
updating the filter weights in a manner to converge to the optimum filter weight. The algorithm
starts by assuming a small weights (zero in most cases), and at each step, by finding the gradient
of the mean square error, the weights are updated. That is, if the MSE-gradient is positive, it
implies, the error would keep increasing positively, if the same weight is used for further
iterations, which means we need to reduce the weights. In the same way, if the gradient is
negative, we need to increase the weights. So, the basic weight update equation is :

where represents the mean-square error. The negative sign indicates that, we need to change
the weights in a direction opposite to that of the gradient slope.

The mean-square error, as a function of filter weights is a quadratic function which means it has
only one extrema, that minimises the mean-square error, which is the optimal weight. The LMS
thus, approaches towards this optimal weights by ascending/descending down the mean-square-
error vs filter weight curve.

Derivation

The idea behind LMS filters is to use steepest descent to find filter weights which
minimize a cost function. We start by defining the cost function as

where is the error at the current sample n and denotes the expected value.
This cost function ( ) is the mean square error, and it is minimized by the LMS. This is
where the LMS gets its name. Applying steepest descent means to take the partial derivatives
with respect to the individual entries of the filter coefficient (weight) vector

where is the gradient operator

Now, is a vector which points towards the steepest ascent of the cost function. To find
the minimum of the cost function we need to take a step in the opposite direction of .
To express that in mathematical terms

where is the step size(adaptation constant). That means we have found a sequential update
algorithm which minimizes the cost function. Unfortunately, this algorithm is not realizable until
we know .

Generally, the expectation above is not computed. Instead, to run the LMS in an online (updating
after each new sample is received) environment, we use an instantaneous estimate of that
expectation. See below.

Simplifications

For most systems the expectation function must be approximated. This can
be done with the following unbiased estimator

where indicates the number of samples we use for that estimate. The simplest case is
For that simple case the update algorithm follows as

Indeed this constitutes the update algorithm for the LMS filter.

LMS algorithm summary

The LMS algorithm for a th order algorithm can be summarized as

Parameters: filter order

step size

Initialisation:

Computation: For

Convergence and stability in the mean

As the LMS algorithm does not use the exact values of the expectations, the weights would never
reach the optimal weights in the absolute sense, but a convergence is possible in mean. That is,
even though the weights may change by small amounts, it changes about the optimal weights.
However, if the variance with which the weights change, is large, convergence in mean would be
misleading. This problem may occur, if the value of step-size is not chosen properly.

If is chosen to be large, the amount with which the weights change depends heavily on the
gradient estimate, and so the weights may change by a large value so that gradient which was
negative at the first instant may now become positive. And at the second instant, the weight may
change in the opposite direction by a large amount because of the negative gradient and would
thus keep oscillating with a large variance about the optimal weights. On the other hand if is
chosen to be too small, time to converge to the optimal weights will be too large.

Thus, an upper bound on is needed which is given as

where is the greatest eigenvalue of the autocorrelation matrix .
If this condition is not fulfilled, the algorithm becomes unstable and diverges.

Maximum convergence speed is achieved when

where is the smallest eigenvalue of R. Given that is less than or equal to this optimum,
the convergence speed is determined by , with a larger value yielding faster convergence.
This means that faster convergence can be achieved when is close to , that is, the
maximum achievable convergence speed depends on the eigenvalue spread of .

A white noise signal has autocorrelation matrix where is the variance of the
signal. In this case all eigenvalues are equal, and the eigenvalue spread is the minimum over all
possible matrices. The common interpretation of this result is therefore that the LMS converges
quickly for white input signals, and slowly for colored input signals, such as processes with low-
pass or high-pass characteristics.

It is important to note that the above upperbound on only enforces stability in the mean, but the
coefficients of can still grow infinitely large, i.e. divergence of the coefficients is still
possible. A more practical bound is

where denotes the trace of . This bound guarantees that the coefficients of do not
diverge (in practice, the value of should not be chosen close to this upper bound, since it is
somewhat optimistic due to approximations and assumptions made in the derivation of the
bound).

Normalised least mean squares filter (NLMS)

The main drawback of the "pure" LMS algorithm is that it is sensitive to the scaling of its input
. This makes it very hard (if not impossible) to choose a learning rate that guarantees
stability of the algorithm (Haykin 2002). The Normalised least mean squares filter (NLMS) is a
variant of the LMS algorithm that solves this problem by normalising with the power of the
input. The NLMS algorithm can be summarised as:
Parameters: filter order

step size

Initialization:

Computation: For

Optimal learning rate

It can be shown that if there is no interference ( ), then the optimal learning rate for
the NLMS algorithm is

and is independent of the input and the real (unknown) impulse response . In the
general case with interference ( ), the optimal learning rate is

The results above assume that the signals and are uncorrelated to each other, which is
generally the case in practice.

Proof

Let the filter misalignment be defined as , we can derive the

expected misalignment for the next sample as:
Let and

Assuming independence, we have:

The optimal learning rate is found at , which leads to:

Wiener filter

In signal processing, the Wiener filter is a filter used to produce an estimate of a desired or
target random process by linear time-invariant filtering of an observed noisy process, assuming
known stationary signal and noise spectra, and additive noise. The Wiener filter minimizes the
mean square error between the estimated random process and the desired process..
Description

The goal of the Wiener filter is to compute a statistical estimate of an unknown signal using a
related signal as an input and filtering that known signal to produce the estimate as an output. For
example, the known signal might consist of an unknown signal of interest that has been
corrupted by additive noise. The Wiener filter can be used to filter out the noise from the
corrupted signal to provide an estimate of the underlying signal of interest. The Wiener filter is
based on a statistical approach, and a more statistical account of the theory is given in the
minimum mean-square error (MMSE) article.

Typical deterministic filters are designed for a desired frequency response. However, the design
of the Wiener filter takes a different approach. One is assumed to have knowledge of the spectral
properties of the original signal and the noise, and one seeks the linear time-invariant filter
whose output would come as close to the original signal as possible. Wiener filters are
characterized by the following:[1]

1. Assumption: signal and (additive) noise are stationary linear stochastic processes with
known spectral characteristics or known autocorrelation and cross-correlation
2. Requirement: the filter must be physically realizable/causal (this requirement can be
dropped, resulting in a non-causal solution)
3. Performance criterion: minimum mean-square error (MMSE)

This filter is frequently used in the process of deconvolution; for this application, see Wiener
deconvolution.

Wiener filter solutions

The Wiener filter problem has solutions for three possible cases: one where a noncausal filter is
acceptable (requiring an infinite amount of both past and future data), the case where a causal
filter is desired (using an infinite amount of past data), and the finite impulse response (FIR) case
where a finite amount of past data is used. The first case is simple to solve but is not suited for
real-time applications. Wiener's main accomplishment was solving the case where the causality
requirement is in effect, and in an appendix of Wiener's book Levinson gave the FIR solution.

Noncausal solution

Where are spectra. Provided that is optimal, then the minimum mean-square error
equation reduces to
and the solution is the inverse two-sided Laplace transform of .

Causal solution

where

consists of the causal part of (that is, that part of this fraction having a
positive time solution under the inverse Laplace transform)
is the causal component of (i.e., the inverse Laplace transform of is
non-zero only for )
is the anti-causal component of (i.e., the inverse Laplace transform of
is non-zero only for )

This general formula is complicated and deserves a more detailed explanation. To write down
the solution in a specific case, one should follow these steps:[2]

1. Start with the spectrum in rational form and factor it into causal and anti-causal
components:

where contains all the zeros and poles in the left half plane (LHP) and contains the zeroes
and poles in the right half plane (RHP). This is called the Wiener–Hopf factorization.

2. Divide by and write out the result as a partial fraction expansion.

3. Select only those terms in this expansion having poles in the LHP. Call these terms
.
4. Divide by . The result is the desired filter transfer function .

Finite impulse response Wiener filter for discrete series

Block diagram view of the FIR Wiener filter for discrete series. An input signal w[n] is
convolved with the Wiener filter g[n] and the result is compared to a reference signal s[n] to
obtain the filtering error e[n].

The causal finite impulse response (FIR) Wiener filter, instead of using some given data matrix
X and output vector Y, finds optimal tap weights by using the statistics of the input and output
signals. It populates the input matrix X with estimates of the auto-correlation of the input signal
(T) and populates the output vector Y with estimates of the cross-correlation between the output
and input signals (V).

In order to derive the coefficients of the Wiener filter, consider the signal w[n] being fed to a
Wiener filter of order N and with coefficients , . The output of the filter
is denoted x[n] which is given by the expression

The residual error is denoted e[n] and is defined as e[n] = x[n] s[n] (see the corresponding
block diagram). The Wiener filter is designed so as to minimize the mean square error (MMSE
criteria) which can be stated concisely as follows:

where denotes the expectation operator. In the general case, the coefficients may be
complex and may be derived for the case where w[n] and s[n] are complex as well. With a
complex signal, the matrix to be solved is a Hermitian Toeplitz matrix, rather than symmetric
Toeplitz matrix. For simplicity, the following considers only the case where all these quantities
are real. The mean square error (MSE) may be rewritten as:
To find the vector which minimizes the expression above, calculate its derivative
with respect to

Assuming that w[n] and s[n] are each stationary and jointly stationary, the sequences
and known respectively as the autocorrelation of w[n] and the cross-correlation
between w[n] and s[n] can be defined as follows:

The derivative of the MSE may therefore be rewritten as (notice that )

Letting the derivative be equal to zero results in

which can be rewritten in matrix form

These equations are known as the Wiener–Hopf equations. The matrix T appearing in the
equation is a symmetric Toeplitz matrix. Under suitable conditions on , these matrices are
known to be positive definite and therefore non-singular yielding a unique solution to the
determination of the Wiener filter coefficient vector, . Furthermore, there exists an
efficient algorithm to solve such Wiener–Hopf equations known as the Levinson-Durbin
algorithm so an explicit inversion of is not required.
Relationship to the least squares filter

The realization of the causal Wiener filter looks a lot like the solution to the least squares
estimate, except in the signal processing domain. The least squares solution, for input matrix
and output vector is

The FIR Wiener filter is related to the least mean squares filter, but minimizing the error
criterion of the latter does not rely on cross-correlations or auto-correlations. Its solution
converges to the Wiener filter solution.

Recursive least squares filter

The Recursive least squares (RLS) adaptive is an algorithm which recursively finds the filter
coefficients that minimize a weighted linear least squares cost function relating to the input
signals. This is in contrast to other algorithms such as the least mean squares (LMS) that aim to
reduce the mean square error. In the derivation of the RLS, the input signals are considered
deterministic, while for the LMS and similar algorithm they are considered stochastic. Compared
to most of its competitors, the RLS exhibits extremely fast convergence. However, this benefit
comes at the cost of high computational complexity.

Motivation

RLS was discovered by Gauss but lay unused or ignored until 1950 when Plackett rediscovered
the original work of Gauss from 1821. In general, the RLS can be used to solve any problem that
can be solved by adaptive filters. For example, suppose that a signal d(n) is transmitted over an
echoey, noisy channel that causes it to be received as

where represents additive noise. We will attempt to recover the desired signal by use
of a -tap FIR filter, :
where is the vector containing the
most recent samples of . Our goal is to estimate the parameters of the filter , and at each
time n we refer to the new least squares estimate by . As time evolves, we would like to
avoid completely redoing the least squares algorithm to find the new estimate for , in
terms of .

The benefit of the RLS algorithm is that there is no need to invert matrices, thereby saving
computational power. Another advantage is that it provides intuition behind such results as the
Kalman filter.

Discussion

The idea behind RLS filters is to minimize a cost function by appropriately selecting the filter
coefficients , updating the filter as new data arrives. The error signal and desired signal
are defined in the negative feedback diagram below:

The error implicitly depends on the filter coefficients through the estimate :

The weighted least squares error function —the cost function we desire to minimize—being a
function of e(n) is therefore also dependent on the filter coefficients:

where is the "forgetting factor" which gives exponentially less weight to older error
samples.
The cost function is minimized by taking the partial derivatives for all entries of the coefficient
vector and setting the results to zero

Next, replace with the definition of the error signal

Rearranging the equation yields

This form can be expressed in terms of matrices

where is the weighted sample covariance matrix for , and is the equivalent
estimate for the cross-covariance between and . Based on this expression we find the
coefficients which minimize the cost function as

This is the main result of the discussion.

Choosing

The smaller is, the smaller contribution of previous samples. This makes the filter more
sensitive to recent samples, which means more fluctuations in the filter co-efficients. The
case is referred to as the growing window RLS algorithm. In practice, is usually chosen
between 0.98 and 1.

Recursive algorithm

The discussion resulted in a single equation to determine a coefficient vector which minimizes
the cost function. In this section we want to derive a recursive solution of the form
where is a correction factor at time . We start the derivation of the recursive
algorithm by expressing the cross covariance in terms of

where is the dimensional data vector

Similarly we express in terms of by

In order to generate the coefficient vector we are interested in the inverse of the deterministic
auto-covariance matrix. For that task the Woodbury matrix identity comes in handy. With

is -by-

is -by-1

is 1-by-

is the 1-by-1 identity matrix

The Woodbury matrix identity follows

To come in line with the standard literature, we define

where the gain vector is

Before we move on, it is necessary to bring into another form

Subtracting the second term on the left side yields

With the recursive definition of the desired form follows

Now we are ready to complete the recursion. As discussed

The second step follows from the recursive definition of . Next we incorporate the
recursive definition of together with the alternate form of and get

With we arrive at the update equation

where is the a priori error. Compare this with the a posteriori

error; the error calculated after the filter is updated:

That means we found the correction factor

This intuitively satisfying result indicates that the correction factor is directly proportional to
both the error and the gain vector, which controls how much sensitivity is desired, through the
weighting factor, .

RLS algorithm summary

The RLS algorithm for a p-th order RLS filter can be summarized as
Parameters: filter order

forgetting factor

value to initialize

Initialization: ,

where is the identity matrix of rank

Computation: For

Note that the recursion for follows an Algebraic Riccati equation and thus draws parallels to
the Kalman filter.

Lattice recursive least squares filter (LRLS)

The Lattice Recursive Least Squares adaptive filter is related to the standard RLS except that it
requires fewer arithmetic operations (order N). It offers additional advantages over conventional
LMS algorithms such as faster convergence rates, modular structure, and insensitivity to
variations in eigenvalue spread of the input correlation matrix. The LRLS algorithm described is
based on a posteriori errors and includes the normalized form. The derivation is similar to the
standard RLS algorithm and is based on the definition of . In the forward prediction case,
we have with the input signal as the most up to date sample. The
backward prediction case is , where i is the index of the sample in the
past we want to predict, and the input signal is the most recent sample.

Parameter Summary

is the forward reflection coefficient

is the backward reflection coefficient

represents the instantaneous a posteriori forward prediction error

represents the instantaneous a posteriori backward prediction error

is the minimum least-squares backward prediction error

is the minimum least-squares forward prediction error

is a conversion factor between a priori and a posteriori errors

are the feedforward multiplier coefficients.

is a small positive constant that can be 0.01

LRLS Algorithm Summary

The algorithm for a LRLS filter can be summarized as

Initialization:

For i = 0,1,...,N

(if x(k) = 0 for k < 0)

 End

Computation:

For k

For i = 0,1,...,N
 Feedforward Filtering

End

End

Normalized lattice recursive least squares filter (NLRLS)

The normalized form of the LRLS has fewer recursions and variables. It can be calculated by
applying a normalization to the internal variables of the algorithm which will keep their
magnitude bounded by one. This is generally not used in real-time applications because of the
number of division and square-root operations which comes with a high computational load.

NLRLS algorithm summary

The algorithm for a NLRLS filter can be summarized as

Initializatio
n:

For i = 0,1,...,N

(if x(k) = d(k) = 0 for k < 0)

 End

Computatio
n:

For k

(Input signal energy)

(Reference signal energy)

For i = 0,1,...,N

Feedforward Filter
End

End