Musical Sound Modeling With Sinusoids Plus Noise

(DRAFT, August 1994)

Xavier Serra Fundació Phonos - Institut Universitari de l’Audiovisual Rambla 31 08002 Barcelona, Spain

1 Introduction

When generating musical sound on a digital computer, it is important to have a good model whose parameters provide a rich source of meaningful sound transformations. Three basic model types are in prevalent use today for musical sound generation: instrument models, spectrum models, and abstract models. Instrument models attempt to parametrize a sound at its source, such as a violin, clarinet, or vocal tract. Spectrum models attempt to parametrize a sound at the basilar membrane of the ear, discarding whatever information the ear seems to discard in the spectrum. Abstract models, such as FM, attempt to provide musically useful parameters in an abstract formula.

This chapter addresses the second category of synthesis techniques: spectrum modeling. The main advantage of this group of techniques is the existence of analysis procedures that extract the synthesis parameters out of real sounds, thus being able to reproduce and modify actual sounds. Our particular approach is based on modeling sounds as stable sinusoids (partials) plus noise (residual component), therefore analyzing sounds with this model and generating new sounds from the analyzed data. The analysis procedure detects partials by studying the time-varying spectral characteristics of a sound and represents them with time-varying sinusoids. These partials are then subtracted from the original sound and the remaining “residual” is represented as a time-varying filtered white noise component. The synthesis procedure is a combination of additive synthesis for the sinusoidal part, and subtractive synthesis for the noise part.

This analysis/synthesis strategy can be used for either generating sounds (synthesis) or transforming pre-existing ones (sound processing). To synthesize sounds we generally want to model an entire timbre family, i.e., an instrument, and that can be done by analyzing single tones and isolated note transitions performed on an instrument, and building a data base that characterizes the whole instrument or any desired timbre family, from which new sounds are synthesized. In the case of the sound processing application the goal is to manipulate any given sound, that is, not being restricted to isolated tones and not requiring a previously built data-base of analyzed data.

Some of the intermediate results from this analysis/synthesis scheme, and some of the techniques developed for it, can also be applied to other music related problems, e.g., sound compression, sound source separation, musical acoustics, music perception, performance analysis,... but a discussion of these topics is beyond the current presentation.

2 Background

Additive synthesis is the original spectrum modeling technique. It is rooted in Fourier’s theorem which states that any periodic waveform can be modeled as a sum of sinusoids at various amplitudes and harmonic frequencies. Additive synthesis was among the first synthesis techniques in computer music. In fact, it was described extensively in the very first article of the very first issue of the Computer Music Journal (Moorer, 1977).

In the early 1970s, Andy Moorer developed a series of analysis programs to support additive synthesis. He first used the “heterodyne filter” to measure the instantaneous amplitude and frequency of individual sinusoids (Moorer, 1973). The heterodyne filter implements a single frequency bin of the Discrete Fourier Transform (DFT), using the rectangular window. The magnitude and phase derivative of the complex numbers produced by the sliding DFT bin provided instantaneous amplitude and frequency estimates. The next implementation (Moorer, 1978) was based on the Digital Phase Vocoder (Portnoff, 1976). In this system, the fast Fourier transform (FFT) was used to provide, effectively, a heterodyne filter at each harmonic of the fundamental frequency. The use of a nonrectangular window gave better isolation among the spectral components.

The main problem with the phase vocoder was that inharmonic sounds, or sounds with timevarying frequency characteristics, were difficult to analyze. The FFT can be regarded as a fixed filter bank or “graphic equalizer”: If the length of the FFT is N, then there are N narrow bandpass filters, slightly overlapping, equally spaced between 0Hz and the sampling rate. In the phase vocoder, the instantaneous amplitude and frequency are computed only for each “channel filter” or “bin.” A consequence of using a fixed-frequency filter bank is that the frequency of each sinusoid is not normally allowed to vary outside the bandwidth of its channel, unless one is willing to combine channels in some fashion which requires extra work. (The channel bandwidth is nominally the sampling rate divided by the FFT length.) Also, the analysis system was really set up for harmonic signals—you could analyze a piano if you had to, but the progressive sharpening of the partials meant that there would be frequencies where a sinusoid would be in the crack between two adjacent FFT bins. This was not an insurmountable condition (the adjacent bins could be combined intelligently to provide accurate amplitude and frequency envelopes), but it was inconvenient and outside the original scope of the analysis framework of the phase vocoder.

The PARSHL program was developed at CCRMA for the purpose of supporting inharmonic and pitch-changing sounds (Smith and Serra, 1987). PARSHL was a simple application of FFT peaktracking technology commonly used in the Navy signal processing community (General Electric, 1977; Wolcin 1980a; 1980b; Smith and Friedlander, 1984). As in the phase vocoder, a series of FFT frames is computed by PARSHL. However, instead of writing out the magnitude and phase derivative of each bin, the FFT is searched for peaks, and the largest peaks are “tracked” from frame to frame. The principal difference in the analysis is the replacement of the phase derivative in each FFT bin by interpolated magnitude peaks across FFT bins. This approach is better suited for analysis of inharmonic sounds and pseudo-harmonic sounds with important frequency variation in time.

Independently at about the same time, Quatieri and McAulay developed a technique similar to PARSHL for analyzing speech (McAulay and Quatieri, 1984; 1986). Both systems were built on top of the short-time Fourier Transform (Allen, 1977).

The PARSHL program worked well for most sounds created by simple physical vibrations or driven periodic oscillations. It went beyond the phase vocoder to support spectral modeling of inharmonic sounds. A problem with PARSHL, however, is that it was unwieldy to represent noiselike signals such as the attack of many instrumental sounds. Using sinusoids to simulate noise is extremely expensive because, in principle, noise consists of sinusoids at every frequency within the band limits. Also, modeling noise with sinusoids does not yield a flexible sound representation useful for music applications. Therefore the next natural step to take in spectral modeling of musical sounds was to represent sinusoids and noise as two separate components (Serra, 1989; Serra and Smith, 1990).

3 The Deterministic plus Stochastic Model

