Quantization in Compressive Sensing: A Signal Processing Approach

TL;DR

This paper investigates how finite-length registers and quantization affect the accuracy of compressive sensing signal reconstruction, providing a mathematical model and validating it with numerical examples.

Contribution

It introduces a unified model for quantization noise in CS and derives an exact formula for expected reconstruction error, considering signal nonsparsity and quantization effects.

Findings

Quantization impacts CS reconstruction accuracy significantly.

Derived an exact formula for expected error in quantized CS.

Validated the theory with numerical examples including noise and approximate sparsity.

Abstract

Influence of the finite-length registers and quantization effects on the reconstruction of sparse and approximately sparse signals is analyzed in this paper. For the nonquantized measurements, the compressive sensing (CS) framework provides highly accurate reconstruction algorithms that produce negligible errors when the reconstruction conditions are met. However, hardware implementations of signal processing algorithms involve the finite-length registers and quantization of the measurements. An analysis of the effects related to the measurements quantization with an arbitrary number of bits is the topic of this paper. A unified mathematical model for the analysis of the quantization noise and the signal nonsparsity on the CS reconstruction is presented. An exact formula for the expected energy of error in the CS-based reconstructed signal is derived. The theory is validated through various numerical examples with quantized measurements, including the cases of approximately sparse signals, noise folding, and floating-point arithmetics.

Tables

Table 1. TABLE I: Variances in the initial estimate for various measurement matrix types for k ≠ k i 𝑘 subscript 𝑘 𝑖 k\neq k_{i}

Measurement matrix	$a_m\left(k\right)$	${\sigma }_{\mu }^2$
Partial DFT	$\frac{1}{\sqrt{M}}{e}^{j2\pi n_m\frac{k}{N}}$	$\frac{N-M}{M\left(N-1\right)}$
Random partial DFT	$\frac{1}{\sqrt{M}}{e}^{j2\pi t_m\frac{k}{N}}$	$1/M$
Equiangular tight frame	$\mu =\sqrt{\frac{N-M}{M\left(N-1\right)}}$	$\frac{N-M}{M\left(N-1\right)}$
Gaussian random	$\sim 𝒩(0,\frac{1}{M})$	$1/M$
Uniform random	$\sim 𝒰(0,\frac{1}{M})$	$1/M$
Bernoulli (binary)	$\sim \pm \frac{1}{\sqrt{M}}$	$1/M$

Equations

$$X(k)=0\text{ for }k\notin \mathbb{K}=\{k_1,k_2,\dots ,k_K\},$$

$$y(m)=\sum _{k=0}^{N-1}{a}_m(k)X(k),m=1,\dots ,M,$$

$$\mathbf{y}=\mathbf{AX},$$

$$y(m)=\sum _{i=1}^{K}X(k_i)a_m(k_i).$$

$$\mathbf{y}=[y(1),y(2),\dots ,y(M){]}^T.$$

$$\text{min }{\Vert \mathbf{X}\Vert }_0\text{ subject to }\mathbf{y}=\mathbf{AX},$$

$$\mathbf{y}={\mathbf{A}}_{MK}{\mathbf{X}}_K.$$

$${\mathbf{X}}_K=({\mathbf{A}}_{MK}^H{\mathbf{A}}_{MK}{)}^{-1}{\mathbf{A}}_{MK}^H\mathbf{y}=\text{pinv}({\mathbf{A}}_{MK})\mathbf{y},$$

$${\mathbf{X}}_0={\mathbf{A}}^H\mathbf{y}={\mathbf{A}}^H\mathbf{A}\mathbf{X},$$

$$X_0(k)=\sum _{m=1}^{M}y(m)a_m^{\ast }(k),$$

$$X_0(k)=\sum _{i=1}^{K}X(k_i)\mu (k_i,k),$$

$$\mu (k_i,k)=\sum _{m=1}^M{a}_m(k_i)a_m^{\ast }(k)$$

$$\mu =\underset{k\ne l}{\max }|\mu (l,k)|,$$

$$E\{X_0(k)\}=\sum _{i=1}^{K}X(k_i)\delta (k-k_i)$$

$$\mathrm{var}\{X_0(k)\}=\sum _{i=1}^K|X(k_i){|}^2\mathrm{var}\{\mu (k_i,k)\}\left(1-\delta (k-k_i)\right),$$

$$\mathrm{var}\{\mu (k_i,k)\}={\sigma }_{\mu }^2.$$

$$1-(K-1)\mu >K\mu ,$$

$${\mathbf{X}}_{0K}=\mathbf{B}{\mathbf{X}}_K={\mathbf{X}}_K+{\mathbf{C}}_K,$$

$${\mathbf{X}}_R=({\mathbf{A}}_K^H{\mathbf{A}}_K{)}^{-1}{\mathbf{A}}_K^H\mathbf{y},$$

$${\mathbf{X}}_R=({\mathbf{A}}_K^H{\mathbf{A}}_K{)}^{-1}{\mathbf{X}}_{0K}.$$

$$({\mathbf{A}}_K^H{\mathbf{A}}_K{)}^{-1}={\mathbf{B}}^{-1}$$

$$\mathbf{y}+\mathbf{\epsilon }=\mathbf{AX},$$

$$\mathrm{E}\{\mathbf{\epsilon }{\mathbf{\epsilon }}^H\}={\sigma }_{\epsilon }^2\mathbf{I}.$$

$$\mathrm{E}\{{\mathbf{X}}_0{\mathbf{X}}_0^H\}-|\mathrm{E}\{{\mathbf{X}}_0\}{|}^2=\mathrm{E}\{{\mathbf{A}}^H\mathbf{\epsilon }{\mathbf{\epsilon }}^H\mathbf{A}\}={\sigma }_{\epsilon }^2\mathbf{I}.$$

