# Wiener (least-squares) inverse filters

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

## Problem 9.22a

Plot cumulative energy as a function of time for wavelets ${\displaystyle A=[1,\;-2,\;3]}$ and ${\displaystyle 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 ${\displaystyle g_{t}}$ to denote the original signal and ${\displaystyle h_{t}}$ the desired signal, the sum of the errors squared is

 {\displaystyle {\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 ${\displaystyle f_{i}}$, we vary the coefficients ${\displaystyle f_{i}}$ to minimize ${\displaystyle E}$. We differentiate ${\displaystyle E}$ with respect to each of the elements ${\displaystyle f_{i}}$ and equate the derivatives to zero. This gives ${\displaystyle (n+1)}$ equations which can be solved for the ${\displaystyle (n+1)}$ elements of the filter. Differentiating ${\displaystyle E}$, we get

{\displaystyle {\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 ${\displaystyle f_{t}*g_{t}}$ with a summation. The derivative now becomes

{\displaystyle {\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 ${\displaystyle f_{t}}$ appears, that is, ${\displaystyle t-k=i}$, so ${\displaystyle k=t-i}$. Substituting this result, we get

{\displaystyle {\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 ${\displaystyle \phi _{gh}}$(${\displaystyle i}$) from equations (9.8a) and (9.8b). Interchanging the order of summation in the right-hand term, letting ${\displaystyle j=t-k}$, and summing over ${\displaystyle j}$, the equation becomes

{\displaystyle {\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

 {\displaystyle {\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 ${\displaystyle A}$ and ${\displaystyle B}$ are ${\displaystyle (1-2z+3z^{2})}$ and ${\displaystyle (3-2z+z^{2})}$, so the roots of ${\displaystyle A}$ are ${\displaystyle {\frac {1}{3}}(1\pm j{\sqrt {2}})}$ and those of ${\displaystyle B}$ are ${\displaystyle (1\pm j{\sqrt {2}})}$, the magnitudes of the roots being ${\displaystyle 1/{\sqrt {3}}}$ and ${\displaystyle {\sqrt {3}}}$. Thus ${\displaystyle A}$ is maximum-phase and ${\displaystyle 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 ${\displaystyle A}$ [1, 5, 14] and for ${\displaystyle 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 ${\displaystyle A}$ and ${\displaystyle B}$.

### Solution

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

{\displaystyle {\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 ${\displaystyle \phi _{gg}(\tau )}$ and ${\displaystyle \phi _{gh}(\tau )}$. Using equations (9.8a) and (9.8e) we get

{\displaystyle {\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

{\displaystyle {\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:

{\displaystyle {\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, ${\displaystyle f_{0}=\Delta _{0}/\Delta =132/1210=0.1091}$, ${\displaystyle f_{1}=88/1210=0.0727}$, ${\displaystyle f_{2}=22/1210=0.0182}$. Applying this filter to ${\displaystyle A}$ gives

{\displaystyle {\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

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

which is far from ${\displaystyle [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 ${\displaystyle A}$ to get ${\displaystyle \delta (t)}$.

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

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

{\displaystyle {\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

{\displaystyle {\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,

${\displaystyle \Delta =1210}$ from calculation for wavelet A;

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

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

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

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

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

{\displaystyle {\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:

${\displaystyle \Delta =1210}$ as before,

{\displaystyle {\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 ${\displaystyle f_{0}=-0.1455}$, ${\displaystyle f_{1}=0.0091}$, ${\displaystyle f_{2}=0.0364}$ and

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

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle A}$ is maximum-phase, but it is much better than when we tried to make ${\displaystyle A}$ into [1,0,0].

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

We have

{\displaystyle {\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

{\displaystyle {\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:

${\displaystyle \Delta =1210}$ as before;

{\displaystyle {\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 ${\displaystyle f_{0}=0}$, ${\displaystyle f_{1}=-0.3182}$, ${\displaystyle f_{2}=0.1818}$, and

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle 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 ${\displaystyle A}$; bottom row, waveform ${\displaystyle B}$.

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

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

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

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

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