Jean-Loïc Le Carrou

François Gautier
email: francois.gautier@univ-lemans.fr

Roland Badeau

Sympathetic modes in the concert harp

Keywords: numbers: 43, 75, Gh

The concert harp is composed of a soundboard, a cavity with sound holes and 47 strings. When one string is plucked, other strings are excited and induce a characteristic 'halo of sound'. This phenomenon, called sympathetic vibrations is due to a coupling between strings via the instrument's body. These sympathetic modes generate the presence of multiple spectral components in each partial of the tone. Resolution of Fourier analysis does not permit their identification. A high resolution Method, called ESPRIT, is used to separate the spectral components which are very close one to another. Some of the measured spectral components in the analysed partials correspond to the response of sympathetic modes. The eigenfrequencies and mode shapes of these modes are investigated using a suitable model of the instrument : this model is based on a waveguide approach in which bending and longitudinal motions of 35 strings connected to an equivalent beam representing the soundboard are described. Identified experimental sympathetic modes are very well captured by the model.

I. INTRODUCTION

The harp is probably one of the oldest string instruments whose origin goes back to the Prehistory where the first men were charmed by the sound produced by their bow's string. This chordophone was first composed of a few strings attached on an arched frame on one side and on a soundboard on the other side. Then, with the increase of the number of strings, a pillar was added between the neck and the soundboard to support the strings' tension. This kind of harp was particularly used in Europe and marks the origin of the current concert harp. Nowadays, the concert harp is composed of 47 strings, from Cb0 (of fundamental frequency 30.9 Hz) to Gb7 (of fundamental frequency 2960 Hz), attached to the soundboard though an eyelet below which they are knotted. The soundboard is designed to withstand the stress imposed by the strings, as for the Camac concert harp used in this study, which is composed of multiple layers of different materials (aluminium, carbon, woods). From an acoustical point of view, the role of the soundboard is to radiate the sound produced by the vibrations of the strings. To some extent, this sound can also be amplified by the soundbox and its five sound holes [START_REF] Le Carrou | Experimental study of A0 and T1 modes of the concert harp[END_REF]. In spite of mechanical and constructional improvements, musicians and harp makers alike are annoyed by the feeling produced by the halo of sound when the instrument is played. Indeed, when one string is plucked, in some tuning configurations some others are also excited by sympathetic vibrations. Although this phenomenon is a fundamental characteristic of the instrument's sound, the instrument maker has to design the harp so that sympathetic vibrations remain reasonable.

The sympathetic couplings between strings lead to phenomena of aftersound and sometimes to beats. These phenomena are due to the fact that some partials of the sound may contain several spectral components whose frequencies are very close one to another. Two kinds of couplings are involved in this situation: couplings between different polarizations of a same string due to the way it is fixed and sympathetic couplings between different strings via the instrument's body [START_REF] Weinreich | Coupled piano strings[END_REF].

Considering the connection point between a single string and the soundboard as a point, the string/structure interaction can be described by a 6 by 6 admittance matrix. In such a description, three translational degrees of freedom and three rotational degrees of freedom are involved [START_REF] Boutillon | Measurement of the three-dimensional admittance at a violin bridge[END_REF]. As a consequence, for one single string, each mode can have up to 6 components, their frequencies being very close one to another in the string's response. However, in practice, 2 components often overshadow the others and correspond to 2 polarizations of the string. This has been shown for the guitar [START_REF] David | Caractérisations acoustiques de structures vibrantes par mise en atmosphère raréfiée. Méthodes d'estimation relatives aux fréquences et amortissements des modes pro-pres. Applications en acoustique musicale[END_REF]. The strong anisotropy of the string/structure interaction (for the out-of-plane and in-plane directions) is responsible for this effect. As a consequence, the decay rates of the two polarizations of a single string are generally different and lead to a double decay effect [START_REF] Weinreich | Coupled piano strings[END_REF].

The case of piano strings, which are grouped in pairs or triplet has often served as a framework for the study of sympathetic couplings. These couplings occur because each group of strings (pair or triplet) is practically tuned in unison. The case of two strings tuned in unison having the same polarization, coupled by a bridge motion is studied in [START_REF] Weinreich | Coupled piano strings[END_REF]: normal modes for this configuration appear in pairs and depend on the strings mistuning and bridge admittance. It is the presence of these two coupled modes that is responsible for phenomena of beats and aftersound. The more complex case of a pair of piano strings tuned in unison, each of them having two polarizations has been studied in [START_REF] Capleton | False beats in coupled piano string unisons[END_REF].

Similar coupling mechanisms exist in string instruments which do not have pairs of strings tuned in unison. This has been illustrated for the American five-string banjo: in [START_REF] Dickey | The structural dynamics of the American five-string banjo[END_REF], this instrument is modelized by an assembly of one-dimensional sub-systems in which waves propagation occur, allowing the time-domain response to be computed. It is shown that when all the strings are incorporated to the model, the decay time is shorter than when only one string is considered. This is explained by the presence of sympathetically driven strings, even if it is difficult to find out how sympathetic vibrations occur with this model.

For the kantele, which is a Finnish plucked five-string instrument, sympathetic vibrations are also identified in [START_REF] Erkut | Acoustical analysis and model-based sound synthesis of the kantele[END_REF]. Through the experimental analysis of the total amount of energy transferred from one string to all the others, it can be shown that the transfer of energy is more pronounced between strings which have simple harmonic relationships.

In a previous paper [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF], it has been shown that the sympathetic phenomenon is due to the presence of particular modes, called sympathetic modes, in the modal basis of the system. These modes have been both theoretically and experimentally identified on a simple academic configuration close to the harp: two strings tuned to the octave and connected to a beam were considered and the modes of this assembly were investigated. The aim of the present paper is to identify the sympathetic modes in the response of a real concert harp.

In the first part, an adequate experimental setup and an analysis based on a 'High Resolution' method are carried out to identify frequency components present in one partial when one string is plucked. Note that 'High Resolution' techniques are suitable for the spectral analysis of signals having very close spectral components, which cannot be separated using a Fourier analysis because of a lack of resolution. In the second part, the sympathetic vibrations are investigated through the use of a physical model of a simplified concert harp. In the final part, a comparison between theoretical and experimental results give an explanation of the origin of the sympathetic phenomenon in the instrument.

II. EXPERIMENTAL INVESTIGATION OF THE HARP'S SYMPATHETIC MODES A. Experimental setup

In order to experimentally highlight the presence of sympathetic modes in the concert harp, the following experimental protocol has been set up: the harp is plucked whereas an accelerometer is used to measure the vibratory signal on the soundboard at point A, as shown in Figure 1. This point is located on the inner surface of the soundboard, between the Db 3 -string and the Cb 3 -string, respectively labelled 30 and 31. The played string is string 31 and the other strings are either damped or stopped during oscillations. Eight experimental configurations presented in Table I have been defined to investigate the characteristics of 2. These signals highlight a fact well-known by the harpists: the stopping of a string plucked does not necessarily imply a fast decrease of the sound level. Indeed, as shown in configurations (3), ( 5), ( 6), [START_REF] Erkut | Acoustical analysis and model-based sound synthesis of the kantele[END_REF] and ( 8), the end signal amplitude is of the same order of magnitude that in the free-strings B. Method used for the extraction of modal parameters

Introduction: limit of the Fourier analysis

The Fourier spectrum of the accelerometer's signal at point A is shown in Figure 3.

