Файл: Ersoy O.K. Diffraction, Fourier optics, and imaging (Wiley, 2006)(ISBN 0471238163)(427s) PEo .pdf

Скачать файл
Заказать решение

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 28.06.2024

Просмотров: 4787

Скачиваний: 1

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

AN ITERATIVE METHOD OF CONTRACTIONS FOR SIGNAL RECOVERY

221

the relation (14.4-11) can be written as

ð14:4-14Þ

ukþp

uk

1 ak a u1 u0

A sequence [an] is called a Cauchy sequence

if,

given any

2

> 0, there is a positive

integer Nð2Þ such that if n; m > Nð2Þ, it follows that

jan

amj <2

ð14:4-15Þ

A basic theorem of analysis is that a sequence converges to a fixed point if it is a Cauchy sequence.

The relation (14.4-14) shows that ½uk& is a Cauchy sequence. It follows that

k

Au

k ¼

Au

¼

u

ð

14:4-16

Þ

lim

!1

where u is the fixed point. In summary, if the mapping by A is strictly nonexpansive, the sequence created by ukþ1 ¼ Auk converges to the unique fixed point u .

14.5 AN ITERATIVE METHOD OF CONTRACTIONS FOR SIGNAL RECOVERY

Assume that there is a measured signal v, which is distorted in some fashion. From v, a more accurate signal is desired to be obtained through iterations. The iterated signal at iteration k is uk, k being an integer for counting iterations. As k increases, uk approaches the recovered signal.

The method to be discussed for this purpose is also known as the method of constrained iterative signal restoration [Schafer]. It is based on the following

iterative equation:

ukþ1 ¼ Fuk

ð14:5-1Þ

where uks are the successive approximations to the recovered signal, and F is an operator. F is obtained both from a distortion operator and a priori knowledge in the form of certain constraints. In general, F is not unique. A simple technique of estimating F is given below.

Prior knowledge can be expressed by a constraint operator:

Cu ! u

ð14:5-2Þ

For example, with a multidimensional signal u, C½u& may be given as

8 u½n&

Cu½n& ¼

<

0

: a

0 u½n& a

u½n& < 0 ð14:5-3Þ u½n& > a

222 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

As another example, if u(t) is an analog signal bandlimited by a system to frequencies below fc, Cu(t) can be written as

1

2 f t

CuðtÞ ¼

ð

uðtÞ

sin p cð

tÞ

dt

ð14:5-4Þ

pðt tÞ

1

Every estimate uk of u should be constrained by C. Cuk is interpreted as an approximation to ukþ1. The result Cuk is further processed with a distortion operator D, if any, to yield DCuk. Interpreting DCuk as an estimate of v, the output errorv is defined as

v ¼ v DCuk

ð14:5-5Þ

The recovered signal error between the (k þ 1)th and kth iterations can be defined as

u ¼ ukþ1 Cuk

ð14:5-6Þ

How does u and v relate? The simplest assumption is that they are proportional, say, u ¼ l v. Then,

Fuk ¼ ukþ1 ¼ Cuk þ u ¼ Cuk þ lðv DCukÞ

ð14:5-7Þ

This can also be written as

Fuk ¼ lv þ Guk

ð14:5-8Þ

where

G ¼ ðI lDÞC

ð14:5-9Þ

and I is the identity operator. The initial approximation to u can be chosen as u0 ¼ lv.

If F is chosen contractive, uk converges to a fixed-point u . Equation (14.5-8) shows that

jFuk Funj ¼ jGuk Gunj

ð14:5-10Þ

Thus, if F is contractive, so is G. With u0 ¼ lv, Eq. (14.5-8) becomes

Fuk ¼ ukþ1 ¼ u0 þ Guk

ð14:5-11Þ
u00k Þ.

ITERATIVE CONSTRAINED DECONVOLUTION

223

The procedure at the kth iteration according to Eq. (14.5-11) is as follows:

1.Apply the constraint operator C on uk, obtaining u0k.

2.Apply the distortion operator D on u0k, obtaining u00k .

3. Compute ukþ1 as u0 þ lðuk

An example of this method is given below for deconvolution.

14.6ITERATIVE CONSTRAINED DECONVOLUTION

Given v and h, deconvolution involves determining u when v is the convolution of u and h:

v ¼ h u

ð14:6-1Þ

When the iterations start with u0 ¼ lv, Eq. (14.5-8) can be written as

ukþ1 ¼ lv þ ðd lhÞ Cuk

ð14:6-2Þ

where the constraint operator is simply the identity operator for the time being. The FT of Eq. (14.6-2) yields

Ukþ1 ¼ lV þ Uk lHUk

ð14:6-3Þ

