# Least-squares prediction and smoothing

Other languages:
English • ‎español
Series Geophysical References Series Digital Imaging and Deconvolution: The ABCs of Seismic Exploration and Processing Enders A. Robinson and Sven Treitel 10 http://dx.doi.org/10.1190/1.9781560801610 9781560801481 SEG Online Store

An important operation is signal extrapolation or, in other words, prediction. Smoothing is yet another important signal operation. Often, the desired quantity is observable only after it has been corrupted or altered by mixture — perhaps additively — with other signals. Such a corrupting signal commonly is called noise, even for cases in which the signal has a physically legitimate presence. For example, reverberations often are regarded as noise, but reverberations are a physical necessity. Such signal-generated noise contains valuable information in its own right, and the trick is to use this information to an advantage whenever possible.

Often it is important to ascertain, in the least-squares sense, what the data would have been like without contamination by noise. Smoothing might be the complete problem to be addressed; alternatively, it might be combined with a prediction problem, which means that we wish to know what the uncontaminated signal will do in the future. Whereas the smoothing problem is clearly distinguishable from the prediction problem, mixed problems involving elements of both are greatly important. Indeed, good smoothing performance usually depends on introduction of a sufficient delay. If the delay is negative, filter performance suffers. On the other hand, the filter becomes a smoothing predictor, which often is a useful tool. Seismic processing must be innovative. Most geophysical data-processing methods are mathematical and physical hybrids and are based on a particular geophysical model.

Before an actual digital filter can be designed, it is necessary to build a model. The time-honored least-squares model requires three signals: the input signal ${\displaystyle x=(\ x_{0},x_{1},x_{2},...)}$, the actual output signal ${\displaystyle y=\left(y_{0}{,\ }y_{\rm {l}}{,\ }y_{2}{,\ .\ ..}\right)}$, and the desired output signal ${\displaystyle z=(\ z_{0},z_{1},z_{2},\ldots )}$. The filter is a finite-impulse-response (FIR) filter. Let such a filter k have a finite number of coefficients, given by

 {\displaystyle {\begin{aligned}k=\left(k_{0}{\ ,\ }k_{1}{\ ,\ }k_{2}{,\dots ,}k_{N-1}\right).\end{aligned}}} (1)

This filter is said to be of length N. The output is the convolution of the filter with the input; that is, ${\displaystyle y=k*x}$. More specifically, the output value ${\displaystyle y_{n}}$ at time n is given by the discrete convolution formula

 {\displaystyle {\begin{aligned}y_{n}=k_{0}x_{n}+k_{\rm {1}}x_{n-1}+\dots +k_{N-{\rm {1}}}x_{n-N+1}.\end{aligned}}} (2)

The input signal and the desired output signal are known. The problem is to find the filter by the least-squares criterion. At time instant n, the value of the output is ${\displaystyle y_{n}}$ and the value of the desired output is ${\displaystyle z_{n}}$. The difference between the desired output and the actual output is the error ${\displaystyle e_{n}=z_{n}-y_{n}}$. The quantity

 {\displaystyle {\begin{aligned}I=\sum _{n}{{\left(z_{n}-y_{n}\right)}^{2}}\end{aligned}}} (3)

is the error energy. The filter f that minimizes the value of the error energy is the optimum filter in the least-squares sense. Setting its partial derivatives with respect to each of the filter coefficients equal to zero minimizes the error energy. The result of that minimization is the set of a system of N linear simultaneous equations known as the normal equations. In their matrix form, the normal equations are

 ${\displaystyle \left[{\begin{array}{l}r_{0}\;\;\;\;\;\;\;r_{1}\;\;\;\;\;\;\ldots \;\;\;r_{N}\\r_{1}\;\;\;\;\;\;\;r_{0}\;\;\;\;\;\ldots \;\;\;r_{N-1}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ldots \\r_{N-1\;}\;\;r_{N-1}\;\;\ldots \;\;\;r_{0}\\\end{array}}\right]\left[{\begin{array}{l}k_{0}\\k_{1}\\\ldots \\k_{N-1}\\\end{array}}\right]=\left[{\begin{array}{l}g_{0}\\g_{1}\\\ldots \\g_{N-1}\\\end{array}}\right],}$ (4)

