Adaptive Filter Notes

Chapter
Adaptive Filters
Introduction
The term adaptive filter implies changing the characteristic of a

filter in some automated fashion to obtain the best possible signal
quality in spite of changing signal/system conditions. Adaptive
filters are usually associated with the broader topic of statistical
signal processing. The operation of signal filtering by definition
implies extracting something desired from a signal containing
both desired and undesired components. With linear FIR and IIR
filters the filter output is obtained as a linear function of the
observation (signal applied) to the input. An optimum linear filter in the minimum mean square sense can be designed to extract
a signal from noise by minimizing the error signal formed by
subtracting the filtered signal from the desired signal. For noisy
signals with time varying statistics, this minimization process is
often done using an adaptive filter.
For statistically stationary inputs this solution is known as a
Wiener filter.1
+
d[n]
e [ n ] Error
Desired
Signal
x[n]
Observation
Signal
Wiener
Filter
y[ n]
MMSE
Estimate
of d[n]
1. Simon Haykin, Adaptive Filter Theory, fourth edition, Prentice Hall, 2002.
ECE 5655/4655 Real-Time DSP
91
Chapter 9 Adaptive Filters
Wiener Filter
An M tap discrete-time Wiener filter is of the form
y[n] =
M1
wm x [ n m ]
(9.1)
m=0
where the w m are referred to as the filter weights

Note: (9.1) tells us that the Wiener filter is just an M-tap
FIR filter
The quality of the filtered or estimated signal y [ n ] is determined from the error sequence e [ n ] = d [ n ] y [ n ]
The weights w m , m = 0, 1, , M 1 are chosen such that
2
E{e [n]} = E{(d[n] y[n]) }
(9.2)
is minimized, that is we obtain the minimum mean-squared

error (MMSE)
The optimal weights are found by setting
2
E { e [ n ] } = 0, m = 0, 1, , M 1
wm
(9.3)
From the orthogonality principle1 we choose the weights

such that the error e [ n ] is orthogonal to the observations
(data), i.e.,
E { x [ n k ] ( d [ n ] y [ n ] ) } = 0, k = 0, 1, , M 1
(9.4)
1. A. Papoulis, Probability, Random Variables, and Stochastic Processes,

third edition, McGraw-Hill, 1991.
92
Wiener Filter
This results in a filter that is optimum in the sense of minimum mean-square error
The resulting system of equations
M1
w m E { x [ n k ]x [ n m ] }
m=0
= E{x[n k](d[n])}
(9.5)
or
M1
w m xx [ m k ] =
m=0
xd [ k ]
(9.6)
for k = 0, 1, , M 1 are known as the Wiener-Hopf or

normal equations
Note: xx [ k ] is the autocorrelation function of x [ n ] and
xd [ k ] is the cross-correlation function between x [ n ] and
d[n]
In matrix form we can write
R xx w o = p xd
(9.7)
where R xx is the M M correlation matrix associated with

x[ n]
...
...
R xx
xx [ 0 ] . . . xx [ M 1 ]
..
=
.
xx [ M + 1 ] . . . xx [ 0 ]
(9.8)
w o is the optimum weight vector given by

93
w o = w o0 w o1 w oM 1
(9.9)
and p xd is the cross-correlation vector given by

p xd = xd [ 0 ] xd [ 1 ] xd [ 1 M ]
(9.10)
The optimal weight vector is given by

1
w o = R xx p xd
(9.11)
As a matter of practice (9.11) can be solved using sample statistics, that is we replace the true statistical auto- and crosscorrelation functions with time averages of the form
N1
1
xx [ k ] ---N
x [ n + k ]x [ n ]
n=0
1
xd [ k ] ---N
N1
x [ n + k ]d [ n ]
n=0
(9.12)
(9.13)
where N is the sample block size
Adaptive Wiener Filter

In an adaptive Wiener filter the error signal e [ n ] is fed back
to the filter weights to adjust them using a steepest-descent
algorithm
94
With respect to the weight vector w , the error e [ n ] can be

