Khanh Mai

Member, IEEE Dominique Pastor

Abdeldjalil Aissa

El Bey

Raphaël Le Bidan

Student Member, IEEE Van-Khanh Mai

Senior Member, IEEE Abdeldjalil Aïssa-El-Bey

Member, IEEE Raphaël Le-Bidan

Robust Estimation of Non-Stationary Noise Power Spectrum for Speech Enhancement

Keywords: Speech enhancement, noise power spectrum estimation, noise reduction, robust statistics

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

N OWADAYS communication electronic support in gen- eral and telephone conversation in particular often take place in noisy and non-stationary environments such as the inside of a car, in the street or inside an airport for example. Hence many research efforts have aimed at improving not only the quality but also the intelligibility of speech. Noise power spectrum estimation is a key issue in designing robust noise reduction methods for speech enhancement. Most noise power spectrum estimation algorithms found in the literature can be classified into four main categories [START_REF] Loizou | Speech enhancement: theory and practice[END_REF], namely histogram-based methods, minimal-tracking algorithms, time-recursive averaging algorithms, and other techniques derived from Maximum-Likelihood (ML) or Bayesian estimation principles, e.g. minimum mean square error (MMSE) methods.

In the first category of algorithm, the noise power spectrum is estimated from the maximum of the histogram in the time-frequency domain of the observed signal power spectrum, the latter being determined by using a first-order smoothing recursion [START_REF] Hirsch | Noise estimation techniques for robust speech recognition[END_REF]. An improvement of this method involves updating the noise power spectrum uniquely on the frames detected as noise-only by a chi-square test [START_REF] Ahmed | A voice activity detector using the chisquare test[END_REF]. However, most of the histogram-based algorithms have the drawback of being relatively complex in terms of computational cost and memory resources [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF].

In the second family of methods, the noise power spectrum is tracked via minimum statistics, according to the reasonable hypothesis that the noise power spectrum level is below that of noisy speech [START_REF] Martin | Spectral subtraction based on minimum statistics[END_REF], [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF]. First, the smoothed noisy speech power spectrum is evaluated by a first-order recursive operation. Then, the noise variance is computed as the statistical minimum of the smoothed power spectrum with a factor of correction. The main difference between the two methods in [START_REF] Martin | Spectral subtraction based on minimum statistics[END_REF] and [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF] lies in the computation of the smoothing parameter used in the first order recursion. In [START_REF] Martin | Spectral subtraction based on minimum statistics[END_REF], the smoothing parameter is chosen empirically, whereas this parameter is derived by minimizing the mean square error between the noise and the smoothed noisy speech power spectrum in [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF]. Minimum-statistics methods require observing the noisy signals on a sufficiently long time interval so as to track speech power instead of noise power. On the other hand, a long time interval is detrimental to the quality of the estimate in case of non stationary noise. A trade-off is thus necessary, leading to a typical time-delay of 1 to 3 seconds in practice. This causes underestimation which decreases in turn the performance of noise reduction algorithms.

Famous methods of the third category include the Minima-Controlled Recursive-Averaging (MCRA) algorithm [START_REF] Cohen | Noise estimation by minima controlled recursive averaging for robust speech enhancement[END_REF] and its many modifications such as the Improved-MCRA (IMCRA) [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF] or the MCRA2 [START_REF] Rangachari | A noise-estimation algorithm for highly non-stationary environments[END_REF] methods. In this class of algorithms, the noise power spectrum in a given frequency bin is estimated by first-order recursive operations where smoothing parameters depend on the conditional speech presence probability in the bin. The main difference between MCRA, MCRA2 and IMCRA lies in the way the speech-presence probability is estimated. MCRA and MCRA2 directly estimate the speech-presence probability frame-by-frame via a smoothing operation whereby, for a given frame, the probability of speech presence is increased when this frame is detected as noisy speech and decreased otherwise. A frame is detected as noisy speech if the ratio of the smoothed noisy speech power spectrum to its local minimum is above a certain threshold, the local minimum being computed by using the minimum-statistics technique proposed in [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF]. Fixed and frequency-dependent thresholds are used in MCRA and MCRA2, respectively. On the other hand, IMCRA derives the speech-presence probability in each bin by a two-step estimation of the speech-absence probability. The first iteration aims at detecting the absence of speech in a given frame, while the second iteration actually estimates the speech-absence probability from the power spectral components in the speech-absence frame. The main disadvantage of these methods is the estimation delay in case of sudden rising noise, this delay being mainly due to the use of the minimum-statistics methods of [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF].

Techniques derived from ML or Bayesian estimation principles overcome the problem of sudden rising noise by estimating the noise power spectrum from the noise periodogram via a statistical criterion. In [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF], [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF], the noise instantaneous power is evaluated by MMSE and then incorporated in a recursive noise power spectrum estimation technique. [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] proposes a simple bias compensation of the noise instantaneous power before estimating the noise power spectrum via the same recursive smoothing and under the same hypotheses as in [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF]. However, the noise instantaneous power estimate in [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] remains biased. In contrast, an unbiased estimator for the noise instantaneous spectrum is obtained in [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] by soft-weighting the noisy speech instantaneous power and the previous noise power spectrum estimate by the conditional probabilities of speech-absence and speech-presence, respectively. The noise power spectrum estimation can also be carried out by recursive ML-Expectation-Maximization [START_REF] Souden | Noise power spectral density tracking: A maximum likelihood perspective[END_REF], similar to MCRA and IMCRA. This approach allows for rapid noise power spectrum estimation and tracking by avoiding the use of minimum-statistics methods.

In this paper, we propose a new approach for noise power spectrum estimation, without requiring any model or any prior knowledge for the probability distributions of the speech signals. Fundamentally, we do not even take into consideration the fact that the signal of interest here is speech. The approach is henceforth called extended-DATE (E-DATE) since it basically extends the d -dimensional amplitude trimmed estimator (DATE), initially proposed in [START_REF] Pastor | Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences[END_REF] for white gaussian noise (WGN), to colored stationary and non-stationary noise. The main principle at the heart of the E-DATE algorithm is the weak-sparseness property of the STFT of noisy signals, according to which the sequence of complex values returned by the STFT in a given timefrequency bin can be modeled as a complex random signal with unknown distribution and whose unknown probability of occurrence in noise does not exceed one half. Noise in each bin is assumed to follow a zero-mean complex gaussian distribution [2, p. 210], so that estimating the noise power spectrum amounts to estimating the noise variance in each bin, the latter being provided by the DATE. The DATE trims the amplitudes in each given bin, after having sorted them by increasing norm. Noise power spectrum estimation by E-DATE is thus similar to and actually extends the quantile-based approach of [START_REF] Stahl | Quantile based noise estimation for spectral subtraction and wiener filtering[END_REF], which relies on assumptions that the weak-sparseness model embraces. More generally, the reader will notice similarities between the proposed method and the state-of-the-art techniques mentioned above. A main difference between the E-DATE approach and standard ones is actually the mathematical justification of the former via the weak-sparseness model, which formalizes more or less standard heuristics in speech processing and yields a reduced number of parameters for more robustness. Although the E-DATE does not rely on minimum-statistics principles or methods, it does however require a time buffer having the same length -typically 80 frames for a frequency sampling rate of 8 kHz -as other popular algorithms.

The paper is organized as follows. In Section II, the main features of the DATE are reviewed. Section III develops the weak-sparseness model for noisy speech. The E-DATE is then introduced in Section IV, following a step-bystep methodology where we successively deal with WGN, stationary noise and non-stationary noise. Two practical implementations of the E-DATE algorithm are then described. The performance of the E-DATE algorithm is evaluated in Section V and compared to state-of-the-art methods in terms of number of parameters and estimation errors. Speech enhancement experimental comparisons using objective as well as pseudo-subjective criteria are also conducted by combining the noise power spectrum estimation methods with a noise reduction system. Conclusions are finally given in Section VI.

II. THE DATE

For the sake of self-completeness, this section presents the DATE in its full generality. Given d -dimensional observations of random signals randomly absent or present in independent and additive WGN, the purpose of the DATE is to estimate the noise standard deviation. Such an estimation may serve to detect the signals or estimate them as in speech denoising. As in [START_REF] Davies | The identification of multiple outliers (with discussion)[END_REF], the DATE addresses the frequently-encountered case where 1) most observations follow the same zero-mean normal distribution with unknown variance, 2) signals of interest have unknown distributions and occurrences in noise. Standard robust scale estimators such as the very popular median absolute deviation (MAD) estimator and the trimmed estimator (Testimator) have performance that degrades significantly when the proportion of signal increases. In contrast, the DATE can still estimate the noise standard deviation when possible signals occur with a probability too large for usual scale estimators to perform well. As indicated by its name, the DATE basically trims the norms of the d -dimensional observations. However, in contrast to the conventional Testimator, which applies to one-dimensional data and fixes the number of outliers to remove, the DATE applies to any dimension and chooses adaptively the number of outliers to discard. It performs the trimming by assuming that the signal norms are above some known lower-bound and that the signal probabilities of occurrence are less than one half. These assumptions bound our lack of prior knowledge about the signals and make it possible to separate signals from noise. Moreover, these assumptions are suitable for signal processing applications where noisy signals are considered as outliers with respect to the noise distribution. They are particularly suitable for observations obtained after sparse transforms capable of representing signals by coefficients that are mostly small except a few ones whose norms are relatively big. In particular, the sequel will exhaustively use the fact that the Fourier transform of speech signals is sparse in a weak sense detailed hereafter.

