Patent Description:
The signal is, in general, a specific type of information representable as a value varying over a signal domain (which may be time and/or space and/or another suitable field), usually in the form of an array of signal values ("data") defined over the signal domain. Generally, the signal will be a technical signal that originated from some technical source, such as a measuring device, recording equipment, imaging system, or other technical signal generator, and is usually encoded as a digital or digitalized signal. A gap in the signal is usually due to some (limited) disturbance, such as a technical fault, noise, interference, which compromises the signal in the gap during or after measurement or recording thereof. One important special case is an audio signal; and the filling of audio signals is also known as audio inpainting. The term "inpainting" has spread to the audio processing community from the image processing field, and is used for filling degraded or missing parts of data, based on its reliable part, and preferably in an unobtrusive way. Audio inpainting is e.g. discussed in the article of Adler et al. [<NUM>], as a signal processing technique to restore gaps, i.e. missing consecutive samples, in an audio signal, while still keeping perceptible audio artifacts as small as possible.

Audio inpainting has attracted a lot of commercial and research interest as it has many important applications. The most common approach to tackle the audio inpainting problem is to exploit prior information about the signal, i.e. to utilize reliable parts of the signal to inpaint the gaps. The reliable parts could be either exploited in the original (time) domain or with respect to some redundant dictionary. In the time domain setting, known inpainting methods (such as those described in the articles Janssen et al. [<NUM>], Oudre [<NUM>], and Etter [<NUM>]) are based on autoregressive (AR) signal modeling; the missing samples are filled by linear prediction using autoregressive coefficients that are learned from the neighborhood of the gap. Due to their excellent performance they are still state-of-the-art.

Other audio inpainting techniques exploit sparsity of the signal. These methods, such as those described in [<NUM>] or in Mokry et al. [<NUM>], formulate the inpainting problem as an optimization task, while taking into account that many real-world audio signals have an (approximately) sparse representation with respect to a suitable dictionary. However, the known sparsity-based inpainting techniques often suffer from a drop of the signal energy within the filled gap, which is the more pronounced the larger the gap is. Moreover, the computing expenses are often considerable.

In view of the above, the invention aims at overcoming the shortcomings of the prior-art reconstructions and inpainting methods. In particular, it is desired to increase the reliability of the restored signal in the gap and provide improved reconstruction performance, while keeping the computational burden moderate. The invention has a primary focus on time-dependent signals, such as audio signals, but is also intended in relation to more general types of input signals, including image and video signals.

The invention is described in independent claims <NUM> and <NUM>. Further optional features of the invention are described in the dependent claims.

In particular, the invention provides a method for signal processing which starts from an input signal dependent on location within a predetermined signal domain in the form of a plurality of signal values associated with respective locations contained in a containing region within the signal domain (the containing region is the range which contains all locations of the signal values and it may be the same as the signal domain or a part of it), where the input signal comprises a gap, and the method comprises the following steps:.

This solution according to the invention proposes a novel approach which is based on modeling of the dictionary to exploit and, if possible, increase sparsity of the signal. The invention starts from the observation that in practice, signals usually have the property of being sparse with respect to a suitable frame representation in the time-frequency framework, such as Short-Time Fourier Transform (STFT), which is also known as Gabor transform, or more generally in the location space/reciprocal space framework. Herein, "sparsity" of an array (i.e. vector or matrix) means that most entries in the array are zero. The measure of sparsity is the number of zero entries, which - except for an inverted sign - corresponds to taking the ℓ<NUM>-norm: the lower the norm, the higher the sparsity. In the context of the invention, the concept of sparsity is relaxed towards using other suitable norms, such as a ℓ<NUM>-norm or ℓ<NUM>,∞-norm. The invention includes a method for modifying an initial dictionary which uses the reliable parts of the input signal around the gap with the goal to obtain a signal representation (a modified dictionary) which has increased sparsity or, in other words, reduced norm. This is achieved by means of a "sparsifying transform" which is an optimization method that "deforms" the elements of the given initial frame (e.g. Gabor vectors) such that the sparsity of the analysis coefficients of the resulting frame is maximized - i.e. its norm is minimized. The deformed frame thus obtained is then used with an inpainting method of known type to fill the input signal in the gap.

Thus, the invention offers the technical effect of an improved method for processing a signal to fill in gaps in the signal, with higher efficiency of the process and reconstruction performance, so as to lower the burden of computational calculations, while increasing the reliability of the restored signal in the gap and reducing computational complexity.

The signal values in the gap are treated as invalid. Depending on the type of signal, invalid values may be represented as zero, or another value which is used to indicate non-reliablity of the value. Also, it is possible that invalid values are identified by an additional marker attributed to the values, indicating the status of the values as "invalid" or "non-reliable" rather than "valid"/"reliable". In many embodiments, however, the actual value of the signal values in the gap is irrelevant because the gap is defined by parameters presented external to the input signal, and the signal values lying in the gap are not used.

The term "signal domain" is used herein to denote the parameter region on which the signal is defined; the term "location" is to be understood broadly and may also include location in time etc..

In a first class of the invention as claimed, the input signal is a time-dependent acoustic signal, and the signal domain is an interval of time. It was found that the method of the invention is particularly suitable for filling audio signal gaps of "medium length", i.e. lengths in the range of <NUM> to <NUM> milliseconds (ms). Evidently, also short gaps (<<NUM>) will be filled efficiently. Generally, the extent of the signal domain will be finite because of practical reasons.

In a second class of the invention as claimed, the input signal is an image signal on a spatial domain, in particular a region defined in a two-dimensional plane, wherein the input signal is an image signal. Another example is a three-dimensional image signal (e.g. for an object of an archaeological reconstruction or of medical imaging) for which the signal domain is a region in two- or three-dimensional space.

The mentioned input signal may additionally be time-dependent, such as a time-dependent image (such as a movie or clip or video), defined on a signal domain which is a product of a time interval and a two- or three-dimensional region, or a subset of this product.

In a further development of the invention, in determining an invertible matrix which optimizes a vector norm of a transformed coefficient vector, a vector norm may be used which is the sum of the moduli of the components of the respective vector (ℓ<NUM>-norm). In a variant, any other variant of the ℓ<NUM>-norm, which is calculated from the moduli of the components of the respective vector, may be used; for instance, with respect to a segmentation of the vector into a number of sub-vectors or cohorts (see below for a specific example), an ℓ<NUM>,∞-norm of the vector, which computes the maximum of the ℓ<NUM>-norms of the sub-vectors in the vector, or an ℓ<NUM>,<NUM>-norm of the vector, which computes the sum-of-squares (i.e. ℓ<NUM>-norm) of the ℓ<NUM>-norms of the sub-vectors over the vector, or another ℓ<NUM>,k-norm. Other suitable vector norms may be the ℓ<NUM>-norm (i.e. the number of vector components which are non-zero) or the ℓ<NUM>,∞-norm (the maximum number of the number of non-zero components in each of the sub-vectors, see eq. (<NUM>) below).

In many advantageous implementations of the method of the invention, the deformation matrix is a unitary matrix, and the invertible matrix used in determining a matrix which optimizes a vector norm of a transformed coefficient vector is also a unitary matrix.

During the calculation of the analysis coefficients for the input signal with respect to the initial frame, it is advantageous to exclude the signal values which are associated with locations within the gap. This can be implemented, for instance, by skipping the respective signal values during the computations (in particular, skipping the values during the summation operations) or by setting the respective signal values to zero in advance.

Advantageously, the initial dictionary signals may be realized as signal functions localized around respective center locations with an associated window length. A common example for this type of signal functions are Gabor frames using windows, such as Von-Hann windows, with a defined window length wg. Another example is a Gauss-type window, where the window length is identified with a suitably defined multiple of the variance σ, e.g. wg = 2σ, or wg = 4σ, or wg = <NUM> (i.e. the full width at half maximum of the Gauss function).

Furthermore, in order to improve avoiding artefacts from the "edge" of the gap, it may be advantageous to further exclude certain dictionary signals (or the corresponding analysis coefficients) and/or more signal values. For instance, also signal values associated with locations in a margin region around the gap may be excluded. This margin region may e.g. comprise guard intervals at the boundaries of the gap, the guard intervals having a width correlating with the window length wg of the initial dictionary signals, namely, a value selected from the range between <NUM> wg and <NUM> wg (inclusive), which includes the particular values of one window length and half the window length.

Furthermore, it may be advantageous to exclude certain initial dictionary signals, namely, those which would "pick up" significant contributions from the signal in the gap, which may be done alternatively or (preferably) in combination with the mentioned exclusion of signal values. More specifically, analysis coefficients may be excluded which correspond to initial dictionary signals which would contain contributions from signal values associated with locations in the gap, preferably locations in the gap plus locations in a margin region around the gap, wherein the margin region comprises guard intervals at the boundaries of the gap, the guard intervals having a width correlating with the window length wg of the initial dictionary signals, namely, a value selected from the range between <NUM> wg and <NUM> wg (inclusive), which includes the particular values of one window length and half the window length. The range of such a contribution may be defined with regard to the window length wg as mentioned above in the context of localized signal functions.

