Convergence analysis of an iterative method for the reconstruction of multi-band signals from their uniform and periodic nonuniform samples.

Abstract

One of the proposed methods for recovery of a band-limited signal from its samples, whether uniform or nonuniform, is the so-called Frame Method (1). In this method the original signal is reconstructed by iterative use of sampling-filtering blocks. Convergence of this method for linear invertible operators has been previously proved. In this paper we show that this method for non-invertible periodic nonuniform samplings as well as non-invertible uniform samples of bandpass (or multi-band) signals will lead to the pseudo-inverse solution. Convergence conditions in case of additive noise will also be discussed.

Key words and phrases: Iterative Method, Periodic Nonuniform Sampling, Pseudo-Inverse Solution, Bandpass and Multiband sampling

2000 AMS Mathematics Subject Classification--94A20, 65F10

1 Introduction

Although uniform sampling at the Nyquist rate for low-pass signals is quite straightforward, the extensions to bandpass and multi-band signals are not trivial. The bandwidth for such signals consists of separate positive frequency intervals with the overall length of B (thus, the Nyquist rate is 2B). Therefore, uniform sampling at the rate 2B is likely to lead to aliasing. The minimum alias-free uniform sampling rate for bandpass and multiband signals is not necessarily 2B [3].

One of the old proposed methods for reducing the sampling rate is the second order sampling [4]. In this method the classical uniform sampling points are substituted by two interlaced uniform sampling sets. The average sampling rate of 2B can be reached by proper selection of the interlacing parameter. Generalization to Nth order samplings (or periodic nonuniform sampling) for bandpass and multi-band signals have been studied by [5, 6, 7, 8]. Periodic nonuniform sampling is also generalized to sampling sets which are unions of shifted lattices ([9, 10]); these sets are not necessarily periodic. Implementations of periodic nonuniform samplings are usually fulfilled by passing the signal through different delaying filters and then uniformly sampling each filter output. The idea of using a general filter bank (instead of delay" filters) has been developed by [11, 12].

The proposed reconstruction methods are mainly based on interpolating functions. Since these functions are bandlimited, they cannot be timelimited. For practical implementations these functions should be truncated. In many cases these truncations result in considerable errors. The alternative methods are iterative approaches. In these methods, by repeated use of a simple but not a perfect reconstruction method, the output gradually converges to the perfect solution [2]. For optimizing the iterative method with respect to the convergence rate, accelerated methods has been proposed [13]. One of the advantages of the iterative methods is that when the sampling scheme is non-invertible, it. still converges while the interpolating methods diverge. We will show that the converging signal is the pseudo-inverse solution. In the next section, we first introduce our specific iterative method and then in Section 3 we will consider its convergence analysis for the case of Hermitian matrices; in this section we consider the effect of additive noise. In Section 4 we will check the results for a special non-invertible case of bandpass sampling where we show how the sampling-filtering block can be modeled with a Hermitian matrix. The extensions to periodic nonuniform sampling for bandpass and muliband signals is given in Section 5. Simulation results and conclusion are given in Sections 6 and 7 respectively.

2 The Iterative Method

One of the methods for reconstruction of a signal from its samples is the iterative method [14, 2]. In this method, by successive use of a crude reconstruction method, we increase the quality (measured by SNR) of the reconstructed signal and we may reach the original signal without error. The block diagram of this method is presented in Fig. 1. As it is shown, to recover the original signal z from its distorted version y, we repeatedly apply lowpass filtering and sampling denoted by the operator D. The mathematical formulation of this method is as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (1)

[FIGURE 1 OMITTED]

In the above equation [lambda] is called the relaxation parameter. Convergence of the iterative method is determined by this parameter. The range of [lambda] thta is required for convergence is usually a continuous interval which includes zero; thus, to avoid divergence, [lambda] is normally set near zero.

3 Convergence Analysis

We shall discuss convergence conditions for the iterative method in this section. The analysis is general for any linear operator D which could be modeled with a Hermitian matrix and is not limited to sampling problems. Let us assume that the input x is an l x 1 vector and the output y is

y = Dx, (2)

where D is an l x l Hermitian matrix. (We will show later that for the convergence of the iterative method, the matrix D should be either Nonnegative Real or Nonpositive Real.) We can also include a noise vector by assuming

y = Dx + n. (3)

The statement of the problem is that given the output y vector, we wish to find the input vector x, i.e., to find the inverse of the system. Since the noise vector or even the existence of the noise is not known, for the purpose of the reconstruction we must assume the distorting operation is defined by (2) while (3) is more likely to hold. Therefore, we implement the iterative method by successive multiplication of the vector by the D matrix. Thus, the reconstructed

vector in each iteration is found by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (4)

We define the error vector and the error matrix as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (5)

Thus, we can rewrite (4) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

[??] [z.sub.k] = - [E.sup.k] x + [lambda]([I.sub.lxl] + E + [E.sup.2] + ... + [E.sup.k-1])n (7)

[??] [x.sub.k] = - ([I.sub.lxl] - [E.sup.k] x + [lambda]([I.sub.lxl] + E + [E.sup.2] + ... + [E.sup.k-1])n. (8)

We are interested in the convergence and uniqueness of the sequence [x.sub.k] when k [right arrow] [infinity]. To find the answers of these questions, we need to evaluate powers of E as shown in (8).

Assuming D is nonnegative definite (Hermitian) with eigen-values [d.sub.1], [d.sub.2], ..., [d.sub.l], we have

0 = [d.sub.1] = [d.sub.2] = ... = [d.sub.c] < [d.sub.c+1] [less than or equal to] [d.sub.c+2] [less than or equal to] ... [less than or equal to] [d.sub.l], (9)

where c is the number of zero eigen-values (0 [less than or equal to] c [less than or equal to] n). Moreover, Hermitian matrices have orthonormal eigen-vectors; let [v.sub.1], [v.sub.2],..., [v.sub.l] be these orthonormal eigen-vectors

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

Therefore:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (11)

and

D = V diag ([d.sub.1], [d.sub.2], ..., [d.sub.l])[V.sup.H], (12)

where diag ([d.sub.1], [d.sub.2], ..., [d.sub.l]) represents a square diagonal matrix with its diagonal elements as [d.sub.l], [d.sub.2], ..., [d.sub.l]. Similarly we can write

E = [I.sub.lxl] - [lambda]D = V diag (1 - [lambda][d.sub.1], 1 - [lambda][d.sub.2], ..., 1 - [lambda][d.sub.l])[V.sup.H]. (13)

We will first study the case with no noise (2) and then we will consider the effect of additive noise as shown in (3).

3.1 The Case Without Additive Noise

Without additive noise, (8) becomes:

[x.sub.k] = ([I.sub.lxl] - [E.sup.k])x. (14)

Hence, the convergence of the algorithm is equivalent to the convergence of [lim.sub.k[right arrow][infinity]] [E.sup.k]. From (13) we have

[E.sup.k] = V diag[((1 - [lambda][d.sub.1]).sup.k], [(1 - [lambda][d.sub.2]).sup.k], ..., [(1 - [lambda][d.sub.l]).sup.k])[V.sup.H]. (15)

To assure convergence, it is required that [absolute value of 1 - [lambda][d.sub.i]] [less than or equal to] 1. Moreover, if for any, i we have 1 - [lambda][d.sub.i] = -1, although [E.sup.k] elements stay bounded, then (14) does not converge. In summary, to assure convergence, we should have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (16)

Since all eigen-values are positive, there exists a nonempty set of k values which satisfy each of these inequalities. If we had both positive and negative eigen-values, no solution for [lambda] would be achieved (equivalent to the previously stated condition of non-negativity or non-positivity of the real part). In the iterative method (4), [lambda] = 0 is a trivial case; thus, we have

convergence [??] 0 < [lambda] < 2/[d.sub.max]. (17)

Prom now on we assume that the above condition is fulfilled by proper choice of [lambda]. If D is an invertible matrix, then it has no zero eigen-values (c = 0 in (9) and, thus, for all i, [absolute value of 1 - [lambda][d.sub.i]] < 1. This means [E.sup.k] [right arrow] 0 as k [right arrow] [infinity] and from (14) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (18)

which means perfect reconstruction (and which was predictable due to the invertibility of D). Now let us assume that D is not invertible and has c zero eigen-values. As a consequence, E has c eigen-values equal to 1 and the rest have absolute values less than 1:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

Now we shall show that this result is also obtained by pseudo-inverse of D:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (23)

Comparing (21) and (23), we get

[x.sub.[infinity]] = [x.sup.+]. (24)

Thus, we prove the important result that, when the matrix D is not invertible, the iterative method converges to the pseudo-inverse solution.

3.2 The Case of Additive Noise

Assuming the additive noise has zero mean and D is invertible, we can write (3) as

y = Dx + n = D(x + [D.sup.-1] n). (25)

Thus, by the same argument as given for the noiseless case, we will have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)

Although the iterative method converges to the above equation, yet D-1 can amplify the noise component, and, hence, the final SNR may not be acceptable.

Now we assume that D is not invertible and has c zero eigen-values and, therefore, (19) holds for powers of E. If we define [e.sub.i] = 1 - [lambda][d.sub.i] (1 [less than or equal to] c [less than or equal to] l), then we will have [e.sub.1] = ... [e.sub.c] = 1 and [absolute value of [e.sub.i]] < 1 for c + 1 [less than or equal to] i [less than or equal to] l. Rewriting (19), we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (27)

Combining (8) with the above result, we get

[??] [x.sub.k] = ([I.sub.lxl] - [E.sup.k])x + [lambda](k [c.summation over (i=1)][v.sub.i][v.sub..sup.H.sub.i] + [l.summation over (i=c+1)] 1 - [e.sup.k.sub.i]/[1 - [e.sub.i][v.sub.i][v.sup.H.sub.i])n. (28)

When [lambda] is properly chosen, [E.sup.k] converges as k [infinity]. Moreover, since [absolute value of [e.sub.i]] < 1 for c + [less than or equal to] i [less than or equal to] l, we know that [e.sup.k.sub.i] [right arrow] 0 as k increases. Hence, the only part in (28) which may diverge is k([[summation].sup.c.sub.i=1][v.sub.i][v.sup.H.sub.i])n:

Convergence [??] ([c.summation over (i=1)] [v.sub.i][v.sup.H.sub.i])n = 0 (29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (30)

On the other hand, we know that {[[v.sub..sub.1]], ..., [v.sub.l]} form an orthonormal basis for our vector space; thus, n [perpendicular to] [v.sub.1], ..., [v.sub.c] implies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

where [??] is an l x 1 vector. Thus, we have shown that the convergence condition requires that n [member of] Range(D). Now if this condition is satisfied, we have convergence of the iterative method:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

We summarize the results in the following theorem:

Theorem 1 The iterative method described in (4) (for a Hermitian distorting matrix D) converges if and only if

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

where [??] is an l x 1 vector and we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

where [d.sub.max] is the maximum eigen-value of the matrix D and [D.sup.+] is its pseudoinverse. For invertible D we have [D.sup.+] = [D.sup.-1] and [??] = [D.sup.-1]n.

4 Iterative Reconstruction from Uniform Samples of Bandpass Signals

In this section we study a special case that is often non-invertible. Let x(t) be a bandpass signal with its Fourier transform X(f) (Fig. 2). If we simply sample this signal uniformly with a rate close to the Nyquist rate, i.e., 2B, the occurrence of aliasing effect for some frequency components is highly probable. Let [f.sub.1], ..., [f.sub.N] be equally spaced frequencies which cover the interval [[f.sub.L], [f.sub.H]. We define the original vector as

[x.sup.(f)] = [[X(-[f.sub.N]), ..., X(-[f.sub.1]), X([f.sub.1]), ..., X([f.sub.N])].sup.T]. (35)

After uniform sampling of x(t), followed by a bandpass filtering, some frequency components are distorted due to the aliasing. We assume the interference to be as shown in Fig. 3. Thus, the k higher frequency components will interfere with the k lower frequency components and the 2N - k middle components remain unchanged; hence, the distortion matrix (shown in (2)) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (36)

It is evident that this matrix, only for k = 0 is invertible and for k > 0, is non-invertible (since the sampling scheme was irreversible). To find the pseudoinverse of the matrix, we must evaluate the eigen-values:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (37)

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (38)

Eigen-values of D can be translated to those of E as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (40)

where [v.sub.1], ..., [v.sub.k] are eigen-vectors of D related to its zero eigen-values. It is easy to check that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (41)

Now that we have the eigenvectors, we can evaluate the final vector produced by the iterative method. For this purpose, we must first calculate the [E.sup.[infinity]] matrix:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (42)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (44)

In other words, using the iterative method, we have only divided the aliased frequencies by 2. Although this result could have been achieved without the use of the iterative method, the strength of this method becomes more evident in more complex cases such as multiband signals. In fact this method, without any knowledge about the aliasing parts, accomplishes the desired results.

5 Periodic Nonuniform Sampling for Bandpass and Multi-band Signals

To avoid high sampling rates for perfect reconstruction of bandpass and multiband signals, Periodic Nonuniform Samplings (PNS) have been suggested. In the Nth order PNS the sampling points {[t.sub.i]} are found by

[t.sub.aN+b] = aT + [k.sub.b] a [member of] Z, b [member of] {0,1, ..., N - 1} 0 [less than or equal to] [k.sub.1] < [k.sub.2] < ... < [k.sub.N-1] < T, (45)

where T is the sampling period and the average sampling rate is N/T. We denote the signal and its Fourier transform by x(t) and X(f), respectively. If we apply the Nth order PNS to the signal x(t), we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (46)

As shown in the above equation, each frequency component of the sampled signal may be affected by more than one component of the original signal. Since the original signal is bandlimited, each frequency component of the sampled signal is the combination of finite components. Let [f.sub.0], [f.sub.1], ..., [f.sub.m] be the initial frequencies which affect the frequency [f.sub.0] of the sampled signal. Thus:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (47)

We can represent the above summations as a matrix equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (48)

Due to the definition of the L function, we know

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (49)

where * denotes the conjugate operation. The above equation shows the fact that the matrix L is Hermitian ([L.sup.H] = L).

Up to now we have only considered m + 1 single frequency components of the continuous Fourier transform. Let us assume that we have sampled the signal X (f) at the sampling points of f [member of] {..., -2/rT, -1/rT, 0, 1/rT, 2/rT, ...}, where r is an arbitrary positive integer. Since the frequency step between these sampling points is an integer multiple of 1/(rT), the frequency components which play a role in aliasing are either all included or all excluded. Moreover, we will discard the out of band components and only focus on the components inside the original band (by means of filtering). The original signal was bandlimited and we have only considered the inband frequency components; thus, we are dealing with finite number of samples and we can consequently define the original and sampled Fourier transform vectors similar to (48). Similar to the above approach, we can show that these vectors are related to each other within a Hermitian matrix. The vectors are estimates of the continuous Fourier transforms and when the parameter r in the frequency sampling part increases, these estimates represent the continuous functions more precisely. We assume that we have chosen a large enough r so that our estimates are suitable representations. Now we are dealing with discrete signals which could be treated via vectors and, as shown here, the sampling-filtering block could be modeled with a Hermitian matrix (distortion matrix).

6 Simulation results

We have verified our theoretical results with MATLAB simulations. Fig. 4 shows a bandpass signals and its reconstructed versions from uniform samples at the Nyquist rate (which are non-invertible) for different iteration numbers. It is evident that for aliased components we have a division by 2 in the frequency domain (compare the sampled signal and the one at the 20th iteration). Fig. 5 shows that after 20 iterations (for [lambda] = 0.9) we have almost reached the final signal and the SNR curve is saturated. The case of noisy samples, which we expect to diverge, is shown in Fig. 6. Although the initial SNR of noisy samples is 25dB, divergence is clear by the reconstructed signal at the 20th iteration. In each iteration, we improve the signal while amplifying the noise. After a finite number of iterations, the pseudo-inverse solution of the signal is reached (neglecting the noise part) but noise amplification is continued; thus, we would expect to have an optimum recovered signal after a finite number of iterations which is given in the 7th iteration with [lambda] = 0.9 (Fig. 7). To include a more general case of PNS, a multiband signal is sampled at times [{n x 50ms, n x 50ms+ 37.5ms, n x 50ms + 41.7ms}.sub.n[epsilon]Z] (Fig. 8). Since SNR curve is saturated before the 20th iteration (Fig. 9), the recovered signal at the 20th iteration (Fig. 8) is the pseudo-inverse solution.

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

7 Conclusion

Uniform sampling at the Nyquist rate for bandpass (or multi-band) signals may not be sufficient to recover the signal exactly. On the other hand, various periodic nonuniform sampling schemes can be used at the Nyquist rate for perfect recovery. In either case, for discrete signals, sampling (uniform or periodic nonuniform) and filtering can be modeled as a Hermitian matrix. For the first case the matrix is usually non-invertible; for the second case the matrix may be invertible. Iterative methods can be used to find the inverse of a Hermitian matrix. Convergence analysis of the iterative methods show that for invertible matrices the iterative method converges to the original signal and, when the matrix is non-invertible, the iterative result is equivalent to the pseudo-inverse of the matrix. For the case of quantized samples or additive noise in the channel, iterative methods converge to the pseudo inverse if a finite number of iterations are used. On the other hand, if infinite number of iterations are used, then the recovered signal converges to the inverse of the system; in this case, depending on the noise structure, the inverse system may amplify or attenuate the noise component.

References

[1] H. G. Feichtinger and K. GrSechenig, Theory and practice of irregular sampling, CRC Press, 1994.

[2] F. Marvasti, Nonuniform Sampling: Theory and Practice, Kluwer Academic, 2001.

[3] N. L. Scott R. G. Vaughan and D. R. White, The theory of bandpass sampling, IEEE Trans, on Signal Process., 39(9), 1973-1984, 1991.

[4] A. Kohlenberg, Exact interpolation of band-limited functions, J.

Appl. Phys., 24(12), 1432-1436, 1953.

[5] T. Kida and T. Kuroda, Relations between the possibility of restoration of bandpass-type band-limited waves by interpolation and arrangement of sampling points, Electron. Commun. Japan, 69(1), 1-10, 1986.

[6] S. C. Scoular and W. J. Fitzgerald, Periodic nonuniform sampling of multiband signals, Signal Process., 28(2), 195-200, 1992.

[7] A. J. Coulson, A Generalization of Nonuniform Bandpass Sampling, IEEE Trans. on Signal Process., 43(3), 694-704, 1995.

[8] Y. P. Lin and P. P. Vaidyanathan, Periodically Nonuniform Sampling of Bandpass Signals, IEEE Trans. Circuits Sys., 45(3), 340-351, 1998.

[9] Hamid Behmard and Adel Faridani, Sampling of bandlimited functions on unions of shifted lattices, Journal of Fourier Analysis and Applications, 8, 2002.

[10] D. Walnut, Nonperiodic sampling of bandlimited functions on unions of rectangular lattices, 1996.

[11] C. Herley and P. W. Wong, Minimum Rate Sampling and Reconstrcution of Signals with Arbitrary Frequency Support, IEEE Trans. on Inf. Theory, 45(5), 1555-1564, 1999.

[12] R. Venkataramani and Y. Bresler, Optimal Sub-Nyquist Nonuniform Sampling and Reconstrcution for multiband signals, IEEE Trans. on Signal Process., 49(10), 2301-2313, 2001.

[13] K. Grochenig, Acceleration of the frame algorithm, IEEE Trans. on Signal Process., 41(12), 3331-3340, 1993.

[14] F. Marvasti, M. Analoui, and M. Gamshadzahi, Recovery of signals from nonuniform samples using iterative methods, IEEE Trans. Signal Proc., 45(3), 340-351, 1991.

Arash Amini

Advanced Communication Research Institute (ACRI) and EE Department

Sharif University of Technology, Tehran, Iran

[e-mail]

Farokh Marvasti

Advanced Conmmnication Research Institute (ACRI) and EE Department

Sharif University of Technology, Tehran, Iran

[e-mail]

(1) It is also called Marvasti's method by Feichtinger and GrSchenig [1] and Grochenig and Strohmer in the sixth chapter of [2].