$$\mathrm{var}\{X_R(k)\}={\sigma }_{\epsilon }^2.$$

$${\Vert {\mathbf{X}}_R-{\mathbf{X}}_K\Vert }_2^2=K{\sigma }_{\epsilon }^2.$$

$${\mathbf{y}}_B={\text{digital}}_B\{\mathbf{AX}\},$$

$${\mathbf{y}}_B={\text{digital}}_B\{\Re \{\mathbf{AX}\}\}+j{\text{digital}}_B\{\Im \{\mathbf{AX}\}\}$$

$$-1\le y(m)<1.$$

$$|e(m)|<\Delta /2,$$

$$\Delta =2^{-B}.$$

Full text

Keywords — compressive sensing, quantization, signal reconstruction, sparse signal processing

I Introduction

Compressive sensing (CS) theory provides a rigorous mathematical framework for the reconstruction of sparse signals, using a reduced set of measurements [1, 2, 7, 3, 4, 5, 6, 8, 9]. Advantages of CS are directly related to the signal transmission and storage efficiency, which is crucial in big data setups. Moreover, the problem of the physical unavailability of measurements, or the problem of a significant signal corruption, are also potentially solvable within the CS framework. Since the establishment of CS, the phenomena related to the reduced sets of measurements and sparse signal reconstruction have been supported by the fundamental theory and well-defined mathematical framework, while the performance of the reconstruction processes have been continuously improved by newly introduced algorithms, often adapted to perform in a particular context, or to solve specific problems [10, 11, 12, 15, 13, 14]. In real applications, many signals are sparse or approximately sparse in a certain transformation domain. This makes the CS applicable in various fields of signal processing [14].

Ideally, the measurements used for the reconstruction should be taken accurately, assuming a very large number of bits in their digital form (providing high precision levels). However, this could be extremely demanding and expensive for hardware implementations [16]. Therefore, in practice, the measurements are quantized, meaning that they are represented using a limited number of bits. Such measurements bring robustness, memory efficiency and simplicity in the corresponding hardware implementation (particularly in sensor design). This paper investigates the quantization influence on the CS reconstruction with a simple yet rigorous characterization of the related phenomena, through the derivation of the corresponding errors. The theory is supported by a relevant theoretical framework and detailed statistical analysis, through extensive numerical experiments.

The most extreme case of quantization is limiting the measurements using only one bit. In previous work [16, 17, 18, 19], one-bit measurements are initially treated as sign constraints, as opposed to the values to be matched in the mean squared sense during the reconstruction process. Quantization to one-bit measurements is suitable for hardware systems since the quantizers do not suffer from dynamic range issues. However, as signs of measurements do not provide amplitude information of the signal, it can be recovered up to a constant scalar factor only. Additionally, more measurements are needed for a successful reconstruction for such systems, exceeding the signal length. In this paper, we focus on the general $B$-bit quantization of available measurements and its effect to the reconstruction accuracy.

If the reconstruction is performed without measurements quantization, the error will be equal to zero or negligible. The quantization step in measurements inevitably introduces reconstruction errors. The effects of quantization in compressive sensing theory has been recently presented [20, 21, 22, 23, 24, 25]. The results mainly include the derivation of quantization error bounds and adaptation of CS algorithms aiming to reduce the distortions related to the quantization [4, 5]. The upper bound of the reconstruction error, for strictly sparse signals, has been derived in [20]. Other reported results are focused on the worst case analysis [21]. Exact asymptotic distortion rate functions have been derived in [21] for scalar quantization, where the reconstruction strategies have been adapted to accommodate quantization errors. An overview of the quantization phenomena in the compressive sensing context is presented in [22]. Therein, the fundamental analysis provides the performance bounds only, with an additional focus on the Sigma-Delta quantization and the related theory. Recently, the effect of quantization on the estimation of sparsity order, support and signals have been studied within a large number of Monte Carlo simulations in [23]. The most frequently used algorithms in compressive sensing are adjusted to the quantization effect in [24]. For the case of one-bit unlimited sampling, a quantization approach using the one-bit modulo samples, [25] shows the bounds of the reconstruction error.

This paper aims to fill the literature gap regarding the exact characterization of the quantization in compressed sensing, by deriving an explicit relation for the mean squared error, instead of the error bounds. The error produced by the quantization of measurements is analyzed from a practical signal processing point of view. Additional to that, the error appearing when approximately sparse signals are reconstructed under the sparsity constraint is also examined. The analysis is expanded to include the effect of the pre-measurements noise in the sparsity domain coefficients, known as the noise folding [26]. The presented theory is unified by exact relations for the expected squared reconstruction errors, derived to take into account all the studied effects. The results are validated using three different reconstruction algorithms. Moreover, we comment on the modifications of the derived relations, required to include the floating-point arithmetics.

The paper is organized as follows. In Section II, basic CS concepts and definitions are briefly presented. Section III introduces a common approach to solve the CS reconstruction problem, including also a brief overview of relevant properties which characterize possible solutions. Section IV puts the quantization within the compressive sensing framework. In Section V, the concept of nonsparse (approximately sparse) signals reconstructed under the sparsity constraint is analyzed, leading to the reconstruction error equation which unifies the studied effects. The theory is expanded, to take into account the noise folding effect, in Section VI, while the Section VII discusses the quantization in floating-point arithmetics. Numerical results verify the presented theory in Section VIII. The paper ends with the concluding remarks.

II Basic Compressive Sensing Definitions

