Solving Time Domain Audio Inverse Problems using Nonnegative Tensor Factorization C ¸agdas ¸Bilen, Alexey Ozerov, and Patrick Pérez Abstract-Nonnegative matrix and tensor factorization (NMF and NTF) are important tools for modeling nonnegative data, which gained increasing popularity in various fields, a significant one of which is audio processing. However there are still many problems in audio processing, for which the NMF (or NTF) model has not been successfully utilized. In this work we propose a new algorithm based on NMF (and NTF) in the short-time Fourier domain for solving a large class of audio inverse problems with missing or corrupted time domain samples. The proposed approach overcomes the difficulty of employing a model in the frequency domain to recover time domain samples with the help of probabilistic modeling. Its performance is demonstrated for the applications such as audio declipping and declicking (never solved with NMF/NTF modeling prior to this work), joint audio declipping/declicking and source separation (never solved with NMF/NTF modeling or any other method prior to this work), compressive sampling recovery and informed source separation (a low complexity encoding scheme that is possible with the proposed approach and has never been proposed prior to this work).

I. INTRODUCTION

Nonnegative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects with nonnegative matrix factorization[END_REF] and nonnegative tensor factorization (NTF) [START_REF] Cichocki | Nonnegative matrix and tensor factorization[END_REF] decompositions have recently found great success in applications to audio modeling, notably for source separation [START_REF] Virtanen | Monaural Sound Source Separation by Nonnegative Matrix Factorization With Temporal Continuity and Sparseness Criteria[END_REF]- [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], compression [START_REF] Nikunen | Object-based audio coding using nonnegative matrix factorization for the spectrogram representation[END_REF], [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF], [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF] and audio inpainting [START_REF] Roux | Computational auditory induction as a missing-data model-fitting problem with Bregman divergence[END_REF]- [START_REF] Simsekli | Score guided audio restoration via generalised coupled tensor factorisation[END_REF]. Within all these applications the power-spectrograms of singlechannel or multichannel audio signals (usually powers of their short-time Fourier transforms (STFT)) are decomposed using NMF or NTF models.

However, these methods address quite poorly the situations when some chunks or samples of audio signals are missing in time domain, as for example in the situations of audio declipping or declicking, as described in a general audio inpainting paper [START_REF] Adler | Audio inpainting[END_REF]. Indeed, the NMF/NTF-based audio inpainting methods [START_REF] Roux | Computational auditory induction as a missing-data model-fitting problem with Bregman divergence[END_REF]- [START_REF] Simsekli | Score guided audio restoration via generalised coupled tensor factorisation[END_REF] assume that the audio data is missing directly in the corresponding time-frequency domain, usually the STFT domain. This is in fact the most convenient situation since the modeling itself is formulated in the STFT domain, and thus it becomes quite easy to take properly into account the missing values. In the case of audio with missing samples in the time domain, one can convert the missing information into an STFT domain formulation by simply assuming that all the STFT frames corresponding to missing time samples are missing in entirety. However, this C ¸. Bilen is with University of Strasbourg, France, e-mail: cagdas@unistra.fr. A. Ozerov and P. Pérez are with Technicolor, France, e-mail: alexey.ozerov@technicolor.com; patrick.perez@technicolor.com.

This work was partially supported by ANR JCJC program MAD (ANR-14-CE27-0002).

will often lead to the loss of a huge amount of available information. Indeed, for example in the case of a clipped audio every STFT frame may be clipped, thus this naive solution would lead to considering the whole signal to be missing, even though there is perhaps only 20 % of the signal that is clipped in the time domain. Another problem of NMF/NTFbased audio inpainting methods [START_REF] Roux | Computational auditory induction as a missing-data model-fitting problem with Bregman divergence[END_REF]- [START_REF] Simsekli | Score guided audio restoration via generalised coupled tensor factorisation[END_REF] which consider fully-missing STFT coefficients is that NMF/NTF models are phase-invariant and thus they only allow estimating the magnitudes of the missing coefficients. As a result, the phase information, which is very important for audio perceptual quality, still needs to be reconstructed somehow. A popular approach by Griffin and Lim [START_REF] Griffin | Signal reconstruction from short-time Fourier transform magnitude[END_REF] is usually used for the phase reconstruction, but it performs quite poorly in many situations. As an alternative, a so-called high resolution NMF (HR-NMF) approach was proposed [START_REF] Badeau | Multichannel high-resolution NMF for modeling convolutive mixtures of non-stationary signals in the timefrequency domain[END_REF], [START_REF] Badeau | Multiplicative updates for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF]. This approach extends the NMF to model temporal dependencies between time-frequency bins, which yields better phase estimates. However, for the moment this approach is quite computationally expensive and it is limited to harmonic sounds. At the same time, when some samples are missing in the time domain and one manages to estimate properly the phase-invariant NMF model and the missing samples from these observations, the resulting phase estimates should be better than those obtained via Griffin and Lim's approach [START_REF] Griffin | Signal reconstruction from short-time Fourier transform magnitude[END_REF], since missing samples in time domain does not mean completely discarding the phase information in the STFT domain.

In this work, we propose a new approach allowing the estimation of lost time domain audio samples of audio sources and/or their mixture via applying a low-rank NMF/NTF model to energy distribution of the signals in the time frequency domain. The proposed method uses Itakura Saito (IS) divergence [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] for measuring how well the given NMF/NTF model parameters are able to estimate the signal variances while using all the information available from all of the known time domain samples from the sources and/or the mixture. The model parameters are estimated using a generalized expectation-maximization (GEM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] and Wiener filtering [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] is used to recover the unknown signals.

Unlike some other approaches that directly apply NMF/NTF model on the STFT coefficient magnitudes or powers, the proposed approach is formulated as a probabilistic Gaussian model on the complex-valued STFT coefficients which enables to estimate the NMF/NTF model in a maximum likelihood (ML) sense directly from the time domain observations, thus avoiding sub-optimally converting this missing information into the STFT domain. Furthermore, thanks to the flexibility of the NMF/NTF representations, the proposed framework can take into account mixtures of several sources, where both the sources and the mixtures can be partially or fully-missing in time domain. Last but not least, when the observed signals are not only lost but also corrupted, such as by noise or quantization, these corruptions can also be taken into account in the proposed approach. Within this general formulation the proposed framework is not limited to audio inpainting, but also becomes useful for different new applications related to audio compression, enhancement and source separation. This work builds on several previous conference/workshop publications [START_REF] Bilen | Audio declipping via nonnegative matrix factorization[END_REF]- [START_REF] Ozerov | Multichannel audio declipping[END_REF] by the authors. Particular instances of the proposed approach for some specific applications have been presented in [START_REF] Bilen | Audio declipping via nonnegative matrix factorization[END_REF]- [START_REF]Compressive sampling-based informed source separation[END_REF] and summarized in [START_REF]Audio inpainting, source separation, audio compression. all with a unified framework based on ntf model[END_REF]. In this paper, we provide a generalized formulation that is applicable but not limited to these specific applications and present more comprehensive and extended experimental results, notably we provide new experiments on compressive sampling recovery. Audio source separation approach presented in [START_REF]Automatic allocation of NTF components for user-guided audio source separation[END_REF] falls within this general formulation as well, but it is not considered here for conciseness. Finally, while we here formulate the framework in the case of single-channel audio mixtures, its extension to the multichannel case is straightforward, which has been demonstrated in [START_REF] Ozerov | Multichannel audio declipping[END_REF] for the declipping application.

The rest of this paper is organized as follows. Some applications for which the proposed approach can be useful, are described in Section II. The problem is formally defined in Section III and the proposed algorithm to solve this problem is described in Section IV. Experiment results for various applications are given in Section V and lastly final remarks and conclusions are presented in Section VI.

II. SOME RELEVANT APPLICATIONS

The proposed approach is useful in (but not limited to) several applications, both new and existing, which are summarized as follows.