When string 31 is plucked, several coupled modes of the system respond to produce the harp sound. This spectrum is composed of partials which are quasi harmonic. Using an appropriate selective filter, a zoom on this spectrum in the vicinity of the first partial and its corresponding time waveform are presented in Figure 4. This figure clearly shows that several sinusoidal components are present in the vibratory signal. Indeed, since free oscillations of the system occur, the response of the instrument to the plucking action measured by the accelerometer, corresponds to the superposition of several modes whose frequencies are very close one to another. In the zoomed area of the spectrum, peaks of every component cannot be clearly separated whereas in the time waveform several sinusoidal components with different damping factors and frequencies appear. This shows the limit of the Fourier analysis. The presence of several vibration modes having close frequencies is a common characteristic for numerous free stringed instruments. Various methods, applied to musical instruments, such as techniques based on the Hilbert transform [START_REF] David | Caractérisations acoustiques de structures vibrantes par mise en atmosphère raréfiée. Méthodes d'estimation relatives aux fréquences et amortissements des modes pro-pres. Applications en acoustique musicale[END_REF][START_REF] Rossi | Instantaneous frequency and short term Fourier transforms: application to piano sounds[END_REF] or High Resolution methods [START_REF] David | Caractérisations acoustiques de structures vibrantes par mise en atmosphère raréfiée. Méthodes d'estimation relatives aux fréquences et amortissements des modes pro-pres. Applications en acoustique musicale[END_REF][START_REF] Laroche | The use of the matrix pencil method for the spectrum analysis of musical signals[END_REF][START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoides modulées. Application aux signaux de musique[END_REF], allow the identification of spectral components.

For the identification of modes with close frequencies in the first partial, we choose to use a High Resolution method: the ESPRIT algorithm (Estimation of Signal Parameters via Rotational Invariance Techniques) [START_REF] Roy | ESPRIT -A subspace rotation approach to estimation of parameters of cisoids in noise[END_REF]. A brief description of this technique and the specificity of its implementation in our context are given in paragraphs II B 2 to II B 4. An application to the harp's signals is then presented in section II C. set of disturbing components, which can be viewed in terms of signal and noise subspaces.

In this framework, the ESPRIT algorithm is based on a particular property of the signal subspace, referred to as the rotational invariance. This property permits to extract the model parameters from the eigenvalues of a so-called spectral matrix, which is obtained from the estimated signal subspace.

The noiseless Exponential Sinusoidal Model (ESM) defines the discrete signal x(t) as a sum of complex exponentials:

x(t) = K k=1 a k e δ k t e j(2πf k +ϕ k ) , t ∈ [0, N -1] (1)
where each frequency

f k ∈ [-1 2 , 1 2 ] is associated to a real magnitude (a k > 0), a phase (ϕ k ∈ [-π, π]
) and a damping or amplification factor (δ k ∈ R). The whole number K is the number of complex exponentials, also called model order, and N is the number of the signal's samples. By defining the complex amplitudes α k = a k e jϕ k and the complex poles

z k = e δ k +j2πf
k , which are supposed to be distinct, the signal model x(t) can be re-written in the following form:

x(t) = K k=1 α k z t k . (2)
For any time t, the data vector

x(t) = [x(t), . . . , x(t + n -1)] T of dimension n > K belongs
to the K-dimensional signal subspace spanned by the Vandermonde matrix

V n =                1 1 • • • 1 z 1 z 2 • • • z K . . . . . . . . . z n-1 1 z n-1 2 • • • z n-1 K                .
It can be noted that this Vandermonde matrix satisfies the following rotational invariance property:

V ↑ = V ↓ D, where D = diag(z 1 . . . z K ), V ↓ is the matrix extracted from V by
deleting the last row, and V ↑ is the matrix extracted from V by deleting the first row.

In practice, the measured signal s(t) is corrupted by an additive noise:

s(t) = x(t) + w(t),
where w(t) is assumed to be white. Although matrix V is unknown, the signal subspace can still be estimated as the principal eigensubspace of the correlation matrix Rss of the measured signal, defined as follows:

Rss = 1 N -n + 1 SS H , (3)
where

S =                s(0) s(1) • • • s(N -n) s(1) s(2) • • • s(N -n -1) . . . . . . . . . . . . s(n -1) s(n) • • • s(N -1)                (4)
is the n × (N -n + 1) Hankel data matrix which involves N successive samples of the signal and the exponent H is the hermitian conjugate. Thus the n × K matrix W formed by the Since the matrices W and V span the same subspace, there exists a K × K non-singular matrix G such that V = W G. It can then be noted that W satisfies an invariance property similar to that of the Vandermonde matrix:

W ↑ = W ↓ Φ, where the K × K matrix Φ = G D G -1
, referred to as the spectral matrix, is similar to the diagonal matrix D. In particular, the eigenvalues of Φ are the complex poles z k .

Finally, the ESPRIT algorithm consists of the following steps:

1. compute the signal subspace basis W by means of an eigenvalue decomposition, 2. compute the spectral matrix Φ by means of the least squares method:

Φ = W + ↓ W ↑ , (5)
where the symbol + denotes the Moore-Penrose pseudo-inverse, 3. estimate the complex poles z k as the eigenvalues Φ.

It is proved that the best performance in terms of statistical efficiency is obtained for a proper dimensioning of the data matrix S: n = N/3 or n = 2N/3 [START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoides modulées. Application aux signaux de musique[END_REF].

In a second stage, the complex amplitudes α k , grouped in a K × 1 vector denoted α, are obtained thanks to the least squares method:

α = [V N ] + s. (6)
In equation ( 6), s denotes the vector containing N successive samples of the signal and

V N is the N × K Vandermonde matrix defined by the poles estimated by the ESPRIT method. Parameters a k and ϕ k are directly deduced as the modulus and phase of the complex amplitudes α k .

Estimation of the number of components

The main difficulty of the method consists in evaluating the number K of components present in the signal. The technique commonly used is the over-estimation of this number and the discrimination of spurious results by means of an indicator such as the components energy or the error between the measured signal and the model. Other more effective methods exist for estimating the model order K such as the ESTimation Error (ESTER) method [START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoides modulées. Application aux signaux de musique[END_REF][START_REF] Badeau | A new perturbation analysis for signal enumeration in rotational invariance techniques[END_REF] also used in the study. It consists in the computation of an inverse error function,

J : p → 1 E(p) 2 2 (7)
where

E(p) = W ↑ (p) -W ↓ (p)Φ(p), (8)
for all possible orders 0 < p < n -1. For determining the value of K, we choose the greatest value of p for which the function J(p) reaches a maximum which is greater than a threshold chosen arbitrarily above the noise level contained in the signal.

ESPRIT method implementation

In order to minimize the computational time and increase the results accuracy [START_REF] Badeau | Méthodes à haute résolution pour l'estimation et le suivi de sinusoides modulées. Application aux signaux de musique[END_REF], the ESPRIT method implementation is carried out according to the following procedure: after centering the studied partial around the null frequency, a finite impulse response (FIR) filter selects the frequency range containing the partial to analyse. The filter is chosen with a linear phase to keep the signal waveform. The filter is known to have a finite transitory response which corresponds to the length of its impulse response. The first filtered signal points belonging to this transitory phase are thus removed from the processing afterwards [START_REF] David | Caractérisations acoustiques de structures vibrantes par mise en atmosphère raréfiée. Méthodes d'estimation relatives aux fréquences et amortissements des modes pro-pres. Applications en acoustique musicale[END_REF].

The filtered and centered signal is highly decimated to limit the computational time of the ESPRIT method. After the estimation of the model order by the ESTER method, as previously explained in section II B 3, the ESPRIT algorithm is then applied. The final validation of the model order is performed using a comparison between the measured and synthesized signals. An illustration of the method implementation is shown in Figure 5.

Note that in the Exponential Sinusoidal Model, the signal is assumed to be complex. For a musical sound, the signal is real and can be written as follows:

x(t) = K k=1 A k e δ k t cos(2πf k t + ϕ k ) (9)
which can be re-written with exponential terms:

x(t) = K k=1 A k 2 e δ k t e j(2πf k t+ϕ k ) + K k=1 A k 2 e δ k t e -j(2πf k t+ϕ k ) (10)
Maximum number of components p Partial to analyse s(t)

Parameters : f k , δ k , a k , φ k
Centering of the partial around the null frequency FIR Filtering around the partial and decimation

Application of the ESPRIT method

Estimation of the number of components K (ESTER method) Because of the filtering in the spectrum around the studied partial, the studied signal corresponds to the first part (with positive frequency) of the equation [START_REF] Laroche | The use of the matrix pencil method for the spectrum analysis of musical signals[END_REF]. In order to find the real amplitude of the measured signal, the amplitude A k obtained by the ESPRIT method has to be multiplied by two.

C. Results

The implemented ESPRIT method, as previously explained, is applied to the signals measured on the concert harp in its final part, between 4s and 8s. The different experimental conditions are described in Table I and the attention is focused on the first partial of string 31. The estimated components found for this partial are gathered in Table II. The model order stretches from 0, for configuration (4), to 4, for configuration (1). Note that for each configuration, coefficients are computed from 5 measurements for estimating a repeatability uncertainty. In Table II, the components are classified in such a way that their amplitude are in descending order from top to bottom. Moreover, they are aligned in order to facilitate the reading of the missing components.

From a general point of view, the repeatability uncertainties calculated for different parameters are extremely small, lower than 1% for the frequencies and around 10% for the damping factors. Nevertheless, for some components the uncertainty can be important as for the last component of the ( 1) configuration. This can happen for components not very present, with a weak energy, and a weak signal-to-noise ratio.

When the instrument is not modified (( 1) and ( 3) to ( 8) configurations) the components' parameters slightly vary, less than 0,5% for the frequencies and, in a maximum of 30% of variation for the damping factors. This result shows that the stopped string does not modify the vibratory behavior of the instrument. Nevertheless, when paper is added to damp all strings except string 31, in configuration (2), the components' damping factors are modified and a slight shift in frequency is caused. The instrument is slightly modified by the added paper but not enough to allow the identification of the components' nature.

Thanks to the developed experimental protocol and the identification method, the components present in the first partial are obtained. The results show that the estimated parameters are stable following the eight experiences. These results are compared in the following section to the modal basis of a simplified harp.

Vibrations of the studied concert harp have been investigated in a previous paper [START_REF] Le Carrou | Experimental study of A0 and T1 modes of the concert harp[END_REF] thanks to an experimental modal analysis. In low frequencies, six consecutive modes have been identified from 24 Hz to 181 Hz. The investigation of the mode shapes show two particular characteristics: mode shapes are symmetric according to the strings plane and for all modes, the bending motion of the soundboard is similar to the first mode shape of a clamped-clamped beam. Thus, although the soundboard is a complex assembly, constituted of a sandwich of several layers of different wood glued together reinforced by a central aluminum bar and by two lateral stiffeners in wood, the soundboard can be described by using an equivalent beam clamped at both ends. The typical first bending mode at 152.2 Hz is shown in Figure 6-A and its corresponding mode shape can schematically be described as an important deflection in the lower two thirds of the soundboard between the two clamped points, one at the pillar level and the other one at the 11-string level. Indeed, on the treble strings, the soundboard is rigidified by the proximity of the edge linked to the soundbox, explaining the absence of movement in low frequencies. Thus, the beam length clamped at both ends is limited to the distance included between the harp's pillar and the 11-string.

According to this experimental result, the vibratory model is composed of a 1m-long beam on which 35 strings are attached as shown in Figure 6-B.

The mechanical properties of the equivalent beam and of each string have to be determined: for the strings, most parameters are directly measured on the harp [START_REF] Carrou | Vibro-acoustique de la harpe de concert (vibro-acoustics of the concert harp)[END_REF] except for the Young's modulus and the density which are supposed to equal the data given in [START_REF] Bell | The physical properties of gut musical instrument strings[END_REF] and [START_REF] Bell | An acoustical investigation of the Concert Harp[END_REF]. Values of the tension are computed from the fundamental frequency of the tones, by considering strings fixed at both ends. For the equivalent beam, the determination of its parameters is more complicated since the geometrical parameters are directly measured on an isolated soundboard, allowing the evaluation of a mean density (ρ = 553 Kg/m 3 ). The area of the cross section A and the second moment of area J are measured at different positions along the axis. Average values of these two parameters are thus retained to characterize the equivalent beam (A eq = 38.3 cm 2 and J eq = 38.9 cm 4 ). Finally, the Young's modulus E of the equivalent beam is determined so that the first bending mode of the beam in the beam-35 strings assembly equals 150 Hz, the eigenfrequency of the fourth mode shown in Figure 6-A (E = 5.9 GPa).

The vibratory model of the harp can thus be considered as an assembly of an equivalent beam connected to 35 strings as described in Figure 6-B. This vibratory model is afterwards named the simplified harp.

B. Harp's sympathetic modes

The modal basis of the simplified harp is computed using the transfer matrix method, which is appropriated for modeling assemblies of one-dimensional sub-structures. This computation is based on a wave guide model for each sub-structure of the beam-strings assembly.

A brief summary of the different steps of the method is proposed in appendix and details are developed in [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF]. Eigenfrequencies are obtained for singularities of a characteristic matrix labeled R R (see equation ( 15) page 31). Each drop of the logarithm of matrix R R 's determinant, defined in equation ( 16), corresponds to an eigenfrequency. [110-160 Hz], showing that modes have very close eigenfrequencies. Among the 151 modes found in the range [0-500 Hz], examples of four modal shapes are presented in Figure 8.

These shapes clearly show that for each mode each substructure does not interact in the same manner. To classify modes of the assembly, we use a criterion called Kinetic Energy Ratio (KER) defined by

KER j (k) = l k 0 ρ k Ψ T j (x)Ψ j (x) dx r lr 0 ρ r Ψ T j (x)Ψ j (x) dx . (11)
In this equation ρ is the mass per unit length of the sub-structure, Ψ j is the mode shape of mode j and x is the generic space variable defined in the Appendix. The KER corresponds to the ratio of kinetic energy of one sub-structure k (of length l k ) divided by the total kinetic energy of the structure. Its value is a percentage and allows us to identify the relative importance of each sub-structure displacement field. In the study, this percentage is rounded to the nearest whole number. Thus, for a given mode, a null value of the KER of a sub-structure indicates that this sub-structure is inactive. Values of KER on each sub-structure are used to classify modes into four groups [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF]: beam modes, string modes, string-string modes, beam-string modes.

In Table III 3) and ( 4)), these two components disappear. In playing configuration, these two sinusoids have close frequencies, separated by only 0.2 Hz, and have very different damping factors and amplitudes. This result thus points out that the string vibrates following its two polarizations. These two polarizations are excited by the harp player in a different way depending on the string plucking, involving different initial amplitudes (showed by the difference of energy of each component). Moreover, as for the piano, one of the polarizations seems to be favored for transmitting its energy to the soundboard, implying a rapidly de-caying component at the plucking moment [START_REF] Weinreich | Coupled piano strings[END_REF]. This result shows that for a good modeling of the strings' modal behavior, the two polarizations of string vibrations have to be taken into account. In the simplified harp, only one polarization is considered, which constitutes a limitation.

When the string is stopped during the oscillations, modes present in the instrument response are selected. In configurations (3), ( 5) and ( 6 Although the model of the instrument suits the identification of sympathetic modes, it can be extended in order to take into account two polarizations per strings. Moreover, from an experimental point of view, the analysis of partials of higher order, where other phenomena such as octave vibrations can be present, could also be performed.

Appendix: determination of the beam-strings assembly modes

In this appendix, we summarize the method for computing normal modes of the beamstrings assembly, illustrated in Figure 9. The detailed method is published in [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF]. In harmonic regime, the vibratory state of each sub-structure (string and beam) of the assembly is described at any point of its neutral line by a state vector X (x) whose compo-nents are kinematic and force variables. Variable x is the spatial coordinate describing the current point on the sub-structure. For the string, the state vector is the following vector:

X s (x) = u s (x) w s (x) N s (x) Q s (x) T , (12)
where u s , w s are longitudinal and transversal displacements and N s , Q s are the longitudinal and transversal forces. In the beam, the state vector is written as follows:

X b (x) = u b (x) w b (x) θ b (x) N b (x) Q b (x) M b (x) T , (13)
where u b , w b and θ b are respectively the longitudinal and transversal displacements and the slope of the beam cross section, and where N b , Q b and M b are respectively the longitudinal and transversal forces and the bending moment. With the transfer matrix method, it can be shown that the vibratory state of each sub-structure at any point x can be expressed in function of the vibratory state at another point x 0 with the following expression:

X (x) = T(x, x 0 )X (x 0 ), (14)
where T(x, x 0 ) is the transfer matrix of the sub-structure (beam or string) between the state vector at point x 0 and the state vector at point x. By using coupling equations at beam-string connection [START_REF] Le Carrou | Modelling of sympathetic string vibrations[END_REF], it can be shown that the force components of state vectors on boundary X f verify the following equation: R R (ω)X f = 0, [START_REF] Bell | The physical properties of gut musical instrument strings[END_REF]

FIG. 1 :

 1 FIG. 1: Schematical representation of the experimental setup.

FIG. 3 : 1 FIG. 4 :

 314 FIG.3: Spectrum of the signal measured at point A when string 31 is plucked.

FIG. 5 :

 5 FIG.5: Summary of the implementation of the ESPRIT method on the studied partial.

11 FIG. 6 :

 116 FIG. 6: (A) Modal shape associated to the fourth mode and description of the vibratory profile in

FIG. 7 :FIG. 8 :

 78 FIG.7: Logarithm of the determinant of the characteristic matrix R R as function of frequency.

  are present in the instrument's response. These sympathetic modes are due to the coupling between different strings via the soundboard. With a waveguide model of the harp, their eigenfrequencies and mode shapes can be accurately determined. This model describes bending and longitudinal motions in the strings, connected to an equivalent beam representing the soundboard and allows the modal basis of the strings-beam assembly to be computed. (b) In the time domain, sympathetic vibrations generate multiple components in the string's partials. Resolution of the Fourier analysis does not permit their identification.This identification is performed using the ESPRIT method, which is a High Resolution method. The determination of the number of elementary spectral components, which is the main difficulty in the implementation of High Resolution methods is successfully performed using the ESTER method.(c) Several components are identified in the first partial of string 31: two of them correspond to two strings modes having different polarizations. Other components correspond to sympathetic modes. In the analysed example, two sympathetic modes are involved. Their number and their eigenfrequencies are very well captured by the proposed model.

FIG. 9 :

 9 FIG. 9: Diagram of a clamped-clamped beam connected to several fixed strings. It includes local

TABLE I :

 I Characteristics of the experimental configurations. String 24 corresponds to Cb 2 of fundamental frequency 246.9 Hz, string 31 to Cb 3 (123.5 Hz), string 35 to Fb 3 (82.4 Hz), string 38 to Cb 2 (61.7 Hz) and string 42 to Fb 1 (41.2 Hz).

	Configuration	Description
	(1)	All strings free to vibrate
	(2)	All strings damped except string 31
	(3)	String 31 stopped during oscillations
	(4)	Strings 24, 31, 35, 38 and 42 stopped during oscillations
	(5)	Strings 31 and 24 stopped during oscillations
	(6)	Strings 31 and 35 stopped during oscillations
	(7)	Strings 31 and 38 stopped during oscillations
	(8)	String 31 and 42 stopped during oscillations
	sympathetic vibrations.
	All measured signals are sampled to 4096 Hz and last 8 seconds. The string is stopped
	by the harp player a few seconds after the plucking of string 31. The temporal signals
	measured following the eight configurations are shown in Figure

TABLE II :

 II Frequencies and damping factors of identified components in the first partial of the accelerometer signals in the eight experimental configurations defined in TableI. The reported uncertainty is an uncertainty with a 95% confidence interval.

					Configurations
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
	123,49±0,01	123,54±0,01
		-0,68±0,02	-0,57±0,01
	123,34±0,01		123,35±0,00		123,34±0,00 123,34±0,00 123,30±0,01
	Components	-0,15±0,03		-0,11±0,01		-0,11±0,01	-0,11±0,01	-0,10±0,01
	123,08±0,00		123,09±0,01		123,08±0,00 123,08±0,00		123,09±0,01
		-0,26±0,01		-0,21±0,01		-0,22±0,03	-0,21±0,02		-0,22±0,07
	123,78±0,02	123,71±0,03
		-0,15±0,12	-0,54±0,06

  , the eigenfrequencies and KER of each sub-structure of modes 26, 27, 28 and 33 are reported. For these four modes, the KER is significant only for four sub-structures: the beam and strings 31, 38 and 42. String 31 (Cb 3 note of fundamental frequency 123.5 Hz) corresponds to the upper octave of string 38 (Cb 2 note at fundamental frequency 61.7 Hz)

TABLE III :

 III Eigenfrequencies and Kinetic Energy Ratio for each sub-structure (strings 42, 38, 35, 31, 24 and the beam). KER expressed in % and rounded to the nearest whole number.

		Frequency			KER (%)
	Mode	(Hz)	string 42 string 38 string 35 string 31 string 24 Beam
	26	122.95	0	99	0	1	0	0
	27	123.29	97	1	0	2	0	0
	28	123.48	1	2	0	97	0	0
	33	149.91	0	0	1	0	0	73
	which is the upper fifth of string 42 (Fb 1 at fundamental frequency 41.2 Hz). The beam's
	KER nearly equals zero for modes 26 to 28, showing that mode shapes are dominated by
	the string's motion since the KER is distributed according to two or three strings, allowing
	the definition of these modes as string-string modes or sympathetic modes [8]. Actually,
	according to the modal superposition principle, if string 31 is plucked, modes 26 to 28 are
	set into vibration and so are strings 38 and 42, generating the phenomenon of sympathetic
	vibrations. Note that for these three sympathetic modes, modifications of the character-
	istics of the equivalent beam only slightly modify the KER distribution, not altering our
	conclusions. Nevertheless, it appears that for other vibratory modes the KER distribution
	can be affected.

TABLE IV :

 IV Comparison between experimental and theoretical results obtained from the vibratory model of the concert harp.

	Modal frequency
	Experimental Theoretical Error
	123.49 Hz	123.48 Hz 0.01 %
	String 31 mode
	123.78 Hz
	Sympathetic mode 31-42 123.30 Hz	123.29 Hz 0.01 %
	Sympathetic mode 31-38 123.09 Hz	122.95 Hz 0.11 %
	IV. DISCUSSION
	When string 31 of the concert harp is free to oscillate (configuration (1) in Table I), two
	sinusoidal components are identified in the first partial. When this string is stopped (config-
	urations (2), (

  ), the same components are present, showing that modes implying strings 24 and 35 do not participate in the first partial of the instrument response. By stopping strings 31 and 38, one vibrating component disappears and by stopping strings 31 and 42, another one is absent from the instrument's response.

	With configurations (7) and (8), we can deduce that strings 38 and 42 participate in the
	sound radiated by the instrument.
	When the harpist plucks string 31, four modes are set into vibrations at the same time:
	two modes involving string 31 (one mode per polarization), one mode involving strings 31-38
	and one mode involving strings 31-42. These last two modes are thus sympathetic modes.
	The stopping of string 31 does not necessary lead to the weakening of these modes, proving
	that the kinetic energy present in strings 42 and 38 is definitely more important than in
	string 31. Moreover, these results show that the sympathetic modes coming from each
	polarization of string 31 are not visible in the response. This fact can be explained by the

weak energies of these modal components, not allowing them to emerge from the noise.

The ESPRIT methodSubspace-based High Resolution methods such as the ESPRIT algorithm are of major interest for estimating mixtures of complex exponentials, because they overcome the spectral resolution limit of the Fourier transform and provide very accurate estimates of the signal parameters. These methods consist in splitting the observations into a set of desired and a

Acknowledgement

The authors acknowledge the instrument maker CAMAC Harps for the lending of the concert harp.

These experimental results can be compared to those obtained from the theory. Experimentally, we found two modes which mostly involve string 31 at 123.59 Hz and at 123.78 Hz.

With the vibratory model, the string mode associated to string 31 is found at 123.48 Hz.

Note that the model takes only one string polarization into account. The agreement between model and measurement is, as expected, very good since the tension value of each string of the simplified harp has been fixed in such a way that the eigenfrequencies of uncoupled strings correspond to those of the real instrument.

V. CONCLUSION

The numerous strings of the concert harp induce sympathetic vibrations, which are responsible for the feeling of halo of sound. This characteristic is important and constitutes a signature of the instrument. In this paper, experimental and theoretical investigations have where R R is defined as a characteristic matrix of the assembly. In the case of free vibrations, vector X f is not equal to zero and the matrix R R is not invertible, implying:

The angular eigenfrequencies ω i of the system are obtained as the solutions of equation [START_REF] Bell | An acoustical investigation of the Concert Harp[END_REF]. For each ω i , relation [START_REF] Bell | The physical properties of gut musical instrument strings[END_REF] provides the components of X f from which we can calculate all state vector components at any point on the assembly.


NMF with time-frequency activations to model non stationary audio events

Romain Hennequin, Student Member, IEEE, Roland Badeau, Member, IEEE, and Bertrand David, Member, IEEE Abstract-Real world sounds often exhibit time-varying spectral shapes, as observed in the spectrogram of a harpsichord tone or that of a transition between two pronounced vowels. Whereas the standard Non-negative Matrix Factorization (NMF) assumes fixed spectral atoms, an extension is proposed where the temporal activations (coefficients of the decomposition on the spectral atom basis) become frequency dependent and follow a timevarying ARMA modeling. This extension can thus be interpreted with the help of a source/filter paradigm and is referred to as source/filter factorization. This factorization leads to an efficient single-atom decomposition for a single audio event with strong spectral variation (but with constant pitch). The new algorithm is tested on real audio data and shows promising results. Index Terms-music information retrieval, non-negative matrix factorization, unsupervised machine learning.

I. INTRODUCTION

T HE decomposition of audio signals in terms of elemen- tary atoms has been a large field of research for years. As we usually encounter very dissimilar audio events (both in their spectral and temporal characteristics), the decomposition on a single basis (such as the Fourier basis) is generally not sufficient to accurately explain the content of a large class of signals. Sparse decomposition techniques [START_REF] Mallat | Matching pursuit with timefrequency dictionnaries[END_REF] use a redundant dictionary of vectors (called atoms) and try to decompose a signal using few of them (much less than the dimension of the space): thus the signal can be accurately decomposed with few elements. When atoms are designed to fit the signal (for instance harmonic atoms for musical signals [START_REF] Gribonval | Harmonic decomposition of audio signals with matching pursuit[END_REF]), these elements become more meaningful, and then a supervised classification can be performed to cluster atoms corresponding to a real event in the signal [START_REF] Leveau | Instrument-specific harmonic atoms for mid-level music representation[END_REF]. These methods are quite powerful and give good results. However, since the dictionary is fixed, it must be designed to fit all possible signals, which is not achievable in practice. Recently methods of data factorization were proposed to simultaneously extract atoms from the signal and provide a decomposition on these atoms. These techniques that we call factorization make use of the natural redundancy of the signal, mimicking human cognition which utilizes this redundancy to understand visual and audio signals: Principal Component Analysis, Independent Component Analysis [START_REF] Comon | Independent component analysis, a new concept[END_REF], [START_REF] Makino | Audio source separation based on independent component analysis[END_REF], sparse coding [START_REF] Samer | Unsupervised analysis of polyphonic music by sparse coding[END_REF] or NMF [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] have been introduced both to reduce the dimensionality and R. Hennequin, R. Badeau and B. David are with the Institut Telecom, Telecom ParisTech, CNRS LTCI, 46, rue Barrault -75634 Paris Cedex 13 -France (email: <forename>.<surname>@telecom-paristech.fr)

The research leading to this paper was supported by the French GIP ANR under contract ANR-06-JCJC-0027-01, DESAM, and by the Quaero Program, funded by OSEO, French State agency for innovation.

Manuscript received December 8, 2009 to explain the whole data set by a few meaningful elementary objects. Thanks to the non-negativity constraint, NMF is able to provide a meaningful representation of the data: applied to musical spectrograms it will hopefully decompose them into elementary notes or impulses. The technique is widely used in audio signal processing, with a number of applications such as automatic music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF], [START_REF] Paulus | Drum transcription with non-negative spectrogram factorization[END_REF], [START_REF] Bertin | Blind signal decompositions for automatic transcription of polyphonic music: NMF and K-SVD on the benchmark[END_REF] and sound source separation [START_REF] Cichocki | New algorithms for nonnegative matrix factorization in applications to blind source separation[END_REF], [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF], [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. However, the standard NMF is shown to be efficient when the elementary components (notes) of the analyzed sound are nearly stationary, i.e. when the envelope of the spectra of these components does not change over time. Nevertheless, in several situations, elementary components can be strongly non stationary. In this article, we will focus on timbral variability, i.e. variability of the spectral shape that we can find in plucked strings sounds or singing voice (sounds of different vowels present greatly dissimilar spectral shapes). However, we will not address pitch variability that is encountered in vibrato or prosody. In case of a noticeable spectral variability, the standard NMF will likely need several nonmeaningful atoms to decompose a single event, which often leads to a necessary post-processing (to cluster the different parts of a single source [START_REF] Fitzgerald | Non-negative tensor factorisation for sound source separation[END_REF]). To overcome this drawback, Smaragdis [START_REF] Smaragdis | Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs[END_REF] proposes a shift-invariant extension of NMF in which time/frequency templates are factorized from the original data: each atom then corresponds to a time-frequency musical event able to include spectral variations over time. This method gives good results, but does not permit any variation between different occurrences of the same event (atom), its duration and spectral content evolution being fixed.

In this paper, an extension of NMF is proposed in which temporal activation becomes frequency dependent: it thus can be interpreted with the help of the classical source/filter paradigm as a source/filter factorization. Our method includes AutoRegressive Moving Average (ARMA) filters estimated from the data, associates a time-varying filter with each source and learns the sources (atoms) in a totally unsupervised way. this method presents some similarity with Durrieu's work [START_REF] Durrieu | An iterative approach to monaural musical mixture de-soloing[END_REF], [START_REF] Durrieu | Singer melody extraction in polyphonic signals using source separation methods[END_REF] in which a source/filter model is used in a NMF framework to extract the main melody of musical pieces. This model permits to efficiently take the strong spectral variations of the human voice into account.

However, our approach is quite different since, in opposition to Durrieu's work, sources are learnt, a time-varying filter is associated with each source and the class of filter we use is more standard.

In section II, we introduce the source/filter decomposition as an extension of NMF. In section III, we derive an iterative algorithm similar to those used for NMF to compute this decomposition. In section IV, we present experiments of source/filter decomposition of the spectrogram of three different sounds, and compare this decomposition to the standard NMF. Conclusions are drawn in section V.

II. MODEL

A. NMF and extension

Given an F × T non-negative matrix V and an integer R such that F R + RT ≪ F T , NMF approximates V by the product of an F × R non-negative matrix W and an R × T non-negative matrix H:

V ≈ V = WH i.e. V f t ≈ Vft = R r=1 w f r h rt ( 1
)
This approximation is generally quantified by a cost function C(V, W, H) to be minimized with respect to (wrt) W and H. A common class of cost functions is designed elementwise:

C(V, W, H) = F f =1 T t=1 d(V f t | Vft )
where d is a scalar divergence (i.e. a function such that ∀(a, b) ∈ R 2 d(a|b) ≥ 0 and d(a|b) = 0 if and only if a = b). Several classes of such divergences have been proposed, for instance Bregman divergences [START_REF] Dhillon | Generalized nonnegative matrix approximations with Bregman divergences[END_REF] and β-divergences [START_REF] Kompass | A generalized divergence measure for nonnegative matrix factorization[END_REF], [START_REF] Cichocki | Csiszars divergences for nonnegative matrix factorization: Family of new algorithms[END_REF]. In this paper, we will focus on the β-divergence which includes usual measures (Euclidean distance, Kullback-Lieber divergence and Itakura-Saito divergence). The β-divergence is defined for β ∈ R\{0, 1} by:

d β (x|y) = 1 β(β -1) (x β + (β -1)y β -βxy β-1 ) (2)
For β = 0 and β = 1, the β-divergence is defined by continuity:

d 0 (x|y) = x y -log x y -1 d 1 (x|y) = x(log x -log y) + (y -x)
For these values, the β-divergence respectively corresponds to Itakura-Saito divergence and Kullback-Leibler divergence, and for β = 2, it corresponds to the Euclidean distance.

One could notice that the singularities in the definition of the β-divergence for β = 0 and β = 1 no longer appear in the partial derivative wrt y: this partial derivative is useful for designing descent methods in order to minimize the cost function C. The first order partial derivative of the β-divergence wrt y is for all β ∈ R:

∂d β (x|y) ∂y = y β-2 (y -x)
When applied to power (squared magnitude) spectrograms, NMF factorizes data into a matrix (or basis) of frequency templates which are the R columns of W and a matrix H whose R rows are the temporal vectors of activations corresponding to each template. For a musical signal made of several notes played by the same instrument, it is hoped that the decomposition leads to spectral templates corresponding to single notes or percussive sounds. H will then display a representation similar to a "piano-roll" (cf. [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF]). This factorization however does not yield an effective representation of a sound presenting a noticeable spectral evolution. For instance a single note of a plucked string instrument most of the time shows high frequency components which decrease faster than low frequency components. This characteristic is not well modeled with a single frequency template. Several templates are needed which results in a less meaningful decomposition: roughly one for low frequency partials and one for high frequency partials. The meaning of each template is lost (a template no longer corresponds to a musical event such as a note).

To address this issue, we propose an extension of NMF where temporal activations become time/frequency activations. The factorization (1) becomes:

V f t ≈ Vft = R r=1 w f r h rt (f ) (3)
where the activation coefficients are now frequency dependent.

To avoid an increase of the problem dimensionality the h rt (f ) coefficients are further parameterized by means of ARMA models (section II-B). Equation ( 3) can be interpreted with the help of the source/filter paradigm: the spectrum of each frame of the signal results from the combination of filtered templates (sources). h rt (f ) corresponds to the time-varying filter associated to the source r. The decomposition thus benefits from the versatility of the source/filter model proved well suited for numerous sound objects.

B. AutoRegressive Moving Average (ARMA) Modeling

h rt (f ) is parameterized following the general ARMA model:

h ARM A rt (f ) = σ 2 rt Q q=0 b q rt e -i2πν f q 2 P p=0 a p rt e -i2πν f p 2 where ν f = f -1 2(F -1)
is the normalized frequency associated to frequency index f ∈ {1, ..., F } (as audio signal are real valued, we only consider frequencies between 0 and the Nyquist frequency). b q rt are the coefficients of the MA part of the filter and a p rt those of the AR part. σ 2 rt is the global gain of the filter: in order to avoid identifiability problems, the first coefficient of all filters is imposed to be equal to 1. For P = Q = 0, h ARM A rt (f ) no longer depends on f and the decomposition corresponds to a standard NMF with temporal activations σ 2 rt . Defining a rt = (a 0 rt , . . . , a P rt ) T and b rt = (b 0 rt , . . . , b Q rt ) T , time/frequency activations can be rewritten as:

h ARM A rt (f ) = σ 2 rt b T rt T(ν f )b rt a T rt U(ν f )a rt
where T(ν) is the (Q + 1) × (Q + 1) Toeplitz matrix with [T(ν)] pq = cos(2πν(p -q)) and U(ν) is similar to T(ν) but of dimension (P +1)×(P +1). MA only and AR only models are included by respectively taking P = 0 and Q = 0. It is worth noting that h ARM A rt (f ) is always non-negative while there exists no non-negativity constraint on b q rt or on a p rt . The parameterized power spectrogram given in equation ( 3) then becomes:

Vft = R r=1 w f r σ 2 rt b T rt T(ν f )b rt a T rt U(ν f )a rt (4)
III. ALGORITHM We choose a general β-divergence cost function:

C(W, A, B, Σ) = F f =1 T t=1 d β (V f t , Vft ) with [W] f r = w f r , [Σ] rt = σ 2 rt , [A] rtp = a p rt and [B] rtq = b q
rt and the expression of d β is given in equation ( 2). The partial derivative of the cost function wrt any variable θ (θ being any coefficient of W, Σ, A or B) is:

∂C(W, A, B, Σ) ∂θ = F f =1 T t=1 V β-2 f t ( Vft -V f t ) ∂ Vft ∂θ (5)
The expression of the gradient of C wrt a vector θ of several coefficients of A or B is the same, replacing the partial derivative by a gradient ∇ θ in [START_REF] Makino | Audio source separation based on independent component analysis[END_REF]. This leads to update rules for a multiplicative gradient descent algorithm similar to those used in [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], [START_REF] Smaragdis | Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs[END_REF]. In such an iterative algorithm, the update rule associated to one of the parameters is obtained from the partial derivative of the cost function wrt this parameter, written as a difference of two positive terms: ∂C ∂θ = G θ -F θ The update rule for θ is then:

θ ← θ × F θ G θ (6)
This rule ensures that θ remains non-negative, becomes constant if the partial derivative is zero and evolves in the opposite direction of the partial derivative (thus in the descent direction).

A. Update of frequency templates:

We derive multiplicative update rules for w f r from the expression of the partial derivative of the cost function with respect to w f r .

The partial derivative of the parameterized spectrogram defined in equation ( 4) with respect to w f0r0 is:

∂ Vft ∂w f0r0 = h ARM A r0t (f 0 )δ f0f
where δ f0f is a Kronecker delta. Then, by replacing this expression in equation ( 5) with θ = w f r0 , we obtain the partial derivative of the cost function with respect to w f0r0 :

∂C(W, A, B, Σ) ∂w f0r0 = N t=1 h ARM A r0t (f 0 ) V β-2 f0t ( Vf0t -V f0t )
This derivative is written as a difference of two positive terms:

G w f 0 r 0 = N t=1 h ARM A r0t (f 0 ) V β-1 f0t and F w f 0 r 0 = N t=1 h ARM A r0t (f 0 ) V β-2 f0t V f0t
Then the update rule of w f0r0 is:

w f0r0 ← w f0r0 F w f 0 r 0 G w f 0 r 0 (7)

B. Update of temporal activation gain:

In the same way as for w f r , we derive multiplicative update rules for σ 2 rt from the expression of the partial derivative of the cost function with respect to σ 2 rt . The partial derivative of the parameterized spectrogram defined in equation ( 4) with respect to σ 2 r0t0 is:

∂ Vft ∂σ 2 r0t0 = w f r0 b T r0t0 T(ν f )b r0t0 a T r0t0 U(ν f )a r0t0 δ t0t
where δ t0t is a Kronecker delta. Then, by substituting this expression into equation ( 5) with θ = σ 2 r0t0 , we obtain the partial derivative of the cost function with respect to σ 2 r0t0 :

∂C(W, A, B, Σ) ∂σ 2 r0t0 = M f =1 w f r0 σ 2 r0t0 h ARM A r0t0 (f ) V β-2 f t0 ( Vft0 -V f t0 )
This derivative is written as a difference of two positive terms:

G σ 2 r 0 t 0 = M f =1 w f r0 σ 2 r0t0 h ARM A r0t0 (f ) V β-1 f t0 and F σ 2 r 0 t 0 = M f =1 w f r0 σ 2 r0t0 h ARM A r0t0 (f ) V β-2 f t0 V f t0
Then the update rule of σ 2 r0t0 is:

σ 2 r0t0 ← σ 2 r0t0 F σ 2 r 0 t 0 G σ 2 r 0 t 0 (8)
We can notice that when Q = 0 and P = 0 (i.e. there is no filter), the update rules of w f r and σ 2 rt are the same as the ones given in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] which are the standard NMF rules for a βdivergence cost function where w f r corresponds to frequency templates and σ 2 rt to temporal activations.

C. Update of filters

The update rules of the coefficients of the filters are derived in a similar way, but the updates are not element-wise, but rather "vector-wise": we derive an update rule for each b rt and for each a rt .

Update of b rt : The gradient of the parameterized spectrogram Vft wrt b r0t0 is:

∇ br 0 t 0 Vft = δ t0t 2w f r0 σ 2 r0t0 a T r0t0 U(ν f )a r0t0 T(ν f )b r0t0
Then, by substituting this expression into equation ( 5) with θ = b r0t0 , we obtain the gradient of the cost function wrt b r0t0 :

∇ br 0 t 0 C = 2 F f =1 w f r0 σ 2 r0t0 V β-2 f t0 ( Vft0 -V f t0 ) a T r0t0 U(ν f )a r0t0 T(ν f )b r0t0 = 2σ 2 r0t0 R r0t0 -R ′ r0t0 b r0t0
where:

R r0t0 = F f =1 w f r0 V β-1 f t0 a T r0t0 U(ν f )a r0t0 T(ν f ) R ′ r0t0 = F f =1 w f r0 V β-2 f t0 V f t0 a T r0t0 U(ν f )a r0t0 T(ν f )
Both matrices R r0t0 and R ′ r0t0 are positive definite under mild assumptions: these matrices are clearly positive semidefinite and it can easily be shown that if there are at least Q+1 different indexes f such that w f r0

V β f t0 = 0 then R r0t0 is non-singular (for R ′ r0t0 , the assumption is very similar). This assumption is always true in practice as long as the frame with index t 0 is not a null vector (i.e. with all samples equal to 0): in this particular case, the decomposition is trivial, and the global gains σ r0t0 should be equal to 0.

Then, we follow the approach given in [START_REF] Badeau | Weighted maximum likelihood autoregressive and moving average spectrum modeling[END_REF] and derive the following update rule for the MA part of the filter:

b r0t0 ← R -1 r0t0 R ′ r0t0 b r0t0 (9)
As R r0t0 and R ′ r0t0 are both non singular, R -1 r0t0 is well defined and b r0t0 is ensured to never be zero.

Update of a rt :

The update rules of a rt are derived in the same way as for b rt . The partial gradient of the parameterized spectrogram Vft with respect to a r0t0 is:

∇ ar 0 t 0 Vft = -2δ tt0 w f r0 h ARM A r0t0 (f ) a T r0t0 U(ν f )a r0t0 U(ν f )a r0t0
Then, by substituting this expression into equation ( 5) with θ = a r0t0 , we obtain the partial gradient of the cost function with respect to a r0t0 :

∇ ar 0 t 0 C(W, A, B, Σ) = 2σ 2 r0t0 S ′ r0t0 -S r0t0 a r0t0
where:

S r0t0 = M f =1 w f r0 V β-1 f t0 b T r0t0 T(ν f )b r0t0 (a T r0t0 U(ν f )a r0t0 ) 2 U(ν f ) and S ′ r0t0 = M f =1 w f r0 V β-2 f t0 V f t0 b T r0t0 T(ν f )b r0t0 (a T r0t0 U(ν f )a r0t0 ) 2 U(ν f )
Both matrices S r0t0 and S ′ r0t0 are positive definite under mild assumptions.

Thus we derive the following update rule for the AR part of the filter:

a r0t0 ← S ′-1 r0t0 S r0t0 a r0t0 (10)

D. Description of the algorithm

The update rules ( 7), ( 8), ( 9) and ( 10) are applied successively to all the coefficients of W, all the coefficients of Σ, all the coefficients of B and all the coefficients of A. Between the updates of each of these matrices (and tensors), the parameterized spectrogram V is recomputed: as for the standard NMF algorithm, this recomputation between each update is necessary to ensure the convergence.

Identifiability: As for the standard NMF, the decomposition (4) which minimizes the cost function is not unique. To cope with identifiability issues, we impose constraints on W, Σ, B and A:

• for all r and t, we impose that b rt and a rt (considered as polynomials) have all their roots inside the unit circle. • for all r, we impose w r = 1 for some norm . . • for all r and t, we impose b 0 rt = 1 and a 0 rt = 1. Thus, at the end of each iteration of our algorithm, we transform b r,t and a r,t by replacing roots outside the unit circle by the conjugate of their inverse and accordingly adapting the gain, normalize each column of W, divide b r,t and a r,t by their first coefficient and update Σ in order not to change Vft by these modifications. All these transformations have no influence on the values of the parameterized spectrogram.

Another choice of filters normalization has been tested: rather than imposing b 0 rt = 1 and a 0 rt = 1, we can impose for all (r, t),

f b T rt T(ν f )b rt = 1 and f a T rt U(ν f )a rt = 1.
It corresponds to a power normalization and then it is more meaningful.

Our algorithm is detailed in Algorithm 1. In the remainder of the article, we will refer to this algorithm by the expression "source/filter factorization".

E. Dimensionality

Since our algorithm is a dimensionality reduction technique like NMF, one should take care of the dimension of the decomposition provided. The dimension of the original data is F T . In the standard NMF with R atoms, the dimension of the parameterized spectrogram is dim W + dim H = R(F + T ) (where dim Φ stands for the number of coefficients of matrix or tensor Φ). With our algorithm, the dimension of the parameters is:

dim W + dim Σ + dim A + dim B = RF + RT (P + Q + 1).
Thus, one should have RF + RT (P + Q + 1) ≪ F T , i.e.

b rt ← R -1 rt R ′ rt b rt a rt ← S ′-1
rt S rt a rt end for bring back roots of all filters inside the unit circle divide the coefficients of all filters by the first one normalize W update appropriately Σ end for R ≪ T and R(P + Q + 1) ≪ F , so P and Q must be small. As in practice F ≤ T , the condition to be respected is

R(P + Q + 1) ≪ F .
One should notice that our decomposition allows a significant reduction of the number of atoms R needed to accurately fit the data when the parts of the sounds present strong spectral variations. Then the total dimension of the parameters obtained with our decomposition remains comparable to the one obtained with standard NMF as will be shown in section IV.

Besides, one should notice that a large number of the coefficients of the filters are useless and therefore do not need to be retained: when the global gain of one of the filter σ rt becomes close to zero, these coefficients (b rt and a rt ) become meaningless and then are useless in the decomposition and can be removed without affecting the values of the parameterized spectrogram.

Finally, in the decomposition (4), all atoms are associated to filters of the same order, but it is also possible to implement a larger model where filters do not have the same characteristics for all atoms. This larger model is not presented in this paper for readability reasons.

F. Computational complexity

The computational complexity of one iteration of source/filter factorization depends on P , Q, R, F and T . The computational complexity of each step of the algorithm is given here:

• Computation of V: O ((P + Q) RF T ) operations.

• Update of W and Σ: O (RF T ) operations each.

• Update of B: O RT (F P + P 3 ) operations.

• Update of A: O RT (F Q + Q 3 ) operations. • Normalization/stabilization: O RT (F + P 3 + Q 3 ) op-
erations. The total complexity of a single iteration of the algorithm is then O ((P + Q) RF T ) operations. With our current implementation in Matlab, 100 iterations of our algorithm applied to a 1025 × 550 spectrogram (corresponding to a 6.5s signal sampled at f s = 22050Hz with 2048-sample-long windows and 75% overlap) with R = 10 atoms, P = 2 and Q = 2 last about 300s (on an Intel R Core TM 2 Duo E8400 @3.00GHz, with Matlab's multithreaded math libraries). In comparison, 100 iterations of standard NMF with the same spectrogram, the same parameters (R = 10, but P = 0 and Q = 0) on the same computer last about 9s: our algorithm appears to be slower than standard NMF. However, this comparison puts at a disadvantage our algorithm, which is designed to work with fewer atoms than standard NMF: in this case our algorithm deals with many more parameters than NMF. If we compare execution times with the same number of parameters, the difference is smaller: for the same spectrogram, 100 iterations of our algorithm with R = 2 atoms, P = 2 and Q = 2 (i.e. with the dimensionality of a standard NMF with R = 10) last about 60s.

About a third of the computation time is due to the computation of the roots of all filters (considered as polynomials) during the stabilization process. Some improvements could be made by considering another regularization method. The inversion of the matrices R rt and S ′ rt (a bit more than 10% of the total computation time) and the computation the frequency response of all filters (slighly less than 10% of the total computation time) are also very costly.

G. Practical implementation and choice of β

We empirically observed the monotonic decrease of the cost function and the convergence of the algorithm for 0 ≤ β ≤ 2 over a large set of tests: this decrease and this convergence are illustrated in particular cases in section IV-D.

However, the algorithm is unstable for 1 ≤ β ≤ 2: some numerical instabilities appear while poles of the filters come close to the unit circle. These instabilities are very common when β becomes close to 2 (Euclidean distance); however, the strong dynamics of audio data is better fitted when β becomes close to 0 as stated in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]. In order to avoid these instabilities, we limit the modulus of poles: the monotonic decrease of the cost function is no longer observed but it permits to avoid non desirable behavior of the decomposition (very resonant order 2 filters only permit to fit one partial).

For the examples in the next section, we chose β = 0.5 since the numerical instabilities were almost negligible and the results were more accurate than with β = 0 (Itakura-Saito divergence).

IV. EXAMPLES

In this section several experiments are presented to show that our algorithm is well adapted to decompose sounds having strong spectral variations. All the spectrograms used in these experiments are power spectrograms obtained from recorded signals by means of a short time Fourier transform (STFT).

Algorithms (standard NMF and source/filter factorization) were initialized with random values (except for the filters which were initially flat) and were run until apparent convergence. The algorithms have been rerun with 100 different initializations in order to maximize the chances to come close to a "good" minimum. Despite these different starting points, the reached solutions were similar in terms of qualitative aspect of the reconstructed spectrograms.

A. Didgeridoo 1) Description of the excerpt: In this section our algorithm is applied to a short didgeridoo excerpt. The didgeridoo is an ethnic wind instrument from northern Australia. It makes a continuous modulated sound produced by the vibrations of the lips. The modulations result from the mouth and throat configuration with the help of which the player is able to control several resonances. Figure 1(a) represents the spectrogram of the excerpt: the sound produced is almost harmonic (with some noise) and a strong moving resonance appears in the spectrogram. We can thus consider that this signal is composed of a single event encompassing spectral variations, and try to decompose it with a single atom (R = 1). The sampling rate of the excerpt is f s = 11025Hz. We chose a 1024-sample-long Hann window with 75% overlap for the STFT.

2) Experiment and results: The spectrogram of the excerpt is decomposed with a standard NMF algorithm for R = 1 atom and R = 5 atoms, and with source/filter factorization for R = 1 atom, with an order 3 AR modeling (Q = 0 and P = 3). Reconstructed spectrograms are respectively represented in figures 1(b), 1(c) and 1(d).

