# Compressive Sensing in Signal Processing: Algorithms and Transform Domain Formulations

Authors: Irena Orović; Vladan Papić; Cornel Ioana; Xiumei Li; Srdjan Stanković
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Compressive sensing has emerged as an area that opens new perspectives in signal acquisition and processing. It appears as an alternative to the traditional sampling theory, endeavoring to reduce the required number of samples for successful signal reconstruction. In practice, compressive sensing aims to provide saving in sensing resources, transmission, and storage capacities and to facilitate signal processing in the circumstances when certain data are unavailable. To that end, compressive sensing relies on the mathematical algorithms solving the problem of data reconstruction from a greatly reduced number of measurements by exploring the properties of sparsity and incoherence. Therefore, this concept includes the optimization procedures aiming to provide the sparsest solution in a suitable representation domain. This work, therefore, offers a survey of the compressive sensing idea and prerequisites, together with the commonly used reconstruction methods. Moreover, the compressive sensing problem formulation is considered in signal processing applications assuming some of the commonly used transformation domains, namely, the Fourier transform domain, the polynomial Fourier transform domain, Hermite transform domain, and combined time-frequency domain.


---

## Body

## 1. Introduction

The fundamental approach for signal reconstruction from its measurements is defined by the Shannon-Nyquist sampling theorem stating that the sampling rate needs to be at least twice the maximal signal frequency. In the discrete case, the number of measurements should be at least equal to the signal length in order to be exactly reconstructed. However, this approach may require large storage space, significant sensing time, heavy power consumption, and large number of sensors. Compressive sensing (CS) is a novel theory that goes beyond the traditional approach [1–4]. It shows that a sparse signal can be reconstructed from much fewer incoherent measurements. The basic assumption in CS approach is that most of the signals in real applications have a concise representation in a certain transform domain where only few of them are significant, while the rest are zero or negligible [5–7]. This requirement is defined as signal sparsity. Another important requirement is the incoherent nature of measurements (observations) in the signal acquisition domain. Therefore, the main objective of CS is to provide an estimate of the original signal from a small number of linear incoherent measurements by exploiting the sparsity property [3, 4].The CS theory covers not only the signal acquisition strategy, but also the signal reconstruction possibilities and different algorithms [8–17]. Several approaches for CS signal reconstruction have been developed and most of them belong to one of three main approaches: convex optimizations [8–11] such as basis pursuit, Dantzig selector, and gradient-based algorithms; greedy algorithms like matching pursuit [14] and orthogonal matching pursuit [15]; and hybrid methods such as compressive sampling matching pursuit [16] and stage-wise OMP [17]. When comparing these algorithms, convex programming provides the best reconstruction accuracy, but at the cost of high computational complexity. The greedy algorithms bring about low computation complexity, while the hybrid methods try to provide a compromise between these two requirements [18].The proposed work provides a survey of the general compressive sensing concept supplemented with the several existing approaches and methods for signal reconstruction, which are briefly explained and summarized in the form of algorithms with the aim of providing the readers with an easier and practical insight into the state of the art in this field. Apart from the standard CS algorithms, a few recent solutions have been included as well. Furthermore, the paper provides an overview of different sparsity domains and the possibilities of employing them in the CS problem formulation. Additional contribution is provided through the examples showing the efficiency of the presented methods in practical applications.The paper is organized as follows. In Section2, a brief review of the general compressive sensing idea is provided together with the conditions for successful signal reconstruction from reduced set of measurements and the signal recovery formulations using minimization approaches. In Section 3, the commonly used CS algorithms are reviewed. The commonly used domains for CS strategy implementation are given in Section 4, while some of the examples in real applications are provided in Section 5. The concluding remarks are given in Section 6.


## 2. Compressive Sensing: A General Overview

### 2.1. Sparsity and Compressibility

Reducing the sampling rate using CS is possible for the case of sparse signals that can be represented by a small number of significant coefficients in an appropriate transform basis. A signal havingK nonzero coefficients is called K-sparse. Assume that signal s exhibits sparsity in certain orthonormal basis Ψ defined by the basis vectors {
ψ
1
,
ψ
2
,
…
,
ψ
N
}. The signal s can be represented using its sparse transform domain vector x as follows:(1)
s
n
=
∑
i
=
1
N
x
i
ψ
i
n
.In matrix notation, the previous relation can be written as(2)
s
=
Ψ
x
.Commonly, the sparsity is measured using the l
0-norm, which represents the cardinality of the support of x: (3)
x
0
=
c
a
r
d
supp
⁡
x
=
K
.In real applications, the signals are usually not strictly sparse but only approximately sparse. Therefore, instead of being sparse, these signals are often called compressible, meaning that the amplitudes of coefficients decrease rapidly when arranged in descending order. For instance, if we consider coefficients x
1
≥
x
2
≥
⋯
≥
x
n, then the magnitude decays with a power law if there exist constants C
1 and q
>
0 satisfying [19](4)
x
i
≤
C
1
i
-
q
,where larger q means faster decay and consequently more compressible signal. The signal compressibility can be quantified using the minimal error between the original and sparsified signal (obtained by keeping only K largest coefficients):(5)
σ
K
x
p
=
min
⁡
x
-
x
K
p
.


### 2.2. Conditions on the CS Matrix: Null Space Property, Restricted Isometry Property, and Incoherence

Instead of acquiring a full set of signal samples of lengthN, in CS scenario, we deal with a quite reduced set of measurements y of length M, where M
<
N. The measurement procedure can be modeled by projections of the signal s onto vectors {
ϕ
1
,
ϕ
2
,
…
,
ϕ
M
} constituting the measurement matrix Φ:(6)
y
=
Φ
s
.Using the sparse transform domain representation of vectors given by (2), we have(7)
y
=
Φ
Ψ
x
=
A
x
,where A will be referred to as CS matrix.In order to define some requirements for the CS matrixA, which are important for successful signal reconstruction, let us introduce the null space of matrix. The null space of CS matrix A contains all vectors x that are mapped to 0:(8)
N
A
=
x
:
A
x
=
0
.In order to provide a unique solution, it is necessary to provide the notion that twoK-sparse vectors x and x
′ do not result in the same measurement vector. In other words, their difference should not be part of the null space of CS matrix A:(9)
A
x
-
x
′
≠
0
.Since the difference between twoK-sparse vectors is at most 2K-sparse, then aK-sparse vector x is uniquely defined if null space of A contains no 2K-sparse vectors. This corresponds to the condition that any 2K columns of A are linearly independent; that is, (10)
s
p
a
r
k
A
>
2
K
,and since s
p
a
r
k
A
∈
[
2
,
M
+
1
], we obtain a lower bound on the number of measurements:(11)
M
>
2
K
.In the case of strictly sparse signals, the spark can provide reliable information about the exact reconstruction possibility. However, in the case of approximately sparse signals, this condition is not sufficient and does not guarantee stable recovery. Hence, there is another property called null space property that measures the concentration of the null space of matrix A. The null space property is satisfied if there is a constant γ
∈
(
0,1
) such that [19](12)
a
S
1
=
γ
a
S
c
1
∀
a
∈
N
A
,for all sets S
⊂
{
1
,
…
,
N
} with cardinalityK and their complements S
c
=
{
1
,
…
,
N
}
∖
S. If the null space property is satisfied, then a strictly K-sparse signal can be perfectly reconstructed by using l
1-minimization. For approximately K-sparse signals, an upper bound of the l
1-minimization error can be defined as follows [20]:(13)
x
-
x
^
1
≤
2
1
+
γ
1
-
γ
σ
K
x
1
,where σ
K
(
x
)
1 is defined in (5) for p
=
1 as the minimal error induced by the best K-sparse approximation.The null space property is necessary and sufficient for establishing guarantees for recovery. A stronger condition is required in the presence of noise (and approximately sparse signals). Hence, in [8], the restricted isometry property (RIP) of CS matrix A has been introduced. The CS matrix A satisfies the RIP property with constant δ
K
∈
0,1 if(14)
1
-
δ
K
x
2
2
≤
A
x
2
2
≤
1
+
δ
K
x
2
2
,for every K-sparse vector x. This property shows how well the distances are preserved by a certain linear transformation. We might now say that if the RIP is satisfied for 2K with δ
2
K
<
1, then there are no twoK-sparse vectors x that can correspond to the same measurement vector y.Finally, the incoherence condition mentioned before, which is also related to the RIP of matrixA, refers to the incoherence of the projection basis Φ and the sparsifying basis Ψ. The mutual coherence can be simply defined by using the combined CS matrix A as follows [21]:(15)
μ
A
=
m
a
x
i
≠
j
,
1
≤
i
,
j
≤
N
A
i
,
A
j
A
i
2
A
j
2
.The mutual coherence is related to the restricted isometry constant using the following bound [22]:(16)
δ
K
=
K
-
1
μ
.


### 2.3. Signal Recovery Using Minimization Approach

The signal recovery problem is defined as the reconstruction of vectorx from the measurements y
=
A
x. This problem can be generally seen as a problem of solving an underdetermined set of linear equations. However, in the circumstances when x is sparse, the problem can be reduced to the following minimization:(17)
m
i
n
x
0
s.t.
y
=
A
x
.The l
0-minimization requires an exhaustive search over all N
K possible sparse combinations, which is computationally intractable. Hence, the l
0-minimization is replaced by convex l
1-minimization, which will provide the sparse result with high probability if the measurement matrix satisfies the previous conditions. The l
1-minimization problem is defined as follows:(18)
min
⁡
x
1
s.t.
y
=
A
x
,and it has been known as the basis pursuit.In the situation when the measurements are corrupted by the noise of levele
:
y
=
Φ
Ψ
x
+
e
=
A
x
+
e and e
2
≤
ε, the reconstruction problem can be defined in a form:(19)
min
⁡
x
1
s.t.
y
-
A
x
2
≤
ε
, called basis pursuit denoising. The error bound for the solution of (19), whereA satisfies the RIP of order 2
K with δ
2
K
<
2
-
1 and y
=
A
x
+
e, is given by(20)
x
-
x
^
2
≤
C
0
σ
K
x
1
+
C
1
2
ε
1
+
2
δ
2
K
,where the constants C
0 and C
1 are defined as [19] (21)
C
0
=
2
1
-
1
-
2
δ
2
K
1
-
1
+
2
δ
2
K
,
C
1
=
2
1
1
-
1
+
2
δ
2
K
.For a particular regularization parameter λ
>
0, the minimization problem (19) can be defined using the unconstrained version as follows:(22)
m
i
n
x
1
2
y
-
A
x
2
2
+
λ
x
1
,which is known as the Lagrangian form of the basis pursuit denoising. These algorithms are commonly solved using primal-dual interior-point methods [22].Another form of basis pursuit denoising is solved using the least absolute shrinkage and selection operator (LASSO), and it is defined as follows:(23)
m
i
n
x
1
2
y
-
A
x
2
2
s.t.
x
1
<
τ
,where τ is a nonnegative real parameter. The convex optimization methods usually require high computational complexity and high numerical precision.When the noise is unbounded, one may apply the convex program based on Dantzig selector (it is assumed that the noise variance isσ
2 per measurement, i.e., the total variance is M
σ
2): (24)
min
⁡
x
1
s.t.
A
∗
y
-
A
x
∞
≤
2
log
⁡
N
σ
,which (for enough measurements) reconstructs a signal with the error bound: (25)
x
-
x
^
2
≤
σ
2
log
⁡
N
K
.The norm l
∞ is infinite norm called also supreme (maximum) norm.Besides thel
1-norm minimization, there exist some approaches using the l
p-norm minimization, with 0
<
p
<
1:(26)
m
i
n
x
p
s.t.
y
-
A
x
2
≤
ε
,or using l
2-norm minimization, in which case the solution is not rigorously sparse enough [23]:(27)
m
i
n
x
2
s.t.
y
-
A
x
2
≤
ε
.