viewed as an M dimensional error surface, that due to the
squared error criterion, is convex cup (a bowl shape)
Mean Square Error
Optimum at [1,1]
w0
w1
Error Surface for M = 2

The filter decorrelates the output error e [ n ] so that signals in
common to both d [ n ] and x [ n ] in a correlation sense are
cancelled
95
A block diagram of this adaptive Wiener (FIR) filter is shown

below
+
d[n]
Desired
Signal
x[n]
Observation
e[ n]
Error
Signal
y[ n]
MMSE
Estimate
of d[n]
Wiener
Filter
Adaptive
Algorithm
Adaptive Wiener FIlter

Least-Mean-Square Adaptation
Ideally the optimal weight solution can be obtained by applying the steepest descent method to the error surface, but since
the true gradient cannot be determined, we use a stochastic
gradient, which is simply the instantaneous estimate of R xx
and p xd from the available data, e.g.,
xx = x [ n ]x T [ n ]
R
(9.14)
p xd = x [ n ]d [ n ]
(9.15)
where
x[n] = x[n] x[n 1] x[n M + 1]
(9.16)
A practical implementation involves estimating the gradient

from the available data using the least-mean-square (LMS)
algorithm
96
The steps to the LMS algorithm, for each new sample at time
n, are:
Filter x [ n ] to produce:
y[n] =
M1
w
[ n ]x [ n ]
m [ n ]x [ n m ] = w
(9.17)
m=0
Form the estimation error:

e[n] = d[n] y[n]
(9.18)
Update the weight vector using step-size parameter :

w
[n + 1] = w
[ n ] + x [ n ]e [ n ]
(9.19)
For algorithm stability, the step-size must be chosen such

that
2
0 < < -------------------------------------tap-input power
(9.20)
where
M1
tap-input power =
E{ x[n k] }
(9.21)
k=0
In theory, (9.20) is equivalent to saying

2
0 < < ---------- max
(9.22)
where max is the maximum eigenvalue of R xx
97
Adaptive Filter Variations1

Prediction
s[n]
d[n]
Delay
x[n]
Adaptive
Filter
y[n]
e[n]
System Identification
s[n]
d[n]
Plant
Adaptive
Filter
x[n]
y[n]
e[n]
Equalization
s[n]
Training
Pattern
Delay
Plant/
Channel
Adaptive
+
Filter
x
[
n
]
+
Noise
y[n]
d[n]
e[n]
1. B. Widrow and S. Stearns, Adaptive Signal Processing, Prentice Hall, New

Jersey, 1985.
98
Adaptive Line Enhancement
Interference Canceling
d[n]
Signal +
Interference
Interference
x[n]
Adaptive
Filter
y[n]
e[n]

A relative of the interference canceling scheme shown above,
is the adaptive line enhancer (ALE)
Here we assume we have a narrowband signal (say a sinusoid) buried in broadband additive noise
x [ n ] = NB [ n ] + BB [ n ]
x[n ]
Adaptive
Filter
e[n]
y[n]
The filter adapts in such a way that a narrow passband forms

around the sinusoid frequency, thereby suppressing much of
the noise and improving the signal-to-noise ratio (SNR) in
y[ n]
99
Example: MATLAB ALE Simulation

A simple MATLAB simulation is constructed using a single
sinusoid at normalized frequency f o = 1 20 plus additive
white Gaussian noise
x [ n ] = A cos [ 2f o n ] + w [ n ]
(9.23)
The SNR is defined as

2
A
SNR = ---------22 w
(9.24)
function [n,x,y,e,wo,F,Wo] = lms_ale(SNR,N,M,mu)

%
[n,x,y,e,wo,F,Wo] = lms_ale(SNR,N,M,mu)
%
%*******LMS ALE Simulation************
%
SNR = Sinusoid SNR in dB
%
N = Number of simulation samples
%
M = FIR Filter length
%
mu = LMS step-size
%
%
n = Index vector
%
x = Noisy input
%
y = Filtered output
%
e = Error sequence
%
wo = Final value of weight vector
%
F = Frequency response axis vector
%
Wo = Frequency response of filter
%**************************************
%
Mark Wickert, 4/27/98
%
%
n
x
x
%
Sinusoid SNR = (A^2/2)/noise_var