Although the didgeridoo is played alone in the analyzed spectrogram, the standard NMF needs many atoms to accurately decompose the power spectrogram. With 1 atom, NMF does not accurately represent the moving resonance (figure 1(b)). With 5 atoms, some spectral variations appear (figure 1(c)), but the resonance trajectory remains a bit unclear. Besides, the signal is not decomposed in a meaningful way (each atom is a part of the sound which has no perceptual meaning) and the dimensionality of the parameters is large (F R + RT = 3290).

In opposition to the standard NMF, source/filter factorization permits to accurately represent the spectral variability of the sound (figure 1(d)) with a single atom, keeping the dimensionality low (F R + T R(Q + 1) = 1093): the moving resonance of the original sound is well tracked, and the total error C is smaller than that of the standard NMF with R = 5. In this case, the decomposition is more efficient and relevant than the standard NMF.

B. Harpsichord 1) Description of the excerpt:

In this section our algorithm is applied to a short harpsichord excerpt, composed of two different notes (C2 and E♭2): first, the C2 is played alone, then the E♭2, and at last, both notes are played simultaneously. The spectrogram of the extract is represented in figure 2(a). As for most of plucked string instruments, high frequency partials of a harpsichord tone decay faster than low frequency partials. This phenomenon clearly occurs in the L-shaped spectrograms of figure 2(a). The sampling rate of the excerpt is f s = 44100Hz. We chose a 2048-sample-long Hann window with 75% overlap for the STFT. The standard NMF needs several atoms per note to accurately decompose the L-shaped power spectrograms: with only 2 atoms (1 per note played), the faster decay of high frequency content does not appear at all (figure 2(b)). With 6 atoms, the attenuation of high frequency partials appears (figure 2(c)), but each atom is a part of a note spectrum and has no real perceptual meaning.