## 2.1. Sparsity and Compressibility

Reducing the sampling rate using CS is possible for the case of sparse signals that can be represented by a small number of significant coefficients in an appropriate transform basis. A signal havingK nonzero coefficients is called K-sparse. Assume that signal s exhibits sparsity in certain orthonormal basis Ψ defined by the basis vectors {
ψ
1
,
ψ
2
,
…
,
ψ
N
}. The signal s can be represented using its sparse transform domain vector x as follows:(1)
s
n
=
∑
i
=
1
N
x
i
ψ
i
n
.In matrix notation, the previous relation can be written as(2)
s
=
Ψ
x
.Commonly, the sparsity is measured using the l
0-norm, which represents the cardinality of the support of x: (3)
x
0
=
c
a
r
d
supp
⁡
x
=
K
.In real applications, the signals are usually not strictly sparse but only approximately sparse. Therefore, instead of being sparse, these signals are often called compressible, meaning that the amplitudes of coefficients decrease rapidly when arranged in descending order. For instance, if we consider coefficients x
1
≥
x
2
≥
⋯
≥
x
n, then the magnitude decays with a power law if there exist constants C
1 and q
>
0 satisfying [19](4)
x
i
≤
C
1
i
-
q
,where larger q means faster decay and consequently more compressible signal. The signal compressibility can be quantified using the minimal error between the original and sparsified signal (obtained by keeping only K largest coefficients):(5)
σ
K
x
p
=
min
⁡
x
-
x
K
p
.


## 2.2. Conditions on the CS Matrix: Null Space Property, Restricted Isometry Property, and Incoherence

Instead of acquiring a full set of signal samples of lengthN, in CS scenario, we deal with a quite reduced set of measurements y of length M, where M
<
N. The measurement procedure can be modeled by projections of the signal s onto vectors {
ϕ
1
,
ϕ
2
,
…
,
ϕ
M
} constituting the measurement matrix Φ:(6)
y
=
Φ
s
.Using the sparse transform domain representation of vectors given by (2), we have(7)
y
=
Φ
Ψ
x
=
A
x
,where A will be referred to as CS matrix.In order to define some requirements for the CS matrixA, which are important for successful signal reconstruction, let us introduce the null space of matrix. The null space of CS matrix A contains all vectors x that are mapped to 0:(8)
N
A
=
x
:
A
x
=
0
.In order to provide a unique solution, it is necessary to provide the notion that twoK-sparse vectors x and x
′ do not result in the same measurement vector. In other words, their difference should not be part of the null space of CS matrix A:(9)
A
x
-
x
′
≠
0
.Since the difference between twoK-sparse vectors is at most 2K-sparse, then aK-sparse vector x is uniquely defined if null space of A contains no 2K-sparse vectors. This corresponds to the condition that any 2K columns of A are linearly independent; that is, (10)
s
p
a
r
k
A
>
2
K
,and since s
p
a
r
k
A
∈
[
2
,
M
+
1
], we obtain a lower bound on the number of measurements:(11)
M
>
2
K
.In the case of strictly sparse signals, the spark can provide reliable information about the exact reconstruction possibility. However, in the case of approximately sparse signals, this condition is not sufficient and does not guarantee stable recovery. Hence, there is another property called null space property that measures the concentration of the null space of matrix A. The null space property is satisfied if there is a constant γ
∈
(
0,1
) such that [19](12)
a
S
1
=
γ
a
S
c
1
∀
a
∈
N
A
,for all sets S
⊂
{
1
,
…
,
N
} with cardinalityK and their complements S
c
=
{
1
,
…
,
N
}
∖
S. If the null space property is satisfied, then a strictly K-sparse signal can be perfectly reconstructed by using l
1-minimization. For approximately K-sparse signals, an upper bound of the l
1-minimization error can be defined as follows [20]:(13)
x
-
x
^
1
≤
2
1
+
γ
1
-
γ
σ
K
x
1
,where σ
K
(
x
)
1 is defined in (5) for p
=
1 as the minimal error induced by the best K-sparse approximation.The null space property is necessary and sufficient for establishing guarantees for recovery. A stronger condition is required in the presence of noise (and approximately sparse signals). Hence, in [8], the restricted isometry property (RIP) of CS matrix A has been introduced. The CS matrix A satisfies the RIP property with constant δ
K
∈
0,1 if(14)
1
-
δ
K
x
2
2
≤
A
x
2
2
≤
1
+
δ
K
x
2
2
,for every K-sparse vector x. This property shows how well the distances are preserved by a certain linear transformation. We might now say that if the RIP is satisfied for 2K with δ
2
K
<
1, then there are no twoK-sparse vectors x that can correspond to the same measurement vector y.Finally, the incoherence condition mentioned before, which is also related to the RIP of matrixA, refers to the incoherence of the projection basis Φ and the sparsifying basis Ψ. The mutual coherence can be simply defined by using the combined CS matrix A as follows [21]:(15)
μ
A
=
m
a
x
i
≠
j
,
1
≤
i
,
j
≤
N
A
i
,
A
j
A
i
2
A
j
2
.The mutual coherence is related to the restricted isometry constant using the following bound [22]:(16)
δ
K
=
K
-
1
μ
.


## 2.3. Signal Recovery Using Minimization Approach

The signal recovery problem is defined as the reconstruction of vectorx from the measurements y
=
A
x. This problem can be generally seen as a problem of solving an underdetermined set of linear equations. However, in the circumstances when x is sparse, the problem can be reduced to the following minimization:(17)
m
i
n
x
0
s.t.
y
=
A
x
.The l
0-minimization requires an exhaustive search over all N
K possible sparse combinations, which is computationally intractable. Hence, the l
0-minimization is replaced by convex l
1-minimization, which will provide the sparse result with high probability if the measurement matrix satisfies the previous conditions. The l
1-minimization problem is defined as follows:(18)
min
⁡
x
1
s.t.
y
=
A
x
,and it has been known as the basis pursuit.In the situation when the measurements are corrupted by the noise of levele
:
y
=
Φ
Ψ
x
+
e
=
A
x
+
e and e
2
≤
ε, the reconstruction problem can be defined in a form:(19)
min
⁡
x
1
s.t.
y
-
A
x
2
≤
ε
, called basis pursuit denoising. The error bound for the solution of (19), whereA satisfies the RIP of order 2
K with δ
2
K
<
2
-
1 and y
=
A
x
+
e, is given by(20)
x
-
x
^
2
≤
C
0
σ
K
x
1
+
C
1
2
ε
1
+
2
δ
2
K
,where the constants C
0 and C
1 are defined as [19] (21)
C
0
=
2
1
-
1
-
2
δ
2
K
1
-
1
+
2
δ
2
K
,
C
1
=
2
1
1
-
1
+
2
δ
2
K
.For a particular regularization parameter λ
>
0, the minimization problem (19) can be defined using the unconstrained version as follows:(22)
m
i
n
x
1
2
y
-
A
x
2
2
+
λ
x
1
,which is known as the Lagrangian form of the basis pursuit denoising. These algorithms are commonly solved using primal-dual interior-point methods [22].Another form of basis pursuit denoising is solved using the least absolute shrinkage and selection operator (LASSO), and it is defined as follows:(23)
m
i
n
x
1
2
y
-
A
x
2
2
s.t.
x
1
<
τ
,where τ is a nonnegative real parameter. The convex optimization methods usually require high computational complexity and high numerical precision.When the noise is unbounded, one may apply the convex program based on Dantzig selector (it is assumed that the noise variance isσ
2 per measurement, i.e., the total variance is M
σ
2): (24)
min
⁡
x
1
s.t.
A
∗
y
-
A
x
∞
≤
2
log
⁡
N
σ
,which (for enough measurements) reconstructs a signal with the error bound: (25)
x
-
x
^
2
≤
σ
2
log
⁡
N
K
.The norm l
∞ is infinite norm called also supreme (maximum) norm.Besides thel
1-norm minimization, there exist some approaches using the l
p-norm minimization, with 0
<
p
<
1:(26)
m
i
n
x
p
s.t.
y
-
A
x
2
≤
ε
,or using l
2-norm minimization, in which case the solution is not rigorously sparse enough [23]:(27)
m
i
n
x
2
s.t.
y
-
A
x
2
≤
ε
.


## 3. Review of Some Signal Reconstruction Algorithms

Thel
1-minimization problems in CS signal reconstruction are usually solved using the convex optimization methods. In addition, there exist greedy methods for sparse signal recovery which allow faster computation compared to l
1-minimization. Greedy algorithms can be divided into two major groups: greedy pursuit methods and thresholding-based methods. In practical applications, the widely used ones are the orthogonal matching pursuit (OMP) and compressive sampling matching pursuit (CoSaMP) from the group of greedy pursuit methods, while from the thresholding group the iterative hard thresholding (IHT) is commonly used due to its simplicity, although it may not be always efficient in providing an exact solution. Some of these algorithms are discussed in detail in this section.


### 3.1. Matching Pursuit

The matching pursuit algorithm has been known for its simplicity and was first introduced in [14]. This is the first algorithm from the class of iterative greedy methods that decomposes a signal into a linear set of basis functions. Through the iterations, this algorithm chooses in a greedy manner the basis functions that best match the signal. Also, in each iteration, the algorithm removes the signal component having the form of the selected basis function and obtains the residual. This procedure is repeated until the norm of the residual becomes lower than a certain predefined threshold value (halting criterion) (Algorithm 1).Algorithm 1:Matching pursuit.
Input:
▸ Measurement matrixA(of size M
×
N)
▸ Measurement vectory(of length M)
Output:
▸ A signal estimatex
^ (of length N)
(1) x
^
=
0, r
0
←
y, i
=
0
(2)while halting criterion falsedo
(3) i
=
i
+
1
(4) ω
i
←
arg
⁡
max
j
=
1
,
…
,
N
⁡
r
i
-
1
,
A
j    (max correlation column)
(5)  x
^
i
←
r
i
-
1
,
ω
i            (new signal estimate)
(6)  r
^
i
←
r
i
-
1
-
ω
i
x
^
i             (update residual)
(7)end while
(8)return
x
^
←
x
^
iThe matching pursuit algorithm however has a slow convergence property and generally does not provide efficiently sparse results.


### 3.2. Orthogonal Matching Pursuit

The orthogonal matching pursuit (OMP) has been introduced [15] as an improved version to overcome the limitations of the matching pursuit algorithm. OMP is based on principle of orthogonalization. It computes the inner product of the residue and the measurement matrix and then selects the index of the maximum correlation column and extracts this column (in each iteration). The extracted columns are included into the selected set of atoms. Then, the OMP performs orthogonal projection over the subspace of previously selected atoms, providing a new approximation vector used to update the residual. Here, the residual is always orthogonal to the columns of the CS matrix, so there will be no columns selected twice and the set of selected columns is increased through the iterations. OMP provides better asymptotic convergence compared to the previous matching pursuit version (Algorithm 2).Algorithm 2:Orthogonal matching pursuit.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA:  A
=
Φ
Ψ
▸ Measurement vectory
Output:
▸ A signal estimatex
^
←
x
K
▸ An approximation to the measurementsyby a
K
▸ A residualr
K
=
y
-
a
K.
▸ A setΩ
K with positions of non-zero elements of x
^.
(1) r
0
←
y, Ω
0
←
⌀, Θ
0
←
[
]
(2) for  i
=
1
,
…
,
K
(3) ω
i
←
arg
⁡
max
j
=
1
,
…
,
N
⁡
r
i
-
1
,
A
j     (max correlation column)
(4) Ω
i
←
Ω
i
-
1
∪
ω
i            (Update set of indices)
(5) Θ
i
←
[
Θ
i
-
1
A
ω
i
]            (Update set of atoms)
(6) x
i
=
arg
⁡
min
x
⁡
r
i
-
1
-
Θ
i
x
2
2
(7) a
i
←
Θ
i
x
i               (new approximation)
(8) r
i
←
y
-
a
i               (update residual)
(9)end for
(10) return
x
K, a
K, r
K, Ω
K


