Multi-rate synchronous optical undersampling of several bandwidth-limited signals

We demonstrate experimentally an optical system for undersampling several bandwidth-limited signals with carrier frequencies that are not known apriori and can be located within a broad frequency region of 0-20 GHz. The system is based on undersampling synchronously at three different rates. The optical undersampling down-converts the entire system bandwidth into a low frequency region called baseband. The synchronous sampling at several rates enables to accurately reconstruct signals even in cases in which different signals overlap in the baseband region of all sampling channels. Reconstruction of three simultaneously generated chirped signals, each with a bandwidth of about 200 MHz, was experimentally demonstrated. © 2010 Optical Society of America OCIS codes: (000.0000) General. References and links 1. P. W. Joudawlkis, J. J. Hargreaves, R. D. Younger, G. W. Titi, J. C. Twichell, “Optical Down-Sampling of WideBand Microwave Signals,” J. Lightwave Technol. 21, 3116−3124 (2004). 2. A. Zeitouny, A. Feldser and M. Horowitz, “Optical sampling of narrowband microwave signals using pulses generated by electroabsorption modulators,” Opt. Commun. 256, 248−255 (2005). 3. J. Kim, M. J. Park, M. H. Perrott and F. X. Kärtner, Photonic subsampling analog-to-digital conversion of microwave signals at 40-GHz with higher than 7-ENOB resolution, Opt. Express 16, 16509−16515 (2008). 4. P. Feng and Y. Bresler, “Spectrum-blind minimum-rate sampling and reconstruction of multiband signals,” in Proc. IEEE Int. Conf. ASSP, vol. 3, pp. 1688−1691, May 1996. 5. R. Venkantaramani and Y. Bresler, ”Optimal sub-Nyquist nonuniform sampling and reconstruction for multiband signals,” IEEE Trans. Signal Process. 49, 2301−2313 (2001). 6. M. Mishali and Y. Eldar, ”Blind multiband signal reconstruction: compressed sensing for analog signals,” IEEE Trans. Signal Process. 57, 993−1009 (2009). 7. A. Feldster, Y. P. Shapira, M. Horowitz, A. Rosenthal, S. Zach and L. Singer, “Optical Under-Sampling and Reconstruction of Several Bandwidth-Limited Signals,” J. Lightwave Technol. 27, 1027−1033 (2009). 8. M. Fleyer, A. Linden, M. Horowitz and A. Rosenthal, ”Multirate Synchronous Sampling of Sparse Multiband Signals,” IEEE Trans. Signal Process. 58, 1144−1156 (2010). 9. J. F. Gravel and J. Wight, “On the Conception and Analysis of a 12-GHz PushPush Phase-Locked DRO,” IEEE Trans. Microwave Theory Tech. 54, 153−159 (2006). 10. G. C. Valley, “Photonic analog-to-digital converters,” Opt. Express 15, 1955−1982 (2007). 11. R. H. Walden, “Analog-to-digital converter survey and analysis,” IEEE J. Sel. Areas Commun. 17, 539−550 (1999). 12. B. Le, T. W. Rondeau, J. H. Reed and C. W. Bosti, “Analog-to-digital converters,” Signal Process. Mag. 69, 69−77 (2005). 13. M. Shinagawa, Y. Akazawa and T. Wakimoto, Jitter Analysis of High-Speed Sampling Systems, IEEE J. SolidState Circuits 25, 220−224 (1990). #127514 $15.00 USD Received 3 May 2010; revised 1 Jul 2010; accepted 1 Jul 2010; published 26 Jul 2010 (C) 2010 OSA 2 August 2010 / Vol. 18, No. 16 / OPTICS EXPRESS 16929 14. M. Rodwell, D. Bloom and K. Weingarten, ”Subpicosecond Laser Timing Stabilization,” IEEE J. Quantum Electron. 25, 817−827 (1989). 15. R. Penrose,”A generalized inverse for matrices,” in Proc. Cambridge Philosophical Society, Cambridge 51, 406−413 (1955).