The DATE basically relies on [1, Theorem 1], which is asymptotic and can be viewed as a method of moments. A detailed presentation of the theoretical background of the DATE is beyond the scope of this paper and the reader is referred to [START_REF] Pastor | Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences[END_REF] for details. However, the following brief heuristic presentation may be convenient for the reader. This heuristic exposure departs from that proposed in [1, Theorem 1], so as to shed different light on the theory behind the DATE.

[U ∈ B ] = 1 if U ∈ B and 1[U ∈ B ] = 0, otherwise.
In addition, Γ is the standard Gamma function and 0 F 1 is the generalized hypergeometric function [15, p. 275]. All the random vectors and variables are henceforth assumed to be defined on the same probability space (Ω, P, E).

Let (Y n ) n∈N be a sequence of d -dimensional random observations such that: In this model, each observation is either noise alone or the sum of some signal and noise. The probability distributions of the signals Λ n are supposed to be unknown. Our purpose is then to estimate σ 0 .

If all the ratios Λ n /σ 0 are known to be above some sufficiently large signal to noise ratio (SNR) ρ, it can be expected that some threshold height σ 0 ξ(ρ) can suitably be chosen to decide with small error probability that Λ n is present (resp. absent) whenever Y n is above (resp. less) σ 0 ξ(ρ). Therefore, most of the non-zero terms in the sum

N n=1 Y n 1[ Y n σ 0 ξ(ρ)] should pertain to noise alone. If the number N n=1 1[ Y n σ 0 ξ(ρ)
] of these non-zero terms is itself large enough, we should have an approximation of the form

N n=1 Y n 1[ Y n σ 0 ξ(ρ)] N n=1 1[ Y n σ 0 ξ(ρ)]
≈ λ σ 0 . Such an approximation can actually be proved asymptotically with the help of some additional assumptions. More precisely, suppose that:

(A1) Λ n , X n and ε n are independent for every n ∈ N; (A2) the set of priors P[ε n = 1] : n ∈ N is upper-bounded by 1/2 and the random variables ε n , n ∈ N, are independent;

(A3) sup n∈N E Λ n 2 < ∞.
These assumptions including (A0) deserve some comments.

To begin with, the independence assumption in (A0) is mainly technical to prove the results stated in [START_REF] Pastor | Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences[END_REF]. In fact, our experimental results below suggest that this assumption is not so constraining in speech processing, where we deal with non-overlapping but not necessarily independent time frames. Assumption (A1) simply means that the two hypotheses for the observation occur independently and that the noise and signal are independent. The model thus assumes prior probabilities of presence and absence through the random variables ε n . However, the impact of these priors is reduced by assuming that the probabilities of presence and absence are actually unknown. The role of Assumption (A2) is then to bound this lack of prior knowledge about the occurrences of the two possible hypotheses that any Y n is supposed to satisfy. Assumption (A3) simply means that the signals Λ n have finite energy.

Under assumptions (A0)-(A3) and with the help of [16, Theorem 1], [1, Theorem 1] then guarantees that σ 0 is the unique positive real number σ such that:

lim ρ→∞ lim sup N →∞ N n=1 Y n 1[ Y n σξ(ρ)] N n=1 1[ Y n σξ(ρ)] -λσ ∞ = 0 (1)
where λ = 2 Γ d +1 2 /Γ d 2 and ξ(ρ) is the unique positive solution in x to the equality 0 F 1 (d /2; ρ 2 x 2 /4) = e ρ 2 /2 . It is thus natural to estimate the noise standard deviation σ 0 by seeking a possibly local minimum of:

N n=1 Y n 1[ Y n σξ(ρ)] N n=1 1[ Y n σξ(ρ)] -λ σ , (2)
when σ ranges over some search interval [σ min , σ max ]. Given a lower bound ρ for the ratios Λ n /σ 0 , the DATE computes the solution in σ to the equality:

N n=1 Y n 1[ Y n σξ(ρ)] N n=1 1[ Y n σξ(ρ)] = λ σ. ( 3
)
Indeed, such a solution trivially minimizes (2). In addition, an application of Bienaymé-Chebyshev's inequality makes it possible to determine the value n min ∈ {1, 2, . . . , N } such that the probability that the number of observations due to noise alone be above n min is larger than or equal to some given probability value Q. The main steps of the DATE are summarized in Algorithm 1, where

Y (1) , Y (2) , . . . , Y (N ) is the sequence Y 1 , Y 2 , . . . , Y N sorted by increasing norm so that Y (1) Y (2) . . . Y (N )
, and where we have defined

M * { Y 1 , Y 2 ,..., Y N } (n) =    1 n n k=1 Y (k) if n = 0 0 if n = 0, (4)
The parameters on which the DATE relies are thus: the dimension d of the observations, the number N of observations and the lower bound ρ for the possible SNRs. The two parameters that directly influence the DATE performance are N and ρ. As recommended in [1, Remark 4], we can use ρ = 4 in practice. Theoretically, N should be large since the theoretical result on which the DATE relies is asymptotic by nature. However, experimental results show that the DATE performance is acceptable when N is above 200. This will Algorithm 1 DATE algorithm for estimation of noise standard deviation Input:

• A finite subsequence {Y 1 , Y 2 , . . . , Y N } of a sequence Y = (Y n ) n∈N of d -dimensional real random vectors satisfying assumptions (A0-A3) above • A lower bound ρ for the SNRs Λ n /σ 0 , n ∈ N • A probability value Q 1 - N 4(N /2-1) 2 Constants: n min = N /2 -N /4(1 -Q), ξ(ρ), λ Output: The estimate σ * {Y 1 ,Y 2 ,...,Y N } of σ 0 Computation of σ * {Y 1 ,Y 2 ,...,Y N } : Sort Y 1 , Y 2 , . . . , Y N by increasing norm so that Y (1) Y (2)
. . . Y (N ) if there exists a smallest integer n in {n min , . . . , N } such that:

Y (n) M * { Y 1 , Y 2 ,..., Y N } (n)/λ ξ(ρ) < Y (n+1) n * = n else n * = n min end if σ * {Y 1 ,Y 2 ,...,Y N } = M * { Y 1 , Y 2 ,..., Y N } (n * )/λ
be confirmed by the application to speech processing of Sections IV and V. Another means to choose the minimal SNR required by the DATE is to resort to the notion of universal threshold [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], as proposed in [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF]. Indeed, the coordinates of all the

N observations Y 1 , Y 2 , . . . ,
(N × d ) = σ 0 2 ln(N × d ) could
then be regarded as the maximum absolute value of these gaussian random variables when N × d is large. Instead of proceeding as in wavelet shrinkage [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] where the universal threshold is utilized to discriminate noisy signal wavelet coefficients from wavelet coefficients of noise alone, the trick proposed in [START_REF] Atto | Detection thresholds for nonparametric estimation[END_REF] and [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF] is to consider λ u (N × d ) as the minimum amplitude that a signal must have to be distinguishable from noise. The minimal SNR can then be defined as

ρ = ρ(N × d ) = λ u (N × d )/σ 0 = 2 ln(N × d ).
It is an interesting fact that the value of ρ(N ×d ) grows rapidly to 4 with N ×d .

In the sequel, we will consider values returned by STFT. The DATE will therefore be applied to sequences of real and complex values, that is, one-and two-dimensional data since complex values can be regarded as 2-dimensional real vectors. It is thus worth recalling the specific values of ξ(ρ) and λ for d = 1 and

d = 2. If d = 1, ξ(ρ) = cosh -1 (e ρ 2 /2 ) = 1 2 ρ + 1 ρ log(1 + 1 -e -ρ 2 ) and λ = 0.7979. If d = 2, ξ(ρ) = I -1 0 (e ρ 2 /
2 )/ρ where I 0 is the zeroth order modified Bessel function of the first kind and λ = 1.2533. Note that 1/λ can be regarded as a bias correction factor, similar to those employed by minimum-statistics approaches.

III. WEAK-SPARSENESS MODEL OF NOISY SPEECH

The main motivation for utilizing the DATE is that noisy speech signals in the time-frequency domain after STFT reasonably satisfy the same type of weak-sparseness model as used to establish [START_REF] Pastor | Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences[END_REF]Theorem 1]. This weak-sparseness model essentially assumes that the noisy speech signal can be represented by a relatively small number of coefficients with large amplitudes. Indeed, let us consider the spectrograms of Figure 1 obtained by STFT of typical examples of clean and noisy speech signals. In the time-frequency domain, speech is composed of a set of time-frequency components or atoms. Most atoms with small amplitudes are masked in the presence of noise. Only the few atoms whose amplitude is above some minimum value remain visible in noise. Clearly, the proportion of these significant atoms does not exceed one half. These remarks lead to the following model for noisy speech STFTs. In the time domain, the observed signal is given by

y(t ) = s(t ) + x(t ) (5)
where s(t ) and x(t ) denotes clean speech and independent additive noise. Note that both are real-valued signals. The signal in the time domain is transformed into the timefrequency domain by STFT since most noise reduction systems operate in this particular transform domain. Hence, all processing is frame-based. Let K be the frame length, or equivalently, the STFT length. The corresponding system model in the time-frequency domain then reads:

Y (m, k) = S(m, k) + X (m, k) (6)
in which m denotes the frame index, k is the frequency-bin index, and S(m, k) (resp. X (m, k)) stands for the STFT component of the speech signal (resp. noise) at time-frequency point (m, k). Following [2, page 210], we model each X (m, k) as a complex Gaussian random variable. Complex values Y (m, k) are manipulated as 2-dimensional real vectors. According to the empirical remarks above, the weaksparseness model first assumes that an atomic speech audio source is either present or absent at any given timefrequency point (m, k). The presence or the absence of this source is modeled by a Bernoulli random variable ε(m, k). This Bernoulli model is tantamount to and justified by the concept of ideal binary masking in the time-frequency domain, as used in audio source separation [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF], [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. The probability of presence is assumed to be less than or equal to 1/2. Thus P ε(m, k) = 1 ≤ 1/2. Second, the atomic audio source must have significant amplitude so as to contribute effectively to the mixture that composes the speech signal. The minimum amplitude that such a source must have will hereafter be denoted by ρ. Let us further denote by Θ(m, k) the underlying atomic audio source. Then, under the previous assumptions, the noisy speech signal at timefrequency point (m, k) can be modeled as: We recognize here the weak-sparseness model [START_REF] Pastor | Wavelet shrinkage: from sparsity and robust testing to smooth adaptation[END_REF] applied to speech processing, in the continuation of [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF]. In summary, our model essentially assumes that the STFT of noisy speech signals satisfies the following three key properties in each time-frequency bin (m, k): (A'1): the presence/absence of speech ε(m, k) and the atomic speech audio source Θ(m, k) are independent, (A'2): the speech-presence probability does not exceed 1/2, (A'3): the instantaneous power of the random clean speech signal is upper-bounded by a finite value.

Y (m, k) = ε(m, k)Θ(m, k) + X (m, k) (7)
Assumptions (A'1-A'3) are adaptations of (A1-A3) to the particular case of noisy speech signals. Regarding (A0), its equivalent form for noisy speech signals is simply Eq. ( 7).

Our purpose is then to estimate the noise power spectrum σ 2 X (m, k) = E |X (m, k)| 2 at any given time-frequency point (m, k). This problem is similar to that addressed in [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF], where the signal of interest was a mixture of audio signals, possibly including speech signals, and where additive noise was stationary, gaussian and white. The DATE was used to estimate the noise power spectrum in [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF] because this estimator does not make prior assumption on the statistical nature of the signals of interest. In the present paper and in contrast to [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF], we do not restrict our attention to WGN and generalize the approach of [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF] to the estimation of colored and possibly non-stationary noise in the presence of speech.

IV. NOISE POWER SPECTRUM ESTIMATION BY E-DATE

In this section, we derive the E-DATE algorithm that will be used for noise power spectrum estimation in all the experiments conducted in Section V. The derivation follows a three-step process, which aims at gradually introducing the modifications required to evolve from the academic WGN model to the much more realistic, but also more challenging, practical case of non-stationary noise. More precisely, we first describe the application of the DATE algorithm to noise power spectrum estimation of noisy speech signals in the time-frequency domain. We extend the DATE to the case of colored stationary gaussian noise, and then discuss the estimation of non-stationary noise. This leads to the E-DATE algorithm, which is specifically designed for noise power spectrum estimation in non-stationary noisy environments, but can be used with stationary noise as well.

In the following, we suppose to be given M noisy speech frames of K samples. The frames are assumed to be nonoverlapping so as to satisfy assumption (A0). The STFTs are normalized by 1/ K .

A. Stationary WGN

In this case, the noise power spectrum is constant and equals σ 2 X over the whole time-frequency plane. Accordingly, and by properties of the (normalized) STFT, each noise sample X (m, k) in the time-frequency domain is a zero-mean circularly-symmetric gaussian complex random variable with variance σ 2 X :

X (m, k) ∼ N c (0, σ 2 X ).
Equivalently, X (m, k) may be viewed as a zero-mean twodimensional real gaussian random vector with covariance matrix (σ 2 X /2)I 2 :

X (m, k) ∼ N 0, (σ 2 X /2)I 2 .
Since the STFT of noisy speech signals is weakly-sparse in the sense of Section III, the M × (K /2 -1) values Y (m, k) for m ∈ {1, 2, . . . , M } and k ∈ {1, 2, . . . , K /2 -1} can be used as inputs of the two-dimensional (d = 2) version of the DATE to provide an estimate σ 2 X of σ 2 X . Note that, in principle, another estimate of σ 2 X could be obtained by applying a one-dimensional (d = 1) DATE on the 2

× M real dataset Y (1, 0), Y (2, 0), . . . , Y (M , 0), Y (1, K /2), Y (2, K /2), . . . , Y (M , K /2).
However, the size of this second dataset is usually much smaller than that of the first one. Thus only the first option is used in practice as it leads to a more reliable estimate. Note also that, due to the Hermitian property of the STFT of real input signals, |Y (m, k)| = |Y (m, K -k)|. Therefore the frequency bins K /2 + 1 to K are not used in the estimation process as they do not bring additional information.

B. Colored stationary noise

For colored stationary noise, the noise power spectrum is no longer constant over the whole time-frequency plane but may vary as a function of frequency. Consequently, each noise sample X (m, k) in a given frequency bin k will now be modeled as a zero-mean complex gaussian random variable with variance σ 2 X (k):

X (m, k) ∼ N c 0, σ 2 X (k) .
Here again, the STFT output sequence Y (m, k) for m = 1, 2, . . . , M is assumed to be weakly-sparse in the sense of Section III so that in each frequency bin k, only a few of these values will have an SNR above ρ and in a proportion that does not exceed 1/2. As a result and as illustrated in Figure 2, the extension to colored stationary noise involves running concurrently K /2 + 1 independent instances of the DATE to estimate σ 2 X (k) in each frequency bin k = 0, 1, 2, . . . , K /2. As discussed earlier, we do not use ) DATE whereas the one dimensional (d = 1) DATE is used for bins 0 and K /2. For colored noise, assumption (A'1) may not always rigorously hold, especially at low frequencies. However, as supported by the experimental results of Section V, this deviation with respect to the underlying theoretical model turns out to be no real issue in practice, thanks to the robust behavior of the DATE, even when the signal presence probability may exceed 1/2 (see [1, Figure 2]).

In contrast to WGN for which the whole time-frequency plane (≈ M K /2 observations) is used to estimate the noise variance σ 2 X , M frames only are available here to estimate σ 2 X (k) in each frequency bin. Clearly a more reliable estimate can be obtained by increasing M , but this increases in return the overall computational cost and may also entail some time-delay. A possible solution is to begin with a first estimate σ 2 X (k) computed over the first M frames, and then to periodically update this estimate as new frames are acquired. For stationary noise, the initial number of frames M need not be very high. Even if the first estimate is not very accurate, it is expected to improve rapidly as new frames enter the estimation process.

C. Extension to non-stationary noise: The E-DATE algorithm

Most practical applications including speech denoising usually face a mix of stationary as well as non-stationary noise. Unlike white or colored stationary noise, the power spectrum of non-stationary noise varies over time and frequency, and, as such, proves to be much more challenging to estimate. Interestingly, non-stationary noise models including car noise, babble noise, exhibition noise and others, usually exhibit some form of local stationarity in time and frequency. In such cases, non-stationary noise can be considered as approximately stationary within short time periods of D consecutive frames, where parameter D has to be defined appropriately for each noise model. This amounts to assuming the existence of a noise power spectrum in this time interval, which is a function of frequency only. The DATE algorithm for colored stationary noise introduced in Section IV-B can then be used to estimate the noise power spectrum within this time window of D frames. This is the basis of the E-DATE algorithm.

Parameter D can be preset once for all or could be optimized for applications where prior knowledge about noise is available. The choice for duration D results from a tradeoff between estimation accuracy, stationarity and practical constraints such as computational cost and time-delay. A large value for D may violate the local stationary property. On the other hand, the number of frames D should be large enough to produce reliable estimates σ 2 X (k). In case D is too small to provide the DATE with a sufficient number of input data, a possible solution consists in grouping several consecutive frequency bins. This is tantamount to assuming that the noise power spectrum is approximately constant over those frequencies. Such a procedure however requires prior knowledge on the noise spectrum properties, which can be irrelevant in practical applications where noise has often unknown type and may evolve across time. For this reason, this solution will not be further studied below.

In summary, the E-DATE algorithm consists in carrying noise power spectrum estimation by running a per-bin instance of the DATE (see Figure 2) on periods of D consecutive non-overlapping frames, where D is chosen so that noise can be considered as approximately stationary within this time interval. Once an estimate of the noise power spectrum has been obtained, it can be used for denoising purpose for instance, but will not be taken into account in the computation of future estimates, as the local power spectrum of non-stationary noise may change significantly from one period of D frames to the next.

Although the E-DATE algorithm was specifically designed for power spectrum estimation of non-stationary noise, it can be used without modification for power spectrum estimation of WGN or colored stationary noise, thereby offering a robust and universal noise power spectrum estimator whose parameters are fixed once for all types of noise considered above. Let us now discuss the practical implementation of the E-DATE algorithm.

D. Practical implementation of the E-DATE algorithm

Two different implementations of the E-DATE algorithm are proposed here.

The first approach is a straightforward block-based implementation of the algorithm described in Section IV-C. It involves estimating the noise power spectrum on each period of D successive non-overlapping frames. This requires storing D frames, calculating the K /2 + 1 estimates σ 2 X (k) using the observations in these D frames, and then waiting for D new non-overlapping frames. The resulting algorithm is called Block-E-DATE (B-E-DATE) and summarized in Algorithm 2, where σ = DATE d ,ρ y 1 , y 2 , . . . , y n denotes the standard deviation estimate σ returned by the d -dimensional DATE with minimal SNR ρ and n real ddimensional inputs y 1 , y 2 , . . . , y n . Estimation of the noise power spectrum over separate periods of D non-overlapping frames reduces the overall algorithm complexity. However, this entails a time-delay of D frames, which must be considered in applications. Consider the particular example of speech denoising illustrated in Figure 3. Noise reduction is performed on a frame-by-frame basis. A new noise power spectrum estimate is provided to the noise reduction system by the B-E-DATE algorithm once every D non-overlapping frames, and then used to denoise each of those D frames. Clearly, denoising cannot start before the first D non-overlapping frames have been recorded. This results in an overall latency of about 1 or 2 seconds for typical sampling rates of 8 and 16 kHz. This delay can then have some impact for speech applications embedded in current mobile devices. It will naturally be lesser in applications such as Active Noise Cancellation (ANC) where frequency rates are much higher.

The delay limitation can be bypassed as follows. First, a standard noise power spectrum tracking method is used to estimate the noise power spectrum during the first D -1 non-overlapping frames. Any of the methods mentioned in the introduction can be used for this purpose. Afterwards, starting from the D th frame onwards, a sliding-window version of the E-DATE algorithm is used to estimate the noise spectrum on a per-frame basis, using the latest recorded D non-overlapping frames. This alternative implementation called Sliding-Window E-DATE (SW-E-DATE) is summarized in Algorithm 3. Its application to speech denoising is illustrated in Figure 4.

The B-E-DATE and the SW-E-DATE algorithm may be viewed as two particular instances of a more general bufferbased algorithm. More precisely, the B-E-DATE algorithm corresponds to the extreme case where the buffer is totally flushed and updated once every D non-overlapping frames. In contrast, the SW-E-DATE algorithm corresponds to the other extreme case where only the oldest frame is discarded in order to store the current one, in a First-In First-Out (FIFO) mode. Clearly, a more general approach between these two extremes consists in partially updating the buffer by renewing only L frames among D. This point has not been further investigated in the present work.

Note finally that the proposed implementations of the E-DATE algorithm are not limited to speech denoising but could find use in any application involving signals corrupted by additive and independent non-stationary noise, and to which the weak-sparseness model locally applies.

Algorithm 2 Block-Extended-DATE (B-E-DATE) algorithm for noise power spectrum estimation

for m ≥ D do if mod (m, D) = 0 m * = m σ X (m * , 0) = DATE 1,ρ Y (m -D + 1, 0), Y (m -D + 2, 0), . . . , Y (m, 0) σ X (m * , K /2) = DATE 1,ρ Y (m -D +1,K /2),Y (m-D +2,K /2),. . . ,Y (m,K /2) for k := 1 to N 2 -1 do σ X (m * , k) = DATE 2,ρ Y (m -D + 1, k), Y (m -D + 2, k), . . . , Y (m, k) σ X (m * , K -k) = σ X (m * , k) end for else σ X (m -D, k) = σ * X (m * , k) end if end for V. PERFORMANCE EVALUATION
Several comparisons and experiments were conducted in order to assess the performance and benefits of the E-DATE noise power spectrum estimator in comparison with other state-of-the-art algorithms. Both the B-E-DATE and the SW-E-DATE implementations were considered in two Algorithm 3 Sliding-Window Extended-DATE (SW-E-DATE) algorithm for noise power spectrum estimation

for m = 1 to the end of signal do if m < D Calculate σ X by an alternative method else σ X (m, 0) = DATE 1,ρ Y (m -D + 1, 0),Y (m -D + 2, 0), . . . ,Y (m, 0) σ X (m, K /2) =DATE 1,ρ Y (m-D +1,K/2,Y (m-D +2,K/2)), . . . ,Y (m,K/2) for k := 1 to K 2 + 1 do σ X (m, k) = DATE 2,ρ Y (m -D + 1, k), Y (m -D + 2, k), . . . , Y (m, k) σ X (m, K -k) = σ X (m, k)
end for end if end for different benchmarks. In subsection V-A, we first compare the number of parameters required by the E-DATE and several classical or more recent noise power spectrum estimators. Then, we compare in subsection V-B the estimation quality of the different algorithms in several distinct noise environments. The combination of the noise power spectrum estimation algorithms with a noise reduction system based on the Log-MMSE algorithm is investigated using the NOIZEUS speech corpus in subsection V-C. Finally, the time-complexity of the E-DATE algorithm is analyzed in subsection V-D.

A. Number of parameters

Table I gives the number of parameters required by the E-DATE as well as by the state-of-the-art noise power spectrum estimation algorithms mentioned in the introduction. Derived from robust statistical signal processing concepts, the E-DATE is the simplest algorithm to configure, with only two parameters to specify, namely the SNR lower bound ρ and the number of frames D. This stands in sharp contrast with other popular approaches such as Minimum Statistics [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF], which involves 7 parameters. In practice, the minimal SNR ρ can be set as explained at the end of Section II so that the only crucial parameter is D. Working with D = 80 nonoverlapping frames of K = 256 samples was found to yield good performance in all the experiments reported here.

B. Noise Estimation Quality

The estimation quality of the noise power spectrum estimation algorithms listed in Table I was evaluated on several noise models using the symmetric segmental logarithmic estimation error measure defined in [START_REF] Hendriks | Noise tracking using DFT domain subspace decompositions[END_REF]. The difference between the estimated noise power spectrum σ 2 X (m, k) and reference noise power spectrum σ 2 X (m, k) is evaluated by

Log Er r = 1 M K M -1 m=0 K -1 k=0 10 log 10 σ 2 X (m, k) σ 2 X (m, k) (8)
where M denotes the total the number of available frames.

For WGN, the theoretical reference noise power spectrum is known and can be substituted to σ 2 X (m, k) in ( 8). This is no longer the case for non-stationary noise involved in the NOIZEUS database. For non stationary noise, the reference noise power spectrum σ 2 X (m, k) is estimated as follows [START_REF] Hendriks | Noise tracking using DFT domain subspace decompositions[END_REF]:

σ 2 X (m, k) = ασ 2 X (m -1, k) + (1 -α)|X (m, k)| 2 , with α = 0.9.
Both the B-E-DATE and the SW-E-DATE implementations of the E-DATE algorithm were evaluated and compared. The SW-E-DATE uses the recently-introduced MMSE2 method [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] as a surrogate algorithm to provide an estimate for the first D-1 frames since, as shown below, this algorithm turns out to offer excellent performance among state-of-the-art noise estimators.

The Log Er r measures obtained with the different noise power spectrum estimators are given in Figure 5. All algorithms have been benchmarked at four SNR levels and against various noise models, namely WGN, auto-regressive (AR) colored stationary noise, and 6 typical non-stationary noise environments.

The results for white and colored stationary noise are given in Figures. 5(a) and 5(b), respectively. The B-E-DATE and SW-E-DATE methods yield the lowest Log Er r error, the best performance being achieved by the B-E-DATE algorithm in WGN. This is no surprise since the underlying DATE algorithm was originally developed for estimating the standard deviation of additive WGN.

For non-stationary noise with slowly-varying noise spectrum like exhibition, car, station or train noise, and depending on the noise level, the B-E-DATE algorithm uniformly obtains either the best score, or comes very close to the best score, as shown in Figures 5(c), 5(d) and 5(e), respectively.

Figures 5(f ), 5(g) and 5(h) present the results obtained with the least favorable types of non-stationary noise. In the case of modulate WGN (resp. babble noise), the SW-E-DATE (resp. B-E-DATE) algorithm yields the smallest Log Er r error. As illustrated in Figure 5(h), the two proposed algorithms are among the best in estimating the very challenging airport noise environment. Their performance closely match those obtained with the state-of-the-art MMSE2 and ML-EM estimators.

C. Performance Evaluation in Speech Enhancement

In complement to the previous study, the performance of the noise power spectrum estimation algorithms listed in Table I have also been evaluated and compared in combination with a noise reduction system. The speech denoising experiments are based on the NOIZEUS database [START_REF] Loizou | Speech enhancement: theory and practice[END_REF], which contains IEEE sentences corrupted by eight types of noise coming from the AURORA noise database, at four SNR levels, namely 0, 5, 10 and 15 dB. The noise reduction algorithm retained for our experiments is the Log-MMSE estimator [START_REF] Ephraim | Speech enhancement using a minimum mean-square error log-spectral amplitude estimator[END_REF]. This method is a standard reference in B-E-DATE SW-E-DATE MARTIN [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF] IMCRA [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF] MCRA2 [START_REF] Rangachari | A noise-estimation algorithm for highly non-stationary environments[END_REF] MMSE1 [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] MMSE2 [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] ML-EM [START_REF] Souden | Noise power spectral density tracking: A maximum likelihood perspective[END_REF] Fig. 6: SNRI with various noise types speech denoising. It can easily be implemented and is known to reduce residual noise without distorting too much the speech signal [START_REF] Loizou | Speech enhancement: theory and practice[END_REF]p. 230,Sec. 7.7].

Two different criteria have been used to compare the different algorithms. The first one is the Signal-to-Noise Ratio Improvement (SNRI) objective criterion standardized in the ITU-T G.160 recommendation for evaluating noise reduction systems [START_REF]Voice Enhancement Devices for Mobile Networks[END_REF]. The SNRI performance obtained with the Log-MMSE combined with the noise power spectrum estimators of Table I are shown in Figure 6 for various noise environments. Note that 4 noise levels were used for each noise type, the final SNRI score being computed as the average score over these 4 levels. We observe that the B-E-DATE and SW-E-DATE yield similar performance measurements and that they outperform all other methods for each type of noise except airport noise. The average SNRI score computed over the 11 noise types and labeled Total at the right of Figure 6 clearly emphasizes the SNRI gain brought by the E-DATE in comparison to other methods.

The second criterion used to assess noise power spectrum estimation in speech enhancement is the composite objective measures proposed in [START_REF] Hu | Evaluation of objective measures for speech enhancement[END_REF] (see also [START_REF] Loizou | Speech enhancement: theory and practice[END_REF]). This criterion introduces three measures C si g , C bak and C ovl that are linear combination of some widely used measures like segmental SNR (segSNR), weighted-slope spectral (WSS), log likelihood ratio (LLR), and perceptual evaluation of speech quality (PESQ): The three measures C si g , C bak and C ovl are designed so as to provide a high correlation with the three usual corresponding subjective measures that are signal distortion (SIG), background intrusiveness (BAK) and Mean Opinion Score (OVRL). We focus here on the C ovl criterion since it has the highest correlation with the real subjective tests. LogErr(dB) MARTIN [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF] MCRA2 [START_REF] Rangachari | A noise-estimation algorithm for highly non-stationary environments[END_REF] IMCRA [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF] MMSE1 [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] MMSE2 [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] ML-EM [START_REF] Souden | Noise power spectral density tracking: A maximum likelihood perspective[END_REF] B-E-DATE SW-E-DATE (a) WGN MARTIN [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF] MCRA2 [START_REF] Rangachari | A noise-estimation algorithm for highly non-stationary environments[END_REF] IMCRA [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF] MMSE1 [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] MMSE2 [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] ML-EM [START_REF] Souden | Noise power spectral density tracking: A maximum likelihood perspective[END_REF] B-E-DATE SW-E-DATE (b) AR noise MARTIN [START_REF]Noise power spectral density estimation based on optimal smoothing and minimum statistics[END_REF] MCRA2 [START_REF] Rangachari | A noise-estimation algorithm for highly non-stationary environments[END_REF] IMCRA [START_REF] Cohen | Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging[END_REF] MMSE1 [START_REF] Yu | A low-complexity noise estimation algorithm based on smoothing of noise power estimation and estimation bias correction[END_REF] MMSE2 [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF] ML-EM [START_REF] Souden | Noise power spectral density tracking: A maximum likelihood perspective[END_REF] B-E-DATE SW-E-DATE (f) modulated WGN Figure 7 shows the C ovl scores obtained with the different noise power spectrum estimators and noise environments. For reference purpose, the C ovl score obtained with noisy speech but without noise reduction is shown in dashed lines in each sub-figure. The good performance of the B-E-DATE and SW-E-DATE are confirmed by the C ovl measures obtained in the case of WGN, AR noise, car noise, station noise and train noise. These results allow us to conclude that the E-DATE approach is well-suited for stationary or slowly varying non-stationary noise. Although not shown here for space limitation, we hasten to mention that very similar trends were observed for the other two criteria C si g and C bak . In the challenging case of airport noise, all the methods in this paper introduce a large signal distortion at 0dB and 5 dB. At 10 and 15 dB, the E-DATE C ovl scores are similar to that obtained by the other methods (see Figure 7(h)). A detailed analysis of the C bak scores in babble and airport noise (see Figure 8) nevertheless reveals that the E-DATE algorithms perform best in terms of background noise reduction. Two final remarks are in order here. First, the B-E-DATE algorithm generally performs better than the SW-E-DATE algorithm. This is particularly evident in Figure 7 and can also be noticed in the other experimental results. This is mainly due to the fact that our implementation of the SW-E-DATE initially resorts to a surrogate algorithm to estimate noise power spectrum during the first D = 80 frames, which has inferior performance to the B-E-DATE. Since these D frames represent a significant part of the total duration of many of the tested utterances, the performance loss incurred by the use of a worse estimator significantly impacts the overall score. Second, in the previous section was evoked the possibility to partially update the buffer by renewing only L frames among D instead of flushing it completely (B-E-DATE), or renewing it only one frame at a time in a FIFO manner (SW-E-DATE). The difference in performance between these two E-DATE implementations suggests that such a partial renewal should not dramatically modify the results. This means that buffer optimization can be performed in practice whenever required by practical constraints, and without significantly impacting the de-

D. Complexity analysis

Tables II and III compare the computational costs of the B-E-DATE and SW-E-DATE implementations, respectively. Each table gives the number of real additions, multiplications, divisions and square roots required to perform the estimate. Both the B-E-DATE and the SW-E-DATE use D frames to compute the noise power spectrum estimate. However computation is performed only once every D frames for the B-E-DATE algorithm, whereas it is performed once per frame in the SW-E-DATE implementation. Hence the number of operations in Table II should be divided by D to allow for a fair per-frame computational cost comparison between the two implementations. For reference purpose, Table IV lists the number of operations required by the MMSE2 estimator of [START_REF] Gerkmann | Unbiased MMSE-based noise power estimation with low complexity and low tracking delay[END_REF]. Inspection of Tables II and IV shows that the B-E-DATE and MMSE2 estimators have similar computational complexity. This is confirmed by execution times of Matlab implementations of these algorithms where the B-E-DATE algorithm is found to have a processing time about 1.53 times that of the MMSE2 algorithm. We also note from Tables II and III that SW-E-DATE requires approximately D/3 times more operations that B-E-DATE. Indeed, B-E-DATE requires 3D multiplications to process D frames at once, whereas SW-E-DATE requires D +2 multiplications per frame. Execution times of Matlab implementations of these algorithms also confirm this ratio.

VI. CONCLUSION

In this paper, we have proposed a novel method to estimate the power spectrum of some non-stationary noise, in applications where a weak-sparse transform makes it possible to represent the signal of interest by a relatively small number of coefficients with significantly large amplitude. The resulting estimator called Extended-DATE (E-DATE) is robust in that it does not use prior knowledge about the signal or the noise except for the weak-spareness property. Compared to other methods in the literature, the E-DATE algorithm has the remarkable advantage of requiring only two parameters to specify. A straightforward block-based implementation of the E-DATE, called B-E-DATE, has first been introduced. This implementation entails an estimation delay, which diminishes as the frequency rate increases. This delay could be reduced by grouping frequency bins. Another solution to shorten this delay involves resorting to a sliding-window implementation called SW-E-DATE, but at the price of a higher computational cost. The B-E-DATE and SW-E-DATE have been benchmarked against various classical and recent noise power spectrum estimation methods in two situations: with and without noise reduction. The experimental results show that the E-DATE estimator generally provides the most accurate noise estimate, and that it outperforms other methods for speech denoising in the presence of various noise types and levels. For its good performance and low complexity, the B-E-DATE should be preferred in practice when frequency rates are high enough to induce acceptable or even negligible time-delay.

Although the present paper focused on noise reduction in speech enhancement systems, it must be emphasized that the E-DATE estimator is not restricted to speech signals and could find other applications in any scenario where noisy signals have a weakly-sparse representation. For many signals of interest, not limited to speech, such a weaklysparse representation can be provided by an appropriate wavelet transform. In this respect, the application of the E-DATE algorithm to audio separation could be considered in continuation of [START_REF] Aziz Sbai | Contribution of statistical tests to sparseness-based blind source separation[END_REF], [START_REF] Socheleau | Robust statistics based noise variance estimation: Application to wideband interception of noncooperative communications[END_REF], [START_REF] Aziz Sbai | Robust underdetermined blind audio source separation of sparse signals in the timefrequency domain[END_REF], [START_REF] Aïssa-El-Bey | Blind identification of sparse simo channels using maximum a posteriori approach[END_REF].

The E-DATE estimator fundamentally relies on the DATE estimator which, as emphasized in [START_REF] Pastor | Robust estimation of noise standard deviation in presence of signals with unknown distributions and occurrences[END_REF], can be regarded as an outlier detector. Consequently the E-DATE can also be used as an outlier detector in each frequency bin. This opens interesting perspectives in voice activity detection based on frequency analysis as well as in the detection and estimation of chirp signals in various types of noise.

Notation:

  In what follows, • is the usual euclidean norm in the space of all d -dimensional real vectors, I d stands for the d × d identity matrix, N (0, σ 2 0 I d ) designates the d -dimensional gaussian distribution with null mean and covariance matrix σ 2 I d and 1[X ∈ B ] stands for the indicator function of the event [U ∈ B ], where U is any random variable and B is any borel set of the real line: 1

( A0 ) 2 0

 A02 The observations Y 1 , Y 2 , . . . , Y n , . . . are mutually independent, Y n = ε n Λ n + X n where X n ∼ N (0, σ I d ) and ε n is Bernoulli distributed with values in {0, 1} for each n ∈ N.

Fig. 1 :

 1 Fig. 1: Spectrograms of clean and noisy speech signals from the NOIZEUS database. The noise source is car noise. No weighting function was used to calculate the STFT.

𝑌 1 , 1 𝑌 1 , 2 Fig. 2 :

 11122 Fig. 2: Principle of noise power spectrum estimation based on the DATE in colored stationary noise

Fig. 3 :

 3 Fig.3: Block E-DATE (B-E-DATE) combined with noise reduction (NR). A single noise power spectrum estimate is calculated every D non-overlapping frames and used to denoise each of these D frames.

Fig. 4 :

 4 Fig.4: Sliding-Window E-DATE (SW-E-DATE) combined with noise reduction. For the first D -1 frames, a surrogate method for noise power spectrum estimation is used in combination with noise reduction. Once D frames are available and upon reception of frame D + , ≥ 0, the SW-E-DATE algorithm provides the NR system with a new estimate of the noise power spectrum computed using the last D frames F +1 , . . . , F +D for denoising of the current frame.

  g = 3.093 -1.029 LLR0.603 PESQ -0.009 WSS C bak = 1.634 + 0.478 PESQ -0.00 WSS + 0.063 segSNR C vol = 1.594 + 0.805 PESQ -0.512 LRR -0.007.WSS

Fig. 5 :

 5 Fig.5: Noise estimation quality comparison of several noise power spectrum estimators at different SNR levels and with different kind of noise.

Fig. 7 :

 7 Fig.7: Speech quality evaluation after speech denoising (C ovl composite criterion).

Fig. 8 :

 8 Fig. 8: Speech quality evaluation after speech denoising (C bak composite criterion).

TABLE I :

 I Number of parameters required by different noise power spectrum estimation algorithms

	Method	MARTIN[7]	IMCRA[5]	MCRA2[9]	MMSE1[10]	MMSE2[11]	ML-ME[12]	E-DATE
	Parameters number	7	10	7	3	5	3	2

TABLE IV :

 IV Computational cost of MMSE2 per new frame and per frequency bin For instance, additional experimental results with airport, babble, station, car and train noises suggest that D can be chosen in the range [50, 80] without really affecting C ovl for SNR > 0dB.

	Addition	Multiplication	Division	Exponent
	12	10	2	1
	noising performance.

TABLE II :

 II Computational cost of B-E-DATE per group of D frames and per frequency bin

		Addition	Multiplication	Division	Square root
	Norm	D	2D	0	D
	Sorting	D log D	0	0	0
	Search n * (worst case)	D(D -1)/2	D	D	0
	Total	D log D + (D + 1)/2	3D	D	D

TABLE III :

 III Computational cost of SW-E-DATE per new frame and per frequency bin

		Addition	Multiplication	Division	Square root
	Norm	1	2	0	1
	Sorting	log D	0	0	0
	Search n * (worst case)	D(D -1)/2	D	D	0
	Total	1 + log D + D(D -1)/2	D + 2	D	1


Ala Aboudib
email: ala.aboudib@telecom-bretagne.eu

Vincent Gripon
email: vincent.gripon@telecom-bretagne.eu

Baptiste Tessiau
email: baptiste.tessiau@ens-rennes.fr

Implementing Relational-Algebraic Operators for Improving Cognitive Abilities in Networks of Neural Cliques

Keywords: Cognitive modeling, artificial neural networks, relational algebra, associative memory, simulated annealing

Associative memories are devices capable of retrieving previously stored messages from parts of their content. They are used in a variety of applications including CPU caches, routers, intrusion detection systems, etc. They are also considered a good model for human memory, motivating the use of neural-based techniques. When it comes to cognition, it is important to provide such devices with the ability to perform complex requests, such as union, intersection, difference, projection and selection. In this paper, we extend a recently introduced associative memory model to perform relational algebra operations. We introduce new algorithms and discuss their performance which provides an insight on how the brain performs some high-level information processing tasks.

I. INTRODUCTION

Associative memories are special types of memories that are capable of high-speed content-based mapping between input queries and outputs. This type of storage differs from classical index-addressable memories in that no explicit address is needed to search for stored data. In order to efficiently provide this associative functionality to memories, different methods of content structuring are required. Artificial neural networks (ANNs) are known to be such an adapted medium to implementing associative memorization. Their design is inspired from neo-cortical neural mechanisms in mammalian brains, believed to be knowledge-associative biological neural networks [START_REF] Anderson | Human associative memory[END_REF].

Many ANN models, differing in topology and functionality, were proposed to act as associative memories. Famous examples include the Perceptron [START_REF] Rosenblatt | The perceptron: a probabilistic model for information storage and organization in the brain[END_REF], self-organizing maps [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF], Hopfield networks [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] and Boltzmann machines [START_REF] Ackley | A learning algorithm for boltzmann machines[END_REF]. A new model was proposed recently by Gripon and Berrou in [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] and generalized by Aliabadi et al. in [7]. This model relies on sparse coded patterns stored as cliques and is based on Hebian learning [START_REF] Hebb | The Organization of Behavior[END_REF]. This sparse coding approach resembles the ones introduced by willshaw [START_REF] Willshaw | Nonholographic associative memory[END_REF] and Palm [START_REF] Palm | On associative memory[END_REF], with an added explicit mapping between stored messages and their representation in the network.

Typically, an associative memory is able to retrieve a previously stored piece of information given partial content, a sort of erasure-retrieval property. However, human ability to handle more complex queries suggests that their representation of information should be able to perform other operations. Relational algebra gives a formal framework to introduce complex operations on tuples of stored messages, including union, intersection, difference, projection and selection. In this paper, we aim at extending the model introduced in [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] to process these operations.

Most of these operations are common as far as human cognition is concerned. For example, selection aims at retrieving the list of all messages that match a given probe (e.g., listing all city names you know that start with the letter 'b'). The union operation can be used for merging date while intersection and difference operation are useful for comparing contents of several memories or memory regions. Our main motivation is to show that existing neural-network-based architectures for associative memories are easily adapted in order to handle these complex queries. There have been several works addressing the relational problem and its biological plausibility such as [START_REF] Blockeel | Using neural networks for relational learning[END_REF] and [START_REF] Hummel | A symbolic-connectionist theory of relational inference and generalization[END_REF], which were more focused on inference and relational learning.

The rest of this paper is organized as follows: in Section II, we introduce the associative memory model extended in this paper. Sections III, IV, V and VI describe how to perform resp. union, intersection, difference and projection using these models. In Section VII, we explain how to handle the more complex selection operator. For this operator, we introduce a novel algorithm using simmulated annealing. Simulation results are provided in Section VIII when independently identically uniformly distributed messages are considered. Section IX is a conclusion.

II. THE MEMORY MODEL AND RELATIONAL ALGEBRA

A. The associative memory model

First, we introduce the associative memory model proposed by Gripon and Berrou in [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] and then extended in [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF]. Let us consider the finite alphabet A made of integers between 1 and . we define a blank character ⊥ ∈ A and A = A ∪ {⊥}. We are interested in storing words made of χ characters over A. We call such a word a message m = m 1 m 2 . . . m χ . Blank characters represent the absence of a character at a given position, such that depending on their number, messages can be regarded as sparse vectors.

An associative memory is a device capable of storing messages and then retrieving them given partial knowledge about some of their nonblank characters. To implement this device, the authors of [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] propose to use a symmetric, binary, χ-partite neural network composed of n = χ • vertices that we shall refer to as units. The authors of [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] explain that units should be considered analogous to cortical microcolumns believed to be the computational building-blocks of the cerebral neo-cortex [13] [14]. This network can be split into χ clusters, each containing the same number of units.

We denote by [n] the set of integers beween 1 and n, let us then index each cluster of the network by an integer in [χ]. We also index units of a given cluster by integers in [ ]. As a result, each unit in the network is uniquely addressed giving a couple (i, j) where i is the index of a cluster and j the index of the unit within this cluster.

We define a function f that maps each message m into a set of couples (i, j) as follows:

f : m → {(i, j) i ∈ [χ], j ∈ [ ] and m i = j} (1)
where m i refers to the i-th character of m.

The network topology is entirely captured by an adjacency matrix W of size χ • such that:

w (i,j)(i ,j ) = 1 if (i, j) and (i , j ) are connected 0 otherwise (2)
To store a message m = m 1 m 2 . . . m χ in such a network, all unit pairs in f (m) are connected pairwise according to 2 forming a clique in the underlying graph. It is worth noticing that cliques corresponding to different messages could sometimes overlap and share connections. Given a set of messages M, we denote by W (M) the adjacency matrix obtained after storage of all messages in M.

A subset of f (m) is called a partial input associated with the message m. The task of an associative memory is then: given a partial input of m ∈ M, retrieve m using W (M).

To retrieve a message from a partial input, an iterative algorithm is performed. Retrieval algorithms and techniques have been discussed in detail in [START_REF] Aboudib | A study of retrieval algorithms of sparse messages in networks of neural cliques[END_REF].

It has been shown in [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] that performance of this structure as an associative memory mainly depends on a density parameter defined as the ratio of the number of connections in the network to the total number of possible connections. Under the hypothesis that messages contain exactly c nonblank characters uniformly distributed over A, density can be approximated by the following equation:

d = 1 -1 - c(c -1) χ(χ -1)l 2 M (3)
In this paper, we aim at extending the functionality of these associative memories, which we shall call Clustered Clique Networks (CCNs), to cover more general problems defined in relational algebra.

B. Connections to relational algebra

In relational algebra, operators are defined on relations (sets of tuples). A set of attributes is associated with each relation. Then, each tuple is defined as a set of instances of these attributes. A CCN in this respect can be viewed as a relation. Each cluster represents an attribute and units within each cluster are instances of that attribute (attribute values). A clique connecting units is equivalent to a tuple. In the following sections, we are going to propose algorithms and methods for implementing some relational operators on relations defined in the form of CCNs. III. UNION Defined in terms of the set theory, the union of a collection of sets S 1 , S 2 , ..., S n is a set S ∪ containing all distinct elements in this collection. It can be described as follows:

S ∪ = S 1 ∪ S 2 ∪ ... ∪ S k (4) = {x|x ∈ S 1 ∨ x ∈ S 2 ∨ ... ∨ x ∈ S k } (5)
where ∪ is the set union operator and ∨ is the Boolean OR function.

In the context of memory storage, union is used to combine the contents of several memories or memory partitions into a single one while avoiding the redundancy resulting from the same data-word being stored multiple times. An example of this is merging the contents of two folders on a computer. This task is straightforward when using a classical indexed memory since messages do not overlap.

Suppose that we have two CCNs with the same dimensions W (M 1 ) and W (M 2 ) that we wish to merge in a single network W ∪ of the same dimensions. We define this operation as follows:

w ∪ (i,j)(i ,j ) =    1 if w(M 1 ) (i,j)(i ,j ) = 1∨ w(M 2 ) (i,j)(i ,j ) = 1 0 otherwise (6)
Following intuitively from this definition is the fact that upon applying the union operation, information is conserved in the new network. That is, if a clique exists in either W (M 1 ) or W (M 2 ), it would also exists in W ∪ . Hence we can rewrite

W ∪ as W ∪ (M) with M = M 1 ∪ M 2 .
The problem here is that combining memories in this fashion can cause a significant growth in the density of W ∪ (M), leading possibly to dramatically low performance in terms of retrieval error rates of stored messages. More formally, if d 1 is the density of W (M 1 ) and d 2 is the density of W (M 2 ), then the density of W ∪ (M) denoted d ∪ is given by:

d ∪ = 1 -(1 -d 1 )(1 -d 2 ) (7)
and thus d ∪ is greater than both d 1 and d 2 . This means that while retrieval error rates may be optimal for W (M 1 ) and W (M 2 ), the network resulting from their union might suffer from some degeneration in performance because of the increased density. Therefore, union should be done only if the resulting network performance is acceptable. The relationship between retrieval error rates and density are presented in details in [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] and [START_REF] Aboudib | A study of retrieval algorithms of sparse messages in networks of neural cliques[END_REF]. According to [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF] all possible connections in the networks should be tested during the union operation. So, given χ(χ-1) 2 2 possible connections [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] in each network, the average-case complexity of such process is Θ(χ 2 2 ).

The phenomenon of degenerated memorization efficiency caused by the increased density is not uncommon in the brain. For example, learning and recalling a new word in a foreign language is typically not a difficult task. However, trying to learn a dozen of new words at the same time might turn out to be much more challenging comprising many memorization errors and confusions and even mixing syllables of different words. New words can be better learned by training and experience, which is partially due to the association of these words with other memories. We suggest that this process, in terms of CCNs is equivalent to adding more clusters to the network such that more units can be added to existing cliques, which lowers its density and increases performance.

IV. INTERSECTION

The intersection among several sets S 1 , S 2 , ..., S k is a set S ∩ containing only those elements that S 1 , S 2 , ..., S k have in common:

S ∩ = S 1 ∩ S 2 ∩ ... ∩ S k (8) = {x|x ∈ S 1 ∧ x ∈ S 2 ∧ ... ∧ x ∈ S k } ( 9
)
where ∩ is the set intersection operator and ∧ is the Boolean AND function. When applied to relations, intersection serves in extracting common tuples between two or more tables having the same number and types of attributes. This is done easily when such a database is stored in an indexed memory. One way to implement intersection between two CCNs is by keeping only those connections that happen to exist in the exact same place in both networks. So, given two CCNs W (M 1 ) and W (M 2 ) of the same type and dimensions, we can define their intersection W ∩ as follows:

w ∩ (i,j)(i ,j ) =    1 if w(M 1 ) (i,j)(i ,j ) = 1∧ w(M 2 ) (i,j)(i ,j ) = 1 0 otherwise (10)
The average-case complexity of this operation is given by Θ(χ 2 2 ) for the same reason as in the union. The density d ∩ of W ∩ as a function of d 1 and d 2 (the densities of W (M 1 ) and W (M 2 ) respectively) is given by:

d ∩ = d 1 .d 2 (11)
We notice from ( 11) that d ∩ is always lower than d 1 and d 2 given that densities have their values in the interval [0,1]. Thus, we guarantee that no density explosion occurs as in union.

Some problems might still occur when performing this operation, as depicted in Figure 1 where four simple identical networks are considered each with a total number of units n = 12 grouped in χ = 3 clusters with a message size of c = χ. Network W (M 1 ) shown in Figure 1a contains two cliques af h and beg and network W (M 2 ) in Figure 1b also contains two cliques af r and beg. We notice that only the clique beg is common between these two networks. By applying the intersection operation as in [START_REF] Palm | On associative memory[END_REF] we obtain W ∩ shown in Figure 1c, which contains the common clique beg as expected but also contains the edge af , which is an undesirable result, because af does not represent a complete message (tuple). We shall call af a residual edge. As a consequence, W ∩ could possibly contain more connections than a an ideal intersection network W ∩ (M) with M = M 1 ∩ M 2 . This increase in density due to residual edges is expected to deteriorate performance [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF].

V. DIFFERENCE The difference of two sets S 1 and S 2 , which can also be called that relative complement of S 2 with respect to S 1 , is a set S \ that contains only those elements of S 1 that are not in S 2 :

S \ = S 1 \ S 2 = {x|x ∈ S 1 and x / ∈ S 2 } ( 12
)
where \ is the set difference operator which is not commutative so that S 1 \ S 2 = S 2 \ S 1 . So, the difference between two database tables is the set of tuples in the first one that do not exist in the other, which is also an easy-to-implement operation in classical memory systems. A simple method of implementing difference between two CNNs W (M 1 ) and W (M 2 ) is by instantiating a new memory W \ containing only connections in W (M 1 ) that do not exist in W (M 2 ). That is:

w \ (i,j)(i ,j ) =    1 if w(M 1 ) (i,j)(i ,j ) = 1∧ w(M 2 ) (i,j)(i ,j ) = 0 0 otherwise (13)
The average-case complexity of this operation is also given by Θ(χ 2 2 ). Intuitively, the density d \ of the new network W \ is always less than or equal to the density d 1 of W (M 1 ). So, if d 1 is well controlled, we would have no problems with the performance of W \ . The density d \ is given by the following relationship:

d \ = d 1 .(1 -d 2 ) (14)
As in the case of intersection, the difference operation defined in [START_REF] Jones | Microcolumns in the cerebral cortex[END_REF] processes data down on the level of individual connections not on the level of cliques. This causes the problem depicted in Figure 1d. In this figure, W \ is the network resulting from applying W (M 1 ) \ W (M 2 ). Ideally, we wish that W \ = W \ (M) with M = M 1 \ M 2 ,i.e., a network that contains only the clique af h. However, according to (13), we only get two edges ah and f h because af is a common edge between W (M 1 ) and W (M 2 ) and thus eliminated by the difference operation. This represents a loss of information since the network W \ no more stores the message corresponding to the clique af h. We shall call this undesirable effect erosion. Actually, no method is yet available for getting an ideal intersection or an ideal difference between CCNs. We consider the methods proposed in this paper as approximations to the real operations.

VI. PROJECTION

In relational algebra, projection is defined as a unary operator Π applied to a tuple R in order to produce a new tuple R Π consisting of k attributes {r 1 , r 2 , ..., r k }, which is a subset of the attributes originally contained in R. This can be written as follows:

R Π = Π r1,r2,...,r k (R) (15)
A network W Π is said to be a projection of a network W on a given set of attributes, if it contains only a subset of the attributes of W ,i.e., W Π contains only a subset of the clusters of W .

Clearly, this operation is very easy to implement in CNNs with a constant time complexity. Moreover, the density d Π of the resulting network is equal to the density of the original network given that connections are uniformly distributed within the latter network:

d Π = d (16)

VII. SELECTION A. Problem definition

In relational algebra, selection or restriction is a unary operator applied to a relation R 1 and returns another relation R 2 . The latter relation contains all tubles in R 1 whose attribute values satisfy a propositional formula ϕ. This can be transcribed as follows:

R 2 = σ ϕ (R 1 ) (17)
We argue in this paper that a mechanism similar to selection might be used by the brain for thinking and memorization. A typical such request would be, for instance, to name all scientific authors one knows whose names start with an 'a'. We aim at using the model proposed in [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] as a substrate for this selection process.

The selection algorithm we shall present here runs continuously, giving multiple (possibly redundant) answers one after another. This appears to us being behaviorally similar to what humans could produce facing a similar query. The process that makes us avoid redundancy is called short-term/working memory. Once its buffer is full or overcrowded, repetitions of the same word/name can occur.

The requirements that our selection algorithm is meant to meet are the following:

1) Determine the sub-graph G of potentially interesting units. 2) In sub-graph G, find all cliques of size c.