The ARMA modeling included in our algorithm leads to a good description of the overall spectrogram shape. 2 atoms (1 per note) are enough to accurately fit the original short time spectrum: each atom is harmonic (figure 3(a)) and corresponds to one note while the decay of high frequency partials is clearly well described by the ARMA modeling (see time/frequency activations h ARM A rt (f ) in figure 3(b)). The dimensionality of the data provided by our algorithm (F R + T R(Q + P + 1) = 5704) remains lower than with a standard NMF with 6 atoms (F R+RT = 9804) and the global error C between the original and the reconstructed spectrogram is approximately the same as the one obtained with the standard NMF with R = 6.

Thus the decomposition provided by source/filter factorization seems to give a more meaningful representation of the given spectrogram than the one obtained with the standard NMF. C. Guitar with wah pedal 1) Description of the excerpt: In this section our algorithm is used to decompose a short extract of electric guitar processed by a wah pedal. The wah pedal (or wah-wah pedal) is a popular guitar effect which consists of a resonant filter, the resonant frequency of which is controlled by means of a foot pedal. This effect is named by emphasizing the resemblance with the human onomatopoeia "Wah". A single note of electric guitar processed by a moving wah pedal presents strong spectral variations and therefore cannot be well represented by a single atom in a standard NMF.

As a wah pedal is well modeled by an AR filter with two complex conjugates poles, we chose to decompose the extract with Q = 0 and P = 2. The chosen extract represented in figure 4(a) is composed of three different notes played successively (the first note is played a second time at the end of the extract). Each note can be viewed as a harmonic pattern which is filtered by a resonant filter, the resonant frequency of which varies between 400Hz and 1200Hz: this resonance clearly appears in the power spectrogram. The sampling rate of the excerpt is f s = 11025Hz. We chose a 1024-sample-long Hann window with 75% overlap for the STFT.