which were derived in detail in Chapter 9.

The known quantities in this system are the autocorrelation ${\displaystyle r_{n}}$ of the input signal and the positive-lag values of crosscorrelation ${\displaystyle g_{n}}$ of the desired output signal with the input signal. The unknown quantities are the values of the filter coefficients. The autocorrelation of the input is

 {\displaystyle {\begin{aligned}r_{n}=\sum _{i}{x_{i+n}}x_{i},\end{aligned}}} (5)

whereas the crosscorrelation between the desired output and the input is

 {\displaystyle {\begin{aligned}g_{n}=\sum _{i}{z_{i+n}}x_{i}.\end{aligned}}} (6)

The solution of the normal equations (equations 4) yields the filter coefficients. The square matrix in the normal equations has the autocorrelation coefficients arranged in Toeplitz form. That is, all the autocorrelation coefficients along any diagonal (the main diagonal or any subdiagonal) are the same. Because of this Toeplitz structure, the well-known Levinson algorithm can be used in the solution of the normal equations (Robinson and Treitel, 2000[1], p. 163-169).

Let us now derive the least-squares prediction filter, and let ${\displaystyle k=}$${\displaystyle \left(k_{0}{,\ }k_{\rm {l}}{\ ,\ }k_{N-{\rm {l}}}\right)}$ denote this prediction filter. Such a filter uses the input signal’s past values to predict that signal’s future values. The input is the signal ${\displaystyle x_{n}}$, and the desired output is the time-advanced version of the input. Let the prediction distance be given by the positive integer ${\displaystyle \alpha }$. The input to the filter is the input ${\displaystyle x_{n}}$ at present time n. The desired output ${\displaystyle z_{n}}$ is the input ${\displaystyle x_{n+\alpha }}$ at the future time ${\displaystyle n+\alpha }$. The prediction filter is designed so that the output ${\displaystyle y_{n}}$ at the present time n is an optimum estimate of the future value ${\displaystyle x_{n+\alpha }}$. If this estimate is denoted by ${\displaystyle {\hat {x}}_{n+\alpha }}$, the filter’s action can be represented by the convolution

 {\displaystyle {\begin{aligned}{\hat {x}}_{n+\alpha }=k_{0}x_{n}+k_{1}x_{n-1}+\ldots +k_{N-1}x_{n-N+1}.\end{aligned}}} (7)

Equation 7 says that the prediction operator acts on an input up to time n and estimates its value at some future time ${\displaystyle n+\alpha }$. Normal equations 4 now can be used.

The crosscorrelation between the desired output and the input is given by the positive lag coefficients of the crosscorrelation between the time-advanced trace and the trace

 {\displaystyle {\begin{aligned}g_{n}=\sum _{i}{x_{i+n+\alpha }}x_{i}=r_{n+\alpha }.\end{aligned}}} (8)

As equation 8 shows, the crosscorrelation between the desired output and the input is equal to the autocorrelation of the input for lags greater than or equal to ${\displaystyle \alpha }$. Thus, the normal equations for the prediction filter can be written as

 ${\displaystyle \left[{\begin{array}{l}r_{0}&r_{1}&.&.&r_{N-1}&\\r_{1}&r_{0}&&&r_{N-2}&\\&.&.&.&\\\\r_{N}&r_{N-1}&.&.&r_{0}&\end{array}}\right]\left[{\begin{array}{l}k_{0}\\k_{1}\\.\\.\\k_{N}\\\end{array}}\right]=\left[{\begin{array}{l}r_{\alpha }\\r_{\varepsilon +1}\\.\\.\\r_{\alpha +N}\\\end{array}}\right].}$ (9)

Their solution yields the coefficients of the optimum prediction filter we seek.

## References

1. Robinson, E. A., and S. Treitel, 2000, Geophysical signal analysis: SEG.