Iane Roulleau - soprano

After graduating in economics and management, Iane Roulleau studied singing at the Guildhall School of Music and Drama in London and then at the Bastille Opera Studio in Paris. She took part in many masterclasses under the tutelage of Suzanne Danco, Rgine Crespin, Michel Snchal and Graziella Sciutti, among others. She has won awards in several French and international competitions: she won first prize in the French contest Voix Nouvelles organised by France Telecom; second prize in the Tunbridge Wells International Young Artists Competition; second prize in the French contest Tournoi des Voix dOr, and the special prize given by the Directeurs de Thtre; first prize in the International UFAM Competition; she was also a semi-finalist in the International Voice Masters Competition of Monte Carlo. Iane has given concerts in Paris at the Bastille Opera Amphitheatre, in Tangiers, Casablanca, London, in Ramatuelle with J.-Ph. Collard and P. Meyer, as well as at many other summer festivals. Since 1995 she has performed regularly in productions of major oratorios such as the Stabat Mater by Pergolese; Exultate Jubilate, C minor and C major Masses, and the Requiem by Mozart; Magnificat by Bach; cantatas and Gloria by Vivaldi; and the Brahms German Requiem. She has sung la Musica (Orfeo - Monteverdi), Constance (Les Dialogues des Carmlites - Poulenc) with the Paris Opera Orchestra; Hro (Batrice et Bndict Berlioz) in Colmar, Mulhouse, and in the Ile de France; Nanetta (Falstaff - Verdi), Clorinda (Cenerentola - Rossini) in Lisbon; Cupidon-Amitta (Der Geduldige Socrates Telemann) conducted by J.C. Malgoire in Warsaw and Magdeburg, broadcast by France Musique in November 1997; Les Noces by Stravinski conducted by M. Janowski with the Radio France Philharmonic Orchestra; Susanna (Le Nozze di Figaro - Mozart) on tour in France; Serpina (La Serva Padrona - Pergolese), and Frasquita (Carmen - Bizet) in Nancy and Marseille Opera Houses and with the Iceland Symphony Orchestra in Reykjavik; Le Martyr de St Sbastien by Debussy conducted by G.Prtre in La Scala, Milan; Elsbeth (Fantasio - Offenbach) in Rennes, Angers, Tours und Nantes Opera Houses, a comic opera broadcast by France Musique in November 2000; the first Lady (The Magic Flute - Mozart) in Paris with the Gruss circus, broadcast on French television in December 2001. In July 2003 she sang in Micromgas by P. Mfano at the Montpellier Festival. In December 2003 she sang Berta (Il Barbiere di Siviglia - Rossini) with the Opra de Chambre de France in Menton. In January/February 2004 she sang in A Midsummer Nights Dream by Mendelssohn with the Lamoureux Orchestra conducted by Y. Sado in the Champs Elyses Theatre, the Caen Theatre and during the French event Les Folles Journes in Nantes. She has recorded the complete work of Germaine Tailleferres Melodies for Cambria with the French baritone Mario Hacquard, and took part in the recording of LArlsienne by Bizet conducted by J.C. Malgoire for Audivis. She recently took part in the recording of Micromegas by P. Mefano with 2E2M Ensemble for Assa. http://www.ianeroulleau.com/ contact e-mail : iroulleau@hotmail.com

PAPERS

