A New Generalized Projection and Its Application to Acceleration of Audio Declipping

Abstract

In convex optimization, it is often inevitable to work with projectors onto convex sets composed with a linear operator. Such a need arises from both the theory and applications, with signal processing being a prominent and broad field where convex optimization has been used recently. In this article, a novel projector is presented, which generalizes previous results in that it admits to work with a broader family of linear transforms when compared with the state of the art but, on the other hand, it is limited to box-type convex sets in the transformed domain. The new projector is described by an explicit formula, which makes it simple to implement and requires a low computational cost. The projector is interpreted within the framework of the so-called proximal splitting theory. The convenience of the new projector is demonstrated on an example from signal processing, where it was possible to speed up the convergence of a signal declipping algorithm by a factor of more than two.

1. Introduction

Proximal algorithms are a modern branch of mathematical optimization, whose principal building blocks are the so-called proximal operators [1,2,3,4]. In applications, the model of the problem often involves the composition of a non-linear function with a linear operator [5,6,7,8,9,10]. When proximal optimization is utilized to find approximate solution iteratively, the proximal operator of such a composition is needed, making the resulting algorithms computationally expensive in general [11]. However, there are special cases in which the proximal operator can be expressed explicitly [2,12,13,14].

In the paper, we introduce and prove a new lemma which provides an explicit formula for the evaluation of the projector onto a convex set, composed with a linear operator. It is motivated by an already known property of proximal operators of functions composed with semi-orthogonal linear operators [2]. Actually, our lemma generalizes the previous result in two directions, since it admits working in a complex vector space and the linear operator is only assumed to satisfy a certain diagonality property. On the other hand, the new projection formula is limited to projectors onto box-type sets. Since projectors represent a special case of proximal operators, our lemma falls into the general concept of proximal minimization; nevertheless, our result is valuable per se, since projections are used in many other contexts and fields of science.

As an example, we apply the novel projection to the problem of signal declipping, i.e., to the restoration of unknown parts of a signal degraded by clipping/saturation, which is a non-linear distortion of an audio signal, described in detail in the experimental part. The declipping model is quite simple, since it is not the goal of this paper to develop a method that outperforms the state of the art, but rather to illustrate the use of the proposed lemma in a technical application arising from the field of audio signal recovery.

In Section 2, the proximal and projection operators are reviewed and the state of the art in this specific area is presented. Next, the novel generalized projection is introduced, including the proof. Section 3 then discusses the result and emphasizes its advantages over the current methods, and also points out its limitations. Based on the new result, an example experiment in the signal processing field is presented in Section 4, where we show that audio declipping benefits from our approach by the computation being accelerated by a factor of at least two.

2. Methods

2.1. Proximal Algorithms

The proximal algorithm is an iterative algorithm that, under certain conditions, provides a sequence of vectors $\left\{{\mathbf{x}}_k\right\}$ converging to the minimizer of a sum of convex functions ${\sum }_i{f}_i\left(\mathbf{x}\right)$. This is achieved by sequential evaluation of the so-called proximal operators associated with each of the summands, $f_i$. The most often used and studied are proximal algorithms that can minimize the sum of two convex functions, $f_1\left(\mathbf{x}\right)+f_2\left(\mathbf{x}\right)$, such as the Forward–backward algorithm or the Douglas–Rachford algorithm [13,15]; however, there exist algorithms for minimizing the sum of arbitrarily many functions [2].

Algorithms that can minimize a sum involving a linear operator have been presented only recently. Chambolle and Pock [5] presented a primal–dual algorithm that minimizes the sum of two functions $f_1\left(\mathbf{x}\right)+f_2\left(L\mathbf{x}\right)$ where neither $f_1$ nor $f_2$ has to be smooth. The Forward–backward based primal—dual (FBB-PD) algorithm [6] can cope with the sum of three functions $f_1\left(\mathbf{x}\right)+f_2\left(\mathbf{x}\right)+f_3\left(L\mathbf{x}\right)$, where $f_2$ is smooth. The most general algorithm published by Condat [7,11] admits the minimization of any finite number of functions with the inclusion of potentially multiple linear operators, $f_1\left(\mathbf{x}\right)+f_2\left(\mathbf{x}\right)+f_3\left(L_3\mathbf{x}\right)+\dots +f_K\left(L_K\mathbf{x}\right)$.

2.2. Proximal Operators

Let $\tilde{\mathbb{R}}$ denote the extended real line, i.e., $\tilde{\mathbb{R}}=\mathbb{R}\cup \{-\infty ,\infty \}$. Proximal operators are the key components in proximal algorithms. Recall that the proximal operator of a convex function $h:{\mathbb{C}}^N\to \tilde{\mathbb{R}}$ maps $\mathbf{z}\in {\mathbb{C}}^N$ to another vector in ${\mathbb{C}}^N$, such that [3,16]

$${\mathrm{prox}}_h\left(\mathbf{z}\right)=\underset{\mathbf{x}\in {\mathbb{C}}^N}{arg~min}\left\{h\left(\mathbf{x}\right)+\frac{1}{2}{\parallel \mathbf{z}-\mathbf{x}\parallel }_2^2\right\}.$$

(1)

The proximal operator is a generalization of a projection operator onto a convex set—this fact is clear if we identify h with the indicator function ${\iota }_K$ of such a set K; in such a case [1],