It is also remarked that the coefficient vector may have various representations: beside a straightforward vector representation, the coefficient vector may have a representation in a matrix-like array ("matrix vector") according to a segmentation into multiple vector segments, where each segment is transformed (deformed) in the same way by the deformation matrix (thus, the transformation operation may be represented as a matrix multiplication of the deformation matrix with the matrix vector). This latter representation is particularly useful where the deformation matrix has a block structure. A specific example of this is the vector representation according to the known real DGT [<NUM>, <NUM>]. The norm of the coefficient vector in the form of a matrix vector can be any suitable norm that is applied to a matrix or to a vector (when the matrix vector is re-ordered into a vector). The segmentation typically accords to a structure in the analysis domain, for instance such that each segment corresponds to a specific localization in the signal domain, grouping together the various dictionary signals associated with different frequencies or reciprocal space channels at that specific localization.

According to a further aspect of the invention, the deformation matrix may have a block structure where the initial frame is divided into "cohorts", preferably a block structure which reflects how the dictionary signals are located in the signal domain. Thus, it may be advantageous to have the initial frame decomposable into a number of cohorts, each cohort containing a plurality of dictionary signals and the cohorts being mutually disjoint (thus, no two cohorts have a dictionary signal in common). For instance, a cohort may correspond to a respective index in the time dimension (or more generally, a location in the signal domain); and the linear transformation according to the invertible matrix (Û) operates on the initial frame (i.e. on the initial dictionary signals) in a manner preserving the cohort structure which defines the division of the initial frames into cohorts. In other words, the matrix only operates within each of the cohorts, but does not combine dictionary signals of different cohorts: this means that the linear combinations "mix" dictionary signals only within a cohort and there is no "mixing" between dictionary signals belonging to different cohorts. Furthermore, the linear transformation (the matrix) preferably operates on the number of cohorts in a manner mutually equivalent between cohorts. This offers the advantage that the matrix can be represented with a largely reduced number of rows/columns, since only the action with respect to one cohort needs to be represented in the matrix coefficients (since the matrix operates uniformly upon the cohorts). In this case, a suitable implementation may provide that the initial dictionary signals are segmented into a matrix vector as mentioned above, or more specifically, arranged according to a matrix wherein the cohorts correspond to columns, and the linear transformation is realized as a matrix operating on the matrix of the initial dictionary signals by matrix multiplication from the left.

Furthermore, each cohort may be associated with a respective location in the signal domain (such as time and/or two- or three-dimensional space). Thus, the initial dictionary signals of each cohort may be signal functions localized around a center location in the signal domain associated with the respective cohort.

Generally, the invention envisages filling (or inpainting) the gap using a sparse reconstruction method. Sparse reconstruction algorithms are well known, such as the known SParse Audio INpainter (SPAIN) method for inpainting of gaps in audio signals. According to another aspect of the invention, a modified sparse reconstruction algorithm may be used in order to achieve a more efficient usage of the increased signal sparsity and benefit from dictionary learning. This aspect suggests that for determining the invertible matrix (denoted U) which optimizes a vector norm of a transformed coefficient vector which is the result of transforming the analysis coefficient vector (denoted v) using the matrix (thus the transformed coefficient vector is Uv; the vector v may be in any suitable representation such as a conventional representation or the mentioned matrix-vector representation; it is reminded that in the vector v some of the analysis coefficients may have been excluded as mentioned above, and likewise for the corresponding matrix) an iterative method may be used according to the steps of.

and repeating steps (b) to (d) until the matrix level value falls below a predetermined lower limit value or a maximum number of iterations is reached. The matrix U thus obtained is then taken as deformation matrix Û.

In order to further reduce calculation expenses, the Hermitian matrix H may be restricted to a Hermitian matrix having at most d non-vanishing off-diagonals, where d is a predetermined positive integer, preferably <NUM> or <NUM>. Furthermore, the Taylor expansion in step (c) may simply be a first-order Taylor expansion, i.e. eiH → I+iH. This is particularly suitable in combination with a block-structure deformation matrix as mentioned above; however, it will be clear that such combination is optional and not required.

The step of defining the gap may be realized simply by accepting a definition of the gap as input to the algorithm, or by determining those signal values which are unreliable algorithmically. Examples of determining signal values as unreliable are.

The above-discussed methods of the invention may conveniently be coded into a computer program. Thus, the invention also includes a computer program product - such as a data storage device or an apparatus employing one of the methods of the invention based on a computer program stored on a memory of the apparatus - which contains computer program code enabling execution of any or all methods of the invention when the computer program code is executed on a computer.

Further suitable aspects and benefits of the invention will become evident from the below description of an embodiment of the invention.

In the following, in order to further demonstrate the present invention, illustrative and non-restrictive embodiments are discussed, as shown in the drawings, which show:.

The detailed discussion of exemplary embodiments of the invention given below is intended to illustrate the basic concepts and further advantageous developments of the invention. It will be evident to the person skilled in the art to freely combine several or all of the features discussed in relation to the embodiment and variants discussed here, as deemed suitable for a specific application of the invention. This also includes the features mentioned in the appended claims, which can be freely combined mutually and with features taken from the discussed embodiments. Throughout this disclosure, terms like "advantageous", "exemplary", or "preferable" indicate elements or dimensions which are particularly suitable (but not essential) to the invention or an embodiment thereof, and may be modified wherever deemed suitable by the skilled person, except where expressly required. It will be appreciated that the invention is not restricted to the exemplary embodiment, which is given for illustrative purpose, merely representing suitable ways to implement the invention. Furthermore, it will be clear that the reference numbers used in the claims are merely intended to improve the readability of the claims, but neither the reference numbers nor the way they are used in the claims are to be construed as limiting the scope of the claims.

<FIG> depicts a flow diagram showing the flow of information during a signal processing procedure according to a first embodiment of the invention, which may be implemented on a computer station by means of a computer program or a suitable combination of cooperating computer program modules.

In a first stage <NUM>, an input signal, denoted x, is provided. Generally, the input signal x is presented in the form of a plurality of signal values associated with respective locations within the signal domain. Thus, in the case where the signal domain is time or another linear field (such as the field of real values), the input signal can be interpreted as a sequence of signal values; for instance, an input signal which is a time-dependent signal may be described as x = {xn, n = <NUM>,. , N - <NUM>} where n is an index corresponding to a sequence of points in time, for instance tn = t<NUM>+nΔt. More generally, the input signal may be seen as an array of the signal values, wherein each signal value can be identified through the (indices of) locations associated with it. Furthermore, gap parameters GPAR are input, in order to determine the gap. This may be done at a stage <NUM> or, in a variant not shown, together with providing the input signal x. The gap parameters specify the gap as a portion (i.e. a subset) of the signal domain by identifying the location of the gap and other gap parameters such as size and/or shape of the gap. Usually, the gap is a simple interval in time or a range in multi-dimensional domain, such as a box-range, a ball, a disk or a range of ellipsoidal shape, but it may be more complex, such as a union of intervals or ranges.

At stage <NUM>, the (initial) dictionary is defined. The dictionary is also referred to as initial frame, and is a set g = {g(p)} of initial dictionary signals g(p), which are defined over the signal domain (or a subset thereof as explained below; p is a symbolic index numbering the dictionary signals). For instance, in relation to time-dependent signals, the initial frame may be chosen as a Gabor frame, which is the set of base functions of Short-Time Fourier Transformation (STFT), where the index p numbers through all combinations of time and frequency channels. Various types of Gabor frames may be used, such as frames using Von-Hann windows. The initial frame g is, as such, suitable for modeling the signal in the containing region by usual analysis, i.e. linear combination of the initial dictionary signals. However, rather than using the initial frame for filling/reconstructing the signal in the gap, a different frame is used which is constructed from the initial frame as explained below.

In most cases the initial dictionary signals g(p) are defined over the entire signal domain. It is remarked, however, that often this is not a necessary condition; rather it is sufficient that they are defined merely on a "containing region" within the signal domain having sufficient extent, provided it contains the gap plus a surrounding region.

In stage <NUM> a learning neighborhood <IMG> is determined. This can be done in a separate module DETNBH or within the computer program implementing the process. The learning neighborhood is usually selected as the containing region minus the gap. It can be advantageous that not only the gap itself is excluded from the neighborhood, but also an additional guard interval is excluded in order to avoid influences from insecure data at, or close to, the edge of the gap.

For instance, the guard interval may be chosen at each boundary of the gap as a stripe having a width which corresponds to a width or half-width of the window of the Gabor frame; more generally suitable values of the width will be in the range of [<NUM> wg , <NUM> wg]. The exclusion of the gap may also be done in the analysis domain, i.e. the space of dictionary functions, by excluding dictionary functions which would contain contributions from the gap, possibly the gap augmented by a guard interval. The parameters for determining the neighborhood <IMG> are derived from the gap parameters GPAR and/or entered as external input parameters NPAR.