Long Interpolation of Audio Signals Using Linear Prediction in Sinusoidal Modeling*
MATHIEU LAGRANGE AND SYLVAIN MARCHAND
(lagrange@labri.fr) (sylvain.marchand@labri.fr)
LaBRI, Universite Bordeaux 1, F-33405 Talence Cedex, France AND JEAN-BERNARD RAULT
(jeanbernard.rault@francetelecom.com)
France Telecom R&D, F-35512 Cesson Sevigne cedex, France
Within the context of sinusoidal modeling, a new method for the interpolation of sinusoidal components is proposed. It is shown that autoregressive modeling of the amplitude and frequency parameters of these components allows us to interpolate missing audio data realistically, especially in the case of musical modulations such as vibrato or tremolo. The problem of phase discontinuity at the gap boundaries is also addressed. Finally, an original algorithm for the interpolation of a missing region of a whole set of sinusoids is presented. Objective and subjective tests show that the quality is improved significantly compared to common sinusoidal and temporal interpolation techniques of missing audio data.
0 INTRODUCTION The sinusoidal model [1], [2] provides a high-quality representation of pseudostationary sounds. Therefore this model is used widely for many musical audio processing purposes such as musical sound processing [3][5] and audio coding [6], [7]. Parameters of the sinusoidal model are extracted from the original sound in a frame-based manner, and a sound that is close to the original one can be synthesized from the extracted parameters. The problem of missing information about sinusoids can occur at both sides of the sinusoidal analysis and synthesis procedure. During the analysis some gaps in the original signal may have been introduced by another module, for example, a module of detection and removal of clicks or transients. During the synthesis, sinusoidal parameters may not be available. For example, in a stream-based audio coding application, some frame packets may be unavailable at the time they are needed for the synthesis. In both cases, information about the sinusoids is available before and after the gap and can be exploited to interpolate the evolution of the partials within the missing region. Let a gap start at frame index n1 and end at frame index n2, corrupting a set of sinusoids S. The aim of the algorithm described in this paper is to interpolate S during the gap. As shown in Fig. 1, the set B represents sinusoids existing before the gap and ending at frame n1. The set A represents sinusoids existing after the gap and beginning at frame n2. Only the sinusoids of these two sets will be considered for the interpolation of the gap.
*Manuscript received 2004 December 6; revised 2005 July 28.
J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
The block diagram in Fig. 1(a) describes the four-step algorithm used to interpolate the missing region. The predicted frequencies and amplitudes in the missing region are computed for each sinusoid of the two sets [Fig. 1(b)]. According to these predicted parameter sets B and A, some sinusoids of B are matched to sinusoids of A. These matched sinusoids then become sinusoids with a missing region [dashed lines in Fig. 1(c)]. This missing region is interpolated using the predicted parameters of the two matched sinusoids. Next, unmatched sinusoids (terminating or beginning with open dots) are extrapolated in the missing region according to their predicted parameters using a specific technique. The interpolated set of sinusoids S is plotted in Fig. 1(c). The remainder of this paper is organized as follows. The sinusoidal model and the limitation of existing interpolation methods are presented in Section 1. The use of autoregressive (AR) modeling for the prediction of the amplitude and frequency parameters of a sinusoid in a missing region is presented in Section 2. Section 3 describes the matching of sinusoids from both sides of the missing region and introduces the use of the predicted parameters to enhance the matching of modulated sinusoids. Next an original method for interpolating the missing parameters of a partial is introduced in Section 4 and is followed by objective and subjective evaluations of this interpolation method. The extrapolation of unmatched sinusoids is presented in Section 5. Finally an algorithm for the interpolation of a whole set of sinusoids in a missing region that makes use of these concepts is compared in Section 6 to known sinusoidal and temporal techniques.

LAGRANGE ET AL.

1 SINUSOIDAL MODELING Sinusoidal modeling aims at representing a sound signal as a sum of sinusoids of given amplitudes, frequencies, and phases. For stationary pseudoperiodic sounds these amplitudes and frequencies evolve slowly and continuously with time, controlling a set of pseudosinusoidal oscillators commonly called partials. (This term will be preferred to sinusoid during the remainder of this paper.) The audio signal s can be calculated from the additive parameters using Eqs. (1) and (2), st =

cos pt

(1) (2)

pt = p0 + 2

f u du
where P is the number of partials and the functions fp, Ap, and p are the instantaneous frequency, amplitude, and phase of the pth partial, respectively. The P triplets (fp, Ap, p) are the parameters of the additive model and represent points in the frequencyamplitude plane at time t. Although potential applications are numerous, few people have paid attention to the interpolation issue. Quatieri and Danisewicz [8] propose an algorithm to interpolate overlapping harmonics for the purpose of separating two speech signals. The amplitude is interpolated linearly, and cubic interpolation is used for the phase. The fre-
quency can be found by the differentiation of the cubic phase polynomial. Although this strategy was originally designed for intraframe parameter interpolation for synthesis purposes [1], this method shows good results for gaps of lengths from 20 to 100 ms during stationary regions of speech sounds. Later on Maher [9] proposed an algorithm to interpolate a whole set of sinusoids with an approximation of missing audio data based on the same principles. This interpolation method based on a polynomial interpolation of the parameters of the partials preserves the harmonic relation among partials together with the envelope of the sound. Yet modulations of the parameters of the partials are not taken into account. For example, the frequency of a partial having natural vibrato is a sinusoid in the timefrequency plane of about 4-Hz frequency. Since the phase polynomial is cubic, the resulting interpolation of the frequency is a quadratic polynomial. A sinusoid is approximated correctly by a quadratic polynomial for less than a quarter of a period. The use of such an interpolation scheme for frequency and phase parameters is limited to segments up to 60 ms. Similarly, if we want to handle natural tremolo, the use of linear interpolation is limited to segments of up to 20 ms. According to Bregman [10] these modulations should be considered, because such modulations play an important role in sound perception: Small fluctuations in frequency occur naturally in the human voice and in musical

