1

Sampling of Pulse Streams: Achieving the Rate of Innovation Kfir Gedalyahu, Ronen Tur and Yonina C. Eldar Department of Electrical Engineering Technion—Israel Institute of Technology, Haifa 32000, Israel E-mail: {kfirge@techunix,ronentur@techunix,yonina@ee}.technion.ac.il.

Abstract—We consider the problem of sampling signals which are comprised of pulse streams. This model belongs to the recently introduced framework of signals with finite rate of innovation. The minimal sampling rate for such signals is the number of degrees of freedom per unit of time, referred to as the rate of innovation. Although sampling of pulse streams was treated in various works, either the rate of innovation was not achieved, or the pulse shape was limited to diracs. In this work we propose multichannel schemes for arbitrary pulse streams, operating at the rate of innovation. The proposed approach is based on modulating the input signal in each channel with a properly chosen waveform, followed by an integrator. We show that the pulses delays and amplitudes can be recovered from the samples using standard spectral estimation tools. The resulting scheme is flexible and exhibits better noise robustness than previous approaches. Index Terms—Finite rate of innovation, sub-Nyquist sampling, time delay estimation.

I. I NTRODUCTION One of the traditional assumptions in sampling theory suggests that in order to perfectly reconstruct a bandlimited analog signal from its samples, it must be sampled at the Nyquist rate, i.e, twice its highest frequency. However, other prior knowledge on the signal, rather than band limitation, can be exploited in order to reduce the rate [1]. A specific prior was explored by Vetterli et al. [2], who considered signals having a finite number of degrees of freedom per unit time, termed by the authors as the finite rate of innovation (FRI) property. It was suggested in [2] that the minimal sampling rate required for this setting, is the rate of innovation (ROI). An interesting model, which belongs to the FRI class, is of streams of pulses. In this model, the signal is defined by the time delays and amplitudes of some known pulse shape. Such signals are prevalent in applications such as bio-imaging, neuronal activity, radar and communications. In [2], a scheme for periodic streams of diracs was proposed, allowing perfect recovery from a number of samples equal to the degrees of freedom of the signal. However in practical applications, finite or infinite streams are typically encountered, rather than periodic streams. For the finite case, various single channel sampling methods were proposed, based on Gaussian kernels [2], or polynomials and exponentials reproducing kernels [3]. However, all these methods are instable for a large number of pulses per unit time. Another sampling scheme was presented in [4], exhibiting superior noise robustness over former approaches. Exploiting the compact support of the sampling kernels in [3], [4], both methods were extended to the infinite case, while the latter exhibits better noise robustness. Although the results in [4] reduced the sampling rate in comparison with previous works, it still cannot achieve the ROI. In [5] the authors proposed a multichannel sampling scheme enabling straightforward reconstruction of the unknown parameters. However, the underlying signal model is limited to a single pulse per sampling period, which is a restrictive requirement. Another multichannel approach comprised of a chain of integrators was proposed in [6]. As we show in simulations, this method exhibits instability in the presence of noise. Although both methods [5], [6] achieve the ROI, they are limited to streams of diracs only. To the

best of our knowledge, a stable and efficient sampling scheme for infinite pulse streams with arbitrary shape is still lacking. Our first contribution is for finite streams of pulses. We design a multichannel sampling scheme, where in each channel the signal is modulated by a parametric waveform and integrated over a finite time interval. We show that by proper selection of the waveform parameters, the signal can be perfectly recovered from a number of samples equal to its degrees of freedom. Extending our results to the infinite setting, we propose a scheme which enables perfect reconstruction of the signal, while sampling at the ROI. In addition, our method treats more general pulse shapes, and is more robust to noise comparing to previous work. We then consider a special case of the infinite model, proposed in [7], of pulse streams with shift-invariant (SI) structure. We show that this additional knowledge on the signal’s structure, can be exploited, in some cases, to further reduce the rate. The remainder of this paper is organized as follows. In Section II we derive a multichannel scheme for finite pulse streams. Section III extends our results to the infinite case, where we show that our approach achieves the ROI. In Section IV we discuss the relation of our results to previous works. Finally, in Section V, we present simulations demonstrating the noise robustness of our method. II. F INITE S TREAM OF P ULSES A. Problem Formulation We consider a finite stream of pulses, defined as x(t) =

L X

al h(t − tl ),

tl ∈ I ⊂ [0, T ), al ∈ C,

(1)

l=1