A sound model assumes certain characteristics of the sound waveform or the sound generation mechanism. In general, every analysis/synthesis system has an underlying model. The sounds produced by musical instruments, or by any physical system, can be modeled as the sum of a set of sinusoids plus a noise residual. The sinusoidal, or deterministic, component normally corresponds to the main modes of vibration of the system. The residual comprises the energy produced by the excitation mechanism which is not transformed by the system into stationary vibrations plus any other energy component that is not sinusoidal in nature. For example, in the the sound of winddriven instruments, the deterministic signal is the result of the self-sustained oscillations produced inside the bore and the residual is a noise signal that is generated by the turbulent streaming that takes place when the air from the player passes through the narrow slit. In the case of bowed strings the stable sinusoids are the result of the main modes of vibration of the strings and the noise is generated by the sliding of the bow against the string, plus by other non-linear behavior of the bow-string-resonator system. This type of separation can also be applied to vocal sounds, percussion instruments and even to non-musical sounds produced in nature.

A deterministic signal is traditionally defined as anything that is not noise (i.e., an analytic signal, or perfectly predictable part, predictable from measurements over any continuous interval). However, in the present discussion the class of deterministic signals considered is restricted to sums of quasi-sinusoidal components (sinusoids with slowly varying amplitude and frequency). Each sinusoid models a narrowband component of the original sound and is described by an amplitude and a frequency function.

A stochastic, or noise, signal is fully described by its its power spectral density which gives the expected signal power versus frequency. When a signal is assumed stochastic, it is not necessary to preserve either the instantaneous phase or the exact magnitude details of individual FFT frames.

Therefore, the input sound s(t) is modeled by,

s(t) = ∑_{r=1}^{R} A_{r}(t) cos[θ_{r}(t)] + e(t)  (1)

where A_{r}(t) and θ_{r}(t) are the instantaneous amplitude and phase of the r^{th} sinusoid, respectively, and e(t) is the noise component at time t (in seconds).

The model assumes that the sinusoids are stable partials of the sound and that each one has a slowly changing amplitude and frequency. The instantaneous phase is then taken to be the integral of the instantaneous frequency ω_{r}(t), and therefore satisfies

θ_{r}(t) = ∫_{0}^{t} ω_{r}(τ)dτ  (2)

where ω(t) is the frequency in radians, and r is the sinusoid number.

By assuming that e(t) is a stochastic signal, it can be described as filtered white noise,

e(t) = ∫_{0}^{t} h(t, τ)u(τ)dτ  (3)

where u(t) is white noise and h(t, τ) is the response of a time varying filter to an impulse at time t. That is, the residual is modeled by the convolution of white noise with a time-varying frequency-shaping filter.

4 General Diagram of the Analysis/Synthesis Process

The deterministic plus stochastic model has many possible implementations and we will present a general one while giving indications on variations that have been proposed. Both the analysis and synthesis are frame based processes with the computation done one frame at a time. Throughout this description we will consider that we have already processed a few frames of the sound and we are ready to compute the next one.

Fig. 1 shows the block diagram for the analysis. First we prepare the next section of the sound to be analyzed by multiplying it with an appropiate analysis window. Its spectrum is obtained by the Fast Fourier Transform (FFT) and the prominent peaks in the magnitude spectrum are detected and incorporated into the existing partial trajectories by means of a peak continuation algorithm. The relevance of this algorithm is that it detects the magnitude, frequency and phase of the partials present in the original sound (the deterministic component). It requires a certain memory of past frames in order to better track the partials in time by making some decisions about the time evolution of the partials. When the sound is pseudo-harmonic a fundamental frequency detection step can improve the analysis by using the fundamental frequency information in the peak continuation algorithm and in the generation of the analysis window for the next frame (pitch-synchronous analysis).

The stochastic component of the current frame is calculated by first generating the deterministic signal with additive synthesis, and then subtracting it from the original waveform in the time domain. This is possible because the phases of the original sound are matched and therefore the shape of the time-domain waveform preserved. The stochastic representation is then obtained by performing a spectral fitting of the residual signal.

In Fig. 2 the block diagram of the synthesis process is shown. The deterministic signal, i.e., the sinusoidal component, results from the magnitude and frequency trajectories, or their transformation, by generating a sine wave for each trajectory (i.e., additive synthesis). This can either be implemented in the time domain with the traditional oscillator bank method or in the frequency domain using the inverse-FFT approach.

The synthesized stochastic signal is the result of generating a noise signal with the time-varying spectral shape obtained in the analysis (i.e., subtractive synthesis). As with the deterministic synthesis, it can be implemented in the time domain by a convolution or in the frequency domain by creating a complex spectrum (i.e., magnitude and phase spectra) for every spectral envelope of the residual and performing an inverse-FFT.

DRAFT, August 1994 7

Figure 1: Block diagram of the analysis process. Figure 2: Block diagram of the synthesis process.

5 Magnitude and Phase Spectra

The computation of the magnitude and phase spectra of the current frame is the first step in the analysis. It is in these spectra that the sinusoids are tracked and the decision takes place as to whether a part of the signal is considered deterministic or noise. The computation of the spectra is carried out by the STFT technique (Allen, 1977; Serra, 1989), discussed elsewhere in this book ([add M. H. Serra reference]).

The control parameters for the STFT (window-length, window-type, FFT-size, and hop-size) have to be set in accordance with the sound to be processed. First of all, a good resolution of the spectrum is needed since the process that tracks the partials has to be able to identify the peaks which correspond to the deterministic component. The time-frequency compromise of the STFT has to be well understood when setting the window-length, window-type, and FFT-size. Also the phase information is particularly important for subtracting the deterministic component to find the residual; we should use an odd-length analysis window and the windowed data should be centered in the FFT-buffer at the origin in order to obtain the phase spectrum without a linear phase trend induced by the window (”zero-phase” windowwing).

Another point concerning the STFT is related to the synthesis part of the system. The synthesis process described here is completely independent from the analysis and does not use the STFT. This implies that the restrictions imposed by the STFT when the inverse transform is also performed, for example, that the analysis windows add to a constant, are unnecessary here. The STFT parameters are more flexible, and we can vary them during the course of one sound, if that is required to improve the detection of partials.

6 Peak Detection

Once the spectrum of the current frame is computed, the next step is to detect its prominent magnitude peaks. Theoretically, a sinusoid that is stable both in amplitude and in frequency (a partial) has a well defined frequency representation: the transform of the analysis window used to compute the Fourier transform. It should be possible to take advantage of this characteristic to DRAFT, August 1994 9