Introduction
Passive electronic warfare (EW) systems include reconnaissance or surveillance equipment that detects and analyzes electromagnetic radiation from radar and communications transmitters in a potential enemy's aircraft, missiles, ships, satellites and ground installations.Since the bandwidth of the radio frequency (RF) spectrum that is currently used for military applications is extremely broad, on the order of tens of GHz, one of the main limitations of current passive EW systems is the analog to digital (A/D) signal conversion.Analog to digital conversion can be performed by multiplying a signal by a train of short pulses.Optics is advantageous for such applications since optical pulse sources can generate a train of high-repetition rate short pulses with a very low jitter as required for sampling signals with a very high carrier frequency.
Although in modern military situation the available frequency spectrum is very broad, at a given time only a small portion of the spectrum is occupied.Therefore, the sampling can be performed below the Nyquist rate (maximum RF signal frequency).Such technique is called undersampling and it has been demonstrated in optical systems [1]- [3].Optical undersampling is obtained by multiplying the RF signal by a train of short optical pulses.This operation downconverts the entire broad spectrum of the original signal into a low frequency region called baseband.In such sampling scheme different parts of the signal spectrum may overlap at the baseband frequency region.This effect, called aliasing deteriorates the baseband spectrum and it prevents the reconstruction of the sampled signal.There is a vast literature on theoretical methods for signals reconstruction from the undersampled data [4]- [6].
In a previous work we have demonstrated theoretically and experimentally a new optical system for undersampling several RF signals [7].This scheme entails gathering samples at several rates with a total sampling rate that is significantly lower than the Nyquist sampling rate.The main advantage of the asynchronous sampling scheme is that it does not require the knowledge of the time offsets between the sampling channels.Hence, the required hardware becomes simple and robust to errors in sampling times and to differences between the sampling channels.However, an accurate reconstruction of the signal spectrum by the asynchronous multi-rate sampling scheme requires that each frequency of the signal support is unaliased in at least one of the sampling channels baseband.
In [8] we have demonstrated theoretically a new sampling scheme that we call Synchronous Multirate Sampling (SMRS).We have shown that synchronization between the sampling channels enables to resolve aliasing and overlaps and to accurately reconstruct signals that could not be reconstructed by the asynchronous multirate sampling scheme [7].The reconstruction of signals in the SMRS scheme is based on solving a system of linear equations that connects the Fourier transforms of the original and baseband signals.Simulation results show that the empirical reconstruction success rate in SMRS scheme is much higher than obtained in asynchronous sampling scheme.For example, the reconstruction success rate of a 4 real-valued signals with a bandwidth of 100 MHz each and a carrier frequency located in the region [0,20] GHz equals about 60% using the asynchronous sampling scheme [7] and more than 98% by using the SMRS scheme with the total sampling rate of 12 GHz.
In this paper we demonstrate experimentally an optical system that implements the SMRS scheme.The system is based on sampling synchronously at three sampling rates.In our experiments we have accurately reconstructed the spectrum of three signals that were generated simultaneously.The signal carrier frequencies of the signals were unknown apriori and the bandwidth of each signal was 100−200 MHz.The bandwidth of the sampling system was about 20 GHz.We have tested about 100 scenarios for different signals carrier frequencies and bandwidths.The SMRS scheme enabled to experimentally reconstruct all the scenarios that were checked.The reconstruction was obtained even for cases in which some down-converted signals overlapped with other signals in all three sampling channels.Such signals can not be reconstructed by the asynchronous sampling scheme.
In order to use the SMRS scheme it should be calibrated and its parameters should be extracted and applied to the signal processing algorithm.We have developed a parametric model of our sampling system and a method to extract the required parameters as described in Section 4 of this manuscript.
The performance of our sampling system is limited by the signal to noise-plus-distortion ratio (SNDR).The SNDR of a signal at a carrier frequency of 11.1 GHz measured with a resolution bandwidth of 100 kHz (that corresponds to a sampling window of 10 µs) equals to about 57 dB.This SNDR value corresponds to a resolution of 9.2 effective number of bits (ENOB).
The rest of the paper is organized as follows: in section 2 we describe the experimental setup, in section 3 we briefly describe the SMRS reconstruction algorithm, in section 4 we describe the method for extraction of the system parameters and in section 5 we describe the experimental results.