2) Experiment and results: As the analyzed sound presents strong spectral variations, the standard NMF needs many atoms to accurately decompose the power spectrogram. Thus one atom no longer corresponds to one note, and the decomposition does not correspond to the analysis that could With source/filter factorization, the strong spectral variations of each note can be accurately represented in the filter activations taking an order 2 AR modeling (Q = 0 and P = 2) as shown in figure 4(d). Then 3 atoms (one for each note) are enough to correctly fit the original spectrogram. Indeed, the global β-divergence error between the original and the reconstructed spectrogram obtained with source/filter factorization is approximately the same as the one obtained with standard NMF with 10 atoms; this β-divergence (obtained with source/filter factorization) is also approximately half that obtained with standard NMF with 3 atoms. Each atom is harmonic and corresponds to one note, and the resonance of the wah pedal clearly appears. The dimensionality of the representation obtained with source/filter factorization remains about half that of NMF with 10 atoms: M R + R(Q + 1)N = 4833. The decomposition provided by our algorithm distinguishes a stationary spectrum representing "average" guitar sounds (contained in W) from the non-stationary effect due to the wah pedal (described by time/frequency activations).

V. CONCLUSION AND FUTURE WORK

In this paper, we proposed a new iterative algorithm which is an extended version of the Non-negative Matrix Factorization based on a source/filter representation. We showed that this representation is particularly suitable to efficiently and meaningfully decompose non stationary sound objects including noticeable spectral variations.