Figure 3: Sound selection and windowing. a. Portion of a violin sound to be used in the analysis of the current frame. b. Hamming window. c. Windowed sound.

DRAFT, August 1994 10

Figure 4: Computing the FFT. a. packing of the sound into the FFT buffer for a zero phase spectrum. b. Magnitude spectrum. c. Phase spectrum.

distinguish partials from other frequency components. However, in practice this is rarely the case since most natural sounds are not perfectly periodic and do not have nicely spaced and clearly defined peaks in the frequency domain. There are interactions between the different components, and the shapes of the spectral peaks cannot be detected without tolerating some mismatch. Only some instrumental sounds (e.g., the steady-state part of an oboe sound) are periodic enough and sufficiently free from prominent noise components that the frequency representation of a stable sinusoid can be recognized easily in the peak detection.

A practical solution is to detect as many peaks as possible and delay the decision of what is a deterministic, or “well behaved” partial to the next step in the analysis: the peak continuation algorithm.

A “peak” is defined as a local maximum in the magnitude spectrum, and the only practical restrictions to be made in the peak search are to have a frequency range and a magnitude threshold. The frequency range can be set to auditory based limits, e.g. 20 − 16000 Hz, and the magnitude threshold can be set to the minimum audible energy, e.g., 4 dB SPL.

Due to the sampled nature of the spectra returned by the FFT, each peak is accurate only to within half a sample. A spectral sample represents a frequency interval of f_{s}/N Hz, where f_{s} is the sampling rate and N is the FFT size. Zero-padding in the time domain increases the number of spectral samples per Hz and thus increases the accuracy of the simple peak detection. However, to obtain frequency accuracy on the level of 0.1% of the distance from the top of an ideal peak to its first zero crossing (in the case of a Rectangular window), the zero-padding factor required is 1000. A more efficient spectral interpolation scheme is to zero-pad only enough so that quadratic (or other simple) spectral interpolation, using only samples immediately surrounding the maximum-magnitude sample, suffices to refine the estimate to 0.1% accuracy.

Even though we cannot rely on the exact shape of the peak to decide whether it is a partial or not, it is sometimes useful to have a measure of how close its shape is to the ideal sinusoidal peak. This measure can be obtained by calculating the difference from the samples of the measured peak to the samples of the analysis window transform centered at the measured frequency and scaled to the measured magnitude. This information, plus the frequency, magnitude, and phase of the peak, DRAFT, August 1994 12

Figure 5: Peak detection. a. Peaks in the magnitude spectrum. b. Peaks in the phase spectrum.

can help in the peak continuation process.

7 Fundamental Frequency Detection

Before continuing a set of peak trajectories through the current frame it is useful to search for a possible fundamental frequency, that is, for periodicity. If it exists, we will have more information to work with, and it will simplify and improve the tracking of partials. This fundamental frequency can also be used to set the length of the analysis window used in the following frame, in order to maintain constant the number of periods to be analyzed at each frame and to get the best time-frequency trade-off possible. This is called a pitch-synchronous STFT analysis.

Given a set of spectral peaks, with magnitude and frequency values for each one, there are many possible fundamental detection strategies (Piszczalski and Galler, 1979; Terhardt, 1982; Hess, 1983; Doval and Rodet, 1993; Maher and Beauchamp, 1994). For this presentation we restrict ourselves to single-source sounds and assume that a fundamental peak or one of its first few partials exists. With these two constraints, plus the fact that there is some number of buffered frames, the algorithm can be quite simple.

The fundamental frequency can be defined as the common divisor of the harmonic series, that best explains the spectral peaks found in the current frame. The first step is to find the possible candidates inside a given range, preferring the lowest succesful candidate. This can be done by stepping through the range by small increments, or by only considering as candidates the frequencies of the spectral peaks and frequencies related to them by simple integer ratios (e.g., 1/2, 1/3, 1/4) that lie inside the range. This last approach simplifies our search enormously.

Once the possible candidates have been chosen we need a way to measure the “goodness” of the resulting harmonic series compared with the actual spectral peaks. A suitable error measure is given by Maher and Beauchamp (Maher and Beauchamp, 1994) which is based on the weighted differences between the measured peaks and the ideal harmonic series. The harmonic series with the smallest error is chosen, and since there is a certain frame memory, the resulting fundamental frequency can be compared with that of previous and future frames before accepting its value.

Once accepted, this information is passed to the peak continuation algorithm and is used to set the window-length in the analysis of the next frame.

8 Peak Continuation

Once the spectral peaks of the current frame have been detected, the peak continuation algorithm adds them to the incoming peak trajectories. The schemes used in PARSHL (Smith and Serra, 1987) and in the sinusoidal model (McAulay and Quatieri, 1984; 1986) find peak trajectories both in the noise and deterministic parts of a waveform, thus obtaining a sinusoidal representation for the whole sound. These schemes are unsuitable when we want the trajectories to follow just the partials. For example, when the partials change in frequency substantially from one frame to the next, these algorithms easily switch from the partial that they were tracking to another one which at that point is closer.

The algorithm described here is intended to track partials in a variety of sounds, even though the behavior of a partial, and therefore the way to track it, varies depending on the sound. Whether we have speech, a harmonic instrumental tone, a gong sound, a sound of an animal, or any other, the time evolution of the component partials will vary. Thus, the algorithm requires some knowledge about the characteristics of the sound that is being analyzed. In the current algorithm there is no attempt to make the process completely automatic and some of the characteristics of the sound are specified through a set of parameters.

The basic idea of the algorithm is that a set of guides advances in time through the spectral peaks, looking for the appropriate ones (according to the specified constraints) and forming trajectories out of them. Thus, a guide is an abstract entity which is used by the algorithm to create the trajectories and the trajectories are the actual result of the peak continuation process. The instantaneous state of the guides, their frequency, magnitude and status, is kept in g_{1}, g_{2}, g_{3}, …, g_{p}, where p is the number of existing guides. These values are continuously updated as the guides are turned on, advanced, and finally turned off.

Before explaining the algorithm, let us describe some of the control parameters used in it.

• maximum-peak-frequency-deviation: The maximum frequency change allowed in a guide from one frame to the next. It is useful to make this parameter a function of relative frequency, i.e., such that the admissible change is bigger for higher frequencies than for lower ones. Thus, the deviation can be defined on a log-scale which is perceptually more meaningful than a linear frequency scale.