Synchronous multirate sampling scheme
Optical undersampling is obtained by modulating a short optical pulse-train by an RF signal [1]− [3].The spectral width of the optical pulses should be sufficiently large to enable the sampling of signals with a high carrier frequency.Undersampling enables to down-convert a broad frequency region into a low frequency region called baseband.The down-converted signal in the baseband region can then be sampled using a conventional electronic analog to digital (A/D) converter.In order to reconstruct signals with unknown carrier frequency (blind reconstruction) we undersample the RF signal using three optical sampling channels that operate with different parameters.
The setup of our system is shown on Fig. 1.The sampling is performed by modulating the amplitudes of three optical pulse-trains with rates F 1 = 3.8 GHz, F 2 = 3.5 GHz and F 3 = 4.0 GHz by an RF signal.The rates F i were chosen to fulfill F i = M i ∆ f , where M i are positive integers and ∆ f = 100 MHz is the bandwidth of a frequency cell in the reconstruction algorithm [8].
The optical pulse-trains were generated by combining optical and electronic devices.Electrical pulses were generated by using a comb generator that was connected to a phase-locked loop dielectric resonator oscillators (PLL-DROs) [7].The electrical pulses at the output of each comb generator were supplied to electro-absorption (EA) modulator that modulated the intensity of a continuous wave (CW) laser.The PLL-DRO enabled generating low-jitter pulses while the EA modulator enabled to shorten the pulse duration from about 50 ps to about 24 ps as required in our system.To compensate losses we used an erbium doped fiber amplifier (EDFA) with a 30 dB gain and electrical amplifiers (AMP) with an electrical gain of about 25 dB.The optical powers and losses in the system are the same as described in our setup for asynchronous sampling [7].
The high frequency signal of a PLL-DRO resonator is stabilized by locking it to a high harmonic of a reference signal via a balanced mixer and a phase detector [9].To synchronize three optical pulse trains we connected three PLL-DROs to the same reference signal.Therefore, three optical pulse trains were locked and synchronized to the same reference signal with a frequency of 100 MHz.The synchronization of three PLL-DROs to a single reference source caused the overlap between the optical pulses of three optical pulse trains every 1/100 MHz = (a) Multirate sampling scheme (b) Optical pulses generator (OPG) Fig. 1.Schematic description of (a) optical sampling system and (b) optical pulse generator (OPG).The synchronous multirate sampling is performed simultaneously at three different rates: F 1 = 3.8 GHz, F 2 = 3.5 GHz and F 3 = 4.0 GHz.A 100 MHz reference signal is used to synchronize between the different channels.OPG consists of a phase-locked loop dielectric resonator oscillators (PLL-DROs), a comb generator, an Electro-Absorption (EA) modulator, and a continuous wave (CW) laser.The optical wavelengths of the OPGs were 1535.04 nm, 1536.61nm and 1544.53 nm.MUX and DEMUX are optical add-drop modules that combine and demultiplex three optical pulse trains.EDFA is an erbium-doped fiber amplifier, MODULATOR is a LiNbO 3 Mach-Zehnder modulator, ϕ is an electrical phase shifter, D is a detector, LPF is an electrical lowpass filter with a bandwidth of 2 GHz and AMP is an electrical amplifier.10 ns.
The optical wavelengths of three optical pulse sources were chosen to be different as in [7].Therefore, the optical pulse trains could be combined without adding a significant loss by using an arrayed-waveguide grating (MUX).Such method is used in optical communication systems that are based on wavelength division multiplexing (WDM) technique.The input RF signal modulated the intensity of the three pulse trains by using a single-drive Mach-Zehnder electrooptic modulator.The V π of the modulator equals 4.6 V and the maximum power of the RF signal at the modulator input was about 0 dBm.The bias voltage of the modulator was adjusted to obtain the minimum second order distortion.Thus, we operated the modulator in its linear operating region.
The signal at the output of the modulator was demultiplexed by another arrayed-waveguide grating (DEMUX).Three optical signals at the output of the optical demultiplexer were detected and electronically sampled by a 3 channel sampling scope that operated at a sampling rate of 5 Gsamples/sec at each channel.The electronic analog to digital conversion was performed simultaneously in all the sampling channels.A fourth A/D sampling channel was used to sample the 100 MHz reference signal.This sampling channel enables to estimate the delay between the electronic A/Ds and the pulse trains as described in the sequel.

Synchronous MRS signal reconstruction algorithm
We briefly describe our signal reconstruction method.A detailed description of the algorithm is given in [8].Our assumption is that the sampled signal is sparse; i.e, it occupies only a small fraction of the frequency bandwidth of the system.Our system bandwidth equals to 20 GHz and the total bandwidth of the signals should not exceed about 2 GHz.
The original RF signal x(t) modulates simultaneously P optical pulse-trains with frequencies where M i are integers.The resulting signals x i (t) contain replicas of the original signal spectrum X( f ) shifted by an integer multiples of F i [2].The baseband frequency region of the i-th channel equals [−F i /2, F i /2].These baseband signals are electronically sampled for each channel synchronously at a rate F s ≥ max i {F i }.
The original signal spectrum X( f ) is connected to the sampled signals spectra X i ( f ) through the system of linear equations [8]: where the vector x (β ) represents the Fourier transform of the sampled signals, the vector x (β ) corresponds to the original signal Fourier transform and Q is a sampling matrix.Vectors x (β ) and x (β ) are given with a resolution ∆ f .The variable 0 ≤ β < ∆ f gives the location within the frequency cells that construct the spectrum vectors.Each column of the matrix Q corresponds to a frequency cell of the original signal spectrum and each row corresponds to a frequency cell of the down-converted signal.In undersampling the number of rows in the matrix Q is smaller than the number of columns.Therefore, a unique reconstruction of the original signal can not be obtained unless more assumptions on the signal are added.Although in modern military environment the available frequency spectrum is heavily used, at a given instant of time only a small portion of the spectrum is occupied.Therefore, in electronic warfare systems it is required to sample and reconstruct multi-band sparse signals.
In the first step of our reconstruction algorithm, possible signals locations are identified.In noisy environment a threshold is applied to identify the baseband signal locations [7].Frequency cells with values above the threshold are identified as a baseband signals.Then, the frequency locations of the original signal that do not contribute to a baseband signals are detected and eliminated from the sampling matrix Q.If the resulting reduced matrix has a full column rank, then the search is complete and the signal can be reconstructed via the pseudo-inverse of the sampling matrix.If the reduced sampling matrix is not full column rank, we assume that the signal obtained after the reduction step is block-sparse.Therefore, we look for a solution that fulfills the reduced system of equations, but also contains the smallest number of bands.The pursuit algorithm that we use to find such solution is a variation of the well-known Orthogonal Matching Pursuit (OMP) modified to search for blocks instead of columns.Numerical simulations indicate that a very high success rate that is close to 100% is obtained when the sum of the sampling rates in all channels is approximately 8 times higher than the total signal frequency support.This success rate is higher than obtained by using previously published methods for undersampling assuming that the number of sampling channels is limited to 3 [8].