In stage <NUM>, a modified dictionary, also referred to as deformed frame ĝ, is generated ("basis optimization"). The deformed frame is a set of modified dictionary signals, i.e. ĝ = {ĝ(p)}; the modified dictionary signals ĝ(p) are invertible linear combinations of the initial dictionary signals, and thus have the same support domain as the initial dictionary signals. The modified dictionary signals are determined using an invertible matrix, denoted Û in <FIG> and referred to as deformation matrix. The deformation matrix is determined as a matrix which increases the sparsity of the input signal with respect to the deformed frame as compared to the sparsity of the input signal with respect to the initial frame (it is noted that increasing the sparsity is the same as reducing the norm of the coefficient vector obtained from the input signal by analysis using the dictionary signals). Suitably, the deformation matrix may be determined by a sub-process which seeks a matrix that optimizes the sparsity of the input signal from a predetermined set of matrices (see below) and, optionally, imposing selectable constraints as explained further below. The basis optimization (i.e. determination of the matrix Û and the deformed frame ĝ) can be done in a separate module BASDEF or within the computer program implementing the process, for instance as a sub-process. Input parameters for the basis optimization, such as the maximal number of iterations in an iterative process, the allowed number of off-diagonals in the matrix, etc., may be entered as external input parameters BPAR. Generally, the deformed frame will be at least partially different from the initial frame; so at least some of the modified dictionary signals will be different from all initial dictionary signals.

Advantageously, the determination of the invertible matrix Û is carried out based on the signal values in the containing region; herein, it is usually advantageous to exclude values as mentioned above, i.e. the signal values lying within the neighborhood <IMG> are used. Thus, first the analysis coefficients for the input signal restricted to the neighborhood are calculated with respect to the initial frame, producing an analysis coefficient vector v = AG x |<IMG>, where the symbol AG stands for the analysis with respect to the initial frame g, and the symbol. |<IMG> indicates restricting the analysis to the neighborhood. The analysis coefficient vector may be in a "standard" vector representation or in another suitable representation, such as the "matrix vector" representation mentioned above according to a suitable segmentation of the analysis domain. Then, starting from the set of invertible matrices or a predetermined subset of the set of invertible matrices, such as the set of unitary matrices, or unitary matrices which are exponentials of a Hermitian having only a small number of non-zero off diagonals, and/or matrices obeying a block-structure, a matrix Û is searched which optimizes a sparsity measure: Preferably, this sparsity measure is determined through a vector norm of a vector obtained by analysis with respect to the frame deformed by the invertible matrix; this is equivalent to applying the invertible matrix to the analysis coefficient vector, i.e. vector Ûv, and determining the vector norm from this transformed coefficient vector Ûv. The vector norm may be suitable chosen from the available vector norms, for instance the ℓ<NUM>-norm ∥ Vv ∥<NUM> (i.e. the number of vector components which are non-zero), or the ℓ<NUM>-norm ∥ Vv ∥<NUM> (the sum of the moduli of the vector components), or the the ℓ<NUM>,∞-norm as defined below in eq. (<NUM>), etc. Once an invertible matrix Û is found that increases the sparsity or even optimizes the sparsity measure, this matrix Û is taken as the desired deformation matrix Û. The deformed frame ĝ is calculated using this deformation matrix Û.

In the next stage <NUM> the reconstruction/filling of the signal is carried out. This is done by reconstructing the signal within the gap from the signal values in the containing region outside the gap (symbolically, x|<IMG>) by means of a sparse reconstruction method using the deformed frame ĝ. The reconstruction method used in this stage <NUM> may be a known sparse reconstruction method, or a modified SPAIN algorithm as disclosed herein. The result is a reconstructed or remodeled output signal x̂ in a final stage <NUM>. The stage <NUM> can be done in a separate module SPSREC or within the computer program implementing the process. Beside parameters derivable from the dictionaries, additional input parameters for the sparse reconstruction, such as the maximal number of iterations, may be entered as external input parameters SPAR.

Optionally, it may be provided that also the set of feasible signals Γx is determined in an additional stage <NUM>, which could be realized in a separate module DETFES. The set of feasible signals Γx is the set of "compatible signals", i.e. any signals which coincide with the input signal in (at least) the neighborhood, and preferably in the containing region (or the signal domain) outside of the gap, or a suitable defined sub-set thereof. This set Γx may be used as additional input to the sparse reconstruction stage (module SPSREC).

In a variant, instead of defining the gap directly, the gap may be determined from the input signal x using a suitable gap definition algorithm; in this case the parameters GPAR may specify parameters used by the gap definition algorithm. Different strategies for identifying a gap may be applicable depending on the specific nature of the underlying gap generating mechanism.

One typical scenario is the mitigation of clipping artefacts in an audio signal. Clipping causes signal degradations due to saturation effects of the amplifier: if the current value (i.e. the instantaneous value) of an audio signal exceeds or falls below certain thresholds, which may be defined through the parameter GPAR, the signal is cut off and set to these thresholds. To obtain the gap, the signal values are compared to the thresholds. A location is identified as a gap element if its corresponding signal value is identical, or close to, one of the thresholds; the set of all locations determined in this way is denoted as gap, its complement comprises all locations declared as reliable.

Another scenario is the compensation of signal packet losses in a communication network (packet loss concealment). For transmission of a signal, such as an audio signal, over a communication network, the signal information is divided into consecutive data packets using well-known coding and decoding techniques. Various impairments during transmission may occur, which will cause checksum errors of individual packets; these packets, in turn, are discarded in the receiver (and are lost). Signal parts which correspond to such lost or discarded audio packets, or more exactly the corresponding locations, will be identified as gap elements; the set of all such elements constitutes the gap.

A further scenario is based on error correcting codes, which are e.g. used to protect signal transmission or signal storage. Error correcting codes usually have the capability to correct errors up to a certain number of failures; if more failures occur, an incorrect code word will be selected. In combination with inpainting and identifying a gap, in case the closest distance (e.g. closest Hamming distance) between the received data word and any code word is in the range of the error correction capability of the used code, the signal part corresponding to this data word is marked as unreliable. This decreases decoding error probabilities in the presence of many failures.

In the following the mathematical background and an implementation of the invention according to another embodiment is discussed in detail, where, as a typical exemplary case, an audio signal is considered for audio inpainting. This exemplary embodiment illustrates the approach according to the invention exploiting the time-frequency sparsity of audio signals. The adaption to other signal types, such as two- or three-dimensional images, is straightforward. The audio signal is provided as a time-domain signal x represented as a sequence of N signal values (also referred to as samples), numbered by a time index n running from <NUM> to N, i.e. n = <NUM>,. , N - <NUM>. Some of the signal values are missing (or unreliable), and the indices of these missing samples are known.

The following notations are used in the present disclosure: Roman letters A, B,. , and a, b,. designate matrices, vectors, and scalars, respectively. <MAT> denotes the largest integer ≤ a. With [ · ] N = [ · modN] we abbreviate the modulo-N operation due to circular indexing. The ith component of the vector u is ui-<NUM>; the element in ith row and jth column of the matrix A is Ai-<NUM>,j-<NUM>. The superscripts T and H and * stand for transposition, Hermitian transposition, and (element-wise) complex conjugation, respectively. The N × N identity matrix is denoted by IN; the M × N all zero matrix is denoted by <NUM>M×N. For a vector u, we write <MAT> for its Euclidean norm and supp(u) for its support. For a set <IMG>, card(<IMG>) denotes its cardinality, and <IMG>(·) is the indicator function on <IMG>, which is <NUM> if its argument is in <IMG> and ∞ otherwise. We use the notation <IMG> to indicate the column submatrix of A consisting of the columns indexed by <IMG>. Similarly, for <MAT> we denote by <IMG> the subvector in <MAT> consisting of the entries of x indexed by <IMG>. Furthermore, <MAT> and ∥ u ∥<NUM> = |u<NUM>| + |u<NUM>| + ··· + |uN-<NUM>| denote the ℓ<NUM>-norm and the ℓ<NUM>-norm of the vector u = [u<NUM>, u<NUM>, ··· , uN-<NUM>]T, respectively. The value <MAT> (with tr(·) denoting the trace of a matrix) is the Froebenius norm of the matrix A; and ∥ A ∥∞,∞ denotes the largest modulus of all entries of A. Re(·) and Im(·) are real and imaginary part of its argument, respectively.

The time-domain signal x is interpreted as <MAT>. The subset of the indices of missing samples is referred to as the gap. The rest of the indices / samples will be considered and called reliable.

It is natural to require that the recovered signal should maintain consistency with the reliable part to the most possible extent. This condition may described formally by introducing a (convex) set Γx as the set of all feasible signals <MAT> where <MAT> is the binary "reliable mask" projection operator keeping the signal samples corresponding to the reliable part, while it sets the others to zero.

