Comparison of threshold-based algorithms for sparse signal recovery
77 th Mediterranean Conference on Embedded Computing MECO’2018, Budva, Montenegro
Comparison of threshold-based algorithms for sparse signal recovery
Tamara Koljenšić, Časlav Labudović
Faculty of Electrical Engineering University of Montenegro Podgorica, Montenegro email: [email protected], [email protected]
Abstract —Intensively growing approach in signal processing and acquisition, the Compressive Sensing approach, allows sparse signals to be recovered from small number of randomly acquired signal coefficients. This paper analyses some of the commonly used threshold-based algorithms for sparse signal reconstruction. Signals satisfy the conditions required by the Compressive Sensing theory. The Orthogonal Matching Pursuit, Iterative Hard Thresholding and Single Iteration Reconstruction algorithms are observed. Comparison in terms of reconstruction error and execution time is performed within the experimental part of the paper.
Keywords-compressive sensing, orthogonal matching pursuit, iterative hard thresholding, single iteration reconstruction I. I NTRODUCTION
The novel theory of Compressive Sensing (CS) is an intensively growing approach in signal processing. Traditional approach for signal reconstruction using acquired samples, follows the Shannon-Nyquist sampling theorem. For signals with high maximal frequency Shannon-Nyquist sampling procedure can result in a large number of samples. Against the traditional perspective in data acquisition, the CS approach allows exact signal reconstruction using small set of randomly acquired samples. The CS offers the possibility to acquire less data then it is commonly done, but still to be able to reconstruct the entire information. This approach attracted much interest in the research community and found wide-ranging applications from astronomy, biology, communications, image and video processing, medicine, to radar. Compressive sensing opens the possibility for simplifying the acquisition devices and apparatus, reducing the number of sensors, acquisition time, and storage capacities. To enable efficient recovery by using the CS approach, the signals have to satisfy certain conditions such as sparsity and incoherence. Sparsity is one of the main requirements that should be satisfied. It implies that the signal in different domains: time, frequency or time-frequency domains has only small number of non-zero coefficients. If the incoherence property is satisfied, the reconstruction using small set of samples is assured. Signal recovery is based on a powerful mathematical algorithms for error minimization. Some of them, like basis pursuit, Dantzig selector and gradient-based algorithms rely on linear programming methods. As much as they’re accurate, they are computationally demanding and not always suitable in practical applications. Many alternative approaches, like greedy algorithms – Matching Pursuit and Orthogonal Matching Pursuit are being developed to reduce the computational complexity. Also, recently proposed threshold based algorithms provide high reconstruction accuracy with low computational complexity: (e.g. Iterative Hard Thresholding - IHT, Iterative Soft Thresholding - IST). In this paper the focus is on the commonly used threshold-based algorithms. The comparison of the several algorithms is done. The sinusoidal multicomponent signal is observed, which can be defined by using the following relation: ( ) , [0, 1] K j knkk x n a e n N , (1) where K denotes number of sinusoidal components in observed signal, while N is signal’s length. Usually, these signals have large number of samples N but small number of non-zero components in frequency domain ( K N ). The reconstruction of such signals using various solutions has been presented in this work. Noiseless signal cases are examined. The organization of this paper is as follows. In Section II, the fundamental concepts of CS theory and signal reconstruction method are given. The overview of the applied CS algorithms is given in the Section III. The CS application to band-limited sparse signals is discussed in Section IV. Concluding remarks are given in Section V. II. B ASICS OF THE C OMPRESSIVE S ENSING APPROACH
In order to reconstruct the signal with high accuracy, traditional methods require sampling frequency to be twice the maximal signal frequency. This leads to a high number of signal samples that are acquired and stored. Beside the large number of samples, another problem is presence of noise, which can lead to missing signal information. New theory, the CS theory, overcomes those limits. It directly senses the data in a th Mediterranean Conference on Embedded Computing MECO’2018, Budva, Montenegro compressed form – i.e. at a lower sampling rate, and allows recovering the intentionally omitted or missing signal samples. Let f be an N -dimensional signal of interest, which is sparse in the transform domain represented with the transformation matrix R N×N Ψ . The sparse representation of f over the basis Ψ is represented by the vector p . Then f can be given by f = Ψp . (2) If Ψ is e.g. the inverse Fourier transform (FT) matrix, then p can be regarded as the frequency domain representation 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 non-zero. If a signal is able to be sparsely represented in a certain domain, the CS technique can be invoked to take only a few linear and non-adaptive measurements. The set of random measurements are selected from signal ( 1) N f × , which can be written by using random measurement matrix ( ) M N Φ as follows: d = Φf (3) From (2) and (3): d = ΦΨp (4) The incoherence is another important condition that basis matrix Ψ and measurement matrix Φ should satisfy to make the CS reconstruction possible. The relation between the number of nonzero samples in the transform domain Ψ and the number of measurements (required for reconstruction) depends on the coherence between the matrices Ψ and Φ . More specifically, a good measurement will pick up a little bit of information of each component in p based on the condition that Φ is incoherent with Ψ. As a result, the extracted information can be maximized by using the minimal number of measurement. III. A LGORITHMS FOR SPARSE RECOVERY
In this section, some of the commonly used threshold-based algorithms are described. The subjects of this analysis are Orthogonal Matching Pursuit (OMP), Iterative Hard Thresholding algorithm (IHT) and Single Iteration Reconstruction algorithm (SIRA). Knowing CS sensing matrix Ω Ψ and measurement vector d , the OMP algorithm approximates signal f as linear combination of columns in Ω . In each iteration, a set of columns is expanded with additional column that correlates best with the residual signal. The algorithm terminates until residual falls below determined threshold. OMP can be summarized as follows: 1. Variables initialization: set the approximation error r = d , the initial solution to 𝐟 and 𝐒 = Ø. 2. Do following steps till the stopping criterion is met: a) arg max n n i n i
S S r Ω , (5) b)
21 1 1 2 arg min , n f n n n f r S f (6) c) n n n n r r - S f , (7) d) n n i n i n n S S r Ω (8) until 𝑛 ≤ 𝐾 , where 𝐾 is number of signal components. The IHT algorithm [12] is an iterative algorithm which uses non-linear operator to reduce the 𝑙 - norm in each iteration. The algorithm solves the following problem: min f d Ωf f (9) From the optimization problem described by(9), the following iterative algorithm is derived. The non-linear operator is denoted as ( ) k H a and sets all but the largest (in terms of magnitude) k elements of a to zero. For given 𝒇 , the algorithm iterate: ( ( )), n n nk H T f f Ω d Ωf (10) Until either 𝑘 > 𝑁 𝑚𝑎𝑥 or n d Ωf 𝐻 𝑘 is defined as follows:
0, ( ) , ik i i i f kH f f f k (11) The Single Iteration Reconstruction Algorithm - SIRA is based on threshold calculation and choice of initial discrete Fourier transform (DFT) components that are above defined threshold. Initial DFT is calculated based on a set of available signal samples. It is shown that components above the threshold correspond to the signal components, while the components below the threshold are considered as noise. We make our analysis on the assumption that some random samples are omitted. The noise appears in signal and variance of it can be modeled as: var ( ... ),1 a k
NM A A AN (12) where M is the number of missing samples, while a N is the number of available samples. This variance will be used in threshold calculation: N T PN (13) The P is the probability that (N-K) components that correspond to noise, are lower than the threshold. The samples positions in the initial DFT that are above the threshold are used for the calculation of the exact DFT signal amplitudes, while the other th Mediterranean Conference on Embedded Computing MECO’2018, Budva, Montenegro frequency positions are filled with zeros. Vector of the initial DFT is formed using the available time domain signal samples, i.e. using the vector of measurements. Let v( M ) denotes the measurement vector. The initial DFT vector V is therefore formed as: ( ) ( ) , 1,..., a N j fa Na
V f v a e f N (14) Positions of the components above the threshold are obtained using following relation: arg{ }. pos V T (15) We found frequency positions, but to find the exact amplitudes on those positions requires solving minimization problem: ( ) ( ) CS CS CS v * X A A A (16) where CS matrix cs A is formed as , ( ) v P pos A , * cs A is Hermitian transpose of the matrix cs A and v P denotes vector of the available signal samples positions. IV. CS RECONSTRUCTION
APPLIED ON SPARSE
BAND-LIMITED
SIGNAL In this section, the three described algorithms, OMP, IHT and SIRA, are tested on the sparse band-limited signal, consisted of 7 components. The components’ magnitudes differ significantly from component to component. The signal is described using the following relation: ... j n N j n N j n N j n Nj n N j n N j n N x A e A e A e A eA e A e A e (17) where component magnitudes are:
A A A A A A A an N is the length of the signal and n ∈ (1, N ). Signal is considered as sparse in the DFT domain. In terms of the reconstruction error and the execution time, the same number of measurements is used: 200 samples (39.06% of the total signal length), 225 (43,94% of the total signal length), 250 (48,82% of the total signal length), 275 (53,71 % of the total signal length) and with M=300 samples (58,59% of the total signal length). Fig. 1a shows time and Fig. 1b shows DFT domain of the original signal. a) b) Figure 1: a) Zoomed time domain of the original, b) Fourier domain of the original signal
Table 1 shows the error and reconstruction time for the obtained measurements. The error is defined as a maximum magnitude difference between the original DFT and the reconstructed one. Table 1: Error and reconstruction time for different numbers of samples
Number of samples M
ALGORITHM
SIRA OMP IHT
TIME(sec)
ERROR TIME(sec)
ERROR TIME(sec)
ERROR TIME(sec)
ERROR TIME(sec)
ERROR From the Table 1 it can be seen that the SIRA algorithm is the fastest with the most accurate reconstruction of the original signal. For M between 200 and 300, the minimum error was obtained using SIRA. Fig.2 shows the different reconstructions for the minimum number of measurements th Mediterranean Conference on Embedded Computing MECO’2018, Budva, Montenegro
Figure 2: Comparison of the different algorithms reconstruction for: SIRA (M=170), IHT (M=34) and OMP (M=200), from the top to the bottom
OMP algorithm was able to guess each position of samples more frequently than SIRA but with less accuracy. IHT could reconstruct the original signal perfectly, but with the information about expected number of components. Fig.3 shows IHT reconstruction when the number of components is larger than predicted.
Figure 3: IHT reconstruction in case of unknown expected components Fig.4 shows error for the minimal measurement numbers needed to successfully reconstruct original signal a) SIRA for M=170; b) IHT for M=34 c) OMP for M=200 Fig. 5a shows the error for all three algorithms for number of measurements range from 200 to 300 and Fig. 5b shows zoomed error for SIRA and IHT a) b) Figure 5: a) Error for the measurements between 200 and 300 b) Zoomed errors for SIRA, IHT V.
CONCLUSION
Performance of the several reconstruction algorithms applied on sparse band-limited signal are considered in the paper. Signal is being reconstructed changing the number of measurements. Comparing maximum error and execution time, the best results were obtained by using SIRA. It is important to acknowledge -11 -4
200 220 240 260 280 300050100150200200 220 240 260 280 30002468 x 10 -6 th Mediterranean Conference on Embedded Computing MECO’2018, Budva, Montenegro that OMP was able to find all samples’ positions in the same number of iterations as SIRA but error was significantly bigger. IHT algorithm demands a priori knowledge of the signal, as one of his inputs is the number of components we are searching for. In case of familiar signal, IHT used 6, 65% of the signal length to reconstruct signal successfully. That's 20% of SIRA's requirements but SIRA can work without previous signal's background. R
EFERENCES [1]
Z. Qin, J. Fan, Y. Liu, Y. Gao, G. Ye Li, “Sparse Representation for Wireless Communications”, https://arxiv.org/abs/1801.08206 [2]
T. Blumensath, M. E. Davies, “Iterative Thresholding for Sparse Approximations”, Journal of Fourier Analysis and Applications, vol. 14, no. 5-6, pp 629-654, December 2008 [3]
S. Stankovic, I. Orovic, E. Sejdic, "Multimedia Signals and Systems: Basic and Advance Algorithms for Signal Processing," Springer-Verlag, New York, 2015. [4]
I. Orović, S. Stanković, “Improved higher order robust distributions based on compressive sensing reconstruction”, IET Signal Processing, 8 (7), 738-748. [5]
E. Candes, J. Romberg, “l1-magic: Recovery of Sparse Signals via Convex Programming”, October 2005. [6]
S. Stankovic, I. Orovic, LJ. Stankovic, “An Automated Signal Reconstruction Method based on Analysis of Compressive Sensed 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 Theory, vol.56, no.9, pp. 4395-4401, Sept. 2010. [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, "Greed 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, LJ. Stankovic, "Combination of Gradient Based and Single Iteration Reconstruction Algorithms for Sparse Signals," 17th IEEE International Conference on Smart Technologies, IEEE EUROCON 2017, 6th-8th July 2017, Ohrid, Macedonia, 2017. [11]
A. Draganic, I. Orovic, N. Lekic, M. Dakovic, S. Stankovic, "Architecture for Single Iteration Reconstruction Algorithm," 4th Mediterranean Conference on Embedded Computing, MECO 2015 [12]
R. Mihajlovic, M. Scekic, A. Draganic, S. Stankovic, "An Analysis of CS Algorithms Efficiency for Sparse Communication Signals Reconstruction," 3rd Mediterranean Conference on Embedded Computing , MECO, 2014 [13]