$$h\left(x\right)={\iota }_K\left(x\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{for}x\in K\hfill \\ \infty \hfill & \mathrm{for}x\notin K,\hfill \end{array}\right.$$

and we get the relation ${\mathrm{proj}}_K={\mathrm{prox}}_{{\iota }_K}$. For a convex h, the map in Equation (1) is uniquely defined [2]. A characterization of proximal operators of convex lower semicontinuous functions is the work of Moreau [16].

In this article, we are interested in the proximal operator of the map $f(L·)=f\circ L$. Let us denote such a composition $h\left(\mathbf{x}\right)=(f\circ L)\left(\mathbf{x}\right)=f\left(L\mathbf{x}\right)$, with $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$. If f is convex, then $f\circ L$ is also convex, which is straightforward to show. Therefore, in line with Equation (1), our goal is to find

$${\mathrm{prox}}_{f\circ L}\left(\mathbf{z}\right)=\underset{\mathbf{x}\in {\mathbb{R}}^N}{arg~min}\left\{f\left(L\mathbf{x}\right)+\frac{1}{2}{\parallel \mathbf{z}-\mathbf{x}\parallel }_2^2\right\}.$$

(2)

Although explicit formulas exist for many particular choices of h in Equation (1) (see, e.g., [2]), the situation gets complicated for the compositions in Equation (2). In addition, note, that when ${\mathrm{prox}}_h$ is available, then ${\mathrm{prox}}_{\nu h}$ is often easy to find for $\nu >0$.

In many practical situations, one can benefit from the following lemma, occurring, for example, in ([2], Table 10.1.x). Let $\mathit{Id}$ denote the identity operator, i.e., $\mathit{Id}\left(\mathbf{x}\right)=\mathbf{x}$; the optional subscript indicates the dimensionality of $\mathbf{x}$.

Lemma 1.

Let L be a linear operator, $L:{\mathbb{R}}^N\to {\mathbb{R}}^M$, such that $L{L}^{\top }=\nu {\mathit{Id}}_M,$$\nu >0$. Then,

$${\mathrm{prox}}_{f\circ L}={\mathit{Id}}_N+{\nu }^{-1}{L}^{\top }\circ ({\mathrm{prox}}_{\nu f}-{\mathit{Id}}_M)\circ L.$$

(3)

$L^{\top }$ denotes the adjoint of L; the symbol ${}^{\top }$ is used here to emphasize that L is assumed to be real. In the finite-dimensional case, $L^{\top }$ can be understood as the transpose of the $M\times N$ matrix identified with L. The property described by the above lemma can be found in [3,12,13,14], exhibiting different styles of proof but limited to real space setting. Equation (3) says that, to evaluate the composite proximal operator, it is sufficient to know the proximal operator of $\nu f$ and to be able to perform L and $L^{\top }$. For an arbitrary $\mathbf{z}\in {\mathbb{R}}^N$, this means

$${\mathrm{prox}}_{f\circ L}\left(\mathbf{z}\right)=\mathbf{z}+{\nu }^{-1}{L}^{\top }({\mathrm{prox}}_{\nu f}\left(L\mathbf{z}\right)-L\mathbf{z}).$$

(4)

The assumption $L{L}^{\top }=\nu \mathit{Id}$ is called “semi-orthogonality” of L in [2]; from the viewpoint of frame theory, such an L corresponds to a synthesis operator of a so-called tight frame with the frame constant $\nu$ [17,18,19,20,21]. Note that as a consequence of the assumption, it automatically holds $M\le N$ due to $M=\mathrm{rank}\left(L{L}^{\top }\right)=\mathrm{rank}\left(L\right)\le N$. Equation (3) drastically reduces the complexity of computations when ${\mathrm{prox}}_f$ is known, since Equation (2) is then solved explicitly when compared with the general case.

Note that a tight frame synthesis operator L acts the same way as an orthonormal operator does, since both satisfy $L{L}^{\top }=\nu \mathit{Id}$ ($\nu =1$ in both the orthonormal case and the case of Parseval tight frame [17]). Expressed in words, the analysis followed by the synthesis reproduces the object. Observe that the opposite $L^{\top }L=\nu \mathit{Id}$ is true if and only if L is orthonormal (unitary) operator ${\mathbb{R}}^N\to {\mathbb{R}}^N$ as a special case of Parseval tight frame with $\nu =1$ and $M=N$.

Example 1.

As an example, the projection onto the set $\left\{\mathbf{z}\right|\parallel L\mathbf{z}-\mathbf{y}{\parallel }_2\le \delta \}$ has to be computed iteratively in general; however, when L is a synthesis operator of a tight frame, such a projection can be evaluated explicitly by

$${\mathrm{proj}}_{\left\{\mathbf{z}\right|\parallel L\mathbf{z}-\mathbf{y}{\parallel }_2\le \delta \}}\left(\mathbf{x}\right)=\mathbf{x}+{\nu }^{-1}{L}^{\top }\left({\mathrm{proj}}_{\left\{\mathbf{z}\right|\parallel \mathbf{z}-\mathbf{y}{\parallel }_2\le \delta \}}\left(L\mathbf{x}\right)-L\mathbf{x}\right),$$

(5)

i.e., using a simple projection onto the ball $\left\{\mathbf{z}\right|\parallel \mathbf{z}-\mathbf{y}{\parallel }_2\le \delta \}$. Note that the set $\left\{\mathbf{z}\right|\parallel L\mathbf{z}-\mathbf{y}{\parallel }_2\le \delta \}$ represents an ellipsoid in ${\mathbb{R}}^N$ when bijective $L:{\mathbb{R}}^N\to {\mathbb{R}}^N$ is considered.

At this point, we recall the pseudoinverse operator, which is used in the rest of the paper. Lemmas 2 and 3 are discussed, for example, in ([22], Section 5.13).

Lemma 2.

Let $L:{\mathbb{C}}^N\to {\mathbb{C}}^M,M\le N$ be a surjective linear operator and $\mathbf{y}\in {\mathbb{R}}^M$. Then, $L^{+}\mathbf{y}$ is the least-norm solution to the system $L\mathbf{x}=\mathbf{y}$, i.e., $\parallel L^{+}{\mathbf{y}\parallel }_2\le {\parallel \mathbf{x}\parallel }_2$, for any $\mathbf{x}$ satisfying $L\mathbf{x}=\mathbf{y}$. The pseudoinverse operator $L^{+}$ is in this surjective case defined through $L^{+}=L^{*}{\left(L{L}^{*}\right)}^{-1}$.

The application of $L^{+}L$ corresponds to projecting onto the range space of $L^{\top }$. Another, more common case of projection is the projection onto an affine subspace. The motivation for the related following lemma is twofold: First, the described projection belongs to the “family” of proximal operators we are concerned with. Second, this result is utilized below in the experiment. To put it into context, this lemma holds for any complex $\mathbf{b}$, but, on the other hand, its application is limited to affine spaces only.

Lemma 3.

Let $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$ be any surjective linear operator and $\mathbf{b}\in {\mathbb{C}}^M$. Then, the projection of a vector $\mathbf{z}\in {\mathbb{C}}^N$ onto the affine space $\Gamma =\{\mathbf{u}\in {\mathbb{C}}^N|L\mathbf{u}=\mathbf{b}\}$ can be expressed as

$${\mathrm{proj}}_{\Gamma }\left(\mathbf{z}\right)=\mathbf{z}+L^{+}\left(\mathbf{b}-L\mathbf{z}\right).$$

(6)

Proof. 

The proof is based upon the observation that the affine subspace $\Gamma$ is parallel to the subspace $\{\mathbf{u}\in {\mathbb{C}}^N|L\mathbf{u}=\mathbf{0}\}$, which is the kernel of the operator L. We thus use the projection of $\mathbf{z}$ onto the kernel of L, which is $(\mathbf{z}-L^{+}L\mathbf{z})$. The demanded projection onto $\Gamma$ is then obtained by the shift between $\Gamma$ and the kernel of L, which is $L^{+}\mathbf{b}$, leading to

$${\mathrm{proj}}_{\Gamma }\left(\mathbf{z}\right)=(\mathbf{z}-L^{+}L\mathbf{z})+L^{+}\mathbf{b},$$

(7)

which equals Equation (6) after rearranging. □

Under the assumption of Lemma 1, the relation in Equation (6) is a special case of Equation (4) in view of $L^{+}={\nu }^{-1}{L}^{\top }$ with f set to the indicator function of set $\left\{\mathbf{b}\right\}$ consisting of this single point. Then, Equation (6) belongs to the family of proximal operators. Without these assumptions, we get a different type of projection. Relaxing the assumptions on $L{L}^{*}$ (where $L^{*}$ denotes the Hermitian transpose or adjoint operator of L) to be an arbitrary diagonal matrix allows us to replace the one-point set by a more general box-type set as shown in the next subsection.

2.3. The New Relation of Projections

First, two concepts known from matrix algebra are adopted. We do this to emphasize that our result could be generalized to an infinite-dimensional (both countable and uncountable) setting, although below we reside in finite-dimensional spaces, where a linear operator uniquely coincides with a matrix. Expansions in separable infinite-dimensional spaces should be interpreted “in norm”, e.g., the notation $\mathbf{x}={\sum }_k{x}_k{\mathbf{e}}_k$ means ${lim}_{n\to \infty }\parallel \mathbf{x}-{\sum }_k^n{x}_k{\mathbf{e}}_k\parallel \to 0$.

Definition 1.

A linear operator D defined on a Hilbert space X is calleddiagonal with respect to a basis ${\left\{{\mathbf{e}}_k\right\}}_{k\in \mathbb{N}}$ if it holds

$$D{\mathbf{e}}_k={\lambda }_k{\mathbf{e}}_k\mathrm{for}\mathrm{all}k\in \mathbb{N},$$

(8)

where ${\lambda }_k$ are complex scalars. In the finite-dimensional complex space, by calling an operator simply diagonal we mean diagonal with respect to the canonical (standard) basis, i.e., to the set $[1,0,\dots ,0],\dots ,[0,\dots ,0,1]$.

The same way as a diagonal matrix does, a bounded diagonal operator performs entry-wise, since $X\ni \mathbf{x}={\sum }_k{x}_k{\mathbf{e}}_k$ and thus $D\mathbf{x}={\sum }_k{x}_{k}D{\mathbf{e}}_k={\sum }_k\left({\lambda }_k{x}_k\right){\mathbf{e}}_k$.

Lemma 4 (the new lemma).

Let $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$ be a surjective linear operator with $M\le N$ and let $L{L}^{*}$ be a diagonal operator. Assume multidimensional interval bounds ${\mathbf{b}}_1,{\mathbf{b}}_2\in {\tilde{\mathbb{R}}}^M$ such that ${\mathbf{b}}_1\le {\mathbf{b}}_2$, with the inequality being interpreted element-wise. Then, the projection of a vector $\mathbf{z}\in {\mathbb{C}}^N$,

$$\begin{array}{cc}\hfill {\mathrm{proj}}_{\Gamma }\left(\mathbf{z}\right)& =\underset{\mathbf{u}\in \Gamma }{arg~min}{\parallel \mathbf{z}-\mathbf{u}\parallel }_2,\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill \mathit{where}\Gamma & =\{\mathbf{u}\in {\mathbb{C}}^N|\Re \left(L\mathbf{u}\right)\in [{\mathbf{b}}_1,{\mathbf{b}}_2],\Im \left(L\mathbf{u}\right)=0\},\hfill \end{array}$$

(9b)

can be evaluated as

$${\mathrm{proj}}_{\Gamma }\left(\mathbf{z}\right)=\mathbf{z}+L^{+}\left({\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(L\mathbf{z}\right)-L\mathbf{z}\right),$$

(10)

where

$${\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(\mathbf{y}\right)=min\left(max({\mathbf{b}}_1,\Re \left(\mathbf{y}\right)),{\mathbf{b}}_2\right),\mathbf{y}\in {\mathbb{C}}^M.$$

(11)

Here, Equation (11) is the projection of the complex vector $\mathbf{y}$ onto a real multidimensional interval $[{\mathbf{b}}_1,{\mathbf{b}}_2]\in {\tilde{\mathbb{R}}}^M$ with min and max functions returning pairwise extremes entry-by-entry.

The proposed lemma says that a projection (following a linear transform L) onto the box-type set $\Gamma$ can be made simpler and faster using projection onto the interval $[{\mathbf{b}}_1,{\mathbf{b}}_2]$, which does not involve L. In Equation (10), the application of $L^{+}$ reduces to entrywise multiplication by the inverse diagonal of $L{L}^{*}$ followed by $L^{*}$ (see Lemma 2). Therefore, the cost of $L^{+}$ will typically be in the order of the cost of $L^{*}$.

Before we prove the main Lemma 4, we present several auxiliary lemmas needed below.

Lemma 5.

Let $L:{\mathbb{C}}^N\to {\mathbb{C}}^M,M\le N$ be a surjective operator, and $L{L}^{*}$ diagonal with respect to a certain basis. Then, ${\left(L^{+}\right)}^{*}{L}^{+}$ is diagonal with respect to the same basis.

Proof. 

The operator L is surjective and therefore its pseudoinverse can be computed as in Lemma 2. Then,

$${\left(L^{+}\right)}^{*}{L}^{+}={\left({\left(L{L}^{*}\right)}^{-1}\right)}^{*}L{L}^{*}{\left(L{L}^{*}\right)}^{-1}={\left({\left(L{L}^{*}\right)}^{-1}\right)}^{*}={\left(L{L}^{*}\right)}^{-1},$$

(12)

which is the inverse of the diagonal operator $L{L}^{*}$, and thus it is diagonal as well. □

Lemma 6.

Given any linear operator $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$ and vector $\mathbf{y}\in {\mathbb{C}}^N$, the vector $L^{+}L\mathbf{y}\in {\mathbb{C}}^N$ is the vector of minimal norm that is mapped by L to $L\mathbf{y}$.

Proof. 

Consider the linear system $L\mathbf{x}=L\mathbf{y}$ with $\mathbf{y}$ fixed. Using the properties of the pseudoinverse (Lemma 2), the minimum-norm solution to such a system is found by pseudoinverting the right hand side, i.e., by $L^{+}\left(L\mathbf{y}\right)$. □

Lemma 7.

Let $W:{\mathbb{C}}^M\to {\mathbb{C}}^N$ be such an operator that $W^{*}W$ is diagonal with respect to a given orthonormal basis $\left\{{\mathbf{e}}_n\right\}$ and let $\mathbf{x},\mathbf{y}\in {\mathbb{C}}^M$. Then, from $|x_m|\le |y_m|\forall m$, where $x_m,y_m$ are the vector coordinates with respect to the same orthonormal basis $\left\{{\mathbf{e}}_n\right\}$, it follows that ${\parallel W\mathbf{x}\parallel }_2\le {\parallel W\mathbf{y}\parallel }_2$.

Proof. 

Since $W^{*}W$ is a diagonal operator with respect to $\left\{{\mathbf{e}}_n\right\}$, according to Definition 1, we have $W^{*}W{\mathbf{e}}_m={\lambda }_m{\mathbf{e}}_m$ for a set of scalars $\left\{{\lambda }_m\right\}$, which in this case are real and nonnegative thanks to the positive semidefiniteness of $W^{*}W$. Now, ${\parallel W\mathbf{x}\parallel }_2^2=\langle W\mathbf{x},W\mathbf{x}\rangle =\langle W^{*}W\mathbf{x},\mathbf{x}\rangle$ using the basic property of the adjoint operator. Thus, for $\mathbf{x}={\sum }_m{x}_m{\mathbf{e}}_m$, it holds that

$$\begin{array}{c}\hfill {\parallel W\mathbf{x}\parallel }_2^2=〈W^{*}W\left(\sum _m{x}_m{\mathbf{e}}_m\right),\sum _m{x}_m{\mathbf{e}}_m〉=\langle \sum _m{x}_m{W}^{*}W{\mathbf{e}}_m,\sum _m{x}_m{\mathbf{e}}_m\rangle \\ \hfill =\langle \sum _m{x}_m{\lambda }_m{\mathbf{e}}_m,\sum _m{x}_m{\mathbf{e}}_m\rangle =\sum _m{\lambda }_m{\left|x_m\right|}^2,\end{array}$$

(13)

where in the last equality we have used the linearity of the inner product and the orthonormality of $\left\{{\mathbf{e}}_m\right\}$. If $|x_m|\le |y_m|$, then also their squares satisfy $|x_m{|}^2\le {\left|y_m\right|}^2$ for each m and therefore ${\parallel W\mathbf{x}\parallel }_2^2={\sum }_m{\lambda }_m|x_m{|}^2\le {\sum }_m{\lambda }_m|y_m{|}^2={\parallel W\mathbf{y}\parallel }_2^2$, since all ${\lambda }_m$ are nonnegative. □

Observe that $W^{*}W$ is symmetric and nonnegatively definite and thus always diagonalizable in the eigenvector orthonormal basis using the spectral theorem. Unfortunately, in applications of Lemma 4, the problem statement is typically not invariant with respect to unitary transforms (see the declipping problem in Section 4). Thus, we have to keep to the standard basis when expanding the projected vector in most cases.

Proof of Lemma 4.

First, we show that Equation (11) is indeed a projector onto $[{\mathbf{b}}_1,{\mathbf{b}}_2]$. Since the set $[{\mathbf{b}}_1,{\mathbf{b}}_2]$ is a multidimensional box (possibly open towards plus or minus infinity), the projection ${\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(\mathbf{y}\right)$ can be evaluated entry-by-entry. Fix m and consider the mth entry ${\mathbf{y}}_m=:y,{\left({\mathbf{b}}_1\right)}_m=:b_1,{\left({\mathbf{b}}_2\right)}_m=:b_2$. In accordance with the definition of projection,

$$\begin{array}{c}\hfill {\mathrm{proj}}_{[b_1,b_2]}\left(y\right)=\underset{b_1\le x\le b_2}{arg~min}\left|x-y\right|=\underset{b_1\le x\le b_2}{arg~min}\sqrt{{(x-\Re \left(y\right))}^2+\Im {\left(y\right)}^2}\\ \hfill =\underset{b_1\le x\le b_2}{arg~min}\left({(x-\Re \left(y\right))}^2+\Im {\left(y\right)}^2\right)=\underset{b_1\le x\le b_2}{arg~min}|x-\Re \left(y\right)|,\end{array}$$

(14)

i.e., the original task reduces to projecting $\Re \left(y\right)$ onto $[b_1,b_2]$. We prove Equation (11) by evaluating the only three possibilities (for a single entry but applied to all of them):

$${\mathrm{proj}}_{[b_1,b_2]}\left(y\right)=min(max(b_1,\Re \left(y\right)),b_2)=\left\{\begin{array}{cc}{b}_1\hfill & \mathrm{for}\Re \left(y\right)<b_1,\hfill \\ b_2\hfill & \mathrm{for}\Re \left(y\right)>b_2,\hfill \\ y\hfill & \mathrm{for}{b}_1\le \Re \left(y\right)\le b_2.\hfill \end{array}\right.$$

(15)

Now, return back to Equation (10). Suppose that $\mathbf{z}\in \Gamma$; then, $\Re \left(L\mathbf{z}\right)\in [{\mathbf{b}}_1,{\mathbf{b}}_2]$; from Equation (11), we have ${\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(L\mathbf{z}\right)=L\mathbf{z}$; and, from Equation (10), we get ${\mathrm{proj}}_{\Gamma }\left(\mathbf{z}\right)=\mathbf{z}$, which is correct. In the complementary case, $\mathbf{z}\notin \Gamma$ and due to L being surjective, there exists at least one $\mathbf{v}\in {\mathbb{C}}^N$ for which

$$L\mathbf{v}=L\mathbf{z}-{\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(L\mathbf{z}\right)=:\mathbf{w}.$$

(16)

Trivially, any such vector $\mathbf{v}$ satisfies $L(\mathbf{z}-\mathbf{v})={\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(L\mathbf{z}\right)\in [{\mathbf{b}}_1,{\mathbf{b}}_2]$, i.e., $\mathbf{z}-\mathbf{v}\in \Gamma$. In line with Equation (9a), we are looking for the shortest $\mathbf{v}$, which is found using the pseudoinverse:

$${\mathbf{v}}^{+}:=L^{+}\mathbf{w}.$$

(17)

Let us pick another vector ${\mathbf{v}}^{\prime }\in {\mathbb{C}}^N$ such that $(\mathbf{z}-{\mathbf{v}}^{\prime })\in \Gamma$, i.e., $L(\mathbf{z}-{\mathbf{v}}^{\prime })\in [{\mathbf{b}}_1,{\mathbf{b}}_2]$. We show that $\parallel {\mathbf{v}}^{\prime }{\parallel }_2\ge {\parallel {\mathbf{v}}^{+}\parallel }_2$. The projection in Equation (11) is performed entrywise and thus it must hold that

$$|w_m|\le |{\left(L{\mathbf{v}}^{\prime }\right)}_m|\mathrm{for}m=1,\dots ,M.$$

(18)

Figure 1 shows an illustration.

From Lemma 5, it follows that ${\left(L^{+}\right)}^{*}{L}^{+}$ is diagonal. Hence, we have by Lemma 7

$$\parallel {\mathbf{v}}^{+}{\parallel }_2=\parallel L^{+}{\mathbf{w}\parallel }_2\le {\parallel L^{+}\left(L{\mathbf{v}}^{\prime }\right)\parallel }_2.$$

(19)

Finally, Lemma 6 states that $\parallel L^{+}\left(L{\mathbf{v}}^{\prime }\right){\parallel }_2\le {\parallel {\mathbf{v}}^{\prime }\parallel }_2$, which together with Equation (19) concludes the proof. □

3. Discussion on the New Result

Our Lemma 4 allows us to cover a broader scope of operators than Lemma 1 does. Specifically, L can be an arbitrary complex operator, and, moreover, $L{L}^{*}$ can be a diagonal operator, not just a multiple of the identity. On the other hand, the generalized relation only holds for a projector onto a box-type set and not for a generic proximal operator.

Example 1 shows that projecting onto an ellipsoid could be performed quickly in the special case of $L{L}^{*}=\nu \mathit{Id}$. In the case of the diagonal $L{L}^{*}$, this is no longer possible.

Remark 1.

Lemma 4 clearly holds when L is the synthesis operator of a tight frame (in such a case, $L{L}^{*}=\nu \mathit{Id}$ is a multiple of the identity). The property is still valid when such a synthesis operator is weighted column-by-column by a diagonal operator D, which is real and positive. This concept actually corresponds to the so-called weighted frames [23,24,25]. Substituting $DL$ for L still meets the assumptions of the theorem, since $\left(DL\right){\left(DL\right)}^{*}=DL{L}^{*}{D}^{*}=\nu D{D}^{*}$, which is again diagonal. Such a particular setup is referred to as the painless non-orthogonal expansion in the literature [26].

Remark 2.

On the contrary, consider that L is an analysis operator, $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$, $M>N$ (i.e., $L^{*}$ is the synthesis operator). First, $L{L}^{*}$ cannot be diagonal by the same argumentation as on page 3 (while actually $L^{*}L$ can be diagonal). Lemma 4 cannot be formulated for such L, even if such a condition in the assumption is omitted; the problem appears early in the proof, in Equation (16). This linear system does not need to have any solution $\mathbf{v}$ for fixed $\mathbf{z}$, and therefore by using Equation (11) we do not need to obtain a vector in Γ.

Remark 3.

It is clear that, if we choose L surjective, but not obeying that $L{L}^{*}$ is diagonal, the proof is ruined in Equation (19), since Lemmas 7 and 5 cannot be used. In such a situation, Equation (10) would provide a feasible solution to the problem in Equation (9), which would not guarantee to be the optimal one.

Remark 4.

Consider the common discrete Fourier transform (DFT) (see [17,27] or [28], Section 3.6). The related synthesis and analysis operators, L and $L^{*}$, respectively, are complex and unitary. Without details, which come in Section 4.2, recall that they have such a special complex-conjugate structure that, when $\mathbf{x}$ is a real signal, $L{L}^{*}\mathbf{x}=\mathbf{x}$ is real again, while $L^{*}\mathbf{x}$ bears the same complex-conjugate structure. In such a case, when the coefficient vector $\mathbf{z}$ is fed into Equation (10), $L\mathbf{z}$ is real and so is its projection ${\mathrm{proj}}_{[{\mathbf{b}}_1,{\mathbf{b}}_2]}\left(L\mathbf{z}\right)$. Thus, $L^{+}$ is applied to a real vector. Since $L{L}^{*}=\mathit{Id}$ for the DFT, it follows that the pseudoinverse $L^{+}=L^{*}=L^{-1}$ bears the structure of $L^{*}$. As a consequence, the resulting projection is also complex-conjugate. Therefore, in theory, taking the real part in Equation (11) is redundant for the DFT. It is shown below that time–frequency operators used in the experiment have this feature as well.

4. Experiment

The so-called clipping is an artifact occurring in signal acquisition or processing, when the dynamic range of a signal is larger than the representation range of a device which is used to record or process the signal. In audio processing, this typically happens when an analog signal is captured with an analog-to-digital converter whose pre-amplifier is set such that the largest signal values are beyond the maximum digital representation level. Clipping is an annoying effect, and it introduces plenty of unwanted higher harmonic components to the signal [29].

The process of restoring as much audio information as possible, based on the degraded signal, is usually called declipping. Such a restoration task is clearly ill-conditioned, since there are infinitely many possible solutions to the problem. Declipping methods thus have to rely on additional information that characterizes the signal.

Besides other types of signal modeling [30,31,32,33,34], sparse priors have recently become very popular. Such a type of regularization assumes that the signal can be well approximated by only a few synthesis coefficients. As the transforms, the cosine (DCT and MDCT) [35,36] or the short-time Fourier transform (STFT) [37,38] are typically used. Since finding truly sparse solutions is NP-hard, various strategies are used to approximate the solution. Two major groups of algorithms can be identified: the greedy algorithms [39,40,41,42,43] and the convex-relaxation algorithms [37,38,44], or algorithms combining both approaches [36].

In our experiment, a convex formulation of the declipping problem was utilized with (synthesis) sparsity in the STFT (also termed Gabor) domain. It is shown that, when the Gabor frame fulfills the requirements of Lemma 4, as a consequence the projection onto the time-domain constraints can be done in a single step, leading to the use of a simple algorithm and speeding up the convergence by a factor of 2–3, compared to the situation when the projection is not used. We show the benefit of our novel projection on a simple declipping model, being aware that more sophisticated models cited above lead to a better quality of reconstruction.

4.1. Problem Formulation

It is assumed that a discrete-time signal ${\mathbf{y}}_{\mathrm{c}}\in {\mathbb{R}}^M$ is given that has been hard-clipped, i.e., any original time-domain sample above the clipping level ${\theta }_{\mathrm{H}}$ has been substituted by ${\theta }_{\mathrm{H}}$ and analogously for the samples below the level ${\theta }_{\mathrm{L}}$. The rest of the samples in ${\mathbf{y}}_{\mathrm{c}}$ have not been altered by clipping, and they are called reliable.

It is furthermore assumed that the unknown original signal $\mathbf{y}$ is sparse or compressible in the STFT domain, formally written as $\mathbf{y}\approx G\mathbf{c}$ with ${\parallel \mathbf{c}\parallel }_0\ll N$, where G is the linear Gabor synthesis operator [17], and where the pseudo-norm ${\parallel \mathbf{c}\parallel }_0$ counts the non-zero elements of $\mathbf{c}$. Characterizing the audio signal with such an operator G represents a suitable choice, since it is close to the kind of time-frequency analysis that takes place in the process of human perception [19,21]. The assumption of sparsity in the STFT domain is most realistic for solo instruments producing tones, for example the violin, double bass, piano, guitar, etc. [45,46] For instruments producing rich spectra as the drums or even for musical bands or orchestras, the number of components in the STFT domain increases and the recovery results get worse. This goes hand in hand with the reasoning that clipping enriches the STFT representation, while the sparsity-based declipping attempts to revert this process and get back the sparse representation.

In this paper, we utilize the ${\ell }_1$-norm as a convex surrogate of the sparsity pseudo-norm [47] and we formulate the optimization problem as

$$\underset{\mathbf{c}}{arg~min}{\parallel \mathbf{c}\parallel }_1\mathrm{subject}\mathrm{to}\left\{\begin{array}{cc}{M}_{\mathrm{R}}G\mathbf{c}=M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}}\hfill & (\mathrm{reliable}\mathrm{samples})\hfill \\ M_{\mathrm{H}}G\mathbf{c}\ge {\theta }_{\mathrm{H}}\hfill & (\mathrm{samples}\mathrm{clipped}\mathrm{from}\mathrm{above})\hfill \\ M_{\mathrm{L}}G\mathbf{c}\le {\theta }_{\mathrm{L}}\hfill & (\mathrm{samples}\mathrm{clipped}\mathrm{from}\mathrm{below}),\hfill \end{array}\right.$$

(20)

where $M_{\mathrm{R}},M_{\mathrm{H}},M_{\mathrm{L}}$ are simple projection operators (masks) selecting only samples related to the reliable and the clipped parts, respectively. To be precise, if we denote $R,H,L$ the disjoint sets coinciding with the three conditions in Equation (20), i.e.,

$$R=\left\{\mathbf{c}\right|M_{\mathrm{R}}G\mathbf{c}=M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}}\},H=\left\{\mathbf{c}\right|M_{\mathrm{H}}G\mathbf{c}\ge {\theta }_{\mathrm{H}}\},L=\left\{\mathbf{c}\right|M_{\mathrm{L}}G\mathbf{c}\le {\theta }_{\mathrm{L}}\},$$

(21)

then $M_{\mathrm{R}},M_{\mathrm{H}},M_{\mathrm{L}}$ select from the complete set of indexes $\{1,\dots ,M\}=R\cup H\cup L$. The formulation in Equation (20) thus looks for coefficients that are approximately sparse and at the same time generate a signal that is consistent with the time-domain constraints.

4.2. The Gabor Operators

As a prelude to Gabor operators, return to the ordinary DFT. The DFT analysis, $F^{*}:{\mathbb{C}}^M\to {\mathbb{C}}^Q$, computes scalar products of the signal with complex exponentials

$$\frac{1}{\sqrt{Q}}{\left[exp\left(\mathrm{i}2\pi 0\frac{q}{Q}\right),exp\left(\mathrm{i}2\pi 1\frac{q}{Q}\right),\dots ,exp\left(\mathrm{i}2\pi (M-1)\frac{q}{Q}\right)\right]}^{\top }$$

(22)

for each frequency $q=0,\dots ,Q-1$. Following Remark 4, we now recall that the DFT operator takes the real signal into the complex domain of dimension Q such that the image has a complex-conjugated structure. To make this statement precise, $F^{*}:{\mathbb{R}}^M\to {\mathbb{K}}^Q$, where we define

$${\mathbb{K}}^Q=\left\{\mathbf{c}\in {\mathbb{C}}^Q|c_0,c_{\frac{Q}{2}}\in \mathbb{R};c_q=\overline{c_{Q-q}},q=1,\dots ,\frac{Q}{2}-1\right\}.$$

(23)

Such a structure follows immediately from Equation (22). In the common DFT case, $Q=M$, which implies that Equation (22) is an orthonormal basis of ${\mathbb{C}}^M$, and thus $F^{*}=F^{-1}$. Nevertheless, the relation between coefficients holds true also in the case $Q>M$ (with Q even) in Equation (22), which is referred to as the overcompleteness in frequency and sometimes as the redundant DFT. In such a situation, the complex exponentials form a Parseval tight frame of ${\mathbb{C}}^M$ [17]. Both cases can be computed using the FFT.

The Gabor transform is a time-frequency transform with an extra structure of the coefficients in the time direction compared to the DFT. The Gabor analysis $G^{*}$ of a signal consists of windowing it and then taking scalar products of such windowed signal with complex harmonics in Equation (22); effectively, $G^{*}$ amounts in computing the FFTs of the windowed signals. The window moves along the signal by means of translation by a constant amount of samples (see below). Generally, $G^{*}:{\mathbb{C}}^M\to {\mathbb{C}}^N$, $N>M$ and the other way round for the synthesis G.

Conditions are known which guarantee that the frame operator $G{G}^{*}$ is invertible. The invertibility is a crucial property since its violation would mean that the synthesis operator G is not surjective. Under even stronger assumptions, the related synthesis operator G moreover bears the same time-frequency structure as $G^{*}$ does; this fact is handy both for the theory and practice, since in such a case, for example, the synthesis can be expressed using windows moving along time direction, as is the case of the analysis. For more details on this wide topic, see, for example, ([18], [Ch. 9]) or ([17], [Ch. 9 and 10]) or [19]. From the computational point of view, systems for which $G{G}^{*}=\nu {\mathit{Id}}_M$ are preferably used in applications (see, e.g., [21,41,48,49,50,51,52]), however let us note that the assumptions of our lemma, in addition, cover an important class of the so-called painless Gabor transforms [26], for which the operator $G{G}^{*}$ is diagonal and the diagonal is not necessarily constant. Gabor Parseval tight frames thus fall into our setup as a special case.

For the Gabor transform, parameters have to be specified such as the window shape and length, time-shift a of the window, and the number of frequency channels Q [19]. The number of time samples M must be a multiple of a [51]; then, the analysis operator is $G^{*}:{\mathbb{R}}^M\to {\mathbb{C}}^{MQ/a}$. Recall now that the DFT analysis is an operator ${\mathbb{R}}^M\to {\mathbb{K}}^Q$; thanks to the just described windowing, it is clear that we can even understand the operator as $G^{*}:{\mathbb{R}}^M\to {\mathbb{K}}^Q\times \dots \times {\mathbb{K}}^Q$, where each of the vectors from ${\mathbb{K}}^Q$ originates in one of the $M/a$ possible window positions. The synthesis $G:{\mathbb{K}}^Q\times \dots \times {\mathbb{K}}^Q\to {\mathbb{R}}^M$ produces a real signal again. It is made clear in Section 4.7 that computations in the algorithms do not tackle the coefficients in such a way that they would fall out of the space ${\mathbb{K}}^Q\times \dots \times {\mathbb{K}}^Q$, making it possible to omit the real-part operator in Equation (11).

4.3. Problem Solution

We use the so-called proximal splitting algorithms to solve Equation (20). Proximal splitting algorithms are introduced in Section 2. First, it is convenient to rewrite Equation (20) to an unconstrained form using indicator functions (defined in Section 2) as

$$\underset{\mathbf{c}\in {\mathbb{C}}^N}{arg~min}{\parallel \mathbf{c}\parallel }_1+{\iota }_R\left(\mathbf{c}\right)+{\iota }_H\left(\mathbf{c}\right)+{\iota }_L\left(\mathbf{c}\right)$$

(24)

where the sets of indexes are defined in Equation (21). Once optimal coefficients are found, signal recovery is simply obtained by their synthesis through G. Equation (24) presents a sum of convex functions, considering the fact that the sets $R,H,L$ are convex. Therefore, Equation (24) is suitable to be solved by a proximal algorithm.

Recall that proximal algorithms rely on evaluating the so-called proximal operators linked to the functions present in the sum. It is known that the proximal operator of ${\lambda \parallel ·\parallel }_1$ is the soft thresholding with threshold $\lambda >0$ [53,54,55] defined such that ${\mathrm{soft}}_{\lambda }\left(\mathbf{z}\right)$ takes the vector argument and performs elementwise mapping

$$z_i\mapsto \mathrm{sgn}\left(z_i\right)·max\left(\right|z_i|-\lambda ,0).$$

(25)

The proximal operator of ${\iota }_K$ is the projection onto the convex set K [1,2], which is discussed in Section 2.2. Performing soft thresholding is a simple operation, but, regarding projections onto the sets $R,H,L$, it becomes more complicated. Explicit formulas for projecting onto these sets are known only for G being an unitary operator or a synthesis operator of a tight frame ($G{G}^{*}=\nu {\mathit{Id}}_M$) defined in real domain (see Lemma 1). Our G, however, maps ${\mathbb{C}}^N\to {\mathbb{R}}^M$ and moreover there are practical Gabor systems for which $G{G}^{*}$ is only diagonal [19,26].

In this experiment, we compared two approaches to overcome this complication. One is to adapt the general Condat algorithm (CA), which allows us to simplify the projections at the cost of an additional application of G and $G^{*}$ in each iteration. The second approach utilizes the novel projection lemma and shows its superiority over the first choice.

4.4. Condat Algorithm

The problem in Equation (24) is equivalent to

$$\underset{\mathbf{c}}{arg~min}{\parallel \mathbf{c}\parallel }_1+{\iota }_R\left(\mathbf{c}\right)+{\iota }_{H^{\prime }}\left(G\mathbf{c}\right)+{\iota }_{L^{\prime }}\left(G\mathbf{c}\right)$$

(26)

where $H^{\prime }=\left\{\mathbf{z}\right|M_{\mathrm{H}}\mathbf{z}\ge {\theta }_{\mathrm{H}}\}$ and $L^{\prime }=\left\{\mathbf{z}\right|M_{\mathrm{L}}\mathbf{z}\le {\theta }_{\mathrm{L}}\}$. Note that projection onto such sets is trivial compared to projections onto H and L, since $H^{\prime }$ and $L^{\prime }$ are time-domain sets for which the closest vector is found by simple elementwise comparison of the vector with clipping thresholds ${\theta }_{\mathrm{H}}$ and ${\theta }_{\mathrm{L}}$, respectively. Regarding the reliable part of the signal, we could also resort to replacing ${\iota }_R\left(\mathbf{c}\right)$ by ${\iota }_{R^{\prime }}\left(G\mathbf{c}\right)$, where projecting onto $R^{\prime }=\left\{\mathbf{z}\right|M_{\mathrm{R}}\mathbf{z}=M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}}\}$ amounts to simple substitution of $\mathbf{z}$ with ${\mathbf{y}}_{\mathrm{c}}$ at the reliable positions. Note, however, that R is an affine set and there exists an explicit formula for the respective projection (see Lemma 3). In our setup, this lemma renders the following corollary.

Corollary 1.

Let G be a linear synthesis operator, $G:{\mathbb{C}}^N\to {\mathbb{C}}^M$ such that $G{G}^{*}$ is diagonal. Let $M_{\mathrm{R}}$ be the “reliable” mask operator, as discussed in Section 4.1. Then, for the projection onto R, we have

$${\mathrm{proj}}_R\left(\mathbf{z}\right)=\mathbf{z}+G^{*}{M}_{\mathrm{R}}^{\top }{\left(M_{\mathrm{R}}G{G}^{*}{M}_{\mathrm{R}}^{\top }\right)}^{-1}(M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}}-M_{\mathrm{R}}G\mathbf{z}).$$

(27)

Proof. 

We use Lemma 3 in the particular case of the diagonal product. The roles of L and $\mathbf{b}$ are played by $M_{\mathrm{R}}G$ and $M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}}$, respectively. The only step to comment on is the determination of the pseudoinverse operator. Recall that G is assumed to be surjective; therefore, $M_{\mathrm{R}}G$ spans the subspace onto which it projects, since $M_{\mathrm{R}}$ is applied elementwise. Therefore, we can use the formula mentioned in Lemma 2 and obtain