Non ideal synchronous multi-rate sampling
In a non-ideal system the implementation of the SMRS reconstruction scheme requires to perform some modifications on the algorithm given in [8].Particularly, in a non-ideal system there is a need to extract parameters of each of the sampling channels.Our experimental results indicate that this extraction should be performed only once over a long time period of few days.
In our system, the RF signal x(t) is undersampled by multiplying the signal with a low-rate pulse train.The resulting signal is then detected, amplified and sampled by an electronic A/D converter.The signal at the input of the i-th (i = 1..3) electronic A/D, x i (t) is given by where * denotes a convolution operator, p i (t), τ i and h i (t) are the pulse shape of the i-th channel, the delay of the i-th pulse train at t = 0 and the impulse response function of the detector and the electrical amplifier of the i-th channel, respectively.Although the electrical signals at the output of the PLL-DROs are synchronized, a time delay τ i must be added due to a slightly different propagation times of the signals in different channels.The time delay τ is added since the electronic A/D starts sampling in a time offset with respect to the pulse trains.In deriving Eq. ( 2) we neglect the transfer function of the modulator.This function is required to obtain the correct spectrum amplitude at each frequency; however, since the same modulator is used to modulate all three pulse trains, it does not affect the stability of the reconstruction.In a non ideal system, Eq. ( 2) can be written in the frequency domain, as shown in Appendix A (Eq. ( 13)): where X( f ), X i ( f ), H i ( f ) and P i ( f ) are the Fourier transform of the original signal x(t), the baseband signal x i (t), the transfer function h i (t) and a pulse profile p i (t), respectively.Equations ( 2) and ( 3) indicate that in order to accurately reconstruct the sampled signal the system parameters τ i , τ, p i (t) and h i (t) should be known.These parameters change the sampling matrix that was given in [8] for an ideal sampling.In Appendix A we develop equation ( 22) that connects vectors of the sampled baseband signals x (β ) and the original signal spectrum x (β ): where the vector x (β ) and the matrix Q are obtained by concatenating vectors W H i (β )x i (β ) and matrices Q i W τ i respectively: where P i ( f ) is the Fourier transform of the i-th channel pulse, and matrices W τ i and W H i (β ) equal: Notice that matrices Q i correspond to a single channel sampling matrix and the matrix Q is the sampling matrix for three sampling channels used in our system.

Measurement of the impulse response h i (t)
The transfer function h i (t) of the i-th sampling channel is determined by the delay of the optical channels between the modulator and the detectors, the detector response function and the response function of the electrical amplifiers and filters.We measured the Fourier transform of h i (t) − H i ( f ) by supplying a sinusoidal wave with a controllable frequency to the sampling system and to additional (fourth) port of the electrical A/D.The frequency of the sinusoidal wave varied between 5 MHz and 2 GHz with a step size of 5 MHz.The transfer function H i ( f ) was measured for each channel i (i = 1 . . .3). Figure 2 shows the amplitude and the phase of H i ( f ) for three channels.The same results were also obtained by using a network analyzer.

Minimizing the relative delays between τ i
In our scheme all three optical pulsed sources are locked to a common reference signal with a frequency of 100 MHz.However, since different optical pulse trains propagate through slightly different fiber lengths before arriving to the modulator, the pulse trains are shifted by a delay τ i .Delays of different optical pulse trains τ i can be made approximately equal by adding a tunable electrical phase shifter (ϕ) to each pulse generator as shown in Fig. 1.The tuning of the phase shifters is performed by measuring the pulse trains at the entrance of the modulator by using a sampling oscilloscope with a 50 GHz bandwidth optical input.Figure 3 shows three optical pulse trains when the phase shifters are adjusted to obtain the best overlap between them.We have verified that the relative delays between the pulse trains did not change over a very long time duration of several days.Therefore, after the calibration process we could assume in our reconstruction algorithm that τ i = 0. Fig. 3. Three optical pulse-trains measured at the modulator input by using a sampling oscilloscope with an optical bandwidth of 50 GHz.Delays between the pulse-trains were adjusted to obtain the best overlap between pulses.

