Psychoacoustically Motivated Audio Declipping Based on Weighted l1 Minimization

TL;DR

This paper introduces a psychoacoustically motivated audio declipping method that uses weighted l1 minimization, improving restoration quality by incorporating hearing thresholds and masking effects.

Contribution

It presents a novel weighting scheme based on psychoacoustic models for sparse audio declipping, enhancing performance over existing methods.

Findings

Improved SDR and PEMO-Q scores with proper weighting.

Outperforms current state-of-the-art declipping techniques.

Low computational complexity maintained.

Abstract

A novel method for audio declipping based on sparsity is presented. The method incorporates psychoacoustic information by weighting the transform coefficients in the $\ell_1$ minimization. Weighting leads to an improved quality of restoration while retaining a low complexity of the algorithm. Three possible constructions of the weights are proposed, based on the absolute threshold of hearing, the global masking threshold and on a quadratic curve. Experiments compare the restoration quality according to the signal-to-distortion ratio (SDR) and PEMO-Q objective difference grade (ODG) and indicate that with correctly chosen weights, the presented method is able to compete, or even outperform, the current state of the art.

Equations

$$y_n=\{\begin{array}{rlrl}& x_n& \text{for}& |x_n|<{\theta }_{\text{c}},\\ & {\theta }_{\text{c}}\cdot \mathrm{sgn}(x_n)& \text{for}& |x_n|\ge {\theta }_{\text{c}},\end{array}$$

$${\arg \min }_{\mathbf{c}}\Vert \mathbf{c}{\Vert }_1\text{ s.t. }\mathbf{c}\in \Gamma ,$$

$$\Gamma =\{\mathbf{c}|M_{\text{R}}D\mathbf{c}=M_{\text{R}}\mathbf{y},M_{\text{H}}D\mathbf{c}\ge {\theta }_{\text{c}},M_{\text{L}}D\mathbf{c}\le -{\theta }_{\text{c}}\},$$

$${\arg \min }_{\mathbf{c}}\Vert \mathbf{c}{\Vert }_1+{\iota }_{\Gamma }(\mathbf{c}),$$

$$\iota_{\Gamma}(\mathbf{c})=\left\{\begin{array}[]{ll}0&\text{for }\mathbf{c}\in\Gamma,\\ +\infty&\text{for }\mathbf{c}\notin\Gamma.\end{array}\right.$$

$${\arg \min }_{\mathbf{c}}\Vert \mathbf{w}\odot \mathbf{c}{\Vert }_1+{\iota }_{\Gamma }(\mathbf{c}).$$

$${\mathrm{proj}}_{\Gamma }(\mathbf{z})=\mathbf{z}-D^{\ast }\left(D\mathbf{z}-\mathrm{proj}(D\mathbf{z})\right),$$

$$\Big{(}\mathrm{proj}(\mathbf{z})\Big{)}_{n}=\left\{\begin{array}[]{ll}y_{n}&\text{for }n\in R,\\ \max(\theta_{\text{c}},z_{n})&\text{for }n\in H,\\ \min(-\theta_{\text{c}},z_{n})&\text{for }n\in L.\end{array}\right.$$

$${\mathrm{soft}}_{\gamma \mathbf{w}}(\mathbf{z})=\mathrm{sgn}(\mathbf{z})\odot \max (|\mathbf{z}|-\gamma \mathbf{w},0).$$

$$T_q(f)=3.64g^{-0.8}-6.5{\mathrm{e}}^{-0.6{\left(g-3.3\right)}^2}+10^{-3}{g}^4,$$

$${\mathbf{w}}_{\text{ATH1}}=(\mathbf{t}-\min (\mathbf{t})+1{)}^{-1},$$

$${\mathbf{w}}_{\text{ATH2}}=-\mathbf{t}+\tau ,$$

$${\mathbf{w}}_{\text{ATH3}}=2\cdot 10^{-5}\cdot 10^{(-\mathbf{t}+\tau )/20},$$

$$\text{SDR}(\mathbf{u},\mathbf{v})=10{\log }_{10}\frac{\Vert \mathbf{u}{\Vert }_2^2}{\Vert \mathbf{u}-\mathbf{v}{\Vert }_2^2}.$$

Full text

I Introduction

The so-called clipping is a non-linear form of signal distortion usually appearing when the signal exceeds its allowed dynamic range. This unwanted phenomenon is most commonly present in audio signals, but clipping can affect all kinds of signal. Clipping has a great negative effect on the perceptual quality of audio signal [1], it reduces the accuracy of automatic voice recognition [2, 3] or it may even potentially damage the audio speaker [4]. To enhance the perceived audio quality and to reduce the negative effects of clipping, it is necessary to perform the recovery of clipped samples. This process is usually termed declipping.

In the past, several declipping methods using various approaches were introduced. The first attempts were based on autoregressive modeling [5] or on the knowledge of the original signal bandwidth [6]. Statistical approaches were exploited in [7, 8] and more recently in [9]. Also, a method based on the Hankel matrix rank minimization was used in [10].

Well-performing methods for the tasks of audio restoration proved to be methods using the sparse prior; for instance methods based on Orthogonal Matching Pursuit (OMP) [11], convex optimization [12, 13], social sparsity [14], or non-convex methods based on the Alternating Directions Method of Multipliers (ADMM) [15, 16, 17, 18, 19]. Besides sparsity-based approaches, methods using non-negative matrix factorization were also adapted for audio declipping [20, 21]. To the best of our knowledge, deep learning techniques have only been used for image declipping [22], not for audio.

All the methods mentioned above (except for [13]) have not exploited any additional information in order to achieve the best possible perceived audio quality of the restored signal, which does not necessarily fully coincide with the physical signal restoration quality. In [13], the authors utilized the effect of simultaneous masking and used the MPEG psychoacoustic model to weight time-frequency coefficients during the restoration process. Such an approach discourages the introduction of distinctively audible signal components (where the masking threshold is low), which are not likely to be present in the original signal, and signal components that are less audible (masking threshold is high) are tolerated to a greater extent. Unfortunately, the algorithm presented in [13] is not fully consistent which means that it allows deviating the reconstructed samples from the reliable samples (i.e., samples that have not been clipped).

In this paper, a fully consistent declipping method based on weighted $\ell_{1}$ minimization is presented. First, the declipping problem is formulated and the minimization task is constructed in Sec. II. The problem is then numerically solved by the Douglas-Rachford algorithm [23], described in Sec. III. Three different approaches incorporating psychoacoustical information in the audio declipping are introduced in Sec. IV. Finally, Sec. V reports on the experiments that have been run and the results obtained.

II Problem formulation

In the case of the so-called hard clipping, which is the subject of interest in this paper, the signal exceeding the prescribed dynamic range $[-\theta_{\text{c}},\theta_{\text{c}}]$ is limited in amplitude such that

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{N}$ denotes the original (clean) signal and $\mathbf{y}\in\mathbb{R}^{N}$ the observed clipped signal. The limiting constant $\theta_{\text{c}}$ is referred to as the clipping threshold. The subscript $n$ in $y_{n}$ refers to the $n$-th sample of signal $\mathbf{y}$, the same as $(\mathbf{y})_{n}$ does.

According to (1), it is possible to divide the signal samples into three disjoint sets $R,H,L$ such that $R\cup H\cup L=\{1,\ldots,N\}$ and, correspondingly, to distinguish the reliable samples (not influenced by clipping), samples clipped from above to the high clipping threshold $\theta_{\text{c}}$ and samples clipped from below to the low clipping threshold ($-\theta_{\text{c}}$), respectively. To select only samples from one specific set, the respective restriction operators $M_{\text{R}}$, $M_{\text{H}}$ and $M_{\text{L}}$ are used.

In the restoration task, it is natural to desire that the recovered signal should not differ from the clipped signal $\mathbf{y}$ at the reliable positions, and at the clipped positions its samples should lie above $\theta_{\text{c}}$ or below $-\theta_{\text{c}}$. Such a task is ill-posed, since the solution is not unique. Therefore, considering some additional information about the signal is crucial. Based on the fact that most of the musical signals are approximately sparse with respect to a certain time-frequency transform, the restoration task can be formulated as finding a signal whose coefficients will be of the highest sparsity. As the true sparsity measure, the $\ell_{0}$-(pseudo)norm can be used; however, for practical reasons, it is usually approximated by various techniques. In this paper, the $\ell_{1}$-norm is used as the “closest” convex surrogate of the $\ell_{0}$-norm. Formally, the minimization problem can be written as

[TABLE]

where $\Gamma\subset\mathbb{C}^{P}$ is a (convex) set of feasible solutions defined as

[TABLE]

where $D:\mathbb{C}^{P}\rightarrow\mathbb{R}^{N}$ is the synthesis operator of a Parseval tight frame, with $P\geq N$ [24, 17, 18]. More specifically, throughout this paper, the Discrete Gabor Transform (DGT) is used in place of the time-frequency transform. This transform is also commonly known as the Short-Time Fourier Transform (STFT). The particular setting of the transform is described in Sec. V.

To find an appropriate algorithm to solve (2), it is convenient to rewrite the problem to an unconstrained form:

[TABLE]

where the hard constraint from (2) is replaced by the indicator function of the set $\Gamma$ defined as

[TABLE]

Problem (4) on its own does not support the incorporation of any additional (psychoacoustical) information. However, it is possible to replace the $\ell_{1}$-norm with the weighted $\ell_{1}$-norm by exploiting a vector of non-negative weights $\mathbf{w}\in\mathbb{R}^{P}$. The minimization task will then attain the following form:

[TABLE]

Here the $\odot$ symbol represents the elementwise product. The larger an element in $\mathbf{w}$ is, the more the respective coefficient in $\mathbf{c}$ gets pushed towards zero in the optimization.

III Algorithmic solution

The Douglas-Rachford algorithm (DRA) [23] is able to find the minimizer of a sum of two convex functions. Problem (6) is exactly an instance of such a setup, since both the summands in the objective function are convex (although not strictly convex). The algorithm is presented in Alg. 1.

The algorithm iterates over two main steps. The first is the projection on $\Gamma$, which corresponds to the proximal operator of $\iota_{\Gamma}$, and the soft thresholding $\mathrm{soft}_{\gamma\mathbf{w}}$ with the vector of thresholds $\gamma\mathbf{w}$, which coincides with the proximal operator of the weighted $\ell_{1}$-norm.

Since only the Parseval tight frames are considered, the projection step can be computed via the explicit formula [18, 25]:

[TABLE]

where $D^{*}:\mathbb{R}^{N}\rightarrow\mathbb{C}^{P}$ is the analysis operator, and the inner projection step is a simple projection in the time domain, which is computed elementwise as:

[TABLE]

The soft thresholding, as the second main step of the DR algorithm, is computed according to

[TABLE]

IV Choice of Weights

The Douglas-Rachford algorithm for audio declipping expects the weighting vector $\mathbf{w}$ at the input. The performance of the algorithm is greatly influenced by this vector. Effectively, weights affect the signal coefficients during the soft thresholding step; recalling (9), the higher weight belongs to a particular coefficient, the more shrunk the coefficient will become. Therefore, it is possible to encourage selecting the important coefficients and, on the other hand, discourage some others (for example, discourage selecting higher frequency components introduced by clipping and thus not likely to be present in the original signal).

Three possible approaches to constructing the vector of weights $\mathbf{w}$ are proposed below.

IV-A Absolute Threshold of Hearing

The human audio perception ranges from 20 Hz to ca 20 kHz. It is also commonly known that the human ear is more sensitive to frequencies around 2 – 5 kHz compared to the higher and lower frequencies [26]. This phenomenon was first characterized in 1933 by Fletcher and Munson as the Equal-loudness contours and the most recent definition from 2003 is listed in the standard ISO226:2003 [27]. The equal-loudness contours indicate the frequency dependency of the sound pressure level of pure tone at a given frequency that is perceived by humans as loud as 1 kHz pure tone associated with the same contour. The minimal loudness at which a harmonic sound is perceived is called the Absolute Threshold of Hearing (ATH) [28] and it can be approximated using the following equation [29]:

[TABLE]

where $g=f/1000$, i.e. it is the frequency in kHz.

The main idea of using the ATH for weighting the coefficients in (6) is to eliminate the negative effects of clipping, especially at frequencies where the human hearing is most sensitive. To do so, large weights should correspond to frequencies with the low respective ATH values and vice versa. There is a plethora of ways how to assign a weight to a corresponding frequency, based on the ATH curve. We examine the following three options:

[TABLE]

where $\mathbf{t}$ represents the vector of the ATH values for equispaced frequencies computed according to (10), and $\tau$ is the parameter setting the maximum value of the ATH in dB. Notice that $\mathbf{w}_{\text{ATH3}}$ is basically $\mathbf{w}_{\text{ATH2}}$ converted from dB(SPL) to the acoustic pressure in Pa. All the operations in (11) are conducted element-wise.

As the last step, the weights are always peak-normalized, so that the highest value of the weights is 1. The normalization does not influence the result, but it affects the speed of convergence and also the proper setting of the parameter $\gamma$ in Alg. 1.

Normalized weights computed from the ATH according to (11) with $\tau=100$, along with the original ATH contour, are shown in Fig. 1.

IV-B Global masking threshold

Apart from the ATH, another phenomenon responsible for human auditory perception is simultaneous masking, where the presence of a certain spectral component with high energy (the masker) masks another spectral component with lower energy (the maskee) [26]. The masker makes the maskee inaudible to humans, although the component is physically present.

Combining this phenomenon with the ATH gives the global masking threshold (GMT), which is a curve that indicates the minimum sound pressure level that a spectral component has to possess in order not to be masked by other spectral components. In other words, spectral components below the GMT will be evaluated as imperceptible.

The information contained in the GMT can be used to focus on restoring perceptually important components of the signal, while tolerating less audible components because they are masked and thus not perceived. Consequently, the weights should be constructed in a similar way to the case of ATH in the previous section, i.e. low values of GMT should produce large weights and vice versa. The same three possibilities shown in (11) were used, only the GMT plays the role of $\mathbf{t}$ now. The resulting vectors of weights are $\mathbf{w}_{\text{GMT1}}$, $\mathbf{w}_{\text{GMT2}}$ and $\mathbf{w}_{\text{GMT3}}$.

The GMT itself is computed by a slightly modified MPEG Psychoacoustic Model 1, where all the maskers are treated as the so-called tonal maskers [30]. In brief, the GMT (for a fixed time frame of the input signal) is obtained as follows: First, the signal is weighted with the Hann window and the power spectral density (PSD) is estimated using the FFT. From the PSD the maskers are identified, and the respective individual masking thresholds are calculated. Finally, the GMT is obtained as a power-additive combination of the ATH and the just described individual thresholds. An example of the GMT, the corresponding PSD estimate and the constructed weights are illustrated in Fig. 2.

Recall that the ground-truth signal is not known in practice. Computing the GMT from the observed (i.e. clipped) signal $\mathbf{y}$ may yield a biased estimate of the significant spectral components, especially for harsh clipping thresholds. Therefore, a way must be found to obtain a more correct GMT for use in (6). We tackle this problem by first computing the declipping problem by simple $\ell_{1}$ minimization without weighting, and we use the recovered signal as the basis for the GMT estimation and generation of the weights in (6).

IV-C Simple parabola

The third option is based on the idea that most of the energy in audio signals is concentrated at lower frequencies and that clipping introduces artificial higher harmonics that were not present in the original signal. Consequently, we consider the option of constructing the weights in such a way that it will suppress the higher harmonics and preserve the low frequencies.

A simple and effective approach to addressing this issue is to weight the coefficients with the identity curve. Better restoration results are obtained when a second-order polynomial is used instead of the identity. Formally, these weights are obtained as $\mathbf{w}_{\text{P}}=\mathbf{m}\odot\mathbf{m}$, where $\mathbf{m}=[1,\dots,M]$, with $M$ denoting the number of frequency channels of the DGT.

V Experiments and Results

The experiments are designed to evaluate and compare the above-presented possibilities of incorporating psychoacoustical information into audio declipping.

A test set of 10 musical audio signals sampled at 44.1 kHz with an approximate length of 7 seconds was thoroughly selected from the EBU SQAM database111https://tech.ebu.ch/publications/sqamcd to be diverse in tonal character and sparsity with respect to the time-frequency transform used. This way, the test set covers excerpts from the glockenspiel to the wind ensemble. Each sound example was first peak-normalized and then clipped, using multiple clipping thresholds $\theta_{\text{c}}=0.1,\dots,0.9$.

Diverse sound excerpts were used with numerous clipping thresholds, and our experience is that to obtain the best declipping results it would be necessary to tune the parameters of the DGT and the parameter $\gamma$ of the DRA manually for each test instance. For simplicity, these parameters were set to a compromise among all the cases. Specifically, Hann window 8192 samples long (i.e. 186 ms) with 75% overlap and 8192 frequency channels were used for the DGT. Such a setting generates a Parseval frame. The restoration algorithm was applied to a whole signal (i.e. the processing was not done frame-by-frame). Inside the DRA, the parameter $\gamma$ was set to 1 and the number of iterations was set strictly to 1000. Such a number provided convergence in each of the test instances.

As a reference, the proposed approaches are compared with the SPADE algorithm [17], which is, to the best of our knowledge, the state-of-the-art among the sparsity-based methods. Unlike the proposed approach, SPADE is designed to process the signal frame-by-frame. On purpose, our SPADE uses the same window shape and length and number of frequency channels. Internal parameters of the SPADE were $r,s=1$ and $\epsilon=0.1$.

V-1 Evaluation using signal-to-distortion-ratio (SDR)

The physical quality of restoration, i.e. how much the recovered signal $\hat{\mathbf{x}}$ is similar to the ground truth $\mathbf{x}$, is evaluated by the $\Delta\text{SDR}$, which is computed as a difference between the SDR values of the restored and the clipped signal, $\Delta\text{SDR}=\text{SDR}(\mathbf{x},\hat{\mathbf{x}})-\text{SDR}(\mathbf{x},\mathbf{y})$, and thus it expresses the SDR improvement in dB. The SDR for two signals $\mathbf{u}$ and $\mathbf{v}$ is computed as

[TABLE]

Fig. 3 illustrates the average $\Delta\text{SDR}$ values for all the proposed choices of weights (see Sec. IV), in comparison with the pure non-weighted version of Algorithm 1 and with the SPADE algorithm. The resulting $\Delta\text{SDR}$ values indicate that weighting with the GMT is a better approach than just a simple ATH curve in all the cases. Also, the best variant of converting the GMT or ATH curves into the actual weight vector $\mathbf{w}$ seems to be the one using the inversion, Eq. (11a). Nevertheless, among all choices, the best results (in terms of $\Delta\text{SDR}$) are obtained by weights derived from a parabola.

Comparing the $\ell_{1}$ model with the SPADE, the weights based on the parabola perform slightly better for mild clipping thresholds ($\theta_{\text{c}}\geq 0.5$) while SPADE performs slightly better for harsher clipping thresholds ($\theta_{\text{c}}<0.5$).

V-2 Evaluation using PEMO-Q

Since the topic of this paper is psychoacoustically motivated, the goal is to recover a signal that would sound the best to the listener; this requirement is not always quite in correspondence with the similarity in terms of the SDR. In other words, the $\Delta\text{SDR}$ is not the right objective measure and thus we use the PEMO-Q [31] evaluator, which should be more consistent with the human auditory system. The PEMO-Q output scale—ODG (the objective difference grade)—ranging from 0 to $-4$ is interpreted as:

[TABLE]

The PEMO-Q ODG values for the experiment are shown in Fig. 4. A comparison of Figs. 3 and 4 confirms the assertion that $\Delta\text{SDR}$ does not correspond quite well to human perception (as modeled by the PEMO-Q). The PEMO-Q results indicate that weights based on the ATH do not help to enhance the quality of restoration compared with the plain $\ell_{1}$ minimization. Weighting based on the GMT, on the contrary, improves the ODG values of restored signals, with the option (11a) being the best out of the three and even outperforming the SPADE algorithm for $\theta_{\text{c}}\geq 0.3$. The best restoration quality, according to PEMO-Q ODG, is delivered by the $\ell_{1}$ minimization weighted with parabola.

The source codes are available at http://www.utko.feec.vutbr.cz/~rajmic/software/declip_psychoacoust.zip.

VI Conclusion

In this paper, a novel sparsity-based approach to audio declipping has been proposed that incorporates psychoacoustic information. The results indicate that weighting helps to improve the obtained quality of restoration not only in SDR but more importantly in the PEMO-Q ODG scale, which better corresponds to human audio perception. The best results were obtained by weighting with the inverse GMT curve and parabola, which even outperform the SPADE algorithm in terms of the PEMO-Q ODG, while being ca $2.5\times$ faster.

As a subject of future research, the combination of the GMT with parabola could be promising since both different options produce good results. We also plan to involve the psychoacoustic weighting into social sparsity based algorithms and into the SPADE algorithm.