Analysis and recovery of multidimensional signals from irregular samples using nonlinear and iterative techniques

SIGNAL PROCESSING Signal Processing 36 (1994) 13-30

ELSEVIER

Analysis and recovery of multidimensional signals from irregular samples using nonlinear and iterative techniques Farokh A . Marvasti*, Chuande Liu t , Gil Adams ', Cornmunications Research Group, Department of Electronics and Electrical Engineering . King's College London, Uni, •e rsitr Strand WC2R 2LS . London, UK

of London .

Received 20 November 1992 : revised 16 April 1993

Abstract A recurring problem in signal processing is reconstruction based upon imperfect information . In this paper we offer two methods . The first one is a simple nonlinear method, based upon spectral analysis, to reconstruct two-dimensional (2-D) band-limited signals from nonuniform samples . The second method is an iterative technique ; by using an intelligent initial condition, we get faster convergence than a simple successive approximation . These methods and analyses are extendable to n-dimensional signals . The feasibility of this procedure is demonstrated through examples of both black-and-white and colour image restoration when some of the pixels are lost . Zusammenfassung Ein hiiufig wiederkehrendes Problem in der Signalverarbeitung is( die Rekonstruktion anhand von nicht perfekter Information . In dieser Arbeit werden zwei Methoden angegeben . Die erste ist eine einfache nichtlineare Methode basierend auf einer Spektralanalyse zur Rekonstruktion von zweidimensionalen (2-D) bandbegrenzten Signalen aus nicht gleichformigen Abtastwerten . Die zweile Methode ist eine iterative Technik : Dutch die Verwendung einer intelligenten Initialbedingung bekommen wir eine schnellere Konvergenz als bet cinfacher sukzessiver Approximation . Diese Methoden and Analysen werden auf n-dimensionale Signale erweitert . Die Wirksamkeit dieser Prozedur wird durch Beispiele von Schwarz weiB- and Farbbild-Rekonstruktionen demonstriert, wenn einige Pixel verlorengegangen sind . Resume Un problcme frequent en traitement du signal est la reconstruction basee sur une information imparfaite . Dans cet article, nous proposons deux methodes . La premiere est une methode non-lineaire simple basee sur ('analyse du spectre pour reconstruire des signaux bi-dimensionnels a bandes limites a partir d echantillons non uniformes . La seconde methode est une technique iterative ; en utilisant une condition initiale intelligente, nous obtenons une convergence plus

-Corresponding author . tPrcsently at Robotics Inc ., Skoki, IL, USA. tEdison Electric Co ., Chicago . IL, USA . 0165-1684/94,/$7-00 L, 1994 Elsevier Science SSDI0165-16840 31E0062-P

B .V . All

rights reserved

P. A . Marvasti et al. Signal Processing 36 (1994 ; 13-30

rapide qu'avec une simple approximation successive . Ces mclhodes et ces analyses peuvent titre etendues aux signaux n dimensions . La faisabilite de cette procedure est demontree par des exemples de restoration d'images noir et blanc on couleur lorsque des pixels sont perdus . a

Key words : Nonuniform samples ; Spectral analysis of irregular samples ; Image recovery from missing samples

l . Introduction

ative methods for the recovery of 1-D signals are given in [10] . The purpose of this paper is to extend the spectral analysis and recovery methods from irregular samples to 2-D signals, and apply the developed techniques to images with missing samples . In the next section, we analyze the spectrum of 2-D irregularly sampled signals, leading to the nonlinear signal recovery method . In Section 3 the nonlinear signal recovery and an iterative reconstruction technique are discussed . In Section 4, experimental results are given for 2-D signals - comparing the nonlinear technique, in terms of signal-to-noise ratios, to the iterative technique and other standard techniques . Visual results for pure 2-D signals (black-and-white) and for color images are presented for inspection.

Irregular sampling is more the rule than the exception for many types of band-limited signals, particularly where the samples deviate from uniform positions only by a finite amount . Examples include estimation of antenna radiation patterns from far fields [15], radio astronomy and image processing [3] to name a few . Another special example is the case of lost pixels in a recording channel or a fading transmission medium [5] . The remaining pixels form a set of nonuniform samples . The spectral analysis and recovery from irregular samples for one-dimensional (1-D) signals is given by the first author in a recent book [8], and application of this analysis for the recovery of speech signals with missing samples is given in [1I] . Iter-

y

2T2 0

0

o

o

0

°

°

0

T2 0

-3T1

°

0

2T

-

0

11

2T1

3T1

'~

0 0

0 o

T2 °

0

0 °

-2T2

Fig. 1 . A general nonuniform sampling position.



F. .4 . Marvasti et al . Signal Processing .36 (1994) 13 30

15

2. Spectral analysis of irregular samples

JP (x,

y) = L 16 (x

- xnm,

V - Ynm)'

(2)

m n

In this section, we would like to analyze the spectrum of nonuniform samples that deviate from uniform positions by a finite amount . The derivation of the spectrum of nonuniform samples leads to a practical method to recover a 2-D signal from the nonuniform samples . We model the set of 2-D nonuniform samples by

f. (x, Y)

(xnm,

Yr .) d(x --

xnm, Y - Ynm),

Now, define the following functions : 91(x, y)

= x - n Tt -

B, (x, Y),

such that 9, (xnm, Ynm) = 9 2 (x, ., y. .) = 0 for any n and m . 0,(x,y) and B 2 (x,y) are any 2-D functions such that

(1)

01(xnm, Ynm)

where f(x,y) is a band-limited 2-D signal . The ir.} is assumed to deviate from the . regular set {xnm, Y uniform set { n T, ,mT 2 } as seen in Figs . 1-3 . Fig . I is the most general case . Fig . 2 represents a nonuniform sampling set that is represented by { .xn,ynm} . Fig . 3 is a special case of a nonuniform sampling set denoted by {xn,ym} . To analyze the spectrum of (1), we first define the comb function as [20]

=

xnm -

02(ynm, ynm)

n T, , mT2 .

b(x -

xnm,Y - Ynm)=d[9t(x,Y),92(x,y)]IJI,

where J is the Jacobian defined by

12

Y21

y21

x-2

(3b)

Eq . (3b) can be interpreted as deviation of nonuniform points from uniform positions . Now, we have the following identity :

y

y22

(3 a)

92(x,Y)=Y - mT2 -62(x,v),

x1

Y2-1

Fig . 2. A special case where the nonuniform samples lie on vertical lines .

(4)



P. A . Marvasti et al. l Signal Processing 36 (1994) 13-30

16

Y

m

r

n

a

Fig . 3 . A

a9i

special case where the nonuniform samples lie on vertical and horizontal lines .

x91092 ra91 0'92 J= ax ax (5) 091 092 ax ay ay ax ' ay ay (It is assumed that J(x,,„, ynm ) # 0 .) Eq . (4) holds by the definition of distribution functions, i .e., f f

notation x, = x - 0 1 (x, y) and 2 Yi = y - 0 (x, y), and using Fourier series expanUsing

092

6(x - x, ., y - yan ) dx dy =

f f

the

sion, we get

fp(x,y)=-. T2

e ,Ili2n/TOx,+p2njT,iy,,

(9)

t i

b (91 , 92) d91 d92

or

=ffa[x - nT, x I J I

- 01(x,y),v - mT2 - 02(x,Y)7

dx dy .

(6) .fp(x, Y)

Substituting (3) in (5), we get a0, 1J1=

(1-

ax)

x0 2

I

ay)

a0,

00 2

ay

ax

(7)

Now by using (3) and (4), the comb function (2) can be written as f,(x,Y) = ~~ JIS(x - nT1 - 01(x,Y),Y - mT2 - 02(x,Y)) . (8)

= T

1 + 21 E cos T2 {

T

t [x - 01(x, Y)]

2nI x I+2E cos-[y-B 2 (x,y)] T2 =1 {

(10)

Eq . (10) is equivalent to the sum of phase modulated signals in two dimensions . Therefore, the Fourier spectrum of fp (x, y) consists of a lowpass component I JI/( T, T2 ) and band-pass or high-pass components around carrier frequencies (2nn/T,,2am/T2 ) as shown in Fig . 4 .



F.A . Maruasli ei al . 1 Signal Processing 36 (1994) 13-30

17

V

J(u,v)

u

Fig 4 . The spectrum of nonuniform samples (comb function) .

The nonuniform samples (1) are derived by multiplying (9) by f (x, y), i .c .,

Notice that for the case of uniform samples, (13) reduces to the following familiar result :

Mx, Y) =f (x, Y)fv(x, Y)

(All, )

= 0(x) = 0,

J=1,

=f (x, y) T T

1

(

f(x .Y)=TI

x

1[+

2m 2 cos -[x - 0, (x,y)) ~,

21d cos -[y-0 2 (x, y)]

x 1+2 1 =,

x

T,

I

Tz

d0,(x)1~ c'02(x,Y)1 1-- a dx J

where 0,(x,y)=01(x) and 0 2 (x,y) remains the same . For the special case represented in Fig . 3, the Jacobian reduces to

IJ1=I[1-0(x)7[1-y(Y)]], where 0, (x, y) = 0(x) and 0 2 (x, y) = q5(y) .

T2

I + 2 I

m=,

Y

2
cos-x

n=,

T,

(11)

For the special case shown in Fig . 2, the Jacobian given in (7) can further be reduced to J1 =

cos

I+2

The spectrum of the above comb function reduces to a periodic function of impulses, and the spectrum of f(x,y) becomes a periodic function of (I/(TIT2 ))Flu, v), the Fourier transform of f(x, y) . The above analysis can be generalized to an n-dimensional signal . A specific application might be in the analysis of a 3-D object or a TV video signal regarded as a 3-D signal (in space and time). To avoid unwieldy equations, we use vector notation, i .e ., the comb function can be written as

(13)

Ib[x - x,7=Y,

dg(x)

dx

6[g(x)7,

(14)



F.A . Marvasri et al. / Signal Processing 36 (1994, 13-30 where x = (x„x2, . . .,xn), xa are irregular samples and dg(x)/dx is the Jacobian defined by 591 OX,

592 5x,

5g,, Sx,

5g, 5g2 5xn OX,

59 5xn

J =

Spectral analysis reveals that the information about the original signal lies in the low-pass component of the ideal nonuniform samples (11) ; i.e ., low-pass filtering the nonuniform samples (11), we get

dg(x) dx

P[f(x,y)]~Zf(x'y) Iji T, T2'

where g=(g1,92, . . .,ga) is similar to (3), viz ., 9t(x,, . . .'xa)=xt-m1T,-0,(x,, . .-,xa), (16) 9a(xi, . . . , xn) = xa - ninT - 6a(x,, . . . . xa),

(17)

The Fourier series expansion of (14) is Y,d[x - x,] = l expJ T m[x - 9(x)] .(18) m „ As before, provided that there is no overlap (similar to Fig. 4), the n-dimensional low-pass filtering yields dg(x) dx The low-pass version of nonuniform samples of an n-dimensional signal f(x) is dg(x) dx

(21)

where P is a low-pass operator. The Jacobian in (21) can be evaluated by low-pass filtering the comb function (10) - as shown in Fig . 4. The result is P[JP(x,y)] s T1T2'

in short g(x) =x-mdiag[T]-0(x) .

3 .1 . A nonlinear recovery method

(22)

The approximation in (21) and (22) is due to the spectral overlap of phase modulated components as shown in Fig. 4. A division of (21) by (22) yields a good approximation of the original signal f(x, y) (see Fig . 5) . The division is possible if the Jacobian function is not equal to zero . This condition puts a limitation on the sampling point deviations . A sufficient condition is that the maximum values of B, (x, y) and 02(x, y) in (3b) - which are always greater than or equal to the nonuniform deviations in (3a) - be less than T, /2n and T2/2n, respectively . For the special case represented in Fig . 3, a sufficient condition is Ixn - nTi I ~ OMa . 6 T1 ym - mT2I 1< 0n+a . \ TZ 7C n

(20)

Therefore, a division of (20) by (19) (similar to Fig . 5) yields f (x) . This spectral analysis leads to a practical reconstruction method, which is our next topic of discussion .

For the proof, see Appendix A . For other sufficient conditions, see [9] . Fig . 5 is an implementation of the above technique, namely the division of (21) by (22) . The input

Nonuniform samples 3 . Practical reconstruction techniques In this section we shall discuss two practical reconstruction methods for image recovery from a set of nonuniform samples, namely a nonlinear method and an iterative technique .

f(x,y)

I J I ftx,vl 2 D LP

2-0 LPF Fig. 5 . The block diagram of the division technique .

F.A . .Mari asmi cr at .

Signal Processing 36 (1994) 1 r 30

to the system represented in Fig . 5 is f (.x, v), which is the set of nonuniform samples . After low-pass filtering in the upper branch, we get the result shown in (21) . 'The blocks shown as a hard limiter and an absolute value in tie lower branch of Fig . 5 are equivalent to taking the impulsive samples with unit area at the nonuniform positions (the comb function f,(x,v) represented in (10)). )''ut :r low-pass filtering, we get (22) . A division of the upper branch by the lower branch is equivalent to dividing (21) by (22).

-pass ~9gyrigq

PS, X ,

,

pp

crcpe pss

P9x, S

3 .2 . An iterative recovery method

An iterative method has been developed to recover the nonuniform samples . The iteration is as follows : Fig . I,

I (x,

Y) =

The hlock diagram

of

the itcrative method .

I'm MX . Y), above iterative technique is an extension to the 1-D case [10] . The proof of the convergence is, however . more involved and is given in Appendix B .

where

f t (X, Y) = k, [PS ./ (x, Y) -- PS .f (x, v)] +

„(x

v),

(23) where P and S are, respectively, the low-pass filtering and the nonuniform sampling process for a 2-D signal . The proof of convergence is given in Appendix B . This iterative technique is shown in Fig . 6, and is explained as follows . The set of nonuniform samples (Sx) is first low-pass filtered and is shown by the block P. The output is denoted by x, . The second iteration is to pass x, through the same nonuniform sampling process S and low-pass filtering P : the result is PS .x, . The error signal .x, - PSx, is added to x, to get the second iteration x,- . Fig . 6 shows a modular implementation where the ill iteration x ; is equal to the error signal . „- PSx t ; „ plus the first iteration x, . Note ; x( that the above iteration is not identical to the iterative methods discussed in [17, 18] . Reference [17] uses smeared pulses as opposed to ideal impulsive samples, and [18] uses the concept of projection onto convex sets (POCS) for the ideal impulsive samples and the iteration is based on sample by sample (for more details see [2 . 21]) . The

4 . Simulation results To simulate the nonuniform samples . we consider the important special case where in a uniform grid of points, some of the pixels are randomly lost . The remaining set of pixels form a set of nonuniform samples . In the following, we simulate the nonlinear method, the iterative method, and compare it to other standard techniques such as lowpass filtering and mean filtering .'

4 .1 . The nonlinear method

the nonlinear method has been tested on a Sun Station for a 256 x 256 band-limited image as shown in Fig . 7(a), sampled at twice the Nyquist rate . To create nonuniform samples of the image,

Median littering was also simulated but will not he discussed here due to its inferior performance compared to mean filtering .



20

F. A. Marvasti et al. Signal Processing 36 (1994) 13-30

Fig. 7. Image recovery using the nonlinear method represented in Fig . 5 : (a) the original 256 x 256 image : (b) the degraded image after 20% pixel loss - smoothed by a low-pass filter, S/N = 8 .9dB ; (c) linear interpolation (mean filtering) . S/N = 25 .5dB ; (d) recovered image using the division technique, S/N = 27 .9 dB .

a random array is generated from IMSL library . After quantizing the random samples, we get a 256x256 array of, for example, about 80% l's and 20% 0's . This array shows a 20% pixel loss and determines the nonuniform positions fp (x, y), and

therefore, by low-pass filtering this array, one gets the Jacobian rJ1, Eq . (10). After multiplying this array (ff (x, v)) by the image, we get the nonuniform samples -f,(x, y) - of the image at the Nyquist rate, Eq . (11). Low-pass filtering the nonuniform samples

The Division Method Fig . 8 . Color image recovery using the nonlinear method represented in Fig . 5: (a) the degraded image after 10% pixel loss ; (b) the image in (a) smoothed by a low-pass filter ; (c) recovered image using the division technique ; (d) the degraded image after 20% pixel loss ; (e) the image in (d) smoothed by a low-pass filter : (f) recovered image using the division technique .

22

F.A . Maruasri ei at. / Signal Processing 36 (1994) 13-30

. ................. .. (b)

Fig . 9 . Image recovery using the nonlinear method at 50% pixel loss: (a) recovery without a hard-limiter/low-pass filter : (b) recovery after a hard-limiter/low-pass filter .

we get Fig. 7(b), which has significant distortion . The division of the low-pass image represented in Fig . 7(b) by Jacobian I J I yields the image shown in Fig . 7(d), which reveals an impressive improvement of 18 dB over simple low-pass filtering . A 2-D linear interpolation (mean filtering) is also shown in Fig . 7(c) for comparison . In mean filtering, we average an array of 3 x 3 pixels to get an estimate of a lost pixel at the center of the array . 2 We observe that the nonlinear technique is about 2 .4 dB better than mean filtering . The nonlinear technique becomes much better than the mean filtering when the sampling rate is increased or, alternatively, if the pixel loss is lower than 20% . Besides black-and-white images, we also perform the simulations on color images for the nonlinear method . For color signals, we need to process each of the three components of color independently . Fig . 8 presents simulation results for the nonlinear method . Fig . 8(a) shows a degraded image after (11),

2 We

tried other array sizes but a 3 x 3 array seems to be the best choice.

10% pixel loss . After low-pass filtering this degraded image, we get Fig . 8(b) . Fig . 8(c) is the recovered image using the nonlinear technique . Significant improvement is observed relative to lowpass filtering . Figs. 8(d-f) are the repetition of the same experiment except at 20% sample loss . Although significant improvement is observed, some white spikes appear in the picture . These white spikes are due to the fact that at higher sample losses, the probability of loss of blocks of pixels increases . The output of the low-pass filter about these blocks is very close to zero. Therefore, a division by zero occurs . To alleviate this problem, we use another hard-limiter at the end of the divider (Fig . 5) followed by a low-pass filter . This additional process tends to smooth the undesirable spikes . Fig . 9 shows the effect of hard-limiting after the nonlinear recovery technique . Fig . 9(a) is the recovered image without hard-limiting at the end . The occurrence of spikes is widespread at 50% pixel loss . Fig . 9(b) is the hard-limited version followed by a low-pass filter . A noticeable improvement is observed . Objective evaluations in terms of S/N versus the percent of sample loss with and



23

F.A . Marcasli et al. / Signal Processing 36 (1994) 13 30

SNR (dB)

A a Nonlinear Method after a Herd limiter and a Low-Pass Filter

a a Is

s

a

so

Percentage of Lost Samples

Fig. 10. Performance of the nonlinear technique with and without a hard-limiter/low-pass filter.

without hard-limiting is shown in Fig . 10. This figure also shows that the proposed nonlinear technique is about 24 dB better than the degraded image at 20% pixel loss . By using a hard-limiter/lowpass filter to smooth out the spikes an additional 8 dB is gained .

4 .2 . The iterative method

For the iterative method, if the initial guess fo (x,y)is properly chosen, the speed of convergence can be increased . In our simulation, we choose either the low-pass version of the nonuniform samples or the mean filtered image as the initial estimate for fo (x,y). Obviously, the latter choice improves the speed of convergence . The simulation results for a 30% lost sample image are shown in Fig . 11 . Fig. 11(a) represents the image degraded after 30% pixel loss - smoothed by a low-pass filter . After two iterations, Fig. 11(b) shows an impressive improvement of about 27 dB (mean filtered esti-

mate is used for fo (x, y)) . Fig. 11(c) shows about 3 dB more improvement after five iterations . Over 4 more dB is gained after 15 iterations (Fig . 11(d)). When the low-pass of nonuniform samples is used as an initial estimate, the improvement in terms of S/ N for 1, 10 and 100 iterations is shown in Fig . 12; a comparison among the iterative technique, the nonlinear technique and mean filtering is also given in Fig . 13 . This figure shows that up to 20% sample loss, 10 iterations outperform the other two techniques . However, beyond 20% sample loss, the nonlinear technique has the best performance. For the same initial condition, Fig. 14 is a bar chart comparing the number of iterations needed to match the S/N of the nonlinear technique and mean filtering . This figure shows that as the percentage of sample losses increases, more iterations are needed to match the performance of the other two techniques . The loss of a large block of adjacent pixels may be fatal for the nonlinear method, however . Under these circumstances, mean filtering and the iterative methods are better alternatives .



24

F.A . Marvasli er at / Signal Processing 36 (1994) 13 30

Fig. 11 . Image recovery using the iterative method represented in Fig. 6 : (a) The degraded image after 30% pixel loss - smoothed by a low-pass filter, S/N = 7.1 dB; (b) recovered image after two iterations, $/N = 31 .87 dB ; (c) recovered image after five iterations, 5/N = 34 .75 dB ; (d) recovered image after 15 iterations, S/N = 38 .9 dB .

15 iterations . The results are surprisingly inThe iterative technique has been tested for distinguishable with the original image (Fig . 15(g)). the color image shown in Fig . 15. Figs . 15(a, c, e) As expected with more sample losses, more are, respectively, the degraded image after 10%, . 15(b,d,f) are the iterations are needed to get the same quality of 20% and 40% pixel loss. Figs image. corresponding recovered images after 5,6 and



25

F.A . Marvasti et al. ! Signal Processing 36 (1994, 13-30

SNR (dB)

A

51

\

10th Imradon

100th Iteration

a \Iteration

Is

s

Degraded Image

0

10

0

70

20

Percentage of Lost Samples

Fig . 12 . Performance of different iterations at different rates of pixel losses .

SNR (dB) 51

20

8

6

0 0

10

30

50

so

70

Percentage of Lost Samples

. Fig . 13 . Comparison of the nonlinear, iterative, and mean filtering methods



FA . Marnasti et aL / Signal Processing 36 (1994) 13-30

26

Number of Iterations needed to match SIN

A 04

Nonlinear Method

55 Mean Filtering

W

21 ai

23 t8

17

12 s 4

0

0

20

I N

m

so

w

70

'

Percentage of " Lost Samples

Fig. 14 . Number of iterations needed to match the performance of mean filtering and the nonlinear methods,

5. Conclusion In this paper we presented a method to describe nonuniform samples of a 2-D signal in the frequency domain . The spectral analysis suggests a simple nonlinear method to reconstruct the signal from its nonuniform samples . The same procedure can be generalized to n-dimensional signals . Since the nonlinear method requires a division, we have to make sure that the denominator is not zero . This requirement puts a constraint on the position of nonuniform points . The error decreases as the set of irregular samples gets closer to the set of uniform sampling set . In terms of pixel losses, significant improvement is observed with the decrease of the percentage of pixel losses . In comparison to mean filtering, the nonlinear technique is always better at random pixel losses. In case a block of adjacent pixels is lost, mean filtering may outperform the nonlinear technique . We also developed a theoretical background to recover images from lost samples using iterative techniques . This technique is based on ideal impul-

live nonuniform samples and the iteration is on the whole image as opposed to pixel by pixel . In comparison, the nonlinear method is easy to implement and performs quite well, particularly when the percentage of the lost samples is low . However, for a high-quality reconstructed image or a high percentage of lost pixels, the iterative method is preferred .

6 . Appendix A We have assumed in the proposed technique (Fig . 5) that I J I in (13) should not be equal to zero . In this appendix, we find a sufficient condition on the maximum values of 0(x a) and 0(y,„) in Eq . (13) . Since 0(x) and 0(y) are any functions that satisfy Eq . (3b), from [6] we know that it is always possible to interpolate band-limited signals 0(x) and ¢(y) (with bandwidths 1/(2T,) and 1/(2T2), respectively) such that 0(x,) = x„ - nT, and 0(y .) = ym - mT2 . According to the Bernstein inequality [12], the slope of a band-limited signal is related to the

The Iteration Method Fig. 15 . Color image recovery using the iterative method represented in Fig . 6 : (a) the degraded image after 10% pixel loss ; (b) recovered image after five iterations ; (c) the degraded image after 20% pixel loss ; (d) recovered image after six iterations ; (e) the degraded image after 40% pixel loss; (t) recovered image after 15 iterations .



28

F. A . Maruasii et al. / Signal Processing 36 (1994) 13-30

From inequalities similar to (A.l) and (A .2), we derive that a sufficient condition would be that the maximum of 0, (x, y) and 0 2 (x, y) be less than T,/(2n) and T 2 /(2n) . The derivation is as follows . According to the Bernstein inequality, we have

01 < 2ltO

0

lma ,(1/2T,)

and

< 2rt0, max (1/2T1 ), a0,

002 < 21t0 2ma ,(l/2T 2) and a02 y

< 2n02max (1/2T2 ) .

For the partial derivatives to be less than 1/2, it is sufficient that the maximum of 0, (x, y) and 0 2 (x, y) be less than T 1 /(2n) and T 2 /(2g) . For the special case shown in Fig . 2, the Jacobian is defined by Eq . (12) . For this special case, we can prove that if d0,/dx and 50 2/3y are less than 1, it is sufficient for the maximum of 0, (x, y) and 02(X, Y) to be less than T l /n and T2/x, respectively .

Fig. 15 . (g) The original signal .

bandwidth of the signal and the maximum amplitude of the signal, i .e ., dO(x) dx

dO(y) dy

7 . Appendix B < 2n0max (2T ~ .

(A .1)

(A .2)

< 2nomax (2 T2

If O ma, and 0 ma , are less than T 1 /n and T2 /n, according to (A .1) and (A.2) 101 and I¢I are less than 1 and hence J in (13) is always positive . Therefore, the nonuniform sampling deviations, O(x„) and o(ym ), should be chosen to be less than or equal to 0ma , and 0 .,, respectively . Note that the condition on the deviations 0(x,) < T l /n and O(y,,) < T2/n do not guarantee that 0 ma , and ,~' ma, are less than T,/7r and T2 /n, respectively . Likewise, we should be able to show that for the Jacobian defined in (7) to be always positive, it is sufficient for the maximum of 0,(x, y) and 0 2 (x, y) to be less than T,/(2n) and T2 /(2sc), respectively . The proof is as follows . If the absolute values of the partial derivatives (7) are less than 1/2, the Jacobian defined in (7) is always positive since

Here we prove that if the sampling set {x„ m ,Yn ,n } is a stable set, the iteration given in (23) will converge to a stable point, which is equal to the original image . ,/A+ i (x,

y) = APS .f (x, y) + (P - )PS)Jk(x, y),

(B .1)

where 2, f(x, y) and fk(x, y) are a convergence constant (relaxation parameter), the original finite energy signal and the kth iteration, respectively . P and S are, respectively, the 2-D band-limiting and ideal nonuniform sampling operators, and PSf(x, y) in (B .1) is the low-pass filtered version of the ideal nonuniform samples, which is known . In the following, we prove that there exists a range for A for which the process will converge to a stable point, which is equal to the desired signal f, i.e .,

km fk (x, y)

= f (X' y) .

PROOF . The iteration will converge: { lim fk 00,

J1=

(

1-Bx l

002

l-

I/ 00, \/ I/a0 2 \I

oy) - vt'y/vdx/

= I(1 - 1/2XI - 1/2) - 1/2 * 1/21 = 0 .

=

f(x,

y)} if P - 1LPS is a contraction . The operator P - APS is a contraction [10] if lf` k +"(x, v)II < IIf "'(x,y) 11

for all k,

(B .2)



F.A . Marra .ti et al ./ Signal Processing 36 (1994) 1330

where IIf (x,y)11 is the L 2 norm of f(x,y) and f (k) (x, y) - fk(x, y) - fk-1(x, y). From (B .1), we have r1 f

I P(f (kl(x, .))-)PS(f(k)(K,y))

II,

(x,Y)1, (B .3)

where 0 < r < 1 . So if (B .3) is satisfied, the theorem is proved . We now show that there is a i. such that (B .3) is true . The left-hand side of (B .3) can be written as

From (B .7) (B .9), we conclude that .[f P(f(k)(x y))PS(f(k)(x,),4)dxdy > AIIf ( "(x,y)'1 2 . (B .10) Therefore, (B .5) is proved by setting k, = A . To prove (B,6), we can write II PS(f "TV" O = J

II

PS( f IkI (X, I,)) II I P( f (1)(X . y)) - .; (kl (x,y) z +) 2 I.PS(f (k) (x,Y))1I 2 = f -

2) f

f P( flk)(x,y))PS(

f PS( f

i z

(kl(x

Y))ILU (k) (xnm, ynm)

x b(x - x .., y -

f (k l(x,y)dxdy

Now, we show that there are positive real such that

29

k,

. (B .4)

= L L

f f [PS (f

Ynm)

tk)(x,

dx dy

y))] ( f

tk) (xmn, Ynm)

and k 2 x 3(x -

f ff P(f (k)(x,y))PS(f(k)(x,y))dxdy > k1 lk) (x,y)11 2

xnm, y - ynm)

dx dy

[PS(f (k)(X, v))] (xnm,

L

ynm)( f

(k

(xnm, ynm)

(B .5)

and

E L

I PS( fl"(_)(, y)) 1 2

< k2

III II f

(x,y)

112

ynm)1

2

x[YI[PS(f (k) (x,Y))

To prove (B .5), we first note that

n

( ~ f O ( f lkl(x. v)) PS( f k)(x y)) dx dy = f f ( f (kI (x . y)) S( f 11) (x, v)) dx dy,

( .f (k) (xnm,

m (B .6)

(xn~ynm)1

2]

L2

J

m

< B If (k) (x,y) I , II PS(f ( k ) (x,

V) II .

(B .11)

(B .7) Hence, from (B .9), we have

since the operator P is self-adjoint (i .e ., f(Px)ydt=,x(Py)dt), P z =P, and Pxk =xk . Now, Eq . (B .7) can be written as x, y)) dx dy f f ( f ( "(x, y)) S ( f ( "((x,

I'PSlf k '(x,Y)) 11 5BIf (k) (x,y)I .

(B.12)

Therefore, hypothesis (B .6) (B .4)-(B .6) . we get

From

is proved .

I P( flk)(x,y))- ).PS( f (k) (x,4) 11 2 < r IIf tkl (x,Y) I 2 , if "k) (xnm,Ynm)] 2 .

(B .8) where r = I + ; 2 k 2

Since the points {xnm,Ynmt are assumed to be a stable sampling set, from the Nikoliskii inequality [1] we have, for all signals) - that are band-limited, ~ nFmf (( m , Ynm)

A <

I lf(x ,

y)

II

z

< B,

- 2Ak, _

In order to satisfy (B .3), r has to be in the range of 0 < r < 1, i .e ., given a particular k, and k z , k has to be in the following region of convergence for the iterative relationship given in (B .1) :

(B .9) 0<4<

where A and B are positive constants . depending only on the bandwidth and (xnm,Ynm).

k, <

kz .

kk1

(B .13)

30

F.A . Maruastiet a(. 7 Signed Processing 36 (1994) 13-30

8. Acknowledgments We would like to thank Mr. David Strong at the R&D Department of Donnelley Inc . for allowing us to use their facilities for our simulation results . 9 . References [I] P .L . Butzer and G . Hinsen, "Two-dimensional nonuniform sampling expansion - An iterative approach . I. Theory of two-dimensional band-limited signals", in : Applicable Analysis . Vol . 32, Gordon and Breach, London . 1989, pp. 53 67. [2] C. Cenkar, H . Feichtinger and M . Hermann, "Iterative algorithms in irregular sampling : A first comparison of methods", in : Conf. Proc. Phoenix Conference on Coinputers and Communications . Scottsdale. Arizona, March 1991 . [3] J .J . Clark, P.D. Lawrence and M .R . Palmer, "A Transformation method for the reconstruction of functions from nonuniformly spaced samples", IEEE Trans . Acoust . Speech Signal Process ., Vol . ASSP-33, No . 4, October 1985. [4] R . Cristescu and G . Marineseu, Applications of the Theory of Distributions, Wiley, New York, 1973 . [5] F . Marvasti. "Spectrum of nonuniform samples", Electronic Lett ., Vol . 20, October 1984. [6] F . Marvasti, A Unified Approach to Zero-Crossings and Nonuniform Sampling of Single and Multidimensional Signals and Systems, NONUNIFORM Publication, Oak Park . IL 60304. [7] F . Marvasti, "Reconstruction of 2-D signals from nonuniform samples or partial information", J . Opt . Soc. Amer., A, January 1989, pp. 52 55 . [8] F . Marvasti, A chapter contribution in : R . Marks II, ed., Advanced Topics in Shannon Sampling and Interpolation Theory, Springer . New York, 1993 . [9] F . Marvasti and T .J . Lee, "Analysis and recovery of sample-and hold and linearly interpolated signals with

irregular samples" . IEEE Trans . Signal Process ., August 1992 . [10] F . Marvasli, M . Analoui and M . Gamshadzahi, "Recovery of signals from nonuniform samples using iterative methods", IEEE Trans . Signal Process ., Vol . 39, No . 4. April 1991, pp . 872 -877 . [11] F . Marvasti, P . Clarkson, M. Dokic, U . Goenchanart and L . Chuande, "Recovery of speech signals with missing samples", IEEE Trans . Signal Process., December 1992 . [121 A . Papoulis, Signal Analysis, McGraw-Hill, New York, 1977 . [13] E . Plotkin . L . Roytman and M .N .S. Swamy. "Nonuniform sampling of band-limited modulated signals", Signal Processing, Vol . 4 . No- 4, July, 1982, pp . 295-303 . [14] E . Plotkin, L .M . Roytman and M .N .S . Swamy, 'Reconstruction of nonuniformly sampled band-limited signals andjitter error reduction" . Signal Processing, Vol . 7, No . 2, October 1984, pp . 151 160 . [l5] Y . Rahmat-Samii and R. Cheung, "Nonuniform sampling techniques for antenna applications", IEEE Trans . Antennas and Propagation, Vol . AP-35, No . 3, March 1987 . [16] H .E . Rowe, Signals and Noise in Communication Systems, Van Nostrand, Princeton, NJ, 1965 . [17] K .D . Sauer and J .P . Allebach, "Iterative reconstruction of band-limited images from nonuniformly spaced samples", IEEE Trans . Circuit and Systems, Vol . CAS-34, No . 12, December 1987 . [I8] S . Yeh and H . Stark, "Iterative and one-step reconstruction from nonuniform samples by convex projections", J . Opt . Sot . Amer ., A, Vol . 7 . No. 3, March 1990, pp . 491-499. [19] J .L. Yen . "On nonuniform sampling of band-limited signals", IRE Trans . Circuit Theory, Vol . CY-3, December 1956, pp . 251 -257 . [20] F. Marvasti, C . Liu and G. Adams, "Spectral analysis of irregular samples of multidimensional signals", IEICE Trans., Vol . E77-A . No . 2. February 1994 . [21] H .G . Feichtinger et al ., Chapter contribution, in : J . Bendetto, ed ., Wavelets, CRC Press, 1994.