Finding a maximum clique in a graph (or equivalently) a minimum cover is a know NP-complete problem. Many algorithms and heuristics were proposed to give approximate solutions to this problem in medium-sized graphs. Examples of these algorithms are [START_REF] Geng | A simple simulated annealing algorithm for the maximum clique problem[END_REF] and [START_REF] Xu | An efficient simulated annealing algorithm for the minimum vertex cover problem[END_REF] that make use of the simulated annealing principle introduced in [START_REF] Brooks | Optimization using simulated annealing[END_REF]] and [START_REF] Černỳ | Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm[END_REF],

which is a probabilistic meta-heuristic method for locating the global optimum of a given function. Another known method widely used in applications such as computational chemistry is the BronKerbosch algorithm [START_REF] Bron | Algorithm 457: finding all cliques of an undirected graph[END_REF], which can efficiently find maximal cliques in an undirected graph.

We propose to use an adaptation of the simulated annealing algorithm proposed by Geng et al. in [START_REF] Geng | A simple simulated annealing algorithm for the maximum clique problem[END_REF]. We also use their same objective function to evaluate our solutions.

B. The proposed selection algorithm

Suppose we have a CCN denoted by W and a partial input message m containing q c known nonblank attributes. The selection operator consists in finding all stored cliques made of a set of units containing f (m). In order to perform this search efficiently, it is sufficient to restrict the search to the subgraph made of only the units connected to units in f (m). We shall call this subgraph G m . We denote its adjacency matrix by A m .