This embodiment employs the known Gabor transform, also known as Short-Time Fourier Transform (STFT), for implementing the initial dictionary signals, from which the sparsified signals, i.e. modified dictionary signals, are constructed as explained below. As is known to the person skilled in the art, the Gabor transform uses a single prototype window function g which is translated in time and frequency-modulated [<NUM>, <NUM>]. The window serves to localize the signal in time. For the discrete Gabor transform (DGT), the translation of the window is characterized by the integer window shift a. The number of modulations (i.e. channels in the transformed "reciprocal space") is denoted by M and will be referred to as the number of frequency channels (in more general implementations, they may also be called wave number channels, reciprocal space channels, or the like). Often it is also required that the signal length N is divisible by a. Then, the system consists of P = MN/a Gabor vectors <MAT> , p = <NUM>,. , P - <NUM>. We will refer to these vectors as the Gabor atoms, and the whole system <MAT> with <MAT> and p = kM + m is referred to as the Gabor dictionary.

Note that the Gabor window g is usually identified with its shorter counterpart, i.e. the interval (time interval, more generally subrange in the signal domain) comprising only those elements of g that are within the smallest interval containing the support of g. The length of this interval is usually much smaller than N and is denoted as window length wg.

We also remark that the set G of Gabor atoms is overcomplete. In fact, there exist suitable combinations of g and the parameters a and M such that the resulting Gabor system forms a frame for <MAT>, i.e. any <MAT> can be represented in a stable way as a linear combination of the.

Gabor vectors [<NUM>, <NUM>]. Although Gabor bases can be constructed, they have undesired properties [<NUM>]. Therefore, overcomplete systems which allow non-unique signal representations are usually preferred. Also the signal is formally treated as a complex vector, although the underlying input signal (e.g. audio signal) is real-valued.

In frame theory, the so-called analysis operator <MAT> produces coefficients from the signal, whereas its adjoint D = AH, the synthesis operator <MAT>, generates a signal from the coefficients. Its composition S = DA is denoted as frame operator. In the context of Gabor frames, we use the subscript notation AG and DG for analysis and synthesis operators, respectively, to emphasize the Gabor structure.

In the majority of implementations of the invention the frames will correspond to the so-called painless case [<NUM>, <NUM>]. Such frames have convenient properties from both theoretical and practical perspective and are typically the ones considered in signal processing. Their analysis, synthesis, and frame operators satisfy the "painless condition" <MAT> with Sn,n > <NUM>, n = <NUM>,. , N - <NUM>, i.e. the frame operator matrix is a diagonal matrix with diagonal elements strictly larger than zero. (For notational simplicity, this disclosure does not distinguish between operators and their matrix representation with respect to the canonical basis throughout. ) Evidently, tight frames, i.e. frames where S = λIN, fall within the class of painless frames. For Gabor frames, wg ≤ M implies (<NUM>). In the following, frames realizing the painless case will be considered primarily, but the generalization to non-painless frames is straightforward.

This section gives a brief review of the SPAIN algorithm, as originally presented in [<NUM>]. The SPAIN algorithm comes in two variants: the first one exploits analysis sparsity, the second one exploits synthesis sparsity (both in the time-frequency domain). Note that the sparse signal processing literature relied for a long time on the so-called synthesis model, where one seeks for a small number of coefficients synthesizing the desired signal [<NUM>, <NUM>-<NUM>]. More recently, the analysis model appeared, where one looks directly for the signal, with the constraint that its coefficients after analysis are sparse [<NUM>, <NUM>]. Generally, however, we can only expect these coefficients to be approximately sparse, since underlying uncertainty principles restrict the maximum degree of sparsity [<NUM>, <NUM>, <NUM>]. If the synthesis/analysis operators are invertible, both approaches are equivalent. As said above, usually they are not invertible. Note also that the synthesis approach is less restrictive, as for redundant systems (i.e. P > N) different coefficients lead to the same signal.

The two variants of SPAIN aim at solving the following optimization tasks, Γx <MAT> <MAT> where (<NUM>) and (<NUM>) present the formulation referred to as the analysis and the synthesis variant, respectively. The abbreviation s. means "such that". In both cases, the minimizing y will be the reconstructed time-domain signal. It is obvious that a brute-force solution of (<NUM>) and (<NUM>) will be infeasible due to its huge computational complexity. As an alternative, SPAIN applies the Alternating Direction Method of Multipliers (ADMM) [<NUM>] adapted to the optimization of the above non-convex problems. The ADMM is able to solve problems of the form <MAT> where <MAT> is a linear operator, and f and g are real convex functions. A reformulation of (<NUM>) yields <MAT> from which an Augmented Lagrangian (in scaled form) is computed. Finally, this Lagrangian is minimized individually with respect to each involved variable in an iterative fashion until a sufficiently accurate solution is obtained. Note that in case of SPAIN, only an approximate solution can be expected, since (<NUM>) and (<NUM>) violate the convexity assumption.

Of special importance in SPAIN is also the segment-wise application of the algorithms, where first the time-domain signal is segmented using overlapping windows, second the algorithms are applied segment-wise using an overcomplete DFT frame, and finally the restored blocks are combined via an overlap-add scheme.

The inventors recognized that when eq. (<NUM>) and (<NUM>) were to be applied to the signal consisting of the gap plus the adjacent parts of length wg before and after the gap as a whole using a Gabor frame of appropriate dimension, the method would fail completely for medium length gaps. This is caused by the regularizer ∥·∥<NUM>, which penalizes the Gabor coefficients globally instead of locally. In particular, the reconstructed signal will be set to zero within the gap because its Gabor coefficients remain to be sparse globally. In order to cope with this behavior, we propose a modification of the original SPAIN algorithm that works without segmentation.

As mentioned above, Gabor systems impose a time-frequency structure. In case of the DGT, the P = MN/a Gabor analysis coefficients can be rearranged into an M × (N/a) matrix (it is assumed that a divides N), whose column and row indices correspond to discrete time and discrete frequency, respectively; this is a specific example of a "matrix vector" representation. Mathematically, we express this operation via the (invertible) mapping <MAT> to denote the mapping between a P-dimensional vector z and a matrix τ(z) having dimension M × (N/a); it will be appreciated that this rearrangement does neither affect the total number nor the values of the coefficients. Note further that for any audio signal <MAT> the complex-valued matrix τ(AGx) is conjugate-symmetric with respect to the frequency/row index. Hence, we do not loose any information if we only keep the rows corresponding to the first <MAT> frequency indices; the remaining rows can be easily reobtained by exploiting conjugate-symmetry. Hereinafter, this restriction operation is expressed via the additional mapping <MAT>; its reversed operation, i.e. the conjugate-symmetric extension, is denoted by <MAT>. For notational simplicity, we use the abbreviation η for the composite mappings, i.e. <MAT> and <MAT>. The mapping <MAT>, x ↦ η(AGx) is often referred to as real DGT [<NUM>, <NUM>].

For each <MAT> let us define the ℓ<NUM>,∞-norm according to: <MAT>.

This ℓ<NUM>,∞-norm uses the matrix vector representation (where each sub-vector xi represents a cohort) and computes the maximum number among the number of non-zero entries in each of the sub-vectors. Instead of (<NUM>), we propose to solve the following two optimization tasks, <MAT> <MAT> where, again, (<NUM>) and (<NUM>) present the formulation referred to as the analysis and the synthesis variant, respectively. In the context of the invention the ADMM technique is adapted to facilitate solving the problems of (<NUM>) and (<NUM>). The choice of ADMM is driven by two main reasons; scalability and fast convergence, which implies that it can be used for large scale problems.

We will first consider the analysis variant (<NUM>). Let us fix a sparsity parameter k and let <MAT>.

In order to apply the ADMM algorithm, we rewrite (<NUM>) with ε = <NUM> as <MAT>.

The Augmented Lagrangian is given as <MAT> leading to the Augmented Lagrangian in scaled form [<NUM>], <MAT>.

The update rules of ADMM now yield, <MAT> <MAT> <MAT>.

Note that the conjugate-symmetry is preserved in each of these steps, provided that the initialization vectors y(<NUM>) and r(<NUM>) are real-valued and conjugate-symmetric, respectively.

The B-update (<NUM>) is solved exactly (cf. <NUM>]) by means of a time-frequency hard-thresholding operator <MAT> defined according to <MAT> where <IMG>(·) is essentially (note: in order to correctly deal with the conjugate-symmetry, the modulus has to be scaled for at most two elements of the vector; i.e. at frequency <NUM> and, for even M, also at frequency M/<NUM>) the conventional hard-thresholding operator for vectors, which preserves the k elements with largest modulus and sets everything else to zero. More specifically, the B-update (<NUM>) is given by <MAT>.

In fact, the y-update (<NUM>) is equivalent to applying the inverse DGT using the canonical dual window to η†(B(i+<NUM>)) - r(i) and then projecting it onto the set of feasible solutions Γx.

The overall algorithm is initialized with a certain pre-defined sparsity parameter k = s, that will be increased in every tth iteration by s, until the stopping criterion is met. A pseudocode reflecting the A-SPAIN algorithm is given as "Algorithm <NUM>" in Table <NUM>; all Tables are given at the end of the description. Note that usually s = t = <NUM>; however, as in the original SPAIN algorithm, we allow for more flexible settings in order to speed up the algorithm (s > <NUM>) or supporting its convergence (t > <NUM>). The input object <MAT> denotes the synthesis operator of the canonical dual Gabor frame [<NUM>]. In practice, of course, only the windows for AG and <MAT> as well as the parameters a and M are passed to the algorithm, and not the full matrices.