In the future, this extended decomposition should be further developed to deal with small pitch variations (like vibrato), for instance using a constant-Q spectrogram like in [START_REF] Smaragdis | Sparse and shift-invariant feature extraction from non-negative data[END_REF], [START_REF] Schmidt | Nonnegative matrix factor 2-D deconvolution for blind single channel source separation[END_REF]. Besides, we plan to introduce harmonic constraints in the basis spectra following [START_REF] Vincent | Harmonic and inharmonic nonnegative matrix factorization for polyphonic pitch transcription[END_REF], [START_REF] Bertin | Enforcing harmonicity and smoothness in bayesian non-negative matrix factorization applied to polyphonic music transcription[END_REF], [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF]. Finally, we plan to investigate the introduction of continuity constraints between filters from one frame to another following the approach given in [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF], [START_REF] Bertin | Enforcing harmonicity and smoothness in bayesian non-negative matrix factorization applied to polyphonic music transcription[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF].

1 Fig. 1 .

 11 Fig. 1. Original power spectrogram of the extract of didgeridoo 1(a) and reconstructed spectrograms 1(b), 1(c) and 1(d)

  Srce/filter factorization, R = 2

Fig. 2 .

 2 Fig. 2. Original power spectrogram of the extract of harpsichord 2(a) and reconstructed spectrograms 2(b), 2(c) and 2(d)