where K is the order of the LP model. We are then looking for a vector aK that minimizes the power of the prediction error, E= xn xn

n=1 N 2

Supposing that a vector aK minimizing the power of the prediction error of the frequencies of Pi is known, the frequencies and amplitudes of Pi are computed by infinite impulse response filtering of the frequencies and amplitudes of Pi [see Fig. 5(b)]. The same strategy is applied to the partial Pj, except that the two extrapolations are done backward [see Fig. 5(c)]. As it will be demonstrated, this extrapolation scheme is able to preserve the modulations of the parameters of the partials in the missing region. However, the predictions of the frequencies of the partials of a harmonic source are computed separately. The proposed prediction scheme also preserves harmonicity provided that the partials of B and A are estimated correctly. Let us consider a set of partials with harmonically related frequencies. The fundamental is denoted by P0 and the harmonics by Pr, with r > 0. The frequencies of the harmonics verify Fr(n) ( r + 1)F0(n). (11)
where li and lj are the lengths of Pi and Pj respectively, and Pk(n) is the triplet of instantaneous parameters of the partial Pk at frame n. Let Pi and Pj denote the predicted amplitude and frequency of the partials of the B and A sets during the missing region, Pi = Pin1 + k, k = 1, , n2 n Pj = Pjn2 k, k = 1, , n2 n f Pkn = kn, Akn, for all k (6) (7) (8)
To predict the evolution of the frequencies of these partials, we consider LP coefficients for each harmonic ar (h) computed using Fr(n) as observations. Because of K Eq. (11) and the scale invariance of LP coefficients [16], we have ar a0. Thus the harmonicity constraint is K K preserved, Fr(n) ( r + 1) F0(n). (12)
where Pk(n) is a couple of instantaneous predicted parameters since the phase will not be predicted, but deduced from the frequency. These parameters should be computed using a relevant method, chosen according to the characteristics of the evolutions of the amplitude and frequency parameters. These evolutions in the timefrequency and timeamplitude planes can be constant, increasing or decreasing exponentially (portamento in the timefrequency plane) or sinusoidal (vibrato in the timefrequency plane and tremolo in the timeamplitude plane).
2.1 Linear Prediction Methods The challenge in linear prediction is to choose a wellsuited method to minimize the error E, given N past samplesconsidered as observationsand the model order K. In this section three methods are described out of many: the autocorrelation method, the covariance method, and the Burg method. Only the method retained is detailed so that it can be implemented easily; the reader is invited to refer to [17], [16] for a complete description of the others. The choice among these three methods is driven by specific constraints: only few observed samples are available and the estimated LP coefficients have to be suitable for extrapolation. The autocorrelation method minimizes the forward prediction error power on an infinite support. In practice the signal is finite. Samples of the x(n) process that are not observed are then set to zero, and observed samples are