Definition: A discrete signal $x(n)$, $n=0,1,\dots,N-1$ is sparse in one of its representation domains $X(k)$ if the number $K$ of nonzero coefficients is much smaller than the total number of samples $N$, that is,

[TABLE]

where $K\ll N$.

Definition: A measurement of a signal is a linear combination of its sparsity domain coefficients $X(k)$,

[TABLE]

or in matrix form

[TABLE]

where $\mathbf{y}$ is an $M\times 1$ ($M$-dimensional) column vector of the measurements $y(m)$, $\mathbf{A}$ is an $(M\times N)$-dimensional measurement matrix with coefficients $a_{m}(k)$ as its elements, and $\mathbf{X}$ is an $N\times 1$ ($N$-dimensional) column sparse vector of coefficients $X(k)$. It is common to normalize the measurement matrix such that the energy of its columns is $1$. In that case, the diagonal elements of the matrix $\mathbf{A}^{H}\mathbf{A}$ are equal to $1$, where $\mathbf{A}^{H}$ denotes a Hermitian transpose of $\mathbf{A}$.

By definition, a measurement of a $K$-sparse signal can be written as

[TABLE]

The compressive sensing theory states that, under certain realistic conditions, it is possible to reconstruct a sparse $N$-dimensional vector $\mathbf{X}$ from a reduced $M$-dimensional set of measurements ($M<N$) belonging to vector $\mathbf{y}$,

[TABLE]

The reconstruction conditions are defined in several forms. The most widely used are the forms based on the restricted isometry property (RIP) and the coherence index [1, 2, 3, 4]. Although providing tight bounds, the RIP based condition is of high calculation complexity. This is the reason why the coherence based relation will be considered in this paper, along with some comments on its probabilistic relaxation.

The reconstruction of a $K$-sparse signal, $\mathbf{X}$ is unique if $K<\left(1+1/\mu\right)/2,$ where the coherence index, $\mu$, is equal to the maximum absolute off-diagonal element of $\mathbf{A}^{H}\mathbf{A}$, assuming its unity diagonal elements.

A simple proof will be provided later.

Formally, compressive sensing aims to solve the optimization problem

[TABLE]

or its corresponding relaxed convex form. Amongst many others, an approach based on matching the components corresponding to the nonzero coefficients, can be used to solve (4). It is further assumed that the CS reconstruction is based on a such methodology. The solution is discussed in the next section, since it will be used to model the quantization noise and other studied effects.

III Problem Solution

To perform the reconstruction, we use an iterative version of the orthogonal matching pursuit algorithm from [10]. Assume that $K$ nonzero values $X(k)$ are detected at positions $k\in\mathbb{K}=\{k_{1},k_{2},\dots,k_{K}\}$. The system of measurement equations becomes

[TABLE]

The system is solved for the nonzero coefficients $X(k)$, $k\in\mathbb{K}$ written in vector form as $\mathbf{X}_{K}$. The matrix $\mathbf{A}_{MK}$ is an $M\times K$ sub-matrix of the $M\times N$ measurement matrix $\mathbf{A}$, where only the columns corresponding to the positions of nonzero elements in $X(k)$ are kept. The solution of the system is

[TABLE]

where $\textrm{pinv}(\mathbf{A}_{MK})$ is the pseudo-inverse of the matrix $\mathbf{A}_{MK}$ and $\mathbf{A}^{H}_{MK}\mathbf{A}_{MK}$ is known as a $K\times K$ Gram matrix of $\mathbf{A}_{MK}$.

Therefore, the problem solution can be split into two steps:

detect the positions of nonzero coefficients, and 2. 2.

apply the reconstruction algorithm at detected positions.

III-A Initial Estimate

Detection of the positions of nonzero coefficients $X(k)$ will be based on the initial estimate concept. An intuitive idea for the initial estimate comes from the fact that the measurements are obtained as linear combinations of the sparsity domain coefficients, with rows of the measurement matrix $\mathbf{A}$ acting as weights. It means that the back-projection of the measurements $\mathbf{y}$ to the measurement matrix $\mathbf{A}$, defined by

[TABLE]

can be used to estimate the positions of nonzero coefficients. For the coefficient at the $k$th position, its initial estimate $X_{0}(k)$ takes the following form:

[TABLE]

or

[TABLE]

where

[TABLE]

are the coefficients of mutual influence (interference) among elements $X(k)$. They are equal to the elements of matrix $\mathbf{A}^{H}\mathbf{A}$, with

[TABLE]

and $\mu(k,k)=1.$ Note that the $\mu$ represents also the coherence index condition.

For various values of $k_{i}$, the off-diagonal elements $\mu(k_{i},k)$ of matrix $\mathbf{A}^{H}\mathbf{A}$ act as random variables, with different distribution for different measurement matrices. For the partial discrete Fourier transform (DFT) matrix, distribution of $\mu(k_{i},k)$ tends to a Gaussian distribution for $1\ll M\ll N$, while for an equiangular tight frame (ETF) measurement matrix, $\mu(k_{i},k)$ takes only the values such that $|\mu(k_{i},k)|=\mu$. Distributions of $\mu(k_{i},k)$ for other measurement matrices can be also easily determined.

The reduced set of measurements (samples) manifests as a noise in the initial estimate, which therefore acts as a random variable, with mean-value and variance given by

[TABLE]

where $\delta(k)=1$ only for $k=0$ and $\delta(k)=0$, elsewhere.

In the analysis of the reconstruction error, we are interested in the variance of random variable $\mu(k_{i},k)$, i.e.

[TABLE]