Fig. 3 .

 3 Fig. 3. Source/filter decomposition (R = 2, Q = 1 and P = 1) of the power spectrogram of the harpsichord excerpt

3 Fig. 4 .

 34 Fig. 4. Original power spectrogram of the extract of electric guitar processed by a wah pedal 4(a) and reconstructed spectrograms 4(b), 4(c) and 4(d)

Figure 4 (

 4 b) represents the power spectrogram reconstructed from the NMF of the original spectrogram with 3 atoms and figure 4(c) with 10 atoms. With 3 atoms, NMF is not able to track the resonance of the wah pedal. With 10 atoms, the resonance appears, but the signal is not well explained (each atom is a part of a note and then has no perceptual meaning) and the dimensionality is higher (F R + RT = 8790).

  The 3 frequency templates (columns of W) obtained are represented in figure 5(a): each template is harmonic with its own fundamental frequency, thus it corresponds to a note (standard NMF with 3 atoms provides similar templates). The time/frequency activations (h ARM A rt (f )) are represented in figure 5(b): the resonance of the Wah pedal appears clearly where the notes are played. Thus the decomposition provided by our algorithm seems to give a more meaningful representation of the given spectrogram.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Source/filter decomposition (R = 3 and P = 2) of the power spectrogram of the Wah processed guitar excerpt


Student Member, IEEE Romain Hennequin

Member, IEEE Bertrand David

Beta-divergence as a subclass of Bregman divergence

Keywords: Beta-divergence, Bregman divergence, nonnegative matrix factorization

In this paper, the β-divergence is shown to be a particular case of Bregman divergence. This result makes it possible to straightforwardly apply theorems about Bregman divergences to β-divergences. This is of interest for numerous applications since these divergences are widely used, for instance in non-negative matrix factorization (NMF).

I. INTRODUCTION

D IVERGENCES are distance-like functions, widely used to assess the similarity between two objects. For instance, Kullback-Liebler (KL) divergence [START_REF] Kullback | On information and sufficiency[END_REF] is used in information theory to compare two probability distributions, and the Itakura-Saito (IS) divergence is used as a measure of the perceptual difference between spectra [START_REF] Itakura | Analysis synthesis telephony based on the maximum likelihood method[END_REF]. Generalized classes of divergences, for instance Bregman divergences, are used in pattern classification and clustering [START_REF] Banerjee | Clustering with Bregman divergences[END_REF]. In non-negative matrix factorization (NMF [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]), divergences are used as cost functions: NMF approximates an F × T non-negative matrix V with the product of two non-negative low-rank matrices:

V ≈ WH,
where the size of W is F × R and the size of H is R × T (with R < F and R < T ).

This approximation is generally quantized with a cost function to be minimized with respect to W and H. This cost function is often an element-wise divergence between V and WH [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF].

Numerous divergences are used as cost functions in NMF. Most common divergences probably are the Euclidean (EUC) distance, the KL divergence (see [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]) and the IS divergence (see [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]).

Several authors proposed generalized divergences which encompass these classical divergences:

• Csiszar's divergence [START_REF] Cichocki | Csiszar's divergences for non-negative matrix factorization: Family of new algorithms[END_REF], which is a generalization of Amari's α-divergence [START_REF] Cichocki | Non-negative tensor factorization using alpha and beta divergences[END_REF]. Both these divergences encompass the KL divergence and its dual.