a) Time Domain Audio Inpainting: The problem of recovering audio samples that are lost or corrupted is often called audio inpainting. We use the term "time domain audio inpainting" to refer to the problems with missing or corrupted time domain audio samples as opposed to the audio inpainting problems with missing samples in the STFT domain, for which NMF/NTF models are already being used prominently. Two specific instances of the time domain audio inpainting problem are called audio declipping and audio declicking, in which one recovers the time domain audio samples that are lost due to clipping and clicking effects caused by audio recording and compression processes. The declipping problem in particular provides additional challenges with respect to the general time domain audio inpainting problem, because it often includes additional constraints for the time domain signal to be estimated. In the recent years, the models based on sparse, cosparse or group-sparse representations in certain dictionaries are shown to be performing best to solve these problems [START_REF] Adler | Audio inpainting[END_REF], [START_REF] Kitić | Consistent iterative hard thresholding for signal declipping[END_REF]- [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF]. Despite the success of modeling audio signals with low rank NMF representations, the time domain audio inpainting problem, especially with the additional constraints in audio declipping, is not trivial to solve with an NMF/NTF model in the time-frequency domain. This is possibly the main reason why these models have not been utilized in time domain audio inpainting problems successfully. The proposed approach can overcome this limitation and provides a new perspective on time domain audio signal recovery with equivalent or better performance than the state of the art.

b) Joint audio inpainting and source separation: The audio source separation is a well known problem for which the NMF/NTF modeling in the time-frequency domain is shown to be quite successful. However in all source separation problems, the audio mixture is assumed to be known perfectly whereas in practice the mixture can also have missing or corrupted samples in time domain. Such a source separation problem with an incomplete and/or corrupted mixture is in fact a new problem that we introduce and address in this work, which, to our knowledge, can not be properly solved by any of the existing source separation approaches in the literature, except by sequentially performing audio inpainting followed by source separation on the reconstructed mixture. The latter sequential approach can be quite suboptimal since neither of these two tasks use all of the information efficiently. The problem of jointly performing the two tasks is for the first time addressed by our proposed approach, which can recover the signal in a way that is more consistent with the multiple source nature of the corrupted mixture while simultaneously estimating the individual sources. c) Compressive Sampling Recovery: Compressive sampling [START_REF] Candès | An introduction to compressive sampling[END_REF] is the theory and application of randomly subsampling a signal that is known to be compressible (e.g., with sparse or low rank representations) in an incoherent domain and making sense of the random samples by using the prior information of compressibility. As our algorithm is well fitted for time domain audio inverse problems, the reconstruction of the randomly sampled audio signals is another field of application for which it can be useful. Looking from the compressive sampling perspective, the compressible characteristics of the audio signals exploited by our algorithm are two fold: i) the significant reduction of the probability space of the possible solutions given the known samples, through the maximum likelihood estimate ii) the further reduction of the possible solutions through the low rank modeling of the NMF/NTF representation in the STFT domain. This application can in fact be seen as another instance of audio inpainting, however we have investigated it separately as the random subsampling changes the characteristic of the problem with respect to the other more typical audio inpainting problems such as audio declipping.

d) Compressive sampling-based informed source separation: Informed source separation (ISS) [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF] is a variant of source separation that is in fact a source compression problem assuming that the mixture is known. The ISS problem can be defined as the problem of encoding multiple audio sources to create a bitstream (also called a side-information) so that the audio from the sources can be recovered given the bit stream and the mixture of the sources. The main difference of ISS from joint compression of multiple audio signals is the assumption that the mixture is available at both encoding and decoding stages. Several ISS methods were proposed [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF],

MIXTURE

Fig. 1: The general framework of the proposed algorithm illustrating recovery of a mixture signal and its sources from a subset of the quantized samples of the sources and/or the mixture.

[31] including those based on the NTF modeling [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF].

In all these approaches the encoding stage is usually more complex and computationally expensive than the decoding stage. The framework proposed in this work can be used to realize a new variant of ISS that we call compressive samplingbased ISS (CS-ISS), where the computational complexity is moved from the encoder to the decoder side by reducing the encoder to simply subsampling the sources in a random and independent fashion and quantizing the samples. The proposed algorithm can then be used to recover the sources at the decoder side given the encoded samples and the mixture, similar to the case of compressive sampling recovery. This new approach is inline with both the compressive sampling [START_REF] Candès | An introduction to compressive sampling[END_REF] paradigm, since the sampling is random and in a sufficiently incoherent domain, and with the distributed source/video coding [START_REF] Xiong | Distributed source coding for sensor networks[END_REF], [START_REF] Girod | Distributed video coding[END_REF], since the posterior source dependencies (the sources are highly correlated a posteriori given the mixture) and the source structure are exploited only at the decoding stage, thus allowing the complexity shift. The CS-ISS also allows independent structures between the encoder and the decoder, i.e., the decoder algorithm can be modified without the need to change the encoder and the encoded bitstream.

III. PROBLEM DEFINITION

Let us consider a single-channel1 mixture that is composed of J sources, among which each of the sources and/or the mixture might be fully, or partially observed and/or corrupted with noise (e.g., quantization noise). For a mixture of length T , the mixture samples, x ′′ t ,2 are measured at a subset Ξ ′′ ⊂ 1, T of the entire time domain. Hence, the measured mixture samples can be represented in terms of in general unknown source samples, s ′′ jt , j ∈ 1, J , as

x ′′ t = J j=1 s ′′ jt + a ′′ t , ∀t ∈ Ξ ′′ , (1)
where a ′′ t represents the noise on the measurement sample due to various effects such as quantization. Furthermore, the individual sources may also be sampled at known subsets of the support Ω ′′ j ⊂ 1, T , j ∈ 1, J to obtain measured source samples, y ′′ jt , such that

y ′′ jt = s ′′ jt + b ′′ jt , ∀t ∈ Ω ′′ j , ∀j ∈ 1, J (2)
where b ′′ jt represents the noise for samples of each source. Lastly, for some problems such as declipping, we may also be given a set of C constraints, Γ ′′ c (s ′′ ), c ∈ 1, C , where each constraint, Γ ′′ c (s ′′ ), is in one of the following forms:

s ′′ jctc ≥ γ ′′ c , s ′′ jctc ≤ γ ′′ c , J j=1 s ′′ jtc ≥ γ ′′ c , J j=1 s ′′ jtc ≤ γ ′′ c (3)
in all of which γ ′′ c is a known constant and t c and j c are known time and source indices respectively.

Generalized Time Domain Audio Inverse Problem: Given all of the above definitions, we define the generalized audio inverse problem in time domain as that of recovering the sources, s ′′ jt ∀t,j (and hence their mixture), given the noisy and incomplete measurements, y ′′ jt ∀t∈Ω ′′ j ,∀j and {x ′′ t } ∀t∈Ξ ′′ , such that the constraints, {Γ ′′ c (s ′′ )} ∀c , are satisfied.

IV. THE PROPOSED APPROACH

A simple illustration of the generalized time domain audio inverse problem is shown in Figure 1, whereas the proposed algorithm in this work to solve this problem is illustrated in Figure 2. The individual steps of the proposed approach are explained in detail through the following subsections.

A. Redefining the Problem for a Frequency Domain Solution

The problem defined in Section III deals with constraints and unknowns in time domain, and as a result solving it with an approach that utilizes STFT domain constraints (such as the NTF model that will be introduced in Section IV-B) can be computationally heavy and even intractable. To rectify this issue, we will introduce the framed-time domain and STFT domain notations, using which, we will define a modified problem that is much easier to handle.

The framed-time domain (or sometimes called windowedtime domain) is the representation of the time domain signal after it is split into (often overlapping) frames of fixed length, F , and multiplied by a fixed windowing function. Assuming that the total number of frames is N , the notations x

′ f n , y ′ jf n , s ′ jf n , a ′ f n , b ′ jf n , Ξ ′ n ⊂ 1, F , Ω ′ jn ⊂ 1,
F represent the framed-time domain counterparts of the time domain notations defined in Section III for the source j ∈ 1, J , the intra-frame index f ∈ 1, F within the frame n ∈ 1, N . The relationships between the framed-time domain variables are similar to that of the time domain counterparts such that3

x ′ f n = J j=1 s ′ jf n + a ′ f n , ∀f ∈ Ξ ′ n , ∀n (4)
y ′ jf n = s ′ jf n + b ′ jf n , ∀f ∈ Ω ′ jn , ∀j, n (5)
We represent the STFT coefficients of the source signals simply by {s jf n } ∀j,f,n . Note that, the STFT coefficients are simply the Fourier transforms of the framed-time domain signals, such that

s jn = [s jf n ] ∀f = Us ′ jn ∀j, n, where s ′ jn [s ′
jf n ] ∀f and U is the normalized Fourier transform matrix satisfying UU H = U H U = I4 .

We now define a modified version of the initial problem using the framed-time domain variables and constraints, all of which can easily be computed from the time domain counterparts. This new definition of the problem has more relaxed conditions from the original problem in the sense that the problem is moved to a larger over-complete domain, and the correlation between the information within different frames is no longer defined. However it should be noted that the solution to this relaxed problem is not necessarily inferior to directly solving the time domain problem, because the models enforced in the STFT domain (such as the NTF model) to solve the relaxed problem often enforces a solution that is consistent among the frames. In the rest of this paper, we shall focus on solving this relaxed problem rather than the initial one.