= 0:N; % actually length N+1
= 1*cos(2*pi*1/20*n); % Here A = 1, Fo/Fs = 1/20
= x + sqrt(1/2/(10^(SNR/10)))*randn(1,N+1);
White Noise -> Delta = 1, so delay x by one sample
910

xd = filter([0 1],1,x);
% Initialize output vector y to zero
y = zeros(1,N+1);
% Initialize error vector e to zero
e = zeros(1,N+1);
% Initialize weight vector to zero
wo = zeros(1,M);
% Initialize filter memory to zero
z = zeros(1,M-1);
% Initialize a vector for holding xd of length M
xdm = zeros(1,M);
for k=1:N+1
% Filter one sample at a time
[y(k),z] = filter(wo,1,x(k),z);
% Form the error sequence
e(k) = x(k) - y(k);
% Update the weight vector
wo = wo + 2*mu*e(k)*xdm;
% Update vector used for correlation with e(k)
xdm = [xd(k) xdm(1:M-1)];
end %end loop on time index k
% Create filter frequency response
[Wo,F] = freqz(wo,1,512,1);
Wo = 20*log10(abs(Wo));
A simulation is run using 1000 samples, SNR = 10 dB, M =

64, and = 0.01 64
[n,x,y,e,wo,F,Wo] = lms_ale(10,1000,64,.01/64);
plot(n,e.^2)
subplot(211)
plot(n,x)
axis([600 1000 -1.5 1.5])
subplot(212)
plot(n,y)
plot(F,Wo)
911
Convergence Occurs
in Here (~275 samples)
ALE x[n] & y[n] for SNR = 10 dB, mu = 0.01/64

1.5
1
x[n]
0.5
0
-0.5
-1
-1.5
600
650
700
750
650
700
750
800
850
900
950
1000
800
850
900
950
1000
y[n]
0.5
0
-0.5
-1
600
Index n
912
C6x Code Examples
ALE Freq. Resp. for SNR = 10 dB, mu = 0.01/64

5
-5
Sinusoid Frequency
-10
|Wo(f)|dB
-15
-20
-25
-30
-35
-40
-45
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Normalized Frequency f
A C version of the above MATLAB code would be very similar except all of the vector operations would have to be
replaced by for loops
TBD
C6x Code Examples

A Two Channel Input Signal + Signal or Signal + Noise Canceller
In this first example a modified version of Chassaings
Adaptnoise.c program is considered
913
In this program the following system is implemented

d[n]
Signal1 +
Signal2 (or noise)
Signal2
x[n]
Adaptive
Filter
y[n]
e[n]
Notice that signal2 is present at both inputs in the exact same

form; in a real application the input x [ n ] would be similar to
signal2 as seen in d [ n ] , but not exactly the same
Floating point C code is shown below:
// Welch, Wright, & Morrow,
// Real-time Digital Signal Processing, 2011
// Modified by Mark Wickert February 2012 to include GPIO ISR start/stop postings
///////////////////////////////////////////////////////////////////////
// Filename: Adaptive_ISRs.c
//
// Synopsis: Interrupt service routine for codec data transmit/receive
//
///////////////////////////////////////////////////////////////////////
#include "DSP_Config.h"
// Function Prototypes
long int rand_int(void);
//
//
//
//
Data is received as 2 16-bit words (left/right) packed into one

32-bit word. The union allows the data to be accessed as a single
entity when transferring to and from the serial port, but still be
able to manipulate the left and right channels independently.
#define LEFT 0
#define RIGHT 1
914
C6x Code Examples

volatile union {
Uint32 UINT;
Int16 Channel[2];
} CodecDataIn, CodecDataOut;
/* add any global variables here */