Extracting parameters of the optical pulsed trains
The classical theory of signal sampling is based on multiplying a signal by a train of impulses with a negligible duration.We use a simple and robust optical pulse source that generates relatively broad pulses (with a duration of about 25 ps) and the electrical bandwidth of our optical sampling system is broad (0-20 GHz).As a result, we had to take into account the pulses shape in our reconstruction algorithm.The non-ideal optical pulses create a different amplitude and phase for each signal replica of the sampled signal as given in Eq. ( 3).This effect is especially important in a system based on SMRS, since in different channels the same signal may be down-converted to baseband from different replicas.
Equation (3) indicates that there is no need to measure the entire pulse spectrum and it is sufficient to measure the pulses spectrum only at frequencies nF i , where n ≥ 0 is an integer number.While the amplitude coefficients of the pulse spectrum P i (nF i ) can be directly measured by using an RF spectrum analyzer, we retrieved the phases of the pulse coefficients P i (nF i ) indirectly by applying an optimization procedure to Eq. ( 4).We supply the system with an RF multi-band signal with known frequency bands that were measured by using RF spectrum analyzer.Since the signal frequencies are known we retain in the matrix Q, given in Eq. ( 4), only columns that correspond to those frequencies.Our system is not ideal and the noise is added to measurements.Therefore, Eq. ( 4) is not fulfilled with equality and we choose to find a least square solution to this equation.This solution x * S is obtained by multiplying both sides of Eq. ( 4) by a pseudo-inverse of the matrix Q S , denoted by Q † S [15]: Q † S depends on pulses parameters P i (nF i ) that we want to retrieve.We estimate these parameters by optimizing Eq. ( 7) to make the calculated baseband signals most consistent with the measured once.We use a square norm error criteria as a cost function for our optimization problem and substitute x * S back to Eq. ( 4).The resulting error equals: The parameters that should be optimized in Eq. ( 8) are P i (nF i ), i = 1 . . . 3 and the initial A/D sampling time τ.
Since the optimization problem ( 8) is highly non-linear and complex we use a MATLAB@ Genetic Algorithm Toolbox to obtain a global minimum for the objective function.After performing the genetic algorithm search, additional search was performed by using a local optimization function to improve the result accuracy.To speed up the optimization we use frequency resolution of ∆ f = 100 MHz.
Our RF source could generate signals in the frequency region between 6 and 12 GHz.Therefore, we could calibrate our system only at that frequencies.The calibration was performed by supplying three signals whose central frequencies were located in the frequency regions (6, 8.75) and (10, 12.25) GHz.The spectral width of each signal was about 150 MHz.This choice of the signal frequencies enables to obtain two pulse coefficients for each pulse train.Therefore, the number of parameters that should be optimized equals 13 − amplitude and phase for each pulses coefficient and the initial A/D sampling time τ.The optimized pulses coefficients amplitudes and phases were P 1 (2F 1 ) = 1.08 exp( j0.53), P 1 (3F 1 ) = 0.93 exp( j0.14) for the 3.8 GHz channel, P 2 (2F 2 ) = 1.20 exp( j0.77), P 2 (3F 2 ) = 0.95 exp( j0.52) for the 3.5 GHz channel, and P 3 (2F 3 ) = 0.65 exp( j0.89), P 3 (3F 3 ) = 0.48 exp( j3.18) for the 4.0 GHz channel.
We have verified that the ratio between the Fourier coefficients obtained by solving the optimization problem are approximately the same (up to 3% difference) as measured directly by using the RF spectrum analyzer.
To extract pulses parameters for the whole system frequency region (0-20 GHz) 6 signals at the following frequency regions may be used: (2, 5.25), (6, 8.75), (10, 12.25), (13.3, 14), (15.75, 17.1), and (19.25, 20) GHz.Assuming that each signal width equals about 100 MHz the resulting sampling matrix Q S will be small and the optimization can be performed in a reasonable time.Additional signals may be added to enhance the optimization accuracy.

Extracting the time offset of the electronic A/D τ
The electronic A/D conversion does not start at the time when the optical-pulse trains overlap.This effect is modeled by adding a time offset τ in Eq. ( 2).The time offset τ changes the sampling matrix in Eq. ( 4) and therefore it should be extracted.The optical pulse trains are synchronized to a reference signal with a frequency of 100 MHz.Therefore, we sampled this reference signal simultaneously with the other three channels by adding another low-rate A/D converter (fourth channel).At the calibration process the signals that are used to extract the parameters of the optical pulses are supplied to the system.The phase of the reference signal was also measured at the beginning of the sampling that is used for the calibration.Then, the value of τ was extracted by solving the optimization problem described in the previous subsection.The extraction of the delay τ should be performed only once.After completing the calibration the value of τ was updated by comparing the phase of the reference 100 MHz signal to the phase of the reference signal that was measured in the calibration process.