Framed-Time Domain Audio Inverse Problem: We define our problem as that of recovering the sources in framedtime domain, {s ′ jf n } ∀j,f,n (or equivalently in STFT domain {s jf n } ∀j,f,n since they are related with a unitary transform) given the noisy and incomplete framed-time measurements, {y ′ jf n } ∀f ∈Ω ′ jn ,∀j,n and {x ′ f n } ∀f ∈Ξ ′ n ,∀n , such that the constraints, {Γ ′ c (s ′ )} ∀c , are satisfied.

B. Applying NTF Model with a GEM Algorithm

In order to make the problem described in Section IV-A easier to solve, we make a number of assumptions:

Assumption 1.
∼ N c (0, σ 2 a,jf n ), b ′ f n ∼ N c (0, σ 2 b,f n ), ∀j, f, n. (6) Assumption 2.
The sources are independently Gaussian distributed: Similarly, the unknown STFT coefficients of the sources, {s jf n } ∀j,f,n , are also independently distributed with zero mean complex valued Gaussian with variance {v jf n } ∀j,f,n , i.e.

s jf n ∼ N c (0, v jf n ), ∀j, f, n. (7)
Even though it is known that the noise in practice (such as quantization noise) is not always Gaussian, modeling the noise as Gaussian is still known to be a good enough approximation that provides significant computational advantage. Similarly the assumption of Gaussian distribution for the sources is also very common in audio community and accepted as a good approximation. It is noted when dealing with non-stationary signals that the assumption of gaussianity in the sources often results in very little loss in the source separation performance with the added benefit of much lower computational requirements [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. However the independence assumption for the source samples is not really accurate even though it greatly simplifies handling the joint probability distribution. Hence in this work we will make another assumption that will help us model the joint relationship of the source samples along the STFT plane and among the sources, with the help of nonnegative tensor factorization. Assumption 3. Variances of the sources form a low rank NTF structure: The tensor of source variances, [v jf n ] j,f,n , can be represented as the sum of few rank-1 nonnegative tensors, i.e.

v jf n = K k=1 q jk w f k h nk , ∀j, f, n (8)
with number of components, K, sufficiently small. This so-called PARAFAC/CANDECOMP [START_REF] Bro | Parafac. tutorial and applications[END_REF] NTF model can be parametrized by θ = {Q, W, H}, such that Q

= [q jk ] j,k ∈ R J×K + , W = [w f k ] f,k ∈ R F ×K + and H = [h nk ] n,k ∈ R N ×K + .
The assumption of a low rank NTF structure on the joint variances of audio sources is well known in the audio source separation community and it is shown to be an accurate model for audio signals in practice [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF], [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF].

We can define the observed mixture vector at frame n, x ′ n , and the observed source vector at frame n for source j, y ′ jn , as

x ′ n x ′ f n ∀f ∈Ξ ′ n ∈ R |Ξ ′ n |×1 , (9)
y ′ jn y ′ jf n ∀f ∈Ω ′ jn ∈ R |Ω ′ jn |×1 . (10)
Hence for each frame we can define the observed data vector, o ′ n , and each unknown source vector, s ′ jn , as

o ′ n y ′ 1n T , • • • , y ′ Jn T , x ′ n T T . ( 11
)
Hence, given the three assumptions above, we propose estimating the NTF model θ in the ML sense as

θ = arg max θ ′ p({o ′ n } ∀n |θ ′ ). (12)
To achieve that we employ a GEM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], while considering as latent data the totality of in general missing source STFT coefficients S = [s jf n ] ∀j,f,n . The algorithm iteratively alternates between an expectation step (E-step) for estimating the posterior power spectra of the signal and a maximization step (M-step) for updating the NTF model parameters. These two main steps can be summarized as follows:

• E-step: Estimate conditional expectations of source power spectra |s jf n | 2 , given the current model θ and the observations:

pjfn = E |s jf n | 2 o ′ n ; θ , ∀j, f, n. (13)
• M-step: Re-estimate NTF model parameters such that the 3-valence tensor of the NTF model approximation, V = [v jf n ] ∀j,f,n , is as close to the 3-valence tensor of estimated source power spectra, P = [p jf n ] ∀j,f,n , as possible with respect to the IS divergence [4]

D IS ( P V) = ∀j,f,n d IS (p jf n v jf n ), (14)
where d IS (x y) = x/ylog(x/y) -1, and pjfn and v jf n are as specified respectively by [START_REF] Adler | Audio inpainting[END_REF] and [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF]. The details of the E-step and the M-step are given in Sections IV-C and IV-D. In certain problems additional steps might also be required to satisfy certain constraints for the time domain signal or the STFT coefficients. It is described in Section IV-E how these additional constraints can be handled by the proposed algorithm. A summary of the overall algorithm is given in Algorithm 1.

C. E-Step: Estimating Posterior Statistics

Following our assumptions of independently Gaussian distributed signals, we can write the posterior distribution of each source frame s jn given the corresponding observed data o ′ n and the NTF model θ (or equivalently V = [v jf n ] ∀j,f,n with v jf n defined in (8)) as s jn |o ′ n ; θ ∼ N c (ŝ jn , Σ sjnsjn ) with ŝjn and Σ sjnsjn being, respectively, posterior mean and posterior covariance matrix of the STFT coefficients, s jn . These terms can be computed respectively by Wiener filtering as [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]

ŝjn = Σ H o ′ n sjn Σ -1 o ′ n o ′ n o ′ n , (15)
Σ sjnsjn = Σ sjnsjn -Σ H o ′ n sjn Σ -1 o ′ n o ′ n Σ o ′ n sjn , (16)
given the definitions of the covariance matrices

Σ o ′ n o ′ n =      Σ y ′ 1n y ′ 1n . . . 0 Σ H x ′ n y ′ 1n . . . . . . . . . . . . 0 . . . Σ y ′ Jn y ′ Jn Σ H x ′ n y ′ Jn Σ x ′ n y ′ 1n . . . Σ x ′ n y ′ Jn Σ x ′ n x ′ n      , (17)
Σ o ′ n sjn = 0 T L1,jn×F , Σ T y ′ jn sjn , 0 T L2,jn×F , Σ T x ′ n sjn T , (18)
Σ y ′ jn y ′ jn = U(Ω ′ jn ) H diag [v jf n ] ∀f U(Ω ′ jn ) + diag σ 2 b,jf n ∀f ∈Ω ′ jn , (19)
Σ x ′ n x ′ n = U(Ξ ′ n ) H diag ∀j v jf n ∀f U(Ξ ′ n ) + diag σ 2 a,f n ∀f ∈Ξ ′ n , (20)
Σ x ′ n y ′ jn = U(Ξ ′ n ) H diag [v jf n ] ∀f U(Ω ′ jn ), (21)
Σ y ′ jn sjn = U(Ω ′ jn ) H diag [v jf n ] ∀f , (22)
Σ x ′ n sjn = U(Ξ ′ n ) H diag [v jf n ] ∀f , (23)
Σ sjnsjn = diag [v jf n ] ∀f , (24)
where

U(Ω ′ jn ) is the F × |Ω ′ jn | matrix of columns from U with index in Ω ′ jn and L 1,jn j-1 l=1 |Ω ′ ln |, L 2,jn J l=j+1 |Ω ′ ln |.
The term diag(x) represents a diagonal matrix with the vector x along the diagonal.

Finally, the posterior power spectra, P = [p jf n ] ∀j,f,n can be computed as

pjfn = E |s jf n | 2 o ′ n ; θ = |ŝ jf n | 2 + Σ sjnsjn (f, f ). ( 25
)

D. M-step: Updating NTF Model Parameters

In order to find the best low rank NTF model estimation, V, of the given posterior power spectra tensor, P, we shall minimize the IS divergence between the two tensors, D IS ( P V) as defined in [START_REF] Griffin | Signal reconstruction from short-time Fourier transform magnitude[END_REF].

A common optimization approach to estimate the model parameters, θ, that minimizes ( 14) is using multiplicative updates (MU) as described in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]. In our case, starting from some initial nonnegative model parameters, the model parameters that minimize ( 14) can be found by applying several iterations of the following updates