Next, considering the synthesis variant (<NUM>), let us fix a sparsity parameter k and rewrite (<NUM>) with ε = <NUM> as <MAT>.

The update rules of ADMM now yield <MAT> <MAT> <MAT>.

Here, the conjugate-symmetry is preserved in each of these steps, provided that both initialization vectors y(<NUM>) and r(<NUM>) are real-valued.

The B-update (<NUM>) is a challenging task to solve; in fact, it is a sparse reconstruction problem utilizing the ℓ<NUM>,∞-norm. However, we will rely on the convenient ADMM behavior that it still converges even when the individual steps are computed only approximately. Adapting the idea in [<NUM>] to our setting, we suggest to apply the time-frequency hard-thresholding operator <MAT> to the analysis coefficients of y(i) - r(i) with respect to the canonical dual frame, or, more formally, <MAT>,.

Here, <MAT> denotes the analysis operator of the canonical dual Gabor frame. Again, the overall algorithm is initialized with a certain pre-defined sparsity parameter k = s, that will be increased in every tth iteration by s, until the stopping criterion is met. A pseudocode reflecting the S-SPAIN algorithm is given as "Algorithm <NUM>" in Table <NUM>.

As pointed out above, all SPAIN variants (including both original and modified versions) rely on the observation that real-world audio signals have approximately sparse representations with respect to Gabor dictionaries. One may ask whether other dictionaries admit even more sparse signal representations and, in turn, improve the audio inpainting capabilities if applied instead of a Gabor dictionary. Clearly, it is not trivial at all to come up with a dictionary which yields better results than the Gabor dictionary.

The invention proposes to learn an optimized dictionary from the reliable signal parts around the gap, aiming at obtaining a more sparse representation within the gap as well, offering that the inpainting performance is enhanced. Of course, this is based on the assumption that the optimum sparsifying dictionary does not change too fast over time. Furthermore, we need some additional reliable signal parts in the neighborhood of the gap to learn the dictionary, which, in case we deal with more than one gap, requires some distance between adjacent gaps. In order to keep the learning effort low, we will not learn it from scratch; instead the invention proposes to "deform" the initial dictionary, i.e. in the present embodiment, a Gabor dictionary.

The starting point is the modified SPAIN algorithm as discussed above. Motivated by (<NUM>), we aim at finding a frame satisfying the painless condition (<NUM>) with analysis operator A such that for the degraded signal x, Ax is as sparse as possible with respect to the ℓ<NUM>,∞-norm in a neighborhood of the gap. More formally, we intend to solve <MAT> where the notation | <IMG> represents the restriction to the neighborhood. Let AG denote the analysis operator of a given Gabor frame satisfying the painless condition (<NUM>). (The analysis vector η(AGx) corresponds to the vector v that was mentioned in the introductory part as denoting the analysis coefficient vector. ) We next "deform" the Gabor frame using a unitary "deformation" operator <MAT> by defining a modified analysis operator <MAT>.

Note that we restrict to unitary deformations, since these preserve the painless condition (<NUM>) due to <MAT>. Furthermore, we are only interested in deformation operators which increase sparsity in frequency direction and which preserve the conjugate-symmetry of the DGT. Therefore, we additionally impose the following structure on W, <MAT> where <MAT> is a unitary matrix (i.e. dimension M × M). It is clear that (<NUM>) represents a well-defined unitary operator. It will be appreciated that (<NUM>) contains a simple matrix product of the matrices V and τ(z) (thus largely reducing the complexity of operation as compared to a general P × P-matrix W acting on a P-dimensional vector z). This operation corresponds to a transformation which only combines ("mixes") the coefficients of the individual columns while excluding a mixing between columns. Moreover all columns are treated in the same manner, i.e. the mixing within each column is the same for all columns. Thus, (<NUM>) denotes a special realization of a matrix block structure where the matrix represents a transformation which preserves the cohort structure.

It is remarked that, as mentioned above, in this embodiment the column index corresponds to the discrete time index; more generally, the cohorts may correspond to a location in the signal domain, while the index within a cohort (i.e. the row index) relates to modulations (frequency/reciprocal space). In other implementations, depending on the respective implementation, the block structure may be finer (e.g. treating segments within a column as cohorts) or coarser (e.g. combining defined groups of columns into cohorts) or more complex (e.g. combining segments of columns, etc.).

Regarding the structure of V, we have to distinguish between M even and M odd. For M even, let <MAT> whereas for M odd, <MAT> with unitary <MAT> and flipping matrix <MAT>.

It is easy to see that for a conjugate-symmetric z, <MAT> with <MAT> respectively. Combining (<NUM>), (<NUM>), and (<NUM>), we obtain <MAT> where <IMG> denotes the set of all unitary matrices Ue/o with the structure described in (<NUM>).

However, the inventors noted that this is a highly non-convex problem, so that we can only hope to find an approximate solution. To this end, it is suggested that the ℓ<NUM>-norms be relaxed to ℓ<NUM>-norms [<NUM>, <NUM>], and furthermore it may be suitable (although not necessary) to replace the max-operation in the definition of the ℓ<NUM>,∞-norm by a summation or other calculation, which is convenient in terms of calculation speed. An advantageous choice is replacing the norm in favor of an ℓ<NUM>-norm or a variant thereof, such as ℓ<NUM>,<NUM>-norm or ℓ<NUM>,∞-norm, where ℓ<NUM>,k-norm stands for the ℓk-norm evaluated for the M'-dimensional vector formed from the norms of the sub-vectors xi; thus, in the following we use an ℓ<NUM>,<NUM>-norm (which is equivalent to the ℓ<NUM>-norm on the entire vector). This results in the following optimization problem: <MAT>.

Finally, it is worthwhile to note that the minimization (<NUM>) is effectively carried out over the set <IMG> of unitary <MAT> matrices U according to (<NUM>), and that the optimum Û is given by <MAT> where <MAT> <MAT>.

Suppose we have found an optimized unitary matrix Û (we will present an algorithm forV̂his task in the next subsection). Inserting it into (<NUM>) and (<NUM>)/(<NUM>) yields Ûe/o and e/o, respectively. From this it is also possible to construct the optimized matrix W of (<NUM>); but in practice this is cumbersome because of the large dimension of W. Therefore, we use a different approach: via (<NUM>) and (<NUM>), we obtain the analysis operator of the deformed frame as, <MAT>.

Likewise, the synthesis operator of its canonical dual frame is given by <MAT>.

The above enables to re-formulate Algorithm <NUM> and Algorithm <NUM> in a way such that the learned sparsity-optimized frame is used instead of the Gabor frame. More specifically, for the analysis variant we replace AG with Â and <MAT> with Ďcd. In the algorithm, the (optimized) matrix Ûe/o is required as an additional input. This algorithm was named A-SPAIN-LEARNED, and a pseudocode reflecting this algorithm is given as "Algorithm <NUM>" in Table <NUM>. A Matlab code program implementing the Algorithm <NUM> is given in Table <NUM>. The Algorithm <NUM> contains calls to functions from the open-source toolbox LTFAT [<NUM>].

For the synthesis variant, an approach is used which optimizes the analysis operator of the canonical dual frame. Then, the synthesis operator of the corresponding original frame is expected to represent any signal in the feasible set Γx with fewer coefficients than the synthesis operator of the given Gabor frame. Again, the algorithm requires the specification of the (optimized) matrix Ûe/o as additional input. This algorithm was named S-SPAIN-LEARNED, and a pseudocode reflecting this algorithm is given as "Algorithm <NUM>" in Table <NUM>. A Matlab code program implementing the Algorithm <NUM> is given in Table <NUM>. The Algorithm <NUM> contains calls to functions from the open-source toolbox LTFAT [<NUM>].

Equation (<NUM>) codes the basis optimization. In order to solve (<NUM>), a basis optimization technique originally developed in a different context of channel estimation [<NUM>, <NUM>] is adapted for the present signal processing tasks. To simplify notation we denote <MAT> and define <MAT> Because the minimization problem (<NUM>) is non-convex (since <IMG> is not a convex set), standard convex optimization techniques will be problematic or cannot be used at all. We therefore propose an approximate iterative algorithm that relies on the following facts: (<NUM>) Every unitary M̃ × M̃ matrix U can be represented in terms of a Hermitian M̃ × M̃ matrix H as U = eiH; and (<NUM>) the matrix exponential U = eiH can be approximated by a Taylor expansion (with respect to H), for instance its first-order Taylor expansion: U ≈ I + iH (in this section we use the shorthand I for the M̃ × M̃ identity matrix IM̃).

