Derivation of Wiener Filter 1

Preliminaries

Before we begin the derivation let us get some definitions in place. Consider the system below where an image v(m, n) of dimension M xN is fed into a filter whose spatial domain representation is h(m, n) and the output is u(m, n). We will frequently use the terms defined in the following subsections in the derivation which begins with section 2

Figure 1: Block diagram representing Image Filtering

1.1

Convolution

In the case of linear and time-invariant systems, convolution with the impulse response of the filter captures the relation between the input and output signals. u(m, n) = v(m, n) ∗ h(m, n) =

=

M X N X m0 =1 n0 =1 M X N X m0 =1

(1)

h(m − m0 , n − n0 )v(m0 , n0 )

(2)

v(m − m0 , n − n0 )h(m0 , n0 )

(3)

n0 =1

Note that * denotes a convolution between two signals not the point wise multiplication. Another concept that we would heavily rely upon is that the FT of a convolution of two is simply the product of the FT of the original signals. Hence, U (ω1 , ω2 ) = V (ω1 , ω2 )H(ω1 , ω2 ) (4) Also there is a relation between the Power Spectral Densities of the input and output which goes as follows Suu (ω1 , ω2 ) = |H(ω1 , ω2 )|2 Svv (ω1 , ω2 ) 1

(5)

1.2

Cross-Correlation Ruv (k, l) = E[u(m, n)v(m − k, n − l)]

(6)

When both the signals in the above equation are the same then this operation is called Autocorrelation. Remark: DSP enthusiasts may note that the above definitions are gonna hold only if the following conditions are satisfied: 1) u and v are jointly stationary, meaning their joint probability density function does not change with time 2) The filter h is space invariant, meaning that if the input image is translated in space, then so shall the output image.

1.3

Relationship between Power Spectral Density and Autocorrelation

Finally the last ingredient to be collected. Presumably you have seen the definition of PSD in terms of the Fourier transform of the signal on my blog. But, that’s not the only way to compute it. PSD can also be written in terms of the Autocorrelation function. It is in fact the Fourier Transform of the Autocorrelation function. Do keep this in mind as we are ready to go ahead with the derivation.

2

Derivation

So we are going to derive the equation in its canonical form, i.e. when the input image also passes through another known filter before the addition of noise. For example, when we try to take a photograph under poor illumination with low shutter speed, then there is a possibility of motion blur. Also the image is noisy because of the limitations of the image acquisition sensor (charge-coupled device in digital cameras). In the world of signal processing the system is represented as follows: In figure 2 h(m, n) is the impulse response or the point spread function (PSF) for the phenomenon of motion blur , meaning that this is the image of a point that would be formed by the camera if there was no noise. While in general there are different techniques for estimating this PSF but in our 2

Figure 2: Block diagram representing the scenario for Wiener Filtering case of noise filtering h(m, n) is simply a delta function since convolution of a signal with a delta function gives the same signal. η(m, n) is the noise which is added to the blurred image. Our job is to find the filter, defined by its impulse response g(m, n), which acts on the noisy signal and produces an output uˆ(m, n) which is as similar as possible to u(m, n), i.e. σe2 = E[{u(m, n) − uˆ(m, n)}2 ]

(7)

Now according to A.K.Jain’s ”Fundamentals of Digital Image Processing” this Mean Squared Error attains its minimum value when the orthogonality condition holds i.e. E[{u(m, n) − uˆ(m, n)}v(m0 , n0 )] = 0

(8) ∞ X

Ruv (m − m0 , n − n0 ) =

g(m − k, n − l)Rvv (k − m0 , l − n0 )

k,l=−∞

(9) If you have some difficulty seeing this then simply write uˆ(m, n) as a convolution of v(m, n) and g(m, n). Then multiply it with v(m, n) and take expectation using the formula of autocorrelation. A little jugglery with the indices would tell us that equation 9 is equivalent to Ruv (m, n) =

∞ X

g(m − k, n − l)Rvv (k, l)

(10)

k,l=−∞

Taking DTFT on both sides of the equation 10 gives Suv (ω1 , ω2 ) = G(ω1 , ω2 )Svv (ω1 , ω2 ) Suv (ω1 , ω2 ) G(ω1 , ω2 ) = Svv (ω1 , ω2 ) 3

(11) (12)

So if we find Suv and Svv then we are done, ain’t we !! Observe that ∞ X

v(m, n) =

h(m − k, n − l)u(k, l) + η(m, n)

(13)

k,l=−∞

Let us find the PSD on both sides of the equation Svv (ω1 , ω2 ) = |H(ω1 , ω2 )|2 Suu (ω1 , ω2 ) + Sηη (ω1 , ω2 )

(14)

Now that Svv is out of the way, let us try to find Suv which is FT of Ruv . Ruv (m − m0 , n − n0 ) = E[u(m, n)v(m0 , n0 )]; ∞ X = h(m0 − k, n0 − l)Ruu (m − k, n − l)

(15) (16)

k,l=−∞

Do the jugglery with the indices again and we get Ruv (m, n) =

∞ X

h(−(m − k), −(n − l))Ruu (k, l)

(17)

k,l=−∞

= h(−m, −n) ∗ Ruu (m, n)

(18)

Take our favorite FT on both sides of this equation and you get Suv (ω1 , ω2 ) = H ∗ (ω1 , ω2 )Suu (ω1 , ω2 )

(19)

Using equations 12, 14 and 19 we get the canonical form of Wiener Filter in frequency domain G(ω1 , ω2 ) =

H ∗ (ω1 , ω2 )Suu (ω1 , ω2 ) |H(ω1 , ω2 )|2 Suu (ω1 , ω2 ) + Sηη (ω1 , ω2 )

(20)

For the purpose of noise filtering we can set H(ω1 , ω2 ) = 1∀ ω1 , ω2 and we get our result...yay!! G(ω1 , ω2 ) =

Suu (ω1 , ω2 ) Suu (ω1 , ω2 ) + Sηη (ω1 , ω2 )

4

(21)