//#define beta 1E-10;//rate of convergence
#define N 64
//# of weights (coefficients)
#define LEFT 0
//left channel
#define RIGHT 1
//right channel
float w[N];
//weights for adapt filter
float delay[N];
//input buffer to adapt filter
short output;
//overall output
short out_type = 1;//output type for slider
float alpha = 10;//alpha from slider
interrupt void Codec_ISR()
///////////////////////////////////////////////////////////////////////
// Purpose:
Codec interface interrupt service routine
//
// Input:
None
//
// Returns:
Nothing
//
// Calls:
CheckForOverrun, ReadCodecData, WriteCodecData
//
// Notes:
None
///////////////////////////////////////////////////////////////////////
{
/* add any local variables here */
WriteDigitalOutputs(1); // Write to GPIO J15, pin 6; begin ISR timing pulse
short i;
float yn=0, E=0, dplusn=0, desired=0, noise=0;
float beta = 1E-11;
beta *= alpha;//scale beta by slider alpha value
if(CheckForOverrun())// overrun error occurred (i.e. halted DSP)

return;
// so serial port is reset to recover
CodecDataIn.UINT = ReadCodecData();// get input data samples
/* add your code starting here */
desired = (float) CodecDataIn.Channel[ LEFT]; //input left chan
noise = (float) CodecDataIn.Channel[ RIGHT]; //input right chan
915

dplusn = desired + noise;
//desired+noise
delay[0] = noise;
//noise as input to adapt FIR
for (i = 0; i < N; i++)
{
yn += (w[i] * delay[i]);
}
//to calculate output of adaptive FIR
E = (desired + noise) - yn;
//"error" signal=(d+n)-yn
//output of adaptive filter
for (i = N-1; i >= 0; i--)

//to update weights and delays
{
w[i] = w[i] + beta*E*delay[i]; //update weights
delay[i] = delay[i-1];
//update delay samples
}
if (out_type == 1)
//if slider in position 1
{
output = (short) E; //error signal as overall output
}
else if (out_type == 2)
{
output = (short) dplusn;//desired+noise
}
CodecDataOut.Channel[ LEFT] = output;//overall output result
CodecDataOut.Channel[RIGHT] = 0;
/* end your code here */
WriteCodecData(CodecDataOut.UINT);// send output data to port
WriteDigitalOutputs(0); // Write to GPIO J15, pin 6; end ISR timing pulse
}
//White noise generator for filter noise testing
long int rand_int(void)
{
static long int a = 100001;
a = (a*125) % 2796203;
return a;
}
///////////////////////////////////////////////////////////////////////
// Filename: main.c
//
916
C6x Code Examples

