Wiener (least-squares) inverse filters

Problems in Exploration Seismology and their Solutions

Series	Geophysical References Series
Title	Problems in Exploration Seismology and their Solutions
Author	Lloyd P. Geldart and Robert E. Sheriff
Chapter	9
Pages	295 - 366
DOI	http://dx.doi.org/10.1190/1.9781560801733
ISBN	ISBN 9781560801153
Store	SEG Online Store

Problem 9.22a

Plot cumulative energy as a function of time for wavelets $A=[1,\;-2,\;3]$ and $B=[3,\;-2,\;1]$

Background

A filter that will change a signal so as to make it as close as possible to a desired signal can be designed in a least-squares sense to minimize the sum of the squares of the “errors” (differences between a desired signal and the filtered signal). Such a filter is called a least-squares filter or Wiener filter.

Using $g_{t}$ to denote the original signal and $h_{t}$ the desired signal, the sum of the errors squared is

${\begin{aligned}E=\mathop {\sum } \limits _{t=0}^{n}(h_{t}-f_{t}*g_{t})^{2}.\end{aligned}}$

(9.22a)

Since all quantities on the right-hand side of equation (9.22a) are fixed except the filter elements $f_{i}$, we vary the coefficients $f_{i}$ to minimize $E$. We differentiate $E$ with respect to each of the elements $f_{i}$ and equate the derivatives to zero. This gives $(n+1)$ equations which can be solved for the $(n+1)$ elements of the filter. Differentiating $E$, we get

${\begin{aligned}{\frac {\partial E}{\partial f_{i}}}=0=\mathop {\sum } \limits _{t=0}^{n}(h_{t}-f_{t}*g_{t}){\frac {\partial (f_{t}*g_{t})}{\partial f_{i}}}\\=\mathop {\sum } \limits _{t=0}^{n}\left(h_{t}-\mathop {\sum } \limits _{t=0}^{n}g_{k}\ f_{t-k}\right){\frac {\partial }{\partial f_{i}}}\left(\mathop {\sum } \limits _{t=0}^{n}g_{k}\ f_{t-k}\right),\end{aligned}}$

using equation (9.2b) to replace $f_{t}*g_{t}$ with a summation. The derivative now becomes

${\begin{aligned}{\frac {\partial }{\partial f_{i}}}\left(\mathop {\sum } \limits _{t=0}^{n}g_{k}\ f_{t-k}\right)=g_{t-i},\end{aligned}}$

since the only nonzero term in the differentiation is that in which $f_{t}$ appears, that is, $t-k=i$, so $k=t-i$. Substituting this result, we get

${\begin{aligned}\mathop {\sum } \limits _{t=0}^{n}\left(h_{t}-\mathop {\sum } \limits _{k}^{}g_{k}\ f_{t-k}\right)g_{t-i}=0=\mathop {\sum } \limits _{t=0}^{n}h_{t}\ g_{t-i}-\mathop {\sum } \limits _{t=0}^{n}\left(\mathop {\sum } \limits _{k}^{}(g_{k}\ f_{t-k})g_{t-i}\right).\end{aligned}}$

The first term is $\phi _{gh}$($i$) from equations (9.8a) and (9.8b). Interchanging the order of summation in the right-hand term, letting $j=t-k$, and summing over $j$, the equation becomes

${\begin{aligned}\mathop {\sum } \limits _{j=0}^{n}\left(f_{j}\mathop {\sum } \limits _{t=0}^{n}g_{t-j}\ g_{t-i}\right)=\mathop {\sum } \limits _{j=0}^{n}f_{j}\phi _{gh}(i-j)\end{aligned}}$

(see problem 9.21a). Thus we arrive at the normal equations

${\begin{aligned}\mathop {\sum } \limits _{j=0}^{n}\phi _{gg}(i-j)f_{j}=\phi _{gh}(i),\qquad i=0,1,2,\ldots ,n.\end{aligned}}$

(9.22b)

Solution

The transforms of $A$ and $B$ are $(1-2z+3z^{2})$ and $(3-2z+z^{2})$, so the roots of $A$ are ${\frac {1}{3}}(1\pm j{\sqrt {2}})$ and those of $B$ are $(1\pm j{\sqrt {2}})$, the magnitudes of the roots being $1/{\sqrt {3}}$ and ${\sqrt {3}}$. Thus $A$ is maximum-phase and $B$ is minimum-phase.

