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

Fig. 1. (a) Block diagram of proposed interpolation algorithm. (b) Left, original set of sinusoids S and right, sets of sinusoids (A and B) used during the interpolation process. (c) Interpolated set of sinusoids S.
892 J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
INTERPOLATION OF AUDIO SIGNALS USING LINEAR PREDICTION
instruments. The fluctuations are not often very large, ranging from less than 1 percent for a clarinet tone to about 1 percent for a voice trying to hold a steady pitch, with larger excursions of as much as 20 percent for the vibrato of a singer. Even the smaller amounts of frequency fluctuation can have potent effects on the perceptual grouping of the components harmonics. Although Bregman only talks about frequency modulations, amplitude modulations are important too. Comments of the experts who performed the listening test presented in Section 4.5 confirmed this assertion. The missing region interpolated using the polynomial scheme was perceived as many simple tones and not as a complex one. As a result, the interpolated part was perceived as artificial. To achieve a more natural interpolation, one needs an interpolation method able to preserve these modulations of the frequency and the amplitude of partials in the missing region. Linear prediction has proven successful for digital audio restoration [11]. Given the AR modeling of parts of the signal before and after the degradation, linearly predicted extrapolations can be added to interpolate the degraded part of the signal (see [12], [13] for further details). Considering that evolutions of the amplitude and frequency parameters of the partials are time signals tooalthough with a much lower sampling ratea similar strategy can be used for the extrapolation and interpolation of amplitudes and frequencies of the partials. 2 PREDICTING EVOLUTION OF THE PARTIALS Let Pi and Pj denote partials of the B and A sets, respectively, Pi = Pin, n = n1 li + 1, , n1 Pj = Pjn, n = n2, , n2 + lj 1 Pkn = fkn, Akn, kn, for all k (3) (4) (5)

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.
where n [n1 + 2, n2] and T is the hop size in seconds. However, a phase discontinuity may occur at the end of the missing region: (n2) (n2). Let e denote the error of the phase extrapolation at n2, (37) e ( n ) (n ).
We satisfy the continuity constraint of phase by spreading the error through the whole missing region. The interpolated phases are then computed during the missing region as follows: n = n + n n1 n2 n1 (38)
where n [n1 + 1,n2] and is chosen to ensure the con tinuity constraint at the end boundary: (n2) (n2).0 Since (n2) is a known 2 modulus, the number of solutions for is infinite. The smallest one is retained, =

e + 2, e e 2, e e, otherwise.
Given the predicted amplitudes and frequencies of the partials from both sides of the missing region, we are able to interpolate reliably the missing region of a partial. The capability of this new interpolation scheme will be evaluated in the remainder of this section, where a synthetized version of the interpolated sinusoidal representation is
Fig. 7. Simulating a missing region.
They were asked to note these four versions using the 100-point Mushra scale. The marks obtained by the two interpolated versions are plotted in Fig. 9. As can be seen in Fig. 9(a), a high-quality interpolation of monophonic signals having vibrato (saxophone tone) or tremolo (vibraphone tone) for missing region sizes close to 1 s is achieved. The audio signals having more complex modulations, such as the singing voice or the orchestra piece, are harder to interpolate, but the LP-based method is a significant improvement [Fig. 9(b)]. 5 EXTRAPOLATION OF UNMATCHED PARTIALS Considering that the matching between partials of B and A is done correctly, unmatched partials of B belong to a note decaying in the missing region and unmatched sinu-
soids of A belong to a note that started in the missing region (see Fig. 1). Let lB be the maximum length of the extrapolation of unmatched partials of B and lA the maximum length of the extrapolation of unmatched partials of A. The extrapolation of unmatched partials Pi or Pj is done according to the predicted parameters Pi or Pj. The predicted frequencies i f are used as is, and the extrapolated phases are comor fj puted using Eq. (36). In general, the amplitudes of the partials have a predictable behavior during the ending of the note (sustain or decay). The predicted amplitude of the unmatched partial of B can then be used safely to detect at which frame the partial should end. The extrapolated amplitude Ai(n) is (n) faded as follows: then the predicted amplitude Ai Ai(n) Ai(n) i(n) (41)
Fig. 8. Objective comparison of LP-based interpolation ( ) and polynomial interpolation ( ) on three sound signals. (a) Saxophone tone with vibrato. (b) Vibraphone. (c) Harpsichord.
900 J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October
with n n1 max Ain1 + lB, 0. in = lB (42)
If the extrapolated amplitude Ai(n) becomes negative at a frame n < n1 + lB, the extrapolated partial Pi ends at frame n 1. As a consequence the partial may end before n1 + lB, as shown in Fig. 1. On the other hand, the amplitude of the partials during an abrupt onset cannot be deduced from the amplitude of the partial during the sustain part of the note. To at least simulate an onset, all unmatched partials of A should begin at the same frame index n2 lA. The extrapolated amplitude Aj(n) is then the predicted amplitude Aj(n) faded as follows: Aj(n) Aj(n) j(n) with n2 n Ajn2 lA. jn = lA (44) (43)

The parameters lB and lA should be chosen according to the targeted application. For interpolating the sinusoidal data lost due to a transmission error, the maximum gap size allowed is generally small due to the limited data buffering capability of the decoder. In this approach the extrapolation should be parameterized to be tolerant to mismatch that occurred during the matching step of the algorithm. This can be done by setting lB lA n2 nto ensure a fade in or out of unmatched sinusoids. During the sinusoidal analysis step or with a digital data restoration application, however, some extra information about the spectral content can be used to estimate the frame index where the unmatched partials of A should start. The parameter lA can be set to an onset index estimate for every gap occurring. 6 INTERPOLATION OF MISSING AUDIO DATA This section compares three methods of prediction for missing audio data using subjective listening tests with the same protocol as the one used in Section 4.5. The temporal method uses 2000 temporal samples from both sides of the
The extrapolated partial Pj starts at the smaller frame index (n + k) > 0, for all k 0. n n2 lA so that Aj
Fig. 9. Results of listening tests comparing polynomial-based method () and LP-based method () for five gap sizes. Symbols means of votes; linesconfidence intervals for each method.
J. Audio Eng. Soc., Vol. 53, No. 10, 2005 October 901
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.

[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.