Adaptive Filter Notes

Chapter

Adaptive Filters

9

Introduction

The term adaptive filter implies changing the characteristic of a filter in some automated fashion to obtain the best possible 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 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 filter.1 Wiener filter. + d[ n ] e [ n ] Error Σ Signal

Desired Signal

x [ n ] Observation

Wiener Filter

y [ n ]

MMSE Estimate of d [n]

Filter Theory, fourth edition, Prentice Hall, 2002. 1.Simon Haykin, Adaptive Filter ECE 5655/4655 Real-Time DSP

9–1

Chapter 9 • Adaptive Filters

Wiener Filter  An M tap tap discrete-time Wiener filter filter is of the form M – 1

y [ n ] =

∑

w m x [ n – m ]

(9.1)

m=0

where the w m are referred to as the filter weights – Note: Note: (9.1) (9.1) tells us that that the Wiener Wiener filter filter is just just an M -tap -tap FIR filter • The quality quality of the filtered filtered or estimated estimated signal signal [ n ] is deterdetermined from the error sequence e [ n ] = d [ n ] – y [ n ] • The The weig weight htss w m , m = 0, 1, …, M – 1 ar aree chos chosen en such such that that 2

2

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

(9.2)

is minimized, that is we obtain the minimum mean-squared error (MMSE) • The The optim optimal al weig weights hts are found found by sett setting ing

∂

2

E { e [ n ] } = 0, m = 0, 1, …, M – 1

∂ wm

(9.3)

• From From the the ort ortho hogo gona nali lity ty pri princ ncip iple le1 we choose the weights such that the error e [ n ] is orthogonal to the observations observations (data), i.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. 9–2

ECE 5655/4655 Real-Time DSP

Wiener Filter

– This results in a filter that is optimum in the sense of minimum mean-square error • The resulting system of equations M – 1

∑

w m E { x [ n – k ] x [ n – m ] } = E { x [ n – k ] ( d [ n ] ) }

(9.5)

m=0

or M – 1

∑

w m φ xx [ m – k ] = φ xd [ – k ]

(9.6)

m=0

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 xx w o

where R xx is the

= p x

(9.7)

× 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


9–3


wo =

w o0 w o1 … w oM – 1

T

(9.9)

and p xd is the cross-correlation vector given by p xd =

φ xd [ 0 ] φ xd [ – 1 ] … φ xd [ 1 – M ]

T

(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 N – 1

1 φ xx [ k ] ≅ ---- ∑ x [ n + k ] x [ n ] N

(9.12)

n=0

N – 1

1 φ xd [ k ] ≅ ---- ∑ x [ n + k ] d [ n ] N

(9.13)

n=0

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

9–4


Adaptive Wiener Filter

• 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)

Optimum at [1,1]

r o r r E e r a u q S n a e M

w1

w0

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


9–5


• 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., T

ˆ xx = x [ n ] x [ n ]

(9.14)

p ˆ xd = x [ n ] d [ n ]

(9.15)

where x[n] =

T

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 9–6


Adaptive Wiener Filter

• The steps to the LMS algorithm, for each new sample at time n, are: – Filter x [ n ] to produce: M – 1

y [ n ] =

∑

T

w ˆ m [ n ] x [ n – m ] = w ˆ [ n ]x[n ]

(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 M – 1

tap-input power =

∑

2

E { x [ n – k ] }

(9.21)

k = 0

• In theory, (9.20) is equivalent to saying 2 0 < µ < -----------

λ ma x

(9.22)

where λ max is the maximum eigenvalue of R xx


9–7


Adaptive Filter Variations1

• Prediction s [ n ]

d[ n]

Adaptive Filter x [ n ] y [ n ]

Delay

+

Σ e[ n]

• System Identification s [ n ]

d[ n]

Plant Adaptive Filter y [ n ]

x [ n ]

+

Σ e[n]

• Equalization s [ n ]

Training Pattern

Delay

Plant/ Channel

Adaptive Σ + Filter x [ n ] +

y [ n ]

d[n]

+

Σ e[n]

Noise

1. B. Widrow and S. Stearns, Adaptive Signal Processing, Prentice Hall, New Jersey, 1985. 9–8


Adaptive Line Enhancement

• Interference Canceling d[n]

Signal + Interference

Interference’

x [ n ]

Adaptive Filter y [ n ]

-

+

Σ e[ n ]

Adaptive Line Enhancement • 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 ]