For the partial DFT matrix, the variance is derived in [6]. For a real-valued ETF measurement matrix, the values $\pm\mu$ are equally probable, producing the variance $\sigma_{\mu}^{2}=\mu^{2}$, where, according to the Welch bound, $\mu^{2}=(N-M)/(M(N-1))$ holds [27, 13]. For the Gaussian measurement matrix, the variance is $\sigma_{\mu}^{2}=1/M$. The same value is obtained for other considered random matrices. The variance $\sigma_{\mu}^{2}$ of $\mu(k_{i},k)$ is presented in Table I for various measurement matrices [6, 13, 15, 14].

III-B Detection of Nonzero Element Positions

The initial estimate can be used as a starting point for an analysis of the reconstruction performance and its outcomes. Potentially, such analysis can lead to the improvements of the reconstruction process. The detection can be done in one step or in an iterative way.

One-step detection: In an ideal case, matrix $\mathbf{A}^{H}\mathbf{A}$ should be such that the initial estimate $\mathbf{X}_{0}$ contains $K$ coefficients higher than the other coefficients. Then by taking the positions of the highest coefficients in (7) as the set $\mathbb{K}$, the signal is simply reconstructed using (6).

Iterative detection: The condition that all $K$ nonzero coefficients in the initial estimate $\mathbf{X}_{0}$ are larger than the coefficient values $X_{0}(k)$ at the original zero-valued positions $k\notin\mathbb{K}$, can be relaxed using an iterative procedure. To find the position of the largest coefficient in $\mathbf{X}$ based on $\mathbf{X}_{0}$, it is sufficient that the coefficient $X_{0}(k)$ has a value larger than the values of the coefficients $X_{0}(k)$ at the original zero-valued coefficient positions $k\notin\mathbb{K}$.

Remark 1:  Solution uniqueness. The worst case for the detection of a nonzero coefficient, with a normalized amplitude $1$, occurs when the remaining $K-1$ coefficients are equally strong (i.e. with unity amplitudes). Then, the influence of other nonzero coefficients to the initial estimate of the considered coefficient may assume its highest possible value. The influence of the $k$th coefficient on the one at the $i$th position is equal to $\mu(k_{i},k)$, given by (10). Its maximum possible absolute value is the coherence index $\mu$. In the worst case, the amplitude of the considered coefficient in the initial estimate is $1-(K-1)\mu$. At the position where there the original coefficient $X(k)$ is zero-valued, in the worst case, the maximum possible contributions $\mu$ of all $K$ coefficients sum up in phase to produce the maximum possible disturbance $K\mu$. The detection of the strongest coefficient is successful if

[TABLE]

producing the well-known coherence condition for the unique reconstruction $K<\left(1+1/\mu\right)/2.$

After the largest coefficient position is found and its value is estimated, this coefficient can be subtracted and the procedure can be continued with the remaining $(K-1)$-sparse signal. If the reconstruction condition is met for the $K$-sparse signal, then it is met for all lower sparsities as well.

The procedure is iteratively repeated for each coefficient. The stopping criterion is that $\mathbf{A}_{MK}\mathbf{X}_{K}=\mathbf{y}$ holds for the estimated positions $\{k_{1},k_{2},\dots,k_{K}\}$ and coefficients $X(k)$.

The method is summarized in Algorithm 1.

Remark 2:  Solution exactness. The coherence index condition guarantees that the positions of the nonzero elements in $\mathbf{{X}}$ will be uniquely determined. Next, we will show that the values of nonzero coefficients will be recovered exactly.

The system of linear equations in (9), for $k\in\mathbb{K}$, can be written in a matrix form as

[TABLE]

where $\mathbf{B}$ is a $K\times K$ matrix with elements $b_{pi}=\mu(k_{p},k_{i})$, $\mathbf{{X}}_{0K}$ is the vector with $K$ elements obtained from the initial estimate as ${X}_{0K}(i)={X}_{0}(k_{i})$, and $\mathbf{X}_{K}$ is the vector with $K$ corresponding coefficients from the original signal. The influence of the other $K-1$ coefficient to the considered coefficient is denoted by $\mathbf{C}_{K}$.

The reconstructed coefficients $\mathbf{X}_{R}$, at the nonzero coefficient positions, are obtained by minimizing $\left\|{\mathbf{y}-\mathbf{A}_{K}\mathbf{X}_{R}}\right\|_{2}^{2}$. They are

[TABLE]

where $\mathbf{A}_{K}$ is a matrix obtained from measurement matrix $\mathbf{A}$ by keeping the columns for $k\in\mathbb{K}$. Since $\mathbf{A}_{K}^{H}\mathbf{y}=\mathbf{{X}}_{0K}$, according to (7), we can rewrite (14) as

[TABLE]

Since $\mathbf{{X}}_{0K}=\mathbf{B}\mathbf{X}_{K}$, the reconstruction is exact if

[TABLE]

holds. Indeed, the elements of matrix $\mathbf{A}_{K}^{H}\mathbf{A}_{K}$ are equal to $\beta_{{i}{j}}=\sum_{n=1}^{M}a^{*}_{n}(k_{i})a_{n}(k_{j})=\mu(k_{j},k_{i})$ meaning that $\mathbf{A}_{K}^{H}\mathbf{A}_{K}=\mathbf{B}$. Therefore, $\mathbf{X}_{R}=\mathbf{X}_{K}$ holds.

The reconstruction algorithm produces correct coefficient values $X(k)$ at the selected positions $k\in\mathbb{K}$. It means that the influence of other $K-1$ coefficients to each coefficient in the initial coefficient estimate ${X}_{0}(k),$ denoted by $C(k),$ is canceled out.