### 3.3. CoSaMP

Compressive sampling matching pursuit (CoSaMP) is an extended version of the OMP algorithm [16].In each iteration, the residual vectorr is correlated with the columns of the CS matrix A, forming a “proxy signal” z. Then, the algorithm selects 2K columns of A corresponding to the 2K largest absolute values of z, where K defines the signal sparsity (number of nonzero components). Namely, all but the largest 2K elements of z are set to zero and the 2K support Ω is obtained by finding the positions of nonzero elements. The indices of the selected columns (2K in total) are then added to the current estimate of the support of the unknown vector. After solving the least squares, a 3K-sparse estimate of the unknown vector is obtained. Then, the 3K-sparse vector is pruned to obtain the K-sparse estimate of the unknown vector (by setting all but the K largest elements to zero). Thus, the estimate of the unknown vector remains K-sparse and all columns that do not correspond to the true signal components are removed, which is an improvement over the OMP. Namely, if OMP selects in some iteration a column that does not correspond to the true signal component, the index of this column will remain in the final signal support and cannot be removed (Algorithm 3).Algorithm 3:CoSaMP.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA: A
=
Φ
Ψ
▸ Measurement vectory
▸K being the signal sparsity
▸ Halting criterion
Output:
▸K-sparse approximation x
^ of the target signal
(1) x
0
←
0, r
←
y, i
←
0
(2) while  halting condition falsedo
(3) i
←
i
+
1
(4) z
←
A
∗
r                   (signal proxy)
(5) Ω
←
s
u
p
p
(
z
2
K
)          (support of best 2K-sparse approximation)
(6) T
←
Ω
∪
s
u
p
p
(
x
i
-
1
)            (merge supports)
(7) x
-
=
arg
⁡
min
x
~
:
supp
(
x
~
)
=
T
⁡
A
x
~
-
y
2
2        (solve least squares)
(8) x
i
←
x
-
K             (Prune: best
K
-sparse approximation)
(9) r
←
y
-
A
x
i                (update current sample)
(10)end while
(11) x
^
←
x
i
(12) return
x
^The CoSaMP can be observed through five crucial steps:(i)
Identification (Line 5). It finds the largest 2s components of the signal proxy.
(ii)
Support Merge (Line 6). It merges the support of the signal proxy with the support of the solution from the previous iteration.
(iii)
Estimation (Line 7). It estimates a solution via least squares where the solution lies within a support T.
(iv)
Pruning (Line 8). It takes the solution estimate and compresses it to the required support.
(v)
Sample Update (Line 9). It updates the “sample,” meaning the residual as the part of signal that has not been approximated.


### 3.4. Iterative Hard Thresholding Algorithm

Another group of algorithms for signal reconstruction from a small set of measurements is based on iterative thresholding. Generally, these algorithms are composed of two main steps. The first one is the optimization of the least squares term, which is done by solving the optimization problem withoutl
1-minimization. The other one is the decreasing of the l
1-norm, which is done by applying the thresholding operator to the magnitude of entries in x. In each iteration, the sparse vector x
i is estimated by the previous version x
i
-
1 using negative gradient of the objective function defined as (28)
J
x
=
1
2
y
-
A
x
2
2
,while the negative gradient is then(29)
-
∇
J
x
=
A
T
y
-
A
x
.Generally, the obtained estimate x
i is sparse, and we need to set all but theK largest components to zero using the thresholding operator. Here, we distinguish two types of thresholding. Hard thresholding [24, 25] sets all but theK largest magnitude values to zero, where the thresholding operator can be written as(30)
H
K
x
=
x
k
,
f
o
r
x
k
>
ξ
0
,
otherwise
,where ξ is theK largest component ofx. The algorithm is summarized in Algorithm 4.Algorithm 4:Iterative hard thresholding.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA: A
=
Φ
Ψ
▸ Measurement vectory
▸K being the signal sparsity
Output:
▸K-sparse approximation x
^
(1) x
0
←
0
(2) whilestopping criterion do
(3) i
←
i
+
1
(4) x
i
←
H
K
x
i
-
1
+
A
T
(
y
-
A
x
i
-
1
)
(5) end while
(6) x
^
←
x
i
(7) return
x
^The stopping criterion for IHT can be a fixed number of iterations or the algorithm terminates when the sparse vector does not change much between consecutive iterations.The soft-thresholding function can be defined as(31)
S
τ
x
k
=
x
k
-
τ
,
x
>
τ
0
,
x
k
<
τ
x
k
+
τ
,
x
<
-
τand it is applied to each element of x.


### 3.5. Iterative Soft Thresholding for Solving LASSO

Let us consider the minimization of the function:(32)
J
x
=
y
-
A
x
2
+
λ
x
1
,as given by the LASSO minimization problem (22). One of the algorithms that has been used for solving (32) is the iterated soft-thresholding algorithm (ISTA), also called the thresholded Landweber algorithm. In order to provide an iterative procedure for solving the considered minimization problem, we use the minimization-maximization [26] approach to minimize J
(
x
). Hence, we should create a function G
k
(
x
) that is equal to J
(
x
) at x
k; otherwise, it upper-bounds J
(
x
). Minimizing a majorizer G
k
(
x
) is easier and avoids solving a system of equations. Hence [27], (33)
G
k
x
=
J
x
+
nonnegative  function  of
x
,
G
k
x
=
J
x
+
x
-
x
k
T
α
I
-
A
T
A
x
-
x
k
=
y
-
A
x
2
+
λ
x
1
+
x
-
x
k
T
α
I
-
A
T
A
x
-
x
k
,and α is such that the added function is nonnegative, meaning that α must be equal to or greater than the maximum eigenvalue of A
T
A: (34)
α
>
max
⁡
e
i
g
A
T
A
.In order to minimize G
k
(
x
), the gradient of G
k
x should be zero:(35)
∂
G
k
x
∂
x
T
=
0
⟹
-
2
A
T
y
+
2
A
T
A
x
+
λ
sign
⁡
x
+
2
α
I
-
A
T
A
x
-
x
k
=
0
.The solution is(36)
x
+
λ
2
α
sign
⁡
x
=
1
α
Α
T
y
-
A
x
k
+
x
k
.The solution of the linear equation(37)
a
+
λ
2
α
sign
⁡
a
=
bis obtained using the soft-thresholding approach:(38)
a
=
S
λ
b
=
b
-
λ
,
b
>
λ
0
,
b
<
λ
b
+
λ
,
b
<
-
λ
.Hence, we may write the soft-thresholding-based solution of (36) as follows:(39)
x
=
S
λ
/
2
α
1
α
Α
T
y
-
A
x
k
+
x
k
.In terms of iteration k
+
1, the previous relation can be rewritten as(40)
x
k
+
1
=
S
λ
/
2
α
1
α
Α
T
y
-
A
x
k
+
x
k
.


### 3.6. Automated Threshold Based Iterative Solution

This algorithm (see [12, 13]) starts from the assumption that the missing samples, modeled by zero values in the signal domain, produce noise in the transform domain. Namely, the zero values at the positions of missing samples can be observed as a result of adding noise to these samples, where the noise at the unavailable positions has values of original signal samples with opposite sign. For instance, let us observe the signals that are sparse in the Fourier transform domain. The variance of the mentioned noise when observed in the discrete Fourier transform (DFT) domain can be calculated as (41)
σ
M
S
2
⟵
M
N
-
M
N
-
1
∑
i
=
1
M
y
i
2
M
,whereM is the number of available andN is the total number of samples and y is a measurement vector. Consequently, using the variance of noise, it is possible to define a threshold to separate signal components from the nonsignal components in the DFT domain.For a desired probabilityP
(
T
), the threshold is derived as (42)
T
=
-
σ
M
S
2
log
⁡
1
-
P
T
1
/
N
.The automated threshold based algorithm, in each iteration, detects a certain number of DFT components above the threshold. A set of positions k
i corresponding to the detected components is selected and the contribution of these components is removed from the signal’s DFT. This will reveal the remaining components that are below the noise level. Further, it is necessary to update the noise variance and threshold value for the new iteration. Since the algorithm detects a set of components in each iteration, it usually needs just a couple of iterations to recover the entire signal. In the case that all signal components are above the noise level in DFT, the component detection and reconstruction are achieved in single iteration.


### 3.7. Adaptive Gradient-Based Algorithm

This algorithm belongs to the group of convex minimization approaches [11]. Unlike other convex minimization approaches, this method starts from some initial values of unavailable samples (initial state) which are changed through the iterations in a way to constantly improve the concentration in sparsity domain. In general, it does not require the signal to be strictly sparse in certain transform domain, which is an advantage over other methods. Particularly, the missing samples in the signal domain can be considered as zero values. In each iteration, the missing samples are changed for +
Δ and for -
Δ. Then, for both changes, the concentration is measured as l
1-norm of the transform domain vectors x
+ (for +
Δ change) and x
- (for -
Δ change), while the gradient is determined using their difference (lines 8, 9, and 10). Finally, the gradient is used to update the values of missing samples. Here, it is important to note that each sample is observed separately in this algorithm and one iteration is finished when all samples are processed.When the algorithm reaches the vicinity of the sparsity measure minimum, the gradient changes direction for almost 180 degrees (line 16), meaning that the stepΔ needs to be decreased (line 17).The precision of the result in this iterative algorithm is estimated based on the change of the result in the last iteration (lines 18 and 19).


## 3.1. Matching Pursuit

The matching pursuit algorithm has been known for its simplicity and was first introduced in [14]. This is the first algorithm from the class of iterative greedy methods that decomposes a signal into a linear set of basis functions. Through the iterations, this algorithm chooses in a greedy manner the basis functions that best match the signal. Also, in each iteration, the algorithm removes the signal component having the form of the selected basis function and obtains the residual. This procedure is repeated until the norm of the residual becomes lower than a certain predefined threshold value (halting criterion) (Algorithm 1).Algorithm 1:Matching pursuit.
Input:
▸ Measurement matrixA(of size M
×
N)
▸ Measurement vectory(of length M)
Output:
▸ A signal estimatex
^ (of length N)
(1) x
^
=
0, r
0
←
y, i
=
0
(2)while halting criterion falsedo
(3) i
=
i
+
1
(4) ω
i
←
arg
⁡
max
j
=
1
,
…
,
N
⁡
r
i
-
1
,
A
j    (max correlation column)
(5)  x
^
i
←
r
i
-
1
,
ω
i            (new signal estimate)
(6)  r
^
i
←
r
i
-
1
-
ω
i
x
^
i             (update residual)
(7)end while
(8)return
x
^
←
x
^
iThe matching pursuit algorithm however has a slow convergence property and generally does not provide efficiently sparse results.


## 3.2. Orthogonal Matching Pursuit