windowed in order to minimize the discontinuity at the boundaries. As a consequence, this method requires N > > K to be effective. Alternatively, the LP coefficients can be estimated on a finite support with the covariance method. This method minimizes the forward prediction error power on a finite support. Since no zeroing of the data is necessary, this method is a good candidate for coefficient estimation of a process having few observed samples. Unfortunately this method should be avoided for data extrapolation because it can lead to filters that are not minimum phase, that is, the estimated poles are not guaranteed to lie within the unit circle. Let ef (n) and eb(n) denote, respectively, the forward and k k backward prediction errors at an intermediate order k, ef n = xn + k
x and y are vectors. Using the Schwarz inequality, it is verified that rk has a magnitude lower than 1. With the Burg method the minimization is done on a finite support and the joint minimization of the forward and backward errors leads to a stable filter. This method is then suitable for data extrapolation with few observed samples. The following algorithm computes the vector a of LP coefficients at order K using the Burg method, ef x eb x a1 for k from 1 to K do efp ef without its first element ebp eb without its last element rk 2ebp efp ebp ebp + efp efp ef efp + rkebp eb ebp + rk efp a a0, a1, , ak, 0 + rk0, ak, ak 1, , a0 end for. 2.2 Linear Prediction Parameters The number of observed samples used to estimate the LP coefficients has to be large enough to be able to extract the signal periodicity, and short enough not to be too constrained by the past evolution. In our system the short-term analysis module uses a sliding timefrequency transform with a hop size of 360 samples on sound signals sampled at CD quality (44.1 kHz). This means that the frequency and amplitude trajectories are sampled at 120 Hz. Since we want to handle natural vibrato with a frequency of about 4 Hz, we need at least 30 samples to get the period of the vibrato. For frequency and amplitude evolutions, since we want to model exponentially increasing or decreasing evolutions (portamento) and sinusoidal evolutions (vibrato, tremolo), the order of the LP model should not be below 2. Most modulations are more complex than the sinusoidal behavior of vibrato or tremolo, thus the order should be set at a higher value. The LP coefficients used to compute the predicted pa rameters Pi and Pj are estimated using the Burg method. This method jointly minimizes the forward and backward prediction errors defined by Eqs. (16) and (17). As a consequence the number of observed samples must be at least twice the model order. In the experiments presented here, N is chosen as the minimum value between 40 and li or lj, respectively, and the model order m is set to the integer value closest to N/2. 3 MATCHING PARTIALS FROM BOTH SIDES OF THE MISSING REGION The first step to interpolate corrupted sinusoidal data in the missing region is to decide which partial of B should be linked to which partial of A to form a unique partial. The problem of matching partials from both sides of the missing region is shown in Fig. 2. The time interval is so long that the evolution of the partials within the missing

(13) (14)

ebn = xn k + k

a hxn k + h.

The Burg method minimizes the average of the forward and backward error power on a finite support in a recursive manner. That is, to obtain ak we minimize 1 k = f + p 2 k k where f = k b = k and a h + rkak1k h, h = 1, 2, , k 1 akh = r k1 h=k k, 1 N k 1 N k (15)

f k n=k Nk1 b k n=0

(16) (17)
rk being the reflection coefficient. By substituting Eq. (18) in Eqs. (16) and (17) we find a recursive expression for the forward and backward errors, ef n = ef n + rkeb n 1 k k1 k1 ebn = eb n 1 + rkef n k k1 k1 where ef n = ebn = xn. (21) (19) (20)
To find rk we differentiate the kth prediction error power with respect to rk, and by setting the derivative to zero we obtain rk = 2

N1 f 2 n=k ek1n

N1 f b n=k ek1n ek1n

. + eb n 12 k1

The minimum-phase property is ensured because the expression of rk is of the form rk 2 xy/(|x|2 + |y|2), where
region has to be taken into account to achieve a good match. We propose that this decision step can be done using predicted information (Pi and Pj) computed using the method introduced in the previous section. This issue is quite similar to the partial tracking problem, but with a much longer time interval between elements to be linked. First a straightforward adaptation of the partial tracking algorithm proposed in [1] is discussed. It will be used in Section 6 for comparison purposes. Couples of partials (Pi, Pj) such that the distance between the last frequency of Pi and the first frequency of Pj is below a given threshold f are matched, |fi(n1) fj(n2)| < f (23)
The normalized Euclidean distance dA (Pi, Pj) between the predicted amplitude is defined similarly. Each couple of partials (Pi, Pj) such that df (Pi, Pj) is below a given threshold f is a candidate for matching. Next these candidates are considered in increasing df distance order. The candidate partials are effectively matched if two criteria involving predicted frequencies and predicted amplitudes are satisfied. These criteria are defined as dfPi,Pj Tf 1 + + f f
where fi(n1) is the last frequency of Pi and fj (n2) is the first frequency of Pj, and f is a threshold parameter in hertz. Yet if the spectrum is changing within the gap interval, this approach may be unsatisfactory, as explained in [9] and shown in Fig. 4(a). Considering that the parameters of the partials have a predictable evolution is useful to match the partials of the two B and A sets more reliably. Unfortunately, considering a simple Euclidean distance between the two predictions in frequency or amplitude may lead to difficulties. If the two predictions vary a lot, the thresholding procedure should be more tolerant than if the two predictions are nearly constant (see Fig. 3). To cope with this problem, a Euclidean distance between the two predictions normalized by the sum of the standard deviation of the two predictions is used to decide whether or not partials from both sides of the missing region should be matched. Let df (Pi, Pj) denote the normalized Euclidean distance f between the predicted frequencies i and j, f dfPi, Pj =

dAPi, Pj Ta 1 + A + A