In summary, the reconstruction algorithm for a coefficient at a position $k\in\mathbb{K}$, works as an identity system to the original signal coefficient in $X_{0}(k),$ eliminating the influence of other coefficient at the same time, Fig. 1.

III-C Noisy Measurements

Assume next that the observations are noisy

[TABLE]

with a zero-mean signal independent noise $\boldsymbol{\varepsilon}$. The noise variance of the assumed additive input noise $\boldsymbol{\varepsilon}$ is $\sigma^{2}_{\varepsilon}$ and the covariance is given by

[TABLE]

Noisy measurements will result in a noisy estimate $\mathbf{X}_{0}=\mathbf{A}^{H}(\mathbf{y}+\boldsymbol{\varepsilon})$. Variance of $X_{0}(k)$ due to the input noise in measurements, is $\sigma^{2}_{{X}_{0}(k)}=\sigma^{2}_{\varepsilon}$ since it has been assumed that the columns of $\mathbf{A}$ have unite energy,

[TABLE]

The noise variance in the reconstructed coefficient is (Remark 1 and Fig. 1)

[TABLE]

Since the noise is independent in each reconstructed coefficient, the total mean squared error (MSE) in $K$ reconstructed coefficients is

[TABLE]

If the partial DFT matrix is formed as a submatrix of the standard inverse DFT matrix (with normalization $1/N$), then we would get $\left\|{{\mathbf{X}_{R}}-{\mathbf{X}_{K}}}\right\|_{2}^{2}=KN^{2}{\sigma_{\varepsilon}}^{2}/M$, as shown, for example, in [15].

IV Quantization Effects

Traditional CS theory does not consider the limitations in the number of bits used for the measurements representation. This can affect the reconstruction performance of the standard CS approaches.

The measurement quantization is particularly important in the hardware implementation context. One-bit measurements are the most extreme case, promising simple, comparator-based hardware devices [16]. The one bit used represents the sign of the sample, i.e. $\mathbf{y}=\text{sign}\{\mathbf{AX}\}$. However, a larger number of samples is required for an accurate reconstruction, which is difficult to achieve using only the sign of a measurement [16, 17].

A more general form of the hardware implementation uses a $B$-bit digital sample of a measurement. We will assume that the measurements are stored into $(B+1)$-bit registers (one sign bit and $B$ bits for the signal absolute value),

[TABLE]

whereas the reconstruction of the coefficients $X(k)$ is done in a more realistic sense for hardware purposes. The requirement for storage is also significantly reduced for such measurements, since the total number of bits is reduced. Note that, for a complex-valued signal $x(n)$, the measurements $\mathbf{y}_{B}$ are also complex, formed as

[TABLE]

where both real and imaginary part of measurements are quantized to $B$-bits.

IV-A Quantization errors

Quantization influences the results of the compressive sensing reconstruction in several ways:

• •

Input signal quantization error, described by an additive quantization noise. This influence can be modeled as a uniform noise with values between the quantization level bounds.

• •

Quantization of the results of arithmetic operations. It depends on the way how the calculations are performed.

• •

Quantization of the coefficients in the algorithm. However, being deterministic for a given measurement matrix, this type of error is commonly neglected from the analysis.

In order to perform an appropriate and exact analysis, some standard assumptions are further made:

• •

The quantization error is a white noise process with a uniform distribution.

• •

The quantization errors are uncorrelated.

• •

The quantization errors are not correlated with the input signal.

The most important source of an error is the quantization of the measurements $y(m)$ and the quantization of the measured sparse signals $X(k)$, referred to as the quantization noise folding. They will be analyzed next.

IV-B Input signal ranges

Assume that registers with $B$ bits, with an additional sign bit, are used and that all measurements are normalized to the range

[TABLE]

The total number of bits in a register is $b=B+1$.

In that case, it is important to notice that the sparse signal coefficients $X(k)$ must be within the range $-\min\{\sqrt{M}/K,1\}<X(k)<\min\{\sqrt{M}/K,1\}$ so that $\mathbf{y}=\mathbf{A}_{MK}\mathbf{X}_{K}$ does not produce a value with amplitude greater or equal to $1$. For the partial DFT matrices, this condition holds in a strict sense, while for the Gaussian matrices it holds in a mean sense (all values whose amplitudes are greater than $1$ are quantized to the closest level with amplitude below $1$). Note that the butterfly schemes for the measurements calculation (as in the quantized FFT algorithms) could extend this bounds for $X(k)$ so that the maximum range $-1<X(k)<1$ can be used.

IV-C Measurements quantization

For the $B$-bit registers, the digital signal values $\mathbf{y}_{B}$ are coded into a binary format. When the signal amplitude is quantized to $B$ bits, the difference in amplitude which produces the quantization is called the quantization error. The quantization error is bounded by

[TABLE]

where $\Delta$ is related to $B$ through

[TABLE]

The quantization error of a signal can be defined as an additive uniform white noise affecting the measurements

[TABLE]

where $\mathbf{e}$ is the quantization error vector with elements $e(m)$. The mean and variance of the quantization noise are calculated as [13]

[TABLE]

Note that, for a complex-valued signal, both real and imaginary parts of samples contribute to the noise. Therefore, in this case, the variance of the quantization noise can be written as

[TABLE]

Considering $\mathbf{y}$ as noisy measurements, the initial estimate will result in a noisy $X_{0}(k)$. Since $X_{0}(k)$ is calculated from (8), with the quantization noise in measurements, its variance will be

[TABLE]

Therefore, the noise variance in the output (reconstructed) coefficients, for the system shown in Fig. 1, is equal to the input noise variance

[TABLE]

