Wiener Filters (Vaseghi - Advanced Digital Signal Processing and Noise Reduction)

Advanced Digital Signal Processing and Noise Reduction, Second Edition.Saeed V. VaseghiCopyright © 2000 John Wiley & Sons LtdISBNs: 0-471-62692-9 (Hardback): 0-470-84162-1 (Electronic)6WIENER FILTERS6.16.26.36.46.56.66.76.8Wiener Filters: Least Square Error EstimationBlock-Data Formulation of the Wiener FilterInterpretation of Wiener Filters as Projection in Vector SpaceAnalysis of the Least Mean Square Error SignalFormulation of Wiener Filters in the Frequency DomainSome Applications of Wiener FiltersThe Choice of Wiener Filter OrderSummaryWiener theory, formulated by Norbert Wiener, forms thefoundation of data-dependent linear least square error filters.Wiener filters play a central role in a wide range of applicationssuch as linear prediction, echo cancellation, signal restoration, channelequalisation and system identification.

The coefficients of a Wiener filterare calculated to minimise the average squared distance between the filteroutput and a desired signal. In its basic form, the Wiener theory assumesthat the signals are stationary processes. However, if the filter coefficientsare periodically recalculated for every block of N signal samples then thefilter adapts itself to the average characteristics of the signals within theblocks and becomes block-adaptive. A block-adaptive (or segmentadaptive) filter can be used for signals such as speech and image that maybe considered almost stationary over a relatively small block of samples.

Inthis chapter, we study Wiener filter theory, and consider alternativemethods of formulation of the Wiener filter problem. We consider theapplication of Wiener filters in channel equalisation, time-delay estimationand additive noise reduction. A case study of the frequency response of aWiener filter, for additive noise reduction, provides useful insight into theoperation of the filter. We also deal with some implementation issues ofWiener filters.179Least Square Error Estimation6.1 Wiener Filters: Least Square Error EstimationWiener formulated the continuous-time, least mean square error, estimationproblem in his classic work on interpolation, extrapolation and smoothingof time series (Wiener 1949). The extension of the Wiener theory fromcontinuous time to discrete time is simple, and of more practical use forimplementation on digital signal processors.

A Wiener filter can be aninfinite-duration impulse response (IIR) filter or a finite-duration impulseresponse (FIR) filter. In general, the formulation of an IIR Wiener filterresults in a set of non-linear equations, whereas the formulation of an FIRWiener filter results in a set of linear equations and has a closed-formsolution. In this chapter, we consider FIR Wiener filters, since they arerelatively simple to compute, inherently stable and more practical.

The maindrawback of FIR filters compared with IIR filters is that they may need alarge number of coefficients to approximate a desired response.Figure 6.1 illustrates a Wiener filter represented by the coefficient vector w.The filter takes as the input a signal y(m), and produces an output signalxˆ (m) , where xˆ (m) is the least mean square error estimate of a desired ortarget signal x(m). The filter input–output relation is given byP −1xˆ ( m ) = ∑ wk y ( m − k )k =0T(6.1)=w ywhere m is the discrete-time index, yT=[y(m), y(m–1), ..., y(m–P–1)] is thefilter input signal, and the parameter vector wT=[w0, w1, ..., wP–1] is theWiener filter coefficient vector.

In Equation (6.1), the filtering operation isexpressed in two alternative and equivalent forms of a convolutional sumand an inner vector product. The Wiener filter error signal, e(m) is definedas the difference between the desired signal x(m) and the filter output signalxˆ (m) :e ( m ) = x ( m ) − xˆ ( m )(6.2)= x ( m) − w T yIn Equation (6.2), for a given input signal y(m) and a desired signal x(m),the filter error e(m) depends on the filter coefficient vector w.180Wiener FiltersInput y(m)y(m–1)–1ww=Ry(m–2)z–1zw0...z–1y(m–P–1)wwP–121–1ryy xyFIR Wiener FilterDesired signalx(m)^x(m)Figure 6.1 Illustration of a Wiener filter structure.To explore the relation between the filter coefficient vector w and theerror signal e(m) we expand Equation (6.2) for N samples of the signalsx(m) and y(m): e(0)   x (0)   y (0)   e (1)   x (1)   y (1) e( 2)  =  x ( 2)  −  y ( 2)     e ( N − 1)   x ( N − 1)   y ( N − 1)  y ( −1)y (0)y (1)y ( −2 )y ( −1)y (0)y ( N − 2)y ( N − 3)y (1 − P ) y(2 − P) y (3 − P ) y ( N − P )  w0  w1  w  2 w  P −1 (6.3)In a compact vector notation this matrix equation may be written ase = x − Yw(6.4)where e is the error vector, x is the desired signal vector, Y is the inputsignal matrix and Yw = xˆ is the Wiener filter output signal vector.

