Multirate Synchronous Sampling of Sparse Multiband Signals

Multirate Sync hronous Sampling of Sparse Multiband Signals Mic hael Fley er, Amir Rosen thal, Alex Linden, and Moshe Horowitz Marc h 13 , 2022 The authors are with the T ec hnion—Isr a el Institute of T ec hnology , Haifa 32000 , Israel (e-mail: mikef @ tx.technion.ac.il; eeamir@tx.technion.ac.il; alin- den@ee.technion.ac.il; horowitz@ee.technion.ac.il). Abstract Recent adv ances in optical systems make them ideal for und ersampling multiband signals that ha ve h igh bandwidth s. In this pap er w e p rop ose a new sc heme for reco nstru cting multiband sparse signals using a small num b er of sampling channels. The scheme, whic h we call synchronous multira te sampling (SMRS), en tails gathering samples synchronously at few different rates whose sum is significantly low er than th e Nyq uist sam- pling rate. The signals are reconstruct ed by solving a sy stem of linear equations. W e h ave demonstrated an accurate an d robust reconstruction of signals using a small num b er of sampling channels that op erate at rel- ative ly hig h rates. Sampling at higher rates increases the signal to noise ratio in samples. The SMRS scheme enables a significant reduction in th e num b er of channels required when the sampling rate increases. W e hav e demonstrated, using only three sampling channels, an accurate sampling and reconstruct ion of 4 real signals (8 bands). The matrices that are u sed to reconstruct the signals in the SMRS scheme a lso hav e lo w cond ition num b ers. This ind icates that the SMRS scheme is robu st to noise in sig- nals. The success of the S MRS sc h eme relies on the assumption that the sampled signals are sparse. As a result most of the sampled sp ectrum ma y b e unaliased in at least one of t he sampling channels. This is in contra st t o multicoset sampling sc h emes in whic h an alias in one c h annel is equiv alent to an alias in all channels. W e hav e demonstrated that th e SMRS scheme obtains similar p erformance using 3 sampling channels and a total sampling rate 8 times t he Landau rate to a n implementation of a multicoset sampling scheme that uses 6 sampling channels with a total sampling rate of 13 times the Landau rate. 1 In tro d uction In many applications of rada rs and communications systems it is de s irable to reconstruct a m ultiband spars e signal from its samples. When the carrier fre- 1 quencies of the signal bands ar e high compared to the ov er all signal width, it is not cost effectiv e and often it is not feasible to sample at the Nyquist rate. It is therefor e desira ble to reconstr uct the sig nal fro m samples taken at ra tes low er than the Nyquist rate. Ther e is a v ast literatur e on reco nstructing signa ls from unders ampled data [2 ] − [5]. Most of metho ds are based on a m ultico set sampling scheme. In a multicoset sampling scheme m low-rate cosets are chosen out o f L cosets of samples o btained from time uniformly distributed sa mples taken at a rate F where F is greater or equal to the Nyquist rate F ny q [3]. In each channel the sampling is offset b y a different predetermined integer mu ltiple of the recipro cal of the rate F . The da ta from the different sampling channels a re then used to reconstruct a signal b y solving a system of linea r equations. Un der certain conditions on the sampling rate and num b er of channels, a prop er choice of the time o ffsets ensur e s that the equations ha ve a unique solutio n in c a se tha t the signal bands are known a prio r i [3], or unknown a prior i [5]. A prev ious w or k has demo nstrated a different sch eme for reconstructing sparse m ultiband signa ls [9]. The scheme, called m ultira te sampling (MRS), ent a ils gathering s amples at P different rates. The num b er P can be s mall and do es not dep end on a ny characteristics of a signa l. The approa ch is not int ended to o btain the minimum sampling r ate. Rather, it is intended to re- construct signals accurately with a very high proba bility a t an ov er all sa mpling rate s ignificantly low er than the Nyquis t r ate under the co nstraint of a small nu mber of channels. The reconstructio n metho d in [9] do es not require synchronization of the different sampling channels. In a ddition to r educing the complexity of the sam- pling hardware, unsync hr onized sampling relaxes the stringent r equirement in m ultico set sampling schemes of a v ery small timing jitter in the s ampling time of the channels. Sim ulatio ns hav e indicated that MRS reconstruction is robust bo th to different s ignal t y p es and to rela tively high no ise. Accurate signal recons truction using the unsy nchronized MRS scheme of [9] requires that ea ch frequency in the s uppo rt of the signal b e unaliased in a t least one of the sampling c hannels. In this pap er w e describ e a new reconstruc tio n algorithm that ov er comes this deficiency by using sy nchronized sampling chan- nels. In our synchronized multirate sampling (SMRS) scheme, the alias ing is resolved, in the same spirit as m ulticos et sampling sc hemes , b y solving a system of linear equations. Ther efore, the SMRS scheme enable s the reconstruction of signals that ca nnot b e r econstructed using MRS scheme. In our SMRS scheme each sampling rate must b e an integer multiple of the same basic frequency . The reconstr uction metho d in our SMRS scheme requires the same frequency resolution in all sampling channels despite the fact that the sampling rate is different in e a ch channel. This requirement is w ell suited to the implementation of the s cheme using th e optical system descr ib ed in [10]. The s ampling is per formed in t wo steps. In the fir st step the entire signal spectrum is do wnconverted into a low frequency region called base ba nd by m ultiplying the signal with a train of short optical puls es [10]. In each of the s ampling channels the pulse rate is different. The frequency downcon verted 2 analog sig nals a re then conv erted in each sampling ch a nnel to digital signals using an A/D electr onic co nv erter that sa mples at the hig hest of the channel rates. The result is that the signal is sampled in a ll channels with the same time r esolution that is deter mined b y the sampling rate of the A/D co nv erter. Alternatively , since each sampling rate must b e an integer m ultiple o f the same basic frequency , a common frequency res o lution can be o btained b y using the same time window for all channels. There is an inherent adv antage to sampling, in each channel, nea r the max- im um sampling rate ma de po ssible by cost and tec hnolo gy . This is beca use sampling at higher r ates increa ses the signal to no ise ratio in sa mpled signals [9]. Our simulation re sults indica te that when the sampling rate in each channel increases, our SMRS scheme requir es s ignificantly fewer sa mpling channels than do es a m ultico set sampling scheme o f [5 ] to obtain comparable reconstruction success. When the s ampling r ate in each channel increas es, the probability that a sparse signa l a liases simultaneously in all sampling channels beco mes very low in an MRS scheme. It is lower than in a multicoset sa mpling scheme in which, bec ause all channels sample at the sa me frequency , an a lia s in one channel is equiv alent to an a lias in all channels. Our n umerical sim ulations indica te that the success r ate of our SMRS scheme is significantly higher than that of a m ul- ticoset sampling schemes o f [5] when the num b er o f sa mpling channels is small (3 in our simulations) a nd the sampling ra te of of each channel is high. Exactly the same data obtained in the SMRS scheme ca n be o btained using a m ulticoset sa mpling sch eme. How ever the multicoset pattern requires many more channels, each o f which samples at a very low ra te. In a n umer ic al example of section 4.1 a 3 c ha nnel multirate sampling pattern is equiv a lent to a 5 8 channel m ultico set sampling pattern. With this m ulticoset pattern the time differe nc e betw ee n t wo consecutive samples is on the order of 1 psec . Such an accuracy cannot b e pr actically a chiev ed. In an SMRS scheme the data is reconstr ucted differently than in a multicoset scheme implementation of [5]. In this multicoset recovery scheme it is ass umed that, in addition to b eing spars e, the spectrum of the signal co nsists of ba nds each of which is narr ow er than the coset s ampling r ate. How ever, o ur SMRS scheme r equires no such a s sumptions o n the or iginating signa l. The recons tr uc- tion o f spa rse sig nals in our SMRS scheme is robust. This is b ecause most of the sampled sp ectrum is unaliased in at leas t o ne of the s ampling c ha nnels. This is in contrast to the equiv alent multicoset sampling scheme with many low rate sampling channels in which the a lias proba bility is very high a nd an alias in one channel is equiv alent to an a lias in all channels. In s uch cases a blind signal recov er y of [5] can b e found only by us ing a pursuit algorithm whereas , in man y cases, the SMRS sc heme doesn’t require a pursuit algor ithm. The signal is re- constructed direc tly by a single matrix inv ersion. By making a very reas onable ph y s ical assumption, we are also able to simplify the reconstruction by r e ducing the n umber of p oss ible sig nal lo cations in a stra ightforw ar d manner . The SMRS s cheme reconstructs signa ls by so lving a sy s tem o f lin ea r equa- tions. Our simulations indicate that linear equa tions in an SMRS scheme are nu mer ically stable in that they ha ve low condition n umbers. This ma kes them 3 rather insensitive to nois e. In the SMRS scheme, when sampling at total sam- pling r ate that is significantly higher than theor etically r equired, the probability that part of the signal spectrum will no t b e aliased in at lea st one o f the sam- pling channels is very high. The unaliased par ts of the sig nal can b e reco vered directly from the sampling channels. O ur sim ulatio n results indicate that the sensitivity of the reconstructed signal to noise added at frequencies of ba seband that are not a liased in at lea st one of sampling channels is very low. 2 Sync hronized MRS In this section w e describe the SMRS scheme. Let F max be an assumed maxi- m um ca rrier frequency and let X ( f ) ∈ L 2 ([0 , F max ]) be the F ourier transfor m of a complex-v alued signa l x ( t ) that is to be reconstruc ted from its samples. Throughout the a nalysis we normalize the F ourie r tra nsform by X ( f ) = Z ∞ −∞ x ( t ) e − j 2 π f t dt. The mo difications r equired to reconstr uct real-v alued s ignals ar e descr ibe d in the App endix. W e a ssume that the signa l x ( t ) to b e sampled, in addition to being bandlimited in [0 , F max ], is multiband; i.e., the supp or t of its s p ec tr um is con tained within a finite disjoin t union of in terv als ( a n , b n ], ea ch of whic h is con taine d in [0 , F max ]. By assumption, max b n ≤ F max . In reco nstructing a signal we do not as sume any a prior i k nowledge of the num b er o r lo cation of the int er v a ls ( a n , b n ]. The only requirement is that the signal be spar se; i.e., that for a signal whose spectral supp or t is co ntained within the N interv a ls ( a i , b i ], P N k =1 b k − a k ≪ F max . In the MRS scheme of [9] a signal is s a mpled at P different sampling r ates F i ( i = 1 . . . P ). If the dela ys of the c hannels are denoted b y ∆ i , the sampled signals x i ( t ) a re given by x i ( t ) = x ( t ) ∞ X n = −∞ δ  t − n F i − ∆ i  (1) where δ ( t ) is a dirac delta “function”. The spectrum X i ( f ) o f the sampled signal in the i th channel sa tisfies X i ( f ) = F i ∞ X n = −∞ X ( f + nF i ) exp[ j 2 π ( f + n F i )∆ i ] . (2) Since the channels ar e synchronized in time, we may set each ∆ i = 0. E quation (2) b ecomes X i ( f ) = F i ∞ X n = −∞ X ( f + nF i ) . (3) 4 It follows from (3) that all the information ab out the i th sampled spec- trum X i ( f ) is contained in the interv al [0 , F i ]. T his interv al is ca lled the i th baseband. T o process the sampled signals it is necessary to dis cretize the fre- quency a xis. W e use the sa me freq uency r esolution (∆ f ) in all of the sampling channels. The sa me frequency reso lution is dir e c tly obtained if the sy stem is implemen ted using an optical sampling system [10]. T o use our reconstruction metho d each sampling ra te should b e chosen to b e an integer m ultiple M i of the basic frequency res olution: F i = M i ∆ f . By using a sa mpling time window with a dur ation T = 1 / ∆ f , the sa me frequency resolution ∆ f is obtained in all of the sampling c hannels . Alternatively , the sa me frequency resolution can be obtained in an optical sy stem b y using A/D con verters with the same rate at th e output of a ll the optical downcon vertors [10]. T o reduce co mputational requirements we ignore r edundant data . Thus, o nly M i ent r ies of the Discrete F ourier tra ns form (DFT) are retained for further p r o cessing. W e repr esent the signals over t he co mmon discretized frequency a xis with the following notations: X i [ k ] = X i ( k ∆ f ) , k = 0 , . . . , M i − 1 , X [ k ] = X ( k ∆ f ) , k = 0 , . . . , M − 1 , where ∆ f is the frequency resolution, M = ⌈ F max / ∆ f ⌉ , X i [ k ] is the sp e c trum of the sampled data fro m the discretized frequencies in the ba seband [0 , F i ] and X [ k ] is the spectr um of the originating signal at the discretized frequencies . Equation (3) ta kes the for m X i [ k ] = F i M − 1 X l =0 X [ l ] ∞ X n = −∞ δ [ l − ( k + nM i )] (4) where δ [ n ] is the Kronecker delta function. Equation (4) can b e w r itten in matrix form a s follows: x i = A i x (5) where x i and x ar e g iven by ( x i ) k +1 = X i [ k ] , 0 ≤ k ≤ M i − 1 , (6) ( x ) k +1 = X [ k ] 0 ≤ k ≤ M − 1 , and A i is a M i × M matrix whose elements are g iven by A i k +1 ,l +1 = F i ∞ X n = −∞ δ [ l − ( k + nM i )] . (7) Each element A i k,l is equal to either F i or 0. This is b eca use there is at most one c o ntribution in the infinite sum o f δ ’s which is ma de when l ≡ k ( mo d M i ). In each r ow of the matrix A i there a re only ⌊ F max /F i ⌋ non ze r o elements. F or each v alue o f i ( i = 1 . . . P ), (5 ) defines a set of linea r equa tions that relate the spe c tr um of the sig nal to the sp ectrum of its sa mples. The vector 5 x in (5) is the sa me for all the P equations b ecause it do esn’t dep end on the sampling. The vector b x and the ma tr ix b A are obtained b y concatenating the v ector s x i and matrices A i as fo llows: b x =      x 1 x 2 . . . x P      , b A =      A 1 A 2 . . . A P      . These for m the sys tem of equa tions b x = b A x . (8 ) The m a trix b A ha s exactly P non-v anishing elemen ts in eac h column that cor- resp ond to the lo c a tions of the sp ectra l r eplica in e ach channel baseband. In case that the signal is r eal-v alued its sp ectr um fulfills X ( f ) = X ( − f ) (9) where a + bj = a − bj is the complex conjugate. The equations fo r reco ns tructing such a sig na l a r e describ ed in the App endix. T o inv er t (8) and calculate the discretized signal sp ectr um ( x ) it is necessa ry that the num b er of rows P P i =1 M i in A b e equal to or lar ger than the num b er of columns M . Defining F total = P P i =1 F i makes this condition equiv alent to the condition F total > F max . (10) The condition o n the sampling ra tes given in (10) is consistent with the r equire- men t that the sampling rate b e g reater than the Nyquist ra te of a g eneral signal whose sp ectra l supp or t is [0 , F max ]. How ever, when sa mpling sparse signals, a n inv ersion o f the matrix may b e p oss ible even if the c ondition (10) is no t fulfilled. Our o b jective is to in vert (8) in the case of spar se s ignals with sampling rates F total < F max . 3 In v ersion A lgorithm In this s ection we descr ibe our inv ersion algo rithm for the SMRS scheme. The purp ose of the algo r ithm is to inv ert (8); i.e., to ca lc ulate the vector x fro m the vector b x . T o in vert the equations with sampling rates low er than those prescrib ed by (10) the as sumption that the sig nal is s parse sho uld b e taken into account. 3.1 Kno wn B and Lo cations In the case in whic h the s ig nal band lo cations ( a n , b n ] a re known (8) can b e simplified easily . All the elemen ts of x that cor resp ond to the frequencies not in 6 the sp ectra l supp ort ∪ n ( a n , b n ] ar e eliminated from (8). All the co lumns of the matrix which corr esp ond to these elements ar e also eliminated. The re duced system o f equations that co rresp onds to (8) is giv en by b x red = b A red x red . (11) A unique so lution exis ts o nly if b A red is full co lumn rank. In this case the inv erse ca n be found using the Mo ore-Penrose pseudo-inverse [12]. In a ma trix of full column rank the num be r of rows must equa l or ex ceed the num b er of columns. Although the entire sp ec tr um is downconv erted to ba seband, we assume that we ar e sampling hig hly spa rse signals. Hence the num b er of non-zer o entries in each x i is significantly s maller than M i in at least one o f the sa mpling channels. The num b er of non-zer o entries in each x i might even b e s maller b ec a use of aliasing. A necessar y condition for a unique inv er se or pseudo-inv er s e o f b A red is that this num b er still b e gr eater than the num b er of non-zero entries of x . This is consistent with a La ndau theor em [1] that s tates tha t one ca nnot reconstruct a signal if the sp ectral density o f the s amples collected fro m all sa mpling channels is less than the sp ectr a l supp o r t of the o r iginating s ignal. The choice o f sampling rates impo ses restrictions on the po s sible v alues of F max for which a n in version of (11) is p os sible. F or the matrix b A red to hav e full column ra nk, it m ust not ha ve any identical columns. Since w e do no t restrict the possible lo cations of the kno wn signal bands, any co mbination of columns of the matrix b A may app ea r in the matrix b A red . Therefore we requir e that b A not hav e an y iden tical c olumns. The matr ix b A is compos ed of P sub-matrices A i whose columns a re per io dic: A i k,l + M i = A i k,l . F or the ma trix b A no t to b e p erio dic, it is required that an y common p erio d of the P sub-ma trices b e larger than M . This condition is met if the leas t common multiple o f the { M i } i is lar ger than M . As a res ult, F max should fulfill F max < l c m ( M 1 , . . . , M P )∆ f , where l cm denotes least common multiple. 3.2 Unkno wn Bands’ Location In the case that the lo catio ns of the bands ( a n , b n ] are no t known a prio ri some additional as s umptions m ust b e ma de . In the m ultico set recovery scheme o f [5] it was assumed that the maximum band width is k nown a prio ri. In o ur SMRS scheme w e do not make any as s umptions on the interv a ls ( a n , b n ] but instead we add ass umptions on a sig nal’s sp ectrum itself. W e assume that, fo r ea ch discretized frequency k ∆ f , any sampled sp ectrum X i ( k ∆ f ) = 0 is due only to lac k of a signal in any o f its replica s and not due to any aliasing. In other words if, f o r any n , X ( f + nF i ) 6 = 0, then X i ( f ) 6 = 0. Another assumption is that there is a unique solution in the case the signal suppo rt is known, i.e., the matrix b A red has a full column rank for a known suppo rt. 7 Applying the first assumption, o ne can detect baseband frequencies in which there is no signal. These frequencies can b e eliminated in the reduction o f (8). W e desc rib e this s imple pro cedure for eliminating frequencies which, a ccording to our assumption, ca nnot b e part of the sp ectral suppor t o f the orig inating signal. The elimination is s imilar to o ne presented in asynchronous-MRS [9]. W e deno te the indicator function χ i [ l ] as follows: χ i [ l ] =  1 , for all l ∈ [0 , M − 1] such that X i ( l ∆ f ) 6 = 0 0 , otherwise . (12) The function X i ( f ) is pe r io dic with perio d F i . Therefor e χ i [ l ] is a p erio dic extension of an indica tor function ov er the baseband f ∈ [0 , F i ). W e define the χ [ l ] as follows: χ [ l ] = P Y i =1 χ i [ l ] , l ∈ [0 , M − 1 ] . (13) The function χ [ l ] equals 1 o ver the intersection of all the upconverted bands of the P sampled signals and it defines the co lumns of the matrix b A that are retained in forming the reduce d matrix b A red . All o ther columns are eliminated and their co rresp o nding elements in the vectors x are a lso eliminated. After the elimination o f the columns from the matrix b A , the matrix rows which corresp ond to zero elements in b x and their cor resp onding elements in the vectors b x a re a lso eliminated. In some cases the function χ [ l ] equals 1 only for frequencies within the spec tral suppo rt of the signal. In such cases the resulting equations are ident ica l to those found in the previous subs e c tion (equation (1 1)). How ever, in o ther cas es, χ [ l ] ma y also equal 1 for frequencies o utside the signal’s true sp ectral supp ort. In s uch cas es the r educed matrix will hav e mo re columns than the matrix in the cas e in which the sp ectral supp ort of the signa l is k nown. As a r esult the inv er sion req uires finding the v alues of more v ariables. Each eliminated zer o ener gy baseba nd comp onent causes elimination of re- sp ective rows and columns. The elimina tion of o ne bas eband entry mea ns that all the frequencies that are downcon verted to that base ba nd en try (the aliasing frequencies) a re a lso eliminated. This is b ecause of our first as sumption that zero entry in the baseband cor resp onds to ze ro entries in all o f the frequency comp onents of the o riginal signal that are do wn-co nv erted to frequency of t he baseband en try . Therefore, elimination of one baseband entry r esults in elimi- nation o f ⌊ F max / min { F i }⌋ to ⌈ F max / max { F i }⌉ co rresp o nding columns. Th us, if the n umber of the zero elements in b x is sufficiently lar ge, the n umber of r ows in the matrix b A red may be la rger than the num b er of columns. If in addition, matrix b A red has a full column rank, the pro blem is either consistent o r ov er determined. In such cases there is a unique inv ers ion to (1 1) which can be found using the Mo ore- Penrose pseudo- inv erse. If the ma trix is not full column rank, the problem is underdeter mined and the inv ersion is not unique. A unique solution in such cases can b e obtained either by increasing the total sa mpling r ate o r b y adding additional ass umptions on the signal. 8 3.3 Ill-p osed cases In ma ny ca ses the matrix in (11) for unknown band loc ations is not full column rank. In these ca ses there are subsets o f columns in the matrix b A red that are linearly dep endent. Using this linear depe ndenc e , a solution to (11) ca n be found. Howev er any solution found can b e used to c onstruct a n infinite num b er of solutio ns to the equation. Th us, there is no unique in versse to (11) and the inv ersion pr oblem is ill-p os e d. T o reco nstruct a signal in the case in which the inv er sion problem is ill- po sed we impos e an additional assumption on the signal. W e assume that in the case the signal supp ort is unknown and the problem is ill-p osed, among all po ssible solutions t o (11), the o riginating s ignal is the one tha t is composed of the minim um num b er of bands. This is the signal we attempt to construct. W e also as s umed e arlier that the r educed matrix which co r resp onds to the c a se in which the signal s upp o rt is known is well-p osed. In case that this assumption is not fulfilled the sampled signal do e s not contain e nough information for solving the problem. Under the thr e e assumptions stated above (the assumption that leads to matrix reduction, the existence of the unique solution t o (11) when the sig nal bands are known, and band-s parsity) the inversion pro blem is reduced to finding the solution of (11) that is comp os ed of the minimum n umber of bands . The problem is NP- complex since we need to tes t every p os sible com binatio n o f bands. The algor ithm describ ed here is o f low er complexity and its purp ose is to find a solution of (11) that is comp osed of the minimum num b er of bands without testing all the com bina tio ns. The resulting algor ithm a ttains a lower success rate but decreases the run-time significantly as compar ed to an NP- complex algorithm. W e do not pro vide the co nditio ns under which the corr ect solutio n is obta ined. Our a lgorithm is ba sed o n the O rth og onal M atching P ursuit ( OM P ) [8]. This algor ithm belo ngs to the category of the ” Greedy Search” alg orithms. The original OMP algo rithm is used to find the sparses t solution x of underdeter- mined equations Ax = b [8] where A is an underdetermined matrix. The sparsest solution is the so lution having th e s mallest norm k x k 0 where k x k 0 is the n umber of non ze r o elements in the vector x . The original OMP algor ithm collects columns of the ma trix A iteratively to construct a reduced ma tr ix A r . A t e a ch iteration n the column o f A whic h is added to A n − 1 r to pro duce a matrix A n r is the column which r esults in the smallest res idua l erro r min x k b − A n r x k 2 2 where for every vector y , k y k 2 2 = P i y 2 i . The iterations a re stopp ed whe n some threshold ǫ is achiev ed. Sufficient conditions ar e giv en for the algor ithm to obtain the co rrect so lution [8 ]. W e denote A = b A red , b = b x red and x = x red . Since w e are seek ing the solution of Ax = b with the smallest n umber of bands and not the sma llest nor m k x k 0 , we modify the OMP algo rithm by instead o f choosing a sing le column as in [8] by selec ting iteratively blocks of possible lo c a tions. The columns of the matrix A fall into J blo cks. Ea ch block is indexed by j and B j contains the 9 index of co lumns of the j th blo ck. Each ind ex set B j ident ifies a p ossible band of the spectr al suppor t of the reconstructed signal. W e star t the iteratio n with the empty set S 0 = ∅ of column indexes, the empty matrix A 0 r , a nd the set B 0 = S J j =1 B j , so that at n th iteratio n the following holds: S n S B n = B 0 . A t th e n th iteration ( n > 1 ) the algo r ithm must decide whic h blo ck to add to A n − 1 r . If the index set B j is c hose n, then S n = S n − 1 ∪ B j and B n = B n − 1 \ B j . The matrix A n r is the matrix whose columns a re selected fro m A acc ording to the index ing se t S n . The blo ck a dded is the one that pro duces the smallest residual error ǫ n = min B j ∈ B n − 1 min x k b − A n r,j x k 2 2 where A n r,j is the matrix obtained by a dding the blo ck indexed by B j to A n − 1 r . The algorithm stops when the threshold ǫ is reached. The threshold ǫ is a very small num b er and reflects up on the finite nu mer ical precis ion o f the c a lculations. The algorithm p erfor med well in our s im ula tions. How ever, there w er e cases in which the support of the reco nstructed s o lution did no t co ntain all the orig- inating bands and cases in which the reconstructed signa l was incor r ect even though all the assumptions on a signals given in sec tion 3.2 were fulfilled. Only after p erforming the last step of the algo rithm is it po ssible to determine the suppo rt of the signal. The a lg orithm failed pr imarily for one of t wo re a sons. One of them was due to the inclusion of a blo ck that reduced the residual erro r on one hand but on the other hand, caused a res ulting matrix A n r to be not full co lumn r ank as hypothesized in our problem (in sectio n 3.2 ). This c an hap- pen, for exa mple, when a blo ck consists o f a cor rect sub-blo ck and erroneo us sub-blo cks. Including any e r roneous sub-blo cks may result in an ill-po sed prob- lem. Another reas o n for failure was a la r ge dynamic range o f the s ignals. When reconstructing suc h signals, co rrect bands may b e ignored by the algorithm in cases that the ener gy within the bands is sig nifica ntly lower than the e ne r gy in other bands. Sufficien t conditions that assure that the alg orithm co nv erges to a unique solution hav e not yet b een determined. 4 Sim ulations results The abilit y of the signal reconstruction alg orithm to reco ver different t yp es of signals was tes ted. In o ne set of sim ulations the abilit y o f the alg orithm to reconstruct multiband complex and rea l-v alued signa ls with different s p e c tral suppo rts, shap es, a nd band widths that w er e not known a pr iori was tested. Additional sim ulations were performed in whic h rea l- v a lued m ultiband signals were co nt a minated by additiv e white noise. Band ca rrier frequencies were chosen from a unifor m distribution ov er the maximum supp o r t: 0 -20 GHz for complex signal and − 2 0 to 20 GHz for real-v alued signals. F or each set o f simulations we counted the mean rate o f ill-p os ed cases in which the mo dified O MP algo rithm was us e d to recover the signa l. Mean times for accurate signal reconstruction were also reco rded. F ailures of the reconstruction were either b ecause o ne of the initial assumptions was not fulfilled or b eca use of the failure o f the mo dified 10 OMP algorithm. In different simulations the width of each band (hence the total bandwidth of the signal) w as v aried. The n umber of bands w a s alwa ys set equal to 4 for complex signals and to 8 for rea l-v alued s ignals (4 p ositive ba nds and 4 negative bands). Ho wev er , the recons truction a lg orithm was unaw are o f this n umber . In all the simulations the frequency resolution was set to 5 MHz. The simulations were per formed on a 2 GHz Cor e2Duo CPU with 2 GB RAM storage in the MA TLAB 7.0 environment (no sp ecial pr ogra mming was p er formed to use b oth cores). 4.1 Ideal m ultiband signals In the case of ideal multiband signals , b ecause of the absence o f ener g y outside of stric tly defined bands, one ca n exp ect a p er fect reconstruction. Accordingly , the a lgorithm w as ev aluated b y a p erfect reconstruction criterion; i.e., a mea n difference betw een the true and the r e c onstructed s ignal s p ec tr um les s than 10 − 10 . Whenever this error was attained, the re construction was deemed to hav e been successful. Otherwis e , it was deemed to hav e failed. The threshold for the mo dified OMP was chosen accordingly ; ǫ = 10 − 20 . Sim ula tions w er e carried o ut for co mplex signals to c ompare the results to those published using the multicoset sampling r ecov ery scheme of [5]. The sampling ra tes w er e chosen to b e 0.9 5, 1 .0 and 1 .05 GHz yielding a total sampling rate F total = 3 . 0 GHz. Different signals w ith 4 ba nds of equal width w er e gener ated. E ach band was chosen to lie within the interv al [0 , 20] GHz. Bo th the real and imaginary sp ectra of the sig na l within each band w ere chosen to b e no rmally distributed. Spec ific a lly , for each freq uency f = k ∆ f in a chosen band, the real and imag- inary c omp onents of X ( f ) were chosen ra ndo mly and indep endent ly fro m a standard nor mal distribution. E ach ba nds’ sp ectra were scaled by a constant α such that ea ch bands’ e nergy was equal to a uniformly generated v a lue E on the in ter v a l [1 , 5 ]; i.e., for sp ecific band, X ( f ) = α [ X r ( f ) + j X im ( f )] , k X ( f ) k 2 = E . These signals were also used to test the mult ico set sampling r econstruction scheme o f [5]. The empirical succes s ra tes were obtained from 100 0 runs, each with a different total bandwidth ( F Landau ). The s uc c e ss rate is shown in Fig. 1. W e ha ve v alidated that the e mpir ical success rate did not significantly change when the num b er of simulation runs was increas e d from 100 0 to 5000 . As is evident from Fig. 1, the empirical success p er centage o f an ideal recon- struction is high when F total /F Landau ≥ 5 . In the SBR2 scheme (downsampling factor L = 199 ) in [5] the e mpir ical suc c e ss rate shows p erfect r e construction for at lea st p = 14 channels. This corresp o nds to F total = p / LT = 1 . 4 GHz and hence F total /F Landau should b e greater than about 3. Although the total sam- pling ra te in [5] is low er tha n in SMRS, the num b er of channels that are used in that scheme is significantly higher compar e d to that used in SMRS wher e o nly 3 channels a re used. 11 Another simulation in [5] shows that fo r low er num ber of c hannels with L = 23 empirical p e r fect reconstruction is achiev ed with at lea st 6 channels and F total /F Landau > 13. In the SMRS scheme empirica l p erfect reconstr uction was obtained using only three channels with a total sampling rate F total /F Landau ≥ 5 . The s ystem par ameters (num ber of sa mpling channels, sampling rates, F max ) that were used in our last simulation ar e the same as those used in our optical sampling exp erimental set up. The fact that the simulation results were obtained in a practical situatio n demonstra tes that our SMRS scheme ca n reconstruc t sparse signa ls per fectly using b oth a fewer n umber of sa mpling channels and with a lower to tal sa mpling r ate tha n are r e q uired by multicoset sampling schemes. W e note that the same data that is obtained in a SMRS pattern can a lways be obtained by a mult ico set sampling pattern since the ratio b etw een each pair of sampling rates is rational. In our example the sampling rate of ea ch coset is equal to 1 /LT = 50 MHz. The n umber of multicoset sampling channels ( p ) is 58 . The time o ffset b etw een the cosets is a multiple of T = 1 399 GH z . The downsampling factor L is 399 GHz / 50 MHz=7980 . Note that since L is not prime, we ar e not g uaranteed to obta in a universal sa mpling patter n [5]. Fig. 1 shows that in o ur scheme 100 p er cent empirica l success was obtained for F total /F Landau > 5. This corresp onds to a maximu m ov erla p in a coset scheme of K = 16 a nd hence p / K = 3 . 6 for p = 58. In compariso n, for p = 1 4 channels and downsampling fac tor of L = 199 in a m ulticoset sampling scheme in [5], th e emp ir ical success rate of 10 0 per cent was obtained for SBR2 algorithm for p/ N = 3 . 5. This ratio is simila r to the equiv alent ratio obtaine d in our scheme (3.6) although the downsampling factor in our equiv alent multicoset scheme is equa l to 7980 . F or appr oximately 95 pe r cent success rate in SB R2 scheme of [5], the ratio p/ N is 2.75. In our scheme, for the same empirical success rate of 95 p erc e nt, F total /F Landau = 3 . 53 and the equiv a lent ra tio of the nu mber of channels to maximal num b er of ov erla ps is p/K = 2 . 4. The mean per centage of ill- p o sed case s is shown in Fig . 2. The figure shows that fo r 4 ≤ F total /F Landau ≤ 10, in most of the t es ted c a ses the in version was ill p osed. Nonetheless, a very high success p ercentage was obtained for these v a lues. This indicates tha t our mo dified OMP algor ithm was very successful in resolving these ca ses. Fig. 3 shows the mean run time as a function of F total /F Landau (constant total sampling rate and v a rying sig nal supp ort). B e c ause matrix inv er sion is the most computatio nally intensiv e op er ation in the a lgorithm, mean r un time decreases with a reduction in signal bandwidth. This is b ecause , with a fixed resolution, matrix size is prop or tional to signa l bandwidth. Moreover, as the ratio F total /F Landau increases, the p o ssible spe c tr al suppo rt obtained at the first step of the reconstructio n increases beyond the increase of the s ig nal bandwidth. W e note that we could reduce the run time without significa nt ly a ffecting the empirical success p er centage b y solving (8) using a low resolution and re- constructing the sig nal us ing a higher resolutio n. The algorithm, mo dified a s explained in the Appendix, w as also tested against rea l-v alued signals. The assumed maximu m frequency F max was se t to 20 GH z. The nu mber of sampling channels w as set to 3 with the sampling 12 3 4 5 6 7 8 9 10 99.4 99.6 99.8 100 F total /F Landau Success pct. Figure 1: E mpirical success p e rcentages for 4- bands co mplex signa ls for differe nt sp ectral suppor ts ( F Landau ) with F Nyquist =20 GHz, total sampling r ate F total = 3 GHz 4 5 6 7 8 9 10 0 20 40 60 80 100 F total /F Landau Ill−conditioned pct. Figure 2: I ll- p o sed cases m ea n percentage fo r 4- bands co mplex signals fo r dif- ferent spectra l supp or ts ( F Landau ) with F Nyquist =20 GHz, to tal sa mpling ra te F total = 3 GHz frequencies c hosen to be F 1 = 3 . 8 GHz, F 2 = 4 . 0 GHz, and F 3 = 4 . 2 GHz re- sulting in a tota l sampling rate, F total = 12 GHz. The sampling frequencies are the same as ar e used in o ur exp eriments ba sed on asy nchronous MRS [9]. The nu mber of bands was set to 8 (4 p ositive and 4 negative frequencies, assuming no carrier frequency so low a s to have the 0 frequency in its s p ec tr um). Each band was chosen to b e of equa l width F Landau / 8. Once a band ( a, b ] was chosen, the sp ectr um o f X ( f ) for f ∈ ( a, b ] w a s determined by the fo llowing fo rmula: X ( f ) = A sin  π ( f − a ) b − a  e j θ . The pha se θ was chosen randomly fro m a uniform distributio n o n [0 , 2 π ] and 13 4 5 6 7 8 9 10 0 0.5 1 1.5 F total /F Landau Time (sec) Figure 3: Mean run times for 4-ba nds complex s ignals for differen t sp ectral suppo rts ( F Landau ) with F Nyquist =20 GHz, tota l sampling rate F total = 3 GHz the amplitude A was chosen ra ndomly from a uniform distribution on [1 , 1 . 2]. Fig. 4 shows the e mpir ical s uccess per centages of the algorithm tested aga inst real v alued signa ls. As is ev ident from the fig ur e, the empir ic al success perc e nt - age is high when F total /F Landau ≥ 8 . W e note that the required s ampling rate is significantly higher in this example than in the complex signa ls simulation. The reason is that in our real ca se ex ample there are twice a s many bands as in the complex case simulations. Hence, after the sampling, an overlap may also o ccur bet ween the negative and the p ositive bands of the r eal signal. W e note that when sampling a r eal s ignal at a sampling r ate F i , it is sufficient to k now the sp e ctrum in a frequency region [0 , F i / 2]. How ever, for re al s ig nals, there is uncertaint y as to whether a sig nal in ba seband is obtained fro m a signal in the po sitive band o r in the nega tive band. The num be r of ill-p osed cases and the mean recovery r un times for the real- v a lued sig nals a re shown in Figs. 5 a nd 6 resp ectively . It can b ee seen that the mean r ate of ill-conditioned cases is muc h low er for r eal-v alued signal simulations than for complex ones. This could be due to the corr elation betw een po sitive and neg ative frequency comp onents of r e al s ig nals. 4.2 Solution Stabilit y The stability of the linear equatio ns used in the r ecov ery s cheme was tes ted via the co ndition num b er fo r the real-v alued signals equations ((21) in the Ap- pendix ). T his c ase is imp o rtant in our exp eriments since we sample r eal s ig- nals [9 ]. The re c onstruction scheme for rea l- v a lued sig nals requires solving tw o systems of linea r equa tio ns; one for the re a l part a nd the other for the ima ginary part. Each of the t wo systems of equations is describ ed by a different matr ix. In each test case we pres ented the maximum v alue of the condition n umbers of the tw o matrices . F or 4-band 20 0-MHz-width r andomly generated signals ( F Landau =1.6 GHz) the conditio n num b er a mong 100 0 r uns was a t most 5 .3. The conditio n num b er s histo gram is shown in Fig. 7. 14 4 5 6 7 8 9 10 95 98 100 F total /F Landau Success pct. Figure 4: Empirical success p ercentages for equal 4-bands real signals for dif- ferent spectra l supp or ts ( F Landau ) with F Nyquist =40 GHz, to tal sa mpling ra te F total = 12 GHz 4 5 6 7 8 9 10 0 20 40 60 80 100 F total /F Landau Ill−conditioned pct. Figure 5: Ill-po sed cases mea n rate for 4-bands r e al signals for different sp ectr al suppo rts ( F Landau ) with F Nyquist =40 GHz, total sa mpling r a te F total = 12 GHz When s ampling a s parse signal, most freq uencies of the sig nal are unaliased in at le ast one o f the sa mpling channels. It can b e easily shown that when a frequency comp onent of the signal is not aliased in a ny of the sampling c han- nels, the rec o nstructed signa l is obtained simply by averaging the c o rresp o nding sampled spectrum at the different sampling channels. The reconstruction of a frequency comp onent that is unaliased in a t leas t one of the sampling channels can also b e easily p erfor med by co pying the cor resp onding unaliased s p ec tr um to the reco nstructed signal. The r efore, for spar se s ignals, the re c o nstruction in the SMRS s cheme is ro bus t and the co nditio n num b er of the matrices is small. Indeed, we have verified that the condition n umber o f the matrix inc r eases as the nu mber of frequencies in the original signal that are alia sed increases. Fig. 8 shows the mean v alues of condition num b ers versus num b er of aliased frequen- cies as calculated for the 4 real sig nals each with a 200 MHz width that are 15 sampled at F total = 12 GHz. T he other pa r ameters of the simulation are the same a s tho se in the simulation that r esulted in Fig. 4. 4.3 Noisy signals The alg o rithm’s per formance was also tested for its abilit y to r econstruct real- v a lued s ignals contaminated by Gaussia n white noise. The presence of no ise demands so me mo dification of the algorithm. One m o dification is in detecting the p ossible ba nds of the supp or t of the originating signal. Bec ause the s pe c tr al suppo rt of white noise is not re s tricted to the sp ectral s uppo rt o f the uncontam- inated signal, the indicator functions in (12) c annot be used. Instead, we a dapt (12) to noisy cas es s imila rly as in [9]. In [9], fo r the indicator function χ i [ l ] to be equal to 1 at any frequency , it w as required that the av era ge energ y of the signal in the neigh b o r ho o d of that frequency be higher than a cer tain threshold. In SMRS we further expa nd each band in χ [ l ] to include additional frequencies that might o therwise be omitted when defining the indicator functions χ i [ l ]. Once the bands are ident ified the matrix equations are co nstructed exactly as in the noiseless ca se. The so lutio n of the linear equations given in (21) is mo dified in the noisy case . Because the added white no ise affects the entire sp ectr um, a signal co ntaminated by white noise can no longer b e considere d m ultiband in the s trict se ns e. Th us one cannot ex pe ct to reco ns truct it p erfectly from samples taken at a total rate low er than the Nyquist rate. Whereas in the ideal noiseless ca se the er ror norm v anishes, with a signal containing noise, one must relent on a perfect reconstructio n and settle for a minim um erro r. In the noisy cas e the s olution to (21) sho uld solve the least square problem min x r,im red k b x r,im red − b A r,im red x r,im red k . When the the matrix b A r,im red is not full column rank, we use the modified OMP algorithm which is adjusted to account for the er rors due to noise. As noted ab ov e, in the no is eless case, o ne can exp ect a p erfect reconstr uction and th us the thre s hold er ror ǫ ca n b e set to 0 or a very small num b er. How ever, with noisy signa ls, some car e must be taken in c ho o sing ǫ . On the one hand, if the ǫ is chosen to o la rge, the a lgorithm may stop befor e a solution is reached. On the other ha nd, if ǫ is c ho s en to o small, the rec o nstructed s ignal may include ba nds that a re not in the origina ting signa l. The pro blem o f to o high threshold is solved by c ha nging the s top criter ion. Instead of stopping the algo rithm when a threshold is attained, we chec k at each iteration whether the blo ck tha t reduces the residual er ror the most ca uses the resulting matr ix to b e rank deficient. When this o ccurs , the iter ation a re s topp ed and the block is not a dded to th e matrix. An additional change is made to the algor ithm when trea ting the blo cks. In the noisele s s case e a ch block corr esp onds to a single ba nd in χ [ l ]. When sampling noisy signa ls, we divide each blo ck into several sub-blo cks. The reason for this division is that, with noisy sig nals, the iden tificatio n of the bands is not accur ate. Ident ified bands may b e wider than the or iginating bands. This is par ticularly true if the thresho ld is chosen small. T his widening may cause the inclusio n of 16 4 5 6 7 8 9 10 0 2 4 6 8 F total /F Landa u Time (sec) Figure 6: Mea n recov er y times for equal 4-bands rea l signals fo r different sp ectral suppo rts (total bandwidth) with F Nyquist =40 GHz, total sampling rate F total = 12 GHz 1 2 3 4 5 0 10 20 30 40 Condition numbers Pct. of runs Figure 7: Condition n umber s histo gram o f the mixing matrix in (21) for 4 ba nds, 200 MHz width rea l sig nals, F total = 1 2 GHz false freq uencies whose corr esp onding columns in b A r,im red are linearly dep endent on the columns c o rresp o nding to the s uppo rt of the o riginating signal. B y using smaller sub-blocks such columns may b e isolated from the rest o f the columns in their co rresp onding band. The recovery scheme w a s tested against real-v alued signals with 8 bands (4 p ositive frequencies ba nds and 4 negative frequencies ba nds ). The signals without noise were genera ted and sampled exactly as in the noiseless simulations of rea l sig na ls. Noise was added r andomly a t each frequency of the pr e -sampled signal according to a nor mal distribution with sta ndard deviation σ = 0 . 0 4; the SNR w as defined by 10 log 10 (1 / ( σ p F max /F 2 )) = 10 . 5 dB. This definition takes int o ac c ount the accumulation of noise in baseband due to sampling. The sampling rates were the same as those in the noisele ss s imu la tions. The indicator functions χ i [ l ] w er e constr ucted using the same par ameters as those used in [9]. 17 40 80 120 1 2 3 4 5 Number of aliased frequencies Avg. condition no. Figure 8: Condition num b er s mean v alue vs . num b er of aliased fr equencies for 4 ba nds, 20 0 MHz width r eal signa ls, F total = 12 GHz 4 5 6 7 8 9 10 88 92 96 100 F total /F Landau Success pct. Figure 9: E mpirical success pe r centages for 4 bands of rea l signals noise that were co ntaminated with a noise with a standa rd deviation of σ = 0 . 04 fo r different sp ectral supp orts of the sig nal ( F Landau ) with F Nyquist =40 GHz, and a total s a mpling r ate F total = 12 GHz Each band in χ [ l ] w as widened by 20 p ercent on each side. The sub-blo cks used in the mo dified OMP had sp ectra l width of 100 MHz. The succe ss was measur ed by the algorithm’s abilit y to a chiev e a low err or l 1 norm be low 2 σ p F max /F 2 = 4 . 47 σ for each re cov ered band. The mean error for each recovered band X rec ( f ) and the true ba nd X ( f ) were ev aluated as follows: 1 | B | Z B | X rec ( f ) − X ( f ) | d f < 2 σ p F max /F 2 where B is the band supp or t. Statistics on recov ering 8 bands 200 MHz width each are ba sed on 1 0000 tests. The sim ulatio n s how ed that, altho ugh the algo rithm’s p erfo rmance in- evitably decrea s ed, it still achiev ed a high empiric a l recov er y r ate (37 failures 18 out of 10000 tests). Additional simulations w er e p erfor med by changing the Landau rates as w a s done in the sim ulatio ns perfor med for the noiseless case. In Fig. 9 the empirical success percentage is pre sented for 1000 sim ulations of noisy s ignals. The results of the simulations ar e similar to those in the noiseless case. When the total s ampling r ate is 8 times higher than the Landa u rate, high success p ercentage w a s achieved. The re c ov ery er ror level dep ended o n the threshold choice. Lower threshold allows mo re accurate recons truction but increases the recov ery t ime. Moreov er , additiona l parameters adjustmen ts are also necessa ry (widening p ercentage, sub-blo cks size). Differen t e rror criter ia are also po ssible. F or example, cho osing l 2 norm instead of l 1 norm and setting the err o r thresho ld to b e 3 . 3 σ as in [9 ] resulted in 9 9 . 5 p ercent empirical success rate in r ecov ering 1.6 GHz La ndau ra te signa ls and 99 . 8 p er cent for 1 .5 GHz Landau rate s ig nals. This is in co ntrast to sim ulatio ns results in Fig. 9 where empirical perfect rec o nstruction w as obtained for those signals und er different criteria. 5 Conclusions In this pap e r w e describ e a synchronized m ultirate sampling scheme for a ccu- rately recons tr ucting spar se m ultiband signals using a small num b er o f sampling channels (3 in our sim ulations) whose total s ampling r ate is significan tly low er than the Nyquist r a te. Although th e same data that is o btained fro m our scheme can b e obtained by a multicoset scheme, such a multicoset scheme requires many mor e channels and a time accura cy that ca nnot b e attained pr actically . Mo reov er , our scheme pro cesses the da ta differently , in a w ay that results in significan tly wider b a se- bands and thus greatly r educes the effects of alia sing. By synch r onizing the sampling channels our sc heme is able to reconstruct signa ls correctly even in most ca ses in which the sig nal is alia s ed in all the sampling channels. The main a dv a ntage of multicoset s ampling schemes is theoretica l. Becaus e in a m ulticoset sa mpling scheme ea ch channel samples at the same frequency , one has a mathematical structure that enables a p erfect re constructions of ideal m ultiband signa ls from samples taken at a tota l sampling rate that is close to a theoretical low er b ound. This bo und is attaine d o nly under sp ecial a ssumptions regar ding the num b er and width of the signa l bands. Mo r eov er, the bound requires the nu mber o f sampling c hannels to b e twice the num b er o f signal bands [5 ]. Hence, in many cases the num ber o f s a mpling channels b eco mes to o high for a practica l implementation. The main adv a ntage of the SMRS is in the use of a small num b er of sa mpling channels that o p er ate a t relatively high ra tes that can reco ns truct accurately sparse signal that consist of several bands. Sampling at higher rates has a fundamen ta l adv antage in that it increases the SNR after sa mpling. Another adv ant a ge is that implemen tation of the SMRS scheme do es not r equire a pr iori knowledge of the maximum width of the s ignal bands. A third adv antage o f our scheme is that, in many cases the signal ca n b e r econstructed by s imple matrix 19 inv ersion rather t ha n thr ough a search a lgorithm as in the m ulticoset re c ov ery scheme of [5]. In the SMRS scheme, when sampling s parse signals, most of the sampled sp ectrum is unaliased in at least one of the sampling channels. Hence, mo st of the spe ctrum can b e recons tr ucted directly from the unaliased parts of the sp ectrum. On the other hand, in multicoset sampling sc hemes an alias in one channel is equiv alent to an alias in all channels. F urthermore, a multicoset scheme downconv erts signals to mu ch low er frequencies, thus increasing the negative effects of alia sing. Our n umer ic al s imulations indicate that the reduced aliasing in o ur SMRS sc heme results in significantly be tter performance ov er a m ultico set sa mpling scheme of [5 ] whos e num b er of sampling channels is small. Also, due to r educed aliasing, it is exp ected that the reconstruction in the SMRS scheme will be r obust when noise is added in the sampling pro ces s. Although we do not hav e a rigo rous criteria for p erfect r e construction, by examining the multicoset patter n that yields the same da ta it migh t be p ossible to obta in necess ary co nditions for a p erfect reconstr uction. Our s cheme also so lves ill-conditioned linear equations b y using a modified OMP algorithm. Whereas we obtained satisfactor y results with it, we do not hav e criter ia for deter mining when the modified OMP algorithm converges to the correct solution. Ho wever, this shortcoming may b e not as significa nt b ecaus e of the p oss ibilit y of using other alg orithms to so lve the equa tions. 6 App endix In this app endix we present the mo dificatio ns to (8) for the real signals r ecov ery . Since the sig na l is r eal-v alued, its sp ectr um fulfills X ( f ) = X ( − f ) (14) where a + bj = a − bj is the c omplex conjugate and a and b ar e r eal nu mber s. It follows fro m (14) and (1), that for ea ch channel index i , a ll the informatio n ab out X i ( f ) is contained in the interv al [0 , F i / 2]. Consequently , it is co nvenien t to choo se the sa mpling frequencies F i such that F i / 2 = ∆ f M i / 2 where M i / 2 is an int eg er. Because the conjugation o p eration a + j b : a + j b 7→ a − j b is not complex linear , (5) needs to b e replaced with tw o sy s tems of equatio ns ; one for the real par t and o ne for the imagina r y part. W e use the following notations to represe nt the s pe ctrum o f the real signals in the discretized frequencies: X i [ k ] = X i ( k ∆ f ) k = −⌊ M i / 2 ⌋ , . . . , ⌊ M i / 2 ⌋ , (15) X [ k ] = X ( k ∆ f ) k = −⌊ M / 2 ⌋ , . . . , ⌊ M / 2 ⌋ . The seque nc e X i [ k ] c o ntains the samples of X i ( f ) in the baseband [ − F i / 2 , F i / 2]. The seq uenc e X [ k ] contains the sa mples of X ( f ) g iven in [ − M ∆ f / 2 , M ∆ f / 2], where M is chosen to fulfill M = ⌈ F Nyquist / ∆ f ⌉ . Equa tio n ( 3 ) now tak es the 20 following form: X i [ k ] = F i ⌊ M / 2 ⌋ X l = −⌊ M / 2 ⌋ X [ l ] ∞ X n = −∞ δ [ l − ( k + nM i )] . (16) Equation (16 ) can b e wr itten in a ma trix fo rm a s x i = A i x (17) where x i and x ar e g iven by ( x i ) k + ⌊ M i / 2 ⌋ +1 = X i [ k ] , − ⌊ M i / 2 ⌋ ≤ k ≤ ⌊ M i / 2 ⌋ , (18) ( x ) k + ⌊ M / 2 ⌋ +1 = X [ k ] , − ⌊ M / 2 ⌋ ≤ k ≤ ⌊ M / 2 ⌋ , and A i is a matrix whose elements a r e given by A i k + ⌊ M i / 2 ⌋ +1 ,l + ⌊ M / 2 ⌋ +1 = F i ∞ X n = −∞ δ [ l − ( k + nM i )] . (19) Note that, since the signal is r eal v alued, all o f its spectra l informa tion is con- tained in the p ositive freq uenc ie s. Each element in A i is equal to either F i or 0. Eq uation (17) for the different sampling channels can b e c o ncatenated as in complex s ignals ca se to yield b x = b Ax . (20) The sp ectrum can b e decomp osed in to its real and imaginary parts. As a result (20) b eco mes b x r = b A r x r , (21) b x im = b A im x im where b x r = Re( b x ) and b x im = Im ( b x ). In addition only comp one nts that cor- resp ond to p ositive freq uencies ar e re ta ined in the vectors b x r and b x im . The elements of the ma trices b A r and b A im are given by b A r k, ⌊ M / 2 ⌋− l +1 = b A k,l +1 + b A k,M − l , l = 0 , . . . , ⌊ M / 2 ⌋ , (22) b A im k, ⌊ M / 2 ⌋− l +1 = b A k,M − l − b A k,l +1 , l = 0 , . . . , ⌊ M / 2 ⌋ . The re c o nstruction is p er formed with (2 1) exactly as in the co mplex cas e. References [1] H. Landau, ”Necess ary densit y conditions for sampling and interp olation of certain entire functions,” Acta Math. , vol. 117, pp. 37 -52, July 1 967. 21 [2] A. Ko hlenberg, ”Exact Interpola tion of Band-limited F unctions,” A ppl. Phys. , vol. 2 4,pp. 1 432- 1436, 1953. [3] R. V enk antaramani and Y. Bres ler, ” Optimal s ub-Nyquist nonunif o rm sa m- pling and reconstr uction for multiband signals,” IEEE T r ans. S ignal Pr o- c ess. , vol. 4 9, pp. 2 301-2 313, O ct. 2001 . [4] Y. M. Lu and M. N. Do, ” A Theory for Sampling Signals from a union of Sub-spaces,” IEEE T r ans. Signal Pr o c ess. , to be published. [5] M. Misha li and Y. Elda r, ”Blind multiband s ignal recostr uction: compress e d sensing for analog s ignals,” arXiv:070 9.1563 , Sept. 2007. [6] Y. P . Lin and P . P . V aidyanathan, ”Perio dically uniform sampling of band- pass signals,” IEEE T r ans. Cir cuits Sys. , vol. 4 5, pp. 3 40-3 51, Mar. 1 9 98. [7] C. Herley and W. W ong, ”Minimum rate sampling and reconstruction of signals with a r bitrary frequency supp or t,” IEEE T r ans. Inform. The ory , vol. 45, pp. 15 55-15 64, Jul. 1999. [8] A.M. Br uckstein, D.L. Donoho and M. Elad, ”F rom Sparse Solutions of Systems of Equations to Sparse Mo deling of Signals and Images ,” SI AM , to be p ublishe d. [9] A. Rosenthal, A. Linden and M. Horowitz, ”Multirate asynchronous sam- pling o f multiband s ignals,” s ubmitted for publication. [10] A. Zeitouny , A. F eldser, and M. Hor owitz, “Optica l sampling o f narr owband microw av e signa ls using pulses generated by electroabs orption mo dulato r s,” Opt. Comm. , v ol. 25 6, pp. 2 48-25 5, Dec. 2 005. [11] I. Stewart a nd D. T all, The F oun dations of Mathematics , O xford, Engla nd: Oxford University Press, 1 977. [12] R. Penrose,”A genera lized inv erse for matrices,” in Pr o c. Cambridge Philo- sophic al So ciety , Ca mbridge, vol. 51, 1 9 55, pp. 40 6-41 3. [13] R. V enk atarama ni a nd Y. Br esler, ”Sampling theorems fo r uniform and per io dic no nuniform MIMO sa mpling of multiband signals ,” IEEE T r ans. Signal Pr o c ess. , vol. 51 , pp. 3152 -316 3, Dec. 20 03. [14] E. J. Cand‘es, J. Romberg and T. T ao, ”Robust uncertaint y pr inciples: exact signal reco nstruction from highly incomplete fre quency information,” IEEE T r ans. Inform. The ory , vol. 5 2, pp. 4 89-50 9. 22