Since only $K$ out of $N$ coefficients are used in the reconstruction, the energy of the reconstruction error is

[TABLE]

where for notation simplicity we have used $\left\|{{\mathbf{X}_{R}}-{\mathbf{X}_{K}}}\right\|_{2}^{2}$ to denote the expected value of the squared norm-two of the vector ${\mathbf{X}_{R}}-{\mathbf{X}_{K}}$. The full and complete notation of the left side of (28) would be $\mathrm{E}\{\left\|{{\mathbf{X}_{R}}-{\mathbf{X}_{K}}}\right\|_{2}^{2}\}$.

IV-D Sparsity to Number of Bits Relation

Based on the previous relations, influence of the quantization with $B$ bits can be related to the sparsity $K$. The error energy in the reconstructed coefficients will remain the same if

[TABLE]

It means that reduction of $B$ to $B-1$ will require sparsity reduction from $K$ to $K/4$. The logarithmic form of the reconstruction error is

[TABLE]

V Nonsparsity Influence

Due to many circumstances, majority of signals in real-world scenarios are only approximately sparse or nonsparse. This means that a signal, in addition to the $K$-sparse large coefficients, has $N-K$ coefficients in the sparsity domain which are small but nonzero. Assume such an approximately sparse (or nonsparse) signal $\mathbf{X}$. The signal is reconstructed under the $K$-sparsity constraint using Algorithm 1, with the reconstruction conditions being satisfied in the compressive sensing sense, thus allowing that the algorithm can detect the $K$ largest coefficients.

The reconstructed signal $\mathbf{X}_{R}$ then has $K$ reconstructed coefficients with amplitudes $X_{R}({k_{1}}),\ X_{R}({k_{2}}),\dots,X_{R}({k_{K}})$. The remaining $N-K$ coefficients which are not reconstructed are treated as a noise in these $K$ largest coefficients. Variance from a nonzero coefficient, according to (13), is $\left|X(k_{l})\right|^{2}(N-M)/(M(N-1)).$ The total energy of noise in the $K$ reconstructed coefficients $\mathbf{X}_{R}$ will be

[TABLE]

where $\mathbf{X}_{K}$ is the sparse version of the original (nonsparse) signal, i.e. a signal with $K$ largest coefficients from $\mathbf{X}$, and others set to zero. Denoting the energy of remaining signal, when the $K$ largest coefficients are removed from the original signal, by

[TABLE]

we get

[TABLE]

For the partial DFT measurement matrix, the result will be

[TABLE]

In the case when the signal is strictly $K$-sparse, i.e. $\mathbf{X=X}_{K}$, and if the reconstruction is performed with non-quantized measurements, the reconstruction would be ideal and the error would be $\left\|\mathbf{X}_{R}\mathbf{-X}_{K}\right\|_{2}^{2}=0$ (or negligible). Since the measurements are quantized by $B$-bits, error of the form (28) will be introduced.

In the case of a nonsparse signal, a general expression is obtained by combining (28) and (32) to get

[TABLE]

This result will be validated by examples in the next section, by calculating the signal-to-noise ratio (SNR) of each result

[TABLE]

and comparing it with the statistical SNR given by

[TABLE]

VI Noise Folding Quantization

Additionally to the nonsparse case, we extend our analysis to include the case when a quantization noise $\mathbf{z}$ exists in the signal coefficients $\mathbf{X}$, prior to taking the measurements [26]. In this case, the measurements are of the form

[TABLE]

which can be rewritten as

[TABLE]

where $\mathbf{v}=\mathbf{e}-\mathbf{Az}$, and the total quantization noise affecting the signal measurements is denoted by $\mathbf{e}$ with covariance $\sigma_{\mathbf{e}}^{2}\mathbf{I}.$ The quantization noise vector $\mathbf{z}$ is random with covariance $\sigma_{\mathbf{z}}^{2}\mathbf{I}$, being independent of $\mathbf{e}$. Therefore, the resulting noise $\mathbf{v}$ is characterized by a covariance matrix

[TABLE]

If the considered measurement matrix $\mathbf{A}$ is formed as a partial Fourier matrix, the relation $\mathbf{AA}^{H}=\frac{N}{M}\mathbf{I}$ holds. The variance of $\mathbf{v}$ is then

[TABLE]

with the covariance matrix $\mathbf{C}=\sigma_{\mathbf{v}}^{2}\mathbf{I}$.

However, for the sparse case, the quantization error is present in only $K$ nonzero elements of $\mathbf{X}$. It means that the noise $\mathbf{A}\mathbf{z}$ variance is $\frac{K}{M}\sigma_{\mathbf{e}}^{2}$ or

[TABLE]

For the nonsparse partial DFT matrix case, $K(N-M)/(M(N-1))\left\|\mathbf{X-X}_{K}\right\|_{2}^{2}$ is added to the right part of the Eq. (41)

[TABLE]

where it is assumed that the quantization of $K$ largest elements in $\mathbf{X}$ is dominant in that part of error. It is shown that all previous relations, for various measurement matrices, can be applied to this case.

VII Floating Point Registers

In floating point registers, the quantization error is modeled as a multiplicative error

[TABLE]

where $\mathbf{e}$ is the quantization error vector with elements $e(n)$. As in the classical digital signal processing, for the analysis of floating point arithmetics, it will be assumed that the sparse coefficients $X(k_{i})$, $i=1,2,\dots,K$, are independent, equally distributed zero-mean random variables, with variance $\sigma^{2}_{X}$. The coefficients $X(k_{i})$ are statistically independent from the measurement matrix $\mathbf{A}$ elements $a_{m}(k)$. The mean value of the quantization error is

[TABLE]