It isassumed that the P initial input signal samples [y(–1), . . ., y(–P–1)] areeither known or set to zero.In Equation (6.3), if the number of signal samples is equal to thenumber of filter coefficients N=P, then we have a square matrix equation,and there is a unique filter solution w, with a zero estimation error e=0, such181Least Square Error Estimationthat xˆ =Yw = x . If N < P then the number of signal samples N isinsufficient to obtain a unique solution for the filter coefficients, in this casethere are an infinite number of solutions with zero estimation error, and thematrix equation is said to be underdetermined. In practice, the number ofsignal samples is much larger than the filter length N>P; in this case, thematrix equation is said to be overdetermined and has a unique solution,usually with a non-zero error. When N>P, the filter coefficients arecalculated to minimise an average error cost function, such as the averageabsolute value of error E [|e(m)|], or the mean square error E [e2(m)], whereE [.] is the expectation operator.

The choice of the error function affects theoptimality and the computational complexity of the solution.In Wiener theory, the objective criterion is the least mean square error(LSE) between the filter output and the desired signal. The least squareerror criterion is optimal for Gaussian distributed signals.

As shown in thefollowings, for FIR filters the LSE criterion leads to a linear and closedform solution. The Wiener filter coefficients are obtained by minimising an2average squared error function E [e (m)] with respect to the filtercoefficient vector w. From Equation (6.2), the mean square estimation erroris given byE[ e 2 ( m )] = E[( x ( m ) − w T y ) 2 ]= E[ x 2 ( m )] − 2 w T E[ yx ( m )] + w T E[ yy T ]w(6.5)= rxx ( 0 ) − 2w T r yx + w T R yy wwhere Ryy=E [y(m)yT(m)] is the autocorrelation matrix of the input signaland rxy=E [x(m)y(m)] is the cross-correlation vector of the input and thedesired signals.

An expanded form of Equation (6.5) can be obtained asP −1P −1P −1k =0k =0j =0E[ e 2 ( m )] = rxx ( 0 ) − 2 ∑ wk ryx ( k ) + ∑ wk∑ w j ryy ( k − j )(6.6)where ryy(k) and ryx(k) are the elements of the autocorrelation matrix Ryyand the cross-correlation vector rxy respectively. From Equation (6.5), themean square error for an FIR filter is a quadratic function of the filtercoefficient vector w and has a single minimum point. For example, for afilter with only two coefficients (w0, w1), the mean square error function is a182Wiener FiltersE [e2]ww1optimalw0Figure 6.2 Mean square error surface for a two-tap FIR filter.bowl-shaped surface, with a single minimum point, as illustrated in Figure6.2.

The least mean square error point corresponds to the minimum errorpower. At this optimal operating point the mean square error surface haszero gradient. From Equation (6.5), the gradient of the mean square errorfunction with respect to the filter coefficient vector is given by∂E[ e 2 ( m )] = − 2E[ x ( m ) y ( m )]+ 2w T E[ y ( m ) y T ( m )]∂w(6.7)T= − 2 r yx + 2 w R yywhere the gradient vector is defined as∂∂∂∂  ∂=,,, ,∂ w  ∂ w0 ∂ w1 ∂ w2∂ wP −1 T(6.8)The minimum mean square error Wiener filter is obtained by settingEquation (6.7) to zero:Ryy w = r yx(6.9)183Least Square Error Estimationor, equivalently,w = R −yy1 r yx(6.10)In an expanded form, the Wiener filter solution Equation (6.10) can bewritten asryy (1)ryy ( 2) w0   ryy (0) ryy (0)ryy (1) w1   ryy (1) w  =  r ( 2)r yy (1)r yy (0) 2   yy   w   r ( P − 1) r ( P − 2) r ( P − 3)yyyy P −1   yy ryy ( P − 1)  −1  ryx (0)   ryy ( P − 2)   ryx (1)  r yy ( P − 3)   ryx ( 2)  ryy ( 0)   r ( P − 1)  yx(6.11)From Equation (6.11), the calculation of the Wiener filter coefficientsrequires the autocorrelation matrix of the input signal and the crosscorrelation vector of the input and the desired signals.In statistical signal processing theory, the correlation values of arandom process are obtained as the averages taken across the ensemble ofdifferent realisations of the process as described in Chapter 3.