• maximum-peak-magnitude-deviation. The maximum relative magnitude change allowed by a guide from one frame to the next.

• peak-contribution-to-guide. The frequency of each guide does not have to correspond to the frequency of the actual partial. It is updated every time it incorporates a new peak. This parameter controls how much the guide frequency changes when a new peak is incorporated and it is useful to confine a guide to a narrow frequency band and to smooth the frequency track.

• fundamental-contribution-to-guide. When a fundamental frequency has been found in the current frame it can be used to refine the value of all the guides. This parameter expresses the contribution of the fundamental frequency to each guide frequency. When it is maximum and peak-contribution-to-guide is 0 the guides are restricted to the harmonic series of the fundamental frequency.

• maximum-sleeping-time. When a guide has not found a continuation peak for a certain number of frames, the guide is killed. This parameter specifies the maximum “inactive” time which translated to the maximum number of frames that the guide can be alive while not finding continuation peaks. If its value is 0 the guide is killed as soon as it does not find a continuation peak. This value is automatically set to the whole duration of the sound in the case of pseudo-harmonic sounds so that guides are not killed during the analysis and they correspond to the harmonic series throughout the duration.

• maximum-length-of-filled-gaps. Given that a certain sleeping time is permitted, we may wish to fill the resulting gaps. This parameter specifies the length of the biggest gap to be filled (a number smaller or equal than maximum-sleeping-time). The gaps are filled by interpolating between the end points in the trajectory.

• minimum-trajectory-length. Once the trajectories have been created and continued through a certain number of frames, it might be useful to delete the short ones in order to maintain only the ones that have tracked stable partials. This parameter sets the minimum length of the trajectories to be accepted.

To describe the peak continuation algorithm let us assume that the guide values at the current frame are g_{1}, g_{2}, g_{3}, …, g_{n}, where n is the number of existing guides. We want to continue the n guides through the peaks of the current frame with values p_{1}, p_{2}, p_{3}, …, p_{m}, thus continuing the corresponding trajectories.

There are three steps in the algorithm: (1) guide advancement, (2) update of the guide values, and (3) start of new guides. Next, these steps are described.

When a fundamental has been found in the current frame the guides are updated according to fundamental-contribution-to-guide, otherwise they keep the values that resulted from the previous frame. Then each guide selects the peak closest in frequency and magnitude to its current value. The rth guide claims peak p_{i} for which |g_{r}.f − p_{i}.f| and |g_{r}.m − p_{i}.m| are a minimum, where g_{r}.f and g_{r}.m are the frequency and magnitude of the guide r and p_{i}.f and p_{i}.m are the frequency and magnitude of the peak r. The relative change in frequency must be less than maximumpeak-frequency-deviation and the relative change in magnitude less than maximum-peak-magnitude-deviation. More weight is given to the frequency difference than to the magnitude difference since the definition of a partial is mainly given by frequency continuity conditions and not by magnitude considerations. Other criteria could be added in this selection process, in particular the phase of the peak and a measure of the “goodness” of the peak. We can select the peak with the phase that makes the most sense given the phase value of the previous peak and also the one that is closest to an ideal sinusoid. These criteria can be used when the other ones do not give us a clear continuation peak. The possible situations are as follows:

1. If a match is found within the maximum frequency and magnitude deviations, the guide is continued (unless there is a conflict to resolve, as described below). The selected peak is incorporated into the corresponding trajectory.

Figure 6: Peak Continuation process.

2. If no match is found, it is assumed that the corresponding trajectory must “turn off” entering the current frame, and its current frequency is matched to itself with zero magnitude. Since the trajectory amplitudes are linearly ramped from one frame to the next, the terminating trajectory ramps to zero over the duration of one hop size. Whether the actual guide is “killed” or not, depends on the maximum-sleeping-time.

3. If a guide claims a peak which has already been claimed by another guide, we give the peak to the guide that is closest in frequency and in magnitude, and the “loser” looks for the next best match. If the current guide loses the conflict, it simply picks the best available non-conflicting peak which is within the maximum-peak-frequency-deviation and maximum-peak-magnitude-deviation. If the current guide wins the conflict, it calls the assignment procedure recursively on behalf of the dislodged guide. When the dislodged guide finds the same peak and wants to claim it, it sees there is a conflict which it loses and moves on. This process is repeated for each guide, solving conflicts recursively, until all possible matches are made.

Once all the existing guides and their trajectories have been continued through the current frame, the guide values are updated. There are two possible situations:

1. If a guide finds a continuation peak, its frequency is updated according to:

ĝ_{r}.f = α(p_{i}.f − g_{r}.f) + g_{r}.f,  αϵ[0, 1]  (4)

where p_{i}.f is the frequency of the peak that the guide has found, and α is the peak-contribution-to-guide. The magnitude of the guide is updated to the magnitude of the peak taken.

2. If a guide has not found a continuation peak for the past maximum-sleeping-time frames, the guide is killed. If it is still within the sleeping-time it keeps the same value. In order to distinguish between the different states of the guides, each one keeps a state variable indicating whether it is alive, sleeping, or killed.

New guides, and therefore new trajectories, are created from the peaks of the current frame that are not incorporated into trajectories by the existing guides. If there are killed or unused guides, a new guide can be started. A guide is created by searching through the “unclaimed” peaks of the frame for the one with the highest magnitude. The frequency value of the selected peak is the frequency of the new guide. The actual trajectory is started in the previous frame where its amplitude value is set to 0 and its frequency value to the same as the current frequency, thus ramping linearly in amplitude to the current frame. In pseudo-harmonic sounds we do not want to start or kill guides in the middle of the sound, even though the trajectories may not be active throughout the whole duration.

The attack portion of most sounds is quite “noisy,” and the search for partials is harder in such a rich spectrum. A useful modification to the analysis is to perform it backwards in time. The tracking process encounters the end of the sound first, and since this is a very stable part in most instrumental sounds, the algorithm finds a very clear definition of the partials. When the guides arrive at the attack, they are already tracking the main partials and can reject non-relevant peaks appropriately, or at least evaluate them with some acquired knowledge.