where (x) is the standard deviation of the vector x, and Tf and Ta are threshold parameters in frequency and amplitude. If these criteria are met for a couple (Pi, Pj), the two partials of the couple are merged in a unique partial Pm, and each couple where Pi or Pj appears is removed from the sorted list. The missing region of the resulting partial Pm is interpolated using the method described in the next section. This process iterates until no satisfactory couple remains. Using this algorithm, the matching is performed even in modulated cases (see Fig. 4) without spurious link in stationary cases [see Fig. 3(b)]. Finally unmatched partials are extrapolated in the missing region using an algorithm described in Section 5. 4 INTERPOLATING THE MISSING INFORMATION WITHIN A PARTIAL Let a couple (Pi, Pj) be represented as a unique partial Pm. The interpolated frequency and amplitude parameters

n21 n=n1+1 f in

jn2 f. n2 n1 1
Fig. 2. Matching partials from both sides of missing region.
J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October 895
of Pm, starting at n1 + 1 and ending at n2 1, are computed by mixing the predicted frequency and amplitude param eters Pi and Pj. The phase continuity at the boundaries of the missing region is then ensured by a method described at the end of this section. 4.1 Frequency Interpolation To compute m(n) given the two predicted frequenf cies i(n) and j(n), a crossfading is carried out by multif f f plying i(n) by a window function w and j(n) by 1 w, f n = w fm
As can be seen in Fig. 5, the forward prediction is of better quality than the backward one, since here Pi is longer than Pj. In general the window function w(t) used to crossfade the two predictions should then be asymmetric in order to favor the prediction done with the largest data set. The symmetric cosine window computed using Eq. (28) is equal to 0.5 in the middle of the missing region, ct = 1 + cos1 + t. 2 (28)
n n1 n n1 n + 1 w fi n2 n1 n2 n1

n. fj (27)

The symmetric crossfading done using this window function is relevant only if the two partials Pi and Pj have the
Fig. 3. Predictions of partials from both sides of missing region. (a) Trombone tone with glissando. (b) Transition between two piano tones.

896 J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
same length. If Pi is three times longer than Pj, the window should reach the 0.5 value at 3/4 of the missing region (see Fig. 6). As a consequence, the window function must fulfill the following constraint: w
This factor is computed according to li and lj, the respective lengths of Pi and Pj. The asymmetric window function is then

li 1 =. li + lj 2

ctr li,lj, ctr lj,li,

li lj otherwise

We propose that such an asymmetric crossfading can be done using an asymmetric factor; rx, y = loglogcx x + y (30)
with t [0, 1]. 4.2 Amplitude Interpolation The amplitude of a partial is often much more modulated than its frequency, as in speech signals. Even if mi-
where log is the Neperian logarithm.
Fig. 4. Results of matching process. (a) Reference method. (b) Proposed method.
J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October 897
cromodulations of the amplitude parameter are preserved, the long-term prediction is not satisfactory. Before the crossfade the amplitude prediction of the partial Pi is constrained to end at a given amplitude equal to the mean amplitude of the partial Pj computed from frame n2 to frame min(n2 + M, n2 + lj 1). The parameter M should be chosen so as to get an energy estimate of the beginning of partial Pj. In the configuration presented in
Section 2.2, M is set to 30. Such a constraint is fulfilled by adding to the predicted amplitude Ai an increment i(n) defined as

n n1 in = n2 n1

minM,lj1 =0

Ajn2 +

minM, lj 1 + 1

Ain2.

Fig. 5. Interpolating frequencies of a partial of a saxophone with vibrato using AR modeling. Forward prediction; backward prediction with LP formalism; predictions crossfaded using an asymmetric window favoring the more reliable predicted samples (those of the forward prediction in this case). (a) Frequencies represented by dots are unavailable.
Fig. 6. Three crossfading windows computed using Eq. (31). From left to right, windows are computed with li/lj {2/1, 3/1, 9/1}.
898 J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
The same strategy is applied to Aj by adding an increment j(n) computed as follows: n2 n jn = n2 n1

minM,li1 =0

minM, li 1 + 1

Ajn1.