For simplicity of representation, we refer to each unit of G m by an integer index k or s where k, s ∈ {0, 1, ..., n -1}, n being the number of units in G m . Our objective now is to find all the cliques in G m that have a size (number of units) of c = c -q.

We consider the following optimization problem: We define ρ as a permutation of units in G m (ρ is an array of size n containing all unit indexes of G m as its elements). For a given ordering of elements in ρ, we consider the first c elements of ρ as indexes of units in G m acting as a potential solution (a clique). So, by permuting ρ's contents we can get a different candidate solution. The objective function used to evaluate these solutions is given by:

F (G m , ρ) = c -2 k=0 c -1 s=k+1 (1 -a ρ[k],ρ[s] ) (18)
where a ρ[k],ρ[s] is an element of the adjacency matrix W . As F (G, ρ) = 0 when a clique is found, the goal is to find permutations that minimize this function.

The algorithm is applied to G m as follows:

Step A single run of this selection algorithm is meant to find one clique ,i.e., one answer. So, in order to output more answers this algorithm should be repeated many times.

VIII. RESULTS

A. Selection

In order to test the proposed selection algorithm, we used a CCN with n = 3240 units grouped in χ = 15 clusters. 15000 randomly generated messages of c = 10 characters were stored giving a network density d of about 0.13. We used the same initial and end temperatures as in [START_REF] Geng | A simple simulated annealing algorithm for the maximum clique problem[END_REF] for the simulated annealing algorithm; 100 and 0.001, respectively. We set the cooling coefficient to 0.996. To construct a selection query message, one already stored message is randomly selected of which 9 characters were erased (set to ⊥) giving an input query message with only one known non-blank character. As described in Section VII, the algorithm we propose might give redundant outputs when executed several times. Figure 2 shows how fast unique outputs were found as a function of the number of executions of the selection algorithm. Curves in the figure are averaged over 100 identical experiments.