Even though U is unitary and I + iH is not, this approximation will be good if ∥ H ∥∞,∞ is small. Because of this condition, it will be advantageous to construct U iteratively: starting with the identity matrix, we perform a small update at each iteration, using the approximation U ≈ I + iH in the optimization criterion, whereas the complete exponential eiH is used for updating U (thus, the iterated U is always unitary). More specifically, at the rth iteration, we consider the following update of the unitary matrix U(r): <MAT> where H(r) is a small Hermitian matrix that remains to be optimized. Note that U(r+<NUM>) is again unitary because both U(r) and eiH(r) are unitary. Ideally, it would be preferred to optimize H(r) according to (<NUM>), i.e. by minimizing <MAT>.

Since this problem is still non-convex, we use the approximation eiH ≈ I + iH (first-order Taylor expansion), and thus the final minimization problem at the rth iteration is <MAT>.

Here, <IMG> denotes the set of all Hermitian M̃ × M̃ matrices H that are small in the sense that ∥ H ∥∞,∞≤ ρr, where ρr is a positive constraint level (a small pr ensures a good accuracy of our approximation U ≈ I + iH and also that eiH(r) is close to I). The problem (<NUM>) is convex and thus can be solved by standard convex optimization techniques [<NUM>].

The next step at the rth iteration is to test whether the cost function is smaller for the new unitary matrix eiH(r) U(r), i.e. whether <MAT>.

In the positive case, we actually perform the update of U(r) and we retain the constraint level ρr for the next iteration: <MAT>.

Otherwise, we reject the update of U(r) and reduce the constraint level ρr: <MAT>.

This construction guarantees the cost function sequence <IMG> ∥ U(r) (ηe/o(AGx))q ∥<NUM>, r = <NUM>,<NUM>,. to be monotonically decreasing.

The above-described iteration process is terminated if ρr falls below a prescribed threshold or if the number of iterations exceeds a certain value. The iteration process is initialized by the identity matrix I (i.e. IM̃), because the "undeformed" Gabor dictionary is known to yield relatively sparse analysis coefficients. We note that efficient algorithms for computing the matrix exponentials eiH(r) are well-known (see e.g. [<NUM>]).

Finally, we would like to emphasize that if we restrict the minimization of (<NUM>) to the set <IMG> of small Hermitian matrices with at most d non-vanishing off-diagonals, we obtain reduced computational complexity, at the costs of a potentially suboptimum variant of the algorithm. In some cases, the choice of such a modification will be mandatory, especially if we deal with large-scale problems. The overall basis optimization algorithm (which also uses an off-diagonal parameter d as input) is represented as a pseudocode as "Algorithm <NUM>" in Table <NUM>. A Matlab code program implementing the Algorithm <NUM> is given in Table <NUM>.

All algorithms <NUM> to <NUM> were implemented and executed in the known Matlab programming language, with Algorithm <NUM> additionally using the convex optimization package CVX, which is a Matlab-based modeling language destined for solving convex programming [<NUM>]. Table <NUM>, Table <NUM> and Table <NUM> show the respective program codes. Table <NUM> shows the program code of a subroutine called by the codes of Table <NUM> and Table <NUM>. Algorithms <NUM> and <NUM> also use functions in the known open-source toolbox LTFAT [<NUM>].

This section deals with how to select the neighborhood <IMG> of the gap, which specifies the (input) training matrix (ηe/o (AGx))|<IMG> for Algorithm <NUM>. It seems to be obvious that it should be as close as possible to the gap since this increases the chance that the learned sparse representation is also a sparse representation for the signal part within the gap. On the other hand, we also expect small "guard" intervals between the gap and training borders to be mandatory to add, in order to avoid that unreliable data from the gap influences our training data. According to the underlying Gabor structure, this guard interval should be at least of length wg (length of Gabor window). Regarding the size of the neighborhood we have to somehow rely on intuition. A larger neighborhood yields a larger training matrix (η(AGx))|<IMG> and, probably, more accurate training results but it also increases the likelihood that signal variations occur within the neighborhood. Such signal variations could yield a dictionary adapted to signal parts which are essentially independent of the signal within the gap. Suppose the signal <MAT> has a gap at the indices n = nB, nB + <NUM>,. Taking into account the guard intervals, the neighborhood is given by <IMG> = <IMG> U <IMG> with <MAT> <MAT> i.e. a part before and a part after the gap. Its overall length is <IMG>.

In the following and with reference to <FIG> and <FIG>, the dictionary learning method according to the invention was tested using a numerical evaluation of an example input signal. In the <FIG> the analysis coefficients of a respective analysis calculation of a signal are depicted as functions of the time and frequency, wherein the shades indicate the moduli (i.e. absolute values) of the coefficients, corresponding to the scale to the right of each figure.

A conventional setup was used which corresponds to the setup of O. Mokry et al. As main performance measure, the signal-to-noise ratio (SNR) is used here, which is defined as <MAT> where xorig and xinp denote original and inpainted signal within the gaps, respectively. Clearly, higher snr values reflect better reconstruction. The average SNR is computed by first computing the particular values of SNR in dB, and then taking the average.

As another quantity for describing the quality reconstruction, the PEMO-Q criterion [<NUM>] may be used, which takes a model of the human auditory system into account, and therefore, is closer to the subjective evaluation than the SNR. The PEMO-Q criterion defines a measured quantity called objective difference grade (ODG), which can be interpreted as the degree of perceptual similarity of xorig and xinp. The ODG attains values from -<NUM> (minus <NUM>, very annoying) up to <NUM> (imperceptible), expressing the effect of audio artifacts in the restored signal.

In the experiments made by the inventors, a collection of ten music recordings sampled at <NUM> were used, which exhibit different levels of sparsity with respect to the original Gabor representation. The signals were obtained from the known EBU SQAM dataset [<NUM>]. In each test instance, the input was a signal containing <NUM> gaps at random positions. The lengths of the gaps ranged from <NUM> up to <NUM>. For fixed lengths, the results over all ten signals consisting of the <NUM> gaps were averaged.

Throughout these experiments a tight Gabor frame was used with the Hann window of length wg = <NUM> samples (approximately <NUM>), window shift a = <NUM> samples and with M = <NUM> frequency channels. We used the fast implementation of Gabor transforms offered by the LTFAT toolbox [<NUM>, <NUM>, <NUM>] in our computations, and we adopted its time-frequency conventions. The calculations were carried out according to the Algorithms <NUM> to <NUM> and the Matlab codes illustrated in the Tables. The code of Algorithm <NUM> is coded and executed in CVX, which is a Matlab-based modeling language destined for solving convex programming [<NUM>]. For the basis optimization algorithm, we set the maximum number of iterations to rmax = <NUM>, the off-diagonal parameter to d = <NUM>, and the remaining parameters to ρstart = <NUM> and ε = <NUM>-<NUM>. Finally, all SPAIN variants used the input parameters s = t = <NUM> and ε = <NUM>.

Referring to <FIG>, one sequence of calculations had the purpose to analyze the sparsity of the analysis coefficients of real-world audio signals with respect to the optimized frame and to compare it with the analysis coefficients using the original Gabor frame. More specifically, <FIG> illustrate analysis coefficients calculated from the signal "a60_piano_schubert" from the EBU SQAM dataset [<NUM>], into which a gap of length <NUM> was introduced by replacing the signal by zero values. <FIG> relate to the signal without gap, <FIG> to the signal with the gap having a length of <NUM>. In <FIG>, the Gabor dictionary was used, in <FIG> the learned dictionary according to the invention. The coefficients within the green rectangles are used for training, corresponding to a neighborhood length parameter of <IMG> = <NUM>wg/a = <NUM>. As expected, the representations with respect to the learned dictionary are more sparse than the ones with respect to the Gabor dictionary. Although the training coefficients have a certain distance from the gap area, the learned dictionary also sparsifies the original signal within the gap area. This shows that the sparsifying dictionary does not change too fast over time.

Referring to <FIG>, another sequence of calculations served to investigate how the size of the learning neighborhood <IMG> around the gap impacts the sparsity of the learned representation. To this end, <FIG> depicts the analysis coefficients of the signal "a25_harp" from the EBU SQAM dataset [<NUM>]. <FIG> show coefficients obtained using Gabor dictionary, <FIG> learned dictionaries. The coefficients within the green rectangles are used for training: the coefficients of <FIG> are obtained from training according to <FIG>, and the coefficients of <FIG> are obtained from training according to <FIG>. Note that the neighborhood length parameter used to obtain <FIG> was <IMG> = <NUM>wg/a = <NUM>, whereas the neighborhood length parameter used to obtain <FIG> was <IMG> = 2wg/a = <NUM>. It is seen that the analysis coefficients with respect to the dictionary obtained by a larger training set are more sparse, especially for the coefficients in the training areas.