where h(t) is a known pulse shape, and {tl , al }L l=1 are the unknown delays and amplitudes. We assume that the number of pulses L is known, and the time-delays {tl }L l=1 are restricted to lie in a finite time interval I ⊂ [0, T ). We further assume that the following condition holds: h(t − tl ) = 0, ∀t ∈ / [0, T ) l = 1 . . . L, (2) i.e., the signal x(t) is confined to the time-window [0, T ). Since the signal is defined by 2L parameters, clearly at least 2L samples are required in order to represent the signal. Our goal is to design a sampling and reconstruction method which perfectly reconstructs the signal x(t) from a minimal number of samples. We first show that our problem can be related to another known signal processing problem, of model-based complex sinusoids (cisoids) parameter estimation [8]. Since x(t) is confined to the interval t ∈ [0, T ), it can be expressed using its Fourier series expansion X 2π x(t) = X[k]ej T kt , t ∈ [0, T ), (3) k∈Z

where X[k] are the corresponding Fourier series coefficients. It can be shown [9] that the coefficients X[k] satisfy µ ¶X L 2π 2π 1 k (4) al ej T ktl , X[k] = H T T l=1

2

1 T

e

RT 0

h

(·)dt

X −⌊ K2 ⌋

i

1 T

+j 2π T ⌊K/2⌋t

s1 (t) =

X

RT 0

(·)dt

c1

−j 2π kt T

s1k e

k∈K

x(t)

x(t)

1 T

RT 0

h

(·)dt

X +⌊ K2 ⌋

1 T

i

−j 2π T ⌊K/2⌋t

sp (t) =

e

X

RT 0

(·)dt

cp

kt −j 2π T

spk e

k∈K

Fig. 1.

Direct multichannel sampling scheme.

Fig. 2.

where H(ω) denotes the continuous-time Fourier transform (CTFT) of h(t). We by K a set of K consecutive indices for which ¡ denote ¢ H 2π k = 6 0, ∀k ∈ K. We assume that such a set exists throughout T the paper. Denote by H the K × K diagonal matrix with kth entry ¡ ¢ 1 H 2π k , and by V(t) the K × L matrix with klth element T T 2π e−j T ktl , where t = {t1 , . . . , tL } is the vector of the unknown delays. In addition denote by a the length-L vector whose lth element is al , and by x the length-K vector whose kth element is X[k]. We may now write (4) in matrix form as x = HV(t)a.

(5)

Since H is invertible by construction we may define y = H−1 x. Writing the kth element of the vector y explicitly we get yk =

L X

al e−j

2π kt l T

.

(6)

l=1

This is a standard problem of finding the frequencies and amplitudes of a sum of L cisoids. The recovery of the cisoids frequencies can be performed using the annihilating filter method [2], or various other techniques [8], as long as K ≥ 2L. Once the time-delays are known, the amplitudes may be obtained via a = V† (t)y, where V† (t) denotes the Moore-Penrose pseudo-inverse of V(t). B. Direct Multichannel Sampling Scheme Our goal now is to design a sampling scheme which allows to obtain the vector x from time-domain samples. For simplicity, we set K to be an odd number, and choose the set K = {−⌊K/2⌋, . . . , ⌊K/2⌋}. We note here that our results can be extended to any set K of consecutive indices. It can be easily seen that the Fourier coefficients, X[k], are directly obtained using the sampling scheme depicted in Fig. 1. The sampling process in each channel involves modulating the signal with a complex exponential, followed by an integrator operating over the interval [0, T ). C. Waveform Modulation Scheme We now extend the scheme presented in the previous subsection, to a more general one, by mixing several Fourier coefficients in each channel. As we will demonstrate, the additional degrees of freedom may be advantageous for practical implementation considerations. We consider p sampling channels. In the ith channel we modulate the signal using the waveform X 2π si (t) = sik e−j T kt , (7)

Multichannel scheme with modulating waveforms.

The resulting scheme is depicted in Fig. 2. We define the p × K matrix S with sik as its ikth element, and c as the p-length sample vector with ith element ci . We may now write (8) in matrix form as c = Sx.

(9)