The energy of a wavelet at any instant is proportional to its amplitude squared. We thus get for the cumulative energy of $A$ [1, 5, 14] and for $B$, [9, 13, 14], as plotted in Figure 9.22a. The cumulative energy of a minimum-phase wavelet at any instant is always greater than that of any other wavelet with the same amplitude spectrum.

Figure 9.22a.  Cumulative energy.

Problem 9.22b

Calculate three-element Wiener inverse filters assuming the desired output is (i) [1, 0, 0] and (ii) [0, 1, 0], then apply the inverse filters to wavelets $A$ and $B$.

Solution

To write the normal equations in equation (9.22b) in explicit form, we give $i$ the values 0, 1, and 2 in succession. The result is

${\begin{aligned}\phi _{gg}(0)f_{0}+\phi _{gg}(-1)f_{1}+\phi _{gg}(-2)f_{2}=\phi _{gh}(0),\\\phi _{gg}(1)f_{0}+\phi _{gg}(0)f_{1}+\phi _{gg}(-1)f_{2}=\phi _{gh}(1),\\\phi _{gg}(2)f_{0}+\phi _{gg}(1)f_{1}+\phi _{gg}(0)f_{2}=\phi _{gh}(2).\\\end{aligned}}$

Shaping wavelet A [1,−2, 3] into [1, 0, 0]

To solve these equations, we need the values of $\phi _{gg}(\tau )$ and $\phi _{gh}(\tau )$. Using equations (9.8a) and (9.8e) we get

${\begin{aligned}\phi _{gg}(0)=1^{2}+2^{2}+3^{2}=14;\quad \phi _{gg}(1)=\phi _{gg}(-1)=-2-6=-8;\\\phi _{gg}(2)=\phi _{gg}(-2)=3;\quad \phi _{gh}(0)=1;\quad \phi _{gh}(1)=0;\quad \phi _{gh}(2)=0.\end{aligned}}$

Substituting these values in the normal equations gives

${\begin{aligned}14\;f_{0}-8f_{1}+3f_{2}&=1,\\-8f_{0}+14f_{1}-8f_{2}&=0,\\3f_{0}-8f_{1}+14f_{2}&=0.\end{aligned}}$

The solution of these equations can be obtained using Cramer’s rule [see Sheriff and Geldart, 1995, equation (15.3b)]. We first calculate the following determinants:

${\begin{aligned}\Delta =\left|{\begin{array}{ccc}{14}{-8}{3}\\{-8}{14}{-8}\\{3}{-8}{14}\end{array}}\right|\\=14(14\times 14-8\times 8)-8(-8\times 3+8\times 14)+3(8\times 8-14\times 3)=1210;\\\Delta _{0}=\left|{\begin{array}{ccc}{1}{-8}{3}\\{0}{14}{-8}\\{0}{-8}{14}\end{array}}\right|=(14\times 14-8\times 8)=132;\\\Delta _{1}=\left|{\begin{array}{ccc}{14}{1}{3}\\{-8}{0}{-8}\\{3}{0}{14}\end{array}}\right|=(-8\times 3+8\times 14)=88;\\\Delta _{2}=\left|{\begin{array}{ccc}{14}{-8}{1}\\{-8}{14}{0}\\{3}{-8}{0}\end{array}}\right|=(8\times 8-14\times 3)=22.\end{aligned}}$

Then, $f_{0}=\Delta _{0}/\Delta =132/1210=0.1091$, $f_{1}=88/1210=0.0727$, $f_{2}=22/1210=0.0182$. Applying this filter to $A$ gives

${\begin{aligned}f_{t}*A_{t}=\left[0.1091,\;0.0727,\;0.0182\left]*\right[1,\;-2,\;3\right]\\=0.1091,0.0727,0.0182\\-0.2182,\;-0.1454,\;-0.0364\\0.3273,0.2181,0,0546\\=0.1091,-0.1455,0.2001,0.1817,0.0546.\end{aligned}}$

Normalizing the wavelet to make the first element equal to +1, we get the wavelet

${\begin{aligned}\left[1,-1.3336,1.8341,1.6654,0.5005\right],\end{aligned}}$

which is far from $[1,\;0,\;0,\;0]$, the rms difference between the two wavelets being 1.43. Thus, it appears that we cannot shape the maximum-phase wavelet $A$ to get $\delta (t)$.

Shaping wavelet $B[3,\;-2,\;1]$ into $[1,\;0,\;0]$