The orthogonal matching pursuit (OMP) has been introduced [15] as an improved version to overcome the limitations of the matching pursuit algorithm. OMP is based on principle of orthogonalization. It computes the inner product of the residue and the measurement matrix and then selects the index of the maximum correlation column and extracts this column (in each iteration). The extracted columns are included into the selected set of atoms. Then, the OMP performs orthogonal projection over the subspace of previously selected atoms, providing a new approximation vector used to update the residual. Here, the residual is always orthogonal to the columns of the CS matrix, so there will be no columns selected twice and the set of selected columns is increased through the iterations. OMP provides better asymptotic convergence compared to the previous matching pursuit version (Algorithm 2).Algorithm 2:Orthogonal matching pursuit.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA:  A
=
Φ
Ψ
▸ Measurement vectory
Output:
▸ A signal estimatex
^
←
x
K
▸ An approximation to the measurementsyby a
K
▸ A residualr
K
=
y
-
a
K.
▸ A setΩ
K with positions of non-zero elements of x
^.
(1) r
0
←
y, Ω
0
←
⌀, Θ
0
←
[
]
(2) for  i
=
1
,
…
,
K
(3) ω
i
←
arg
⁡
max
j
=
1
,
…
,
N
⁡
r
i
-
1
,
A
j     (max correlation column)
(4) Ω
i
←
Ω
i
-
1
∪
ω
i            (Update set of indices)
(5) Θ
i
←
[
Θ
i
-
1
A
ω
i
]            (Update set of atoms)
(6) x
i
=
arg
⁡
min
x
⁡
r
i
-
1
-
Θ
i
x
2
2
(7) a
i
←
Θ
i
x
i               (new approximation)
(8) r
i
←
y
-
a
i               (update residual)
(9)end for
(10) return
x
K, a
K, r
K, Ω
K


## 3.3. CoSaMP

Compressive sampling matching pursuit (CoSaMP) is an extended version of the OMP algorithm [16].In each iteration, the residual vectorr is correlated with the columns of the CS matrix A, forming a “proxy signal” z. Then, the algorithm selects 2K columns of A corresponding to the 2K largest absolute values of z, where K defines the signal sparsity (number of nonzero components). Namely, all but the largest 2K elements of z are set to zero and the 2K support Ω is obtained by finding the positions of nonzero elements. The indices of the selected columns (2K in total) are then added to the current estimate of the support of the unknown vector. After solving the least squares, a 3K-sparse estimate of the unknown vector is obtained. Then, the 3K-sparse vector is pruned to obtain the K-sparse estimate of the unknown vector (by setting all but the K largest elements to zero). Thus, the estimate of the unknown vector remains K-sparse and all columns that do not correspond to the true signal components are removed, which is an improvement over the OMP. Namely, if OMP selects in some iteration a column that does not correspond to the true signal component, the index of this column will remain in the final signal support and cannot be removed (Algorithm 3).Algorithm 3:CoSaMP.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA: A
=
Φ
Ψ
▸ Measurement vectory
▸K being the signal sparsity
▸ Halting criterion
Output:
▸K-sparse approximation x
^ of the target signal
(1) x
0
←
0, r
←
y, i
←
0
(2) while  halting condition falsedo
(3) i
←
i
+
1
(4) z
←
A
∗
r                   (signal proxy)
(5) Ω
←
s
u
p
p
(
z
2
K
)          (support of best 2K-sparse approximation)
(6) T
←
Ω
∪
s
u
p
p
(
x
i
-
1
)            (merge supports)
(7) x
-
=
arg
⁡
min
x
~
:
supp
(
x
~
)
=
T
⁡
A
x
~
-
y
2
2        (solve least squares)
(8) x
i
←
x
-
K             (Prune: best
K
-sparse approximation)
(9) r
←
y
-
A
x
i                (update current sample)
(10)end while
(11) x
^
←
x
i
(12) return
x
^The CoSaMP can be observed through five crucial steps:(i)
Identification (Line 5). It finds the largest 2s components of the signal proxy.
(ii)
Support Merge (Line 6). It merges the support of the signal proxy with the support of the solution from the previous iteration.
(iii)
Estimation (Line 7). It estimates a solution via least squares where the solution lies within a support T.
(iv)
Pruning (Line 8). It takes the solution estimate and compresses it to the required support.
(v)
Sample Update (Line 9). It updates the “sample,” meaning the residual as the part of signal that has not been approximated.


## 3.4. Iterative Hard Thresholding Algorithm

Another group of algorithms for signal reconstruction from a small set of measurements is based on iterative thresholding. Generally, these algorithms are composed of two main steps. The first one is the optimization of the least squares term, which is done by solving the optimization problem withoutl
1-minimization. The other one is the decreasing of the l
1-norm, which is done by applying the thresholding operator to the magnitude of entries in x. In each iteration, the sparse vector x
i is estimated by the previous version x
i
-
1 using negative gradient of the objective function defined as (28)
J
x
=
1
2
y
-
A
x
2
2
,while the negative gradient is then(29)
-
∇
J
x
=
A
T
y
-
A
x
.Generally, the obtained estimate x
i is sparse, and we need to set all but theK largest components to zero using the thresholding operator. Here, we distinguish two types of thresholding. Hard thresholding [24, 25] sets all but theK largest magnitude values to zero, where the thresholding operator can be written as(30)
H
K
x
=
x
k
,
f
o
r
x
k
>
ξ
0
,
otherwise
,where ξ is theK largest component ofx. The algorithm is summarized in Algorithm 4.Algorithm 4:Iterative hard thresholding.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸ CS matrixA: A
=
Φ
Ψ
▸ Measurement vectory
▸K being the signal sparsity
Output:
▸K-sparse approximation x
^
(1) x
0
←
0
(2) whilestopping criterion do
(3) i
←
i
+
1
(4) x
i
←
H
K
x
i
-
1
+
A
T
(
y
-
A
x
i
-
1
)
(5) end while
(6) x
^
←
x
i
(7) return
x
^The stopping criterion for IHT can be a fixed number of iterations or the algorithm terminates when the sparse vector does not change much between consecutive iterations.The soft-thresholding function can be defined as(31)
S
τ
x
k
=
x
k
-
τ
,
x
>
τ
0
,
x
k
<
τ
x
k
+
τ
,
x
<
-
τand it is applied to each element of x.


## 3.5. Iterative Soft Thresholding for Solving LASSO

Let us consider the minimization of the function:(32)
J
x
=
y
-
A
x
2
+
λ
x
1
,as given by the LASSO minimization problem (22). One of the algorithms that has been used for solving (32) is the iterated soft-thresholding algorithm (ISTA), also called the thresholded Landweber algorithm. In order to provide an iterative procedure for solving the considered minimization problem, we use the minimization-maximization [26] approach to minimize J
(
x
). Hence, we should create a function G
k
(
x
) that is equal to J
(
x
) at x
k; otherwise, it upper-bounds J
(
x
). Minimizing a majorizer G
k
(
x
) is easier and avoids solving a system of equations. Hence [27], (33)
G
k
x
=
J
x
+
nonnegative  function  of
x
,
G
k
x
=
J
x
+
x
-
x
k
T
α
I
-
A
T
A
x
-
x
k
=
y
-
A
x
2
+
λ
x
1
+
x
-
x
k
T
α
I
-
A
T
A
x
-
x
k
,and α is such that the added function is nonnegative, meaning that α must be equal to or greater than the maximum eigenvalue of A
T
A: (34)
α
>
max
⁡
e
i
g
A
T
A
.In order to minimize G
k
(
x
), the gradient of G
k
x should be zero:(35)
∂
G
k
x
∂
x
T
=
0
⟹
-
2
A
T
y
+
2
A
T
A
x
+
λ
sign
⁡
x
+
2
α
I
-
A
T
A
x
-
x
k
=
0
.The solution is(36)
x
+
λ
2
α
sign
⁡
x
=
1
α
Α
T
y
-
A
x
k
+
x
k
.The solution of the linear equation(37)
a
+
λ
2
α
sign
⁡
a
=
bis obtained using the soft-thresholding approach:(38)
a
=
S
λ
b
=
b
-
λ
,
b
>
λ
0
,
b
<
λ
b
+
λ
,
b
<
-
λ
.Hence, we may write the soft-thresholding-based solution of (36) as follows:(39)
x
=
S
λ
/
2
α
1
α
Α
T
y
-
A
x
k
+
x
k
.In terms of iteration k
+
1, the previous relation can be rewritten as(40)
x
k
+
1
=
S
λ
/
2
α
1
α
Α
T
y
-
A
x
k
+
x
k
.


## 3.6. Automated Threshold Based Iterative Solution

This algorithm (see [12, 13]) starts from the assumption that the missing samples, modeled by zero values in the signal domain, produce noise in the transform domain. Namely, the zero values at the positions of missing samples can be observed as a result of adding noise to these samples, where the noise at the unavailable positions has values of original signal samples with opposite sign. For instance, let us observe the signals that are sparse in the Fourier transform domain. The variance of the mentioned noise when observed in the discrete Fourier transform (DFT) domain can be calculated as (41)
σ
M
S
2
⟵
M
N
-
M
N
-
1
∑
i
=
1
M
y
i
2
M
,whereM is the number of available andN is the total number of samples and y is a measurement vector. Consequently, using the variance of noise, it is possible to define a threshold to separate signal components from the nonsignal components in the DFT domain.For a desired probabilityP
(
T
), the threshold is derived as (42)
T
=
-
σ
M
S
2
log
⁡
1
-
P
T
1
/
N
.The automated threshold based algorithm, in each iteration, detects a certain number of DFT components above the threshold. A set of positions k
i corresponding to the detected components is selected and the contribution of these components is removed from the signal’s DFT. This will reveal the remaining components that are below the noise level. Further, it is necessary to update the noise variance and threshold value for the new iteration. Since the algorithm detects a set of components in each iteration, it usually needs just a couple of iterations to recover the entire signal. In the case that all signal components are above the noise level in DFT, the component detection and reconstruction are achieved in single iteration.


## 3.7. Adaptive Gradient-Based Algorithm

This algorithm belongs to the group of convex minimization approaches [11]. Unlike other convex minimization approaches, this method starts from some initial values of unavailable samples (initial state) which are changed through the iterations in a way to constantly improve the concentration in sparsity domain. In general, it does not require the signal to be strictly sparse in certain transform domain, which is an advantage over other methods. Particularly, the missing samples in the signal domain can be considered as zero values. In each iteration, the missing samples are changed for +
Δ and for -
Δ. Then, for both changes, the concentration is measured as l
1-norm of the transform domain vectors x
+ (for +
Δ change) and x
- (for -
Δ change), while the gradient is determined using their difference (lines 8, 9, and 10). Finally, the gradient is used to update the values of missing samples. Here, it is important to note that each sample is observed separately in this algorithm and one iteration is finished when all samples are processed.When the algorithm reaches the vicinity of the sparsity measure minimum, the gradient changes direction for almost 180 degrees (line 16), meaning that the stepΔ needs to be decreased (line 17).The precision of the result in this iterative algorithm is estimated based on the change of the result in the last iteration (lines 18 and 19).


## 4. Exploring Different Transform Domains for CS Signal Reconstruction

### 4.1. CS in the Standard Fourier Transform Domain

As the simplest case of CS scenario, we will assume a multicomponent signal that consists ofK sinusoids that are sparse in the DFT domain:(43)
s
n
=
∑
i
=
1
K
A
i
e
j
2
π
k
i
n
/
N
,where the sparsity level is K
≪
N, where N is total signal length, while A
i and k
i denote amplitudes and frequencies of signal components, respectively. Sinces is sparse in the DFT domain, then we can write(44)
s
=
Ψ
F
=
F
-
1
F
,where F is a vector of DFT coefficients where at mostK coefficients are nonzero, while F
-
1 is the inverse Fourier transform matrix of size N
×
N.Ifs is a signal in compressive sensing application, then only the random measurements y
⊂
s are available and these are defined by the set ofM positions {
n
1
,
n
2
,
n
M
}. Therefore, the measurement process can be modeled by matrix Φ:(45)
y
=
Φ
F
-
1
F
=
A
F
.The CS matrix A represents a partial random inverse Fourier transform matrix obtained by omitting rows from F
-
1 that corresponds to unavailable samples positions:(46)
A
=
ψ
0
n
1
ψ
1
n
1
⋯
ψ
N
-
1
n
1
ψ
0
n
2
ψ
1
n
2
⋯
ψ
N
-
1
n
2
⋮
⋮
⋱
⋮
ψ
0
n
M
ψ
1
n
M
⋯
ψ
N
-
1
n
M
,with the Fourier basis functions: (47)
ψ
k
n
=
e
j
2
π
k
n
/
N
.The CS problem formulation is now given as (48)
min
⁡
F
1
s.t.
y
=
A
F
.