compared to the signal synthetized from the original sinusoidal representation. 4.4 Objective Evaluation We simulate a missing region in the sinusoidal representation S by deleting parameters of the partials existing before and after the missing region. The other partials are left as they are, as illustrated in Fig. 7. Missing parameters of partials are then interpolated during the missing region using the polynomial or the LP-based interpolation scheme. In all the experiments reported here the interpolation scheme described in [1] is used for the intraframe interpolation of the parameters of the partials. The amplitude is interpolated linearly and phases are computed using a maximally smooth cubic polynomial. The reconstruction signal-to-noise ratio (R-SNR) is used to evaluate the performance of the algorithm tested, R-SNR = 10 log10
The corrected amplitudes are then asymmetrically crossfaded to provide the interpolated amplitude, Amn = w n n1 Ain + in n2 n1 n n1 + 1w Ajn + jn. n2 n1
4.3 Phase Interpolation Using the interpolation strategy described in [8], the phase of a partial Pm is interpolated using a maximally smooth cubic polynomial having four constraints at the boundaries fm(n1), m(n1) and fm(n2), m(n2). The interpolated frequencies are then obtained by phase differentiation. Inversely we propose to integrate the predicted frequencies m using the trapezoidal method. A phase increment f defined below is added to each phase in the missing region to ensure phase continuity at boundaries. Let us denote by (n) the unwrapped phase at frame n, [(n) ( n) mod 2]. The subscript m is omitted for convenience. In a first approximation the missing phases may be computed from f n1 + 1 = n1 + T f n1 + n1 + 1 n = n1 + 1 + T (35) (36)

M1 m=0

xm xm2

f( )1 f()]

where x(m) is the original temporal signal and x(m) the synthesized signal of the sinusoidal representation interpolated using one of the two tested interpolation strategies. For every gap size the result plotted in Fig. 8 is the mean R-SNR for every position of the gap. The LP-based interpolation is designed for musical modulation management (vibrato, tremolo) and therefore performs better for the saxophone tone or the vibraphone tone and performs as well as the polynomial method in the stationary case, such as for the harpsichord tone. 4.5 Subjective Evaluation The two methods are compared by a subjective test performed at France Telecom R&D with ten experts in audio processing. Four audio signals were used: a saxophone tone, a vibraphone tone, a soprano female voice, and an orchestra piece. It tests the interpolation for gap sizes from 80 to 820 ms. For every gap size and audio file, the experts were asked to listen to the original set of partials synthesized as an explicit reference signal. After this first listening, they were asked to note four versions, one with no interpolation performed, one with interpolation performed using the polynomial approach, one with the interpolation performed using the LP approach, and the original set of partials synthesized as a hidden reference.

region to estimate two sets of 1000 LP coefficients using the Burg method. The predictions obtained by filtering are crossfaded using the window computed using Eq. (28). The two other methods are based on sinusoidal modeling. First, two sets of partials (B and A) are extracted using the sinusoidal analysis technique described in [15]. The interpolated set of partials S is computed using one of the two sinusoidal schemes and synthesized. With the polynomial method the matching of partials is done according to Eq. (23) with f 04 Hz. The missing phases and frequencies are computed using the maximally smooth cubic phase polynomial whereas the amplitude is interpolated linearly. Extrapolated parameters of unmatched partials are computed using the algorithm detailed in Section 5, considering constant frequencies and amplitudes as predicted parameters. With the proposed method the matching is done using the algorithm described in Section 3 with Tf 5.0 and Ta .1.0 Interpolation of the missing parameters of the partials uses the method described in Section 4, and the extrapolated parameters of the unmatched partials are computed using the algorithm detailed in Section 5 with lB n2 nand lA ( n2 n1 1)/2. Five audio signals are used: a violin tone with vibrato, a piano tone, an orchestra piece, a gong tone, and the recording of two female soprano voices. The gap can be at a sustained or at a transitional segment of the sound. LP-based temporal interpolation has proven successful for the interpolation of up to thousands of samples from CD-quality audio signals without audible distortion [13], [18]. The interpolation quality of longer gaps depends on the characteristics of the signal. If it consists of stationary partials like in the piano tone, the attenuation phenomenon is lightly pronounced [see Fig. 10 (a), left]. Yet if the interpolated signal has roughly the same number of partialsaround tenwith vibrato, the attenuation is very pronounced [see Fig. 10 (a), right]. This attenuation problem explains why the marks obtained by this method range from 30 to 50 when the parameters of the partials are modulated (see Fig. 11). The sinusoidal model can be used to cope with this attenuation problem [see Fig. 10(b), (c)]. The sinusoidal interpolation scheme based on polynomial interpolation outperforms the temporal method for gap sizes up to 320 ms (see Fig. 11). In counterpart, all kinds of modulations disappear. This effect is perceived by the listeners as a freezing of the sound throughout the interpolated region. For larger gaps linear interpolation gives an artificial interpolation, rated poorly by the listeners. The rating can be even worse than the one obtained by the temporal method. This is the case for the interpolation of a 820-ms gap of the violin tone (see Fig. 11). The proposed method keeps the advantages of the two previous methods while avoiding some of their disadvantages. Use of a sinusoidal model avoids the problem of attenuation as long as long-gap interpolation can be achieved. In addition AR modeling of the parameters of the partials is useful to preserve the modulations important to perception. The gong tone and the two soprano voices have partials with small-range modulations. The violin tone with vibrato has a larger range of frequency modulations.

