From left: Original image, blurred image, image
deblurred using Wiener deconvolution.
In
mathematics , Wiener deconvolution is an application of the
Wiener filter to the
noise problems inherent in
deconvolution . It works in the
frequency domain , attempting to minimize the impact of deconvolved noise at frequencies which have a poor
signal-to-noise ratio .
The Wiener deconvolution method has widespread use in
image deconvolution applications, as the frequency spectrum of most visual images is fairly well behaved and may be estimated easily.
Wiener deconvolution is named after
Norbert Wiener .
Definition
Given a system:
y
(
t
)
=
(
h
∗
x
)
(
t
)
+
n
(
t
)
{\displaystyle \ y(t)=(h*x)(t)+n(t)}
where
∗
{\displaystyle *}
denotes
convolution and:
x
(
t
)
{\displaystyle \ x(t)}
is some original signal (unknown) at time
t
{\displaystyle \ t}
.
h
(
t
)
{\displaystyle \ h(t)}
is the known
impulse response of a
linear time-invariant system
n
(
t
)
{\displaystyle \ n(t)}
is some unknown additive noise,
independent of
x
(
t
)
{\displaystyle \ x(t)}
y
(
t
)
{\displaystyle \ y(t)}
is our observed signal
Our goal is to find some
g
(
t
)
{\displaystyle \ g(t)}
so that we can estimate
x
(
t
)
{\displaystyle \ x(t)}
as follows:
x
^
(
t
)
=
(
g
∗
y
)
(
t
)
{\displaystyle \ {\hat {x}}(t)=(g*y)(t)}
where
x
^
(
t
)
{\displaystyle \ {\hat {x}}(t)}
is an estimate of
x
(
t
)
{\displaystyle \ x(t)}
that minimizes the
mean square error
ϵ
(
t
)
=
E
|
x
(
t
)
−
x
^
(
t
)
|
2
{\displaystyle \ \epsilon (t)=\mathbb {E} \left|x(t)-{\hat {x}}(t)\right|^{2}}
,
with
E
{\displaystyle \ \mathbb {E} }
denoting the
expectation .
The Wiener deconvolution filter provides such a
g
(
t
)
{\displaystyle \ g(t)}
. The filter is most easily described in the
frequency domain :
G
(
f
)
=
H
∗
(
f
)
S
(
f
)
|
H
(
f
)
|
2
S
(
f
)
+
N
(
f
)
{\displaystyle \ G(f)={\frac {H^{*}(f)S(f)}{|H(f)|^{2}S(f)+N(f)}}}
where:
G
(
f
)
{\displaystyle \ G(f)}
and
H
(
f
)
{\displaystyle \ H(f)}
are the
Fourier transforms of
g
(
t
)
{\displaystyle \ g(t)}
and
h
(
t
)
{\displaystyle \ h(t)}
,
S
(
f
)
=
E
|
X
(
f
)
|
2
{\displaystyle \ S(f)=\mathbb {E} |X(f)|^{2}}
is the mean
power spectral density of the original signal
x
(
t
)
{\displaystyle \ x(t)}
,
N
(
f
)
=
E
|
V
(
f
)
|
2
{\displaystyle \ N(f)=\mathbb {E} |V(f)|^{2}}
is the mean power spectral density of the noise
n
(
t
)
{\displaystyle \ n(t)}
,
X
(
f
)
{\displaystyle X(f)}
,
Y
(
f
)
{\displaystyle Y(f)}
, and
V
(
f
)
{\displaystyle V(f)}
are the Fourier transforms of
x
(
t
)
{\displaystyle x(t)}
, and
y
(
t
)
{\displaystyle y(t)}
, and
n
(
t
)
{\displaystyle n(t)}
, respectively,
the superscript
∗
{\displaystyle {}^{*}}
denotes
complex conjugation .
The filtering operation may either be carried out in the time-domain, as above, or in the frequency domain:
X
^
(
f
)
=
G
(
f
)
Y
(
f
)
{\displaystyle \ {\hat {X}}(f)=G(f)Y(f)}
and then performing an
inverse Fourier transform on
X
^
(
f
)
{\displaystyle \ {\hat {X}}(f)}
to obtain
x
^
(
t
)
{\displaystyle \ {\hat {x}}(t)}
.
Note that in the case of images, the arguments
t
{\displaystyle \ t}
and
f
{\displaystyle \ f}
above become two-dimensional; however the result is the same.
Interpretation
The operation of the Wiener filter becomes apparent when the filter equation above is rewritten:
G
(
f
)
=
1
H
(
f
)
1
1
+
1
/
(
|
H
(
f
)
|
2
S
N
R
(
f
)
)
{\displaystyle {\begin{aligned}G(f)&={\frac {1}{H(f)}}\left[{\frac {1}{1+1/(|H(f)|^{2}\mathrm {SNR} (f))}}\right]\end{aligned}}}
Here,
1
/
H
(
f
)
{\displaystyle \ 1/H(f)}
is the inverse of the original system,
S
N
R
(
f
)
=
S
(
f
)
/
N
(
f
)
{\displaystyle \ \mathrm {SNR} (f)=S(f)/N(f)}
is the
signal-to-noise ratio , and
|
H
(
f
)
|
2
S
N
R
(
f
)
{\displaystyle \ |H(f)|^{2}\mathrm {SNR} (f)}
is the ratio of the pure filtered signal to noise spectral density. When there is zero noise (i.e. infinite signal-to-noise), the term inside the square brackets equals 1, which means that the Wiener filter is simply the inverse of the system, as we might expect. However, as the noise at certain frequencies increases, the signal-to-noise ratio drops, so the term inside the square brackets also drops. This means that the Wiener filter attenuates frequencies according to their filtered signal-to-noise ratio.
The Wiener filter equation above requires us to know the spectral content of a typical image, and also that of the noise. Often, we do not have access to these exact quantities, but we may be in a situation where good estimates can be made. For instance, in the case of photographic images, the signal (the original image) typically has strong low frequencies and weak high frequencies, while in many cases the noise content will be relatively flat with frequency.
Derivation
As mentioned above, we want to produce an estimate of the original signal that minimizes the mean square error, which may be expressed:
ϵ
(
f
)
=
E
|
X
(
f
)
−
X
^
(
f
)
|
2
{\displaystyle \ \epsilon (f)=\mathbb {E} \left|X(f)-{\hat {X}}(f)\right|^{2}}
.
The equivalence to the previous definition of
ϵ
{\displaystyle \epsilon }
, can be derived using
Plancherel theorem or
Parseval's theorem for the
Fourier transform .
If we substitute in the expression for
X
^
(
f
)
{\displaystyle \ {\hat {X}}(f)}
, the above can be rearranged to
ϵ
(
f
)
=
E
|
X
(
f
)
−
G
(
f
)
Y
(
f
)
|
2
=
E
|
X
(
f
)
−
G
(
f
)
H
(
f
)
X
(
f
)
+
V
(
f
)
|
2
=
E
|
1
−
G
(
f
)
H
(
f
)
X
(
f
)
−
G
(
f
)
V
(
f
)
|
2
{\displaystyle {\begin{aligned}\epsilon (f)&=\mathbb {E} \left|X(f)-G(f)Y(f)\right|^{2}\\&=\mathbb {E} \left|X(f)-G(f)\left[H(f)X(f)+V(f)\right]\right|^{2}\\&=\mathbb {E} {\big |}\left[1-G(f)H(f)\right]X(f)-G(f)V(f){\big |}^{2}\end{aligned}}}
If we expand the quadratic, we get the following:
ϵ
(
f
)
=
1
−
G
(
f
)
H
(
f
)
1
−
G
(
f
)
H
(
f
)
∗
E
|
X
(
f
)
|
2
−
1
−
G
(
f
)
H
(
f
)
G
∗
(
f
)
E
{
X
(
f
)
V
∗
(
f
)
}
−
G
(
f
)
1
−
G
(
f
)
H
(
f
)
∗
E
{
V
(
f
)
X
∗
(
f
)
}
+
G
(
f
)
G
∗
(
f
)
E
|
V
(
f
)
|
2
{\displaystyle {\begin{aligned}\epsilon (f)&={\Big [}1-G(f)H(f){\Big ]}{\Big [}1-G(f)H(f){\Big ]}^{*}\,\mathbb {E} |X(f)|^{2}\\&{}-{\Big [}1-G(f)H(f){\Big ]}G^{*}(f)\,\mathbb {E} {\Big \{}X(f)V^{*}(f){\Big \}}\\&{}-G(f){\Big [}1-G(f)H(f){\Big ]}^{*}\,\mathbb {E} {\Big \{}V(f)X^{*}(f){\Big \}}\\&{}+G(f)G^{*}(f)\,\mathbb {E} |V(f)|^{2}\end{aligned}}}
However, we are assuming that the noise is independent of the signal, therefore:
E
{
X
(
f
)
V
∗
(
f
)
}
=
E
{
V
(
f
)
X
∗
(
f
)
}
=
0
{\displaystyle \ \mathbb {E} {\Big \{}X(f)V^{*}(f){\Big \}}=\mathbb {E} {\Big \{}V(f)X^{*}(f){\Big \}}=0}
Substituting the power spectral densities
S
(
f
)
{\displaystyle \ S(f)}
and
N
(
f
)
{\displaystyle \ N(f)}
, we have:
ϵ
(
f
)
=
1
−
G
(
f
)
H
(
f
)
1
−
G
(
f
)
H
(
f
)
∗
S
(
f
)
+
G
(
f
)
G
∗
(
f
)
N
(
f
)
{\displaystyle \epsilon (f)={\Big [}1-G(f)H(f){\Big ]}{\Big [}1-G(f)H(f){\Big ]}^{*}S(f)+G(f)G^{*}(f)N(f)}
To find the minimum error value, we calculate the
Wirtinger derivative with respect to
G
(
f
)
{\displaystyle \ G(f)}
and set it equal to zero.
d
ϵ
(
f
)
d
G
(
f
)
=
0
⇒
G
∗
(
f
)
N
(
f
)
−
H
(
f
)
1
−
G
(
f
)
H
(
f
)
∗
S
(
f
)
=
0
{\displaystyle \ {\frac {d\epsilon (f)}{dG(f)}}=0\Rightarrow G^{*}(f)N(f)-H(f){\Big [}1-G(f)H(f){\Big ]}^{*}S(f)=0}
This final equality can be rearranged to give the Wiener filter.
See also
References
Rafael Gonzalez, Richard Woods, and Steven Eddins. Digital Image Processing Using Matlab . Prentice Hall, 2003.
External links