Compressed Sensing of Analog Signals in Shift-Invariant Spaces

1 Compressed Sensing of Analog Signals in Shift-In v ariant Spaces Y onina C. Eldar Abstract A traditional assumption under lying most data con verter s is that the signal should be sampled at a rate exceeding twice the highest frequen cy . This stateme nt is based o n a worst-case scenar io in which th e signa l occup ies the entire av ailable b andwidth . In practice , m any sign als are sparse so tha t only part of the b andwidth is used. I n this paper, we develop meth ods fo r low-rate sampling o f continuo us-time sparse signals in shift-inv arian t ( SI) spaces, generated by m kernels with period T . W e mo del spar sity by treating the case in wh ich on ly k out o f the m generato rs are active, h owe ver, we d o not know which k are chosen . W e sh ow how to sample such signals at a rate much lower th an m /T , which is the minim al samplin g rate witho ut explo iting sparsity . Ou r ap proach co mbines ideas from analog sampling in a subspace with a recently developed block diagr am th at converts an in finite set o f sparse equations to a finite coun terpart. Using these two c ompon ents we formu late our pr oblem within the fram ew o rk of finite compressed sensing (CS) an d then rely on algorithm s developed in that context. The disting uishing feature of ou r results is that in c ontrast to stan dard CS, which treats finite-leng th vectors, we consider sampling of ana log signals fo r which no und erlying finite-dim ensional model exists. The pr oposed framework allows to extend much of the rec ent literature o n CS to the ana log domain. I . I N T RO D U C T I O N Digital applications hav e developed rapidly over the las t few decad es. Signal process ing in the discrete domain inherently relies on sa mpling a continuous-time signa l to ob tain a discrete-time rep resentation. The traditional assumption u nderlying most a nalog-to-digital c on verters is tha t the samples must b e ac quired at the Shanno n- Nyquist rate, corresponding to twice the h ighest frequency [1], [2]. Although the bandlimited as sumption is often app roximately met, many sign als can be more adequ ately modeled in alternative bases other than the Fourier basis [3], [4], or pos sess further s tructure in the Fourier domain. Re search Department of Electrical Engineering , T echnion—Israel Institute of T echnology , Haifa 32000, Israel. Phone: +972-4 -8293256 , fax: +97 2- 4-829575 7, E-mail: yonina@ee.techn ion.ac.il. This work was supported in part by the Israel S cience Foundation under Grant no. 1081/07 and b y th e Euro pean Comm ission in the frame work of t he F P7 Network of Excellence in Wireless COMmunications NEW COM++ (contract no. 216715). in sampling theory over the past two de cades has s ubstantially enlarged the class of sampling problems that can be treated ef ficiently and reliably . Th is resu lted in many new sampling theories w hich acco mmodate mo re gen eral signal sets as well a s various linear an d nonlinear d istortions [5], [6 ], [4], [7], [8], [9], [10], [11]. A signal class that plays an impo rtant role in sampling theory are signals in sh ift-in v ariant (SI) spac es. Such functions can b e expres sed as linear combinations of shifts of a set of g enerators with pe riod T [12], [13], [14], [15], [16]. This mode l encompas ses many signals used in commun ication an d signal processing . For example, the s et of bandlimited functions is SI with a single ge nerator . Other examp les include splines [4], [17] and pu lse amplitude modulation in communications . Us ing multiple gen erators, a larger set of s ignals can be described su ch as mu ltiband func tions [18], [19], [20], [21], [22], [23]. Sa mpling theories similar to the Shann on theorem can b e developed for this s ignal c lass, which allo ws to sample and recons truct suc h functions using a broad variety o f filters. Any s ignal x ( t ) in a SI spac e ge nerated by m functions shifted with period T can be pe rfectly rec overed from m sampling seq uence s, obtained by filtering x ( t ) with a ba nk of m filters a nd uniformly sampling their o utputs at times nT . The overall sampling rate of such a system is m/T . In Section II we sh ow explicitly how to re cover x ( t ) from the se s amples b y an a ppropriate filter ban k. If the signal is generated b y k o ut o f the m gene rators, then as long as the c hosen subs et is known, it suffices to sample at a rate of k /T correspo nding to uniform s amples with p eriod T a t the output of k filters. Howe ver , a more difficult que stion is whether the rate can be redu ced if we know that only k of the generators are acti ve, but we do not k now in a dvance whic h ones . Since in principle x ( t ) may be comp rised of any of the generators, it may see m at first that the rate canno t be lower than m/T . This que stion is a spec ial c ase of sampling a signal in a union o f sub space s [24], [25], [26]. In our problem, x ( t ) lies in one of the sub space s expressed by k gene rators, howe ver we do not know which subspac e is c hosen . Necess ary a nd s ufficient con ditions where d eri ved in [24 ], [25] to ens ure that a sampling operator over such a union is in vertible. In our se tting this reduce s to the requirement that the sampling rate is at least 2 k /T . However , no c oncrete sa mpling metho ds where given that ens ure e f ficient and stable recovery , and n o rec overy algorithm was provided from a given set o f samp les. Finite-dimensional unions wh ere treated in [26], for which s table recovery methods where d ev eloped. Anothe r special case o f sampling on a union of sp aces that has been studied extensively is the p roblem unde rlying the field of compresse d sens ing (CS). In this s etting, the goal is to rec over a len gth m vector x from p < m linea r mea surements, whe re x is known to b e k -sparse in some basis [27], [28]. Many stable and ef fi cient recovery algorithms have been proposed to rec over x in this setting [29], [30], [31], [32], [28], [33], [26]. A fundamen tal difference betwee n o ur problem and mainstream CS papers is that we aim to sample and reconstruct continuous signals, while CS focuses o n recovery of finite vectors. The methods d ev eloped in the context of CS rely on the finite na ture o f the problem an d cannot be immediately ado pted to infinite-dimensional settings without disc retization or heuristics. Our goal is to direc tly red uce the ana log sampling rate, without first requiring the Nyquist-rate samp les and then applying finite-dimension al CS techn iques. Several attempts to extend CS ideas to the an alog domain were developed in a s et of conference s p apers [34], [35]. Howev e r , in both pap ers an underlying d iscrete model was as sumed which enab led immediate application of known CS tec hniques. A n alternati ve a nalog framew o rk is the work on finite rate of innov a tion [36], [37], in which x ( t ) is mod eled a s a finite linear combination o f sh ifted diracs (some extensions are given to other generators as well). The algorithms developed in this con text exploit the similarity between the giv en problem and spectral estimation, and again rely on finite dimension al methods. In con trast, the mo del we treat in this paper is inherently infinite dimens ional as it in volves a n infinite sequ ence of s amples from wh ich we would like to recover an analog signal with infinitely many parame ters. In such a s etting the measureme nt matrix in stand ard CS is replaced by a more gene ral linea r o perator . It is therefore no longe r c lear how to cho ose su ch an o perator to ens ure stability . Furthermore, even if a stable ope rator can be implemente d, it will result in infi nitely many compres sed samples . As stan dard CS a lgorithms operate o n finite-dimens ional optimization prob lems, they cannot be applied to infinite dimensional s equenc es. In ou r previous work, we considered a s parse an alog sampling problem in which the signal x ( t ) has a multiband structure, s o tha t its Fourier transform cons ists of at most N bands , each of width limited by B [21], [22 ], [23 ]. Explicit sub-Nyquis t sampling and reconstruction schemes we re dev eloped in [21], [22], [23] that ensure pe rfect recovery of multiband signa ls at the minimal p ossible rate, without requiring k nowledge of the ba nd locations. The proposed algorithms rely on a set of operations grouped under a block named continuous-to-finite (CTF). The CTF , which is further developed in [38], esse ntially trans forms the c ontinuous recon struction problem into a finite dimensional eq uiv alent, without discretization or he uristics. Th e resulting problem is formulated within the framew ork of CS, and thus c an be solved e f ficiently using known tractable a lgorithms. The sa mpling methods used in [21], [23] for blind multiband s ampling are tailored to that specific setting, and are not applicable to the more general mode l we con sider here. Our goa l in this paper is to ca pitalize on the key elements from [21], [38], [23 ] that enab le CS of multiband signals and extend them to the mo re general SI setting by combining res ults from standard s ampling theo ry a nd CS. Althou gh the ide as we present are rooted in our previous work, the ir a pplication to mo re ge neral a nalog CS is no t immediately obvious . T o extend our work, it is cruc ial to setup the more gen eral p roblem treated here in a pa rticular way . Therefore, a large part of the pape r is focus ed on the problem setup, and reformulation of previously de ri ved results. W e then show explicitly how signals in a SI union crea ted by m ge nerators with period T , can be sa mpled and stably recovered at a rate much lower tha n m/T . Specifically , if k out o f m generators a re active, the n it is sufficient to use 2 k ≤ p < m uniform sequen ces at rate 1 /T , where p is determined by the requirements of standard CS. The paper is organized a s follows. In Sec tion II we p rovide backg round material on Ny quist-rate s ampling in SI space s. Although most of these results are known in the literature, we review them here since our interpretation of the recovery method is e ssential in treating the sp arse setting. The sparse SI mod el is presented in Sec tion III. In this se ction we also review the ma in elements of CS needed for the development of our a lgorithm, and elab orate more on the ess ential difficulty in extend ing them to the analog setting. Th e d if feren ce b etween s ampling in genera l SI sp aces and our p revious work [21] is also h ighlighted. In Section V we present o ur s trategy for CS o f SI an alog signals. Some exa mples of our framew ork are discu ssed in Sec tion VI. I I . B AC K G RO U N D : S A M P L I N G I N S I S P AC E S T rad itional sa mpling theory d eals with the problem of recovering an unknown function x ( t ) ∈ L 2 from its uniform samp les at times t = nT where T is the sampling period. More generally , the signal may b e pre-filtered prior to sampling with a filter s ∗ ( − t ) [4], [7], [11] , [39], [16], where ( · ) ∗ denotes the complex conjugate, a s illustrated in the le ft-hand side of Fig. 1. The s amples c [ n ] can be represented as the inner products c [ n ] = x ( t ) ✲ s ∗ ( − t ) ✲ ✑ ✑ ❄ t = nT ✲ φ − 1 S A ( e j ω ) ✲ ❧ × ✲ a ( t ) ✲ x ( t ) c [ n ] d [ n ] P n ∈ Z δ ( t − nT ) ✻ Fig. 1. Non -ideal sampling and reconstruction. h s ( t − n T ) , x ( t ) i whe re h s ( t ) , x ( t ) i = Z ∞ t = −∞ s ∗ ( t ) x ( t ) dt. (1) In order to rec over x ( t ) from these s amples it is typically ass umed that x ( t ) lies in an appropriate subsp ace A of L 2 . A common choice of subsp ace is a SI subs pace ge nerated by a sing le ge nerator a ( t ) . Any x ( t ) ∈ A ha s the form x ( t ) = X n ∈ Z d [ n ] a ( t − nT ) , (2) for some ge nerator a ( t ) and s ampling pe riod T . Note that d [ n ] are n ot nec essarily pointwise samples o f the signal. If   φ S A ( e j ω )   > α > 0 , a.e. ω , (3) where we defined φ S A ( e j ω ) = 1 T X k ∈ Z S ∗  ω T − 2 π T k  A  ω T − 2 π T k  , (4) then x ( t ) ca n be p erfectly recons tructed from the sa mples c [ n ] in Fig. 1 [39], [6]. T he function φ S A ( e j ω ) is the discrete-time Fourier transform (DTF T) of the samp led cross-correlation se quence : r sa [ n ] = h s ( t − n T ) , a ( t ) i . (5) T o empha size the fact that the DTFT is 2 π -periodic we use the notation Φ ( e j ω ) . Re covery is obtained by filtering the s amples c [ n ] by a disc rete-time filter with frequency response H ( e j ω ) = 1 φ S A ( e j ω ) , (6) follo wed by modulation by an impulse train with p eriod T and fi ltering with a n an alog filter a ( t ) . The overall sampling and rec onstruction s cheme is illustrated in Fig. 1 . E vidently , SI sub spaces allow to retain the basic flavor of the Shannon sampling the orem in wh ich sampling and recovery are implemented by filtering o perations. In this pa per we co nsider more gene ral SI space s, generated by m functions a ℓ ( t ) ∈ L 2 , 1 ≤ ℓ ≤ m . A finitely-generated SI subspac e in L 2 is de fined as [12], [13], [15] A = ( x ( t ) = m X ℓ =1 X n ∈ Z d ℓ [ n ] a ℓ ( t − nT ) : d ℓ [ n ] ∈ ℓ 2 ) . (7) The functions a ℓ ( t ) are referred to as the generators of A . In the Fourier d omain, we c an represent any x ( t ) ∈ A as X ( ω ) = m X ℓ =1 D ℓ ( e j ωT ) A ℓ ( ω ) , (8) where D ℓ ( e j ω ) = X n ∈ Z d ℓ [ n ] e j ωn , (9) is the DT FT of d ℓ [ n ] . Througho ut the pap er , we u se upp er- case letters to denote Fourier transforms: X ( ω ) is the continuous -time Fourier trans form of the function x ( t ) , and C ( e j ω ) is the DTFT o f the sequen ce c [ n ] . In orde r to guarantee a unique s table repres entation of a ny signal in A by coe f ficients d ℓ [ n ] , the generators a ℓ ( t ) are typically chosen to form a Rie sz bas is for L 2 . This means that there exists constants α > 0 and β < ∞ such that α k d k 2 ≤      m X ℓ =1 X n ∈ Z d ℓ [ n ] a ℓ ( t − nT )      2 ≤ β k d k 2 , (10) where k d k 2 = P m ℓ =1 P n ∈ Z | d ℓ [ n ] | 2 , and the norm in the middle term is the standard L 2 norm. By taking Fourier transforms in (10) it follows that the gen erators a ℓ ( t ) form a Riesz b asis if and only if [13] α I  M AA ( e j ω )  β I , a.e. ω , (11) where M AA ( e j ω ) =       φ A 1 A 1 ( e j ω ) . . . φ A 1 A m ( e j ω ) . . . . . . . . . φ A m A 1 ( e j ω ) . . . φ A m A m ( e j ω )       . (12) Here φ A i A ℓ ( e j ω ) is de fined by (4) with A i ( e j ω ) , A ℓ ( e j ω ) replacing S ( e j ω ) , A ( e j ω ) . Throug hout the paper we assume that (11) is satisfied. Since x ( t ) lies in a space ge nerated by m func tions, it makes sense to sa mple it with m filters s ℓ ( t ) , a s in the left-hand side of Fig. 2. The sa mples are given b y ✲ s ∗ m ( − t ) ✲ ✑ ✑ ✑ ❄ t = nT ✲ ✲ ♥ × ✲ a m ( t ) ✲ c m [ n ] P n ∈ Z δ ( t − nT ) ✻ d m [ n ] . . . . . . ✲ s ∗ 1 ( − t ) ✲ ✑ ✑ ✑ ❄ t = nT ✲ ✲ ♥ × ✲ a 1 ( t ) ✲ c 1 [ n ] P n ∈ Z δ ( t − nT ) ✻ d 1 [ n ] x ( t ) ✲ M − 1 S A ( e j ω ) ♥ ✲ x ( t ) ✻ ❄ Fig. 2. Sampling and reconstruction in shift-i n variant spaces. c ℓ [ n ] = h s ℓ ( t − nT ) , x ( t ) i , 1 ≤ ℓ ≤ m. (13) The following prop osition provides a simple Fourier-domain relationship between the samp les c ℓ [ n ] of (13) and the expa nsion coe f ficients d ℓ [ n ] of (7): Pr opo sition 1: Let c ℓ [ n ] = h s ℓ ( t − nT ) , x ( t ) i , 1 ≤ ℓ ≤ m be a s et of m seque nces obtained by filtering x ( t ) of (7) with m filters s ∗ ℓ ( − t ) and s ampling their outputs at times nT , as depicted in the left-hand side of Fig. 2. Denote by c ( e j ω ) , d ( e j ω ) the vectors with ℓ th elements C ℓ ( e j ω ) , D ℓ ( e j ω ) , respe ctiv e ly . Then c ( e j ω ) = M S A ( e j ω ) d ( e j ω ) , (14) where M S A ( e j ω ) =       φ S 1 A 1 ( e j ω ) . . . φ S 1 A m ( e j ω ) . . . . . . . . . φ S m A 1 ( e j ω ) . . . φ S m A m ( e j ω )       , (15) and φ S i A ℓ ( e j ω ) are de fined by (4). Pr oof: The p roof follows immediately by taking the Fourier transform of (13): C ℓ ( e j ω ) = 1 T X k ∈ Z S ∗ ℓ  ω T − 2 π T k  X  ω T − 2 π T k  = 1 T m X i =1 D i ( e j ω ) X k ∈ Z S ∗ ℓ  ω T − 2 π T k  A i  ω T − 2 π T k  , (16) where we used (8) a nd the fact that D ( e j ω ) is 2 π -periodic. In vector from, (16) reduces to (14 ). Proposition 1 can be us ed to recover x ( t ) from the giv en samples as long as M S A ( e j ω ) is in vertible a.e. in ω . Under this cond ition, the expa nsion coefficients can be computed as d ( e j ω ) = M − 1 S A ( e j ω ) c ( e j ω ) . Given d ℓ [ n ] , x ( t ) is formed by modulating eac h coe f ficient seq uence by a p eriodic impuls e train P n ∈ Z δ ( t − nT ) with period T , and filtering with the corresp onding a nalog filter a ℓ ( t ) . In order to ensure stable recovery we require α I  M S A ( e j ω )  β I a.e. In particular , we may choo se s ℓ ( t ) = a ℓ ( t ) d ue to (11). The resulting s ampling and reconstruction sche me is dep icted in Fig. 2. The approach of Fig. 2 results in m sequenc es of s amples, eac h at rate 1 /T , leading to an average s ampling rate of m/T . Note that from (7 ) it follows that in each time step T , x ( t ) c ontains m new pa rameters, so that the signal has m d egrees of freedom over ev ery interval of length T . Therefore this sampling strategy ha s the intuiti ve property that it requires o ne samp le for eac h degree of freedom. I I I . U N I O N O F S H I F T - I N V A R I A N T S U B S P A C E S Evidently , when subsp ace information is av a ilable, perfect reco nstruction from linear samples is often ach iev able. Furthermore, re covery is pos sible using a simple filter bank. A more interesting sc enario is whe n x ( t ) lies in a union of SI subs paces of the form [26] x ( t ) ∈ [ | ℓ | = k A ℓ , (17) where the n otation | ℓ | = k mea ns a union (or s um) over at most k elemen ts. Here we con sider the case in wh ich the un ion is over k out of m possible subs paces {A ℓ , 1 ≤ ℓ ≤ m } , where A ℓ is ge nerated by a ℓ ( t ) . Thus , x ( t ) = X | ℓ | = k X n ∈ Z d ℓ [ n ] a ℓ ( t − nT ) , (18) so that on ly k out of the m sequenc es d ℓ [ n ] in the s um (18) are no t ide ntically zero. No te that (18) no longer defines a sub space . Indeed, while the un ion c ontains signals of the form P n d i [ n ] a i ( t − nT ) for two distinct values of i , it doe s not include their linear c ombinations. In principle, if we know wh ich k seque nces are non-zero, the n x ( t ) can be recovered from s amples at the ou tput of k filters using the scheme of Fig. 2. The resulting average s ampling rate is k /T s ince we hav e k sequen ces, each at rate 1 /T . Alternativ ely , even without knowledge of the active subspa ces, we c an rec over x ( t ) from s amples at the o utput of m filters resulting in a sampling rate of m/T . Although this strategy does not require knowledge of the activ e subspac es, the price is an increase in sampling rate. In [24], [25], the authors dev eloped nec essary and sufficient con ditions for a samp ling o perator to b e in vertible over a union of subspa ces. Spec ializing the res ults to our problem implies that a minimal rate of at least 2 k/T is needed in order to ens ure that the re is a uniqu e SI signal co nsistent with the samples . Thus, the fact that we do not know the exact subs pace lead s to an increa se of at least a factor two in the minimal rate. However , no co ncrete methods were p rovided to recons truct the original sign al from its samples. Furthermore, althoug h conditions for in vertibility we re provided, these do not nece ssarily imply that a stable and efficient recovery is poss ible at the minimal rate. Our goa l is to develop algorithms for rec overing x ( t ) from a set of 2 k ≤ p < m sa mpling sequen ces, obtained by sampling the outputs of p filters at rate 1 /T . Before developing ou r sampling scheme, we fi rst explain the dif ficulty in address ing this problem and its relation to prior work. A. Compressed Sensing A sp ecial cas e of a union of s ubspa ces that has been treated extensiv ely is CS of finite vectors [27], [28]. In this setup, the problem is to recover a fin ite-dimensional vector x of length m from p < m linear measuremen ts y where y = Mx , (19) for s ome ma trix M of size p × m . S ince the equ ations (19) are u nderdetermined, more information is needed in order to rec over x . The prior as sumed in the CS literature is that x = Φ α where Φ is an m × m in vertible ma trix, and α is k -sparse , so that it ha s a t most k non-zero e lements. Th is prior can be viewed as a union of subspac es where eac h subspa ce is span ned by k columns of Φ . A sufficient condition for the uniquene ss o f a k -sparse solution α to the e quations y = MΦ α is that A = MΦ has a Krusk al-rank of at least 2 k [32], [40]. The Krusk al-rank is the maximal n umber q such that every s et of q columns of A is linearly ind epende nt [41]. This unique α ca n be rec overed by solving the optimization problem [27]: min α k α k 0 s . t . y = A α , (20) where the pseudo -norm ℓ 0 counts the numbe r of n on-zero entries. Therefore, if we a re n ot conc erned with stability and compu tational c omplexity , then 2 k measurements a re enou gh to recover α exac tly . Since (20) is known to be NP-ha rd [27],[28], s everal alternative algorithms ha ve been propo sed in the literature tha t have p olynomial complexity . T wo prominent approac hes are to replace the ℓ 0 norm by the c on vex ℓ 1 norm, and the orthogo nal matching pursuit algorithm [27], [28]. For a given sparsity le vel k , the se techniques are guaranteed to recover the true s parse vector α as long as ce rtain conditions on A are s atisfied, such as the restricted isometry property [28], [42]. The efficient methods propose d to recover x all require a nu mber o f meas urements p that is larger than 2 k , however still considerab ly sma ller than m . For example, if A is cho sen as p random ro ws from the Fourier transform matrix, then the ℓ 1 program will recover α with overwhelming probability as long a s p ≥ ck log m whe re c is a co nstant. Other ch oices for A are random matrices consisting of Gaus sian or Bernoulli random vari ables [33], [43]. In these c ases, on the o rder of k log( m/k ) mea surements are neces sary in order to be ab le to recover α efficiently with high probability . These results have also been generalized to the mult iple-measuremen t vector (MMV) model in which the problem is to recover a matrix X from matrix measureme nts Y = AX , where X has at most k non-zero ro ws. Here again, if the Kruskal-rank of A is a t lea st 2 k , then there is a unique X cons istent with Y . Th is uniqu e solution c an b e obtained by solving the c ombinatorial prob lem min X | I ( X ) | , s . t . Y = AX , (21) where I ( X ) is the set of indices corresponding to the non-zero rows of X [40]. V arious efficient algorithms that coincide with (21) unde r certain cond itions on A hav e also bee n propos ed for this p roblem [44],[40],[38]. B. Compressed Sensing of Analog Signals Our problem is similar in spirit to finite CS: we would like to s ense a spars e sign al using fewer me asuremen ts than req uired without the s parsity a ssumption. Howe ver , the fundamental difference between the two stems from the fact that our problem is defined over a n infinite-dimensional s pace o f continuous functions. As we now s how , trying to represent it in the same form as CS by replacing the finite matrices by appropriate operators, raise s several difficulti es that precludes direct application of CS-type results. T o s ee this, suppose we rep resent x ( t ) in terms of a sparse expansion, by de fining an infin ite-dimensional operator Φ( t ) co rresponding to the c oncatena tion of the func tions a ℓ ( t − nT ) , a nd an infinite se quence α ∈ ℓ 2 which cons ists of the c oncatena tion of the sequ ences d ℓ [ n ] . W e may the n write x ( t ) = Φ( t ) α which rese mbles the finite expansion x = Φ α . Since d ℓ [ n ] is identically ze ro for several values o f ℓ , α will contain many zero elements. Next, we can de fine a meas urement operator M ( t ) so that the mea surements are given by y = Aα where A = M ( t )Φ( t ) . In ana logy to the fin ite setting, the recovery properties of α s hould depen d on A . Howe ver , immediate application of CS ideas to this op erator equation is impossible. As we have see n, the ab ility to rec over α in the fin ite s etting depend s on its s parsity . In o ur ca se, the spa rsity of α is a lways infinite. Furthermore, a practica l way to ensure stable recovery with high probability in c on ventional CS is to draw the elements of A at rand om, with the number of rows of A p roportional to the sparsity . In the operator se tting, we can not clearly define the dimensions of A or d raw its elements at random. Even if we can develop c onditions on A suc h that the measure ment s equenc e y = Aα un iquely dete rmines α , it is s till not c lear how to recover α from y . The immediate extension of bas is pursuit to this con text would be: min α k α k 1 s . t . y = Aα. (22) Although (22) is a c on vex problem, it is defined over infinitely many variables, with infinitely many constraints. Con vex programming techn iques such a s semi-infinite prog ramming a nd generalized s emi-infinite programming, allow only for infinite con straints while the optimization vari able mu st be finite. Therefore, (22) canno t be solved using stan dard optimization tools as in finite-dimensional CS. This discuss ion raises three important que stions we nee d to address in orde r to ad apt CS results to the analog setting: 1) How do we ch oose an a nalog sampling o perator? 2) Can we introduce structure into the s ampling ope rator and still preserve stability? 3) How do we so lve the resulting infin ite-dimensional recovery p roblem? C. Pr eviou s W ork on Analog Compressed Se nsing In o ur p revious work [21], [22], [23 ], we treated a special case of analog CS. S pecifically , we conside red blind multiband sampling in which the goal is to sample a bandlimited signal x ( t ) whose frequ ency res ponse consists of at mo st N ba nds of length limited by B , with un known supp ort. In Section VI we show that this model can be formulated as a union of SI subspac es (18). In order to sample and recover x ( t ) at rates much lower than Nyquist, we proposed two types of s ampling method s: In [21] we c onsidered multicoset sa mpling while in [23] we modu lated the signal by a periodic function, prior to standard lo w-pa ss sampling. Bo th s trategies lead to p sequen ces of low-rate uniform samples which, in the Fourier domain, ca n be related to the unknown x ( t ) via an infinite measureme nt vector (IMV) model [38], [21]. This mo del is a n exten sion of the MMV p roblem to the c ase in which the goal is to recover infinitely ma ny unknown vectors that share a joint sp arsity pattern. Using the IMV techniques developed in [21], [38], x ( t ) can then be rec overed by solving a finite-dimensional con vex optimization problem. W e elab orate more on the IMV mod el below , as it is a key ingredient in ou r propos ed sa mpling strategy . The s ampling method d escribed above is tailored to the multiband mo del, an d exploits the fact that the spectrum has many intervals that are identically zero. Applying this approach to the general SI setting will not lead to pe rfect recovery . In order to extend ou r previous results, we the refore ne ed to reveal the key ideas that allow CS of analog signals, rather than analyzing a spe cific set of s ampling equa tions a s in [21], [22], [23]. Further study of blind multiband sampling suggests two key elements en abling analog CS: 1) Fourier domain a nalysis of the sequenc es of samples . 2) Choosing the sampling functions such that we obtain a n IMV mode l. Our app roach is to capitalize on these two compon ents and extend them to the mo del (18). T o develop an a nalog CS s ystem, we design p < m sampling filters s i ( t ) that enable p erfect recovery of x ( t ) . In view of our previous discuss ion, our task ca n be rephrased as dete rmining sampling filters suc h tha t in the Fourier domain, the resu lting sa mples can be des cribed b y an IMV s ystem. In the next sec tion we revie w the key elements of the IMV problem. W e then s how how to app ropriately choos e the sampling filters for the general mode l (18). I V . I N FI N I T E M E A S U R E M E N T M O D E L In the IMV model the goal is to recover a se t of unk nown vectors x ( λ ) from meas urement vectors y ( λ ) = Ax ( λ ) , λ ∈ Λ , (23) where Λ is a s et whose cardina lity can b e infinite. In particular , Λ may be uncoun table, suc h as the frequencies ω ∈ [ − π , π ) . Th e k -spa rse IMV model as sumes that the vectors { x ( λ ) } , which we deno te for brevity by x (Λ) , share a joint spa rsity pattern, that is, the no n-zero elemen ts are sup ported on a fixed location set of size k [38]. As in the finite-cas e it is easy to see that if σ ( A ) ≥ 2 k , wh ere σ ( A ) is the Kruskal-rank of A , then x (Λ) is the unique k -sp arse solution of (23) [38]. The major difficulty with the IMV model is how to rec over the solution set x (Λ) from the infinitely many equations (23). One subop timal strategy is to con vert the problem into an MMV by solving (23) only over a finite set of values λ . Howe ver , clearly this strategy cannot g uarantee perfect rec overy . Instead, the ap proach in [38] is to recover x (Λ) in two steps. F irst, we fin d the su pport se t S of x (Λ) , and then reconstruct x (Λ) from the data y (Λ) and k nowledge o f S . Once S is found, the secon d step is straightforward. T o see this, no te that us ing S , (23 ) can be written a s y ( λ ) = A S x S ( λ ) , λ ∈ Λ , (24) where A S denotes the matrix con taining the columns of A whose indices be long to S , a nd x S ( λ ) is the vector consisting of entries of x ( λ ) in locations S . Since x (Λ) is k -sparse, | S | ≤ k . Therefore, the columns of A S are linearly independ ent (beca use σ ( A ) ≥ 2 k ), implying that A † S A S = I , where A † S =  A H S A S  − 1 A H S is the pseudo -in verse of A S and ( · ) H denotes the He rmitian conjuga te. Multiplying (24) by A † S on the left gives x S ( λ ) = A † S y ( λ ) , λ ∈ Λ . (25) The comp onents in x ( λ ) not su pported on S are all zero. Therefore (25) a llo ws for exact recovery of x (Λ) once the fin ite set S is correc tly identified. It rema ins to determine S efficiently . In [38] it was shown that S can be fou nd exactly by solving a finite MMV . The steps used to formulate this MMV are grouped under a block referred to as the continuous -to-finite (CTF) block. The esse ntial idea is that every finite c ollection of vectors spann ing the s ubspac e span( y (Λ)) contains sufficient information to recover S , as incorporated in the follo wing theo rem [38]: Theorem 1: Suppo se that σ ( A ) ≥ 2 k , an d let V be a matrix with column span equal to span( y (Λ)) . The n, the linear sy stem V = AU (26) has a unique k -spa rse solution U whose s upport is eq ual S . The advantage of The orem 1 is that it allows to avoid the infi nite structure of (23) an d ins tead find the finite set S by solving the single MMV s ystem of (26). The add itional requirement of The orem 1 is to cons truct a matrix V having column s pan equal to span( y (Λ)) . Th e following p roposition, p roven in [38], sugges ts such a proce dure. T o this end , we as sume that x (Λ) is pie cewise co ntinuous in λ . y (Λ) Find a frame for y (Λ) Reconstruct the set S = I ( x (Λ)) Q = R λ ∈ Λ y ( λ ) y H ( λ ) dλ Q = VV H V S = I ( ¯ U ) S Prop os ition 2 Theorem 1 Solv e MMV V = AU for sparsest solution m atrix ¯ U Fig. 3. The fundamental stages for the recov ery of the non-zero location set S in an IMV model using only one finite-dimensional problem. Pr opo sition 2: If the integral Q = Z λ ∈ Λ y ( λ ) y H ( λ ) dλ, (27) exists, then every matrix V satisfying Q = VV H has a column span equal to span( y (Λ)) . Fig. 3, taken from [38], summarizes the reduction steps tha t follow from Theo rem 1 a nd Proposition 2. Note, that ea ch bloc k in the figu re ca n be replaced by another s et of ope rations having an equiv alen t func tionality . In particular , the c omputation of the matrix Q o f Prop osition 2 can b e avoided if a lternati ve methods a re e mployed for the construction o f a frame V for span( y (Λ)) . In the figure, I ind icates the joint sup port set of the correspond ing vectors. V . C O M P R E S S E D S E N S I N G O F S I S I G N A L S W e now combine the ide as of Sections II and IV in order to develop efficient s ampling strategies for a union of s ubspac es of the form (18). Our app roach con sists of filtering x ( t ) with p < m filters s i ( t ) , and un iformly sampling their outputs at rate 1 /T . The de sign of s i ( t ) , 1 ≤ i ≤ p relies on two ing redients: 1) A matrix A cho sen s uch that it solves a discrete CS problem in the dimensions m (vector length) and k (sparsity). 2) A set of functions h i ( t ) , 1 ≤ i ≤ m wh ich c an b e used to sample a nd reconstruct the entire set of generators a i ( t ) , 1 ≤ i ≤ m , namely s uch that M H A ( e j ω ) is stably in vertible a.e . The matrix A is de termined by co nsidering a finite-dimensiona l CS problem where we would like to recover a k -sparse vector x of length m from p measuremen ts y = Ax . The value of p c an be chose n to guarantee exact recovery with combinatorial optimization, in which c ase p ≥ 2 k , or to lead to efficient recovery (poss ibly only with h igh p robability) req uiring p > 2 k . W e show b elow that the s ame A cho sen for this discrete problem c an be used for analog CS. The functions h i ( t ) are chosen so that they ca n be u sed to recover x ( t ) . Howe ver , since there a re m s uch functions, this results in more meas urements than actually needed . W e de ri ve the p roposed sampling scheme in three s teps: First, we cons ider the problem of compressively measuring the vector se quence d [ n ] , whose ℓ th compo nent is gi ven by d ℓ [ n ] , wh ere only k out o f the m seque nces d ℓ [ n ] a re non-ze ro. W e show that this ca n be ac complished by using the matrix A above a nd IMV recovery theory . In the second step, we obtain the vector seq uence d [ n ] from the gi ven s ignal x ( t ) using a n appropriate filter bank of m an alog filters, and sa mpling their outpu ts. Finally , we merge the first two steps to arri ve a t a bank o f p < m analog filters tha t can compressively sample x ( t ) directly . These steps are detailed in the 3 ensuing subs ections. A. Union o f Discrete Se quence s W e begin by treating the problem of sampling and recovering the sequenc e d [ n ] . This can be acco mplished by using the IMV model introduced in S ection IV. Indeed, su ppose we measu re d [ n ] with a s ize p × m matrix A , that allows for CS of k -spa rse vectors of leng th m . Th en, for e ach n , y [ n ] = Ad [ n ] , n ∈ Z . (28) The sys tem of (28) is an IMV mode l: For every n the vector d [ n ] is k -sparse. Fu rthermore, the infinite set o f vectors { d [ n ] , n ∈ Z } has a joint sparsity p attern s ince at most k of the s equen ces d ℓ [ n ] are non -zero. As we described in Sec tion IV, suc h a system of e quations c an be solved by trans forming it into an equiv alent MMV , whose recovery prop erties a re determined by those of A . Since A was designe d such that CS tec hniques will work, we are gu aranteed that d [ n ] can be perfectly recovered for each n (or recovered with high probability). The reconstruction algorithm is depicted in Fig. 3. Note that in this case the integral in c omputing Q becomes a s um over n : Q = P n ∈ Z y [ n ] y H [ n ] (we assume here that the su m exists). Instead of solving (28) we may also cons ider the Freque ncy-domain set of equa tions: y ( e j ω ) = Ad ( e j ω ) , 0 ≤ ω < 2 π , (29) where y ( e j ω ) , d ( e j ω ) are the vectors whose compo nents are the frequency response s Y ℓ ( e j ω ) , D ℓ ( e j ω ) . In principle, we may ap ply the CTF b lock of Fig. 3 to e ither representations, dep ending on which ch oice off ers a simpler method for determining a basis V for the range of { y (Λ) } . When d esigning the measurements (28), the only freedom we have is in ch oosing A . T o gene ralize the class of sensing ope rators we note that d [ n ] c an also be recovered from y ( e j ω ) = W ( e j ω ) Ad ( e j ω ) , 0 ≤ ω < 2 π , (30) for any in vertible p × p matrix W ( e j ω ) with elements W iℓ ( e j ω ) . T he me asuremen ts of (30) c an be ob tained d irectly in the time domain as y i [ n ] = p X ℓ =1 w iℓ [ n ] ∗ m X r =1 A ℓr d r [ n ] ! , 1 ≤ i ≤ p , (31) where w iℓ [ n ] is the in verse transform of W iℓ ( e j ω ) , and ∗ de notes the c on volution o perator . T o recover d ℓ [ n ] from y ( e j ω ) , we note that the modified mea surements ˜ y ( e j ω ) = W − 1 ( e j ω ) y ( e j ω ) obey an IMV model: ˜ y ( e j ω ) = Ad ( e j ω ) , 0 ≤ ω < 2 π . (32) Therefore, the CTF block can be a pplied to ˜ y ( e j ω ) . As in (31), we may use the CTF in the time doma in b y n oting that ˜ y i [ n ] = p X ℓ =1 b iℓ [ n ] ∗ y ℓ [ n ] , (33) where b iℓ [ n ] is the in verse DTF T of B iℓ ( e j ω ) , with B ( e j ω ) = W − 1 ( e j ω ) . The extra freedom offered b y choo sing an arbitrary in vertible matrix W ( e j ω ) in (30) will be useful when we discuss an alog sampling, a s different choices lea d to dif feren t samp ling functions. In Section VI we will see an example in which a proper s election of W ( e j ω ) leads to ana log sampling functions that are ea sy to implement. B. Biorthogonal Expans ion The previous se ction estab lished that given the ability to samp le the m se quence s d ℓ [ n ] we can recover them exactly from p < m discrete-time seque nces obtaine d via (30) or (31). Re construction is performed by applying the CTF b lock to the modified measuremen ts e ither in the frequency doma in (32) or in the time doma in (33). The drawback is that we do no t have a ccess to d ℓ [ n ] but rather we are given x ( t ) . In Fig. 2 and Section II we have se en that the seq uences d ℓ [ n ] can be ob tained by sampling x ( t ) with a set of functions h ℓ ( t ) for which M H A ( e j ω ) of (15) is s tability in vertible, and then filtering the sampled sequen ces with a mu ltichannel disc rete-time filter M − 1 H A ( e j ω ) . Thus, we can first apply this front-end to x ( t ) , wh ich will produce the seq uence of vectors d [ n ] . W e ca n the n use the results of the p revious subs ection in order to s ense these seque nces efficiently . The resulting measureme nt seq uences y ℓ [ n ] are depicted in Fig. 4, whe re A is a p × m matrix s atisfying the requiremen ts of CS in the appropriate dimension s, and W ( e j ω ) is a size p × p filter b ank that is in vertible a.e. Combining the analog filters h ℓ ( t ) with the discrete-time multichannel filter M − 1 H A ( e j ω ) , we can express d ℓ [ n ] as d ℓ [ n ] = h v ℓ ( t − nT ) , x ( t ) i , 1 ≤ ℓ ≤ m, n ∈ Z , (34) where v ( ω ) = M −∗ H A ( e j ωT ) h ( ω ) . (35) Here v ( ω ) , h ( ω ) are the vectors with ℓ th elemen ts V ℓ ( ω ) , H ℓ ( ω ) and ( · ) −∗ denotes the con jugate of the in verse. ✲ h ∗ m ( − t ) ✲✑ ✑ ❄ t = nT ✲ ✲ y p [ n ] . . . . . . . . . ✲ h ∗ 1 ( − t ) ✲ ✑ ✑ ❄ t = nT ✲ ✲ y 1 [ n ] x ( t ) ✲ M − 1 H A ( e j ω ) ✲ d m [ n ] ✲ d 1 [ n ] A W ( e j ω ) Fig. 4. Ana log compressed sampling with arbitrary filters h i ( t ) . The inn er produ cts in (34) can be obtained by filtering x ( t ) with the bank o f filters v ∗ ℓ ( − t ) , an d un iformly samp ling the ou tputs at times nT . T o see that (34) holds, let c ℓ [ n ] be the samp les res ulting from filtering x ( t ) with the m filters v ℓ ( t ) and uniformly sampling their outputs at rate 1 /T . From Propo sition 1, c ( e j ω ) = M V A ( e j ω ) d ( e j ω ) . (36) Therefore, to establish (34) we need to show tha t M V A ( e j ω ) = I . Now , from (15), [ M V A ( e j ω )] iℓ = 1 T X k ∈ Z V ∗ i  ω T − 2 π T k  A ℓ  ω T − 2 π T k  = 1 T m X r =1 [ M − 1 H A ( e j ω )] ir X k ∈ Z H ∗ r  ω T − 2 π T k  A ℓ  ω T − 2 π T k  = [ M − 1 H A ( e j ω )] i [ M H A ( e j ω )] ℓ = I iℓ , (37) where [ Q ] ir is the ir th eleme nt of the matrix Q , and [ Q ] i , [ Q ] i are the i th row an d c olumn respec ti vely of Q . Therefore, as required, M V A ( e j ω ) = I . The functions { v ℓ ( t − nT ) } have the prope rty that they are b iorthogonal to { a ℓ ( t − nT ) } , tha t is h v ℓ ( t − nT ) , a i ( t − r T ) i = δ ℓi δ nr , (38) where δ ℓi = 1 if ℓ = i , and 0 otherwise. Th is follows from the f act that in the Fourier domain, (38) is equ i valent to M V A ( e j ω ) = I . (39) Evidently , we can c onstruct a set of b iorthogonal fun ctions from any s et of functions h ℓ ( t ) for whic h M H A ( e j ω ) is stably in vertible, via (35). Note that the biorthogonal vectors in the space A are uniqu e. This follows from the fact that if two sets { v 1 i ( t ) } , { v 2 i ( t ) } sa tisfy (38), then h g ℓ ( t − nT ) , a i ( t − r T ) i = 0 , for all ℓ, i, n, r , (40) where g i ( t ) = v 1 i ( t ) − v 2 i ( t ) . Since { a i ( t − nT ) } span A , (40) implies that g i ( t − mT ) lies in A ⊥ for any i, m . Howe ver , if both v 1 i ( t ) and v 2 i ( t ) are in A , then so is g i ( t − mT ) , from which we conclud e that g i ( t − mT ) = 0 . Thus, a s long as we start with a se t of functions h i ( t ) tha t sp an A , the sa mpling functions v i ( t ) re sulting from (35) will be the sa me. Howe ver , their implementation in h ardware is dif fe rent, since h i ( t ) represen ts an analog filter while M H A ( e j ω ) is a discrete-time filter bank. Th erefore, d if feren t cho ices of h i ( t ) lead to dis tinct ana log filters. C. CS of Analog Signals Although the sampling sc heme of Fig. 4 resu lts in compress ed measu rements y ℓ [ n ] , they are s till obta ined by an analog front-end that operates a t the h igh rate m/T . Howe ver , o ur g oal is to red uce the rate at the an alog front-end. This can be ea sily acc omplished by moving the discre te filters M − 1 H A ( e j ω ) , A W ( e j ω ) back to the an alog domain. In this way , the co mpressed meas urement sequenc es y ℓ [ n ] can b e obtained d irectly from x ( t ) , by filtering x ( t ) with p filters s ℓ ( t ) and uniformly sampling their outputs at times nT , leading to a system with sampling rate p/T . An explicit expression for the resulting s ampling functions is given in the follo w ing theo rem. Theorem 2: Let the compres sed measuremen ts y ℓ [ n ] , 1 ≤ ℓ ≤ p be the output of the hy brid filter ban k in Fig. 4. Then { y ℓ [ n ] } c an be obtained by filtering x ( t ) of (18) with p filters { s ∗ ℓ ( − t ) } a nd sampling the outputs at rate 1 /T , where s ( ω ) = W ∗ ( e j ωT ) A ∗ v ( ω ) = W ∗ ( e j ωT ) A ∗ M −∗ H A ( e j ωT ) h ( ω ) . (41) Here s ( ω ) , h ( ω ) are the vectors with ℓ th e lements S ℓ ( ω ) , H ℓ ( ω ) respec ti vely , and the compon ents V ℓ ( ω ) of v ( ω ) = M −∗ H A ( e j ωT ) h ( ω ) a re F o urier transforms of generators v ℓ ( t ) suc h that { v ℓ ( t − nT ) } a re biorthogonal to { a ℓ ( t − nT ) } . In the time domain, s i ( t ) = m X ℓ =1 p X r =1 X n ∈ Z w ∗ ir [ − n ] A ∗ r ℓ v ℓ ( t − nT ) , (42) where w ir [ n ] is the in verse trans form of [ W ( e j ω )] ir and v i ( t ) = m X ℓ =1 X n ∈ Z ψ ∗ iℓ [ − n ] h ℓ ( t − nT ) , (43) where ψ iℓ [ n ] is the in verse trans form of [ M − 1 H A ( e j ω )] iℓ . Pr oof: Suppos e tha t x ( t ) is filtered by the p filters s i ( t ) and then unifor mly sampled at nT . From Proposition 1, the s amples can be expressed in the Fourier do main as c ( e j ω ) = M S A ( e j ω ) d ( e j ω ) . (44) In order to prove the theorem w e need to show tha t M S A ( e j ω ) = W ( e j ω ) A . Let B ( e j ω ) = W ∗ ( e j ω ) A ∗ , (45) so that s ( ω ) = B ( e j ωT ) v ( ω ) . Then , [ M S A ( e j ω )] iℓ = 1 T X k ∈ Z S ∗ i  ω T − 2 π T k  A ℓ  ω T − 2 π T k  = 1 T m X r =1 B ∗ ir ( e j ω ) X k ∈ Z V ∗ r  ω T − 2 π T k  A ℓ  ω T − 2 π T k  = [ B ∗ ( e j ω )] i [ M V A ( e j ω )] ℓ , (46) where [ Q ] i , [ Q ] i are the i th row a nd column respectively of the matrix Q . Th e first equality follows from the fact that B ( e j ω ) is 2 π periodic. From (46), M S A ( e j ω ) = B ∗ ( e j ω ) M V A ( e j ω ) = W ( e j ω ) A , (47) where we used the fact tha t M V A ( e j ω ) = I due to the biorthogon ality property . Finally , if s ( ω ) = B ( e j ωT ) v ( ω ) , the n s i ( t ) = m X ℓ =1 X n ∈ Z b iℓ [ n ] v ℓ ( t − nT ) , (48 ) where b iℓ [ n ] is the in verse DTFT of [ B iℓ ( e j ω )] iℓ . Using (45 ) together with the fact that the in verse transform of Q ∗ iℓ ( e j ω ) is q ∗ iℓ [ − n ] , res ults in (42). The relation (43) follows from the same conside rations. Theorem 2 is the main res ult which allows for compres siv e samp ling of analog signa ls. Specifica lly , starting from any matrix A that satisfie s the CS requ irements o f finite vectors, and a set of sampling functions h i ( t ) for which M H A ( e j ω ) is in vertible, we ca n c reate a multitude of sa mpling fun ctions s i ( t ) to compres siv ely sample the underlying a nalog signal x ( t ) . The sens ing is performed b y fi ltering x ( t ) with the p < m correspo nding filters, and samp ling their outpu ts at rate 1 /T . Recons truction from the compressed measureme nts y i [ n ] , 1 ≤ i ≤ p is obtained by applying the CT F block of F ig. 3 in order to rec over the s equenc es d i [ n ] . The original s ignal x ( t ) is then con structed by modulating appropriate impulse trains and filtering with a i ( t ) , as depicted in Fig. 5. ✲ s ∗ p ( − t ) ✲✑ ✑ ❄ t = nT ✲ ✲ ❧ × ✲ a m ( t ) ✲ y p [ n ] P n ∈ Z δ ( t − nT ) ✻ d m [ n ] . . . . . . ✲ s ∗ 1 ( − t ) ✲✑ ✑ ❄ t = nT ✲ ✲ ❧ × ✲ a 1 ( t ) ✲ y 1 [ n ] P n ∈ Z δ ( t − nT ) ✻ d 1 [ n ] x ( t ) ✲ CTF ❧ ✲ x ( t ) ✻ ❄ Fig. 5. Compressed sensing of analog signals. The sampling functions s i ( t ) are obtained by combining the blocks i n Fig. 4 and are given in Theorem 2 . As a final comment, we note that we ma y a dd an inv e rtible diagonal matrix Z ( e j ω ) p rior to multiplication by A . Inde ed, in this case the measuremen ts are giv e n by y ( e j ω ) = W ( e j ω ) AZ ( e j ω ) d ( e j ω ) = A ˜ d ( e j ω ) , (49) where ˜ d ( e j ω ) has the same sparsity profile as d ( e j ω ) . Therefore, ˜ d ( e j ω ) can be recovered using the CTF block. In order to reconstruct x ( t ) , we first filter each of the non-zero seque nces ˜ d i [ n ] with the con volutional in verse of Z i ( e j ω ) . In this section we discuss ed the basic elements that allow recovery from c ompressed a nalog signals: we first use a biorthogonal s ampling set in order to access the c oefficient sequenc es, an d then employ a c on ventional CS mixing ma trix to compres s the measureme nts. Re covery is poss ible by using an IMV model and ap plying the CTF block of Fig. 3 e ither in time or in frequency . In practical a pplications we have the freedom to choos e W ( e j ω ) and A so that we e nd up with an alog sampling functions that are easy to implement. T wo examp les a re cons idered in the next s ection. V I . E X A M P L E S A. P e riodic Sparsity Suppose that we are gi ven a s ignal x ( t ) that li es in a SI s ubspa ce gen erated by a ( t ) so that x ( t ) = P n ∈ Z d [ n ] a ( t − nT ′ ) . Th e coe f ficients d [ n ] have a pe riodic s parsity pattern: Out of e ach c onsec uti ve group of m c oefficients, there are at most k non-zero values, in a given pattern. For examp le, s uppose that m = 7 , k = 2 and the sparsity profile is S = { 1 , 4 } . Then d [ n ] c an be no n-zero only at indices n = 1 + 7 ℓ or n = 4 + 7 ℓ for s ome integer ℓ . Decompo sing d [ n ] into blocks d [ n ] of length m , the sparsity pa ttern of x ( t ) implies tha t the vectors { d [ n ] } are jointly k -sparse. Since x ( t ) lies in a SI sub space spanne d by a single generator , we can sample it b y first prefiltering with the filter Q ( ω ) = 1 φ H A ( e j ωT ′ ) H ( ω ) , (50) where h ( t ) is any function su ch that φ H A ( e j ω ) d efined by (4) is non-zero a .e. on ω , and the n sampling the output at rate 1 /T ′ , as in Fig. 1. W ith this choic e, the samples c [ n ] = h q ( t − nT ′ ) , x ( t ) i are equal to the unknown coefficients d [ n ] . W e may the n us e stand ard CS techn iques to compres siv e ly sample d [ n ] . For examp le, we can sample d [ n ] se quentially by con sidering blocks d [ n ] of length m , and using a s tandard CS matrix A designe d to sample a k -sparse vec tor of length- m . Alternatively , we c an exploit the joint sp arsity by combining sev eral b locks and sampling them toge ther using MMV technique s, o r applying the IMV me thod of Section IV. Howev er , the se approach es still require an ana log samp ling rate o f 1 /T ′ . Thus, the rate red uction is only in discrete time, whereas the a nalog sampling rate remains at the Nyquist rate. Since many of the coefficients are ze ro, we would like to directly reduce the an alog rate so as not to ac quire the zero values of d [ n ] , rathe r than acquiring them first, and then compre ssing in disc rete time. T o this e nd, we note that our problem may be viewed as a special case of the general model (18) with a i ( t ) = a ( t − ( i − 1)) , 1 ≤ i ≤ m a nd d i [ n ] = d [ i + nT ] , 1 ≤ i ≤ m with T = mT ′ . Therefore, the rate can be reduced by using the tools of Section V. From Theo rem 2 , we first need to cons truct a se t of fun ctions { v ℓ ( t − nT ) } that are biorthogona l to { a ℓ ( t − nT ) } . It is e asy to s ee that v ℓ ( t ) = q ( t − ( ℓ − 1)) , 1 ≤ ℓ ≤ m , (51) with Q ( ω ) given by (50) cons titute a biorthogon al set. Inde ed, with this choice [ M V A ( e j ω )] iℓ = 1 T e j ( i − ℓ ) ω/T · X k ∈ Z e − j ( i − ℓ )2 πk /T Q ∗  ω T − 2 π T k  A  ω T − 2 π T k  = 1 T e j ( i − ℓ ) ω/T T − 1 X r =0 e − j ( i − ℓ )2 πr /T G ( e j ( ω − 2 π r ) /T ) , (52) where we defined G ( e j ω ) = 1 T X k ∈ Z Q ∗ ( ω − 2 π k ) A ( ω − 2 π k ) . (53) From (50), G ( e j ω ) = 1 . Combining this with the relation 1 T T − 1 X r =0 e − j 2 πr m/T = δ [ m ] , (54) it follows from (52) that M V A ( e j ω ) = I . W e now use Theorem 2 to c onclude that any sampling fun ctions of the form s ( ω ) = W ∗ ( e j ωT ) A ∗ v ( ω ) with V ℓ ( ω ) = Q ( ω ) e − j ( ℓ − 1) ω and Q ( ω ) given by (50) can be us ed to co mpressively sample the sign al at a rate p/T = p/ ( mT ′ ) < 1 /T ′ . In particular , giv en a matrix A of size p × m that s atisfies the CS re quirements, we ma y choose s ampling functions s i ( t ) = m X ℓ =1 A ∗ iℓ q ( t − ( ℓ − 1)) , 1 ≤ i ≤ p. (55) In this strategy , each sa mple is equal to a linear comb ination of s everal values d [ n ] , in contrast to the high-rate method in which each s ample is exactly e qual d [ n ] . As a special c ase, suppose that T ′ = 1 , a ( t ) = 1 on the interv al [0 , 1] and is zero otherwise. Thus, x ( t ) is piecewise constant over intervals of leng th 1 . Choosing h ( t ) = a ( t ) in (50) it follows that q ( t ) = 1 . T his is becaus e r aa [ n ] = h a ( t − n ) , a ( t ) i = δ [ n ] s o that φ AA ( e j ω ) = 1 . Therefore, v ℓ ( t ) = a ( t − ( ℓ − 1)) a re biorthogonal to a ℓ ( t ) . One way to a cquire the co efficients d [ n ] is to filter the s ignal x ( t ) with v ( t ) = a ( t ) and sa mple the output at rate 1 . This correspo nds to integrating x ( t ) over intervals of length on e. Since x ( t ) ha s a con stant value d [ n ] over the n th interval, the o utput will indeed be the se quence d [ n ] . T o reduc e the rate, we may instead u se the p < m sampling functions (55) and s ample the output at rate 1 /m . Th is is equiv alent to fi rst multiplying x ( t ) b y p periodic s equenc es with pe riod m . Each seq uence is piecewise con stant over intervals of length 1 with values A iℓ , 1 ≤ ℓ ≤ m . The continuo us output is the n integrated over intervals of length m to produc e the samples y ℓ [ n ] . Applying the CTF block to these measureme nts allows to recover d [ n ] . Although this sp ecial c ase is rather simple, it highlights the main idea . Furthermore, the same technique s c an b e use d even when the g enerator a ( t ) has infinite len gth. B. Multiband Sa mpling Consider next the multiband s ampling problem [21], [23] in which we h av e a complex signal that consis ts of at mos t N frequency band s, ea ch of length no larger than B . In a ddition, the signa l is b andlimited to 2 π /T . If the ba nd locations are kn own, the n we can recover such a sign al from n onuniform s amples at an average rate o f N B / (2 π ) which is typically muc h smaller than the Ny quist rate 2 π /T [18], [19 ], [20]. When the b and locations are unknown, the problem is much more comp licated. In [21] it was shown that the minimum sampling rate for such signa ls is N B /π . Furthermore, explicit algorithms were developed which achieve this rate. Here we illustrate how this problem can be formulated within our framework. Di viding the frequency interval [0 , 2 π /T ) into m sections , eac h of equal length 2 π / ( mT ) , it follows that if m ≤ 2 π / ( B T ) then each band co mprising the sign al is containe d in no more than 2 intervals. S ince there a re N bands, this implies that at most 2 N sections contain ene r gy . T o fit this problem into our gen eral model, let A i ( ω ) = √ mT A  ω − 2 π ( i − 1) mT  , 1 ≤ i ≤ m, (56) where A ( ω ) is a low-pass filter (LPF) on [0 , 2 π / ( mT )) . Th us, A i ( ω ) describes the supp ort of the i th interv a l. Since any multiband signal x ( t ) is supported in the frequency domain over at most 2 N se ctions, x ( t ) can be written as X ( ω ) = m X i =1 A i ( ω ) D i ( ω ) , (57) for some D i ( ω ) supported on the i th interv a l, where at most 2 N func tions are nonze ro. Since the s upport of D i ( ω ) has leng th 2 π / ( mT ) , it can be written as a Fourier series D i ( ω ) = X n ∈ Z d i [ n ] e − j ωnmT △ = D i ( e j ωmT ) . (58) Thus, our signal fits the general mode l (18), whe re there a re at mo st 2 N seque nces d i [ n ] that are nonze ro. W e now use our gen eral results to obtain samp ling fun ctions that can be u sed to sample an d rec over such signals at rates lower than Nyquist. On e poss ibility is to choose h i ( t ) = a i ( t ) . S ince the functions a i ( t ) a re orthonormal (as is evident by cons idering the freque ncy domain represe ntation), we have that M H A ( e j ω ) = M AA ( e j ω ) = I . Consequ ently , the res ulting sampling func tions are s i ( t ) = m X ℓ =1 A ∗ iℓ a ℓ ( t ) . (59) In the Fourier doma in, S i ( ω ) is ba ndlimited to 2 π /T and piecewise constant with v alues √ mT A ∗ iℓ over interv als of length 2 π/ ( mT ) . Alternati ve s ampling functions are thos e used in [21]: s i ( t ) = δ ( t − c i T ) , 1 ≤ i ≤ p, (60) where { c i } a re p distinct integer values in the range 1 ≤ c i ≤ m . Since x ( t ) is band limited, samp ling with the filters (60) is equiv alent to using the bandlimited functions S i ( ω ) = e − j c i ω T , 0 ≤ ω ≤ 2 π T . (61) T o show that these filters c an be obtained from o ur general framework incorporated in T heorem 2, we ne ed to choose a p × m matrix A and an in vertible p × p matrix W ( e j ω ) suc h that G i ( ω ) = S i ( ω ) where g ( ω ) = W ∗ ( e j ωmT ) A ∗ v ( ω ) , (62) and v ( ω ) represe nts a biorthogonal s et. In our setting, we can ch oose V i ( ω ) = A i ( ω ) d ue to the o rthogonality of a i ( t ) . Let A be the matrix cons isting of the ro ws c i , 1 ≤ i ≤ p of the m × m Fourier matrix A iℓ = 1 √ m e j 2 π ( ℓ − 1) c i /m , (63) and ch oose W ( e j ω ) a s a diag onal matrix with i th d iagonal element W i ( e j ω ) = 1 √ T e j c i ω /m , 0 ≤ ω < 2 π . (64) From (62), G i ( ω ) = W ∗ i ( e j ωmT ) m X ℓ =1 A ∗ iℓ A ℓ ( ω ) . (65) Since A ℓ ( ω ) is equal to √ mT over the ℓ th interv al [( ℓ − 1)2 π/ ( mT ) , ℓ 2 π /mT ) and 0 otherwise, P m ℓ =1 A ∗ iℓ A ℓ ( ω ) is p iecewise cons tant with values equal to √ mT A ∗ iℓ on interv als of length 2 π / ( mT ) . In a ddition, on the ℓ th interval, W ∗ i ( e j ωmT ) = 1 √ T e − j c i T ( ω − ( ℓ − 1)2 π / ( mT )) = ( √ m/ √ T ) e − j c i ω T A iℓ . (66) Consequ ently , o n this interval, G i ( ω ) is equal to √ mT W ∗ i ( e j ωmT ) A ∗ iℓ = e − j c i ω T . Since this express ion doe s no t depend on ℓ , G i ( ω ) = S i ( ω ) for 0 ≤ ω ≤ 2 π /T . From our ge neral results, in orde r to recover the original signal x ( t ) we need to app ly the CTF to the modified measureme nts ˜ y ( e j ω ) = W − 1 ( e j ω ) y ( e j ω ) . Sinc e W ( e j ω ) is diagonal, the D TFT of the i th seque nce ˜ y i [ n ] is gi ven by ˜ Y i ( e j ω ) = 1 √ T e − j c i ω /m Y i ( e j ω ) , 0 ≤ ω ≤ 2 π . (67) This c orrespond s to a scaled non-integer delay c i /m of y i [ n ] . Su ch a de lay can be realized by first upsamp ling the seque nce y i [ n ] by factor of m , lo w -pass filtering with a LPF with c ut-of f π /m , shifting the resulting se quenc e by c i , and then down sampling by m . This coinc ides with the approac h s uggeste d in [21] for applying the CTF directly in the time domain. Here we see that this proces sing follo ws directly from our gen eral framew ork. W e have shown that a particular ch oice of A and W ( e j ω ) results in the s ampling strategy of [21]. Alternativ e selections can lea d to a variety of diff erent sampling functions for the s ame p roblem. The add ed value in this context is that in [21] there is no discuss ion on what type of sampling me thods lea d to stable recovery . The framew ork we developed in this paper can be applied in this specific se tting to su ggest more general types of stable sa mpling and rec overy strategies. V I I . C O N C L U S I O N W e developed a ge neral frame work to treat sampling of sparse analog signals. W e focus ed on signals in a SI space generated by m kernels, where o nly k out of the m generators are active. The difficulty arises from the fact that we do no t know in advance which k are chos en. Our a pproach was based on merging ide as from stand ard an alog sampling, with results from the emerging fie ld of CS. The latter foc uses on sens ing finite-dimension al vec tors that have a sparsity s tructure in some trans form d omain. Although our problem is inherently infin ite-dimensional, we showed that by us ing the notion of b iorthogonal s ampling sets a nd the recently developed CTF block [38], [21], we ca n conv ert our p roblem to a finite-dimensional counterpa rt that ta kes on the form of an MMV , a proble m which has been treated pre viously in the CS literature. In this pa per , we focused o n s ampling u sing a bank of analog fi lters. An interesting future direction to pursue is to extend these ideas to other sa mpling architectures that may be easier to implement in ha rdware. As a final note, mos t of the lit erature to date on the exciting field o f CS ha s focuse d on sens ing of finite- dimensional vectors. On the other hand , traditional sa mpling theo ry focuses on infi nite-dimensional continuous- time signals. It is our hope that this work ca n serve a s a s tep in the direction of merging thes e two important a reas in sa mpling, leading to a more general notion o f comp ressiv e s ampling. V I I I . A C K N O W L E D G E M E N T The autho r would lik e to thank Mosh e Mishali for many fruitful discus sions, and the re viewers for us eful comments on the manus cript which helped improve the exposition. R E F E R E N C E S [1] C. E. S hannon, “Communications in the presence of noise, ” Pr oc. IRE , vol. 37, pp. 10–21, Jan 1949. [2] H. Nyquist, “Certain topics in telegraph transmission theory , ” EE T rans. , vol. 47, pp. 617–644, Jan. 1928. [3] I. Daubechies, “T he wave let transform, time-frequency localization and signal analysis, ” I EEE T rans. Inform. Theory , vol. 36, pp. 961–10 05, Sep. 1990. [4] M. Unser , “Sampling—50 years after Shannon, ” IEEE P r oc. , vol. 88, pp. 569–587, Apr . 2000. [5] A. Aldroubi and M. Unser , “Sampling procedures i n function spaces and asymptotic equi valence wit h Shannon’ s sampling theory , ” Numer . Funct. Anal. Optimiz. , vol. 15, pp. 1–21, Feb. 1994. [6] M. Unser and A. Aldroubi, “ A general sampling theory for nonideal acquisition devices, ” IEEE T rans. Signal Proce ssing , vol. 42, no. 11, pp. 2915–29 25, Nov . 1994. [7] P . P . V aidyanathan, “Generalizations of the sampling theorem: Sev en decades after Nyquist, ” IEEE Tr ans. C ir cuit Syst. I , vol. 48, no. 9, pp. 1094–1109, Sep. 2001. [8] Y . C. E ldar , “Sampling and reconstruction in arbitrary spaces and oblique dual frame vectors, ” J . F ourier Analys. Appl. , vol. 1, no. 9, pp. 77–96, Jan. 2003. [9] Y . C. Eldar and T . G. Dvorkind , “ A minimum squared-error f rame work for generalized sampling, ” IE EE T rans. Signal Pr ocessing , vol. 54, no. 6, pp. 2155–21 67, Jun. 2006. [10] T . G. Dvorkin d, Y . C. El dar , and E . Matusiak, “ Nonlinear and non-ideal sampling: Theory and methods, ” IEEE T rans. Signal Pr ocessing , vol. 56, no. 12, pp. 5874– 5890, Dec. 2008. [11] Y . C. Eldar and T . Michaeli, “Beyond bandlimited sampling: Nonlinearities, smoothness and sparsity , ” to appear in IE EE Signal Proc . Maga zine . [12] C. de Boor , R. DeV ore, and A. Ron, “The structure of fi nitely generated shift-inv ariant spaces in L 2 ( R d ) , ” J. Funct. Anal , vol. 119, no. 1, pp. 37–78 , 1994. [13] J. S. Geronimo, D. P . Hardin, and P . R. Massopu st, “Fr actal functions and wav elet expansions based on se veral scaling f unctions, ” J ournal of Appr oximation Theory , vol. 78, no. 3, pp. 373–40 1, 1994. [14] O. Christansen and Y . C. Eldar, “Oblique dual frames and shift-in v ariant spaces, ” Appl. Comp. Harm. Anal. , vol. 17, no. 1, pp. 48–68, 2004. [15] O. Christensen and Y . C. Eldar , “Generalized shift-in v ariant systems and frames for subspaces, ” J. F ourier Analys. Appl. , vol. 11, pp. 299–31 3, 2005. [16] A. Aldroubi and K. Gr ¨ ochenig, “Non-uniform sampling and reconstruction in shift-inv ariant spaces, ” Siam Revie w, vol. 43, pp. 585–62 0, 2001. [17] I. J. Schoenbe rg, Car dinal Spline I nterpolation , Philadelphia, P A: SIAM, 1973. [18] Y .-P . L in and P . P . V aidyanathan, “Periodically nonuniform sampling of bandpass signals, ” IEEE T rans. C ir cuits Syst. II , vol. 45, no. 3, pp. 340–351, Mar . 1998. [19] C. Herley and P . W . W ong, “Minimum rate sampling and reconstruction of si gnals with arbitrary frequency support, ” IEEE T rans. Inform. Theory , vol. 45, no. 5, pp. 1555–1 564, July 1999. [20] R. V enkataramani and Y . Bresler, “Perfect reconstruction formulas and bounds on aliasing error in sub-nyq uist nonuniform sampling of multiband signals, ” IEEE T rans. I nform. Theory , vol. 46, no. 6, pp. 2173– 2183, Sep. 2000. [21] M. Mishali and Y . C. E ldar , “Bl ind multi-band signal reconstruction: Compressed sensing for analog signals, ” IEEE T rans. Signal Pr ocess. , vol. 57, pp. 993–1009 , Mar . 2009. [22] M. Mishali and Y . C. Eldar, “Spectrum-blind reconstruction of multi-band signals, ” in Proc. Int. Conf. A coust., Speech , Signal Pr ocessing (ICASSP-2008), (Las V e gas, USA) , April 2008, pp. 3365–336 8. [23] M. Mishali and Y . C . Eldar, “From theory to practice: Sub-Nyquist sampling of sparse wideband analog signals, ” arXiv 0902.429 1; submitted to IEEE Selcted T opics on Signal Pr ocess. , 2009. [24] Y . M. Lu and M. N. Do, “ A theory for sampling signals from a union of subspaces, ” IEE E T rans. Signal Pr ocessing , vol. 56, no. 6, pp. 2334–2345, 2008. [25] T . Blumensath and M. E . Davies, “Sampling t heorems f or signals from the union of fi nite-dimensional linear subspaces, ” IEEE Tr ans. Inform. Theory , to appear . [26] Y . C. E ldar and M. Mishali, “Robust recovery of signals fr om a union of subspaces, ” submitted to IEEE T rans. Inform. Theory , July 2008. [27] D. L. Donoho, “Compressed sensing, ” IEEE Tr ans. on Inf. Theory , vol. 52, no. 4, pp. 1289–1306 , Apr 2006. [28] E. J. Cand ` es, J. Romberg, and T . T ao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, ” IEEE T rans. Inform. T heory , vol. 52, no. 2, pp. 489–509 , F eb . 2006. [29] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries, ” IEE E T rans. Signal P r ocessing , vol. 41, no. 12, pp. 3397–3415, Dec. 1993. [30] S. S. Chen, D. L. Donoho, and M. A. Saunders, “ A tomic decomposition by basis pursuit, ” SIAM J. Scientific C omputing , vol. 20, no. 1, pp. 33–61, 1999. [31] I. F . Gorodnitsky , J. S. George, and B. D. R ao, “Neuromagnetic source i maging with FOCUS S: A recursiv e weighted minimum norm algorithm, ” J. Electroen cephalo g. Cl inical Neur ophysiol. , vol. 95, no. 4, pp. 231–251, Oct. 1995. [32] D. L. Donoho and M. Elad, “Maximal sparsity representation via ℓ 1 minimization, ” Pr oc. Natl. Acad. Sci. , vol. 100, pp. 2197–2202 , Mar . 2003. [33] E. J. Cand ` es and T . T ao, “Decoding by l inear programming, ” IEEE T rans. Inform. Theory , vol. 51, no. 12, pp. 4203–42 15, Dec. 2005. [34] J. A. Trop p, M. B. W akin, M. F . Duarte, D. Baron, and R . G. Baraniuk, “Random filt ers for compressi ve sampling and reconstruction, ” in Pr oc. IEEE International Confer ence on Acoustics, Speec h and Signal Pr ocessing ICA SSP 2006 , May 2006, vol. 3. [35] J. N. Laska, S. Kirolos, M. F . Duarte, T . S. Ragheb, R. G. Baraniuk, and Y . Massoud , “Theory and implementation of an analog- to-information con verter using random demodulation, ” i n P r oc. IEEE International Symposium on Circu its and Systems ISCAS 2007 , May 2007, pp. 1959–1 962. [36] M. V etterli , P . Marziliano, and T . Blu, “Sampling signals with finite rate of innov ation, ” IEEE Tr ans. Signal Processing , vol. 50, pp. 1417–1 428, June 2002. [37] P .L. Dragotti, M. V etterli, and T . Blu, “Sampling moments and reconstructing signals of finite r ate of innov ation: Shannon meets Strang-Fix, ” IEEE T rans. Sig. Proc. , vol. 55, no. 5, pp. 1741–1757, May 2007. [38] M. Mishali and Y . C. Eldar, “Reduce and boost: Recovering arbitrary sets of jointly sparse vectors, ” IEE E T rans. Signal Proce ss. , vol. 56, no. 10, pp. 4692–4 702, Oct. 2008. [39] Y . C. Eldar and O. Christansen, “Characterization of oblique dual frame pairs, ” J. Applied Signal Proc essing , pp. 1–11, 2006, Ar ticle ID 92674. [40] J. Chen and X. Huo, “Theoretical results on sparse representations of multiple-measurement vectors, ” IEEE T rans. Signal Pr ocessing , vol. 54, no. 12, pp. 4634–4 643, Dec. 2006. [41] J. B. Kr uskal, “Three-way arrays: Rank and uniqueness of trili near decompositions, with application to arithmetic complexity and statistics, ” Linear Alg. Its Applic. , vol. 18, no. 2, pp. 95–138, 1977. [42] E. J. Cand ` es, “The restricted i sometry property and its implications for compressed sensing, ” C. R . Acad. Sci. P aris, Ser . I , vol. 346, pp. 589–592, 2008. [43] E. J. Cand ` es and J. Romberg, “Sparsity and i ncoherence i n compressiv e sampling, ” In verse P r ob . , vol. 23, no. 3, pp. 969–985, 2007. [44] S. F . Cotter , B. D. Rao, K. E ngan, and K. Kreutz-Delgado, “Sparse solutions to linear in verse problems with multiple measurement vectors, ” IEE E T rans. Signal Proce ssing , vol. 53, no. 7, pp. 2477–2488 , July 2005.