### 4.2. CS in the Polynomial Fourier Transform Domain

In this section, we consider the possibility of CS reconstruction of polynomial phase signals [28]. Observe the multicomponent polynomial phase signal vectors, with elementss(n): (49)
s
n
=
∑
i
=
1
K
s
i
n
=
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
,where the polynomial coefficients are assumed to be bounded integers. The assumption is that the signal can be considered as sparse in the polynomial Fourier transform (PFT) domain. Hence, the discrete PFT form is given by(50)
X
k
1
i
,
…
,
k
p
i
=
∑
n
=
0
N
-
1
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
×
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!
e
-
j
2
π
/
N
n
k
1
. If we choose a set of parameters (
k
2
i
,
k
3
i
,
…
,
k
p
i
) that match the polynomial phase coefficients (
a
2
i
,
a
3
i
,
…
,
a
p
i
)
(51)
k
2
i
,
k
3
i
,
…
,
k
p
i
=
a
2
i
,
a
3
i
,
…
,
a
p
i
,then the ith signal component is demodulated and becomes a sinusoid: r
i
e
j
(
2
π
/
N
)
n
a
1
i which is dominant in the PFT. In other words, the spectrum is highly concentrated at k
=
a
1
i. Therefore, we might say that if (51) is satisfied then the PFT is compressible with the dominant ith component. Note that the sparsity (compressibility) in the PFT domain is observed with respect to the single demodulated component.In order to define the CS problem in the PFT domain, instead of the signals itself, we will consider a modified signal form obtained after the multiplication with the exponential term:(52)
ξ
n
=
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!
,given in the PFT definition (50). The new signal form is obtained as(53)
z
n
=
s
n
ξ
n
=
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!or in the vector form:(54)
z
=
s
ξ
.Then, the PFT definition given by (50) can be rewritten in the vector form as follows:(55)
X
=
F
z
,where F is the discrete Fourier transform matrix (N
×
N).For a chosen set of parameters(
k
2
i
,
k
3
i
,
…
,
k
p
i
) in ξ that is equal to the set (
a
2
i
,
a
3
i
,
…
,
a
n
i
) in s, X
=
X
i is characterized by one dominant sinusoidal component at the frequency a
1
i.


#### 4.2.1. CS Scenario

Now, assume thatz is incomplete in the compressive sensing sense, and instead of z we are dealing with M available measurements defined by vector y
z, and (56)
y
z
=
A
z
X
,where A
z is the partial random Fourier matrix obtained by omitting rows from F that correspond to the unavailable samples. When (
k
2
i
,
k
3
i
,
…
,
k
p
i
)
=
(
a
2
i
,
a
3
i
,
…
,
a
p
i
), thenX can be observed as a demodulated version of the ith signal component X
i, having the dominant ith component in the spectrum with the support k
1
i. The rest of the components in spectrum are much lower than X
i and could be observed as noise. Hence, we may write the minimization problem in the form(57)
min
⁡
X
i
1
s.t.
y
z
-
A
z
X
i
2
<
T
,whereT is certain threshold. The dominant components can be detected using the threshold based algorithm (Algorithm 5) described in the previous section. Using an iterative procedure, one may change the values of parameters k
2
i
,
…
,
k
p
i between k
m
i
n and k
m
a
x, ∀
i. The algorithm will detect the support k
1
i when the set (
k
2
i
,
k
3
i
,
…
,
k
p
i
) matches the set (
a
2
i
,
a
3
i
,
…
,
a
p
i
); otherwise, no component support is detected. Hence, as a result of this phase, we have identified the sets of signal phase parameters: k
i
=
k
1
i
,
k
2
i
,
…
,
k
n
i
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).Algorithm 5:Automated threshold based iterative algorithm.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸N
M
=
n
1
,
n
2
,
…
,
n
M, M –number of available samples, N-original signal length
▸ MatrixA: A
=
Φ
Ψ (A is obtained as a partial random Fourier matrix: A = Φ
Ψ
=
Ψ
{
ϑ
},
ϑ
=
{
ψ
n
1
,
ψ
n
2
,
…
,
ψ
n
M
}
▸ Measurement vectory
=
Φ
f
Output:
(1) i
=
1, N
M
=
{
n
1
,
n
2
,
…
,
n
M
}
(2) p
=
⌀.
(3) P
←
0.99                  (Set desired probability
P)
(4) σ
M
S
2
←
M
(
(
N
-
M
)
/
(
N
-
1
)
)
∑
i
=
1
M
(
y
(
i
)
2
/
M
)      (Calculate variance)
(5) T
i
←
-
σ
M
S
2
log
⁡
(
1
-
P
(
T
)
1
/
N
)           (Calculate threshold)
(6) X
i
←
A
-
1
y                  (Calculate initialDFTof
y)
(7) k
i
←
arg
⁡
X
i
>
T
i
/
N            (Find comp. in
X
i  above the T
i)
(8) Update
p
←
p
∪
k
i
(9) A
c
s
i
←
A
(
k
i
)                   (CS matrix)
(1
0) X
~
i
←
A
c
s
i
T
A
c
s
i
-
1
A
c
s
i
T
y
(11) for
∀
p
∈
p:
(1
2) y
←
y
-
(
M
/
N
)
X
i
p
exp
⁡
(
j
2
π
p
N
M
/
N
);            (Update
y)
(13)Update
σ
M
S
2
←
(
M
(
N
-
M
)
/
(
N
-
1
)
)
∑
y
2
/
M
(14)If
∑
y
2
/
M
<
δ
break; Else
(15)Set
i
=
i
+
1, and go to (5).
(16)return
X
~
i,  k


#### 4.2.2. Reconstruction of Exact Components Amplitudes

The next phase is reconstruction of components amplitudes. Denote the set of measurements positions by(
n
1
,
n
2
,
…
,
n
M
), such that [
y
(
1
)
,
y
(
2
)
,
…
,
y
(
M
)
]
=
[
s
(
n
1
)
,
s
(
n
2
)
,
…
,
s
(
n
M
)
], while the detected sets of signal phase parameters are k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).In order to calculate the exact amplitudesr
1
,
r
2
,
…
,
r
K of signal components, we observe the set of equations in the form(58)
s
n
1
⋮
s
n
M
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
N
k
11
+
⋯
+
n
N
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
N
p
k
p
K
r
1
⋮
r
Kor in other words we have another system of equations given by (59)
y
=
A
R
,where R
=
r
1
,
…
,
r
K
T contains the desired K signal amplitudes. The matrix A of size N
×
K is based on the PFT: (60)
A
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
M
k
11
+
⋯
+
n
M
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
M
p
k
p
K
.The rows of A correspond to positions of measurements (
n
1
,
n
2
,
…
,
n
M
), and columns correspond to phase parameters k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
), for i
=
1
,
…
,
K. The solution of the observed problem can be obtained in the least squares sense as follows:(61)
X
=
A
∗
A
-
1
A
∗
y
.The resulting reconstructed signal is obtained as(62)
s
n
=
r
1
e
-
j
2
π
/
N
n
a
11
+
⋯
+
n
p
a
p
1
/
p
!
+
⋯
+
r
K
e
-
j
2
π
/
N
n
a
1
K
+
⋯
+
n
p
a
p
K
/
p
!
.


### 4.3. CS in the Hermite Transform Domain

The Hermite expansion usingN Hermite functions can be written in terms of the Hermite transform matrix. First, let us define the Hermite transform matrix H (of size N
×
N) [29, 30]:(63)
H
=
1
N
ψ
0
1
ψ
N
-
1
1
2
ψ
0
2
ψ
N
-
1
2
2
⋯
ψ
0
N
ψ
N
-
1
N
2
ψ
1
1
ψ
N
-
1
1
2
ψ
1
2
ψ
N
-
1
2
2
⋯
ψ
1
N
ψ
N
-
1
N
2
⋮
⋮
⋱
⋮
ψ
N
-
1
1
ψ
N
-
1
1
2
ψ
N
-
1
2
ψ
N
-
1
2
2
⋯
ψ
N
-
1
N
ψ
N
-
1
N
2
.Here, the pth order Hermite basis function is defined using the pth order Hermite polynomial as follows: (64)
ψ
p
n
=
σ
2
p
p
!
π
-
1
/
2
e
-
n
2
/
2
H
p
n
σ
,where σ is the scaling factor used to “stretch” or “compress” Hermite functions, in order to match the signal.The Hermite functions are usually calculated using the following fast recursive realization:(65)
ψ
0
n
=
1
π
4
e
-
n
2
/
2
σ
2
,
ψ
1
n
=
2
n
π
4
e
-
n
2
/
2
σ
2
,
ψ
i
n
=
n
2
i
Ψ
i
-
1
n
-
i
-
1
i
Ψ
i
-
2
n
,
∀
i
≥
2
.The Hermite transform coefficients can be calculated in the matrix form as follows: (66)
c
=
H
s
,where the vector of Hermite transform coefficients is c
=
c
0
,
c
1
,
…
,
c
N
-
1
T, while the vector of signal samples is s
=
[
s
(
1
)
,
s
(
2
)
,
…
,
s
(
N
)
]
T.Following the Gauss-Hermite approximation, the inverse matrixΨ
=
H
-
1 contains N Hermite functions and it is given by(67)
Ψ
=
ψ
1
1
ψ
2
1
⋯
ψ
N
1
ψ
1
2
ψ
2
2
⋯
ψ
N
2
⋮
⋮
⋱
⋮
ψ
1
N
ψ
2
N
⋯
ψ
N
N
.Now, in the context of CS, let us assume that a signal s isK-sparse in the Hermite transform domain, meaning that it can be represented by a small set of nonzero Hermite coefficientsc: (68)
s
n
=
∑
p
=
k
1
k
K
c
p
ψ
p
n
,where c
p is the pth element from the vector c of Hermite expansion coefficients. In the matrix form we may write (69)
s
=
Ψ
c
,such that most of the coefficients in c are zero values. Furthermore, assume that only the compressive sensing is done using random selection ofM signal values [
n
1
,
n
2
,
…
,
n
M
]. Then, the CS matrix A can be defined as random partial inverse Hermite transform matrix: (70)
A
=
ψ
0
n
1
ψ
1
n
1
⋯
ψ
N
-
1
n
1
ψ
0
n
2
ψ
1
n
2
⋯
ψ
N
-
1
n
2
⋮
⋮
⋱
⋮
ψ
0
n
M
ψ
1
n
M
⋯
ψ
N
-
1
n
M
.For a vector of available measurements vector y
=
(
s
(
n
1
)
,
s
(
n
2
)
,
…
,
s
(
n
M
)
), the linear system of equations (undetermined system ofM linear equations and N unknowns) can be written as(71)
y
=
A
c
.Since c has only K nonzero components, the reconstruction problem can be defined as(72)
m
i
n
c
1
s.t.
y
=
A
c
.If we identify the signal support in the Hermite transform domain by the set of indices [
k
1
,
k
2
,
…
,
k
K
], then the problem can be solved in the least squares sense as(73)
c
^
=
Θ
c
s
∗
Θ
c
s
-
1
Θ
c
s
∗
y
,where (
∗
) denotes the conjugate transpose operation, while(74)
Θ
c
s
=
ψ
k
1
n
1
ψ
k
1
n
2
…
ψ
k
1
n
M
ψ
k
2
n
1
ψ
k
2
n
2
⋯
ψ
k
2
n
M
⋮
⋮
⋱
⋮
ψ
k
K
n
1
ψ
k
K
n
2
⋯
ψ
k
K
n
M
,and c
^
=
c
k
1
,
c
k
2
,
…
,
c
k
K
T. The estimated signal s can be now obtained according to (68).