$$\begin{array}{cc}\hfill {\left(M_{\mathrm{R}}G\right)}^{+}& ={\left(M_{\mathrm{R}}G\right)}^{*}{\left(M_{\mathrm{R}}G{G}^{*}{M}_{\mathrm{R}}^{*}\right)}^{-1}=G^{*}{M}_{\mathrm{R}}^{\top }{\left(M_{\mathrm{R}}G{G}^{*}{M}_{\mathrm{R}}^{\top }\right)}^{-1}.\hfill \end{array}$$

(28)

Equation (27) follows immediately. □

Equation (27) employs application of one synthesis and one analysis. Note that $\left(M_{\mathrm{R}}G{G}^{*}{M}_{\mathrm{R}}^{\top }\right)$ is a square diagonal matrix whose size equals the number of reliable samples, and, as such, its inversion is simple and can be precomputed.

The complete Condat algorithm for declipping is described in Algorithm 1. The sufficient condition on the step sizes $\tau$ and $\sigma$ ensuring convergence of the algorithm involves norms of the linear operators $L_j$ of the function to minimize: $\tau \sigma \parallel {\sum }_j{L}_j^{*}{L}_j\parallel \le 1$ (see [7,11]). In our case, checking the problem in Equation (26) reveals that it involves an identity operator and twice the synthesis G, and it can be shown that our particular convergence criterion becomes $\tau \sigma \parallel \mathit{Id}+G^{*}G+G^{*}G\parallel \le 1$, where $\parallel ·\parallel$ is the operator norm. With the help of the inequality $\parallel \mathit{Id}+G^{*}G+G^{*}G\parallel \le \parallel \mathit{Id}\parallel +\parallel G^{*}G\parallel +\parallel G^{*}G\parallel$, we find a (suboptimal) bound for the step sizes. Since $G{G}^{*}$ is diagonal and it holds $\parallel L^{*}L\parallel =\parallel L{L}^{*}\parallel$ for any operator, we have $\parallel G^{*}G\parallel =max\mathrm{diag}|G{G}^{*}|=:\mu$. Thus, $\tau \sigma \le 1/(1+2\mu )$ ensures convergence; for Parseval frames, it simplifies to $\tau \sigma \le 1/3$. However, the choice of particular $\sigma$ and $\tau$ influences the convergence speed. Regarding $\rho$, which plays the role of an extrapolation parameter, setting it between 1 and 2 usually provides the fastest convergence in practice [11].

