Overcomplete Frame Thresholding for Acoustic Scene Analysis
Romain Cosentino Randall Balestriero Richard Baraniuk Ankit Patel
eess.AS cs.SD stat.ML







[
    Gane Samb Lo ^1,2,3
    December 30, 2023  File: 
=============================




In this work, we derive a generic overcomplete frame thresholding scheme based on risk minimization. Overcomplete frames being favored for analysis tasks such as classification, regression or anomaly detection, we provide a way to leverage those optimal representations in real-world applications through the use of thresholding. We validate the method on a large scale bird activity detection task via the scattering network architecture performed by means of continuous wavelets, known for being an adequate dictionary in audio environments.






§ INTRODUCTION



From anomaly prediction in medical data to image recognition, the input signal is often corrupted with noise <cit.>. A classical approach in signal processing is based on the following analysis-synthesis framework: the data are analyzed (change of basis), denoised in their new basis, and then either fed to a machine learning classifier or synthesized (inverse the change of basis) before being processed by another algorithm. There is a wild number of analysis-synthesis operator, each of them being well defined mathematically by their specific properties. In this work, we will focus on overcomplete frame operators. As a matter of fact, those operators have been applied and studied in a wide range of areas such as wireless communications <cit.>, signal processing <cit.>, quantum computing <cit.>, functional analysis <cit.>. In fact, they yield interesting properties such as robustness <cit.>, stability and compactness <cit.>. 
A simple intuition behind those properties can be trivially described via the human vocabulary analogy. Consider the atoms of a dictionary as the vocabulary of a person. The higher the cardinal of the dictionary is, the higher the person will be precise and concise in his saying. As such, the analysis of his discourse will be eased thanks to its sparsity. In signal processing, sparse representations provide a fruitful pre-processing step allowing simple machine classifiers to produce high accuracy results. 
In this work, we will focus on overcomplete frames yielding a time-frequency representation as having an efficient representation of a signal is a synonym of extracting both its time and frequency content <cit.>. Especially, wavelet frames are known to provide a robust time-frequency representation for non-stationary signals as it is localized, both, in time and frequency. 

In most cases, wavelets can be used for three different tasks. First, for compression purposes: in fact, it has been demonstrated that wavelet transforms provide sparse representations particularly adapted for compression of natural signals with non-stationary components <cit.>. Secondly, for denoising purposes. Due to the sparsity of the induced representation and the ability to derive analytical solutions, optimized denoising solution with reconstruction guarantees have been developed <cit.>. Finally, the last task corresponds to analysis in which one aims at detecting subtle patterns of the input signal by using over-complete dictionaries. The analysis task can be considered today as crucial, especially when considering the need for pre-processing for current machine learning tasks. Unfortunately, thresholding methods for overcomplete frames have been tackled for special cases.




Therefore, we propose to remedy this absence by deriving a general analytical thresholding scheme for overcomplete frames principled via empirical risk minimization being computationally efficient and with the ability to leverage overcomplete frameworks within a corrupted environment. 

The paper is organized as follows: <ref> and <ref> are devoted to the related work and contribution of the paper,  the section <ref> shows the theoretical results including the derivation of the over-complete hard-thresholding technique, derivation of error bounds and presentation of computational complexity and pseudo-code. Finally, in section <ref> we provide experiments demonstrating the ability of the proposed scheme via the introduction of a Sparse Deep Scattering Network (SDSN), an extension of the Deep Scattering Network (DSN) <ref>. The task comparing the DSN and SDSN consists of a bird detection task using the Freefield1010[http://machine-listening.eecs.qmul.ac.uk/bird-audio-detection-challenge/]  audio scenes dataset in <ref>. The appendix is divided into two parts, Appendix <ref> provides both, the pre-requisite and details about the wavelets families used in our experiments;  In Appendix <ref>, we propose a review of the orthogonal denoising framework proposed by Donoho et. al. as well as the  mathematical details and proofs for the over-complete thresholding non-linearity that we introduce.



 §.§ Related Work


During the last decades, the properties of wavelets have been exploited to perform denoising on a wide range of signals <cit.>. In fact, the decomposition of a signal by mean of a wavelet transform concentrates the signal of interest in few coefficients having high energy, while the noise is uniformly propagated throughout all the time-frequency components <cit.>. This property is the direct consequence of the number of vanishing moment of the considered wavelet. In fact, it is defined as their orthogonality with respect to a certain order p of regular signals y  ∈𝒞^p. This property makes them an edge-detector basis (in the broad sense of the term) capable to provide sparse representations.
Donoho et. al. via a statistical signal processing approach to the problem, were able to provide powerful orthogonal basis denoising techniques holding theoretical guarantees <cit.>. However, orthogonal wavelet decomposition provides poor time-frequency resolution and are more prone to instabilities.
Since overcomplete frames provide extremely precise and often sparse representations, different approaches have been taken to tackle special cases of an overcomplete dictionary. This ability to fit the signal such closely should allow for more precise and robust denoising capacities.
An extension to the work of Donoho et. al., one step ahead toward overcomplete frames, has been achieved by Berkner et. al. in <cit.> using Translation Invariant Discrete Wavelet Transform and Biorthogonal-DWT. Then, using independent Component Analysis, Marusic et. al. <cit.> were able to tackle denoising in the case of multiple DWT. Others are exploiting the tools available in sparse coding-based models to perform denoising in learned "overcomplete frames" <cit.>. Therefore as these work rely on specific structures of overcomplete dictionaries that are often based on orthogonal filters drastically limit their applications. 

 


 §.§ Contributions


The contribution of the paper is the derivation of a general analytical thresholding technique for overcomplete frames. This technique based on empirical risk minimization will be studied by providing error bounds, computational complexity as well as pseudo-code allowing one to implement the method. We demonstrate the generality and robustness of our approach by mean of redundant continuous wavelet dictionaries. First, we propose the visualization of the thresholded scalograms as well as the evaluation of the sparsity before and after denoising. Secondly, we propose quantitative experiments validating the proposed scheme. Those later are achieved via introducing our thresholding technique into the Deep Scattering Network (DSN), which leads to the development of the Sparse Deep Scattering Network (SDNS). We compare the performance of the SDNS and the DSN on a bird activity detection task. As we will see, bringing denoising into this state-of-the-art pre-processing technique allows to further increase its linearisation capability and improves the accuracy of the classifier for a corrupted dataset composed of extremely diverse audio scenes. 



§ OVERCOMPLETE FRAMES DENOISING


In the case of analysis-synthesis framework with an overcomplete dictionary, the first question that arises is whether the mapping from the domain to the co-domain preserves distinctness. For sake of generality, we will assume that such operator has a nontrivial kernel, thus the application is not injective. Therefore, the dual of the dictionary is achieved via the pseudo-inverse transform. Notice that, Berkner et. al. in <cit.> proposed the use of the Moore-Pensore inverse to build the reconstruction dictionary of Bi-orthogonal wavelets and TIDWT dictionaries. 
Let's assume the observed signal, denoted by y, is corrupted with white noise such that y= x +ϵ where x is the signal of interest and ϵ∼𝒩(0,σ^2).
We now denote by W ∈ℂ^(n × (J × Q)) × n the matrix composed by the the filters at each time and scale (i.e: explicit convolution matrix of the operator 𝒲) such that ∀ x ∈ℝ^n, μ(x) = Wx is the wavelet transform of x. In other world, for each time position and each scale, the filters vectors ψ_i^T ∈ℂ^n,  ∀ i ∈{1,...,n×(J× Q)} are the row vectors of W. Then, we denote by W^†∈ℂ^ n × ((J × Q) × n ) the generalized inverse of W such that  W^†W=I. Each column of W^†, are the inverse filters ψ^†_i ∈ℂ^n,  ∀ i ∈{1,...,n×(J× Q)}. The estimate of x is given by  x̂_W,D(y) = W^†D^SWy where D^S is a diagonal binary operator such that,
 
    D^S_i,i= δ_i = {
     1     if    i ∈ S  
     
     0     if    i ∈ U 
    .

with U and S denoting respectively the set of selected and unselected wavelet coefficients. We also define D^U such that I = D^U + D^S. Our thesholding estimate is based on the mean-square error between the signal of interest x and the synthesized thresholded observed signal y: 

    ℛ^⋆(x,W)   =  min_δ𝔼 x-x̂_W,D(y) ^2 =  min_δ𝔼 W^† (Wx-D^SWy ^2,

 which can be rewritten in the following explicit form,

    ℛ^⋆(x,W) =    min_δ∑_i,j=1^n × (J × Q)μ_iμ_jψ^†^T_iψ^†_j 1_{δ_i = 0, δ_j = 0 }
        + σ^2∑_i,j=1^n × (J × Q)  ( ψ^†^T_iψ^†_j ) ψ_i^Tψ_j 1_{δ_i =1,δ_j =1 },

where μ =  Wx are the co-domain's coefficients. It is clear that the correlation inherent to the redundant information contained in the dictionary implies an intractable factorial optimization problem as the ideal risk is now dependent on all the possible pairs in the frequency axis. Since this ideal risk requires the intervention of an oracle, we develop an upper bound to the ideal risk such that we benefit an explicit equation for the thresholding operator that is adapted to any over-complete transformation and possesses a tractable and analytic expression. We propose to use a min-max formulation as an upper bound of this ideal risk explicitly derived in Appendix <ref>. 
The upper-bound on the optimal risk is denoted by ℛ_up and defined as,

    ℛ_up(x,W) = ∑_k =1^n(J*Q+1)min(ℛ_up^U(x)[k],ℛ_up^S[k]),

where we denote by ℛ_up^U the upper bound error term corresponding to unselected coefficients:

    ℛ_up^U(x)[k] =∑_j =1^n × (J × Q) | μ_k(x)μ_j(x)    ψ^†^T_kψ^†_j |,

which means that if we do not select a coefficient (e.g: μ_k), the risk is equal to a weighted value of this coefficient. The weight are all the coefficient (μ_j, ∀ j) weighted by the correlation between the reconstruction filter  (ψ^†_k) and each of the reconstruction filters (ψ^†_j , ∀ j). 
Now, let's ℛ_up^S[k] be the upper bound error term corresponding to the selected coefficients:

    ℛ_up^S[k] = σ^2∑_j=1^n × (J × Q) |  (ψ^†^T_kψ^†_j) ψ_k^Tψ_j   |.

In the case where the coefficient is selected, the error corresponds to the energy of the noise propagated by the redundant frame. As a matter of fact, the noise is weighted by the correlation of, both the filter ψ^T_k and the reconstruction filters ψ^†_k, with all filters ψ_j^T and reconstruction filters ψ^†_j, ∀ j.
Now, one way to evaluate this upper-bound, is to assume an orthogonal basis, and to compare it with the optimal risk in the orthogonal case which leads to the following proposition.

Assuming orthogonal filter matrix W_O, the upper bound ideal risk coincides with the orthogonal ideal risk:

    ℛ_up(x,W_O) = ℛ_O(x,W_O).


The proof is derived in <ref>. Therefore, even with the upper-bound, our optimal risk, in the case of bijective transform, can still be used and the optimal solution can be recovered. 
In order to apply the ideal risk derive, ones needs an oracle decision regarding the signal of interest. In real application, the signal of interest x is unknown. We thus propose the following empirical risk:

    ℛ̃(y,W) =  ∑_k =1^n × (J × Q)min(ℛ_up^U(y)[k],ℛ_up^S[k]).

This risk corresponds to the empirical version of the ideal risk where the observed signal y is evaluate in the left part of the minimization function. In order to compare this empirical risk with the ideal version, we propose the bound analysis with respect to the optimal risks in the two possible extreme cases:

In the case where D^S=I, the empirical risk coincides with the upper bound ideal risk:

    ℛ̃(y,W) =   ℛ_up(x,W).




In the case where D^U=I, the following bound shows the distance between the empirical risk and the upper bound ideal risk:

    ℛ̃(y,W) ≤ℛ_up(x,W) + ∑_k,j=1^n × (J × Q) C[k,j] × |  ψ^†^T_kψ^†_j |, a.s.

where, 

    C[k,j] =      | μ_k(x)  | ψ_j _1 σ√(2/π) +   | μ_j(x)  |ψ_k _1 σ√(2/π)  + σ^2 (1- 2/π).


Refer to <ref> for proofs.
As the empirical risk introduces the noise in the left part of the risk expression, this term represents the propagation of this noise throughout the decomposition. Naturally, we from the proposition <ref>, that the more redundant the dictionary is, the less the bound is tight. 



 §.§ Algorithm and Complexity

In this section we present practical guide and discussion concerning the application of the method for large scale tasks.
In term of complexity, by using the proposed local method, we are asymptotically linear in time hence number of signals to treat and quadratic in the number of filters present in the over-complete filter-bank. However, the latter can be ”pushed out” of the computations by performing it once for all signals. In fact, most of the costly terms to compute are constant for all treated signals (the pseudo inverse of the filters correlation). As a result the most demanding operation remains in the computation of the correlation matrix to be done for each column of the scalogram. Again this is quadratic in the number of filters but linear in time and with computations completely independent from one another allowing very efficient parallelization or map-reduce schemes. 





§ EXPERIMENTAL VALIDATION WITH CONTINUOUS WAVELET TRANSFORMS:
BIRD ACTIVITY DETECTION TASK

We propose to validate the two contributions over a large-scale audio dataset. 
The data set consists of 7000 field recording signals from  freefield1010[http://machine-listening.eecs.qmul.ac.uk/bird-audio-detection-challenge/] collected via the Freesound[https://arxiv.org/abs/1309.5275] project. This collection represent a wide range of audio scenes such as birdsong, city, nature, people, train, voice, water...
The focus in this paper is the bird audio detection task that can be formally defined as a binary classification task, where each label corresponds to the presence or absence of birds. Each signal is 10sec. long, and has been sampled at 44.1Khz. The evaluation of the results is performed via the Area Under Curve metric on 33% of the data. The experiments are repeated 50 times.
The total audio length of this dataset is thus of slightly less than 20 hours of audio recordings. To put in comparison, it is about 20 × larger than CIFAR10 in term of number of scalar values in the dataset. The results comparing our algorithm to the DSN are achieved with two different wavelet and are provided in Table <ref>.  For all the experiments, the octave and quality parameters of the layers are J1=5,Q1=8,J2=4,Q2=1. As the feature of interests are birds songs, only high frequency content requires high resolution, the thresholding is applied per window of 2^16 representing ≈ 1.5 sec.


 §.§ Continuous Wavelet Transform

We first denote the collection of scaling factors by, 

    λ = {λ_j, j= 1,...,J × Q },  λ_j=2^(j-1)/Q,
 
The hyper-parameters are, J and Q which represents respectively the number of octave to decompose and  the quality coefficients a.k.a the number of wavelets per octave.
Based on the collection Λ, the wavelet filter-bank can be derived by scaling the mother wavelet ψ_0 leading to the filter-bank denoted as

    𝒲[ψ_0,λ]=[     ψ_λ_1;         …; ψ_λ_J × Q;           ],  ψ_λ_j(t)=1/√(λ_j)ψ(t/λ_j),

where we omit the upper-script over the λ parameter to avoid redundancy. Also, we denote this filter-bank as 𝒲[ψ_0,λ]:=𝒲. 
We define the continuous wavelet transform as,

    y ⋆𝒲= [  y ⋆𝒲(1,.);          …; y ⋆𝒲(JQ,.) ].





 §.§ Qualitative Analysis

we first provide some examples of denoising over different audio scenes.
In the Fig. <ref>, we present the visualization and evaluation of the sparsity within the wavelet domain for several signals. The evaluation of the sparsity is achieved by the following sparsity ratio:

    α =  1- x _0/n × (J× Q),

where . _0 denotes the l_0 norm. It is clear that the closer the ratio is to 1, the sparser the representation is. On the other hand, if the ratio equals zeros, none of the coefficients within the representation are null. 


        
    

Given the fact that the data set contains recording via mobile application in diverses audio scenes, it is clear that many signals contain energy within all frequency bands. As such, the wavelet decomposition, even if non-orthogonal, is not sparse in most cases. Our denoising technique has the ability to remove the background noise and conserves only the foreground activity of each signal. In the special case of noise activity, one can see in Fig.  <ref>, that our algorithm conserves the bird song's chirp while removing almost all the other information. 




 §.§ Quantitative Analysis: Sparse Deep Scattering Network


We now propose a quantitative analysis via solving the supervised classification task by using a variant of the standard scattering network including our proposed denoising scheme.
The Deep Scattering Network, first developed in <cit.> and first successfully applied in <cit.> is based on a cascade of linear and non-linear operators on the input signal. The linear transformation is a wavelet transform, and the nonlinear transformation is a complex modulus. At each layer, the convolution of the representation with a scaling function leads to the so-called scattering coefficients.
This network is stable (Lipschitz-continuous) and suitable for machine learning tasks as it removes spatio-temporal nuisances by building space/time-invariant features.
The translation invariance property is provided by the scaling function that acts as an averaging operator on each layer of the transform leading to an exponential decay of the scattering coefficients <cit.>. Since the continuous wavelet transform increases the number dimension, the complex modulus is used as its contractive property reduces the variance of the projected space <cit.>. 


The first layer of a DSN corresponds to standard scalogram, the second considers each frequency dimension of this scalogram as the input signal of the second layer. All the features are computed by applying a low-pass filter on each of these representations. In the SDSN, a non-linear thresholding operator is applied on each representation before to process both, the next layer as well as the scattering coefficient of the current layer.





  §.§.§ Sparse Deep Scattering Network


We denote by U^(1) the application of the filter-bank of the first layer onto a signal y corresponding to the first layer of the scattering network. The output of this first layer and as previously mentioned consists of a 2D tensor of shape (J^(1)Q^(1),N) with N the length of the input signal.
We omit here boundary conditions, sub-sampling, and consider a constant shape of N throughout the representations.
We thus obtain,

    U^(1)[y]=|y ⋆𝒲^(1)|,

where |.| operator corresponds to an element-wise complex modulus application.
As previously mentionned, we define the convolution operation between those two objects as,

    y ⋆𝒲^(1)= [          y ⋆𝒲^(1)(1,.);                      …; y ⋆𝒲^(1)(J^(1)Q^(1),.) ].

Then, we define by 𝒯 the thresholding operator minimizing the empirical risk,

    𝒯 = min_δℛ̃(y,W^(1)).

This non-linearity is applied to the latent representation U^1 producing the following sparse representation: 𝒯[U^(1)[y]]:=U^(1)_T[y].
From this sparse representation, the scattering coefficients can be computed leading to 

    S^(1)  = U^(1)_T[y] ⋆ϕ^(1),

ϕ^(1) being a ”low-pass” filter. The application of the latter is performed for each row of the representation.
This application of a low frequency band-pass filter allows for symmetries invariances, inversely proportional to the cut-off frequency of ψ.
From this, the second layer is computed. 
We now denote the second layer representation as U^(2)[U_T^(1)[y]]  being
 a 3D tensor. 
 It defined as

    U^(2)[U_T^(1)[y]](j_2,.,.)= | U^(1)_T[y] ⋆𝒲^(2)(j_2)|.

Then, the thresholding operator is applied to each new representation U^(2)[U_T^(1)[y]](.,j_1,.) independently. The final sparse representation is thus denoted by  U_T^(2)[U_T^(1)[y]].
Given this 3D representation, the second layer scattering coefficients are defined as

    S^(2)=U_T^(2)[U_T^(1)[y]] ⋆ϕ^(2),

We present an illustration of the network computation in Fig.  <ref>.

As can be seen in the proposed example, while the first layer provides time-frequency information, the second layer characterizes transients as demonstrated in <cit.>. With this extended framework, we now dive into the problem of thresholding redundant frame, cases where the quality factor, Q, is greater than 1 which are in practice needed to bring enough frequency precision. 
We provide in Fig.  <ref> illustration showing the effect of each SDSN layers. It is clear that even if the first layer representation has been denoised, the second layer's latent representation, by decomposing each of the U^(1)_T frequential dimension, still contains inherent noise. This noise is removed as it can be observed in the second thresholded layer representation U^2_T. 


As we will demonstrate below, our method as well as each contribution taken independently and jointly lead to significant increase in the final performance. We compare our results against the DSN. In all cases, the scattering coefficients are then fed into a random forest <cit.> with parameters[n_estimator: 100, min_samples_split:  130 ,class_weights:'balanced_subsample']  based on the sklearn library <cit.>.






  §.§.§ Controlled Denoising Experiment

We based our implementation on <cit.> leveraging the Fourier based computations of localized filter in the frequency domain. For all input signals we perform a renormalization to have unitary energy. 
We propose complementary experiments in this section to further demonstrate the need for thresholding even in the context of stable representations offered by the scattering network. To do so we simulate an environment for which we combine denoised and noised samples in both the training and test sets with different mixture proportions.
Thus we propose in table <ref> the described experiment achieved via the Gammatone wavelet only: for both the SDSN and DSN models. 
We evaluate the performance of the classifier for different mixture configuration as we can see in the table <ref>. The row corresponds to the percentage of noisy data in the training set, and the columns are the percentage of noisy data in the test set. It is clear that the percentage of denoised data in each set is one minus the percentage of noisy data. Then, in the table <ref>, we show the difference between the mean accuracy of the SDSN-gammatone with respect to each mixture mean accuracy. 


The configuration depicted in the diagonal of the tables <ref> and <ref> corresponds to the ”ideal” case of identical proportion of noise in the train and test set. However for many practical applications more volatile regimes can be witnessed. Those correspond either to the upper diagonal part for which the proportion of noisy samples is more important in the test set than the train set. Conversely the lower diagonal corresponds to the symmetric case. We can see that even in the ideal setting, introduction of the SDSN brings reasonable robustness to the classifier pipeline. More importantly for the off diagonal cases, a significant gain is achieved via the SDSN.






plain


§ BUILDING A DEEP SCATTERING NETWORK




 §.§ Continous Wavelet Transform



"By oscillating it resembles a wave, but by being localized it is a wavelet".
 Yves Meyer 


Wavelets were first introduced for high resolution seismology <cit.> <cit.> and then developed theoretically by Meyer et al. <cit.>.
Formally, wavelet is a function ψ∈𝕃^2 such that:

    ∫ψ(t) dt = 0,

it is normalized such that ψ_𝕃^2 = 1. 
Often, wavelets can be categorized in two distinct groups, the discrete wavelets and the continuous ones. The discrete wavelets are constructed based on a system of linear equation that represent the atom's property. These wavelet, when scaled in a dyadic fashion form an orthonormal basis.
On the other hand, the continuous wavelets have an explicit formulation and build an over-complete dictionary when successively scaled. In this work, we focus on the continuous wavelets as they provide a more complete tool for analysis of signals. 
In order to perform a time-frequency transform of a signal, we first build a filter bank based on the mother wavelet. This wavelet is names the mother wavelet since it will be dilated and translated in order to create the filters that will constitute the filter bank. Notice that wavelets have a constant-Q property, thereby the ratio bandwidth to center frequency of the children wavelets are identical to the one of the mother. Then, the more the wavelet atom is high frequency the more it will be localized in time.
The usual dilation parameters follows a geometric progression and belongs to the following set:

    Λ = { 2^(j-1)/Q, j= 1,...,J × Q  }.
 
Where the integers J and Q denote respectively the number of octaves, and the number of wavelets per octave.
Having selected a geometric progression ensemble, the dilated version of the mother wavelet in the time are computed as follows:   
    
    ψ_λ(t)=1/λψ ( t/λ ),    ∀λ∈Λ,

and can be calculated in the Fourier domain as follows:

    ψ̂_λ(ω) = ψ̂(λω),  ∀λ∈Λ.

Notice that in practice the wavelets are computed in the Fourier domain as the wavelet transform will be based on a convolution operation which can be achieved with more efficiency in the frequency plane. 
By construction, the children wavelets have the same properties than the mother one. As a result, in the Fourier domain:

    ψ̂_λ = 0,  ∀λ∈Λ.

Thus, to create a filter bank that cover all the frequency support, one needs a function that captures the low frequencies contents. The function is called the scaling function and satisfies the following criteria:

    ∫ϕ(t) dt = 1.










 §.§ Wavelet Families


Among the continuous wavelets, different families exist. Each posses different properties, such as bandwidth, center frequency. This section is dedicated to the development of the families that are important for the analysis of diverse signals.


  §.§.§ The Morlet wavelet

The Morlet wavelet (Fig.  <ref>) is built by modulating a complex exponential and a Gaussian window defined in the time domain by,

    ψ^M(t) = π^-1/4 e^iω_0 te^-t^2/2,

where ω_0 defines the frequency plane.
In the frequency domain, we denote it by ψ̂^M(t),

    ψ̂^M(ω) = π^-1/4e^-(ω-ω_0)^2/2, ∀ω∈ℝ^⋆_+,

thus, it is clear that ω_0 defines the center frequency of the mother wavelet.
With associated frequency center and standard deviation denoted respectively by ω_c^λ_i and Δ^λ_iω, ∀ j ∈{0,...,J*Q-1 } are:

    ω_c^λ_i   = ω_0/λ_i, 
    Δ^λ_iω   = 1/2 λ_i^2.

Notice that for the admissibility criteria ω_0 = 6, however one can impose that zeros-mean condition facilely in the Fourier domain. Usually, this parameter is assign to the control of the center frequency of the mother wavelet. Then, we are able to vary the parameter ω_0 in order to have different support of Morlet wavelet. 

The Morlet wavelet, is optimal from the uncertainty principle point of view <cit.>. The uncertainty principle, when given a time-frequency atoms, is the area of the rectangle of its joint time-frequency resolution. In the case of wavelet, given the fact that their ratio bandwidth to center frequency are equal, it implies that this area is equal for the mother wavelets and its scaled versions. As a result, because of its time-frequency versatility this wavelet is wildly used for biological signals such as bio-acoustic <cit.>, seismic traces <cit.>, EEG <cit.> data.



  §.§.§ The Gammatone wavelet

The Gammatone wavelet is a complex-valued wavelet that has been developed by <cit.> via a transformation of the real-valued Gammatone auditory filter which provides a good approximation of the basilar membrane filter <cit.>. Because of its origin and properties, this wavelet has been successfully applied for classification of acoustic scene <cit.>.






The Gammatone wavelet (Fig.  <ref>) is defined in the time domain by,

    ψ^G(t) =  ( 2 π(i-σ)t^m-1+(m-1)t^m-2 )e^-2piσ te^2 pi i t,

and in the frequency domain by,

    ψ̂^G(ω) =  i ω (m-1)!/ (σ +i(ω-σ)  )^m.

A precise work on this wavelet achieved by V. Lostalnen in <cit.> allows us to have an explicit formulation of the parameter σ such that the wavelet can be scaled while respecting the admissibility criteria:

    σ^2 =  r^2/m (1 - r^2/m) m^2 ξ^2/2 ( √(1+ B^2/(1 - r^2/m)^2 m^2 ξ^2) -1  ),

where ξ is the center frequency and B is the bandwidth parameter. Notice that B=(1-2^-1/Q) ξ with ξ= 2 π/1+ 2^1/Q induce a quasi orthogonal filter bank.
The associated frequency center and standard deviation denoted respectively by ω_c^λ_i and Δ^λ_iω, ∀ j ∈{0,...,J*Q-1 } are thus:

    ω_c^λ_i   = ξ, 
    Δ^λ_iω   = B .

For this wavelet, thanks to the derivation in <cit.>, we can manually select for each order m the center frequency and bandwidth of the mother wavelet, which ease the filter bank design.

An important property that is directly related to the auditory response system is the asymmetric envelope, thereby the Gammatone wavelet is not invariant to time reversal to the contrary of the Morlet wavelet that behaves as a Gaussian function. Thus, for tasks such as sound classifications, this wavelet provides an efficient filter that will be prone to perceive the sound's attacks. Beside this suitable property for specific analysis, this wavelet is near optimal with respect to the uncertainty principle. Notice that, when m →∞ it yields the Gabor wavelet <cit.>. Another interesting property of this wavelet is the causality, by taking into account only the previous and present information, there is no bias implied by some future information and thus it is suitable for real-time signal analysis.



§ THRESHOLDING




 §.§ Thresholding in an orthogonal basis

Assuming that the observed signal y, is corrupted with white noise, 
 
    y = x+ ϵ,

 where ϵ is a vector of i.i.d centered normal distributions 𝒩(0,σ^2). 
 Now we define the estimate of x by x̂_W,D such that:
 
    x̂_W,D(y) = W^TD^SWy

 where W denotes the orthogonal basis and D^S is a diagonal binary operator such that,
 
    D^S_i,i= δ_i = {
     1     if    i ∈ S  
     
     0     if    i ∈ U 
    .
,
 where U and S denote respectively the set of selected and unselected wavelet coefficients. We also define D^U such that I = D^U + D^S.
 This estimate corresponds to a thresholding operation in the new basis and the inverse transform of this truncated representation. 

 We define the denoising problem as the solution of the following mean-square error:
 
    ℛ_o^⋆(x,W)   =  min_δ𝔼 x-x̂_W,D(y) ^2 
        =   min_δ𝔼 W^T (Wx-D^SWy ^2 
        =  min_δ𝔼D^U Wx - D^SW(x+ϵ) ^2 
        =  min_δ D^U Wx  ^2 + σ^2 tr (D^S WW^T D^S) 
        =  min_δ∑_i^n [Wx]_i^2 1_{δ_i = 0 } + σ^2 1_{δ_i =1}
        =  ∑_i^nmin ([Wx]_i^2,σ^2).

 Therefore, the optimal D^S^⋆ and D^U^⋆ given by the following δ values:

    δ_i =  1_ | [Wx]^2_i  |>σ^2.



 §.§ Upper-bound Non-orthogonal Risk & Empirical Risk




    x̂_W,D(y) = W^†D^SWy


    ℛ^⋆(x,W)   =  min_δ𝔼 x-x̂_W,D(y) ^2 
        =   min_δ𝔼 W^† (Wx-D^SWy ^2 
        =  min_δ𝔼W^† (D^U Wx - D^SWϵ) ^2 
        =  min_δ W^⋆D^U Wx  ^2 + σ^2 tr (W^TD^S W^†W^†^T D^SW),

Developing the previous expression and denoting by μ =  Wx the wavelet coefficient vector, we have:

    ℛ^⋆(x,W) =    min_δ∑_i,j=1^n × (J × Q)μ_iμ_jψ^†^T_iψ^†_j 1_{δ_i = 0, δ_j = 0 }
        + σ^2∑_i,j=1^n × (J × Q)  ( ψ^†^T_iψ^†_j ) ψ_i^Tψ_j 1_{δ_i =1,δ_j =1 }.

we first use the triangular inequality,


    ℛ^⋆(x,W)     ≤min_δ∑_i,j=1^n × (J × Q) |     μ_iμ_jψ^†^T_iψ^†_j | 1_{δ_i = 0, δ_j = 0 }
        +  σ^2∑_i,j=1^n × (J × Q) | (ψ^†^T_iψ^†_j) ψ_i^Tψ_j   | 1_{δ_i =1,δ_j =1 }

Now let's,

    ℛ^U[i,j]=    |    μ_iμ_j ( ψ^†^T_iψ^†_j)    |,

and,

    ℛ^S[i,j]= σ^2 |(ψ^†^T_iψ^†_j) ψ_i^Tψ_j   |.

Then, based on the following min-max formulation, we obtain an upper bound of the ideal risk, that, when minimized will approximate the ideal risk in the overcomplete case:

    ℛ^⋆(x,W)     ≤∑_k =1^n × (J × Q)min_δ_k∑_j=1^n × (J × Q)max_δ_j, j≠ k (1_{δ_k = 0, δ_j = 0 }ℛ^U[k,j] + 1_{δ_k =1,δ_j =1 }ℛ^S[k,j] ) 
       ≤∑_k =1^n × (J × Q)min_δ_k∑_j=1^n × (J × Q) (max_δ_j, j≠ k 1_{δ_k = 0, δ_j = 0 }ℛ^U[k,j]  +  max_δ_j, j≠ k  1_{δ_k =1,δ_j =1 }ℛ^S[k,j] ) 
        =  ∑_k =1^n × (J × Q)min_δ_k   1_{δ_k =0 } [ ∑_j =1^n × (J × Q) |   μ_kμ_jψ^†^T_kψ^†_j | ] 
         + 1_{δ_k =1 } [ σ^2∑_j=1^n × (J × Q) | (ψ^†^T_kψ^†_j) ψ_k^Tψ_j   |  ].

Now, let's denote by ℛ_up^U the error term corresponding to unselected coefficients:

    ℛ_up^U[k] =∑_j =1^n × (J × Q) | μ_kμ_jψ^†^T_kψ^†_j |,

and by ℛ_up^S for the selected ones:

    ℛ_up^S[k] = σ^2∑_j=1^n × (J × Q) |  (ψ^†^T_kψ^†_j) ψ_k^Tψ_j   |.

we have that,

    ℛ_up(x,W)     =  ∑_k =1^n × (J × Q)min_δ_k   1_{δ_k =0 }ℛ_up^U[k]  + 1_{δ_k =1 }ℛ_up^S[k]  
        =  ∑_k =1^n × (J × Q)min(ℛ_up^U[k],ℛ_up^S[k]).




 §.§ Comparison Upper Bound Ideal Risk with Orthogonal Ideal Risk


Proposition 1.
The comparison of this upper bound risk given an orthogonal dictionary and the one derived in the orthogonal case is as follows:

 If the basis is orthogonal, we have,

    (ψ^†^T_kψ^†_j) =  1,    k=j 
     0,

and,

    ψ_k^Tψ_j =   1,    k=j 
     0,

Therefore, the upper-bound derived recovers the ideal risk in the orthogonal case.




 §.§ Comparison Upper Bound Ideal Risk with Empirical Risk


Proposition 2.





If D^S=I, the empirical risk is equal to:

    ℛ̃(y,W) = σ^2∑_k,j=1^n × (J × Q) |  (ψ^†^T_kψ^†_j) ψ_k^Tψ_j   |.
   
and the upper bound risk is:   

    ℛ_up(x,W) = σ^2∑_k,j=1^n × (J × Q) |  (ψ^†^T_kψ^†_j) ψ_k^Tψ_j   |.

Thus both coincide as this restriction on the support of the risk makes it independent of both x and y.
Proposition 3.
In the case where D^U=I,

    ℛ̃(y,W)    = ∑_k,j=1^n × (J × Q) | μ_k(y) μ_j(y)  ψ^†^T_kψ^†_j | 
        = ∑_k,j=1^n × (J × Q) | ( μ_k(x) μ_j(x) +  μ_k(x) μ_j(ϵ) +  μ_k(ϵ) μ_k(x) + μ_k(ϵ) μ_j(ϵ) )  | ×
          | ψ^†^T_kψ^†_j |,

by the triangular inequality, we have that:

    ℛ̃(y,W)≤   ∑_k,j=1^n × (J × Q) (  | μ_k(x) μ_j(x)  | + |  μ_k(x) μ_j(ϵ)  | +   | μ_k(ϵ) μ_k(x) | +  | μ_k(ϵ) μ_j(ϵ) )  | ×
            |  ψ^†^T_kψ^†_j |,

by the monotony of expectation and the Fubini theorem, we have almost surely:


    ℛ̃(y,W)≤ℛ_up(x,W) +  ∑_k,j=1^n × (J × Q) C[k,j] × |   ψ^†^T_kψ^†_j |      a.s.,

where C[k,j] is equals to,

    C[k,j] =     | μ_k(x)  | ψ_j _1 σ√(2/π) +   | μ_j(x)  |ψ_k _1 σ√(2/π)  + σ^2 (1- 2/π).