An interesting property of the resulting curve is that most unique answers are obtained during earlier executions of the algorithm. So, referring back to Figure 2, about 45 unique answers out of 47 possible ones are found by the 250 th execution, 40 answers by the 125 th execution and about 30 answers by the 60 th . In other words, about 89% of answers were obtained when only 50% of total executions were achieved and 67% of answers after 24% of total executions. This is a natural result for a selection by replacement experiment where all answers have an equal chance of being chosen at each execution.

We suggest that such result bears some qualitative resemblance to the way human beings memorize lists of mental objects where it is common that the last few items turn out to be more difficult and time consuming to recall because of the distraction caused by redundant answers coming to mind and other phenomena.

B. Intersection and difference

A comparison among average complexities of some operators when applied to CCNs and two other known data structures (ordered lists and binary search trees) storing sparse messages of the form m = m 1 m 2 . . . m χ is provided in Table I. An interesting observation is the fact that the order of complexity of operators using a CCN is close to that of a binary search tree given that setting χ is preferable in practice for a higher network capacity [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF]. However, the complexity of union, intersection and difference of ordered lists is lower by a factor of χ than that of CCNs while the complexity of insertion is 2 /χ higher. IX. CONCLUSION AND FUTURE WORK In this paper, we have introduced some methods for applying certain algebraic-relational operators on a new class of neural-network-based associative memories we call CCNs. We argued that the process of recalling a list of items (which can also be mapped to more general memorization tasks) in the brain can be behaviorally assimilated to the selection operation known to relational algebra. We proposed an algorithm for implementing this process using the principle of simulated annealing. Then, we showed that the results we got have some resemblance to what might be obtained by a human subject in terms of redundancy.