The variance is

[TABLE]

For $y(m)=\sum_{i=1}^{K}X(k_{i})a_{m}(k_{i})$ we get

[TABLE]

for all measurement matrices with normalized energy columns, when their elements $a_{m}(k)$ are equally distributed. This means that the quantization noise $y(m)e(m)$ has the variance $\sigma^{2}_{X}\sigma^{2}_{\mathbf{e}}K/M$ and we can write (Remark 1 and Fig. 1)

[TABLE]

with

[TABLE]

All formulas, in various considered scenarios, can now be rewritten, including the cases of nonsparse signals and noise folding.

For example, if the measurements are normalized such that $\mathrm{E}\{y^{2}(m)\}=\sigma^{2}_{X}K/M=1$ then $\left\|{{\mathbf{X}_{R}}-{\mathbf{X}_{K}}}\right\|_{2}^{2}=K\sigma^{2}_{\mathbf{e}},$ that is the floating-point arithmetics produces the same results as the fixed-point arithmetics. However, if the range of measurement values is lower, for example, $\mathrm{E}\{y^{2}(m)\}=\sigma^{2}_{X}K/M=1/10$, then the floating-point arithmetics will produce ten times lower error, $\left\|\mathbf{X}_{R}-\mathbf{X}_{K}\right\|_{2}^{2}=K\sigma^{2}_{\mathbf{e}}/10$.

VIII Numerical Results

Example 1: One realization of a sparse and nonsparse signal will be considered as an illustration of the reconstruction.

a) Consider an $N=256$-dimensional signal of sparsity $K=10$, whose $M=N/2$ available measurements are stored into registers with $B=6$ bits. The measurements matrix is a partial DFT matrix with randomly selected $M$ out of $N$ rows from the full DFT matrix, with columns being energy normalized. The the sparsity domain coefficients are assumed in the form

[TABLE]

where $\nu(p)$ is a random variable with uniform distribution from [math] to $0.4$.

Since this signal is sparse, the reconstruction error is defined by (28). The SNR is defined by (35) with $\left\|\mathbf{X-X}_{K}\right\|_{2}^{2}=0$. The original and the reconstructed signals are shown in Fig. 2(top). The statistical SNR is $SNR_{st}=42.35$ dB and $SNR_{th}=42.56$ dB.

b) The signal from a), with $K=10$ significant coefficients, is considered here. However, we will also assume that the remaining $N-K$ coefficients are small but not zero-valued,

[TABLE]

with $\nu(p)$ being a random variable with uniform distribution from [math] to $0.4$ as in a) and $N=256$. The number of bits in the registers where the measurements are stored is $B=6$. The original signal and the signal reconstructed under sparsity constraint, using $M=N/2$ measurements, are shown in Fig. 2(bottom). The statistical SNR is $SNR_{st}=33.33$ dB and the theoretical is $SNR_{th}=33.68$ dB.

Example 2: Statistical analysis of the sparse signal reconstruction whose form is given in (46) is performed in this example. Random uniform changes of coefficient amplitudes $\nu(p)$ are assumed from [math] to $0.2$. The numbers of quantized measurements $M=192$, $M=170$, and $M=128$ are considered. Typical cases for the measurements quantization to $B\in\{4,6,8,10,12,14,16,18,20,24\}$ bits are analyzed.

Signals with sparsity levels $K\in\{5,10,15,20,25,30\}$ are considered. The average statistical signal-to-nose ration $SNR_{st}$ and theoretical signal-to-noise ration $SNR_{th}$ values in 100 realizations are presented in Fig. 3(a)-(c). Black dots represent the statistical results, $SNR_{st}$, and the dash-dot lines show the theoretical results, $SNR_{th}$. The agreement is high.

For nonsparse signals we used the model in (47). Random changes of coefficient amplitudes $\nu(p)$ are assumed from [math] to $0.2$, while the amplitudes of the coefficients $X(k)$ for $k_{p}\notin\mathbb{K}$ are of the form $X(k_{p})=\exp(-p/(8K))$ in order to reduce its influence to the quantization level. With such amplitudes of the nonsparse coefficients, the quantization error dominates in the reconstruction up to $B=14$, while the nonsparse energy is dominant for $B\geq 16$, as can be seen in Fig. 3(d)-(f). Statistics is again in full agreement with the theoretical results.

Finally, the noise folding effect is included, taking into account that the input coefficients $X(k)$ are quantized, in addition to the quantization of measurements $y(m)$. Since the folding part of quantization error is multiplied by $K/M\ll 1$ in (41), the results do not differ from those presented in Fig. 3(a)-(c). In order to test the influence of noise folding we assumed that the quantized input coefficients $X(k)$ contain an additional noisy. An input additive complex-valued i.i.d. Gaussian noise with variance $\sigma_{z}=0.0001$ is added to these coefficients. This noise is of such a level that it does not influence quantization error for $B<14$. However, for $B\geq 14$, it becomes larger than the quantization error and its influence becomes dominant. The results with the quantization and noise folding, with additional noise, are shown in Fig. 3(g)-(i).

Example 3: The statistical analysis is extended to other forms of the measurement matrices, namely the ETF, the Gaussian, and the uniform random matrix. All three forms of the signal and quatization error are considered here with $M=128$ measurements. Sparse and nonsparse signals described in Example 2 are used in the analysis. The reconstruction error with various number of bits $B=\{4,6,8,10,12,14,16,18,20,24\}$ used in the quantization, and various assumed sparsity levels $K$ are shown in Fig. 4(a)-(c). The results for nonsparse signals, reconstructed with assumed sparsity, are presented in Fig. 4(d)-(f). The noise folding is analyzed for a reduced number of bits in the quantization of $X(k)$ and presented in Fig. 4(g)-(i).