### 4.4. CS in the Time-Frequency Domain

Starting from the definition of the short-time Fourier transform(75)
STFT
n
,
k
=
∑
m
=
0
M
-
1
x
n
+
m
e
-
j
2
π
m
k
/
M
,we can define the following STFT vector and signal vector: (76)
S
T
F
T
M
n
=
STFT
n
,
0
,
…
,
STFT
n
,
M
-
1
T
,
x
n
=
x
n
,
x
n
+
1
,
…
,
x
n
+
M
-
1
T
,respectively. The STFT at an instantn can be now rewritten in the matrix form:(77)
S
T
F
T
M
n
=
F
M
x
n
,where F
M is the standard Fourier transform matrix of size M
×
M with elements(78)
F
M
m
,
k
=
e
-
j
2
π
k
m
/
M
.If we consider the nonoverlapping signal segments for the STFT calculation, then n takes the values from the set {
0
,
M
,
2
M
,
…
,
N
-
M
}, where M is the window length while N is the total signal length. The STFT vector, for n
∈
{
0
,
M
,
2
M
,
…
,
N
-
M
}, is given by(79)
S
T
F
T
=
S
T
F
T
M
0
S
T
F
T
M
M
⋮
S
T
F
T
M
N
-
M
.Consequently, (77) can be extended to the entire STFT for all considered values of n
∈
{
0
,
M
,
2
M
,
…
,
N
-
M
} as follows [31, 32]:(80)
S
T
F
T
M
0
S
T
F
T
M
M
⋮
S
T
F
T
M
N
-
M
=
F
M
0
M
⋯
0
M
0
M
F
M
⋯
0
M
⋯
⋯
⋯
⋯
0
M
0
M
⋯
F
M
x
0
x
M
⋮
x
N
-
M
,or in the compact form(81)
S
T
F
T
=
F
x
.In the case of time-varying windows with lengths {
M
1
,
M
2
,
…
,
M
Q
}, the previous system of equations becomes(82)
S
T
F
T
M
0
0
S
T
F
T
M
1
M
0
⋮
S
T
F
T
M
Q
N
-
M
Q
=
F
M
0
0
⋯
0
0
M
F
M
1
⋯
0
⋯
⋯
⋯
⋯
0
0
⋯
F
M
Q
x
0
x
M
⋮
x
N
-
M
.Note that the signal vector is x
=
[
x
(
0
)
T
,
x
(
M
)
T
,
…
,
x
(
N
-
M
)
T
]
T
=
[
x
(
0
)
,
x
(
1
)
,
…
,
x
(
N
-
1
)
]
T and can be expressed using the N-point inverse Fourier transform as follows:(83)
x
=
F
N
-
1
X
,whereX is a sparse vector of DFT coefficients. By combining (81) and (83), we obtain [31](84)
S
T
F
T
=
F
F
N
-
1
X
.Furthermore, if we assume thatSTFT vector is not completely available but has only a small percent of available samples, then we are dealing with the CS problem in the form(85)
y
S
T
F
T
=
A
X
,where y
S
T
F
T is a vector of available measurements fromSTFT, while A is the CS matrix obtained as partial combined transform matrix Ψ
=
F
F
N
-
1 with rows corresponding to the available samples y
S
T
F
T. In order to reconstruct the entireSTFT, we can define the minimization problem in the form(86)
m
i
n
X
1
s.t.
y
S
T
F
T
=
A
X
,which can be solved using any of the algorithms provided in the previous section.


## 4.1. CS in the Standard Fourier Transform Domain

As the simplest case of CS scenario, we will assume a multicomponent signal that consists ofK sinusoids that are sparse in the DFT domain:(43)
s
n
=
∑
i
=
1
K
A
i
e
j
2
π
k
i
n
/
N
,where the sparsity level is K
≪
N, where N is total signal length, while A
i and k
i denote amplitudes and frequencies of signal components, respectively. Sinces is sparse in the DFT domain, then we can write(44)
s
=
Ψ
F
=
F
-
1
F
,where F is a vector of DFT coefficients where at mostK coefficients are nonzero, while F
-
1 is the inverse Fourier transform matrix of size N
×
N.Ifs is a signal in compressive sensing application, then only the random measurements y
⊂
s are available and these are defined by the set ofM positions {
n
1
,
n
2
,
n
M
}. Therefore, the measurement process can be modeled by matrix Φ:(45)
y
=
Φ
F
-
1
F
=
A
F
.The CS matrix A represents a partial random inverse Fourier transform matrix obtained by omitting rows from F
-
1 that corresponds to unavailable samples positions:(46)
A
=
ψ
0
n
1
ψ
1
n
1
⋯
ψ
N
-
1
n
1
ψ
0
n
2
ψ
1
n
2
⋯
ψ
N
-
1
n
2
⋮
⋮
⋱
⋮
ψ
0
n
M
ψ
1
n
M
⋯
ψ
N
-
1
n
M
,with the Fourier basis functions: (47)
ψ
k
n
=
e
j
2
π
k
n
/
N
.The CS problem formulation is now given as (48)
min
⁡
F
1
s.t.
y
=
A
F
.


## 4.2. CS in the Polynomial Fourier Transform Domain

In this section, we consider the possibility of CS reconstruction of polynomial phase signals [28]. Observe the multicomponent polynomial phase signal vectors, with elementss(n): (49)
s
n
=
∑
i
=
1
K
s
i
n
=
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
,where the polynomial coefficients are assumed to be bounded integers. The assumption is that the signal can be considered as sparse in the polynomial Fourier transform (PFT) domain. Hence, the discrete PFT form is given by(50)
X
k
1
i
,
…
,
k
p
i
=
∑
n
=
0
N
-
1
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
×
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!
e
-
j
2
π
/
N
n
k
1
. If we choose a set of parameters (
k
2
i
,
k
3
i
,
…
,
k
p
i
) that match the polynomial phase coefficients (
a
2
i
,
a
3
i
,
…
,
a
p
i
)
(51)
k
2
i
,
k
3
i
,
…
,
k
p
i
=
a
2
i
,
a
3
i
,
…
,
a
p
i
,then the ith signal component is demodulated and becomes a sinusoid: r
i
e
j
(
2
π
/
N
)
n
a
1
i which is dominant in the PFT. In other words, the spectrum is highly concentrated at k
=
a
1
i. Therefore, we might say that if (51) is satisfied then the PFT is compressible with the dominant ith component. Note that the sparsity (compressibility) in the PFT domain is observed with respect to the single demodulated component.In order to define the CS problem in the PFT domain, instead of the signals itself, we will consider a modified signal form obtained after the multiplication with the exponential term:(52)
ξ
n
=
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!
,given in the PFT definition (50). The new signal form is obtained as(53)
z
n
=
s
n
ξ
n
=
∑
i
=
1
K
r
i
e
j
2
π
/
N
n
a
1
i
+
n
2
a
2
i
/
2
+
⋯
+
n
p
a
p
i
/
p
!
e
-
j
2
π
/
N
n
2
k
2
i
/
2
+
⋯
+
n
p
k
p
i
/
p
!or in the vector form:(54)
z
=
s
ξ
.Then, the PFT definition given by (50) can be rewritten in the vector form as follows:(55)
X
=
F
z
,where F is the discrete Fourier transform matrix (N
×
N).For a chosen set of parameters(
k
2
i
,
k
3
i
,
…
,
k
p
i
) in ξ that is equal to the set (
a
2
i
,
a
3
i
,
…
,
a
n
i
) in s, X
=
X
i is characterized by one dominant sinusoidal component at the frequency a
1
i.


### 4.2.1. CS Scenario

Now, assume thatz is incomplete in the compressive sensing sense, and instead of z we are dealing with M available measurements defined by vector y
z, and (56)
y
z
=
A
z
X
,where A
z is the partial random Fourier matrix obtained by omitting rows from F that correspond to the unavailable samples. When (
k
2
i
,
k
3
i
,
…
,
k
p
i
)
=
(
a
2
i
,
a
3
i
,
…
,
a
p
i
), thenX can be observed as a demodulated version of the ith signal component X
i, having the dominant ith component in the spectrum with the support k
1
i. The rest of the components in spectrum are much lower than X
i and could be observed as noise. Hence, we may write the minimization problem in the form(57)
min
⁡
X
i
1
s.t.
y
z
-
A
z
X
i
2
<
T
,whereT is certain threshold. The dominant components can be detected using the threshold based algorithm (Algorithm 5) described in the previous section. Using an iterative procedure, one may change the values of parameters k
2
i
,
…
,
k
p
i between k
m
i
n and k
m
a
x, ∀
i. The algorithm will detect the support k
1
i when the set (
k
2
i
,
k
3
i
,
…
,
k
p
i
) matches the set (
a
2
i
,
a
3
i
,
…
,
a
p
i
); otherwise, no component support is detected. Hence, as a result of this phase, we have identified the sets of signal phase parameters: k
i
=
k
1
i
,
k
2
i
,
…
,
k
n
i
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).Algorithm 5:Automated threshold based iterative algorithm.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸N
M
=
n
1
,
n
2
,
…
,
n
M, M –number of available samples, N-original signal length
▸ MatrixA: A
=
Φ
Ψ (A is obtained as a partial random Fourier matrix: A = Φ
Ψ
=
Ψ
{
ϑ
},
ϑ
=
{
ψ
n
1
,
ψ
n
2
,
…
,
ψ
n
M
}
▸ Measurement vectory
=
Φ
f
Output:
(1) i
=
1, N
M
=
{
n
1
,
n
2
,
…
,
n
M
}
(2) p
=
⌀.
(3) P
←
0.99                  (Set desired probability
P)
(4) σ
M
S
2
←
M
(
(
N
-
M
)
/
(
N
-
1
)
)
∑
i
=
1
M
(
y
(
i
)
2
/
M
)      (Calculate variance)
(5) T
i
←
-
σ
M
S
2
log
⁡
(
1
-
P
(
T
)
1
/
N
)           (Calculate threshold)
(6) X
i
←
A
-
1
y                  (Calculate initialDFTof
y)
(7) k
i
←
arg
⁡
X
i
>
T
i
/
N            (Find comp. in
X
i  above the T
i)
(8) Update
p
←
p
∪
k
i
(9) A
c
s
i
←
A
(
k
i
)                   (CS matrix)
(1
0) X
~
i
←
A
c
s
i
T
A
c
s
i
-
1
A
c
s
i
T
y
(11) for
∀
p
∈
p:
(1
2) y
←
y
-
(
M
/
N
)
X
i
p
exp
⁡
(
j
2
π
p
N
M
/
N
);            (Update
y)
(13)Update
σ
M
S
2
←
(
M
(
N
-
M
)
/
(
N
-
1
)
)
∑
y
2
/
M
(14)If
∑
y
2
/
M
<
δ
break; Else
(15)Set
i
=
i
+
1, and go to (5).
(16)return
X
~
i,  k


### 4.2.2. Reconstruction of Exact Components Amplitudes

The next phase is reconstruction of components amplitudes. Denote the set of measurements positions by(
n
1
,
n
2
,
…
,
n
M
), such that [
y
(
1
)
,
y
(
2
)
,
…
,
y
(
M
)
]
=
[
s
(
n
1
)
,
s
(
n
2
)
,
…
,
s
(
n
M
)
], while the detected sets of signal phase parameters are k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).In order to calculate the exact amplitudesr
1
,
r
2
,
…
,
r
K of signal components, we observe the set of equations in the form(58)
s
n
1
⋮
s
n
M
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
N
k
11
+
⋯
+
n
N
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
N
p
k
p
K
r
1
⋮
r
Kor in other words we have another system of equations given by (59)
y
=
A
R
,where R
=
r
1
,
…
,
r
K
T contains the desired K signal amplitudes. The matrix A of size N
×
K is based on the PFT: (60)
A
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
M
k
11
+
⋯
+
n
M
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
M
p
k
p
K
.The rows of A correspond to positions of measurements (
n
1
,
n
2
,
…
,
n
M
), and columns correspond to phase parameters k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
), for i
=
1
,
…
,
K. The solution of the observed problem can be obtained in the least squares sense as follows:(61)
X
=
A
∗
A
-
1
A
∗
y
.The resulting reconstructed signal is obtained as(62)
s
n
=
r
1
e
-
j
2
π
/
N
n
a
11
+
⋯
+
n
p
a
p
1
/
p
!
+
⋯
+
r
K
e
-
j
2
π
/
N
n
a
1
K
+
⋯
+
n
p
a
p
K
/
p
!
.