Further evaluations were made to assess how the size of the training neighborhood impacts the inpainting performance. To that end, we considered the same two cases as in the previous subsection, i.e. two different neighborhoods, one with length parameter <IMG> = <NUM>wg/a = <NUM> and the other one with length parameter <IMG> = <NUM>wg/a = <NUM>. It was found that a version with smaller learning neighborhood outperforms a corresponding version with larger neighborhood for almost all gap lengths in terms of SNR and for all gap length in terms of ODG. This behavior is observed for both analysis and synthesis variants. It seems that for the considered scenarios, the benefit of only using learning coefficients close to the gap outweighs potential other disadvantages of small training sets. Hence, for all further comparisons, we restricted to the training neighborhood with length parameter <IMG> = <NUM>wg/a = <NUM>. Furthermore, it was found that the method according to the invention produced distinctly better outcomes of restored signals, as compared to various other existing audio inpainting methods, in particular in terms of ODG values, over the whole gap length range (<NUM> to <NUM>).

The experimental results demonstrated that the proposed dictionary learning approach yields large performance gains, in particular in combination with the modified analysis SPAIN variant, and significantly outperforms known inpainting approaches. It is emphasized that the presented dictionary learning framework will be applicable to other sparsity-based inpainting techniques apart from SPAIN as well. Furthermore, it is expected that suitable adaptations to e.g. (re-)weighted ℓ<NUM>-minimization methods will also increase their performance. Moreover, it is emphasized that it is possible to replace the Gabor frame by any other frame (e.g. non-stationary Gabor frames) satisfying the painless condition as long as this frame offers a localization in the time-frequency structure, or more generally, in the location space/reciprocal space structure for other types of input signals.