As long as S has full column rank, where p ≥ K is a necessary condition, we can obtain x out of the samples by x = S† c. The direct sampling scheme presented earlier, is a special case of this more general scheme, with p = K and S = I, the identity matrix. We summarize our results in the following theorem. Theorem 1. Consider a finite stream of pulses given by (1) and assume that condition (2) is satisfied. Using the multichannel sampling scheme depicted in Fig. 2, the signal x(t) can be perfectly reconstructed from the samples {ci }pi=1 , as long as p ≥ K ≥ 2L and the matrix S in (9) is left invertible. We note here that in contrast to [2], [3], the proposed approach is stable in the presence of noise. In addition, as we discuss in Section IV-A, the scheme in [4] can be viewed as a special case of our method, and therefore this work generalizes the results of [4]. We now give two examples for useful modulating waveforms. 1) Cosine and Sine waveforms: Assume p = K. We take the first ⌊K/2⌋ waveforms as cos( 2π kt), the next ⌊K/2⌋ waveforms T are taken as sin( 2π kt) and the last one will be the constant function T 1. Clearly, each of these waveforms is of the required form, since each sine and cosine can be expressed as a sum of two complex exponentials. It can be easily verified, that this selection yields an invertible matrix S. The advantage of this choice, is in the fact that the waveforms are real, in contrast to complex exponentials required in the direct scheme. 2) Periodic Waveforms: The second example is based on the fact that every periodic function can be expressed using its Fourier series expansion. Passing such a periodic signal through a low pass filter (LPF), will reject its high frequency components. The resulting signal will have the required form of (7). Proper selection of periodic functions can yield a left invertible mixing matrix S. One example for such a selection, is taking the waveforms as periodic streams of rectangular pulses modulated by ±1 [9], [10]. For this setting, it can be shown that a single waveform is sufficient, where each channel uses a delayed version of the common waveform [9]. This reduces the design complexity, since a single waveform generator replaces the numerous oscillators required in the direct scheme of Fig. 1.

k∈K

where sik are weighting coefficients. The modulation operation is then followed by integration over the interval [0, T ). The sample at the ouput of the ith channel is given by X ci = sik X[k]. (8) k∈K

III. I NFINITE S TREAM OF P ULSES A. General Model We now consider an infinite stream of pulses defined by X x(t) = al h(t − tl ), tl ∈ R, al ∈ C. l∈Z

(10)

3

1Z (·)dt T Im s1 (t) =

X

2π

s1k e−j T

t = mT

The use of the results from [7], which treated a similar set of equations, gives us the following condition for a unique recovery of the delays from (13):

c1[m]

kt

k∈K

x(t) 1Z (·)dt T Im sp (t) =

X

2π

spk e−j T

cp[m]

Multichannel scheme for infinite stream of pulses.

We assume that there are no more than L pulses in any interval Im , [(m − 1)T, mT ] , m ∈ Z and that within each interval condition (2) holds. The maximal number of degrees of freedom per unit time, also known as the ROI [2], is 2L/T . We now present a multichannel sampling and reconstruction scheme which operates at this rate. Consider an extension of the sampling scheme presented in Section II-C, where we sample every T seconds. Upon each sample we reset the integrators, so that the mth sample results from an integral over the interval Im . The resulting sampling scheme is depicted in Fig. 3. Since the mth sample is influenced by the interval Im only, the infinite problem may be reduced into a sequence of finite streams of pulses which we have already solved in the previous section. We state our result in the following theorem. Theorem 2. Consider an infinite stream of pulses given by (10). Assume that there are no more than L pulses within any interval Im , [(m − 1)T, mT ] , m ∈ Z, and that condition (2) holds for all intervals Im , m ∈ Z. Consider the multichannel sampling scheme depicted in Fig. 3. Then, the signal x(t) can be perfectly reconstructed from the samples {ci [m]}pi=1 , m ∈ Z as long as p ≥ K ≥ 2L, and the matrix S in (9) is left invertible. It should be emphasized that Theorem 2 presents the first sampling scheme for general pulse shapes, operating at the ROI. Furthermore, as we show in simulations, our method is more stable than previous approaches. B. Stream of Pulses with Shift-Invariant Structure We now focus on a special case of the infinite model (10), where the signal has an additional SI structure. This structure is expressed by the fact that in each period T , the delays are constant relative to the beginning of the period. Such a signal can be described as x(t) =

L XX

aℓ [m]h(t − tℓ − mT ),

η = dim (span ({a[m], m ∈ Z}))

(15)

tℓ ∈ I ⊂ [0, T ),

(11)

where aℓ [n] denotes the ℓth pulse amplitude in the mth period. This model was first proposed in [7]; we discuss the relation to this work in Section IV-B. Our aim is to show that the special SI structure can be exploited in order to reduce the sampling rate, and to improve the noise robustness of the recovery algorithm. To this end, we define the vector c[m] as the length-p vector containing the samples of the mth period. We then define the K-length vectors y[m] as y[m] = (SH)−1 c[m],