The peak continuation algorithm presented is only one approach to the peak continuation problem. The creation of trajectories from the spectral peaks is compatible with very different strategies and algorithms; for example, hidden Markov models have been applied (Garcia, 1992; Depalle, Garcia and Rodet, 1993). An N Markov model provides a probability distribution for a parameter in the current frame as a function of its value across the past N frames. With a hidden Markov model we are able to optimize groups of trajectories according to a defined criteria, such as frequency continuity. This type of approach might be very valuable for tracking partials in polyphonic sounds and complex inharmonic tones. In particular, the notion of “momentum” is introduced, helping to properly resolve crossing fundamental frequencies.

9 Deterministic Subtraction

The output of the peak continuation algorithm is a set of peak trajectories updated for the current frame. From these trajectories a series of sinusoids can be synthesized which reproduce the instantaneous phase and amplitude of the partials of the original sound. Thus, it is possible to subtract the synthesized sinusoids from the original sound and obtain a residual which is substantially free of the deterministic part.

One frame of the deterministic part of the sound, s(m), is generated by

d(m) = ∑_{r=1}^{R} Â_{r} cos[mω̂_{r} + φ̂_{r}], m = 0, 1, 2, …, S − 1 (5)

where R is the number of trajectories present in the current frame and S is the length of the frame. To avoid “clicks” at the frame boundaries, the parameters (Â_{r}, ω̂_{r}, φ̂_{r}) are smoothly interpolated from frame to frame.

Let (Â_{r}^{(l−1)}, ω̂_{r}^{(l−1)}, φ̂_{r}^{(l−1)}) and (Â_{r}^{l}, ω̂_{r}^{l}, φ̂_{r}^{l}) denote the sets of parameters at frames l − 1 and l for the rth frequency trajectory (we will simplify the notation by omitting the subscript r). These parameters are taken to represent the state of the signal at time 0 (the left endpoint) of the frame.

The instantaneous amplitude Â(m) is easily obtained by linear interpolation,

Â(m) = Â^{l−1} + \frac{(Â^{l} − Â^{l−1})}{S}m  (6)

where m = 0, 1, …, S − 1 is the time sample into the lth frame.

Frequency and phase values are tied together (frequency is the phase derivative), and both control the instantaneous phase θ̂(m), defined as

θ̂(m) = mω̂ + φ̂  (7)

Given that four variables affect the instantaneous phase: ω̂^{(l−1)}, φ̂^{(l−1)}, ω̂^{l}, and φ̂^{l}, we need three degrees of freedom for its control, but linear interpolation gives only one. Therefore, we need a cubic polynomial as an interpolation function,

θ̂(m) = ζ + κm + ηm^{2} + ιm^{3}  (8)

It is unnecessary to go into the details of solving this equation since they are described by McAulay and Quatieri (McAulay and Quatieri, 1986). The result is

θ̂(m) = φ̂^{(l−1)} + ω̂^{(l−1)}m + ηm^{2} + ιm^{3}  (9)

where η and ι are calculated using the end conditions at the frame boundaries,

η = \frac{3}{S^{2}} (φ̂^{l} − φ̂^{l−1} − ω̂^{l−1}S + 2πM) − \frac{1}{S} (ω̂^{l} − ω̂^{l−1})
ι = −\frac{2}{S^{3}} (φ̂^{l} − φ̂^{l−1} − ω̂^{l−1}S + 2πM) + \frac{1}{S^{2}} (ω̂^{l} − ω̂^{l−1})  (10)

This gives a set of interpolating functions depending on the value of M, among which we select the maximally smooth function. This is done by choosing M to be the integer closest to x, where x is

x = \frac{1}{2π} [(φ̂^{l−1} + ω̂^{l−1}S − φ̂^{l}) + \frac{S}{2} (ω̂^{l} − ω̂^{l−1})]  (11)

Finally, the synthesis equation for frame l becomes

d^{l}(m) = ∑_{r=1}^{R^{l}} Â_{r}^{l}(m) cos[θ̂_{r}^{l}(m)]  (12)

which goes smoothly from the previous to the current frame with each sinusoid accounting for both the rapid phase changes (frequency) and the slowly varying phase changes.

The synthesized deterministic component can be subtracted from the original sound in the time domain by

e(n) = w(n) ∗ (s(n) − d(n)) n = 0, 1, …, N  (13)

where e(n) is the residual, w(n) a smoothing window, s(n) the original sound, d(n) the deterministic component, and N the length of the window. It is desirable to set N smaller than the windowlength used in the deterministic analysis in order to improve the time resolution of the residual signal. While in the deterministic analysis the window-length was chosen large enough to obtain a good partial separation in the frequency domain, in the deterministic subtraction we are looking for especially good time resolution. This is particularly important in the attacks of percussion instruments, and we will discuss this issue from that perspective even though it is also needed wherever there are fast dynamic changes.

It is easy to observe that the attacks of most percussion instruments are mainly composed of what we call “stochastic components”. The attack has a very broad band of energy and it is only after this attack that the main modes of vibration of the system (the deterministic component) become important. In the stochastic component we do not require the frequency resolution that is needed in the deterministic analysis, but instead we want a better time resolution. The deterministic component cannot maintain the sharpness of the attacks, because even though we can have a very small hop-size in the analysis, we also need a long window-length, and this length will determine the time resolution. Since it is this deterministic signal that we subtract from the original sound, the resulting residual signal has the sharp attacks also smeared. In order to improve the stochastic analysis we can “fix” this residual so that the sharpness of the attacks of the original sound are preserved. This is accomplished by defining a small window (a good choice is two times the hop-size) in the deterministic subtraction and comparing the resulting residual with the original waveform. Whenever the residual has a greater energy than the original waveform, its amplitude is re-scaled to the magnitude of the original. Then the stochastic analysis is performed on this scaled residual. Thus, the smaller we choose the hop-size (which is the same for all steps in the analysis) the better time resolution we will get in the residual.

Tests on this residual can be performed to check whether the deterministic plus stochastic decomposition has been successful (Serra, 1994a). Ideally the resulting residual should be as close as possible to a stochastic signal. Since the autocorrelation function of white noise is an impulse, a measure of correlation relative to total power could be a good measure of how close we are to white noise,

c = \frac{∑_{l=1}^{L−1} |r(l)|}{(L − 1)r(0)}  (14)

where r(l) is the autocorrelation estimate for L lags of the residual, and c will be close to 0 when the signal is stochastic. A problem with this measure is that it does not behave well when partials DRAFT, August 1994 22