However inmany practical situations there are only one or two finite-durationrealisations of the signals x(m) and y(m). In such cases, assuming the signalsare correlation-ergodic, we can use time averages instead of ensembleaverages. For a signal record of length N samples, the time-averagedcorrelation values are computed asryy (k ) =1NN −1∑ y ( m) y ( m + k )(6.12)m =0Note from Equation (6.11) that the autocorrelation matrix Ryy has a highlyregular Toeplitz structure.

A Toeplitz matrix has constant elements alongthe left–right diagonals of the matrix. Furthermore, the correlation matrix isalso symmetric about the main diagonal elements. There are a number ofefficient methods for solving the linear matrix Equation (6.11), includingthe Cholesky decomposition, the singular value decomposition and the QRdecomposition methods.184Wiener Filters6.2 Block-Data Formulation of the Wiener FilterIn this section we consider an alternative formulation of a Wiener filter for ablock of N samples of the input signal [y(0), y(1), ..., y(N–1)] and thedesired signal [x(0), x(1), ..., x(N–1)]. The set of N linear equationsdescribing the Wiener filter input/output relation can be written in matrixform as xˆ (0)   y (0)  xˆ (1)   y (1) xˆ (2)   y ( 2)=  xˆ ( N − 2)   y ( N − 2)  xˆ ( N − 1)   y ( N − 1)y ( −1)y ( −2)y (0)y (−1)y (1)y ( 0)y ( N − 3)y ( N − 4)y ( N − 2)y ( N − 3)y (2 − P )y (3 − P )y (4 − P )y ( N − P)y ( N + 1 − P)y (1 − P)   w0 y (2 − P )   w1 y (3 − P)   w2 y ( N − 1 − P ) wP − 2 y ( N − P )   wP −1 (6.13)Equation (6.13) can be rewritten in compact matrix notation asxˆ =Y w(6.14)The Wiener filter error is the difference between the desired signal and thefilter output defined ase = x − xˆ(6.15)= x − YwThe energy of the error vector, that is the sum of the squared elements ofthe error vector, is given by the inner vector product ase T e = ( x − Y w) T ( x −Y w)= x T x − x TYw − w TY T x + w TY TYw(6.16)The gradient of the squared error function with respect to the Wiener filtercoefficients is obtained by differentiating Equation (6.16):∂ e Te= − 2 x TY + 2 w TY TY∂w(6.17)185Block-Data FormulationThe Wiener filter coefficients are obtained by setting the gradient of thesquared error function of Equation (6.17) to zero, this yields(Y TY )w = Y T xor(w = Y TY)−1Y T x(6.18)(6.19)Note that the matrix YTY is a time-averaged estimate of the autocorrelationmatrix of the filter input signal Ryy, and that the vector YTx is a timeaveraged estimate of rxy the cross-correlation vector of the input and thedesired signals.

Theoretically, the Wiener filter is obtained fromminimisation of the squared error across the ensemble of differentrealisations of a process as described in the previous section. For acorrelation-ergodic process, as the signal length N approaches infinity theblock-data Wiener filter of Equation (6.19) approaches the Wiener filter ofEquation (6.10):( )lim  w = Y TYN →∞ −1 TY x  = R −yy1rxy(6.20)Since the least square error method described in this section requires ablock of N samples of the input and the desired signals, it is also referred toas the block least square (BLS) error estimation method.