Regularized Sampling of Multiband Signals
Abstract
This paper presents a regularized sampling method for multiband signals, that makes it possible to approach the Landau limit, while keeping the sensitivity to noise at a low level. The method is based on bandlimited windowing, followed by trigonometric approximation in consecutive time intervals. The key point is that the trigonometric approximation “inherits” the multiband property, that is, its coefficients are formed by bursts of nonzero elements corresponding to the multiband components. It is shown that this method can be well combined with the recently proposed synchronous multirate sampling (SMRS) scheme, given that the resulting linear system is sparse and formed by ones and zeroes. The proposed method allows one to trade sampling efficiency for noise sensitivity, and is specially well suited for bounded signals with unbounded energy like those in communications, navigation, audio systems, etc. Besides, it is also applicable to finite energy signals and periodic bandlimited signals (trigonometric polynomials). The paper includes a subspace method for blindly estimating the support of the multiband signal as well as its components, and the results are validated through several numerical examples.
I Introduction
Sampling is the operation that makes the discrete processing of continuous signals possible. The basic tool for this operation is Shannon’s Sampling Theorem, which states that a signal can be reconstructed without error from its samples taken at a rate equal to twice its maximum frequency (Nyquist rate), [1]. In some situations, however, the signal is multiband, in the sense that its spectral support is composed of a finite set of disjoint intervals. In this case, sampling at the Nyquist rate can be very inefficient. The design of alternative sampling schemes for this kind of signals has been an active research topic for decades. It was first demonstrated by H. J. Landau in [2] that if the signal is sampled nonuniformly, then the average sampling rate is lower bounded by the measure of its spectral support, (Landau lower bound). Later, it was shown in [3] that there exist nonuniform sampling patterns that approach this bound as much as desired. The sampling scheme in this reference was the socalled “multicoset sampling”, which consists in selecting a sampling rate above the Nyquist rate, and then choosing a periodic nonuniform subsequence of it. The subsequent investigations on this topic have mainly focused on devising low complexity and stable (low sensitivity) implementations of multicoset sampling, [4, 5, 6, 7]. Recently, compressed sensing techniques have been applied to multicoset sampling in [8, 9], so as to detect the band positions, as well as to reconstruct each of the components. An alternative method to discretize multiband signals is the random demodulator in [10, 11], in which the complexity of the analog processing is increased, with the aim of employing only analogtodigital converters of moderate input bandwidth. Finally, it is worth mentioning the synchronous multirate sampling (SMRS) scheme in [12]. Here, the multiband signal is sampled in several uniform grids with different rates, in order to produce a bunch of lowpass signals, in which the spectrum of the multiband signal appears aliased. The key idea is that this aliasing may not affect all frequencies of the original spectral support for all sampling rates. This fact is exploited in [12] by properly selecting the number of sampling frequencies, so as to achieve the reconstruction of the multiband signal. This technique has been further developed in a second paper, which also includes the application of compressed sensing [13].
An assumption that is common to these approaches is that of finite energy: the multiband components are assumed to be in the Lebesgue space , and the machinery of Fourier analysis for this space is applied extensively. This way to proceed simplifies the analysis, because signals in have a proper spectrum function, and the wellknown properties of the Fourier transform apply. Nevertheless, there is a mismatch between this model and the situation encountered in multitude applications, (communications, audio, navigation, etc). Quite often, it is not possible to sample a time interval that contains most of the signal’s energy due to its long duration. This limitation has negative effects, that can be readily noticed in the singleband case by analyzing the Sampling Theorem. Given a finite energy signal whose spectrum lies in , the Sampling Theorem provides the representation
(1) 
where is an integer and is any time shift following . Yet in practice, the infinite set of samples is not available, and then (1) must be truncated at some index
(2) 
This last formula is well known for its poor accuracy, [14]. The problem lies in that (1) converge only because the samples converge to zero as , due to the fact that has finite energy. So, if the energy of is mostly outside the sampling interval , the accuracy is poor due to the sinc tails corresponding to the neglected samples.
This problem is well known in the signal processing field, and it can be solved by regularizing the interpolation formula in (2), using various filter design techniques, [15]. A typical method consists in assuming that there is some excess sampling bandwidth, i.e, that the spectrum of lies in with (and not ), and then multiplying the sinc samples in (2) by a set of weights ,
(3) 
This regularization converts the illbehaved interpolator in (2) into a wellbehaved one, and this can be readily seen in the error trends: while the interpolation error of (2) decreases only as , [14], there are sets of coefficients for which the error of (3) decreases exponentially as , [20]. The price paid for the regularization is that it is not possible to represent signals using the full sampling bandwidth anymore, and it is also not possible to implement filters with sharp transitions, since this would produce convergence problems similar to those of the sinc series.
By analogy with the single band case, the multiband problem is currently in the stage equivalent to that of the sinc series in (1) and its truncated version in (2); i.e, some sampling schemes have been identified that approach the Landau rate, and allow the recovery of the multiband signal and its components, but there is no regularization procedure available. This lack of regularization can be noticed in analytical devices, like the slicing of the spectrum using sharp transitions, and in the use of discrete sequences with full spectrum, like coset sequences.
The purpose of this paper is to present a regularized sampling scheme for multiband signals. Relative to the existing methods in the literature, the one in this paper provides a flexible interface between the multiband signal and its discrete representation, that allows one to trade sampling efficiency for noise sensitivity. In the scheme, the reconstruction procedure is specified in the form of weighted trigonometric polynomials, which are able to interpolate the multiband signal and its components in consecutive time intervals.
The starting point in the derivation of the scheme is the boundness assumption which substitutes the usual finite energy assumption, i.e, the multiband components are assumed to have bounded time amplitude, but their energy can be unlimited. To make no assumption about the energy is meaningful since, for long signals, it usually lies outside the sampling interval and can be arbitrarily large. And to work assuming bounded components reflects clear physical and engineering constraints: signals have bounded peak power due to restrictions like the maximum transmitted power, or due to systems like automatic gain controls. With the new assumption, the multiband components have a lowpass equivalent whose real and imaginary parts belong to the Bernstein space , [16, chapter 6], i. e, they can be described as entire functions of exponential type that are real and bounded on the real axis. This new description has the drawback that the basic tool in spectral analysis, namely the Fourier transform (Plancherel theorem), is not directly usable. So, basic analytic procedures in the literature on multiband sampling [3, 12, 13, 17, 10, 11], like the division of the spectrum in disjoint intervals, the use of discrete sequences with full spectrum (like coset sequences), and the use of lowpass filters with sharp transitions are not valid on bounded bandlimited signals. However, these tools can be substituted by a simple yet powerful analytical device, as shown in this paper. If is the multiband signal, the device consists in multiplying by a bandlimited window , that approximately “selects” a given finite interval, in the sense that its timedomain content is concentrated in it. The key result is that there exist windows , for which the product is a trigonometric polynomial with negligible error in the interval selected by . Besides, the trigonometric polynomial “inherits” the multiband structure of , i.e, its coefficients are formed by bursts of values corresponding to the bands of . Sec. III is dedicated to this analytical device.
After this, the design of a sampling scheme is addressed in Sec. IV, where it is designed for sparse trigonometric polynomials, since the product can be viewed as a polynomial of this type with negligible error. In this section, it is shown that a finite version of the SMRS scheme produces a sparse linear system, in which the sensitivity to noise can be reduced by increasing the number of samples. Afterward, Sec. V shows how the finite sampling grid for can exactly match an infinite SMRS scheme for . This makes it possible to interpolate and its components at any . The interpolation error is then analyzed in Sec. VI. Finally, the blind sampling problem is studied in Sec. VII, where the MUSIC algorithm from direction of arrival estimation is adopted.
Since the notation in the paper is extensive, the basic symbols and operators are described in the next subsection, and there is a list of symbols at the end of the paper in Table I. The next section sketches the sampling method, in order to give a broad view of the concepts involved. The novel aspects of the paper are, fundamentally, the contents of Secs. III to V and Sec. VII.
Ia Notation
Several conventions have been adopted in the paper to simplify the notation:

Definitions are performed using the operator “”.

denotes the interval
(4) for arbitrary and .

The symbols , , and denote sets of distinct indices (integers). In and the subscripts ’’ and ’’ remind one of the meaning of these sets: thus, and are the index sets of and , respectively, after windowing using .

denotes the number of elements of a set .

The operator “” denotes periodization with period , that is, for a signal , it is
(5) 
Vectors and matrices are written in bold font and in lower and upper case, respectively, (, ).

denotes the identity matrix.

denotes the element of matrix ; denotes the th element of vector ; and the th column of a matrix .
Ii Sketch of the proposed method
Consider a multiband signal formed by bounded bandlimited components ,
(6) 
where the band of is , and these bands are known and appear in increasing order, , . The th component can be viewed as a signal whose baseband version
(7) 
has real and imaginary parts in the space . The spectral support of is the set
(8) 
The multiband sampling problem can be posed in a finite time interval, centered at an arbitrary instant , , , in which a set of samples is taken at distinct instants , but assuming that the components must be interpolated only inside a shorter interval , . The objective is then to find a set of functions , such that the components and can be interpolated using
(9) 
and
(10) 
for in . If there is a method to solve this last problem, then it can be applied repeatedly in overlapping intervals with integer , so as to interpolate and the components at any . So, for example, the value of any of the components in a regular grid of instants , with integer and arbitrary , could be obtained by evaluating (9) in the intersection of this grid with , and repeating the procedure for consecutive values of .
In this setting, consider a bandlimited window with spectrum lying in , , that is approximately time limited to the interval . Besides, assume that it is bounded between one and a minimum value greater than zero in , that is, in , . Note that is not the typical window in filter design, that sharply selects a given time interval and nulls the rest of the time axis. Actually, it is the dual case: the spectrum of sharply selects the band , and is very small but not strictly zero outside . If is smaller than the minimum separation among bands of ,
(11) 
then the product is also a multiband signal with components ,
(12) 
and spectral support
(13) 
Thus, is another multiband signal in which the bands of have been expanded, but are still disjoint due to (11). Now, the multiband sampling problem can be posed in (12) with sampling instants and sample values . If there is a satisfactory solution for this last problem, then the original multiband signal and its components can be readily interpolated in , simply by dividing by .
In the windowed model in (12), the signal is approximately timelimited. This allows one to apply the Sampling Theorem, but with the time and frequency domains switched, i.e, is repeated periodically with period in the time domain, and this operation samples the frequency domain. Here, the window acts as a timedomain “lowpass filter” with nonuniform response. The result of this operation is a trigonometric polynomial, with coefficients equal to samples of the spectrum of , taken at the frequencies with integer , that lie in (13). If denotes the integer indices of the frequencies inside the band of ,
(14) 
and denotes the union of these sets,
(15) 
then this approximate sampling theorem says that is the polynomial
(16) 
for in , where the are unknown samples of the spectrum of . Besides, the components can be interpolated using
(17) 
for in .
As will be shown in this paper, the accuracy of this approximation depends on the window . Actually, it is shown in the next section that the interpolation error is bounded by , where is a bound on the amplitude of , and a bound on the summation
(18) 
The fact is that there exist windows, as the one presented in this paper, for which grows with , and decreases exponentially as with trend roughly equal to . Thus, there is a moderate value of for fixed such that (17) can be regarded as an equality, for any fixed numerical precision. There are several methods to construct this kind of window, [18, 19, 20]. The one proposed in this paper in Sec. IIIA is based on the properties of the Fourier transform of the KaiserBessel function. Fig. 1 shows this kind of interpolation for a baseband BPSK signal using the window in this paper. The bold curve is the windowed signal, that is selected for in (18). Fig. 2 shows the interpolation error, which is below 200 dB roughly in . (For details, see Sec. VIIIA.)
Coming back to the formula in (16), if is sufficiently small, then can be regarded as a trigonometric polynomial. This implies that the windowed problem in Eq. (12) can be cast in a generic setting, as that of computing the coefficients of a sparse trigonometric polynomial from its samples at the abscissas . If is an index set, then this kind of polynomial can be expressed as
(19) 
The polynomial on the right in (16) has this form with and . In this generic setting, the sampling problem reduces to solving the linear system
(20) 
If this system has full column rank, then the coefficients can be obtained using
(21) 
where denotes the pseudoinverse of (20). Then, the solution of the sampling problem for is obtained by substituting this last equation into (19). This yields
(22) 
where
(23) 
Besides, since any “component” of is specified by a subset , , its associated sampling formula is obtained by replacing with in (22),
(24) 
Finally, the solution for a sparse polynomial in (22) can be applied to the windowed model in (12), if is identified with each of the sets , and with the full set of spectral samples . So, if (22) is used with sample values and index set , and then the effect of the window in (12) is removed by dividing by , the result is an interpolation formula for
(25) 
for in . And the same can be done for interpolating , but with index set ,
(26) 
The last two formulas solve the initial problem in (9), with functions given by
(27) 
This is in broad terms the sampling method proposed in this paper for multiband signals, and the next three sections are dedicated to clarify its fundamental aspects. The windowing and spectral sampling are explained in the next section. The key points in it are the approximate truncation of the multiband signal using , and the dual sampling theorem, which is demonstrated by means of the Poisson’s summation formula. Also, a specific window is selected in Subsection IIIA. Then, Sec. IV studies the selection of the instants so that the linear system in (20) has low sensitivity to perturbations. It turns out that a finite SMRS scheme produces a sparse linear system of ones and zeroes, in which it is possible to reduce the sensitivity by slightly oversampling. Next, Sec. VII shows that the finite SMRS scheme in Sec. IV can be perfectly integrated into an infinite SMRS scheme for the initial multiband signal. This means that the finite scheme in an interval with arbitrary employs samples lying in the infinite scheme exclusively.
Iii An approximate dual sampling theorem: approximation by means of trigonometric polynomials
Recall the multiband signal in (6) with components , and let denote specific bounds for them,
(28) 
Also assume that it is necessary to approximate the value of for some or all , , in an interval for a specific . To perform this approximation, take a bandlimited window in with spectral support contained in , that is bracketed in as for a specific , and that additionally fulfills the concentration property in Eq. (18) with . Next, form the new signal
(29) 
This product expands each of the bands of , so that the spectral support of is contained in the set (13). Besides, this windowing affects all components of equally, i.e, also consists of components ,
(30) 
where
(31) 
If is smaller than the minimum separation among bands [Eq. (11)], it turns out that is also a multiband signal (its bands are disjoint). The inequalities in (18) and (28) imply that and can be well approximated by their respective periodic versions, and , defined using (5). It is
(32) 
and
(33) 
Now, the signals and are bandlimited and periodic. So, by the Poisson’s summation formula [16, Sec. 2.3], they are sparse trigonometric polynomials, whose coefficients correspond to nonzero samples of the spectra of and at frequencies of the form , respectively. The sets of indices of these frequencies, denoted and , have already been defined in (14) and (15). Thus, and can be expressed as
(34) 
and
(35) 
where is the (unknown) value of the spectrum of at frequency . So, the disjoint spectral intervals of translate into disjoint bursts of coefficients in . Note that and have proper spectrum functions, because they are bandlimited signals that lie in . This is not true for either or .
Finally, the desired formulas are obtained from (34) and (32). For in , it is
(36) 
with error bounded by , and
(37) 
with error bounded by .
A specific window is presented in the next subsection.
Iiia Selection of a window function
As shown in several applications [20, 21, 22, 23], a bandlimited window with excellent concentration properties can be constructed from the Fourier transform of the KaiserBessel function. For the interpolation problem in this paper, the proposed window is
(38) 
where
(39) 
In (38), the second sinc function in the numerator is the Fourier transform of the KaiserBessel function, and the concentration of is roughly proportional to the denominator in (38), which increases exponentially as when . The first sinc function in (38) is necessary so as to damp the tails of , given that it would otherwise not belong to . A bound for the inequality in (18) is derived in Ap. A for this window, which is well approximated by
(40) 
Here, is the value of in (38) that achieves the corresponding in (40). Note the exponential rate in this last expression for . As a consequence, the model accuracy can be increased arbitrarily by increasing . For it is and for it is . The actual bound on the interpolation error is given by , where depends on . A value for can be obtained numerically using the definition
(41) 
The error bound versus is plotted in Fig. 3 for . In this figure and were selected using (40). Note that any practical accuracy can be achieved by increasing the product .
Iv Low sensitivity sampling of a sparse trigonometric polynomial
The generic sampling problem for sparse trigonometric polynomials, as posed in Sec. II, Eqs. (19) and (20), does not impose any constraint on the distribution of the instants . However, the condition number of (20) can vary wildly with it, as can be easily checked numerically. The distribution in which the sampling points are equally spaced and cover is interesting, because it exploits the sparseness of the polynomial, i.e, it produces equations in which most coefficients are zero.
To see this point, assume is uniformly sampled at
(42) 
for , , where is an arbitrary instant, and define the coefficients
(43) 
Next, compute the DFT of the sequence (divided by ), defining a corresponding set of coefficients :
(44) 
This derivation shows that can be computed either from the samples using its definition,
(45) 
or from the coefficients using
(46) 
So, if the samples are known, Eq. (45) allows one to compute the with the best possible conditioning, given that DFT matrices have condition number equal to one. And Eq. (46) shows that is just the sum of some of the coefficients (those with index congruent with modulo ). In some cases, the sum in Eq. (46) may contain a single element, and then one of the coefficients would be revealed by , i.e, with and . Usually, however there are several summands, but if the are known, then (46) is in any case a sparse linear system of ones and zeroes.
This argument suggests that a regular grid like (42) is a good choice, in order to obtain a well conditioned linear system. But the selection of is problematic. On the one hand, if is larger than or equal to the difference between the largest and the smallest element in , then the sum in (46) has a single summand, and the reveal the coefficients completely. But the oversampling would be too large in this case. On the other hand, if is smaller than that difference, the system in (46) could be underdetermined. The solution proposed in the sequel consists in adopting a finite version of the SMRS scheme in [13], in which several integer moduli are used instead of a single one.
Let us state first the solution of the interpolation problem for this scheme, and then study its sensitivity. For a fixed set of moduli , let denote a finite SMRS scheme with initial instant ,
(47) 
and let denote the scaled DFT for modulo as in (44),
(48) 
Then, Eq. (46) for all moduli can be written as
(49) 
, . This is a linear system with unknowns . If it has full column rank, then its solution can be written is terms of its pseudoinverse, denoted , i.e,
(50) 
Now, Eqs. (43) and (19) allow one to write in terms of the coefficients ,
(51) 
Finally, if Eqs. (48) and (50) are substituted into this equation, and the summations are then reordered, the result is
(52) 
where
(53) 
The last two equations specify the solution to the interpolation problem using the finite SMRS scheme. Note that (52) can be used to reconstruct any component of , simply by restricting the index set in , i.e, if , then
(54) 
Also, note that here is the specific reconstruction function in the SMRS scheme with initial instant , while in Sec. II, Eq. (22), is a reconstruction function for a generic sampling set and fixed .
The selection of the moduli can be done now so as to make this method viable, i.e, so that the linear system in (49) has full column rank, and the sensitivity to perturbations of (52) is low. Let us define first the sensitivity measure
(55) 
and its supremum in ,
(56) 
is the standard deviation factor due to any perturbations in the samples . So, if these samples are contaminated by perturbations of zero mean and variance , then the value of computed using (54) has variance .
The following method gives a suitable collection of moduli :