Figure 7: Deterministic subtraction. a. original sound. b. deterministic synthesis. c. residual sound.

are still left in the signal; for example, it does not always decrease as we progressively subtract partials from a sound. A simpler and sometimes better indication of the quality of the residual is to measure the energy of the residual as a percentage of the total sound energy. Even though a problem with this measure is that it cannot distinguish subtracting partials from subtracting noise, and its value will always decrease as long as we subtract some energy, it is still a good measure for choosing the best analysis parameters.

This sound decomposition is useful in itself for a number of applications. The deterministic component is a set of partials, and the residual includes noise and very unstable components of the sound. This technique has been used by Schumacher and Chafe (Chafe, 1990; Schumacher and Chafe, 1990) to study bow noise in string instruments and breath noise in wind instruments. In general, this decomposition strategy can give a lot of insight into the makeup of sounds.

The residual component is the part of the instrumental sounds that the existing synthesis techniques have a harder time reproducing, and it is especially important during the attack. A practical application would be to add these residuals to synthesized sounds in order to make them more realistic. Since these residuals remain largely invariant throughout most of the instrumental range, only a few residuals would be necessary to cover all the sounds of a single instrument.

10 Computation of the Stochastic Part

One of the underlying assumptions of the current model is that the residual is a stochastic signal. Such an assumption implies that the residual is fully described by its amplitude and its general frequency characteristics. It is unnecessary to keep either the instantaneous phase or the exact spectral shape information. Based on this, a frame of the stochastic residual can be completely characterized by a filter, i.e., this filter encodes the amplitude and general frequency characteristics of the residual. The representation of the residual for the overall sound will be a sequence of these filters, i.e., a time-varying filter. Initially, there is a set of filter coefficients for each frame.

The filter design problem is generally solved by performing some sort of curve fitting in the magnitude spectrum of the current frame (Strawn, 1980; Sedgewick, 1988). Standard techniques are: spline interpolation (Cox, 1971), the method of least squares (Sedgewick, 1988), or straight line approximations (Phillips, 1968). For our purpose a simple line-segment approximation to the log-magnitude spectrum is accurate enough and gives the desired flexibility.

One way to carry out the line-segment approximation is to step through the magnitude spectrum and find local maxima in each of several defined sections, thus giving equally spaced points in the spectrum that are connected by straight lines to create the spectral envelope. The accuracy of the fit is given by the number of points, and that can be set depending on the sound complexity. Other options are to have unequally spaced points, for example, logarithmically spaced, or spaced according to perceptual criteria.

Another practical alternative is to use a type of least squares approximation called linear predictive coding, LPC (Makhoul, 1975; Markel and Gray, 1976). LPC is a popular technique used in speech research for fitting an nth-order polynomial to a magnitude spectrum. For our purposes, the line-segment approach is more flexible than LPC, and even though LPC results in less analysis points, the flexibility is considered more important.

11 Representation of the Analysis Data

The output of the deterministic analysis is a set of amplitude and frequency functions with one breakpoint for every frame. From these functions a series of sinusoids can be synthesized which reproduce the deterministic part of the sound. The phase trajectories are not kept because they are unnecessary in the final synthesis, they are perceptually irrelevant in most cases, and they make it harder to perform modifications.

It is important to mention that we have found some situations in which the preservation of phases has made a difference in the quality of the resynthesis. These are: badly analyzed sounds, very low instrumental tones, and some vocal sounds. In the case of badly analyzed sounds, some of the trajectories may actually be tracking non-deterministic parts of the signal, in which case the phase of the corresponding peaks is important to recover the noisy characteristics of the signal. In the case when the analyzed sound has a very low fundamental, maybe lower than 30 Hz, and the DRAFT, August 1994 25

Figure 8: Stochastic approximation from the sound in Fig. 7. a. Original spectrum. b. Residual spectrum and its line-segment approximation.

partials are phase-locked to the fundamental, the period is perceived as a pulse and the phase of the partials is required to maintain this perceptual effect. Also in the case of some vocal sounds, the higher partials have a high degree of modulation that cannot be completely recovered from the frequency and magnitude information of the partials, but that seems to be maintained when we add the phases of the peaks.

The resulting amplitude and frequency functions can be further processed to achieve a data reduction of the representation or to smooth the functions. A data reduction strategy is to perform a line-segment approximation on each function, thus reducing the number of breakpoints (Grey, 1975; Strawn, 1980). However, for the purpose of easy manipulation of the representation it is useful to have equally spaced points along each function, and thus it may be better to keep one breakpoint per frame as returned by the analysis, unless data reduction is a priority. Another alternative for data reduction is to combine groups of similar functions into a single one, thus reducing the number of functions.

The stochastic analysis returns an envelope for every frame. These envelopes can either be interpreted as a series of envelopes or frequency-shaping filters, one per frame, or as time-varying, equally spaced band-pass filters, each one centered at each breakpoint. It is convenient to normalize these envelopes by dividing them by their average magnitude so that we can control the spectral shape of the noise independently of its time-varying magnitude.

12 Modifications of the Analysis Data

The representation resulting from the analysis is very suitable for modification purposes, permiting a great number of sound transformations. For example, time-scale modifications are accomplished by resampling the analysis points in time and results in slowing down or speeding up the sound while maintaining pitch and formant structure. Due to the stochastic and deterministic separation, this representation is more successful in time-scale modifications than others spectral representations. With it, the noise part of the sound remains “noise” no matter how much the sound is stretched.

In the deterministic representation each function pair, amplitude and frequency, accounts for a DRAFT, August 1994 27

Figure 9: Analysis representation of a vocal sound. a. Deterministic frequencies. b. Deterministic magnitudes.

DRAFT, August 1994 28

Figure 10: Analysis representation of a vocal sound. a. Stochastic magnitude. b. Stochastic coefficients.

partial of the original sound. The manipulation of these functions is easy and musically intuitive. All kinds of frequency and magnitude transformations are possible. For example, the partials can be transposed in frequency, with different values for every partial and varying during the sound. It is also possible to decouple the sinusoidal frequencies from their amplitude, obtaining effects such as changing pitch while maintaining formant structure.

The stochastic representation is modified by changing the shape of each of the envelopes and the time-varying magnitude. Changing the envelope shape corresponds to a filtering of the stochastic signal. Their manipulation is much simpler and more intuitive than the manipulation of a set of all-pole filters, such as those resulting from an LPC analysis.