Algorithm 1: Condat algorithm (CA) adapted to solving Equation (26).

In practice, we terminate the main loop when a proper convergence criterion is fulfilled (discussed below). Usually, when the algorithm stops, the solution might happen to lie slightly out of the feasible set. Therefore, we normally perform a final projection of the current solution onto the feasible set. (This also applies to the second algorithm below.)

The advantage of the Condat algorithm is its universality; it would work for audio declipping with arbitrary (meaningful) linear operator G. In the following, however, we exploit the special properties of G, which are assumed.

4.5. Douglas–Rachford Algorithm

Now, we approach the problem in Equation (24) directly. Notice that it could be rewritten as a sum of two convex functions

$$\underset{\mathbf{c}}{arg~min}{\parallel \mathbf{c}\parallel }_1+{\iota }_K\left(\mathbf{c}\right)$$

(29)

with $K=\left\{\mathbf{z}\right|M_{\mathrm{R}}G\mathbf{z}=M_{\mathrm{R}}{\mathbf{y}}_{\mathrm{c}},M_{\mathrm{H}}G\mathbf{z}\ge {\theta }_{\mathrm{H}},M_{\mathrm{L}}G\mathbf{z}\le {\theta }_{\mathrm{L}}\}$ gathering all the declipping constraints into a single set. Due to the form of two convex summands, the simple Douglas–Rachford (DR) algorithm [2,13] can be applied to solve Equation (29), if the projector onto K is available—to develop such a projector, we use our new lemma.