• Bregman divergence [START_REF] Bregman | The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex programming[END_REF], [START_REF] Inderjit | Generalized nonnegative matrix approximations with Bregman divergences[END_REF], which encompasses the EUC distance, the KL divergence and the IS divergence.

The authors are with the Institut Telecom, Telecom ParisTech, CNRS LTCI, 75634 Paris Cedex 13, France (e-mail: romain.hennequin@telecomparistech.fr; bertrand.david@telecom-paristech.fr; roland.badeau@telecomparistech.fr).

The research leading to this paper was supported by the CONTINT program of the French National Research Agency (ANR), as a part of the DReaM project (ANR-09-CORD-006-03), and by the Quaero Programme, funded by OSEO, French State agency for innovation.

• β-divergence, introduced in [START_REF] Eguchi | Robustifying maximum likelihood estimation[END_REF] and studied as a cost function for NMF in [START_REF] Kompass | A generalized divergence measure for nonnegative matrix factorization[END_REF] which also encompasses the EUC distance, the KL divergence and the IS divergence. NMF is widely used in numerous areas such as image processing [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF], text mining [START_REF] Pauca | Text mining using non-negative matrix factorizations[END_REF], email surveillance [START_REF] Berry | Email surveillance using nonnegative matrix factorization[END_REF], spectroscopy [START_REF] Gobinet | Application of non-negative matrix factorization to fluorescence spectroscopy[END_REF] and audio processing [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF], [START_REF] Paulus | Drum transcription with nonnegative spectrogram factorization[END_REF].

In this paper, we give a proof that the β-divergence is actually a subclass of Bregman divergence. This result is supposed to be known in a certain community [START_REF] Nielsen | The dual voronoi diagrams with respect to representational bregman divergences[END_REF], [START_REF]Families of Alpha-Beta, Gamma-Divergences: Flexible, and Robust Measures of Similarities. Andrzej cichocki and shun-ichi amari[END_REF], but we propose to give a full demonstration of it in the wide framework of element-wise divergences and we present application to illustrate its interest. This result indeed permits to immediately particularize properties derived for Bregman divergence to β-divergence.

II. DIVERGENCE

In this section, we define the concept of divergence, element-wise divergence, and the particular case of Bregman divergence and β-divergence.

A. Definition

Divergences are distance-like functions which measure the separation between two elements. Definition 2.1: Let S be a set. A divergence on S is a function D : S × S → R satisfying:

∀(p, q) ∈ S × S D(p||q) ≥ 0, and D(p||q) = 0 iff p = q.
As a distance, a divergence should be non-negative and separable. However, a divergence does not necessarily satisfy the triangle inequality and the symmetry axiom of a distance. In order to avoid the confusion with distances, the notation D(p||q) is often used instead of the classical distance notation D(p, q).

B. Bregman divergence

Definition 2.2: Let S be a convex subset of a Hilbert space and Φ : S → R a continuously differentiable strictly convex function. The Bregman divergence [3] D Φ : S × S → R + (where R + is the set of non-negative real numbers) is defined as:

D Φ (x||y) = Φ(x) -Φ(y) -x -y, ∇Φ(y) ,
where ∇Φ(y) stands for the gradient of Φ evaluated at y and ., . is the standard Hermitian dot product.

The value of Bregman divergence D Φ (x||y) can be viewed as the difference between the function Φ(x) and its first order Taylor series at y. Thus, adding an affine form to Φ does not change D Φ .

III. ELEMENT-WISE DIVERGENCE

A. Definition

In this section, S = R N + or S = (R + \{0}) N . On such sets, one can define element-wise divergences: a divergence on R N + (resp. (R + \{0}) N ) is called element-wise if there exists a divergence d on R + (resp. R + \{0}) such that:

∀x = (x 1 , ..., x n ), ∀y = (y 1 , ..., y n ) D(x||y) = N n=1 d(x n |y n ).

B. Element-wise Bregman divergence

Element-wise Bregman divergences are a subclass of Bregman divergences for which Φ is the sum of N scalar, continuously differentiable and strictly convex element-wise functions:

∀x = (x 1 , x 2 , . . . , x N ) ∈ S Φ(x) = N n=1 φ(x i ).
d β (x|y) =      x y -log( x y ) -1 β = 0 x(log x -log y) + (y -x) β = 1 x β +(β-1)y β -βxy β-1 β(β-1) β ∈ R\{0, 1}.
One should notice that the previous definition of β-divergence is continuous with respect to β in the sense that:

∀β 0 ∈ R, ∀x, y ∈ R + \{0} d β0 (x, y) = lim β→β0 d β (x, y),
particularly for β 0 = 0 and β 0 = 1.

From this divergence on R + \{0}, one can define an element-wise β-divergence on (R + \{0}) N :

D β (x||y) = N n=1 d β (x|y).

IV. β-DIVERGENCE AS A BREGMAN DIVERGENCE

In this section, we show that the Bregman divergence encompasses the β-divergence in a natural way.

For β ∈ R, let φ β : R + \{0} → R be the function defined as:

∀x ∈ R + \{0} φ β (x) =    -log x + x -1 β = 0 x log x -x + 1 β = 1 x β β(β-1) -x β-1 + 1 β otherwise.
As shown in Appendix A, this definition is continuous with respect to β in the sense that:

∀β 0 ∈ R, ∀x ∈ R + \{0} lim β→β0 φ β (x) = φ β0 (x).
For all β ∈ R, φ β is smooth on R + \{0} and its second derivative is:

φ ′′ β (x) = x β-2 . ( 1
)
Thus φ β is strictly convex and one can define the Bregman divergence D φ β associated to φ β :

D φ β (x||y) = N n=1 φ β (x n ) -φ β (y n ) -φ ′ β (y n )(x n -y n ).
Straightforward calculations (see Appendix B) show that for all β ∈ R, D φ β = D β is a β-divergence. Thus the Bregman divergence encompasses β-divergence.

V. APPLICATIONS

In this section, we present examples showing how our result can particularize properties of the Bregman divergence to the β-divergence, in order to illustrate its potential fields of application.

A. Non-negative matrix factorization

The multiplicative update rule of H for minimizing a Bregman divergence D φ cost function given in [START_REF] Inderjit | Generalized nonnegative matrix approximations with Bregman divergences[END_REF] is:

H ← H. W T (φ ′′ (WH).V) W T (φ ′′ (WH).(WH)) .
The product ".", the fraction bar and φ ′′ are element-wise operations on the corresponding matrices. We can directly derive the (already well-known [START_REF] Cichocki | Csiszar's divergences for non-negative matrix factorization: Family of new algorithms[END_REF]) update rule of H for a β-divergence D β cost function using (1): 1) ) .

H ← H. W T ((WH) .(β-2) .V) W T ((WH) .(β-
This illustrates the interest of deriving general properties about the Bregman divergence instead of the β-divergence.

B. Right type centroid

The right type centroid is used in clustering as a "center" of a point cloud with respect to an asymmetric divergence: the right type centroid can thus be thought as an average typical point of a set.

Definition 5.1: Given a divergence D, the right type centroid of a finite set X = {x 1 , x 2 ..., x n } ⊂ S is defined as:

c D right = argmin

VI. CONCLUSION

In this paper, we presented a proof that the general class of Bregman divergence encompasses the β-divergence in a natural way. This results permits to straightforwardly apply theorems about the Bregman divergence to the β-divergence. As the latter is widely used in methods such as NMF, which has applications in numerous areas (signal processing, clustering, data mining, spectroscopy), the field of application of this result is quite wide.

RÉFÉRENCES APPENDIX

A. Continuity of φ β with respect to β

With the little-o notation, one can write as β → 0:

φ β (x) = e β log x β(β -1) - x β -1 + 1 β = 1 + β log x + o(β) β(β -1) - x β -1 + β -1 β(β -1) = 1 + log x -x β -1 + o(1).
Then:

lim β→0 φ β (x) = -log x + x -1.
In a similar way, one can write as β → 1:

  Then D Φ (x||y) = N i=1 d φ (x i |y i ) where d φ (x|y) = φ(x)φ(y)φ ′ (y)(xy) and thus, the divergence is element-wise. For element-wise Bregman divergences, we can equivalently denote the divergence D Φ or D φ . C. β-divergence Definition 3.1: Let β ∈ R. The β-divergence on R + \{0} is defined by:

Theorem 5 . 1 :

 51 When D = D β is a β divergence, c D β right is unique, independent of β and is equal to µ = 1 n n i=1 x i . Proof: It was shown in [1] that, when D = D Φ is a Bregman divergence, c DΦ right is unique, independent of Φ and is equal to µ = 1 n n i=1 x i .As β-divergence is a subclass of Bregman divergence, the proof is straightforward.

- 1 - 1 =

 11 φ β (x) = xe (β-1) log x β(β -1) β + 1 + (β -1) log x + o(β -1)) β(βx) = x log xx + 1.B. Equivalence between the Bregman divergence and the βdivergenceFor β ∈ R\{0, 1}:d φ β (x|y) = x β β(β -1) -x β -1 -y β β(β -1) + y β -1 -y β-1 β -1 -1 β -1 (xy) = 1 β(β -1) (x β + (β -1)y ββxy β-1 ) = d β (x|y).It is straightforward to check that the equality d φ β (x|y) = d β (x|y) also holds for β ∈ {0, 1}:d φ0 (x|y) = -log x + x -(-log y + y) -(-1 y + 1)(xy) = -log x + log y + (xy) + x y (xy) d 0 (x|y),d φ1 (x|y) = x log xx + 1 -(y log yy + 1) -log y(xy) = x(log x -log y) + (yx) = d 1 (x|y).