Interesting effects are accomplished by changing the relative amplitude of the two components, thus emphasizing one or the other at different moments in time. Eventhough we have to realize that the characterization of a single sound by two different representations may cause problems. When different transformations are applied to each representation it is easy to create a sound in which the two components, deterministic and stochastic, do not fuse into a single entity. This may be desirable for some musical applications, but in general it is avoided, and requires some practical experimentation with the actual representations.

The best analysis/synthesis computation is generally considered the one that results in the best perceptual identity with respect to the original sound. Once this is accomplished, transformations are performed on the corresponding representation. For musical applications, however, this may not be always desirable. Very interesting effects result from purposely setting the analysis parameters “wrong”. We may, for example, set the parameters such that the deterministic analysis only captures partials in a specific frequency range, leaving the rest to be considered stochastic. The result is a sound with a much stronger noise component.

Very useful sound modifications are the result of interpolating the data from two or more sounds (Serra, 1994b). By controlling how the interpolation process is done on the different parts of the representation and in time, a large number of new sounds will result. This type of sound processing has been traditionally called cross-synthesis, eventhough a more appropiate term would be sound hybridization. With our spectral modeling method we can actually explore the timbre

Figure 11: Extraction of musical parameters from the analysis representation of a single note of an instrument.

space created by a set of sounds and define paths to go from one sound to another.

Even though this representation is powerful and many musically useful transformations are possible, we can still go further in the direction of a musically powerful representation based on analysis. The goal is to be able to control the perceptually relevant musical parameters of a sound, and the current representation is still far from that. Steps in this direction consist on extracting parameters, such as spectral shape, vibrato, overall amplitude and frequency evolution, from the current representation. These parameters can be extracted, modified and added back into the analysis data before the synthesis is done without any degradation of the resulting sound. This process is easily implemented when the input sound is a single note, in which case the musical parametrization can be quite complete. Figure 10 shows a block diagram of the steps that should be done, however a discussion of the details involved in each of the steps is beyond this presentation.

13 Deterministic Synthesis

The determnistic component is generated with additive synthesis, similar to the sinusoidal synthesis that was part of the analysis, with the difference that now the phase trajectories are discarded.

The instantaneous amplitude Â(m) of a particular partial is obtained by linear interpolation,

Â(m) = Â^{l−1} + \frac{(Â^{l} − Â^{l−1})}{S}m  (15)

where m = 0, 1, …, S − 1 is the time sample in the lth frame.

The instantaneous phase is taken to be the integral of the instantaneous frequency, where the instantaneous radian frequency ω̂(m) is obtained by linear interpolation,

ω̂(m) = ω̂^{l−1} + \frac{(ω̂^{l} − ω̂^{l−1})}{S}m  (16)

and the instantaneous phase for the rth trajectory is

θ̂_{r}(m) = θ̂_{r}(l − 1) + ω̂_{r}(m)  (17)

Finally, the synthesis equation becomes

d^{l}(m) = ∑_{r=1}^{R^{l}} Â_{r}^{l}(m) cos[θ̂_{r}^{l}(m)]  (18)

where Â(m) and θ̂(m) are the calculated instantaneous amplitude and phase.

A very efficient implementation of additive synthesis is based on the inverse-FFT (Rodet and Depalle, 1992). Eventhough this approach looses some of the flexibility of the traditional oscillator bank implementation, specially the instantaneous control of frequency and magnitude, the gain in speed is significant. This gain is based on the fact that a sinusoid in the frequency domain is a sinc-type function and on these functions not all the samples carry the same weight. To generate a sinusoid in the spectral domain it is sufficient to calculate the samples of the main lobe of the sinctype function with the appropiate magnitude, frequency and phase values. We can then synthesize as many sinusoids as we want by adding these main lobes in the FFT buffer and performing an inverse-FFT to obtain the resulting time-domain signal. By an overlap-add process we then get the time-varying characteristics of the sound.

The synthesis frame rate is completely independent of the analysis one. In the implementation using the inverse-FFT we want to have a high frame rate, so that there is no need to interpolate the frequencies and magnitudes inside a frame.

14 Stochastic Synthesis

The synthesis of the stochastic component can be understood as the generation of a noise signal that has the frequency and amplitude characteristics described by the spectral envelopes of the stochastic representation. The intuitive operation is to filter white noise with these frequency envelopes, that is, performing a time-varying filtering of white noise, which is generally implemented by the timedomain convolution of white noise with the impulse response of the filter. But in practice, the easiest and most flexible implementation is to generate the stochastic signal by an inverse-FFT of the spectral envelopes. As in the deterministic synthesis, we can then get the time-varying characteristics of the stochastic signal by an overlap-add process.

Before the inverse-FFT is performed, a complex spectrum (i.e., magnitude and phase spectra), has to be obtained from each frequency envelope. The magnitude spectrum is generated by linear interpolating the spectral envelope to a curve of length N/2, where N is the FFT-size and multiplying it by the average magnitude that was extracted in the analysis. There is no phase information in the stochastic representation, but since the phase spectrum of noise is a random signal, the phase spectrum can be created with a random number generator. To avoid a periodicity at the frame rate, different values are generated at every frame.

15 Conclusions

Modeling sounds by their time-varying spectral characteristics is a well known powerful tool. But only a few of the possible approaches are musically useful. Our discussion has been focused on the musical goal of getting a general and intuitive sound representation based on analysis, from which we can manipulate musical parameters while maintaining the perceptual identity with the original sound when no transformations are made. The sinusoids plus noise, or deterministic plus stochastic, model gives us a powerful starting point in this direction with many musical applications and possibilities for further development. We have presented the basic concepts involved in obtaining this representation and we have discussed ways to transform and synthesize sounds from the analyzed data.

An interesting direction in which to continue the work on spectral modeling using the sinusoids plus noise model is to go beyond the representation of single sounds and towards the modeling of entire timbre families. Such as all the sounds generated with an acoustic instrument. Therefore, being able to represent their common characteristics by a common set of data and keeping separate only the sound differences. As part of this process it is also important to model the articulation between notes, so that we can generate expressive phrasing based also on analysis. The result is a very powerful synthesis technique that has both the sound identity properties of sampling and the flexibility of FM.

A software package that implements most of the techniques presented in this chapter is publicly available through the Internet at the ftp server: ccrma-ftp.stanford.edu, under the name sms.tar.Z.