Experimental results
In our experiments we have sampled an RF source that could generate three signals simultaneously.Each signal had a bandwidth of 100−200 MHz and a carrier frequency between 6 to 12 GHz.The limitation on the carrier frequency is caused only due to our RF-signal generator and not due to the bandwidth of the sampling system.The sufficient sampling rate for the electronic A/D converter in each channel equals to the frequency of the corresponding pulse train.To simplify the hardware, the sampling rate of all the electronic A/D converters was set to 5 Gsamples/sec.This sampling rate is higher than theoretically required and therefore redundant frequency components at frequencies that are higher than the repetition rate of the corresponding pulse train were discarded.As a result, the effective total sampling rate equals 11.3 Gsample/sec.The total frequency support of the RF signals used in our experiment approximately equals 2 • 3 • 150 MHz = 900 MHz (the support of a real signal equals twice the bandwidth of its positive frequency components).The ratio between the total sampling rate to the signal support equals about 12.5.This ratio is higher than the minimum ratio of about 8 that was required in the numerical simulations [8].We note that since in our experiments we used the same reconstruction algorithm as in our numerical simulations, we expect that the minimum ratio that is required in the experiments will be about 8.The minimal theoretical sampling rate that is required for a blind multiband signal reconstruction equals twice the signal support.
The sampling window duration equals 8.2 µs that corresponds to a frequency resolution of 122 kHz.To reduce software runtime the frequency cell width was chosen to be ∆ f = 4 GHz/840 = 4.7619 MHz.This is consistent with our sampling rates requirement: After inverting the matrix and performing the signal reconstruction the frequency resolution of the reconstructed signal was increased to that of the sampled data (122 kHz) by setting β n = n∆ f /39, 0 ≤ n < 39, in Eq. ( 22) (Appendix A) where n is an integer number.To decrease the reconstruction error and to increase the robustness of the system to variations in system parameters we apply a simple algorithm that enhances the consistency of the sampling equations.The algorithm is described in Appendix B.
Figure 4 shows the reconstructed spectrum and a close up on each of the reconstructed signals.Figure 5 shows the baseband spectra measured by the A/D at the output of the sampling channels.The baseband spectra are calculated from the reconstructed signals (blue curves) using Eq. ( 4) and compared to the spectra of the measured baseband signals (red dashed curves).In this example, the signal with the highest carrier frequency (around 11.25 GHz) overlaps with two other signals at the baseband of all the sampling channels.We note that the signal frequency locations were not known before the reconstruction.An accurate reconstruction is obtained although the signal around 11.25 GHz overlaps after down-conversion with other signals and it can not be directly reconstructed from the baseband spectrum of any single channel.Therefore, such spectrum can not be reconstructed by using asynchronous multi-rate sampling scheme [7] with same sampling rates.In such sampling scheme each frequency of the signal should be unaliased in at least one of the sampling channels.In order to verify the reconstruction algorithm results we have compared the reconstructed spectrum to that measured by an RF spectrum analyzer as shown on Fig. 6.These results indicate that an excellent quantitative agreement is obtained between two spectra.The spectrum measurement by using the RF spectrum analyzer required many RF pulses.On the other hand, our system requires only a single RF pulse for measuring the entire sparse spectrum and hence it can be used in a real-time applications.
In cases where no overlap between different signals occurs in one of the sampling channels, the reconstruction could be also verified by up-converting the baseband signal in that channel according to Eq. ( 3). Figure 7 shows the reconstructed spectrum and a zoom on each signal (blue curve) in such case.Figure 8 shows the baseband spectra that are calculated from the reconstructed signal (blue curve) and compared to the spectrum of the measured baseband signals (red dashed curves).The down-converted signals in the experiment do not overlap at the 3.8 GHz channel.Therefore, the original spectrum could be directly retrieved from the baseband spectrum of that channel.The up-conversion was performed according to the signal frequencies that were measured by the RF spectrum analyzer.The up-converted spectrum which is shown by the green curve in Fig. 7 is in excellent quantitative agreement with the blind reconstruction of the spectrum (blue curve).
In our experiments we have detected and reconstructed more than 100 different spectra.An accurate reconstruction was obtained for all spectra that were checked.The worst run-time of    Since no aliasing occurs at the baseband of the 3.8 GHz channel, the original spectrum can be also directly calculated from the measured spectrum of that channel.An excellent agreement is obtained between the reconstructed and directly calculated spectra.  .Baseband spectrum of a signal at the 3.8 GHz sampling channel measured using the RF spectrum analyzer (green curve) and compared to the Fourier transform of the sampled baseband signal (blue curve).The input signal frequency equals 11.1 GHz and its power equals 6.5 dBm.The resolution bandwidth of the spectrum analyzer equals 100 kHz.The acquisition window of the sampling equals 8.19 µs that corresponds to a frequency resolution of about 122 kHz. the signal reconstruction algorithm was about 0.6 seconds when running in Matlab on Intel Pentium Core2 Duo 2.9 GHz machine with a 2 GBytes RAM Memory.Most of the computational effort was spent on calculating the Fourier transform of the sampled signals and on enhancing the consistency by the algorithm described in Appendix B. We notice, however that in all our experiments the signals were reconstructed without applying any pursuit algorithm, since the reduced matrix was full-rank for all cases.The run times can be dramatically improved by using a digital signal processor (DSP).