where the capital letters indicate the FTs of the corresponding lowercase-lettered variables in Eq. (14.6-2)

Eq. (14.6-3) is a first-order difference equation in the index k. It can be written as

Ukþ1

ð1 lHÞUk ¼ lV

ð14:6-4Þ

whose solution is

V

lH&kþ1u½k&

Uk ¼

½1

ð14:6-5Þ

H

where u½k& is the 1-D unit step sequence. As k ! 1, the solution becomes

U1 ¼

V

ð14:6-6Þ

H

if

j1 lHj < 1

ð14:6-7Þ

Equation (14.6-6) is the same result obtained with the inverse filter. The advantage of the interactive procedure is that stopping the iterations after a finite number of steps may yield a better result than the actual inverse filter.

224 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

Equation (14.6-7) shows that the procedure will break down if H equals 0 at some frequencies. Equation (14.6-7) can also be written as

Re½H& > 0

ð14:6-8Þ

In order to avoid testing the validity of Eq. (14.6-8) at all frequencies, a better procedure involves convolving both sides of Eq. (14.6-1) with inverted h (h ½ m& in 1-D and h ½ m; n& in 2-D with discrete signals). Then, we get

v0 ¼ x h0

ð14:6-9Þ

where h0 is a new impulse response sequence whose transfer function is given by

H0 ¼ jHj2

ð14:6-10Þ

and v0 is the convolution of v with the inverted h .

In solving Eq. (14.6-9) by the iterative method, Eq. (14.6-8), namely, Re½H0& > 0 is automatically satisfied provided that H 6¼ 0 at any frequency.

Since the iterative solution converges toward the inverse filter solution, it may also end up in an undesirable behavior. For example, a unique result may not be obtained if H ! 0. The use of a constraint operator often circumvents this problem. For example, u may be assumed to be positive over a finite support. In 2-D, the

constraint operator for this purpose is given by

8

uk½m; n& M1 m < M2; N1

n

N2

Cuk m; n

uk m; n 0

14:6-11

½

:

otherwise½ &

ð

Þ

& ¼ <

0

Letting uk0 ¼ Cuk, Eq. (14.6-2) in the DTFT domain becomes

U0

¼ l

V

þ

U0

l

HU0

ð

14:6-12

Þ

kþ1

k

k

The iterations defined by Eq. (14.6-12) can be done by using the DFT. However, u0k ¼ Cuk is to be implemented in the signal domain. The resulting block diagram is shown in Figure 14.4.

Figure 14.4. The block diagram for iterative constrained deconvolution utilizing the Fourier transform.

METHOD OF PROJECTIONS

225

In many applications the impulse response is shift variant. The iterative method discussed above can still be used in the time or space domain, but the algebraic equations in the transform domain are no longer valid.

EXAMPLE 14.2 A 2-D Gaussian blurring impulse response sequence is given by

h

m; n

& ¼

exp

m2 þ n2

½

100

The image sequence is assumed to be

u½m; n& ¼ ½dðm 24Þdðn 32Þ þ dðm 34Þdðn 32Þ&

(a)Determine the blurred image v ¼ u h. Visualize h and v.

(b)Compute the deconvolved image with l ¼ 2, 65 iterations, and a positivity constraint.

Solution: (a) The convolution of u and h is given by [Schafer, Mersereau, Richards]

v m; n

exp

ðm 24Þ2 þ ðn 32Þ2

exp

ðm 34Þ2 þ ðn

32Þ2

& ¼

100

# þ

100

#

½

"

"

Plots of h and v are shown in Figure 14.5(a) and (b), respectively.

(b)Figure 14.5(c) shows the deconvolution results after 65 iterations.

14.7METHOD OF PROJECTIONS

A subset of contractions is projections. An operator P is called a projection operator or projector if the mapping Pu satisfies

j

Pu

u

j ¼

v S0

j

v

u

j

ð

14

:

7-1

Þ

inf

2

where u 2 S, and S0 is a subset of S. In the context of signal recovery, S0 is the subset to which the desired signal belongs to. g ¼ Pu is called the projection of u onto S0. This is interpreted as the projection operator P generating a vector g 2 S0 which is closest to u 2 S. G is unique if S0 is a convex set. A very brief review of a convex set is given below.

Let u1 and u2 belong to S0. S0 is convex, if for all l, 0 l 1, lu1 þ ð1 lÞu2 also belongs to S0. A convex set and a nonconvex set are visualized in Figure 14.6.

It can be shown that if S0 is convex, then P is contractive [Youla]. Hence, it is desirable that S0 is indeed convex. If so, the resulting method is called projections onto convex sets. We will consider this case first.

Смотрите также файлы