## 4.2.1. CS Scenario

Now, assume thatz is incomplete in the compressive sensing sense, and instead of z we are dealing with M available measurements defined by vector y
z, and (56)
y
z
=
A
z
X
,where A
z is the partial random Fourier matrix obtained by omitting rows from F that correspond to the unavailable samples. When (
k
2
i
,
k
3
i
,
…
,
k
p
i
)
=
(
a
2
i
,
a
3
i
,
…
,
a
p
i
), thenX can be observed as a demodulated version of the ith signal component X
i, having the dominant ith component in the spectrum with the support k
1
i. The rest of the components in spectrum are much lower than X
i and could be observed as noise. Hence, we may write the minimization problem in the form(57)
min
⁡
X
i
1
s.t.
y
z
-
A
z
X
i
2
<
T
,whereT is certain threshold. The dominant components can be detected using the threshold based algorithm (Algorithm 5) described in the previous section. Using an iterative procedure, one may change the values of parameters k
2
i
,
…
,
k
p
i between k
m
i
n and k
m
a
x, ∀
i. The algorithm will detect the support k
1
i when the set (
k
2
i
,
k
3
i
,
…
,
k
p
i
) matches the set (
a
2
i
,
a
3
i
,
…
,
a
p
i
); otherwise, no component support is detected. Hence, as a result of this phase, we have identified the sets of signal phase parameters: k
i
=
k
1
i
,
k
2
i
,
…
,
k
n
i
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).Algorithm 5:Automated threshold based iterative algorithm.
Input:
▸ Transform matrixΨ, Measurement matrix Φ
▸N
M
=
n
1
,
n
2
,
…
,
n
M, M –number of available samples, N-original signal length
▸ MatrixA: A
=
Φ
Ψ (A is obtained as a partial random Fourier matrix: A = Φ
Ψ
=
Ψ
{
ϑ
},
ϑ
=
{
ψ
n
1
,
ψ
n
2
,
…
,
ψ
n
M
}
▸ Measurement vectory
=
Φ
f
Output:
(1) i
=
1, N
M
=
{
n
1
,
n
2
,
…
,
n
M
}
(2) p
=
⌀.
(3) P
←
0.99                  (Set desired probability
P)
(4) σ
M
S
2
←
M
(
(
N
-
M
)
/
(
N
-
1
)
)
∑
i
=
1
M
(
y
(
i
)
2
/
M
)      (Calculate variance)
(5) T
i
←
-
σ
M
S
2
log
⁡
(
1
-
P
(
T
)
1
/
N
)           (Calculate threshold)
(6) X
i
←
A
-
1
y                  (Calculate initialDFTof
y)
(7) k
i
←
arg
⁡
X
i
>
T
i
/
N            (Find comp. in
X
i  above the T
i)
(8) Update
p
←
p
∪
k
i
(9) A
c
s
i
←
A
(
k
i
)                   (CS matrix)
(1
0) X
~
i
←
A
c
s
i
T
A
c
s
i
-
1
A
c
s
i
T
y
(11) for
∀
p
∈
p:
(1
2) y
←
y
-
(
M
/
N
)
X
i
p
exp
⁡
(
j
2
π
p
N
M
/
N
);            (Update
y)
(13)Update
σ
M
S
2
←
(
M
(
N
-
M
)
/
(
N
-
1
)
)
∑
y
2
/
M
(14)If
∑
y
2
/
M
<
δ
break; Else
(15)Set
i
=
i
+
1, and go to (5).
(16)return
X
~
i,  k


## 4.2.2. Reconstruction of Exact Components Amplitudes

The next phase is reconstruction of components amplitudes. Denote the set of measurements positions by(
n
1
,
n
2
,
…
,
n
M
), such that [
y
(
1
)
,
y
(
2
)
,
…
,
y
(
M
)
]
=
[
s
(
n
1
)
,
s
(
n
2
)
,
…
,
s
(
n
M
)
], while the detected sets of signal phase parameters are k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
).In order to calculate the exact amplitudesr
1
,
r
2
,
…
,
r
K of signal components, we observe the set of equations in the form(58)
s
n
1
⋮
s
n
M
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
N
k
11
+
⋯
+
n
N
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
N
p
k
p
K
r
1
⋮
r
Kor in other words we have another system of equations given by (59)
y
=
A
R
,where R
=
r
1
,
…
,
r
K
T contains the desired K signal amplitudes. The matrix A of size N
×
K is based on the PFT: (60)
A
=
e
j
2
π
/
N
n
1
k
11
+
⋯
+
n
1
p
k
p
1
⋯
e
j
2
π
/
N
n
1
k
1
K
+
⋯
+
n
1
p
k
p
K
e
j
2
π
/
N
n
2
k
11
+
⋯
+
n
2
p
k
p
1
⋯
e
j
2
π
/
N
n
2
k
1
K
+
⋯
+
n
2
p
k
p
K
⋮
⋮
e
j
2
π
/
N
n
M
k
11
+
⋯
+
n
M
p
k
p
1
⋯
e
j
2
π
/
N
n
M
k
1
K
+
⋯
+
n
M
p
k
p
K
.The rows of A correspond to positions of measurements (
n
1
,
n
2
,
…
,
n
M
), and columns correspond to phase parameters k
i
=
(
k
1
i
,
k
2
i
,
…
,
k
n
i
)
=
(
a
1
i
,
a
2
i
,
…
,
a
n
i
), for i
=
1
,
…
,
K. The solution of the observed problem can be obtained in the least squares sense as follows:(61)
X
=
A
∗
A
-
1
A
∗
y
.The resulting reconstructed signal is obtained as(62)
s
n
=
r
1
e
-
j
2
π
/
N
n
a
11
+
⋯
+
n
p
a
p
1
/
p
!
+
⋯
+
r
K
e
-
j
2
π
/
N
n
a
1
K
+
⋯
+
n
p
a
p
K
/
p
!
.


## 4.3. CS in the Hermite Transform Domain

The Hermite expansion usingN Hermite functions can be written in terms of the Hermite transform matrix. First, let us define the Hermite transform matrix H (of size N
×
N) [29, 30]:(63)
H
=
1
N
ψ
0
1
ψ
N
-
1
1
2
ψ
0
2
ψ
N
-
1
2
2
⋯
ψ
0
N
ψ
N
-
1
N
2
ψ
1
1
ψ
N
-
1
1
2
ψ
1
2
ψ
N
-
1
2
2
⋯
ψ
1
N
ψ
N
-
1
N
2
⋮
⋮
⋱
⋮
ψ
N
-
1
1
ψ
N
-
1
1
2
ψ
N
-
1
2
ψ
N
-
1
2
2
⋯
ψ
N
-
1
N
ψ
N
-
1
N
2
.Here, the pth order Hermite basis function is defined using the pth order Hermite polynomial as follows: (64)
ψ
p
n
=
σ
2
p
p
!
π
-
1
/
2
e
-
n
2
/
2
H
p
n
σ
,where σ is the scaling factor used to “stretch” or “compress” Hermite functions, in order to match the signal.The Hermite functions are usually calculated using the following fast recursive realization:(65)
ψ
0
n
=
1
π
4
e
-
n
2
/
2
σ
2
,
ψ
1
n
=
2
n
π
4
e
-
n
2
/
2
σ
2
,
ψ
i
n
=
n
2
i
Ψ
i
-
1
n
-
i
-
1
i
Ψ
i
-
2
n
,
∀
i
≥
2
.The Hermite transform coefficients can be calculated in the matrix form as follows: (66)
c
=
H
s
,where the vector of Hermite transform coefficients is c
=
c
0
,
c
1
,
…
,
c
N
-
1
T, while the vector of signal samples is s
=
[
s
(
1
)
,
s
(
2
)
,
…
,
s
(
N
)
]
T.Following the Gauss-Hermite approximation, the inverse matrixΨ
=
H
-
1 contains N Hermite functions and it is given by(67)
Ψ
=
ψ
1
1
ψ
2
1
⋯
ψ
N
1
ψ
1
2
ψ
2
2
⋯
ψ
N
2
⋮
⋮
⋱
⋮
ψ
1
N
ψ
2
N
⋯
ψ
N
N
.Now, in the context of CS, let us assume that a signal s isK-sparse in the Hermite transform domain, meaning that it can be represented by a small set of nonzero Hermite coefficientsc: (68)
s
n
=
∑
p
=
k
1
k
K
c
p
ψ
p
n
,where c
p is the pth element from the vector c of Hermite expansion coefficients. In the matrix form we may write (69)
s
=
Ψ
c
,such that most of the coefficients in c are zero values. Furthermore, assume that only the compressive sensing is done using random selection ofM signal values [
n
1
,
n
2
,
…
,
n
M
]. Then, the CS matrix A can be defined as random partial inverse Hermite transform matrix: (70)
A
=
ψ
0
n
1
ψ
1
n
1
⋯
ψ
N
-
1
n
1
ψ
0
n
2
ψ
1
n
2
⋯
ψ
N
-
1
n
2
⋮
⋮
⋱
⋮
ψ
0
n
M
ψ
1
n
M
⋯
ψ
N
-
1
n
M
.For a vector of available measurements vector y
=
(
s
(
n
1
)
,
s
(
n
2
)
,
…
,
s
(
n
M
)
), the linear system of equations (undetermined system ofM linear equations and N unknowns) can be written as(71)
y
=
A
c
.Since c has only K nonzero components, the reconstruction problem can be defined as(72)
m
i
n
c
1
s.t.
y
=
A
c
.If we identify the signal support in the Hermite transform domain by the set of indices [
k
1
,
k
2
,
…
,
k
K
], then the problem can be solved in the least squares sense as(73)
c
^
=
Θ
c
s
∗
Θ
c
s
-
1
Θ
c
s
∗
y
,where (
∗
) denotes the conjugate transpose operation, while(74)
Θ
c
s
=
ψ
k
1
n
1
ψ
k
1
n
2
…
ψ
k
1
n
M
ψ
k
2
n
1
ψ
k
2
n
2
⋯
ψ
k
2
n
M
⋮
⋮
⋱
⋮
ψ
k
K
n
1
ψ
k
K
n
2
⋯
ψ
k
K
n
M
,and c
^
=
c
k
1
,
c
k
2
,
…
,
c
k
K
T. The estimated signal s can be now obtained according to (68).