Using wavelet $B$, we get the following values for $\phi _{gg}(\tau )$ and $\phi _{gh}(\tau )$:

${\begin{aligned}\phi _{gg}(0)=14;\phi _{gg}(1)=-8,\phi _{gg}(2)=3;\phi _{gh}(0)=3;\phi _{gh}(1)=0;\phi _{gh}(2)=0.\end{aligned}}$

The normal equations now become

${\begin{aligned}14\;f_{0}-8f_{1}+3f_{2}&=3,\\-8f_{0}+14f_{1}-8f_{2}&=0,\\3f_{0}-8f_{1}+14f_{2}&=0.\end{aligned}}$

Proceeding as before,

$\Delta =1210$ from calculation for wavelet A;

${\begin{aligned}\Delta _{0}=\left|{\begin{array}{ccc}{3}&{-8}&{3}\\{0}&{14}&{-8}\\{0}&{-8}&{14}\end{array}}\right|=3(14\times 14-8\times 8)=396;\\\Delta _{1}=\left|{\begin{array}{ccc}{14}&{3}&{3}\\{-8}&{0}&{-8}\\{3}&{0}&{14}\end{array}}\right|=-3(-8\times 14+8\times 3)=264;\\\Delta _{2}=\left|{\begin{array}{ccc}{14}&{-8}&{3}\\{-8}&{14}&{0}\\{3}&{-8}&{0}\end{array}}\right|=3(8\times 8-14\times 3)=66.\end{aligned}}$

The solution is

${\begin{aligned}f_{0}=396/1210=0.3273,f_{1}=264/1210=0.2182,f_{2}=66/1210=0.0545.\end{aligned}}$

Applying the filter [0.3273, 0.2182, 0.0545] gives

${\begin{aligned}f_{t}*B_{t}&=\left[0.3273,\;0.2182,\;0.0545\left]*\right[3,\;-2,\;1\right]\\&=0.9819,0.6546,0.1635\\&\qquad \qquad -0.6546,\;-0.4364,\;-0.1090\\&\qquad \qquad \qquad \qquad 0.3273,0.2182,0.0545\\&0.9819,0.,\quad 0.0544,0.1092,0.0545.\end{aligned}}$

Normalizing the first value to 1 gives

${\begin{aligned}f_{t}*g_{t}=[1,0,0.056,0.114,0.056],\end{aligned}}$

which is close to the desired wavelet of $[1,\;0,\;0]$ beause $B$ is minimum-phase. The rms difference from the desired output is 0.069.

Shaping wavelet $A[1,\;-2,\;3]$ into $[0,\;1,\;0]$

We use the values of $\phi _{gg}(\tau )$ from the previous calculations: $\phi _{gg}(0)=14$, $\phi _{gg}(1)=-8,\,\phi _{gg}(2)=3$; and $\phi _{gh}(0)=-2$, $\phi _{gh}(1)=1$, $\phi _{gh}(2)=0$. The normal equations thus become

${\begin{aligned}14\;f_{0}-8f_{1}+3f_{2}=-2,\\-8f_{0}+14f_{1}-8f_{2}=1,\\3f_{0}-8f_{1}+14f_{2}=0.\end{aligned}}$

We next calculate the determinants:

$\Delta =1210$ as before,

${\begin{aligned}\Delta _{0}=\left|{\begin{array}{ccc}{-2}&{-8}&{3}\\{1}&{14}&{-8}\\{0}&{-8}&{14}\end{array}}\right|=-176,\\\Delta _{1}=\left|{\begin{array}{ccc}{14}&{-2}&{3}\\{-8}&{1}&{-8}\\{3}&{0}&{14}\end{array}}\right|=11\\\Delta _{2}=\left|{\begin{array}{ccc}{14}&{-8}&{-2}\\{-8}&{14}&{1}\\{3}&{-8}&{0}\end{array}}\right|=44.\end{aligned}}$

Thus $f_{0}=-0.1455$, $f_{1}=0.0091$, $f_{2}=0.0364$ and

${\begin{aligned}f_{t}*A_{t}=\left[-0.1455,\;0.0091,\;0.0364\left]*\right[1,\;-2,\;3\right]=\end{aligned}}$

${\begin{aligned}-0.1455,0.0091,0.0364\\0.2910,-0.0182,-0.0728\\-0.4365,0.02730.1092\\-0.1455,0.3001,-0.4183,-0.0455,0.1092.\end{aligned}}$