denotes the dimension of the minimal subspace containing the vector set {a[m], m ∈ Z}. This condition implies that in some cases K and the number of sampling channels p, can be reduced beyond the lower limit 2L for the general model, depending on the value of η. The recovery of the delays from (13) can be accomplished using methods such as ESPRIT [11] or MUSIC [12], which use the combined information from all the vectors y[m]. When η = L, these algorithms can be applied directly on (13) to recover the delays, using p ≥ K ≥ L + 1 sampling channels. When η < L, an additional smoothing stage [13] is required prior to the use of these methods. Therefore, when the amplitudes of the pulses vary sufficiently from period to period, the common information about the delays can be utilized in order to reduce the sampling rate to (L + 1)/T . Moreover, the approach presented here can improve the delays estimation compared with the one used for the general model, since it uses the mutual information between periods, rather than recovering the delays for each period T separately. We can interpret this result by calculating the ROI of (11). In the time interval [−mT, mT ], there are L unknown delays, and 2Lm unknown amplitudes, therefore the ROI is given by 2mL + L L = . (16) 2mT T The sampling rate required by our method for this case, is close to that value and equals (L + 1)/T . We summarize our results in the following theorem. ρ = lim

m→∞

Theorem 3. Consider a stream of pulses with SI structure, given by (11), for which condition (2) holds. Using the multichannel sampling scheme depicted in Fig. 3, the signal x(t) can be perfectly reconstructed from the samples {ci [m]}pi=1 , m ∈ Z, as long as the matrix S in (9) is left invertible and ( L + 1 when η = L p≥K≥ 2L when η < L. IV. R ELATED S AMPLING S ETUPS A. Single-Channel Sampling with the SoS Filter

m∈Z ℓ=1

(12)

where H is given by (5). It can be easily shown that, under condition (2), these vectors satisfy y[m] = V(t)a[m],

(14)

where

kt

k∈K

Fig. 3.

K ≥ 2L − η + 1,

t = mT

(13)

where a[m] denotes the L-length vector whose ℓth element is given by aℓ [m].

The work in [4] considered single-channel sampling schemes of pulse streams. The proposed approach is based on a special filter, comprised of a sum of sinc (SoS) functions in the frequency domain. This filter can be expressed in the time domain as µ ¶X t g(t) = rect bk ej2πkt/T , (17) T k∈K

where {bk }k∈K are nonzero coefficients. Assuming a periodic input signal, by passing the signal through the SoS filter and sampling its output uniformly at a rate of T /p, we obtain the following samples X j 2π kn c[n] = bk X[k]e p , n = 0, . . . , p − 1. (18) k∈K

Using the parametric matrix V(t), defined in (5), with the vector ts = {0, T /p, . . . , (p − 1)T /p}, and defining the diagonal matrix B with kth diagonal element bk , (18) can be written in matrix form as c = V(−ts )Bx.

(19)

4

60 40 20

MSE [dB]

Therefore, this scheme can be viewed as a special case of our mixing multichannel sampling method, with a mixing matrix S = V(−ts )B. In that case, the p samples taken over one period in [4], are equal to the samples at the output of the p channels in our scheme. The sampling scheme in [4] was later extended to non periodic signals, by using an r-fold periodic continuation of the SoS filter, where r is a parameter depending on the support of the pulse h(t). However, this approach does not achieve the ROI for infinite streams of pulses.

0 −20 −40 −60 tones rectangular integrators

−80 −100 0

B. Multichannel Schemes for Shift-Invariant Pulse Streams The work in [7] treated the signal model (11) discussed in Section III-B, and proposed the scheme depicted in Fig. 4. In each sampling channel, the input signal is filtered by a band-limited sampling kernel s∗ℓ (−t) followed by a uniform sampler operating at a rate of 1/T . After sampling, a properly designed digital correction filter bank W−1 (eωT ) is applied on the sequences. It was then shown that the ESPRIT algorithm can be applied on the corrected samples, in order to recover the unknown delays. The sampling rate of this scheme is similar to the one presented here, i.e, it generally requires a sampling rate of 2L/T , where for certain signals it can be reduced to (L + 1)/T . s∗1 (−t) x (t)

d1 [n]

..



W−1 ejωT

s∗p (−t)

Fig. 4.

c1 [n] t = nT

cp [n] t = nT



.. dp [n]

Proposed sampling scheme in [7].

The advantage of the scheme in [7], is that it does not require condition (2), and therefore applies also to pulses with infinite time support. However, for pulses satisfying (2) our approach has several advantages. First, the equivalent stage for the digital correction in [7], is the appliciation of the matrix (SH)−1 on the samples. This correction stage can be viewed as a one-tap digital filter bank, in contrast to the filter W−1 (eωT ) which generally has a larger number of taps, and therefore the filtering operation has higher computational complexity. An additional benefit of our method, is that the technique of [7] requires collecting an infinite number of samples, even in the case where the input signal contains a finite number of periods. Moreover, consider an infinite signal which we are interested only in a certain finite time interval of it. The method in [7] does not allow processing of the interval of interest alone, since the samples are influenced by all the periods outside this interval as well. In contrast, our method can collect samples only from the relevant periods, since we integrate over finite time intervals. V. S IMULATIONS We now examine the performance of our method in the presence of noise, when working at the ROI. Since the only method which operates at the same rate, without additional constraints on the pulses locations, is the one presented in [6], we compare our results to the ones achieved by this method. As modulating waveforms we use the two options discussed in Section II-C: cosine and sine waveforms (tones), and filtered rectangular alternating pulses (rectangular). For all methods, we used only p = 5 sampling channels. We focus on one period of the input signal, which consists of 2 closely spaced diracs with delays t = [0.256T, 0.38T ]T , and amplitudes a = [1, 0.8]T . For each scheme the samples are obtained by analytic calculation: successive integration of the input signal for

10

20

30

40

50

60

70

80

90

SNR [dB]

Fig. 5. Performance in the presence of noise: delays MSE for two supported schemes vs. integrators approach [6].

[6], and integration of the modulated input signal in each channel for our approach. Then, Gaussian white noise with a given signal to noise ratio (SNR) is added to the samples. The estimation error of the time delays as a function of SNR is depicted in Fig. 5, for all methods. Evidently, our approach outperforms the integrators method [6] in terms of noise robustness, for both modulating waveforms. VI. C ONCLUSION In this work we proposed new sampling schemes for pulse streams. These methods allow perfect recovery of the sampled signals, while operating at the ROI. In contrast to previous works treating dirac streams only, our approach extends to more general pulse shapes. As we demonstrate by simulation, our method exhibits better noise robustness than the one presented in [6]. R EFERENCES [1] Y. C. Eldar and T. Michaeli, “Beyond bandlimited sampling,” IEEE Signal Process. Mag., vol. 26, no. 3, pp. 48–68, May 2009. [2] M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE Trans. Signal Process., vol. 50, no. 6, pp. 1417– 1428, Jun 2002. [3] P. L. Dragotti, M. Vetterli, and T. Blu, “Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets strang-fix,” IEEE Trans. Signal Process., vol. 55, no. 5, pp. 1741–1757, May 2007. [4] R. Tur, Y. C. Eldar, and Z. Friedman, “Low rate sampling of pulse streams with application to ultrasound imaging,” submitted to IEEE Trans. Signal Process. [5] C. Seelamantula and M. Unser, “A generalized sampling method for finite-rate-of-innovation-signal reconstruction,” IEEE Signal Process. Lett., vol. 15, pp. 813–816, 2008. [6] J. Kusuma and V. Goyal, “Multichannel sampling of parametric signals with a successive approximation property,” IEEE Int. Conf. Image Process. (ICIP2006), pp. 1265 –1268, Oct 2006. [7] K. Gedalyahu and Y. C. Eldar, “Time delay estimation from low rate samples: A union of subspaces approach,” IEEE Trans. Signal Process., vol. 58, no. 6, pp. 3017–3031, Jun 2010. [8] P. Stoica and R. Moses, Introduction to Spectral Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1997. [9] R. Tur, K. Gedalyahu, and Y. C. Eldar, “Multichannel sampling of pulse streams at the rate of innovation,” submitted to IEEE Trans. Signal Process. [10] M. Mishali and Y. C. Eldar, “From theory to practice: Sub-Nyquist sampling of sparse wideband analog signals,” IEEE J. Sel. Topics Signal Process., vol. 4, no. 2, pp. 375–391, Apr 2010. [11] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984–995, Jul 1989. [12] R. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 276–280, Mar 1986. [13] T.-J. Shan, M. Wax, and T. Kailath, “On spatial smoothing for directionof-arrival estimation of coherent signals,” IEEE Trans. on Acoustics, Speech and Signal Processing, vol. 33, no. 4, pp. 806–811, Aug 1985.