Comparison of threshold-based algorithms for sparse signal recovery

7 th Mediterranean Conferenc e on Embedde d Compu ting MECO’20 18, Budva, M ontenegro Compar is on of threshold-based algorithms for sparse signal recovery Tamara Kol j en šić, Časlav Labudov ić Faculty of Electrical E ngineering University of Montene gro Podgorica, Montenegro email: koljensict@ gmail.com, caslavlabudovic@ gmail.com Abstract — Intensively growing approach in signal processing and acquisition, the Compressive Sensing ap proach, allows sparse signals to be recovered from small number of r andomly acquired signal coefficients . This paper analyses s ome of the c ommonly u sed threshold-based algorithms for sparse signal reconstruction . Signals satisfy the co nditions requi red by the Compressive Sensing theory . The Orthogonal Matching Pursuit , Iterative Hard Thresholding and Single Iteration R econstruction algorith ms are observed. Comparison in terms of reconstruction error and execution ti me is performed within the ex perimental part of the paper. Keywords-compressive sensing, orthogonal matching p ursuit, iterative hard thresholding, single iteration reconstruction I. I NTRODUCTI ON The novel theory of Compressiv e Sensing (CS) is an intensively growing approach in signal processing. Tradit ional approach for signal rec onstruction using acqui red samples, follow s th e Shan non-Nyquist sampling theo rem. Fo r s ignals with high maximal frequency Shannon-Nyquist samplin g procedur e can result in a large number of samples . Against the traditi onal perspective in data acquisiti on, the CS approach allow s exact signal reconst ruction using small set o f rando mly acquired sam ples. The CS o ffers the possibility to acquire less data then it is comm only done, but st ill to b e able to re construct the entire information . This approach attra cted much interest in the research com munity and found wide-ranging applications from astron omy, biology , commun ications, image and vi deo processing, medicine , to radar . Compressive sen sing opens the possibili ty for simplif ying the ac quisition devices and apparatu s, reducing the number of sens ors, acquisition time, an d st orage capaciti es. To enable efficient recovery by using the CS approach , the signals have t o sa tisfy certain c ondition s such as s parsity and incoheren ce. Sparsity is on e of th e main requirements that should be satisfie d. It implies that the signal in different domains: tim e, frequency or time-frequency d omains has only small num ber o f non-zero coefficients. If the incohere nce property is satisfi ed, the reconstruction using sm all set o f samples is assure d. Signal recove ry is based on a powerful m athematical algorithm s for err or minim ization. Some of th em, lik e basis pursuit, Dantzig selector and gradi ent-based algor ithms rely on linear p rogramm ing methods. As much as they’ re accu rate, they are computat ionally deman ding and not alw a ys suitable in practical applicat ions. Many al ternative approach es, lik e gree dy algorithm s – Matchin g Pursuit and Orthogon al Matching Pursuit are being devel oped to reduc e the computational complexity . Also, recently proposed thresh old based algori thms provide hig h reconstruct ion a ccuracy w ith low computati onal co mplexi ty: (e.g. Iterat ive Hard Thresholdin g - IHT, Iterativ e Soft Threshol ding - IST) . In this paper the focus is on the comm only used threshold- based algorith ms. The comparison o f the several algorithm s is done. The s inusoida l multicom ponent signal is observed, w hich can be defined by using the f ollowin g relation: 2 1 ( ) , [0, 1 ] K j kn k k x n a e n N        , (1) where K deno tes nu mber of sinus oidal components in obse rved signal, while N is signal’s length . Usu ally, these sig nals have large number of sam ples N but small number o f non-zero components in frequency domain ( K N ). The reconstructi on of such signals using various soluti ons has been presente d in this work. Noisel ess sign al cases are exam ined. The organization of this pape r is as follow s. In Section II, the fundam ental concepts of CS theory and signal reconstructi on method are giv en. The overview of the ap plied CS alg orithm s is given in the Section III. The CS application to b and-limited sparse signals is discussed in Section IV. Concluding remark s are given in Sect ion V. II. B ASI CS OF THE C OMPRESSIV E S ENSING APP ROACH In order to reconstruct t he signal w ith high ac curac y, traditional methods req uire sa mpling frequenc y to be t wice the maximal sig nal freque ncy. Th is leads to a high number of sig nal samples t hat ar e acquired and stored. B eside the large number of samples, a nother problem is p resence of noise, w hich can lead to missing si gnal infor mation. Ne w theory, the CS theo ry, overcomes those limits. It directly senses the data in a 7 th Mediterranean Conferenc e on Embedde d Compu ting MECO’20 18, Budva, M ontenegro compressed for m – i.e. at a lower sampling rate, and allows recovering the intentionally o mitted or missing signal samples . Let f be an N -dimensional signal of interest, which is sparse in the tran sform do main repr esented with the tr ansformation matrix R  N× N Ψ . T he sparse representation of f over the basis Ψ is represented by the vector p . Then f can b e given by f= Ψp . (2) If Ψ i s e.g. the inv erse F ourier tr ansform (F T) matr ix, then p can be regar ded as the frequen cy domain represent ation of the time domain signal, f . Signal f is said to be K -sparse in the Ψ domain if there are only () K K N  out of N coefficients in p that are n on-zero. If a signal is able to be s parsely repres ented in a certain domain, the CS techn ique can be invoked to take only a few lin ear and n on-adaptiv e measurem ents. The set o f ra ndom measure ments are selected fro m s ignal ( 1 ) N f× , w hich can b e w ritte n by using random measurement matrix () MN  Φ as follows: d= Φf (3) From (2) and (3 ): d= Φ Ψ p (4) The in coherence i s a nother important co ndition that basis matrix Ψ and measurement matrix Φ should satisfy to make the CS reconstructio n po ssible. The relation bet ween the number of nonzero sa mples in the tra nsform do main Ψ a nd the number of m eas urements (required for reconstruction) depends on the coherence b etween the matrices Ψ and Φ . More specifically, a good measurement will pick up a little b it of information of each componen t in p based on the condition that Φ is incohere nt with Ψ . As a result, the extracted in formation can be maximized by using the minimal number of measurement. III. A LGO RITHMS FOR SPA RSE RECOVERY In this section , some of the co mmonly used threshol d-based algorithm s are described. T he subjects o f this analy sis are Orthogona l Matching Pursuit (OMP), Iterative Hard Threshol ding alg orithm (I HT) and Single Iteration Reconstruc tion algo rithm (S I RA ). Knowin g CS sensing matrix  Ω Ψ  and measurem ent vector d , the OMP algorith m approxim ates signal f as linear combinati on o f columns in Ω . In each iterat ion, a set of columns is expande d w ith additional colum n that correlates best with the residual signal. The algorith m terminates until residual falls below determ ined threshold. OMP can be summ arized as follow s: 1. Var iables in itializ ation: set the approxim ation error 0 r = d , the initia l soluti on to   and    Ø. 2. Do foll owing steps till the st opping crit erion is m et: a) 1 1, ar g ma x n n i n i      S S r Ω , (5) b) 2 1 1 1 2 ar g min , n f n n n     f r S f (6) c) 11 n n n n   r r - S f , (7) d) 11 1 and arg max | n n i n i nn          S S r Ω (8) until    , wh ere  is num ber of signal components. The IHT algorithm [12] is an ite rative alg orithm w hich uses non-linea r operator t o r educe the   - n orm in each iteration. The algorithm solves the foll owing pr oblem: 0 m in f   d Ωf f (9) From the opt imizati on problem described by (9 ), the follow ing iterative algorithm is derived. The n on-linear ope rator is denoted as () k Ha and sets all but the largest (in term s of magn itude) k elements of a to zero. For given   , the algorithm iterate: 1 ( ( )) , n n n k H     T ff Ω d Ω f ( 10 ) Until eithe r     or 2 n   d Ωf   is define d as follow s: 0, () , i ki ii fk Hf f f k        ( 11 ) The Single Iterati on Reconstr uction Algorith m - SIRA is based on threshold calcula tion an d ch oice of ini tial discr ete Fourier transform (DFT) components that are above defined threshol d. Initial DFT is calculate d based on a set o f availa ble signal s amples. It is show n th at com ponents abov e the thresh old correspon d t o the signal com ponents, while the components below the threshol d are c onsider ed as no ise. We make our analysis on the assumption that some random samples are omitted. The noise appears in signal and varian ce o f it can be modeled as: 2 2 2 12 v ar ( ... ), 1 a k N M A A A N      ( 12 ) where M is the number of missing samples, while a N is the number of available samples. T his varia nce will be used in threshold calculation : 2 1/ 2 10 1 ( va r lo g ( 1 ) N TP N    ( 13 ) The P is th e probability that (N-K) components that correspond to noise, are lower than the threshold. T he samples positions in the initia l DFT that are above the threshold are used for the calculation of the exact DFT signal amplitudes, while the other 7 th Mediterranean Conferenc e on Embedde d Compu ting MECO’20 18, Budva, M ontenegro frequency p ositions are filled with zeros. Vector of the initi al DFT is form ed using the available time domain signal sam ples, i.e. using the vector of measurements. Let v( 1 M  ) denotes the measurement vector. The initial DFT vector V is theref ore formed as: 2/ 1 ( ) ( ) , 1 , ..., a N j fa N a V f v a e f N      ( 14 ) Positions of the co mponents abo ve the thres hold are obtained using following relation: ar g{ }. po s V T  ( 15 ) We found frequency positions, but to find the exact amplitudes on those positions requires solv ing m inimization pro blem: 1 ( ) ( ) CS CS CS v   * X A A A ( 16 ) where CS matrix cs A is form ed as , () v P po s A , * cs A is Hermitian transpose o f the matrix cs A and v P denotes vector o f the available signal sam ples posi tions. IV. CS RECONSTRUCT ION APPLIED ON SPARSE BAND-LIMIT ED SIGNAL In t his section, the three described algorithms, OMP, I H T and SIRA, are tested on th e sparse band-limited signal, consisted of 7 components. The co mponents’ magnitudes dif fer significantly fro m co mponent to component. The signal is described using the following relation: 2 32 / 2 38 / 2 130 / 2 148 / 1 2 3 4 2 272 / 2 415 / 2 435 / 567 ... j n N j n N j n N j n N j n N j n N j n N x A e A e A e A e A e A e A e             ( 17 ) where co mponent magnitudes are: 1 2 3 4 5 6 7 3.5 ; 3 ; 1.75 ; 2.5 ; 3.75 ; 2.3 ; =3.3 ; A A A A A A A       an N is the length o f the signal and n  (1 , N ). Sig nal is considered as sp arse in the DFT domain. In ter ms of the r econstruction err or and the execution time , the same number of measure ments is u sed: 200 samples (39.06% of the total signal le ngth), 225 ( 43,94% of the total signal le ngth), 250 (48 ,82% of the total sig nal le ngth), 275 (53,71 % of the total signal length) and with M=3 00 samples (58,59% of the total signal length). Fig. 1a shows ti me and Fig. 1b shows DFT d omain of the original s ignal. a) b) Figure 1: a) Zoomed time domain of the original, b) Fourier domain of t he original signal Table 1 shows the error and reco nstruction time for the obtained measure ments. T he error is d efined as a maximum magnitude difference bet ween the original DFT and the reconstructed one. Table 1: Er ror and reconstruction ti me for different numbers of samples Number of samples M ALGORITHM SIRA OMP IHT 200 TIME(sec) 0. 01 39 45 0.0 984 06 0.0 552 30 ERROR 11 3 .5 3 9 1 1 0   15 5 .12 99 6 6 .6 9 8 1 0   225 TIME(sec) 0.0 135 77 0.0 764 67 0.0 670 02 ERROR 11 3 .5 3 9 1 1 0   81.3 499 6 5 .8 89 2 1 0   250 TIME(sec) 0. 01 65 28 0. 09 56 28 0.0 686 34 ERROR 11 3 .5 3 9 1 1 0   70.1 819 6 7 .0 7 3 1 0   275 TIME(sec) 0. 01 40 68 0. 08 99 88 0. 06 33 23 ERROR 11 3 .5 3 9 1 1 0   9.7435 6 6 .2 6 3 1 0   300 TIME(sec) 0. 01 88 51 0. 10 39 83 0. 05 92 31 ERROR 11 3 .5 3 9 1 1 0   8.6813 6 6 .1 3 3 1 0   From the T able 1 it can be see n that the SIR A a lgorith m is the fastest with the most acc urate reconstruction o f the original signal. For M between 200 and 300, the minimum err or was obtained using SIRA. Fig.2 shows t he d ifferent reconstructio ns for the minimu m number of meas urements 0 10 20 30 40 50 60 70 80 90 100 0 5 10 15 20 0 50 100 150 200 250 300 350 400 450 500 0 500 1,000 1,500 2000 7 th Mediterranean Conferenc e on Embedde d Compu ting MECO’20 18, Budva, M ontenegro Figure 2: Compariso n of the differe nt algorithms recons truction for: SIRA (M=170), I HT (M=34) and OMP (M=200), from the top to the bottom OMP algorithm was ab le to guess each positio n of sampl es more frequently tha n SIRA but with less ac curacy. IHT could reconstruct the original si gnal perfectly , but with the information about expected number of compo nents. Fig.3 sho ws IHT reconstruction when the number of components is lar ger than predicted . Figure 3 : IHT re construction in case of unknow n expected components Fig.4 shows error for the minimal measureme nt numbers needed t o succ essful ly reconstruct origi nal signal a) SIRA for M=17 0; b) IHT for M= 34 c) OMP for M= 200 Fig. 5a shows the error for all three algorithms for number of measurements range fro m 2 00 to 300 and Fig. 5b shows zoomed error for SIRA and I HT a) b) Figure 5: a) Er ror for the measurements between 200 and 300 b) Zoom ed error s for SIRA, IHT V. CONCLUSION Performance of the several reconstruction algorith ms applied on sparse band -limited sign al ar e co nsidered in the pap er. Signal is being reconstructed changing the number of measurements. Comparing maximum error and execution time, the best results were obtained by usi ng SI RA. It is important to ackno wledge 0 50 100 150 200 250 300 350 400 450 500 0 500 1000 1500 2000 0 50 100 150 200 250 300 350 400 450 500 0 500 1000 1500 2000 0 50 100 150 200 250 300 350 400 450 500 0 500 1000 1500 2000 0 20 40 60 80 100 120 140 160 180 200 0 500 1000 1500 2000 0 100 200 300 400 500 -2 -1 0 1 2 3 4 x 10 -11 0 100 200 300 400 500 -1 -0 .5 0 0.5 1 x 10 -4 0 100 200 300 400 500 -4 00 -3 00 -2 00 -1 00 0 100 200 220 240 260 280 300 0 50 100 150 200 200 220 240 260 280 300 0 2 4 6 8 x 10 -6 7 th Mediterranean Conferenc e on Embedde d Compu ting MECO’20 18, Budva, M ontenegro that OMP was able to fi nd all sa mples’ po sitions in the sa me number of iterations as SIRA but error w as significantly bigger. IHT alg orithm d emands a priori knowledge of the si gnal, as one of h is inputs is the number of components we are searching for. In case of familiar signal, IHT used 6 , 65% of the signal lengt h to reconstruct sig nal succes sfully. That's 20 % of SI RA's requirements but SIR A can work witho ut previous si gnal's background. R EFERENCES [1] Z. Qin, J. Fan, Y. Liu, Y. Gao, G. Ye Li, “Sparse Representation for Wirele ss Communications”, https://ar xiv.org/abs/1801. 08206 [2] T . Blumensath, M. E. Davies, “Iterative Thresholding for Sparse Approximatio ns”, Journal of Fourier Analysis and Applications, vol. 14, no. 5-6, pp 629-6 54, December 2008 [3] S. Stankovic, I. Orovic, E. Sejdic, "Multimedia Signals and Systems: Basic and Advance Algorithms for Signal Processing," Springer-Verlag, New Yo rk, 2015. [4] I. Orović, S. S tanković, “Improved higher o rder robust distributions based on compressive sensing reconstruction”, IET Signal Processing, 8 (7), 738-748. [5] E. Candes, J. Romberg, “l1 -magic: R ecovery of Sparse Signals via Convex Pro gramming”, October 2005. [6] S. Stankov ic, I . Orovic, LJ. Stankovic, “A n Automated Signal Reconstruction Method based on Analysis of C ompressive S ensed Signals in Noisy Environment,” Signal Processing, vol . 104, Nov 2014, pp. 43 - 50, 2014. [7] M. A. Davenport, M. B. Wakin, " Analysis of Orthogonal Matching Pursuit Using the Restricted Isometry Property," IEEE Transactions on Information T heory, vol .56, no.9, pp. 4395-4401, Sept. 201 0. [8] LJ. Stankovic, S. Stankovic, M . Amin, "Missing Samples Analysis in Signals for Applications to L-estimation and Compressive Sensing," Signal Processing, vol. 94, Jan 2014, pp. 401-408, 2014. [9] J. A. Tropp, "Gree d is good: algorithmic results for sparse approximation," IEEE Transactions on Information Theory, vol.50, no.10, pp.2231-2242 , Oct. 2004. [10] S. Stankovic, S. Vujovic, I . Orovic, M. Dakovic, L J. Stankovic, "Combination o f G radient Base d and Single I t eration Reco nstruction Algorithms for Sparse Signals," 17th IEEE I nternational Conference on Smart Technologies, I EEE EUROCON 2017, 6th -8th July 2017, Ohrid, Macedonia, 20 17. [11] A. Draganic, I . Orovic, N. L ekic, M. D akovic, S. Stankovic, "A rchitecture for Single Iteration Reconstruction Algorithm," 4th Mediterranean Conference o n Embedded Comp uting, MECO 20 15 [12] R. Mihajlovic, M . Sce kic, A . Drag anic, S. Stankov ic, "An A nalysis of CS Algorithms Efficiency for Sparse Communication Signals Reconstruction," 3rd Me diterranean Confer ence on Embedded Computing , ME CO, 2014 [13] A. Draganic, I. Orovic, S. Stankovic, "On some common compressive sensing recove ry algorithms a nd a pplications - Revie w paper," Facta Universitatis, Series: El ectronics and Energe tics, accepted for publication, 20 17