Wiener (least-squares) inverse filters

From SEG Wiki
Jump to navigation Jump to search
ADVERTISEMENT

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): \begin{align} 14\; f_{0} -8f_{1} +3f_{2} =-2,\\ -8f_{0} +14f_{1} -8f_{2} =3,\\ 3f_{0} -8f_{1} +14f_{2} =0. \end{align}

We next calculate the following determinants:

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): \Delta =1210 as before;

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): \begin{align} \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{align}

We now find that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): f_{0} =0 , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): f_{1} =-0.3182 , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): f_{2} =0.1818 , and

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): \begin{align} 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{align}

Normalizing the second value to 1 gives

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): \begin{align} f_{t} *B_{t} =[0,\; 1,\; -0.0964,\; -0.0476,\; 0.1904], \end{align}

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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): B .

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

To make Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): A into [1, 0, 0]: rms error = 1.426;

To make Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): B into [1, 0, 0]: rms error = 0.069;

To make Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): A into [0, 1, 0]: rms error = 0.764;

To make Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): 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

External links

find literature about
Wiener (least-squares) inverse filters
Retrieved from "http://wiki.seg.org/w/index.php?title=Wiener_(least-squares)_inverse_filters&oldid=133092"