// Synopsis: Main program file for demonstration code
//
///////////////////////////////////////////////////////////////////////
int main()
{
// call StartUp for application specific code
// defined in each application directory
StartUp();
// main stalls here, interrupts drive operation
while(1) {
;
}
}
///////////////////////////////////////////////////////////////////////
// Filename: StartUp.c
//
// Synopsis: Adaptive filter initialize
//
///////////////////////////////////////////////////////////////////////
#define N 64
//# of weights (coefficients)
extern float w[N];//weights for adapt filter
extern float delay[N];//input buffer to adapt filter
void StartUp()
{
InitDigitalOutputs();
// initialize filter weights and delay states
DSP_Init();
int T = 0;
for (T = 0; T < N; T++)
{
w[T] = 0.0;
delay[T] = 0.0;
}
//init variables
//init weights of adaptive FIR
//init buffer for delay sampless
917
The left channel input (desired) is added to the right channel input (noise) and forms the LMS input d [ n ] (dplusn)
The right channel input by itself forms the input x [ n ]
(noise)
A direct form FIR filter is implemented using floats in the
first for loop
In the second for loop the filter tap weights w[n] are
updated using coefficient beta which can be adjusted using
a GEL file; the filter history in delay[n] is also updated in
this loop
GEL Slider to toggle
between dplusn
and the error E
GEL slider to set the

value of coefficient
beta between 1E-12
and 1E-10
A second GEL slider controls which filter output goes to the

left codec channel
This adaptive filter was tested with two sinusoid inputs, one
at 1 kHz as the desired signal and one at 2.5 kHz the noise
signal
Analog outputs were captured using the PC sound card and
918
C6x Code Examples
spectrum analyzed using MATLAB

Desired 1.0 kHz Sinusoid + 2.5 kHz Noise
Power Spectrum Magnitude (dB)
30
20
10
0
10
20
30
40
50
60
500
1000
1500
2000
2500
Frequency
3000
3500
4000
Error Signal Output: Noise Removed

30
20
10
0
10
20
30
40
50
60
500
1000
1500
2000
2500
Frequency
3000
3500
4000
919
Adaptive Line Enhancer

In this first example a modified version of Chassaings
Adaptpredict.c program is considered
In this program the following system is implemented
x [ n ] = NB [ n ] + BB [ n ]
x[n ]
Adaptive
Filter
e[n]
y[n]
Floating point C code is shown below:

//Adaptpredict_AIC23.c Adaptive predictor to cancel interference
//using interrupts
// Enhanced version of Chassaing program
// by Mark Wickert Spring 2009
//**********************************************************
#include "DSK6713_AIC23.h"
//set sampling rate {8,16,24,32,44.1,48,96}
Uint32 fs=DSK6713_AIC23_FREQ_8KHZ;
#define DSK6713_AIC23_INPUT_MIC 0x0015
#define DSK6713_AIC23_INPUT_LINE 0x0011
Uint16 inputsource=DSK6713_AIC23_INPUT_LINE; // select input
//Uint16 inputsource=DSK6713_AIC23_INPUT_MIC; // select input
//***********************************************************
// For left and right channel processing
// left and right help in packed 32 bit word
#define LEFT 0
#define RIGHT 1
union {Uint32 combo; short channel[2];} AIC23_data;
920
C6x Code Examples

// Program globals
#define N 64
//# of coefficients of adapt FIR
#define Ns 10 //length splusn buffer
float splusn[Ns];
//buffer wideband signal+interference
float w[N];
//buffer for weights of adapt FIR
float delay[N];
//buffer for input to adapt FIR
Uint32 out_type = 1;//output type for slider
float alpha = 10;//alpha from slider
interrupt void c_int11()
//interrupt service routine
{
short i;
float yn, E;//yn=out adapt FIR, error signal
float wb_signal;//wideband desired signal
float noise;
//external interference
short output;//output to codec
float beta = 1E-13;//rate of convergence
beta *= alpha;//scale beta by slider alpha value
AIC23_data.combo = input_sample(); //in l&r as 32-bit
wb_signal = (float) AIC23_data.channel[LEFT]; //desired on l
//noise = (float) AIC23_data.channel[RIGHT];
//noise on r
noise = 0.1*((short) rand_int());//intern noise
splusn[0] = wb_signal + noise; //wb signal+interference
delay[0] = splusn[3];
//delayed input to adapt FIR
yn = 0;
//init output of adapt FIR
for (i = 0; i < N; i++)
yn += (w[i] * delay[i]);
E = splusn[0] - yn;
//output of adapt FIR filter

//(wb+noise)-out adapt FIR
for (i = N-1; i >= 0; i--)

{
w[i] = w[i]+(beta*E*delay[i]);//update wgts of FIR
delay[i+1] = delay[i];
//update buf delay samps
if (i < Ns-1) {
splusn[i+1] = splusn[i]; //update buf corr wb
}
921

}
//out_type = 2;
/*
if (out_type == 1)//if slider in position 1
{
output = ((short)E);//WB signal
//output = ((short)wb_signal);//WB signal
}
else if (out_type == 2)
{
output =((short)yn);//NB signal
//output =((short)wb_signal);//NB signal
}*/
output =((short)yn);
output_sample(output);//output resultto left chan.
return;
}
void main()
//main function
{
int T = 0;
for (T = 0; T < N; T++)
//init variables
{
w[T] = 0.0;
//init weights of adaptive FIR
delay[T] = 0.0;
//init buffer for delay samples
if (T < Ns)
{
splusn[T] = 0;
//init wideband+interference
}
}
comm_intr();
//init DSK, codec, McBSP
while(1);
//infinite loop
}
//16 bit word random number generator
int rand_int(void)
{
static int a = 100001;
a = (a*125) % 2796203;
922
C6x Code Examples

return a;
}
The inputs are assumed to be a wideband signal, e.g., noise,

and a narrowband signal, e.g., a sinusoid
The inputs arrive at the left and right codec inputs and are
summed in code as wb_signal + noise
An input delay buffer of Ns samples is created in the array
splusn[], here Ns = 10
splusn[3] serves as x [ n ] , so here = 3
A GEL file allows the left channel of the codec to output
either the wideband signal at setting 1 or the narrowband signal at setting two
The adaptation coefficient is set to be 10
12
923
A 1.8 kHz sinusoid is input from a function generator and

wideband noise is input from another function generator, the
frequency response of the converged filter weights w[n] is
displayed in CCS
The corresponding spectra of y [ n ] and e [ n ] are captured via

the PC sound card
The Wideband Signal Output
30
20
10
0
10
20
30
40
50
60
924
500
1000
1500
2000
2500
Frequency
3000
3500
4000
C6x Code Examples
The Narrowband Signal Output
30
20
10
0
10
20
30
40
50
60
500
1000
1500
2000
2500
Frequency
3000
3500
4000
925
926

Adaptive Filter Notes

Uploaded by

Copyright:

Available Formats

Adaptive Filter Notes

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Adaptive Filter Notes

Uploaded by

Copyright:

Available Formats

Chapter

The term adaptive filter implies changing the characteristic of a

Chapter 9 Adaptive Filters

where the w m are referred to as the filter weights

E{e [n]} = E{(d[n] y[n]) }

is minimized, that is we obtain the minimum mean-squared

From the orthogonality principle1 we choose the weights

1. A. Papoulis, Probability, Random Variables, and Stochastic Processes,

ECE 5655/4655 Real-Time DSP

for k = 0, 1, , M 1 are known as the Wiener-Hopf or

where R xx is the M M correlation matrix associated with

w o is the optimum weight vector given by

Chapter 9 Adaptive Filters

and p xd is the cross-correlation vector given by

The optimal weight vector is given by

where N is the sample block size

Adaptive Wiener Filter

ECE 5655/4655 Real-Time DSP

Adaptive Wiener Filter

With respect to the weight vector w , the error e [ n ] can be

Mean Square Error

Error Surface for M = 2

ECE 5655/4655 Real-Time DSP

Chapter 9 Adaptive Filters

A block diagram of this adaptive Wiener (FIR) filter is shown

Adaptive Wiener FIlter

A practical implementation involves estimating the gradient

ECE 5655/4655 Real-Time DSP

Adaptive Wiener Filter

Form the estimation error:

Update the weight vector using step-size parameter :

For algorithm stability, the step-size must be chosen such

In theory, (9.20) is equivalent to saying

where max is the maximum eigenvalue of R xx

ECE 5655/4655 Real-Time DSP

Chapter 9 Adaptive Filters

Adaptive Filter Variations1

1. B. Widrow and S. Stearns, Adaptive Signal Processing, Prentice Hall, New

ECE 5655/4655 Real-Time DSP

Adaptive Line Enhancement

Adaptive Line Enhancement

The filter adapts in such a way that a narrow passband forms

Chapter 9 Adaptive Filters

Example: MATLAB ALE Simulation

The SNR is defined as

function [n,x,y,e,wo,F,Wo] = lms_ale(SNR,N,M,mu)

Sinusoid SNR = (A^2/2)/noise_var

ECE 5655/4655 Real-Time DSP

Adaptive Line Enhancement

A simulation is run using 1000 samples, SNR = 10 dB, M =

ECE 5655/4655 Real-Time DSP

Chapter 9 Adaptive Filters

ALE x[n] & y[n] for SNR = 10 dB, mu = 0.01/64

ECE 5655/4655 Real-Time DSP

C6x Code Examples

ALE Freq. Resp. for SNR = 10 dB, mu = 0.01/64

C6x Code Examples

ECE 5655/4655 Real-Time DSP

Chapter 9 Adaptive Filters