Example 4: The analysis of quantization effects is done with the assumption that the quantization errors are uncorrelated. This condition is met for all previously considered matrices. However, in the case of Bernoulli measurement matrix and a small signal sparsity this condition does not hold, meaning that we can not expect quite accurate estimation of the statistical error using the previous formulas. To explain this effect, we will start with the simplest case of signal whose sparsity is $K=1$. The measurements are $y(m)=a_{m}(k_{1})X(k_{1})$. For all previously considered matrices $y(m)$ and $y(n)$ are different for $m\neq n$ and the quantization errors are independent. However, for the Bernoulli measurement matrix we have $y(m)=\pm X(k_{1})/\sqrt{M}$. These measurements will produce only two possible quantization errors for all $m=1,2,\dots,M$. It means that $M/2$ errors in the initial estimate will sum up in phase, producing the mean square error with variance $\operatorname{var}\{X_{R}(k)\}=\frac{M}{2}\sigma^{2}_{\mathbf{e}}$. This is significantly higher than $\operatorname{var}\{X_{R}(k)\}=\sigma^{2}_{\mathbf{e}}$ in other cases. For $K=2$, we get the measurements $y(m)=(\pm X(k_{1})\pm X(k_{2}))/\sqrt{M}$, producing four possible values for $y(m)$ and only four possible values of the quantization error. For large $K$, the number of possible levels increases and the result for the variance converges to the one for uncorrelated quantization errors $\operatorname{var}\{X_{R}(k)\}=\sigma^{2}_{\mathbf{e}}$, obtained under the assumption that all $M$ measurements $y(m)$ are different. The results for the Bernoulli measurement matrix, with the described correstion for small $K$, are shown in Fig. 5.

Example 5: The result in previous four examples are obtained based on the OMP reconstruction method (Algorithm 1), which is also used for the derivation of the theoretical results. Here we will show that we may expect similar results for other reconstruction methods as far as the reconstruction conditions are met. The simulation for the reconstruction with the partial DFT measurement matrix with sparse and nonsparse signals, including noise folding, are repeated with the iterative hard thresholding (IHT) reconstruction method, given in Algorithm 2 [12, 24]. The theoretical and statistical errors are shown in Fig. 6(a)-(c), showing high agreement between the statistical and theoretical results.

Example 6: In this example, the reconstruction of sparse and nonsparse signals, including noise folding, is performed using the Gaussian measurement matrix and the Bayesian-based method, [28, 29], summarized in Algorithm 3. The theoretical and statistical errors for the signals from Example 1 are shown in Fig. 7(a)-(c). Again, a high agreement between the statistical and theoretical results is obtained.

IX Conclusions

The effect of quantization noise to the reconstruction of signals under sparsity constraint is analyzed in this paper. If the measurements are not quantized, the reconstruction would be ideal and the error will be negligible. However, the quantization is a very useful step since the hardware realization of systems cannot store the exact values of samples. We derived the exact error of the reconstruction due to quantization noise. The cases when a signal is not strictly sparse is analyzed, as well as the noise folding effect. The reconstruction performance is validated on numerical examples, and compared to the statistical error calculation concluding a high agreement between them.

Influence of quantization to the reconstruction condition

Quantization noise can be included in the coherence index based relation for the reconstruction. The worst case amplitude of the considered normalized coefficient in the initial estimate is $1-(K-1)\mu-\nu\Delta/2$ where

[TABLE]

This inequality follows from the relation between the norm-two and the norm-one of a vector. For the partial DFT matrix, the random partial DFT matrix, and the Bernoulli matrix the equality $\nu=\sqrt{M}$ holds. Following the same reasoning as in Remark III-B, we may conclude that at a position where the the ordinal coefficient $X(k)$ is zero-valued, in the worst case, the maximum possible disturbance is $K\mu+\sqrt{M}\Delta/2$. The detection of the strongest coefficient position is always successful if $1-(K-1)\mu-\sqrt{M}\Delta/2>K\mu+\sqrt{M}\Delta/2$, producing the condition for reconstruction

[TABLE]

Influence of the quantization error to the uniqueness condition will be negligible if $\sqrt{M}\Delta\ll 1$ holds.

Coherence index based values guarantee exact reconstruction, however, they are pessimistic. More practical relations can be obtained by considering the probabilistic analysis [6, 30]. The resulting disturbance in the initial estimate at the position $k\in\mathbb{K}$ due to the other coefficients and quantization noise behaves as a Gaussian random variable $\mathcal{N}(1,(K-1)\sigma_{\mu}^{2}+\sigma^{2}_{\mathbf{e}})$, for $K\gg 1$ . The initial estimate at $k\notin\mathbb{K}$ behaves as $\mathcal{N}(0,K\sigma_{\mu}^{2}+\sigma^{2}_{\mathbf{e}})$. Probabilistic analysis may provide approximative relations among $N$, $M$, and $K$, for a given probability. We have performed the statistical analysis with various measuremet matrices. The results of this analysis lead to the conclusion that for high probabilities of reconstruction we may neglect the quantization effect influence to the reconstruction condition for $B\geq 4$.

Note that for large sparsities $K$ we have found that the reconstruction probability can be improved by increasing the upper limit for iterations in Algorithm 1 for a few percents, with respect to the expected sparsity $K$. After the iterations are completed, the expected sparsity $K$ is used in the final reconstruction. This solves the problem that the iterative reconstruction in Algorithm 1 cannot produce the exact result if it misses one of the nonzero coefficient positions during the iterative process for large $K$.