For all these sounds the ratings go from 90 to 70 in a regular decay for gap sizes from 320 to 820 ms. The soprano voices can even be interpolated during 1.6 s with a good mark. The partials extracted from the orchestra piece have complex modulations because they represent harmonics of several notes and noise. The prediction capability is then lower than in the previous cases, but a fair quality can be achieved for gap sizes up to 450 ms. If the gap occurs during a transition, some important information is lost and the quality of interpolation is lower, (see Fig. 12). In this case the temporal scheme seems to be better appreciated, probably because of the attenuation effect that simulates a fade in/fade out centered at the middle of the gap. Concerning sinusoidal interpolation schemes, the quality is improved by the use of the matching algorithm presented in Section 3, which is useful to avoid a mismatch of partials of different tones. 7 CONCLUSION In this paper an enhanced method is proposed for the interpolation of audio signals based on linear prediction in sinusoidal modeling. It is shown that AR modeling of the parameters of the partials allows those partials to be interpolated reliably. Partials having simple modulations such as vibrato or tremolo allow high-quality interpolation for gap sizes up to 1 s. More complex modulations are harder to interpolate, but the proposed method shows a significant improvement over the polynomial method. Since these modulations are important to perception [10], the sinusoidal interpolation of missing audio data is more realistic. The listening tests showed that the proposed method provides fair interpolation for complex polyphonic signals for gap sizes up to 450 ms and good interpolation for monophonic modulated tones for gap sizes up to 1600 ms. 8 REFERENCES [1] R. J. McAulay and T. F. Quatieri, Speech Analysis/ Synthesis Based on a Sinusoidal Representation, IEEE Trans. Acoust., Speech, Signal Process., vol. 34, pp. 744754 (1986). [2] J. O. Smith and X. Serra, An Analysis/Synthesis Program for Non-Harmonic Sounds Based on a Sinusoidal Representation, in Proc. Int. Computer Music Conf. (ICMC) (Computer Music Assoc., San Francisco, CA, 1987). [3] K. R. Fitz and L. Haken, Sinusoidal Modeling and Manipulation Using Lemur, Computer Music J., vol. 20, no. 4, pp. 4459 (Winter 1996). [4] X. Serra, Musical Sound Modeling with Sinusoids plus Noise, in Musical Signal Processing, Studies on New Music Research ser. (Swets & Zeitlinger, Lisse, The Netherlands, 1997) pp. 91122. [5] S. Marchand and R. Strandh, InSpect and ReSpect: Spectral Modeling, Analysis and Real-Time Synthesis Software Tools for Researchers and Composers, in Proc. Int. Computer Music Conf. (ICMC) (International Computer Music Assoc., Beijing China, 1999, Oct.), pp. 341344.