The invention is not restricted to signals on a one-dimensional domain, but can also be used for filling or reconstructing of signals on a two- or higher-dimensional domain with gaps (such as a "hole" in an image). The inpainting of signals "living" in a two- or more-dimensional domain according to the invention is made according to the same basic concepts as presented above, and in the following discussion of a further embodiment, only the differences to the above embodiment for audio inpainting are discussed to clarify suitable adaptations. Thus, in the following the mathematical background and an implementation of the invention according to a further embodiment is discussed in detail, where, as a typical exemplary case, a grayscale image signal is considered for image inpainting. The generalization to more than two dimensions, e.g. processing of a three-dimensional image signal, will be evident to the skilled person from the following discussion. Also, the extension to colored images (or other types of vector information) is straightforward by considering each color component (or vector component) separately. We assume that the grayscale image signal is provided as a two-dimensional (2D) signal x represented as a collection of N<NUM>N<NUM> signal values (also referred to as samples), numbered by a location (position) index (n<NUM>, n<NUM>) where n<NUM> = <NUM>,. , N<NUM> - <NUM> and n<NUM> = <NUM>,. , N<NUM> - <NUM>. Each signal value represents the (light) intensity of an image pixel at a certain position using, e.g., the convention <NUM> ≤ x ≤ <NUM> (<NUM>. black; <NUM>. Furthermore, the dynamic range is usually limited, so that only a finite number of quantized values are considered within this <NUM>-<NUM>-range. We assume that some of the signal values are missing (or unreliable), and the position indices of these missing samples are known. The subset of the indices of missing samples is referred to as the gap. Typically, the gap is a 2D sub-region within {<NUM>,. , N<NUM> - <NUM>} × {<NUM>,. , N<NUM> - <NUM>}, of, e.g., circular or rectangular shape. However, more general shapes are also possible.

It is natural to require that the recovered signal should maintain consistency with the reliable part to the most possible extent. This condition may be described formally by introducing a set Γx as the set of all feasible signals <MAT> where MR: [<NUM>,<NUM>]N<NUM>×N<NUM> → [<NUM>,<NUM>]N<NUM>×N<NUM> is the binary "reliable mask" projection operator keeping the signal samples corresponding to the reliable part, while it sets the others to zero.

This embodiment employs the known 2D Gabor transform, for implementing the initial dictionary signals, from which the modified dictionary signals, are constructed. As is known to the person skilled in the art, the 2D Gabor transform uses a single prototype 2D window function g which is 2D-translated in location and 2D-frequency-modulated [<NUM>, <NUM>]. The window serves to localize the signal with respect to location (position). For the discrete 2D-Gabor transform (DGT), the translation of the window is characterized by the integer window shifts a<NUM> and a<NUM>. The total number of 2D-modulations (i.e. channels in the transformed "reciprocal space") is given by M<NUM>M<NUM>, where M<NUM> and M<NUM> denote the number of individual 1D-modulations in the respective subdirections, and will be referred to as the number of frequency channels (in more general implementations, they may also be called wave number channels, reciprocal space channels, or the like). Often it is also required that the N<NUM> and N<NUM> are divisible by a<NUM> and a<NUM>, respectively. Then, the system consists of P = M<NUM>M<NUM>N<NUM>N<NUM>/(a<NUM>a<NUM>) Gabor matrices <MAT>, p = <NUM>,. , P - <NUM>. We will refer to these matrices as the Gabor atoms, and the whole system <MAT> <MAT> with <MAT> is referred to as the Gabor dictionary.

Often, the prototype 2D window function g is chosen as product, i.e., gn<NUM>,n<NUM> = (g<NUM>)n<NUM> (g<NUM>)n<NUM>, of individual 1D window functions g<NUM> and g<NUM>, such as 1D Von-Hann windows and 1D Gauss-type windows (as used in the audio inpainting example described above).

We also remark that the set G of Gabor atoms is overcomplete. In fact, there exist suitable combinations of g and the parameters a<NUM>, a<NUM>, M<NUM>, and M<NUM> such that the resulting Gabor system forms a frame for <MAT>, i.e., any <MAT> can be represented in a stable way as a linear combination of the Gabor matrices [<NUM>, <NUM>]. Note that the signal is formally treated as a complex matrix, although the underlying input signal x (e.g. grayscale image) is real-valued with <NUM> ≤ x ≤ <NUM>. Let us denote by <MAT> with <MAT>, where indN<NUM>×N<NUM>: {<NUM>,. , N<NUM> - <NUM>} × {<NUM>,. , N<NUM> - <NUM>} → N<NUM>N<NUM>, (n<NUM>, n<NUM>) ↦ n<NUM>N<NUM> + n<NUM>, the row-wise stacking operator (any other stacking operator, e.g., column-wise stacking operator, would be equally suitable, as long as different stacking operators are not mixed). It follows that the system <MAT> forms a frame for <MAT> and any vectorized grayscale image <MAT> can be represented in a stable way as a linear combination of the vectors in Gvec. Moreover, we can associate an analysis operator <MAT>, which produces coefficients from vectorized image signals, and its adjoint D = AH, the synthesis operator <MAT>, which generates vectorized 2D signals from the coefficients. Its composition S = DA is denoted as frame operator. Again, we use the subscript notation AG and DG for analysis and synthesis operators, respect-tively, to emphasize the underlying 2D Gabor structure. In the majority of implementations of the invention the frames will correspond to the so-called painless case [<NUM>, <NUM>], for which the frame operator matrix S is a diagonal matrix with diagonal elements strictly larger than zero.

Also in this embodiment, a segmentation into cohorts may be used. We partition the frame consisting of the P = M<NUM>M<NUM>N<NUM>N<NUM>/(a<NUM>a<NUM>) frame vectors vecN<NUM>×N<NUM>{g(p)} into N<NUM>N<NUM>/(a<NUM>a<NUM>) cohorts C(k<NUM>,k<NUM>) = {vecN<NUM>×N<NUM>{(k<NUM>,m<NUM>,k<NUM>,m<NUM>)}: m<NUM> = <NUM>,. , M<NUM> - <NUM>, m<NUM> = <NUM>,. , M<NUM> - <NUM>} indexed by k<NUM> = <NUM>,. , (N<NUM>/a<NUM>) - <NUM>, k<NUM> = <NUM>,. , (N<NUM>/a<NUM>) - <NUM>. Note that each cohort C(k<NUM>,k<NUM>) consists of M<NUM>M<NUM> frame vectors which are all "localized" around the vectorized location (position) indN<NUM>×N<NUM>(a<NUM>k<NUM>, a<NUM>k<NUM>). (The matrices in <MAT> are localized around (a<NUM>k<NUM>, a<NUM>k<NUM>).

It follows that the P analysis coefficients of a vectorized image can be rearranged into an M<NUM>M<NUM> × N<NUM>N<NUM>/(a<NUM>a<NUM>) matrix, whose column and row indices correspond to location (position) and channel number in the transformed "reciprocal space", respectively; this is a specific example of a "matrix vector" representation. Mathematically, we express this operation via the (invertible) mapping <MAT> to denote the mapping between a P-dimensional vector z and a matrix τ(z) having dimension M<NUM>M<NUM> × N<NUM>N<NUM>/(a<NUM>a<NUM>); it will be appreciated that this rearrangement does neither affect the total number nor the values of the coefficients. Note further that for any vectorized image signal <MAT> the complex-valued matrix τ(AGx) satisfies certain conjugate-symmetry relations due to the real-valuedness of x which could be further exploited by means of a restriction operation expressed via the additional mapping <MAT> (see also the audio inpainting example described above). Its reversed operation, i.e. the conjugate-symmetric extension, is denoted by <MAT>. For notational simplicity, we use the abbreviation η for the composite mappings, i.e. <MAT> and <MAT>. We emphasize that the use of the restriction operation is optional; if no restriction operation is used, we set <MAT>, and <MAT>.

This ℓ<NUM>,∞-norm uses the matrix vector representation (where each sub-vector xi represents a cohort) and computes the maximum number among the number of non-zero entries in each of the sub-vectors. We propose to solve the following two optimization tasks, <MAT> <MAT> where, again, (7a) and (8a) present the formulation referred to as the analysis and the synthesis variant, respectively. In both cases the y that (jointly with B) achieves the minimum will be the reconstructed vectorized image signal; the reconstructed 2D image signal is then given by <MAT>. In the context of the invention, the ADMM technique is adapted to facilitate solving the problems of (7a) and (8a). Due to the introduced vectorized representations, the actual implementation is essentially the same as in the audio inpainting example described above; potential remaining minor adaptations are evident and straightforward.

The starting point is the sparse image inpainting formulation as discussed above. We aim at finding a frame satisfying the painless condition (<NUM>) with analysis operator A such that for the degraded vectorized image signal vecN<NUM>×N<NUM>{x} , A vecN<NUM>×N<NUM>{x} is as sparse as possible with respect to the ℓ<NUM>,∞-norm in a neighborhood of the gap. More formally, we intend to solve <MAT> where the notation | <IMG> represents the restriction to the neighborhood. Let AG denote the analysis operator of a given Gabor frame satisfying the painless condition (<NUM>). (The analysis vector η(AG vecN<NUM>×N<NUM>{x}) corresponds to the vector v that was mentioned in the introductory part as denoting the analysis coefficient vector. ) We next "deform" the Gabor frame using a unitary "deformation" operator <MAT> by defining a modified analysis operator <MAT>.

Note that we restrict to unitary deformations, since these preserve the painless condition (<NUM>) due to <MAT>. Furthermore, we are only interested in deformation operators which increase sparsity in frequency direction and which, in case the restriction operation σ(·) is applied, preserve the mentioned conjugate-symmetry relations of the DGT (for real-valued signals). Therefore, we additionally impose the following structure on W, <MAT> where <MAT> is a unitary matrix (i.e. dimension M<NUM>M<NUM> × M<NUM>M<NUM>). It is clear that (<NUM>) represents a well-defined unitary operator. It will be appreciated that (17a) contains a simple matrix product of the matrices V and τ(z) (thus largely reducing the complexity of operation as compared to a general P × P-matrix W acting on a P-dimensional vector z). This operation corresponds to a transformation which only combines ("mixes") the coefficients of the individual columns while excluding a mixing between columns. Moreover all columns are treated in the same manner, i.e. the mixing within each column is the same for all columns. Thus, (17a) denotes a special realization of a matrix block structure where the matrix represents a transformation which preserves the cohort structure.

It is remarked that, as mentioned above, in this embodiment the column index corresponds to the location (position) index; the index within a cohort (i.e. the row index) relates to 2D modulations (channel number in the transformed "reciprocal space"). In other implementations, depending on the respective implementation, the block structure may be finer (e.g. treating segments within a column as cohorts) or coarser (e.g. combining defined groups of columns into cohorts) or more complex (e.g. combining segments of columns, etc.).

For ease of understanding, we continue with the case that does not apply the restriction operation σ(·), so that M' = M<NUM>M<NUM>, <MAT>, and <MAT>. Combining (16a), (17a), and (15a), we obtain <MAT> where <IMG> denotes the set of all unitary M' × M' matrices V. However, the inventors noted that this is a highly non-convex problem, so that we can only hope to find an approximate solution. To this end, it is suggested that the -ℓ<NUM>-norms be relaxed to ℓ<NUM>-norms [<NUM>, <NUM>], and furthermore it may be suitable (although not necessary) to replace the max-operation in the definition of the ℓ<NUM>,∞-norm by a summation or other calculation, which is convenient in terms of calculation speed. An advantageous choice is replacing the norm in favor of an ℓ<NUM>-norm or a variant thereof, such as ℓ<NUM>,<NUM>-norm or ℓ<NUM>,∞-norm, where ℓ<NUM>,k-norm stands for the ℓk-norm evaluated for the M'-dimensional vector formed from the norms of the sub-vectors xi; thus, in the following we use an ℓ<NUM>,<NUM>-norm (which is equivalent to the ℓ<NUM>-norm on the entire vector). This results in the following optimization problem: <MAT>.

Suppose we have found an optimized unitary matrix V̂. From this it is also possible to construct the optimized matrix W of (17a); but in practice this is cumbersome because of the large dimension of W (as we had with the one-dimensional case above for audio inpainting). Therefore, we use a different approach: via (17a) and (16a), we obtain the analysis operator of the deformed frame as, <MAT>.

The above enables to re-formulate (7a) and (8a) in a way such that the learned sparsity-optimized frame is used instead of the Gabor frame. Note that optimization problem (23a) can be solved by the same basis optimization technique as used in the audio inpainting example described above.

Finally, this section deals with how to select the neighborhood <IMG> of the gap, which specifies the (input) training matrix (η(AGvecN<NUM>×N<NUM>{x})) |<IMG> used in (22a). It seems to be obvious that it should be as close as possible to the gap since this increases the chance that the learned sparse representation is also a sparse representation for the signal part within the gap. On the other hand, we also expect small "guard" areas between the gap and training borders to be mandatory to add, in order to avoid that unreliable data from the gap influences our training data. Regarding the size of the neighborhood we have to somehow rely on intuition. A larger neighborhood yields a larger training matrix (η(AGvecN<NUM>×N<NUM>{x})) |<IMG> and, probably, more accurate training results but it also increases the likelihood that variations of the image signal occur within the neighborhood. Such variations could yield a dictionary adapted to image parts which are essentially independent of the image within the gap. To illustrate convenient choices of the neighborhood, we will consider circular or rectangular shaped gaps in the following. The extension to more general shapes is evident by deforming and/or combining suitable base shapes.

Let us first assume that the image x ∈ [<NUM>,<NUM>]N<NUM>×N<NUM> has a circular shaped gap around a center point (c<NUM>,c<NUM>) ∈ {<NUM>,. , N<NUM> - <NUM>} × {<NUM>,. , N<NUM> - <NUM>} with radius r, i.e., that image samples are missing at positions <MAT>.

Then, it is convenient to choose a neighborhood with the shape of an annulus, i.e., <MAT> where rguard ≥ <NUM> takes into account the guard area and R determines the size of the neighborhood. Note that (η(AGvecN<NUM>×N<NUM>{x})) |<IMG> consists precisely of those columns of η(AGvecN<NUM>×N<NUM>{x}) that correspond to the cohorts specified by the neighborhood <IMG>.

Next, let us assume that the image x ∈ [<NUM>,<NUM>]N<NUM>×N<NUM> has a rectangular shaped gap, i.e., that that image samples are missing at positions <MAT>.

Then, it is convenient to choose a neighborhood of the following ("frame-like") shape, <MAT> with <MAT> <MAT> <MAT> <MAT> where <MAT> and <MAT> take into account the guard area and L<NUM> and L<NUM> determine the size of the neighborhood. Again, note that (η(AGveN<NUM>×N<NUM>{x})) | consists precisely of those columns of η(AGvecN<NUM>×N<NUM>{x}) that correspond to the cohorts specified by the neighborhood <IMG>.

Claim 1:
Computer-implemented method for signal processing,
comprising the steps of
providing an input signal (x) dependent on location within a predetermined signal domain in the form of a plurality of signal values associated with respective locations contained in a containing region within the signal domain, wherein the input signal is a time-dependent acoustic signal or an image signal, and the signal domain is an interval of time or a region defined in a two-dimensional plane, respectively,
defining a first portion within the containing region, said first portion being referred to as gap, wherein signal values associated with locations in the gap are assumed to be invalid,
providing a first set of dictionary signals (g) defined on the containing region, said first set being referred to as initial frame and its elements as initial dictionary signals, wherein the initial frame is usable for modeling the signal within the containing region by linear combination of the initial dictionary signals, the computer-implemented method being characterised by further comprising the steps of:
determining a second set of dictionary signals (ĝ) defined on the containing region, said second set being referred to as deformed frame and its elements as modified dictionary signals, wherein the modified dictionary signals are invertible linear combinations of the initial dictionary signals, namely, linear combinations according to an invertible matrix (Û, V̂, Ŵ) referred to as deformation matrix, wherein at least some of the modified dictionary signals are different from all initial dictionary signals,
wherein the deformation matrix (Û, V̂, W) is determined based on the signal values in the containing region outside the gap, by
calculating analysis coefficients for the input signal with respect to the initial frame, to obtain an analysis coefficient vector (v),
determining an invertible matrix (U) which optimizes a vector norm of a transformed coefficient vector (Uv) which is the result of transforming the analysis coefficient vector using the matrix, and
using the invertible matrix (U) thus obtained as the deformation matrix (Û, V̂, W), and
reconstructing the signal within the gap from the signal values in the containing region outside the gap by means of a sparse reconstruction method using the deformed frame (ĝ) or the deformation matrix (Û, V̂).