Measuring the system performance and the sampling jitter
The system performance can be characterized by the effective number of bits (ENOB).For a sinusoidal wave the ENOB is approximately connected to the signal-to-noise-plus-distortion ratio (SNDR) by [10]- [12]: ENOB = (SNDR(dB) − 1.76)/6.02.(9) To measure the SNDR for baseband signals we provided the system a sinusoidal wave with a frequency of 11.1 GHz.The frequency of the signal was chosen in a way that ensures that the second order distortion at 22.2 GHz is within the bandwidth of our system.The downconverted signal at the 3.8 GHz channel was sampled by the electronic A/D and it was also measured using the RF spectrum analyzer.When the signal power was increased to 6.5 dBm the second harmonic distortions were slightly higher than the noise floor.Figure 9 shows the spectrum of the down-converted signal measured by an RF spectrum analyzer with a resolution bandwidth of 100 kHz.Such resolution corresponds to the acquisition window on the order of 10 µs.The SNDR depends on the resolution bandwidth.Assuming a window duration of 10 µs the SNDR for a 100 kHz resolution bandwidth equals about 57 dB.This result corresponds to 9.2 effective bits.Assuming a window duration of 0.1 µs the resolution bandwidth equals 10 MHz and the SNDR equals 41.7 dB which corresponds to about 6.6 effective bits.One of the causes of noise added by the down-conversion is the timing jitter of the optical pulse-trains.
Assuming that the input signal is a sinusoidal wave with a frequency f c , the signal-to-noise ratio of the baseband signals is theoretically limited by the aperture jitter [12,13]:  We used Agilent E5052B Signal Source Analyzer for phase noise measurements.Figure 10 shows the single side-band phase noise spectrum L( f ) of three optical pulse-trains.
To estimate the standard deviation of the jitter we integrated the power spectrum L( f ) over the specified bandwidth [13]: where f 0 and f 1 define the integration region.The initial frequency f 0 is connected to the sampling window duration, T W , by f 0 = 1/4T W [14].In our scheme the sampling window duration equals 8.19 µs and therefore the initial integration frequency is f 0 = 30.5 kHz.The jitter in a frequency region of 31 kHz − 1 MHz equals 11.5, 13.1 and 17.1 fs for the 3.8, 3.5 and 4.0 GHz sampling channels, respectively.Using the measured phase noise and Eq.(10) we can estimate the SNR that is caused by the jitter.Assuming a sinusoidal input wave with a frequency of 20 GHz the jitter induced SNR equals 56.8 dB, 55.7 dB and 53.6 dB for the 3.8, 3.5 and 4.0 GHz channel, respectively.These results imply that for a sampling window duration of about 10 µs and a signal with a 20 GHz carrier frequency, the noise induced by the jitter is expected to be similar to the noise level measured in our system for the 11.1 GHz sine wave signal.By assuming a shorter window duration (larger resolution bandwidth) the amplitude noise is expected to be the dominant noise in the system.

Conclusions
We have demonstrated experimentally a new system for sampling and reconstruction of several bandwidth-limited signals with unknown carrier frequencies that can be located in a broad frequency region of 0-20 GHz.The system is based on a synchronous undersampling at several sampling rates.Accurate reconstruction was obtained even when the spectrum of some downconverted signals was deteriorated in all the sampling channels due to aliasing.Such spectra could not be reconstructed by using the asynchronous multi-rate sampling scheme [7].The optical pulse trains were generated by combining an electrical comb generator and an electroabsorption modulator.The locking of the sampling channels was obtained by providing the same reference signal to the PLL-DROs that operated at different frequencies.The reconstruction of signals that alias in part of the sampling channels is based on solving a set of linear equations.Therefore, the system should be calibrated and the parameters of each sampling channel should be extracted once in a long time (few days).We have developed methods for extracting the required system parameters.
We have demonstrated experimentally the reconstruction of three signals that were generated simultaneously.The bandwidth of each signal was 100-200 MHz and their carrier frequencies were located in a broad frequency region (6)(7)(8)(9)(10)(11)(12).The frequency locations of the signal were not known to the spectrum reconstruction algorithm.The system could accurately reconstruct signals in all the scenarios we tested (more than 100).It could also solve cases when one of the signals overlapped with other signals in the baseband region of all three sampling channels.The bandwidth of our system is about 20 GHz and the ENOB is about 9 bits, assuming a sampling window duration of 10 µs that corresponds to a resolution bandwidth of 100 kHz.