16 References

Allen, J.B. 1977. “Short Term Spectral Analysis, Synthesis, and Modification by Discrete Fourier Transform.” IEEE Transactions on Acoustics, Speech and Signal Processing 25(3):235–238.

Chafe, C. 1990. “Pulsed Noise in Self-sustained Oscillations of Musical Instruments.” Proceedings of the IEEE Int. Conf on Acoustics, Speech, and Signal Processing, Alburquerque, NM, 1990.

Cox, M. G. 1971. “An algorithm for approximating convex functions by means of first-degree splines.” Computer Journal, vol. 14, pp. 272–275.

Depalle, Ph., G. Garcia and X. Rodet. 1993. “Analysis of Sound for Additive Synthesis: Tracking of Partials Using Hidden Markov Models.” Proceedings of the 1993 International Computer Music Conference. San Francisco: Computer Music Association.

Doval, B., and X. Rodet. 1993. “Fundamental frequency estimation and tracking using maximum DRAFT, August 1994 34

likelihood harmonic matching and HMMs.” Proceedings of the ICASSP ’93, 221–224. Garcia G. 1992. “Analyse des Signaux Sonores en Termes de Partiels et de Bruit. Extraction

Automatique des Trajets Fréquentiels par des Modèles de Markov Cachés.” Mémoire de DEA en Automatique et Traitement du Signal, Orsay, 1992.

General Electric Co. 1977. “ADEC Subroutine Description.” Heavy Military Electronics Dept.,

Syracuse NY, 13201, June 1977. Grey, J.M. 1975. An Exploration of Musical Timbre. Ph.D. Dissertation, Stanford University. Hess, W. 1983. Pitch Determination of Speech Signals New York: Springer-Verlag. Maher, R. C. and J. W. Beauchamp. 1994. “Fundamental Frequency Estimation of Musical Signals

using a two-way Mismatch Procedure.” Journal of the Acoustical Society of America 95(4):2254– 2263.

Makhoul, J. 1975. “Linear Prediction: A Tutorial Review.” Proceedings of the IEEE 63:561–580. Markel, J.D. and A.H. Gray. 1976. Linear Prediction of Speech. New York: Springer-Verlag. McAulay, R.J. and T.F. Quatieri. 1984. “Magnitude-only Reconstruction using a Sinusoidal Speech

Model.” Proceedings of the 1984 IEEE International Conference on Acoustics, Speech and Signal Processing. New-York: IEEE Press.

McAulay, R.J. and T.F. Quatieri. 1986. “Speech Analysis/Synthesis based on a Sinusoidal Repre-

sentation.” IEEE Transactions on Acoustics, Speech and Signal Processing 34(4):744–754. Moorer, J. A. 1973. “The Hetrodyne Filter as a Tool for Analysis of Transient Waveforms.”

Memo AIM-208, Stanford Artificial Intelligence Laboratory, Computer Science Dept., Stanford University.

Moorer, J. A. 1977. “Signal Processing Aspects of Computer Music–A Survey.” Computer Music

Journal 1(1):4–37. Moorer, J. A. 1978. “The Use of the Phase Vocoder in Computer Music Applications.” Journal of

the Audio Engineering Society 26(1/2):42–45.

Phillips, G. M. 1968. “Algorithms for Piecewise Straight Line Approximation.” Computer Journal, vol. 11, pp. 211–212.

Piszczalski, M. and B. A. Galler. 1979. “Predicting Musical Pitch from Component Frequency Ratios.” Journal of the Acoustical Society of America 66(3):710–720.

Portnoff, M.R. 1976. “Implementation of the Digital Phase Vocoder Using the Fast Fourier Transform.” IEEE Transactions on Acoustics, Speech and Signal Processing 24(3):243–248.

Rodet, X. and P. Depalle. 1992. “Spectral Envelopes and Inverse FFT Synthesis.” 93rd Convention of the Audio Engineering Society. San Francisco, October 1992.

Schumacher, R. T., and C. Chafe. 1990. “Detection of Aperiodicity in Nearly Periodic Signals.” Proceedings of the IEEE Int. Conf on Acoustics, Speech, and Signal Processing, Alburquerque, NM, 1990.

Sedgewick, R. 1988. Algorithms Reading, Massachusetts: Addison-Wesley.

Serra, X. 1989. A System for Sound Analysis/Transformation/Synthesis based on a Deterministic plus Stochastic Decomposition. Ph.D. Dissertation, Stanford University.

Serra, X. and J. Smith. 1990.“Spectral Modeling Synthesis: A Sound Analysis/Synthesis System based on a Deterministic plus Stochastic Decomposition.” Computer Music Journal 14(4):12–24.

Serra, X. 1994a. “Residual Minimization in a Musical Signal Model based on a Deterministic plus Stochastic Decomposition.” Journal of the Acoustical Society of America 95(5-2):2958–2959.

Serra, X. 1994b. “Sound Hybridization Techniques based on a Deterministic plus Stochastic Decomposition Model.” Proceedings of the 1994 International Computer Music Conference. San Francisco: Computer Music Association.

Smith, J.O. and B. Friedlander. 1984. “High Resolution Spectrum Analysis Programs.” TM no. 5466-05, Systems Control Technology, Palo Alto CA, April 1984.

Smith, J.O. and X. Serra. 1987. “PARSHL: An Analysis/Synthesis Program for Non-Harmonic Sounds based on a Sinusoidal Representation.” Proceedings of the 1987 International Computer Music Conference. San Francisco: Computer Music Association.

Strawn, J. 1980. “Approximation and Syntactic Analysis of Amplitude and Frequency Functions for Digital Sound Synthesis.” Computer Music Journal 4(3):3–24.

Terhardt, E., G. Stoll, and M. Seewann. 1982. “Algorithm for Extraction of Pitch and Pitch Salience from Complex Tonal Signals.” Journal of the Acoustical Society of America 71(3):679– 688.

Wolcin J.J. 1980. “Maximum A Posteriori Line Extraction: A Computer Program.” USC TM no. 801042, March 20, 1980.

Wolcin J.J. 1980. “Maximum A Posteriori Estimation of Narrowband Signal Parameters.” NUSC TM no. 791115, June 21, 1979, Journal of the Acoustical Society of America 68(1):174–178.