q jk ← q jk f,n w f k h nk pjfn v -2 jf n f,n w f k h nk v -1 jf n , (26)
w f k ← w f k n,j q jk h nk pjfn v -2 jf n n,j h nk q jk v -1 jf n , (27)
h nk ← h nk f,j q jk w f k pjfn v -2 jf n f,j w f k q jk v -1 jf n . ( 28
)
In the beginning of the proposed GEM algorithm, the model parameters can be initialized randomly with nonnegative values. In the following iterations however, the update of the model parameters can be always applied starting from the current model parameters (instead of randomly initializing them before MU iterations each time P is updated). This not only reduces the number of iterations needed to apply the MU updates, but also helps avoiding local minima, since it is known that optimization based on MU rules is susceptible to getting stuck in local minima when initialized too far away from the solution.

E. Applying Additional Constraints

In many practical audio inverse problems, there may be additional knowledge on the signal to be estimated apart from the observed samples. We shall consider mainly two complementary types of knowledge on the signal to be treated, which provide:

i Constraints on STFT coefficients, such as specific harmonics being active, some frequency or time bins being silent, or simply the frequency response being symmetric (time domain signal being real valued); ii Constraints on framed-time domain samples, such as the constraints, {Γ ′ c (s ′ )} ∀c , that were defined earlier. Since our approach already constrains the STFT coefficients, enforcing additional constraints on STFT coefficients is usually as straightforward as modifying the model parameters, θ, accordingly during the MU iterations or simply initializing them in a specific way. For instance, the symmetry in frequency (hence being real valued in time) can be enforced if the matrix W is updated to be always symmetrical along the frequency axis. Similarly if the harmonics in the sources are known a priori, W can be initialized with a specific dictionary and then may never be updated to enforce using only these specific harmonics. Another example is, when certain entries of the matrices W, H and Q are known to be zero, they can be simply initialized to be zero and these zero values will be automatically enforced in the following MU iterations. Lastly, in certain applications, it is even possible to change the model to enforce additional structures on the matrices W, H and Q, such as sparsity by small modifications on the MU equations [START_REF] Le Roux | Sparse NMF? half-baked or well done?[END_REF].

Dealing with constraints on framed-time domain samples, unlike the constrains on the STFT coefficients, is not straightforward. When a framed-time domain sample is known to be clipped or quantized, the original value of this sample is known to be above (below) a certain threshold or to lay within a certain interval, and the resulting posterior probability distribution of the sample is no longer Gaussian. As a result, estimating the posterior power spectrum with this modified probability distribution is not as simple as described in Section IV-C. To overcome this problem, we propose to estimate the posterior power spectrum by computing the posterior mean, ŝjn and the posterior covariance, Σ sjnsjn , as described in Section IV-C, but then projecting them so as to satisfy the time domain constraints to obtain modified statistics, sjn and Σ sjnsjn respectively. As a result, the modified posterior power spectrum (to be used as the input for the NTF model update in M-step) is obtained as

pjfn = |s jf n | 2 + Σ sjnsjn (f, f ). ( 29
)
We define several approaches to compute the aforementioned modified statistics to satisfy the constraints in framed-time domain samples:

1) Unconstrained: The simplest way to perform the estimation is to ignore completely the constraints, treating the problem as a more generic audio inpainting in time domain. Hence during the iterations, the "constrained" signal is taken simply as the estimated signal, i.e. sjn = ŝjn , ∀n, j, as is the posterior covariance matrix, Σ sjnsjn = Σ sjnsjn , ∀n, j. 2) Ignored projection: Another simple way to proceed is to ignore the constraint during the iterative estimation process and to enforce it at the end as a post-processing of the estimated signal. In this case, the signal is treated the same way as in the unconstrained case during the iterations. 3) Signal projection: A more advanced approach is to update the estimated signal at each iteration so that the magnitude obeys the constraints. As an example, let us suppose we have a constraint in the form s ′ jcfcnc ≥ γ ′ c and it is not satisfied by the estimated posterior mean, i.e. ŝ′ jcfcnc < γ ′ c . We can simply set s′ jcfcnc = γ ′ c and s′ jf n = ŝ′ jf n for the rest of the support (and sjn = Us ′ jn ). Formally we can define,

{s ′ jf n } ∀j,f,n = argmin {z ′ jf n } ∀j,f,n ∀j,f,n |z ′ jf n -ŝ′ jf n | 2 s.t. {Γ ′ c (z ′ )} ∀c (30)
Note that this approach does not update the posterior covariance matrix, i.e. Σ sjnsjn = Σ sjnsjn , ∀n, j. 4) Covariance projection: In order to update the posterior mean and the posterior covariance matrix in a consistent manner, we can re-compute the posterior mean and the posterior covariance by ( 15) and ( 16) respectively, by treating the projected signal samples in [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF] at the support Ω ′ m,jn {f |s ′ jf n = ŝ′ jf n } as observed values for the current iteration. If the resulting estimation of the sources violates the time domain constraints on additional indices, those samples are also projected to obey the constraints and treated as observed. This process is repeated until a posterior mean, sjn , and a posterior covariance, Σ sjnsjn , that are consistent with all the time domain constraints, are obtained. Note that in addition to updating the pos- until convergence criteria met 9: end procedure terior covariance matrix, this approach also updates the entire posterior mean (or estimated signal) and not just the posterior mean at the indices of violated constraints.

V. EXPERIMENT RESULTS

The proposed algorithm is adapted to solve a number of audio inverse problems, which are listed in Section II. For each of these problems, we performed a set of experiments on various audio examples and compared the performance to that of known state of the art algorithms.

In the experiments below all the audio signals are sampled at 16 kHz, and the STFT within the various instances of the proposed algorithm is computed using a half-overlapping sine window of 1024 samples (64 ms).

A. Audio Declipping

In the audio declipping problem, we consider an audio signal with no known source information (as such it is modeled as a single source, J = 1) that is clipped to a known threshold of magnitude τ > 0. Thus the signal is accurately known for a subset of the support, Ξ ′′ , where signal magnitude is smaller than τ . For the remaining support, Ξ′′ = 1, T \ Ξ ′′ , the signal is unknown but obeys the time domain constraints of the form,

s ′′ tc ≥ τ, or s ′′ tc ≤ -τ. ( 31
)
We also assume that there is no observation noise, i.e., σ 2 a,f n = σ 2 b,jf n = 0, ∀j, f, n, in [START_REF] Nikunen | Object-based audio coding using nonnegative matrix factorization for the spectrogram representation[END_REF]. In [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF], various state of the art audio declipping algorithms are compared based on the experiments performed on music and speech examples. We have repeated these experiments using our approach with the same methodology and the datasets as reported in [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF] and provided an overall comparison of our algorithm to the other approaches. The experiment procedure can be summarized as follows; 10 music and 10 speech signals, each of length of 4 seconds, are scaled to have maximum magnitude of 1 in time domain, and then artificially clipped at eight different clipping thresholds (uniformly spaced from 0.2 to 0.9). The proposed algorithm is tested with four different methods to handle the clipping constraints as described in Section IV-E, namely Unconstrained (NMF-U), Ignored Projection (NMF-IP), Signal Projection (NMF-SP) and Covariance Projection (NMF-CP). The music signals are declipped with 20 NMF components (K = 20), while 28 components are used for speech signals (K = 28). The proposed GEM algorithm is run for 50 iterations. The performance of the proposed algorithm is compared to five state of the art methods: iterative hard-thresholding (HT) [START_REF] Kitić | Consistent iterative hard thresholding for signal declipping[END_REF], cosparsity (Cosp) [START_REF]Sparsity and cosparsity for audio declipping: a flexible nonconvex approach[END_REF], orthogonal matching pursuit (OMP) [START_REF] Adler | Audio inpainting[END_REF], social sparsity with empirical Wiener operator (SS-EW) and social sparsity with posterior empirical Wiener operator (SS-PEW) [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF].

The performance metric that is used to compare the algorithms is the improvement of the signal to noise ratio (computed only on the clipped regions) with respect to the clipped signal, SNR m , that is computed as [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF]:

SNR m = 10 log 10 ∀t∈ Ξ′′ |x ′′ o,t | 2 ∀t∈ Ξ′′ |x ′′ o,t -x ′′ e,t | 2 , (32)
where x ′′ o,t is the original time domain signal sample and x ′′ e,t is the estimated signal sample. Finally, the performance is measured in terms of the SNR m improvement, which is the difference between the SNR m computed on the estimated signal and the SNR m computed on the clipped signal.