We have also demonstrated that CCNs can be used as classic data-structures by approximating operators such as union, intersection, difference and projection. We saw that the implementation of union and its related density explosion problem raised the question as to how the brain organizes information with high correlation or high density. Two possible mechanisms the brain might be using are forgetting rarely used "data" (by the decay of synaptic weights) and tagging pieces of correlated data with different contextual information. Similar mechanisms might be integrated in CCNs by allowing connections to have real values with a decay parameter and by providing contextual tagging in the form of CCNs existing on a separated level of a hierarchy of networks. Another possible solution is to design networks with dynamic sizes to prevent exceeding a maximum allowed density.

Figure 1 .

 1 Figure 1. Intersection and Difference between two CCNs. W ∩ shown in 1c is the network resulting from the intersection of W (M 1 ) and W (M 2 ). W \ shown in 1d is the difference of W (M 1 ) from W (M 2 ).

  of the selection algorithm Number of unique outputs found Accumulated number of unique outputs found Target number of unique outputs

Figure 2 .

 2 Figure 2. Average accumulated number of unique outputs as a function of the number of iterations in a network of n = 3240, χ = 15, c = 10, 15000 stored messages. The input message used has only one non-blank character.

  1: Parameter initialization. Initial temperature T 1 , end temperature T 2 , current temperature t = T 1 and cooling coefficient α. ], ρ[1], ..., ρ[c -1]} than ρ[w] has, then the index w is rejected and we go back to step 3. but if w has already been rejected for more than 8n times consecutively as defined in [16], then ρ[w] and ρ[u] are swapped and a new permutation ρ is obtained. Step 4: Compute F (G m , ρ ) and terminate if it evaluates to zero. Step 5: Compute∆F = F (G m , ρ ) -F (G m , ρ).If ∆F 0, then accept the new permutation ρ by setting ρ = ρ . Otherwise accept ρ with probability p = e

	Set the initial permutation as ρ[k] = k, k ∈ {0, 1, ..., n -1}. Random setting of permutation is also possible. Step 2: Compute F (G m , ρ) and terminate if it evaluates to zero. Step 3: Randomly choose two integer indexes u and w of ρ such that u ∈ {0, 1, ..., c -1} and w ∈ {c , c + 1, ..., n -1}. Condition 1: If ρ[w] has more or the same number of connec-tions with {ρ[0], ρ[1], ..., ρ[c -1]} than ρ[u] has, then ρ[w] and ρ[u] are swapped and thus a new permutation ρ is obtained. Condition 2: If ρ[u] has more connections with {ρ[0-∆F t . Step 6: Update current temperature as t = αt. If t < T 2 terminate the algorithm; otherwise, go back to step 3.

TABLE I .

 I COMPLEXITY OF SOME RELATIONAL OPERATORS IN SEVERAL TYPES OF DATA STRUCTURES.

		CCN	ordered list	binary search tree
	Insertion(Storing)	Θ(χ 2 )	Θ(χ 2 )	Θ(χ log( ))
	Union	Θ(χ 2 2 )	Θ(χ 2 )	Θ(χ 2 log( ))
	Intersection	Θ(χ 2 2 )	Θ(χ 2 )	Θ(χ 2 log( ))
	Difference	Θ(χ 2 2 )	Θ(χ 2 )	Θ(χ 2 log( ))

ACKNOWLEDGMENT This work was supported by the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) / ERC grant agreement n 290901.