First, define the “lower” and “upper” bounding vectors ${\mathbf{b}}_{\mathrm{L}},{\mathbf{b}}_{\mathrm{H}}\in {\tilde{\mathbb{R}}}^M$ such that for $m=1,\dots ,M$,

$${\left({\mathbf{b}}_{\mathrm{L}}\right)}_m=\left\{\begin{array}{cc}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m\hfill & \mathrm{for}{\theta }_{\mathrm{L}}<{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m<{\theta }_{\mathrm{H}},\hfill \\ {\theta }_{\mathrm{H}}\hfill & \mathrm{for}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m={\theta }_{\mathrm{H}},\hfill \\ -\infty \hfill & \mathrm{for}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m={\theta }_{\mathrm{L}},\hfill \end{array}\right.{\left({\mathbf{b}}_{\mathrm{H}}\right)}_m=\left\{\begin{array}{cc}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m\hfill & \mathrm{for}{\theta }_{\mathrm{L}}<{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m<{\theta }_{\mathrm{H}},\hfill \\ \infty \hfill & \mathrm{for}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m={\theta }_{\mathrm{H}},\hfill \\ {\theta }_{\mathrm{L}}\hfill & \mathrm{for}{\left({\mathbf{y}}_{\mathrm{c}}\right)}_m={\theta }_{\mathrm{L}}.\hfill \end{array}\right.$$

(30)

Note that defining the multidimensional box $[{\mathbf{b}}_{\mathrm{L}},{\mathbf{b}}_{\mathrm{H}}]$ this way matches the set of feasible solutions K. Specifically, we have $K=\left\{\mathbf{c}\right|{\mathbf{b}}_{\mathrm{L}}\le G\mathbf{c}\le {\mathbf{b}}_{\mathrm{H}}\}$. Due to Lemma 4, the projection onto K is realized via

$${\mathrm{proj}}_K\left(\mathbf{c}\right)=\mathbf{c}+G^{+}\left({\mathrm{proj}}_{[{\mathbf{b}}_{\mathrm{L}},{\mathbf{b}}_{\mathrm{H}}]}\left(G\mathbf{c}\right)-G\mathbf{c}\right),$$

(31)

where

$${\mathrm{proj}}_{[{\mathbf{b}}_{\mathrm{L}},{\mathbf{b}}_{\mathrm{H}}]}\left(G\mathbf{c}\right)=min\left(max({\mathbf{b}}_{\mathrm{L}},\Re \left(G\mathbf{c}\right)),{\mathbf{b}}_{\mathrm{H}}\right).$$

(32)

The Douglas–Rachford algorithm is presented in Algorithm 2. It always converges; nevertheless, the parameter $\gamma$ is responsible for the convergence speed.

Algorithm 2: Douglas–Rachford algorithm (DR) solving Equation (29)

4.6. Comparison of the Algorithms

If converged, both presented algorithms produce the same solution when the objective function is strictly convex ([56], Part B1). In our case, Equation (24) is not strictly convex due to the ${\ell }_1$-norm, and therefore we can hardly hope for the same minimizers output from CA and DR. The objective function, i.e., the ${\ell }_1$-norm of coefficients, should coincide, however.

Regarding the computational cost, the CA (see Algorithm 1) requires, per iteration:

• Sparsifying step: one soft thresholding, which is performed elementwise, and thus it is $\mathcal{O}\left(N\right)$, and one analysis $G^{*}$, which is $\mathcal{O}(NlogN)$
• Reliable part: one synthesis G and one analysis $G^{*}$, both $\mathcal{O}(NlogN)$
• Each of the clipped parts: one synthesis, $\mathcal{O}(NlogN)$, and one elementwise projection, $\mathcal{O}\left(N\right)$.

We neglect simple addition of vectors and multiplication of a vector with a scalar, and it is assumed that, where applicable, parts of the formulas are precomputed. It is clear that the overall cost is dominated by the transforms, namely every iteration requires three applications of G and two of $G^{*}$.

On the other hand, the DR (see Algorithm 2) requires:

• Sparsifying step: one soft thresholding, which is $\mathcal{O}\left(N\right)$
• Projection onto K: one synthesis G and one pseudoinverse $G^{+}$, which is in the order of $G^{*}$, i.e., $\mathcal{O}(NlogN)$ in our particular setup; projection that is performed elementwise, $\mathcal{O}\left(N\right)$.

Again, the overall cost is dominated by the transforms, but we can see that the DR is much less demanding per iteration. Notice that this does not automatically mean that it will converge more quickly. We provide an experiment on this issue later.

The CA is a primal–dual method ([1], Ch. 5) and needs more memory than the DR does since variables in the dual space have to be kept, compared to the DR. The ratio of the memory requirements of the two algorithms depends also on the redundancy of the transform $G^{*}$, i.e., on the ratio $M/N$.

The practical advantage of the DR lies also in the fact that only one step size has to be hand-tuned instead of two in the CA.

4.7. Redundancy of the Real-Part Operator

The operator $\Re (·)$ can be omitted in Equations (11) and (32) if its argument is real. As mentioned several times above, if the coefficient vector $\mathbf{c}\in ({\mathbb{K}}^Q\times \dots \times {\mathbb{K}}^Q)\subset ({\mathbb{C}}^Q\times \dots \times {\mathbb{C}}^Q)$, then the synthesis $G\mathbf{c}$ is real. It remains to check whether there exists an operation in any of the above algorithms that would make the coefficients fall out of this subspace: The soft thresholding is harmless in this regard. In the consolidated projection in Equation (31), $G^{+}=G^{*}{\left(G{G}^{*}\right)}^{-1}$ has the same complex-conjugate structure as $G^{*}$ since $G{G}^{*}$ is a real diagonal matrix, according to the assumption. In the projection in Equation (27) onto the reliable set, it is straightforward to see that, if the input coefficients belong to ${\mathbb{K}}^Q\times \dots \times {\mathbb{K}}^Q$, then the same holds for the output, since $\left(M_{\mathrm{R}}G{G}^{*}{M}_{\mathrm{R}}^{\top }\right)$ is a real diagonal matrix.

4.8. Results

For the following experiments, five audio excerpts sampled at 16 kHz with length of approximately 5 s were chosen. Because the goal of the experiments was the proof of concept demonstrating the convenience of the proposed projector in a practical example, rather than developing a novel method for audio declipping outperforming the current state of the art, these audio excerpts are fully sufficient for this purpose. Used signals differ in tonal content and sparsity with respect to the used time–frequency representation to cover a wide spectrum of audio signals; they consists of recordings of acoustic guitar, double bass, artificial signature tone, speech and string orchestra.

As a preprocessing step, the signals were peak-normalized and then artificially clipped to multiple thresholds ${\theta }_{\mathrm{c}}\in \{0.3,0.5,0.7\}$. The clipping was considered to be symmetric, i.e., ${\theta }_{\mathrm{c}}={\theta }_{\mathrm{H}}=-{\theta }_{\mathrm{L}}$.

Regarding the used transform, the common setting for audio declipping on 16 kHz sampled signals was adopted, i.e., the Gabor transform (STFT) with 1024-sample long Hann window (corresponding to 64 ms) and a 75% overlap. In all cases, the number of frequency channels was the same as the window length, i.e., $Q=1024$, except for the experiment shown in Figure 5, where the number of frequency channels was twice as many as the window length, specifically 2048 channels. In both settings, it holds that $G{G}^{*}=\mathit{Id}$.

It has been noticed that the convergence performance of Algorithms 1 and 2 is influenced by the setting of parameters. These parameters were manually tuned such that both declipping algorithms performed well for all testing audio excerpts and all clipping thresholds. Specifically, the parameters of the CA were set to $\tau =0.5$ and $\sigma =0.666$ and the DR parameter $\gamma$ was set to 1.

The results were evaluated using $\Delta \mathrm{SDR}$, the signal-to-distortion ratio improvement, which is the difference between SDR of the restored signal and SDR of the clipped signal. Formally,

$$\Delta \mathrm{SDR}=\mathrm{SDR}(\mathbf{y},\widehat{\mathbf{y}})-\mathrm{SDR}(\mathbf{y},{\mathbf{y}}_{\mathrm{c}}),$$

(33)

where $\mathbf{y}$ denotes the original (undistorted, ground-truth) signal, $\widehat{\mathbf{y}}$ represents the restored signal and ${\mathbf{y}}_{\mathrm{c}}$ represents the clipped signal, while the SDR is the standard ratio computed as

$$\mathrm{SDR}(\mathbf{u},\mathbf{v})=10{log}_{10}\frac{{\parallel \mathbf{u}\parallel }_2^2}{{\parallel \mathbf{u}-\mathbf{v}\parallel }_2^2}\left[\mathrm{dB}\right].$$

(34)

Figure 2 and Figure 3 are designed to demonstrate how the Condat algorithm and the Douglas–Rachford algorithm act over time, therefore the development of the $\Delta \mathrm{SDR}$ and the objective function (${\ell }_1$-norm of the coefficients $\mathbf{c}$) for both algorithms are presented in these figures. Whereas Figure 2 displays only data obtained from reconstructing the “acoustic guitar” excerpt clipped at the level of 0.3, Figure 3 uses the average of all five audio signals and three tested clipping thresholds. It can be concluded from the figures that the objective function somewhat correlates with the $\Delta \mathrm{SDR}$ curve and both cases show that both algorithms converge to the solution with the same final objective function value, but the DR algorithm converges more quickly, i.e., the curves level up more quickly, reaching slightly higher $\Delta \mathrm{SDR}$ value. This is caused most likely by the projection step, where, in the CA, there are three individual projections (onto R, H and L) combined in a sum, thus it is not as accurate as in the DR case, where the projection step in Equation (31) is exploited.

The horizontal axis was converted to time in seconds in the figures, but they also offer an overview of the number of iterations, since the iterations of CA and DR do not consume the same amount of time. For this reason, apart from the $\Delta \mathrm{SDR}$ and objective function graphs, the (+ and ×) markers emphasizing every 100th iteration are also displayed in the plots. Note that, although both algorithms ran for exactly 1000 iterations, since the time axis is limited to 26 s, only 570 iterations of the Condat algorithm are shown in Figure 2, and 581 iterations in Figure 3.

The average computational times for all sound excerpts and all computed clipping thresholds with a fixed number of 1000 iterations were 44.73 s for the CA and 23.82 s for the DR.

More detailed time comparison of both algorithms can be seen in Figure 4 and Figure 5 where the time ratio (time of DR divided by the time of CA) is plotted for every sound excerpt and clipping threshold ${\theta }_{\mathrm{c}}$ combination. Here, both algorithms ran until the objective function differed by 0.1% from the state of full convergence, the assessment of which was determined to be reached after exactly 3000 iterations. With the average value of the plotted time ratios being ca. 0.53, it is possible to claim that the DR algorithm is about twice as fast as the CA.

Figure 5 presents the same time ratios as Figure 4 with only one difference—the number of frequency channels of the Gabor transform is set to 2048 instead of 1024. The average value of the time ratio is now around 0.42. It can be assumed that, the more DGT coefficients we work with (longer windows, bigger overlaps or more frequency channels), the more significant time difference between the two algorithms there would be.

The algorithms were implemented and tested in MATLAB R2017a using the LTFAT toolbox [51,57] and, supporting the idea of reproducible research, the codes are available at http://www.utko.feec.vutbr.cz/~rajmic/software/project_accel_declip.zip in supplementary. All experiments ran on a PC with Intel i7-3770 3.4 GHz, 24 GB RAM and OS Windows 10 Pro (version 1809).

4.9. Other Applications

The new lemma is useful in situations where a problem to be solved requires projecting onto the box-type feasible set $\Gamma$ as defined in Equation (), under the condition that $L:{\mathbb{C}}^N\to {\mathbb{C}}^M$ with $M\le N$ and $L{L}^{*}$ is a diagonal operator. Such a setup corresponds to the “painless” operator L as the operator that is used to synthesize the signal from its coefficients.

As the natural extension of the audio declipping problem, the image desaturation (i.e., recovering images that have been “clipped” to the maximum value due to a low dynamical range of a digital system) can benefit from the lemma as well. In a similar fashion, our projection can be used to keep the processed signal (or image) in a prescribed dynamic range, as present in the example in [11]. Dequantization is another example: each quantized signal sample has its origin in one of the non-overlapping intervals and our lemma is suitable to find signal coefficients that are consistent with the quantized observation [58]. As the last example, we name the signal inpainting problem, where a portion of samples is completely missing, while the others are considered reliable [35,50,59]. In this problem, nothing is directly known about the missing samples, thus the set $\Gamma$ is defined such that the lower and upper bounds are set to minus and plus infinity, respectively; with this definition, projection lemma can be used here as well.

5. Conclusions

A new projector onto a box-shaped convex set in the transformed domain is presented in this article. The projector forms a generalization of previous results and its use was demonstrated on a simple audio signal declipping task.