The average performance of all the algorithms for declipping of music and speech signals is represented on Figure 3. It can be seen from the overall results that the proposed algorithm with the covariance projection (NMF-CP) has almost identical performance with the social sparsity based methods (SS-EW and SS-PEW) proposed in [START_REF] Siedenburg | Audio declipping with social sparsity[END_REF] while outperforming others. It can be also seen in the results that the model based algorithms (social sparsity and the NMF model) significantly outperform the methods relying on just sparsity (OMP and HT).

Regarding the effect of clipping constraints, the first thing to notice is that the performance of NMF-U with respect to NMF-IP (and NMF-SP) shows that simple constraints on the signal magnitude can noticeably improve the performance especially for music signals, hence they should not be ignored when possible. NMF-IP and NMF-SP are shown to have almost identical performance, even though the latter applies the constraints on the posterior mean of the signal at every iteration and the former simply applies a post processing to the final result. This observation combined with the superior performance of NMF-CP compared to the other methods demonstrates the importance of updating the posterior power spectrum more accurately for the success of the NMF-based methods.

Even though the performance is not better than the social sparsity approaches at first glance, the proposed algorithm has room for improvements in various aspects:

• NMF model can be easily extended to other, more structured NMF-like models such as source-excitation model or harmonic NMF [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. As shown in [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] Fig. 3: The average performance of all the audio declipping algorithms as a function of the clipping threshold (lower threshold corresponds to more severe clipping).

• It is shown in the results that the performance of our method depends significantly one the way the clipping constraint is handled. Therefore an alternative, more accurate computation of the posterior power spectrum might also improve the results further, whereas in dictionary based methods there is no approximation for the clipping constraints, hence performance improvement in this regard is not possible. It should be noted that dealing with constraints in both time domain and STFT domain comes at a computational cost in the Wiener filtering stage of the proposed algorithm. Luckily, this step is independent for each frame of the signal and hence can be easily parallelized, e.g., using graphical processing units (GPUs), to get a significant speed-up. On the other hand, estimating the signal independently within each window comes with the disadvantage that the estimation is not possible when there are no observed samples within a window. In practice, however, the loss of an entire window due to clipping is not probable for natural audio signals when the window size is chosen properly and the clipping threshold is not extremely low.

B. Joint Audio Declipping and Source Separation

The problem of joint audio declipping and source separation naturally arises when one would like to perform source separation on a mixture that is degraded due to clipping effects. Additionally one can also consider the case when source separation is not really needed, but a multi-source model is still employed to improve the performance of declipping when dealing with mixtures of different sources.

The global setup of our modeling to handle this problem is the same as the one for declipping in Section V-A, except that J > 1 sources are considered instead of just one. In order to assess the performance of declipping and source separation using the proposed algorithm, 5 different music mixtures 5 , each composed of 3 sources (bass, drums and vocals), are considered under 3 different clipping conditions. For each mixture with a maximum magnitude of 1 in time domain, 3 clipping levels at the thresholds of 0.2 (heavy clipping), 0.5 (moderate clipping) and 0.8 (light clipping) are considered, resulting in a total of 15 mixtures with different clipping levels. Each mixture is reconstructed by joint declipping and source separation, sequential declipping and source separation and only source separation ignoring the clipping artefacts. The proposed GEM algorithm (run for 100 iterations) has been used for all the reconstructions6 with K = 15 components. Inline with [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for userguided audio source separation[END_REF] and so as to inject some information about the sources to be separated, the sources in the mixtures are artificially silenced during a percentage of the total time. An example of the activation periods of the sources and corresponding indices set to zero in H during NTF model estimation are shown in Figure 4. Similarly Q is simply a J ×K matrix with a single 1 on each column and zeros everywhere else, which describes the assignment of the components to the sources.

The results of the simulations can be seen in Figure 5. Signal to noise ratio on the clipped support (SNR m ) for the declipped mixture is shown to demonstrate the declipping performance while signal to distortion ratio (SDR) is shown to demonstrate the source separation performance as described in [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF].

The results in Figure 5 show that when the clipping is severe, joint optimization is almost always preferable since it provides improvement on both the quality of the mixture and the quality of the separated sources with respect to source separation without declipping. This is as opposed to the sequential approach which provides comparable quality improvement in the mixture at the expense of the performance in source separation. In fact, for heavy clipping the declipping in sequential approach often reduced the performance of source separation noticeably with respect to separation without declipping. As the clipping gets lighter, the performance of sequential method approaches to that of joint method, and finally performs slightly better for light clipping. The joint optimization, however, still has few drawbacks which could be improved upon. The declipping in the sequential approach is performed with 15 components without any restrictions whereas the joint optimization is performed with the additional limitation that each source uses 5 components independently. Hence it is not possible that two sources share a common component in the joint optimization. This can be overcome by devising better methods to inject the prior information regarding the sources, see, e.g., [START_REF]Automatic allocation of NTF components for user-guided audio source separation[END_REF]. It should be also noted that the sequential optimization is approximately twice as fast as joint optimization due to handling much less complicated problems in either phases of the sequential processing. The fact that the Wiener filtering stage is independent for each window and can be parallelized to provide significant speed improvements, can be helpful to overcome this problem in the future.

C. Compressive Sampling Recovery

In order to demonstrate the ability of the proposed approach to reconstruct randomly subsampled signals, we have randomly subsampled a typical music signal of 4 seconds at different average rates (percentage of samples at 2, 4, 8, 16, 32), and then reconstructed with our algorithm in a similar fashion to the experiments in Section V-A, but without any clipping constraints (hence J = 1 and σ 2 a,f n = σ 2 b,jf n = 0, ∀j, f, n). The reconstruction is performed with different number of components (K = 2, 8, 24, 32, 48, 72) in order to observe the sensitivity of the results to the parameter K. In order to provide a reference for the reconstruction capability of our algorithm, the results with shape preserving piecewise cubic interpolation are also provided 7 .

The reconstruction results can be observed in Figure 6. The first thing to notice is that the reconstruction results with the proposed algorithm (solid lines) are significantly better than the results with simple interpolation (dashed lines) as expected. Another noticeable behaviour in the results is that once the number of components, K, is sufficiently large, the reconstruction performance does not seem to suffer. This behaviour is unlike what we have observed for other problems such as declipping, for which the choice of number of components is an important factor for obtaining best performance. Looking 7 For the interpolation, the interp1() function of Matlab 2016a is used with phcip method, which gave the best results among the available interpolation methods. more closely to the estimated NMF components, we have seen that the maximum likelihood estimate combined with random sampling already provided a strong prior for signal estimation and the benefit from low rank model was minimal in this case. Hence, as long as the number of components are chosen sufficiently large, the estimated variances, V, are effectively independent of K.

D. Informed Source Separation

The goal of ISS is essentially to compress multiple audio sources so that they can be efficiently decoded when the mixture of the sources is available. Our proposed algorithm can also be used for this purpose when combined with random subsampling and quantization. Thanks to the efficiency of the signal estimation with the proposed approach, the sources can be quantized and randomly subsampled at the compression stage, and from these samples, the proposed algorithm can estimate the sources at the decoding stage. A summary of this scheme is shown in Figure 7. We call this approach compressive sampling-based informed source separation (CS-ISS) because it exploits the principles of compressive sampling through random sampling of a tight frame transform domain (namely STFT domain) while having a compact representation in another domain (namely low rank NTF representation of the power spectra of the sources).

In order to assess the performance of our approach, three (J = 3) 11-second long sources of a music recording are encoded and then decoded using the proposed CS-ISS with different levels of quantization (16 bits, 11 bits, 6 bits and 1 bit) and different raw sampling bitrates8 per source (0.64, 1.28, 2.56, 5.12 and 10.24 kbps/source). Since uniform quantization is used, the noise variance in time domain is is pre-defined and known during both encoding and decoding. The quantized samples are truncated and compressed using an arithmetic encoder with a zero mean Gaussian distribution assumption. At the decoder side, following the arithmetic decoder, the sources are decoded from the quantized samples using 50 iterations of the GEM algorithm with the number of components fixed at K = 18, i.e. in average 6 components per source. The quality of the reconstructed samples is measured with SDR as described in [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. The resulting encoded bitrates and SDR of decoded signals are presented in Table I along with the percentage of the encoded samples in parentheses.

σ 2 = ∆ 2 /
Note that the compressed rates in Table I differ from the corresponding raw bitrates due to the variable performance of the entropy coding stage, which is expected.

The performance of CS-ISS is compared to a classical ISS approach with a more complicated encoder and a simpler decoder presented in [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF], as well as much better performing coding based approach proposed in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF]. Both the classical ISS and coding-based ISS algorithms are used with NTF model quantization and encoding in a similar fashion as in the experiments described by [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF], i.e., NTF coefficients are uniformly quantized in logarithmic domain, quantization step sizes of different NTF matrices are computed using equations ( 31)-( 33) from [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] and the indices are encoded using an arithmetic coder based on a two-state Gaussian mixture model (GMM) (see Fig. 5 of [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF]). The approach is evaluated for different quantization step sizes and different numbers of NTF components, i.e., ∆ = 2 -2 , 2 -1.5 , 2 -1 , . . . , 2 4 and K = 4, 6, . . . , 30. The results are generated with 250 iterations of model update. The performance of CS-ISS and the earlier approaches are shown in Figure 8 in which CS-ISS clearly outperforms the classical ISS approach and is on par with coding-based ISS approach, even though both of these approaches can use an optimized number of components as opposed to our decoder which uses a fixed number of components (the encoder is very simple and does not compute or transmit this value). The performance difference with classical ISS is due to the high efficiency achieved by the CS-ISS decoder thanks to the incoherency of random sampled time domain and of maximum likelihood estimation along with low rank NTF model. Also, the classical ISS approach [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF] is unable to perform beyond an SDR of 10 dBs due to the lack of additional information about STFT phase as explained in [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF]. The results indicate that the rate distortion performance exhibits a similar behaviour as to the coding-based ISS algorithm. It should be reminded that the Fig. 8: The rate-distortion performance of CS-ISS using different quantization levels of the encoded samples. The performance of the ISS algorithm from [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF] and the codingbased ISS algorithm from [START_REF] Ozerov | Coding-based informed source separation: Nonnegative tensor factorization approach[END_REF] are also shown for comparison.

proposed approach distinguishes itself by its low complexity encoder and hence can still be advantageous against other ISS approaches with better or seemingly equivalent rate distortion performance. The performance of CS-ISS in Table I and Figure 8 indicates that different levels of quantization may be preferable in different rates. Even though neither 16 bits nor 1 bit quantization seem well performing, the performance indicates that 16 bits quantization may be superior to other schemes when a much higher bitrate is available. Coarser quantization such as 1 bit, on the other hand, had very poor performance in the experiments. The choice of quantization can be performed in the encoder with a simple look up table as a reference. One must also note that even though the encoder in CS-ISS is very simple, the proposed decoder is significantly high complexity, typically higher than the encoders of traditional ISS methods. However, this can also be overcome by exploiting the independence of Wiener filtering among the frames in the proposed decoder with parallel processing, e.g., using graphical processing units (GPUs).

VI. CONCLUSIONS

In this paper, we have presented a novel approach for time domain signal estimation in the maximum likelihood manner which can be used for applying the low rank NTF model to the power spectrum of the signal in many types of problems that were never before solved using the NMF/NTF model. The presented algorithm is demonstrated to be very effective for several different audio inverse problems while providing multiple advantages to other existing methods. For the audio declipping problem, clipped sections of music and speech signals are restored using the proposed approach as well as state of the art methods, and the proposed algorithm is shown to be highly competitive while providing complementary advantages. For audio source separation and mixture declipping, the proposed algorithm is shown to be capable of jointly solving these two separate problems which was not possible with any other method in the literature. Joint handling of these problems is also demonstrated to be more effective than sequentially approaching each problem in case of severe distortions. The proposed algorithm is also shown to be highly effective for the reconstruction of randomly subsampled signals such as in the case of compressive sampling approaches. This advantage of our algorithm is further utilised for the problem of informed source separation, to create a compression scheme which uses the principles of compressive sampling and distributed coding For this application, the proposed algorithm is not only shown to achieve compression performance equivalent to that of the state of the art, but also shown to have unique advantages, specifically having a very simple encoder as well as the decoding stage being independent of the encoding stage.

The NMF and NTF representations are gaining a lot of popularity in signal modelling community and we see the algorithm presented in this paper to be a step towards the application of these models to a wider class of signal estimation problems. Even though the provided examples in this paper are from audio inverse problems, the proposed algorithm is by no means limited to audio applications, but any application for which a low rank NMF/NTF model is an accurate representation for the power spectrum.

We consider several improvements and extensions to the proposed algorithm as future work. An extension to multichannel audio is an interesting step for dealing with real world audio problems. Furthermore, adapting the proposed algorithm for imaging problems with multiple additive components, such as imaging through transparent and reflective surfaces, is another intriguing direction.

Fig. 2 :

 2 Fig. 2: Illustration of the proposed algorithm to solve the generalized time domain audio inverse problem.

  Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (a) Declipping performance for clipping level 0.2. Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (b) Source separation performance for clipping level 0.2. Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (c) Declipping performance for clipping level 0.5. Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (d) Source separation performance for clipping level 0.5. Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (e) Declipping performance for clipping level 0.8. Sep. Sequential Declip./S-Sep. Joint Declip./S-Sep. (Proposed) (f) Source separation performance for clipping level 0.8.

Fig. 5 :

 5 Fig. 5: The declipping and source separation performance of joint optimization compared to sequential.

Fig. 4 :Fig. 6 :

 46 Fig.4: Depiction of the experiment set-up for the injection of information on different source characteristics. The information on the known silent time durations of each source (represented in different colors) is directly utilized by setting the corresponding coefficients in H to 0. Note that the matrices H and W are formed by concatenating the Hs for each source and W s for each source (depicted above) along the component dimension respectively.

Fig. 7 :

 7 Fig. 7: The encoding and decoding processes for the compressive sensing-based informed source separation.

  -16 bits (proposed) CS-ISS -11 bits (proposed) CS-ISS -6 bits (proposed) CS-ISS -1 bit (proposed)

  The noise is independently Gaussian distributed with known variance: The noise time samples for the observations, {a ′ jf n } ∀j,f,n and {b ′ f n } ∀f,n , are independently distributed with zero mean Gaussian with known variances, {σ 2

a,jf n } ∀j,f,n and {σ 2 b,f n } ∀f,n respectively, i.e. a ′ jf n

  Estimate ŝjn , Σ sjnsjn , ∀n, j, given θ, o ′

	Algorithm 1 GEM algorithm for solving Time Domain Audio
	Inverse Problems with NTF model
	1: procedure RESTORE-AUDIO-WNTF
	2:	Initialize nonnegative θ = {W, H, Q} randomly
	3:	repeat
	4:	E-step : jn
		∀n, j	⊲ see § IV-C
	5:	Time domain constraints : Estimate sjn , Σ sjnsjn ,
		∀n, j and P given {Γ ′ c } ∀c	⊲ see § IV-E
	6:	M-step : Update θ given P	⊲ see § IV-D
	7:	STFT domain constraints : Modify {W, H, Q}
		as needed	⊲ see § IV-E
	8:

  [START_REF] Simsekli | Score guided audio restoration via generalised coupled tensor factorisation[END_REF] where ∆ is the quantization step size (henceσ 2 b,jf n = ω 2 f ∆ 2 /12, where ω 2 f are the framing (or STFT) window coefficients. The mixture is available in entirety at the decoder, therefore the noise variance of the mixture is zero (σ 2 a,f n = 0). It is assumed that the random sampling pattern 73%) 1.75 / 13.30 dB (1.45%) 3.50 / 19.47 dB (2.91%) 7.00 / 24.66 dB (5.82%) 6 bits 0.27 / 4.17 dB (0.67%) 0.54 / 7.62 dB (1.33%) 1.08 / 12.09 dB (2.67%) 2.18 / 14.55 dB (5.33%) 4.37 / 16.55 dB (10.67%) 1 bit 0.64 / -5.06 dB (4.00%) 1.28 / -2.57 dB (8.00%) 2.56 / 1.08 dB (16.00%) 5.12 / 1.59 dB (32.00%) 10.24 / 1.56 dB (64.00%)

	Bits per			Raw rate (kbps / source)
	Sample	0.64	1.28	2.56	5.12	10.24
			Compressed Rate / SDR (% of Samples Kept)
	16 bits 0.50 / -1.64 dB (0.25%) 1.00 / 4.28 dB (0.50%)	2.00 / 9.54 dB (1.00%)	4.01 / 16.17 dB (2.00%) 8.00 / 21.87 dB (4.00%)
	11 bits 0.43 / 1.30 dB (0.36%)	0.87 / 6.54 dB (0.

TABLE I :

 I The final bitrates (in kbps per source) after the entropy coding stage of CS-ISS with corresponding SDR (in dBs) for different (uniform) quantization levels and different raw bitrates before entropy coding. The percentage of the samples kept is also provided for each case in parentheses. Results corresponding to the best rate-distortion compromise are in bold.

For sake of simplicity, we only consider the single-channel case here. The proposed algorithm in this paper can be readily extended to multichannel case in a similar way as it is done in[START_REF] Ozerov | Multichannel audio declipping[END_REF] for the declipping application.

Throughout this paper the time domain signals will be denoted by letters with two primes, e.g., x ′′ , the framed-time domain signals by letters with one prime, e.g., x ′ , and complex-valued STFT coefficients by letters with no prime, e.g., x.

From this point on, we shall use simply ∀n to denote ∀n ∈ 1, N , ∀f to denote ∀f ∈ 1, F and ∀j to denote ∀j ∈ 1, J , unless a subset of these sets is specified, e.g. Ξ ′ n .

x T and x H represent the non-conjugate transpose and the conjugate transpose of the vector (or matrix) x respectively.

The mixtures are taken from the "professionally produced music recordings" task dataset of SiSEC 2015 source separation evaluation campaign (https://sisec.inria.fr/sisec-2015/).

For declipping only, the algorithm is used with a single source (as in Sec. V-A), and for source separation only, the algorithm is used with the observed support set being the entire time axis.

The raw sampling bitrate is defined as the bitrate before the entropy encoding step.


Dongguang Li
email: d.li@ecu.edu.au

Measurements of Luminescence from Cleaved Silicon

This paper outlines the results from experiments performed to gain further information about the structure and properties of cleaved silicon surfaces, using vacuum cleavage luminescence detection methods.

The experiments involved detecting the luminescence produced by cleaving thin silicon plates within a high vacuum, by a process of converting the luminescence to an amplified electrical signal.

INTRODUCTION

It was discovered that single-crystal silicon, Si, when cleaved in a vacuum, emitted light in at least two wavelength regions at low intensities with duration of hundreds of microseconds [1] [2]. This paper outlines the results from experiments performed to gain further information about the structure and properties of cleaved silicon surfaces. The experiments involved detecting the luminescence produced by cleaving thin silicon plates within a high vacuum, by a process of converting the luminescence to an amplified electrical signal. The experiments were based on the assumption that surface cleavage and reconstruction may cause electrons to become excited, and the resulting recombination process which involves the combining of an electron in the conduction band with a hole in the valence band, would result in an emission of energy that would be detectable. This hypothesis was supported by simple calculations that predicted that every broken atomic bond on a silicon surface should radiate one photon, thus generating a strong detectable emission signal. Assuming the energy would be approximate the bond breaking energy of 2.1 eV, the corresponding wavelength was calculated to be within the visible region of the spectrum at 590 nm. Silicon has an indirect bulk band gap therefore band edge transitions are indirect and mostly non-radiative. The radiative changes in bond-energies during the surfaceatom reconstruction process had been observed and was found to be much less than anticipated, falling instead within the infrared region.

CLEAVED SURFACE STRUCTURE OF SILICON

Silicon is a prominent material used in semiconductor devices therefore it is important to understand its surface structure and properties. Although the bulk structure is well known, the complex rearrangements of the first few layers of atoms of clean surface structures have been more difficult to determine. Clean surfaces have generally been obtained by using techniques such as high temperature thermal flashing, or ion sputtering followed by thermal annealing in an ultra high vacuum (UHV), but an alternative method has been to use direct sample cleaving.

Haneman [3] outlines various methods used to investigate the microstructures of clean cleaved silicon surfaces, including scanning tunneling microscopy (STM) and lowenergy electron diffraction (LEED). Data obtained from these techniques appear to be consistent with buckled dimer structures for the (100) 2x1, 2x2 and 4x2 surfaces, and with the dimer-adatom-stacking fault (DAS) model for the annealed (111) 7x7 surface [3]. However other low-index surface structures of silicon are not as well understood. In particular, the structure of silicon surfaces produced by cleavage in ultrahigh vacuum shows a (2x1) cell, which usually converts upon annealing to an intermediate (5x5) structure at about 370 o C, depending upon the annealing time. This structure then converts to the (7x7) form at around 550 -600 o C. Measurements of angle-resolved photoemission [4] and polarized optical reflectivity [5] and absorption [6] were consistent with the surface chain structure model proposed by Pandey [7] to explain the large (0.8 eV) dispersion deduced from the photoemission in the Γ-J direction of the two-dimensional Brillouin zone. This model appeared to be consistent with the STM measurements and with the ion scattering data after buckling modifications were introduced, but it wasn't supported by the LEED measurements [8]. Upon close examination, the Pandey model was found to depend upon an assumed calculation of the presence of very strong shear forces existing within the first few layers for which no origin could be found in the cleavage process. The experimental results gave indication of new bond overlaps forming while sub-surface bonds were rupturing.

Therefore, the Three-Bond Scission (TBS) Model, which possesses the surface chain features required to satisfy the optical data, is proposed. The TBS model assumes that the cleavage passes through the three-bond plane, not through the single-bond (111) plane.

EXPERIMENTAL DESIGN

As shown in Fig. 1, the specimen was supported in a specimen holder at the center of curvature of an aluminium vacuum chamber, with a detector located 5-8 mm above the specimen. After cleaving, the resultant luminescence signal detected was amplified by a low noise amplifier then fed into a digital storage adaptor (DSA) that was triggered by the emission with a set pretrigger period. By incorporating several choices of time base and sensitivities, the signal was captured and directly plotted using an x-y recorder. A cathode ray oscilloscope (CRO) was used for display, and signals were stored and processed on a computer system (CMP).

The system was calibrated by testing with a pulse stroboscope (flash) and a light emitting diode (LED) driven by a waveform generator which imitated one shot light signal.

Obtaining High Vacuum

Assuming a unit sticking coefficient, whereby the probability that a molecule striking the surface will stick is 100%, contamination of the sample by a layer of gas forming a compact monolayer (approx. 8 x 10 14 molec/cm 2 ) occurs at a pressure of 10 -3 Torr in 2.2 ms and at 10 -6 Torr in 2.2 s within an unrestricted supply of ambient gas. These times are extended in systems with volumes of a few litres for surfaces of the order of a few cm -2 .

Cleavage luminescence usually has the maximum duration of about 1 ms and it is possible to carry out cleavage experiments at high vacuum pressures of 10 -3 to 10 -6 Torr before contamination happens, provided sticking coefficients significantly less than 1 are involved as is usually the case with semiconductor surfaces. A high vacuum condition was achieved by using three different pumps (Rotary, Diffusion and Ion), with the pressure measurable by the current drawn from the ion pump power supply and by an ionization gauge.

Installation of Specimen & Detecting Unit

Cleavage was performed in vacuum which required both specimen and detector to be sealed hermetically in the vacuum chamber, with the specimen fixed on a sample holder combined with a cleaver in the bottom of the vacuum chamber. In room-temperature detectors, an electrical feedthrough mounted on the side of the vacuum chamber could be used for the leading-out of wires of the detector. But the larger liquid-nitrogen-cooled-detectors required a perspex disk of 160cm diameter to replace the chamber cover as a detector holder, with a hole drilled in the center of the holder to fit the outer diameter of the detector case, so that the detector could be located exactly over the sample to be cleaved. This sealing method was used between the detector case and the detector holder as well as between the detector holder and the vacuum chamber. To overcome the mechanical vibration transferred by the metal case which could seriously distort the true signal and reduce the ratio of signal-tonoise, a small size room temperature detector was suspended by fine wires from the Perspex holder over the sample, which effectively obstructed the vibration from the cleaver.

Cleavers

The cleaver design had to reduce the mechanical vibration to a minimum level when the sample was being cleaved. Cleavage can be achieved by bending, pulling or pulling-bending, and several cleaver types were used to allow for several experimental requirements as crack behaviour depended heavily on the cleavage method used.

The block bending cleaver was found to be the best for reducing mechanical vibration. Other types of cleavers used included the paddle bending cleaver, the Wilson seal pulling cleaver, the linear motion feedthrough (LMF) pulling cleaver (Fig. 2) and the paddle pulling-bending cleaver (Fig. 3). The two or more time constants, T 1 and T 2, observed for the Si bulk recombination may be related to band structure, different photon cross sections or the nature of

Dongguang Li Measurements of Luminescence from Cleaved Silicon

the surface as well (tunneling processes). A possible explanation is that after electrons are excited by cleavage from valence band to conduction band, the radiative recombinations occur through several different paths rather than only one process, which results in several overlapped luminescence signals with different time constants. Experimentally, only a filter with a peak wavelength was 400 nm could effectively block the Si signal, so the conclusion was drawn that the Si signal A has a wide wavelength range of 460-1040 nm (or 1.1-2.7 eV).

Signal B

The luminescence wavelength from the second cleavage signal, B, corresponds to the recombination of cleavageexcited carriers occurring across the surface state gap. Signal B is harder to detect because the longer wavelength photovoltaic or photoconductive detectors have less detectivity by nearly 3 orders of magnitude than Si photovoltaic detectors. The requirement of keeping the detectors at liquid nitrogen temperature also meant that spurious long-time signals (seconds) were radiated to the detector, caused by vibration of the room temperature parts during cleavage, but these were easily detected and recognized, as indicated in Fig. 6. Problems of background radiation were overcome by using optical interference filters (coarse and segment) mounted directly onto the sapphire window that sealed the vacuum enclosure over the liquid-nitrogen detector. Our experiments then detected the B radiation and showed that the durations were of length of several hundred microseconds, with the average length 200 µs.

A precise value for the wavelength range of Signal B was obtained from an extensive series of experiments with optical interference filters placed over the liquid-nitrogencooled InSb photovoltaic detector.

To determine the energy ranges accurately, the transmittance of Segment 2 with the mask in place was measured with a Cary 17 spectrophotometer. The peak transmittance was 20% with the FWHM (full-width at half maximum) points occurring at 0.266 and 0.253 eV. The points of 10% of peak transmittance occurred at 0.268 and 0.248 eV. Therefore the energy range of Signal B was determined from this to be 0.26 + 0.01 eV.

Dongguang Li Measurements of Luminescence from Cleaved Silicon

DISCUSSION

The silicon bulk band gap radiation A is explained by assuming that electrons are excited by the cleavage process into the conduction band whence they decay. Since Si has a high density of surface states, there is a strong surface field which would drive conduction electrons into the interior so that their luminescent decay should not be affected by the presence of any nonradiative paths introduced by contamination at the surface. However, the experiments show that the presence of air quashes the radiation, forcing a conclusion that the excited electrons are held at the surface, which could be regarded as a small potential well. The excited electrons would also occupy surface states. From direct measurements of surface reflectivity [16] and surface absorption [6] on cleaved surfaces with polarized light, it is known that there is a band gap between the surface states with a peak optical-absorption coefficient at approximately 0.45 eV. Hence the observed emission B in the range 0.25-0.27 eV could well be due to transitions across the surface-state gap. Since the duration of B is about 200 µs, it is inferred that the radiative recombination must be indirect. In addition to surface state emission, the role of defects has to be considered, since some cleaved Si surfaces display properties that have been ascribed to these (usually interpreted as steps) although other imperfections are also possible. The energy states of such "bad cleaves" have been characterized by optical measurements [16] [17], electron energy loss spectroscopy [18][19], pumped photoemission [20][21] and scanning tunneling spectroscopy [22], among others. These agreed in general that there exists a band of defect states with an onset absorption energy at room temperature at about 0.2 eV and a peak at about 0.35 eV. In the case of cleaved surfaces subjected to intense laser irradiation, measurements by Boker et al [20] [21] found short decay times of a fraction of a nanosecond. Very significantly, the lifetimes were shorter for the cases of poorer cleaves.

The position of the defect band with respect to the band edges can be estimated from knowing that the Fermi level is at 0.4 eV above the bulk valence band edge [23], and that the upper surface state band has a peak density at about the Fermi level for the very heavily doped Si (ntype, 8x10 18 /cm -3 ) [14] and about 0.2 eV higher for normally doped material where the band is known to be much less occupied. The lower surface state band has a direct gap from the upper of 0.45 eV. These figures place the top of the lower surface state band at about the top of the bulk valence band. The defect band has a peak at about 0.35 eV above this, i.e. close to the Fermi level. It is obvious that the energy of 0.26 eV does not fit very well with the peak energies of the defect states (0.35 eV), but could correspond to the onset region of the band, reported as about 0.2 eV [19]. However, if we adopt this interpretation, it must be explained how it is possible to have such a long emission time (200 µs) from localized defect states. One would have to regard them as trap states with an energy barrier against transition to a state from which fast emission might be possible. From the discussion above of the position of the defect band, the only states that are readily accessible are the upper surface states. According to the measurements of Bokor et al [20][21] the decay durations of electrons in the upper band are about 2 ns even on defect surfaces. However, in those experiments, the surface was excited into photoemission by a photon pulse of 60 ps duration and power approximately 10 7 watts cm -2 . These conditions differ greatly from ours in that the surface is highly excited with a large supply of holes in the valence bonds and of phonons, which are missing in our case. The extra holes and phonons can account for a greatly enhanced transition probability from the upper to the lower surface band. The conditions used by us, in which excitation is due in effect to the energy supplied by lattice vibrations initiated by the bond rupture, will be more characteristic of unexcited Si.

We also note that the result (200 µs) is reasonable for an indirect surface state gap in that long lifetimes are characteristic of device-grade Si, the bulk indirect-gap lifetime being milliseconds. As discussed previously, the data can thus be interpreted as arising, at least in part, from conduction-band electrons dropping to the surface-state conduction-band minimum and there making transitions to the surface-state valence band. The lowest-energy transitions are indirect, requiring phonons of almost the entire zone-width momentum, and thus accounting for the long observed durations. The data thus supports the feature of the surface state theories in that there is an indirect gap. This feature could not be detected by optical absorption measurements on surface states since it would be too weak.

From the above considerations, the Signal B is explained as due to electrons decaying from the bottom of the surface state conduction band. They reach that band either directly or via defects. In the latter case, the long decay times could be due to trapping times. In the former, they are due to the indirect nature of the surface state gap. In either case, the measured value of 0.26 + 0.01 eV would provide a determination of the minimum gap. With respect to the signal C, there are two related questions: what is the origin of luminescence C and why is it apparently more readily detected during pulling than during bending? At this stage, it is ascribed to recombination at defects on cleaved surfaces since these could have short recombination times. It is known that cleaved surfaces can have, for example, paramagnetic defects in the shape of small partial splits under step edges [24]. Various point defects, mainly vacancies, appear to have been seen in scanning tunneling microscopy study [22][25]. There is as yet, insufficient information to enable one to correlate any particular defect with the mode of cleavage but these results suggest that there is some correlation and that pulling and bending cleavage give noticeably different luminescence results. A microscopic search for topographic differences on many cleavage surfaces formed by the different cleavage techniques failed to find general differences. All surfaces had areas of tear marks, river patterns, and large mirror-like regions. Therefore, we cannot at this stage correlate atomic-scale differences with micron-scale features.

Fig. 1 .

 1 Fig. 1. Block diagram of experimental system

Fig 2 .

 2 Fig 2. LMF Pulling cleaver

Fig 4 . 1 =

 41 Fig 4. Types of detectors used. Front L-R: Si (36mm 2 ), Ge (9mm 2 ), Si (1mm 2 ) Back L-R: Ge (20 mm 2 ); MCT (4mm 2 ) ; Si(100mm 2 ); InSb (13mm 2 )4. DETACTED SIGNALS

Fig 5 .

 5 Fig 5. Cleavage luminescence from Si -two time constants

Fig 6 .

 6 Fig 6. Cleavage luminescence from Si on MCT detector.

4. 3 .

 3 Signal C Signal C was simultaneously found with the usual signal A when faster amplifiers (amplifier B or C) were used, but it is completely different from A in having a duration of approximately 15 µs rather than about 300 µs for A. Its true duration may in fact be less, as the progress of the crack through the Si wafer would cause continuous emission of light for the duration of the crack. To detect Signal C, it was necessary to use a very fast amplifier with a sufficiently short response time, although it meant a loss in signal-to-noise performance for any longer signals. A three-resistor-network feedback unit was used to replace the original 10 MΩ single feedback resistor of the first stage of the amplifier. Signal C was only detectable when cleavage was initiated by longitudinal pulling (i.e. by using the Wilson seal pulling cleaver or LMF pulling cleaver). Optical filter experiments were carried out to determine the exact wavelength range of the signal C. Among seven optical narrow band-pass filters used, four passed signal C and three blocked it. Experimental results obtained indicated that signal C from cleaved Si also has a wide wavelength spread of 660-1040 nm (or 1.1 -1.88 eV), compared to signal A which has an even wider wavelength range (0.46 -1.1 µm), a characteristic in addition to the difference in duration, which also distinguished signal C from signal A.

Fig 7 .

 7 Fig 7. Cleavage luminescence from Si (using pulling cleavers)