+

Σ -

–∆

z

Adaptive Filter x [ n – ∆ ]

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 ] ECE 5655/4655 Real-Time DSP

9–9


Example: MATLAB ALE Simulation

• A simple MATLAB simulation is constructed using a single sinusoid at normalized frequency o = 1 ⁄ 20 plus additive white Gaussian noise x [ n ] = A cos [ 2 π f o n ] + w [ n ]

(9.23)

• The SNR is defined as 2

A SNR = ---------2 2 σ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

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

9–10


Adaptive Line Enhancement 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)


9–11


Convergence Occurs in Here (~275 samples)

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

] n [ x

0.5 0

-0.5 -1 -1.5 600

650

700

750

800

850

900

950

1000

650

700

750

800

850

900

950

1000

1

0.5

] n [ y

0

-0.5

-1 600

Index n

9–12


C6x Code Examples

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

0

-5

Sinusoid Frequency

-10

-15

B d | ) f (-20 o W | -25

-30

-35

-40

-45 0

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 Chassaing’s Adaptnoise.c program is considered


9–13


• 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

9–14


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) //left channel #define LEFT 0 //right channel #define RIGHT 1 float w[N]; //weights for adapt filter //input buffer to adapt filter float delay[N]; 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


9–15

Chapter 9 • Adaptive Filters 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

//to update weights and delays for (i = N-1; i >= 0; i--) { 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; } // Welch, Wright, & Morrow, // Real-time Digital Signal Processing, 2011 /////////////////////////////////////////////////////////////////////// // Filename: main.c //

9–16


C6x Code Examples // Synopsis: Main program file for demonstration code // /////////////////////////////////////////////////////////////////////// #include "DSP_Config.h" int main() { // call StartUp for application specific code // defined in each application directory StartUp();

// main stalls here, interrupts drive operation while(1) { ; } } // Welch, Wright, & Morrow, // Real-time Digital Signal Processing, 2011 /////////////////////////////////////////////////////////////////////// // Filename: StartUp.c // // Synopsis: Adaptive filter initialize // /////////////////////////////////////////////////////////////////////// #include "DSP_Config.h" #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

}


9–17


• 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

9–18


C6x Code Examples

spectrum analyzed using MATLAB

Desired 1.0 kHz Sinusoid + 2.5 kHz Noise 30 ) B 20 d ( e 10 d u t i n 0 g a M−10 m u r −20 t c e p−30 S r e−40 w o P−50

−60

0

500

1000

1500

2000 Fre uenc

2500

3000

3500

4000

Error Signal Output: Noise Removed 30 ) B 20 d ( e 10 d u t i n 0 g a M−10 m u r −20 t c e p−30 S r e−40 w o P−50

−60

0

500

1000


1500

2000 Fre uenc

2500

3000

3500

4000

9–19


Adaptive Line Enhancer

• In this first example a modified version of Chassaing’s Adaptpredict.c program is considered • In this program the following system is implemented x [ n ] = NB [ n ] + BB [ n ]

+

Σ -

–∆

z

Adaptive Filter x [ n – ∆ ]

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;

9–20


C6x Code Examples // Program globals #define N 64 #define Ns 10

//# of coefficients of adapt FIR //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 }


9–21

Chapter 9 • Adaptive Filters } //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(); while(1);

//init DSK, codec, McBSP //infinite loop

} //16 bit word random number generator int rand_int(void) { static int a = 100001; a = (a*125) % 2796203;

9–22


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

9–23


• 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 ) B 20 d ( e 10 d u t i n 0 g a M−10 m u r −20 t c e p−30 S r e−40 w o P−50

−60

9–24

0

500

1000

1500

2000 Fre uenc

2500

3000

3500

4000


C6x Code Examples

The Narrowband Signal Output 30 ) B 20 d ( e 10 d u t i n 0 g a M−10 m u r −20 t c e p−30 S r e−40 w o P−50

−60

0

500

1000


1500

2000 Fre uenc

2500

3000

3500

4000

9–25


9–26


Adaptive Filter Notes

Recommend Documents