Appendix A
In this appendix we develop equations that connect the original signal to the down-converted signals in the non-ideal sampling system.
Equation (2) gives the connection between the continuous time original signal x(t) and the down-converted signal x i (t).
The connection between the spectrum of the original signal X( f ) and the down-converted signals X i ( f ) is given by and δ ( f ) is the Dirac delta function.By defining X( f ) = X( f ) exp ( j2π f τ) we obtain where f is a baseband frequency 0 ≤ f < F i .
To obtain matrix equations similar to ideal case analyzed in [8] we define an integer k and scalar β , 0 ≤ β < ∆ f , so that for any 0 ≤ f ≤ F i , f = k∆ f + β .Equation ( 13) can be written as We define vectors and matrices Using these definitions we rewrite Eq. ( 14) by By using the same frequency resolution ∆ f for all the sampling channels we are able to construct a system of linear equations that connect the original and the down-converted spectra.By defining M = F max /∆ f to be the number of cells in the support of the original signal X ( f ) Eq. ( 17) becomes Equation ( 18) can be written in a matrix-vector form.We define an M i × M matrix Q i whose elements are given by: Vectors x i (β ) and x (β ) are given by By substituting Eq. ( 19) and (20) to Eq. ( 18) we obtain the system of linear equation for i-th channel: For P sampling channels a single system of linear equations can be obtained: where the vector x (β ) and the matrix Q are obtained by concatenating vectors W H i (β )x i (β ) and matrices Q i W τ i as follows: . .

Appendix B
To enhance the consistency of the sampling equations we apply the following simple method after the reconstruction of the signal.In the first step we extract the support of the signal by using the reconstruction method as in [8].Let S be a set of columns corresponding to reconstructed signals frequency locations.We denote by Q S a sampling matrix with columns corresponding to these locations.This matrix is constructed by retaining corresponding columns in a matrix Q given in Eq. (22).Equation ( 22) connects the original signal to the Fourier transforms of the sampled signals.However, since our sampling model contains only the main effects and since the noise is present in the system, Eq. ( 22) is not fulfilled with equality.The solution we obtained minimizes the square norm of the error between the Fourier transform of the measured sampled signals and the baseband signal that are calculated from the reconstructed signal: x (β ) − Q S x S (β ).The solution vector x S is found by multiplying x (β ) by a pseudo-inverse matrix Q † S [15].In order to further reduce the error we define a set of coefficients α i , i = 1, 2, 3 . . .P that multiply each sampling channel.These coefficients compensate errors in the baseband spectrum.Assuming α i b i (β ) be a compensated baseband signal vector, where b i (β ) = W H i (β )x i (β ) as defined in Appendix A. Our goal is to find the parameters that enhance the consistency of our sampling model by minimizing the least squares error.Using our consistency criteria we obtain the following minimization problem: min x(β ), α 1 ,...,α P , . . .
where we optimize Eq. ( 23) for parameters α i and a reconstructed signal x S simultaneously.Let M = Q S Q † S − I.We define matrices M i with the number of columns equal to the length of the baseband vectors b i such that: i.e stacking the matrices M i in a row provides the matrix M. It is easy to verify that the minimization problem (23) can be converted to the following eigenvalue problem: where V is equal to a stacked matrix

Fig. 2 .
Fig. 2. Transfer function H i ( f ) of three sampling channels measured in the frequency region 5 MHz−2 GHz with a resolution of 5 MHz.

#Fig. 4 .
Fig. 4. Reconstructed spectrum of three signals and a zoom on each signal.The signal located around 11.25 GHz overlaps with other signals in the baseband region of all channels.

Fig. 5 .
Fig. 5. Baseband spectra of three sampling channels calculated from the reconstructed signal (blue curves) and compared to the spectrum of the measured baseband signals.

Fig. 6 .Fig. 7 .
Fig. 6.Comparison between the reconstructed spectrum (blue curve) and the spectrum measured by the RF spectrum analyzer (red dashed curve).

Fig. 8 .Fig. 9
Fig. 8. Baseband spectra of three sampling channels calculated from the reconstructed signal (blue curves) and compared to the spectrum of the measured baseband signals.

[M 1 b 1
|M 2 b 2 . . .|M P b P ].The eigenvector a corresponding to a smallest eigenvalue λ contains coefficients α i that are used as a gain for each channel sequence.The reconstructed signal x S is obtained by multiplying pseudo-inverse matrix Q † S .