Normalizing the second value to 1 gives

${\begin{aligned}f_{t}*A_{t}=[-0.485,\;1,\;-1.394,\;-0.152,\;0.364],\end{aligned}}$

and the rms difference from the desired output is 0.764. The result is poor because the input wavelet $A$ is maximum-phase, but it is much better than when we tried to make $A$ into [1,0,0].

Shaping wavelet $B[3,\;-2,\;1]$ into $[0,\;1,\;0]$

We have

${\begin{aligned}\phi _{gg}(0)=14,\;\phi _{gg}(1)=-8,\;\phi _{gg}(2)=3;\phi _{gh}(0)=-2,\\\phi _{gh}(1)=3,\;\phi _{gh}(2)=0.\end{aligned}}$

The normal equations are

${\begin{aligned}14\;f_{0}-8f_{1}+3f_{2}=-2,\\-8f_{0}+14f_{1}-8f_{2}=3,\\3f_{0}-8f_{1}+14f_{2}=0.\end{aligned}}$

We next calculate the following determinants:

$\Delta =1210$ as before;

${\begin{aligned}\Delta _{0}=\left|{\begin{array}{ccc}{-2}&{-8}&{3}\\{3}&{14}&{-8}\\{0}&{-8}&{14}\end{array}}\right|=0,\\\Delta _{1}=\left|{\begin{array}{ccc}{14}&{-2}&{3}\\{-8}&{3}&{-8}\\{3}&{0}&{14}\end{array}}\right|=385,\\\Delta _{2}=\left|{\begin{array}{ccc}{14}&{-8}&{-2}\\{-8}&{14}&{3}\\{3}&{-8}&{0}\end{array}}\right|=220.\end{aligned}}$

We now find that $f_{0}=0$, $f_{1}=-0.3182$, $f_{2}=0.1818$, and

${\begin{aligned}f_{t}*B_{t}=[0,0.3182,0.18181,0.0182]*[3,\;-2,\;1]\\=0,0.9546,0.5454\\0,\;-0.6364,\;-0.3636,\\0,0.3182,0.1818\\=0,0.9546,-0.0910,-0.0454,0.1818.\end{aligned}}$

Normalizing the second value to 1 gives

${\begin{aligned}f_{t}*B_{t}=[0,\;1,\;-0.0964,\;-0.0476,\;0.1904],\end{aligned}}$

and the rms difference from the desired output is 0.109. The result is not as good as for the previous case, but it is much better than for the wavelet $A$. Wiener filtering works best when input and desired wavelets have the same phase.

Figure 9.22b.  Original and Wiener-filtered waveforms.; first value = 1 Top row, waveform $A$; bottom row, waveform $B$.

The original and filtered wavelets are shown in Figure 9.22b. Reviewing the errors:

To make $A$ into [1, 0, 0]: rms error = 1.426;

To make $B$ into [1, 0, 0]: rms error = 0.069;

To make $A$ into [0, 1, 0]: rms error = 0.764;

To make $B$ into [0, 1, 0]: rms error = 0.109.

Continue reading

Previous section	Next section
Autocorrelation	Interpreting stacking velocity
Previous chapter	Next chapter
Reflection field methods	Geologic interpretation of reflection data
Table of Contents (book)

Also in this chapter

• Fourier series
• Space-domain convolution and vibroseis acquisition
• Fourier transforms of the unit impulse and boxcar
• Extension of the sampling theorem
• Alias filters
• The convolutional model
• Water reverberation filter
• Calculating crosscorrelation and autocorrelation
• Digital calculations
• Semblance
• Convolution and correlation calculations
• Properties of minimum-phase wavelets
• Phase of composite wavelets
• Tuning and waveshape
• Making a wavelet minimum-phase
• Zero-phase filtering of a minimum-phase wavelet
• Deconvolution methods
• Calculation of inverse filters
• Inverse filter to remove ghosting and recursive filtering
• Ghosting as a notch filter
• Autocorrelation
• Wiener (least-squares) inverse filters
• Interpreting stacking velocity
• Effect of local high-velocity body
• Apparent-velocity filtering
• Complex-trace analysis
• Kirchhoff migration
• Using an upward-traveling coordinate system
• Finite-difference migration
• Effect of migration on fault interpretation
• Derivative and integral operators
• Effects of normal-moveout (NMO) removal
• Weighted least-squares

External links

find literature about Wiener (least-squares) inverse filters