## 4.4. CS in the Time-Frequency Domain

Starting from the definition of the short-time Fourier transform(75)
STFT
n
,
k
=
∑
m
=
0
M
-
1
x
n
+
m
e
-
j
2
π
m
k
/
M
,we can define the following STFT vector and signal vector: (76)
S
T
F
T
M
n
=
STFT
n
,
0
,
…
,
STFT
n
,
M
-
1
T
,
x
n
=
x
n
,
x
n
+
1
,
…
,
x
n
+
M
-
1
T
,respectively. The STFT at an instantn can be now rewritten in the matrix form:(77)
S
T
F
T
M
n
=
F
M
x
n
,where F
M is the standard Fourier transform matrix of size M
×
M with elements(78)
F
M
m
,
k
=
e
-
j
2
π
k
m
/
M
.If we consider the nonoverlapping signal segments for the STFT calculation, then n takes the values from the set {
0
,
M
,
2
M
,
…
,
N
-
M
}, where M is the window length while N is the total signal length. The STFT vector, for n
∈
{
0
,
M
,
2
M
,
…
,
N
-
M
}, is given by(79)
S
T
F
T
=
S
T
F
T
M
0
S
T
F
T
M
M
⋮
S
T
F
T
M
N
-
M
.Consequently, (77) can be extended to the entire STFT for all considered values of n
∈
{
0
,
M
,
2
M
,
…
,
N
-
M
} as follows [31, 32]:(80)
S
T
F
T
M
0
S
T
F
T
M
M
⋮
S
T
F
T
M
N
-
M
=
F
M
0
M
⋯
0
M
0
M
F
M
⋯
0
M
⋯
⋯
⋯
⋯
0
M
0
M
⋯
F
M
x
0
x
M
⋮
x
N
-
M
,or in the compact form(81)
S
T
F
T
=
F
x
.In the case of time-varying windows with lengths {
M
1
,
M
2
,
…
,
M
Q
}, the previous system of equations becomes(82)
S
T
F
T
M
0
0
S
T
F
T
M
1
M
0
⋮
S
T
F
T
M
Q
N
-
M
Q
=
F
M
0
0
⋯
0
0
M
F
M
1
⋯
0
⋯
⋯
⋯
⋯
0
0
⋯
F
M
Q
x
0
x
M
⋮
x
N
-
M
.Note that the signal vector is x
=
[
x
(
0
)
T
,
x
(
M
)
T
,
…
,
x
(
N
-
M
)
T
]
T
=
[
x
(
0
)
,
x
(
1
)
,
…
,
x
(
N
-
1
)
]
T and can be expressed using the N-point inverse Fourier transform as follows:(83)
x
=
F
N
-
1
X
,whereX is a sparse vector of DFT coefficients. By combining (81) and (83), we obtain [31](84)
S
T
F
T
=
F
F
N
-
1
X
.Furthermore, if we assume thatSTFT vector is not completely available but has only a small percent of available samples, then we are dealing with the CS problem in the form(85)
y
S
T
F
T
=
A
X
,where y
S
T
F
T is a vector of available measurements fromSTFT, while A is the CS matrix obtained as partial combined transform matrix Ψ
=
F
F
N
-
1 with rows corresponding to the available samples y
S
T
F
T. In order to reconstruct the entireSTFT, we can define the minimization problem in the form(86)
m
i
n
X
1
s.t.
y
S
T
F
T
=
A
X
,which can be solved using any of the algorithms provided in the previous section.


## 5. Experimental Evaluation

Example 1.
Let us observe the QRS complex extracted from the real ECG signal. It has been known that these types of signals are sparse in the Hermite transform (HT) domain [33, 34]. As an illustration, we present the Hermite transform and discrete Fourier transform (DFT) in Figure 1, from which it can be observed that HT provides sparse representation while DFT is dense. The analyzed signal is obtained from the MIT-ECG Compression Test Database [35], originally sampled uniformly with T
=
1
/
250
[
s
] being the sampling period, with amplitude gain 1/400. Extracted QRS is of length N
=
51, centered at the R peak. Note that, in order to provide the sparsest HT representation, the scaling factor σ should be set to the optimal value [33, 34], which is 6.8038 in this example, that is, σ
Δ
t
=
0.027 in seconds, proportional to the scaling factor 0.017 presented in [34] for signal with 27 samples, when it is scaled with the signal lengths ratio. Also, the QRS signal is resampled at the zeros of Hermite polynomial using sinc interpolation functions as described in [33].
The QRS signal is subject to compressed sensing approach, meaning that it is represented by 50% of available measurements. The available measurements and the corresponding HT (which is not sparse in the CS conditions) are shown in Figure2.
In order to reconstruct the entire QRS signal from available measurements, two of the presented algorithms are employed, namely, the OMP and gradient-based reconstruction algorithm.
The reconstruction results are shown in Figure3 in both time and transform domain. It can be observed that both algorithms provide quite successful reconstruction performance, which is certainly better in the case of gradient-based approach. In order to measure and compare the reconstruction performance of these two algorithms, the mean square error (MSE) is calculated, showing that the gradient-based algorithm in the considered case provides MSE = 13.85, much lower than MSE = 31.64 for the OMP algorithm.Figure 1
(a) Hermite transform coefficients of the selected QRS complex (σ
=
6.8038) and (b) Fourier transform coefficients (absolute values).
(a)
(b)Figure 2
(a) Measurements of the selected QRS complex (50% of the total samples number) and (b) Hermite transform of measurements vector.
(a)
(b)Figure 3
Results of signal reconstruction using gradient-based algorithms: (a) original signal (solid line) and reconstructed signal (dash line) and (b) Hermite coefficients of the original signal (solid line) and the reconstructed signal (squares). Reconstruction results using OMP algorithm: (c) original signal (solid) and reconstructed signal (dash) in time domain and (d) Hermite coefficients of original signal (solid line) and reconstructed signal (squares).
(a)
(b)
(c)
(d)Example 2.
In this example, the reconstruction of images using adaptive gradient-based algorithm (Algorithm6) is considered. The images are generally not strictly sparse in any domain and as such could be quite demanding for reconstruction. For that reason, as stated in Section 3.7, the adaptive gradient-based algorithm is employed, since it does not require the strict sparsity condition. As a suitable transform domain representation for natural image, the two-dimensional DCT (
I
→
2
D
D
C
T
) is used. We might say that the image can be observed as approximately sparse in the 2DDCT domain.
The image with missing pixels is divided into16
×
16 blocks, where 40% of the pixels are missing in each block (missing pixels are denoted in white in Figures 4(a), 4(c), and 4(e)). Hence, the algorithm is applied on a block-by-block basis where the measurement vector y is created for each observed image block using the available block pixels. The value of gradient step Δ
=
128 is used. The resulting reconstructed images are shown in Figures 4(b), 4(d), and 4(f).
The quality of the reconstructed image is measured using the structural similarity index (SSIM) [36] between the original image (with full set of pixels) and reconstructed image. The SSIM ranges between 0 and 1, 1 being identical and 0 being not similar. It is designed to improve the traditional metrics such as peak signal-to-noise ratio (PSNR) in order to be more consistent with the perception of the human visual system. The achieved values of SSIM for the test images in Figure 4 are SSIM = 0.9, SSIM = 0.94, and SSIM = 0.95.Algorithm 6:Adaptive gradient-based algorithm.
Input:
▸ Available samples positionN
M
=
n
1
,
n
2
,
…
,
n
M
▸ Missing samples positionN
∖
N
M, where N
=
{
1,2
,
…
,
N
}
▸ Available samplesf
n, n
∈
N
M
▸ Measurement vector:y
n
=
f
(
n
), for n
∈
N
M; 0, for n
∈
N
∖
N
M
Output:
(1) Set i
=
0, y
(
0
)
(
n
)
←
y
(
n
), Δ
←
max
⁡
y
(
0
)
(
n
)
(2) repeat
(3) Set  y
p
(
n
)
=
y
(
i
)
(
n
)
(4) repeat
(5)  i
←
i
+
1
(6)  for  n
s
←
0  to  N
-
1  do
(7)   if
n
s
∈
N
∖
N
M
(8)    x
+
(
k
)
←
I
y
(
i
)
(
n
)
+
Δ
δ
(
n
-
n
s
)
(9)    x
-
(
k
)
←
I
y
(
i
)
(
n
)
-
Δ
δ
(
n
-
n
s
)
(10)     G
(
i
)
(
n
s
)
←
(
1
/
N
)
x
+
1
-
x
-
1
(11)   else
G
(
i
)
(
n
s
)
←
0
(12)   end if
(13)  y
(
i
+
1
)
(
n
s
)
←
y
(
i
)
(
n
s
)
-
G
(
i
)
(
n
s
)
(14)  end for
(15)  β
i
=
arcos
⁡
(
G
i
-
1
G
i
/
G
i
-
1
2
2
G
i
2
2
)
(16) until
β
i
<
170
°
(17) Δ
←
Δ
/
10;
(18) Ξ
=
10
log
10
⁡
(
∑
n
∈
N
∖
N
M
y
p
(
n
)
-
y
(
i
)
(
n
)
2
/
∑
n
∈
N
∖
N
M
y
(
i
)
(
n
)
2
)
(19) until
Ξ
<
Ξ
m
a
x
(20)return
y
(
i
)
(
n
)
Reconstructed signalx
R
(
n
)
=
y
(
i
)
(
n
)Figure 4
((a), (c), and (e)) Available measurements/pixels (missing values are denoted by white pixels). ((b), (d), and (f)) Reconstructed images.
(a)
(b)
(c)
(d)
(e)
(f)Example 3.
Let us assume that a radar signal consists of stationary (sinusoidal) components belonging to the rigid body and nonstationary components belonging to the disturbances caused by rotating target parts and noise. The Fourier transform and the STFT of the considered radar signal type are given in Figure5.
The total signal length is 4096 samples. The STFT is calculated for the window widthw
=
64 samples (STFT is of size 64
×
64). The aim of this experiment is to extract the rigid body components which are presented as a set of stationary sinusoids.
The STFT values of the original full signal along columns and then the 75% of the largest values (possibly belonging to rotating components and noise) are set to zero or in other words discarded from the STFT plane. Also, a small percentage of the smallest values along the sorted columns are also set to zero with the aim of discarding small values belonging to noise. The remaining STFT values are unsorted and placed to their original positions in the STFT plane. These values represent a small amount of available measurements in the STFT domain and we are faced with the CS problem. In Figure6, we show the available measurements from the STFT after discarding all unwanted components.
In order to reconstruct the stationary rigid body components from the available measurements, the CS matrix is defined as partial combined transform matrix with rows corresponding to the available samples in the STFT (according to the procedure described in Section4.4). The STFT of the reconstructed stationary components is given in Figure 6(a), while the corresponding DFT is given in Figure 6(b).Figure 5
(a) STFT of the observed signal; (b) DFT of the observed signal.
(a)
(b)Figure 6
(a) STFT of the observed signal; (b) DFT of the observed signal.
(a)
(b)


## 6. Conclusion

The paper reviews the fundamental concepts of compressive sensing theory, comprising the main requirements, conditions, and common optimization problem formulations. Several algorithms for signal reconstruction from incomplete measurements are summarized, allowing an insight into the diversity of approaches, the related complexity, assumptions, and efficiency. As a special contribution, the paper presents some interesting compressive sensing formulations for different types of signals and sparsity domains, such as the Fourier transform domain, polynomial Fourier transform domain, Hermite transform, and time-frequency domain. Particularly, this part aims to reveal the possibilities and perspectives of using compressive sensing in different signal processing scenarios.



---

*Source: 7616393-2016-10-25.xml*