Select positive integers and following and then set , .

Check whether the linear system
(57) is overdetermined (or determined). If not increase either or .

Finally, the sensitivity can be reduced by adding other which are selected so as to minimize , or if the performance for a specific components must be improved.
V Interpolation of a bounded multiband signal using the SMRS scheme
As shown in Sec. III, Eq. (37), the sampling problem for a bounded multiband signal can be reduced to a problem of the same type, but involving a sparse trigonometric polynomial. Specifically, it was shown that there is a polynomial such that
(58) 
where lies in . This formula specifies how the multiband signal can be interpolated if the polynomial is known, and how the polynomial can be obtained from samples of .
The objective of this section is to show that the finite SMRS scheme in the previous section of the form
(59) 
can be integrated into the infinite SMRS scheme
(60) 
More precisely, there is a such that there is a onetoone correspondence between the samples of in the scheme (60) that lie in , and those of in the scheme (59). To demonstrate this, assume that it is necessary to sample at one of the instants in (59). Let denote the unique integer such that lies in ,
(61) 
From (58) and recalling that is periodic, it is
(62) 
The argument of in this formula belongs to the scheme in (60) for any pair of indices if . So, for this , it is
(63) 
This is the onetoone correspondence between the samples of in (59) and those of in the intersection of (60) with .
Now, it is possible to write an interpolation formula for . For this it is enough to “sample” using (63), then obtain the polynomial approximately using the formula in (54), and finally undoing the windowing [dividing by ]. To simplify the notation, define first the instants
(64) 
From (63), the samples of are
(65) 
If the sampling formula in (54) is applied to with these sample values, it is
(66) 
And, finally, (58) gives the desired formula for ,
(67) 
The formula for is the same but with replaced with .
Vi Error analysis
In the error analysis in the sequel, the sampling scheme is denoted using a single index as instead of two, so as to simplify the notation. This also affects the function in (53), that will be denoted . Also, the error terms will be written inside square braces for readability.
Let us recall the sampling and interpolation process on in the last three sections, but keeping track of the different errors, so as to produce an error bound. First, the multiband signal can be approximated using the polynomial , [Eq. (37)],
(68) 
Next, can be reconstructed using the sampling formula in (52) with , , and samples ,