[6] H. Purnhagen and N. Meine, HILNThe MPEG-4 Parametric Audio Coding Tools, in Proc. IEEE Int. Symp. on Circuits and Systems (ISCAS 2000), vol. 3 (2000 May), pp. 201204. [7] B. den Brinker, E. Schuijers, and W. Oomen, Parametric Coding for High-Quality Audio, presented at the 112th Convention of the Audio Engineering Society, J.
Audio. Eng. Soc. (Abstracts), vol. 50, p. 510 (2002 June), convention paper 5554. [8] T. F. Quatieri and R. G. Danisewicz, An Approach to Co-channel Talker Interference Suppression Using a Sinusoidal Model for Speech, IEEE Trans. Acoust., Speech, Signal Process., vol. 38, pp. 5669 (1990 Jan). [9] R. C. Maher, A Method for Extrapolation of Miss-
Fig. 10. Temporal representations of piano tone and violin tone with vibrato, interpolated using three methods tested during 820 ms. Two vertical lines fix boundaries of missing region; two symmetric lines inside this region approximate envelope of original sound. (a) Temporal interpolation. (b) Polynomial-based sinusoidal interpolation. (c) LP-based sinusoidal interpolation.
J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October 903
ing Digital Audio Data, J. Audio Eng. Soc. (Engineering Reports), vol. 42, pp. 350357 (1994 May). [10] A. S. Bregman, Auditory Scene Analysis: The Perceptual Organization of Sound (MIT Press, Cambridge, MA, 1990). [11] A. J. E. M. Janssen, R. N. J. Veldhius, and L. B.
Vries, Adaptive Interpolation of Discrete-Time Signals that Can Be Modeled as Autoregressive Processes, IEEE Trans. Acoust., Speech, Signal Process., vol. 34, pp. 317330 (1986). [12] W. Etter, Restoration of a Discrete-Time Signal Segment by Interpolation Based on the Left-Sided and
Fig. 11. Results of listening tests comparing polynomial-based method (), LP-based method (), and temporal method () for three gap sizes. Symbolsmeans of votes, linesconfidence intervals for each method.
Fig. 12. Results of the listening tests comparing polynomial-based method (), LP-based method (), and temporal method () for three gap sizes on a transitional segment of piano. Symbolsmeans of votes, linesconfidence intervals for each method.
904 J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
Right-Sided Autoregressive Parameters, IEEE Trans. Acoust., Speech, Signal Process, vol. 44, pp. 11241135 (1996). [13] I. Kauppinen, J. Kauppinen, and P. Saarinen, A Method for Long Extrapolation of Audio Signals, J. Audio Eng. Soc., vol. 49, pp. 11671180 (2001 Dec.). [14] M. Lagrange, S. Marchand, M. Raspaud, and J. B. Rault, Enhanced Partial Tracking Using Linear Prediction, in Proc. Digital Audio Effects (DAFx) Conf. (Queen Mary, University of London, 2003 Sept.), pp. 141146. [15] M. Lagrange, S. Marchand, and J. B. Rault, Using Linear Prediction to Enhance the Tracking of Partials, in

Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP), vol. 4 (2004 May), pp. 241244. [16] S. M. Kay, Autoregressive Spectral Estimation: Methods, in Modern Spectral Estimation, Signal Processing ser. (Prentice-Hall, Englewood Cliffs, NJ, 1988), pp. 228231. [17] J. Makhoul, Linear Prediction: A Tutorial Review, Proc. IEEE, vol. 63, pp. 561580 (1975 Nov.). [18] I. Kauppinen and K. Roth, Audio Signal ExtrapolationTheory and Applications, in Proc. Digital Audio Effects (DAFx) Conf. (University of the Federal Armed Forces, Hamburg, Germany, 2002 Sept.), pp. 105110.

THE AUTHORS

M. Lagrange

S. Marchand

J.-B. Rault
Mathieu Lagrange was born in Caen, France, in 1978. He studied computer science at the University of Rennes 1, France, where he obtained his masters degree in 2000. He received a postgraduate diploma with a focus on spectral sound synthesis from the University of Bordeaux 1, Talence, France. Dr. Lagrange carried out research on sound analysis and coding at the France Telecom Laboratories in partnership with the LaBRI (computer science laboratory), University of Bordeaux 1, there he received a Ph.D. degree in 2004. He is particularly involved in spectral sound analysis, audio restoration, and auditory scene analysis. He is a member of SCRIME (Studio de Creation et de Recherche en Informa tique et Musique Electroacoustique) at the University.
Dr. Marchand was appointed associate professor at the LaBRI (computer science laboratory), University of Bordeaux 1, in 2001. He is particularly involved in spectral sound analysis, transformation, and synthesis. He is a member of SCRIME (Studio de Creation et de Recherche en Informatique et Musique Electroacoustique) at the University.
Sylvain Marchand was born in Pessac near Bordeaux, France, in 1972. He studied computer science at the University of Bordeaux 1, Talence, France. He obtained his masters degree in 1995 and a postgraduate diploma in algorithmics the following year. In the meantime he carried out research in computer music and sound modeling. He received a Ph.D. degree in 2000.
Jean-Bernard Rault received a Ph.D. degree in signal processing and telecommunications from the University of Rennes, France, in 1987. Dr. Rault then joined the CCETT in Rennes, France, to collaborate on the European project Eureka 147 (DAB) in the area of digital audio compression. From 1990 to 1992 he spent two years at Thomson-LER, there he was involved in the Multicarrier Digital Modulation studies. Since 1993 he has been a France Telecom representative with ISO/MPEG and participates in the development of the MPEG Audio coding standards. He has also been involved in several European projects (MoMuSys, Cinenet, Nadib, Song, Ardor) to contribute to audio-related work packages.