# MULTI-RIPPLE LOSS FILTER FOR WAVEGUIDE PIANO SYNTHESIS

Skip other details (including permanent urls, DOI, citation information)This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License. Please contact mpub-help@umich.edu to use this work in a way not covered by the license. :

For more information, read Michigan Publishing's access and usage policy.

Page 00000001 MULTI-RIPPLE LOSS FILTER FOR WAVEGUIDE PIANO SYNTHESIS Jukka Rauhala, Heidi-Maria Lehtonen, and Vesa VSlimSki Helsinki University of Technology Laboratory of Acoustics and Audio Signal Processing P.O. Box 3000, FI-02015 TKK, Espoo, Finland ABSTRACT The multi-ripple filter is a new approach to model losses in digital waveguide synthesis. Stringed musical instruments produce decaying sounds due to physical losses, which need to be modeled. In digital waveguide string model, the losses are modeled with a loss filter. Another function for the loss filter is that it is needed to make the model stable. In piano synthesis, a high-order loss filter is needed, because of the large variation of partial decay times, which correlate to the audible sound. The multiripple loss filter is a combination of a single one-pole IIR filter and an arbitrary-order sparse FIR filter. It is able to powerfully model the differences in the partial decay rates, especially with higher orders. The multiripple loss filter can be efficiently designed by using the frequency-sampling method. Finally, a comparison between the multi-ripple loss filter and some common filters is presented showing that the multi-ripple loss filter performs well with a small computational cost. 1. INTRODUCTION Physical modeling of musical instruments is a very popular research field at the moment. One essential task in the physical modeling is to simulate the physical losses occurring in string instruments. The piano is one of the most popular musical instruments. It has a characteristic sound, which changes tonally according to time from the key press. There are a number of partials, which seem to dominate the sound after a few seconds. This must be simulated in the sound synthesis, if the goal is to produce realistic piano sounds. Digital waveguide synthesis is one of the most popular physical modeling techniques at the moment [1], [2]. In digital waveguide synthesis, the losses are modeled with a loss filter, which is included in the feedback loop together with the delay line. The goal is to match the decay rates with the magnitude response, and to have a linear or nearly linear phase response. The existing solutions use common filter design methods to design loss filters for digital waveguide synthesis. For example, a linear prediction approach is proved to produce fairly good results [1]. VSlimSki et al. proposed the "rippleOloss filter, which included a onepole filter and a feedforward comb filter [3]. The main idea of the ripple filter is to design the one-pole filter to simulate the general trend, and then to use the comb filter to tune a single partial gain accurately. The multi-ripple filter is proposed as an extension of the ripple filter, which had only one feedforward factor. The multi-ripple filter includes a one-pole IIR filter and a sparse FIR filter in cascade. In the multi-ripple filter, the only restriction for the number of feedforward paths (which is also the order of the FIR filter) is that it must be equal or smaller than the length of the delay line. When the order is increased, the quality of the produced decay rates improves. Computationally a single multiplication and one addition per sample are needed, when the order is increased by one. On the other hand, the multi-ripple filter does not need any extra memory apart from the one delay for the one-pole filter, since it can share the memory with the delay line. One of the advantages of the multi-ripple filter is that it can be designed by using the frequency-sampling method [4]. In the frequency-sampling method, the fact that the magnitude response of a single FIR comb filter is a cosine function is used to design the FIR filter coefficients with the help of the Fourier series. Moreover, a very efficient way to calculate the Fourier series coefficients is to use the fast Fourier transform (FFT). Hence, once the data for the FFT is prepared carefully, we can easily obtain the parameters for the multi-ripple filter by using the FFT of the data. 2. LOSS FILTER IN DIGITAL WAVEGUIDE SYNTHESIS Digital waveguide synthesis is a physical modeling technique, which is based on the traveling-wave solution of the wave equation [1]. Figure 1 shows a simplified digital waveguide synthesis string model, which consists of a loss filter, a delay line, a dispersion filter, and a tuning filter [5]. x(n) _Ay y(n) Tuning Dispersion Delay line Loss filter filter filter Figure 1. A basic digital waveguide string model. The loss filter design has three important goals: accurate modeling of the losses, stability, and linear phase delay response. First, the magnitude response of the filter should correspond accurately to the desired decay rates of the partials. Some researchers suggest that it is important to put more effort to model the highest gain peaks accurately [6], but it is still unclear, whether it is

Page 00000002 perceptually more important to model the highest peaks or the differences accurately. Second, the magnitude response of the loss filter should remain smaller than one in order to make the string model stable. Third, the filter should have a nearly constant phase delay response, because additional phase delay will affect the tuning of the partials. In designing the loss filter, we can take advantage of the fact that it is possible to form a closed-form formula between the decay time constant and the corresponding loss filter gain at the same frequency: gk = exp(-l/(fozk)), where gk is the gain of the loss filter, Tk is the decay time constant, fo is the fundamental frequency, and k is the partial number. 3. MULTI-RIPPLE LOSS FILTER design the feedforward filter for the whole frequency range. The transfer function of the combined multi-ripple filter and delay line can be written as N 1+ rnz-Rn + z-L H(z)= b-=1az -, (1) 1+az-1 where b is the gain factor, a is the one-pole filter parameter, N is the filter order, rn is the ripple gain, and Rn is the ripple delay length. To ensure stability, it is required that lal < 1. The ripple delay length can be denoted as Rn = round((1 - rrate,n )L), (2) where L is the length of the delay line without the multiripple filter, and rrate,n is the rate of the ripple. 3.1. Filter structure VSlimSki et al. suggested a feedforward comb filter, "ripple filterq to implement a loss filter in the digital waveguide synthesis [3]. However, we have found out that it is not suitable for piano synthesis, since it is not able to produce enough variation in the magnitude response. Hence, we propose a new extension for the ripple filter, the multi-ripple filter. Figure 2 shows the structure of the generalized Nthorder multi-ripple filter combined with the delay line. The filter includes two filters in cascade: a one-pole filter and a feedforward comb filter. The order of the feedforward filter can be anything up to the length of the delay line, but in practice we have found out that approximately filter order of 4 produces satisfactory results. 3.2. Filter design The frequency-sampling filter design method takes advantage of the property of the FIR filter that its magnitude response is a weighted sum of cosine functions [4], which are the Fourier series of an even function. Hence, the filter parameters for the desired magnitude response can be obtained by determining the Fourier series of the magnitude response, if the magnitude response is an even function. The feedforward part of the multi-ripple filter is a sparse FIR filter, which can be designed by using the frequency-sampling method. The magnitude response of the filter can be written by using the Fourier series as Hc(f)= ao +Y an cos, n=l NP) (3) where f is the frequency where the response is calculated, N, is the number of data points where the Fourier series are calculated by using a discrete transform, such as FFT, N is the filter order, and q, is the index number of the filter coefficient. When comparing to (1) and (2), we can notice that n rrate'n = N (4) Figure 2. Structure of the Nth-order multi-ripple filter combined with the delay line, b is the gain factor, a is the one-pole filter parameter, rlE rN are the ripple gain parameters, and RIE RN are the ripple delay lengths. The multi-ripple filter could be designed without the one-pole filter. However, the one-pole filter has a great advantage as it provides stability without much effort. We can select a number of partials for designing the feedforward filter and trust that the one-pole filter ensures stability on the higher frequencies, as long as the ripple gains are small and the one-pole filter has a negative coefficient. Otherwise, we would need to Additionally, we can determine that an/ao = r,, and b = (g + ao)(1 + a), since we need to get rid of the constant term. Hence, we can obtain the filter parameters from the Fourier series. 4. DESIGN STEPS Next we describe how to design the multi-ripple loss filter to match a specification obtained from a recorded piano tone.

Page 00000003 4.1. One-pole filter design The first step is to estimate the loop gain g. We do this by taking the first ten partials and by calculating the average gain. Then, we need to determine the one-pole filter parameter a. One easy way to do this is to use the least square error method. 4.2. Choose data points The next step is to calculate the difference between the desired loop gain, and the filter including the loop gain and the one-pole filter parameter (we need to calculate the b loop factor as b - g(1 + a)). Out of all partials, we need to choose the anchor points for the loss filter. We suggest that there are three types of important anchor points in the loop gain: amplitude maxima, loop gain maxima, and local maxima/minima. First, we choose the partial numbers of the two highest amplitude peaks. Second, we choose the partial numbers of the five highest loop gain values. Third, we take the first 20 partials, fit a polynomial of order 4 with Matlab's (Qolyfit' function, and differentiate the polynomial with Matlab's (Qolyder' function [7]. Then, we find the zeros in the derivative, which are the local minima and maxima, and find the corresponding partial numbers. Finally, we need to ensure that each partial number exists only once, and that the first partial is included. An example of the point selection can be seen in Figure 3. 0.994 0.98 - D 0.97 - - S0.96 - - S0.95 - 0.934 -09 - *-Desired gain O Amplitude maximum 0.92 Gain maximum 0 Local maximum/minimum 0.911111 0 10 20 30 40 50 Partial number Figure 3. An example of selected partials for Fourier series calculation. 4.3. Prepare the data for Fourier series calculation To determine the Fourier series out of the chosen data points, we need to generate a data vector. Since we have selected only a number of data points, we need to create data to fill the missing points. The goal is to find values, which generate the smallest error in the coefficients (for example, filling the gaps with zeros would create a large error). One good way is interpolation. Since most likely all the chosen data points are among the first 20 partials, we need to define some additional data points for the interpolation in order to set the general trend for the tail. We can do this by fitting a second-order polynomial on the original data, and by choosing corresponding value from point kmax+ 10, which is the maximum partial number in the chosen points. One requirement for the use of Fourier series in this case is that the input data needs to be an even function. Hence, we need to select a mirror point kmirror, which is chosen to be kmax+26, and mirror the data points so that the data is symmetric around the mirror point. Then, we can interpolate by using the chosen data points. We do this with the help of Matlab's Onterpl' function, which uses shape-preserving piecewise cubic interpolation [7]. Figure 4 shows an example of interpolated data and selected data points. 20 30 Partial number 50 Figure 4. The data that is mirrored and used to calculate the Fourier series coefficients. 4.4. Determine the Fourier series coefficients The next step is to determine the Fourier series coefficients, as seen in (3). The fastest way to do it is to compute the FFT of the data with 2(kmirror -1) points. 4.5. Determine the comb filter parameters The last task is to determine the actual filter parameters out of the Fourier series coefficients. First, we need to decide the feedforward filter order. It can be done by using a graph, which shows the squared error for each filter order, as given in Figure 5. We can find a point in the graph, where the desired performance meets the computational requirements. Fx 10-4 &-- - - 0 - --- _____ 0 10 20 30 40 50 Filter order Figure 5. The squared error of the selected data points (partials) with different filter order values.

Page 00000004 Second, we take the N highest coefficients a,, where N is the decided filter order and a, is the Fourier series coefficient n > 0. Third, we can mark the ripple gains rn to be the selected coefficient values an/ao. Fourth, we can calculate the ripple delay values Rn by using (2) and (4). Finally, the loop gain b needs to be modified by using formula b = (g + ao)(1 + a). 5. RESULTS We present a design example where we approximate the loss filter specification using the proposed filter and standard, non-sparse filters. The data for this example has been analyzed from a piano tone with a low fundamental frequency (30.9 Hz). The required loop gain for the 50 lowest harmonics is shown in Fig. 6(a) with small dots. Figure 6(b) shows the corresponding decay time, i.e., the time it takes for each harmonic to decay 8.7 dB. The characteristics of the multi-ripple filter are shown with solid line in Fig. 6. This filter consists of a one-pole filter and five feedforward comb filters. It can be seen that the filter magnitude response follows accurately the specification at the lowest few harmonics and at the highest magnitude peaks around 200 Hz and 600 Hz. The approximation is less accurate at frequencies above 1000 Hz. For comparison, we present two alternative filter designs. They were both designed using Matlab's Onvfreqz' function with the optimal error-weighting function proposed in [6]. The dash-dot line in Fig. 6 corresponds to an IIR filter with a first-order denominator and fourth-order numerator, which has the same number of multipliers as the multi-ripple filter. Obviously, this standard IIR filter cannot follow accurately the peaks of the magnitude specification. The dashed line shows the characteristics of an IIR filter with a firstorder denominator and a numerator (FIR part) of order 201. This filter is practically as good as the multi-ripple filter for this data, but the computational cost to implement it is over 30 times larger than that of the multiripple filter. This comparison shows that the power of expression of the new structure is excellent for this purpose, since it can perform well with such a small computational cost. 6. CONCLUSION A new filter structure for waveguide synthesis of string instruments was proposed. The filter consists of the cascade of a one-pole filter and a sparse FIR filter. Together they can approximate well the fluctuating magnitude specification that appears in loss filter design for waveguide piano synthesis, for example. A comparison against non-sparse loss filters designed using standard design tools with an optimal error-weighting function shows that the new structure performs well with a small computational cost. (a) 1, i---------------------------------- 0.98. 0.96 (9 0.94 0.92 - - - - 4,2 S,3 02 D1 0 200 400 600 800 1000 1200 1400 161 (b) Frequency (Hz) %. ' 0 200 400 600 800 1000 1200 1400 1600 Frequency (Hz) Figure 6. (a) The gain specification (circles) and the magnitude responses of three example filter designs: the multi-ripple (solid), the high-order conventional (dashed), and low-order conventional (dash-dot) filter. (b) The corresponding decay time data. 7. ACKNOWLEDGMENTS This work was supported by the Academy of Finland (project no. 104934). 8. REFERENCES [1] J. O. Smith III, "Techniques for Digital Filter Design and System Identification with Application to the Violin,OPh.D. thesis, Dept. of Music, Stanford Univ., Stanford, CA, 1983. [2] M. Karjalainen, V. VSlimSki, and T. Tolonen, "Plucked-string models: from the KarplusStrong algorithm to digital waveguides and beyond,O Computer Music J., vol. 22, no. 3, pp. 17-32, 1998. [3] V. VSlimSki, H. Penttinen, J. Knif, M. Laurson, and C. Erkut, "Sound synthesis of the harpsichord using a computationally efficient physical model,OEURASIP J. Applied Signal Processing, vol. 4, no. 7, pp. 934-948, 2004. [4] T. W. Parks and C. S. Burrus. Digital Filter Design. John Wiley & Sons, Inc., 1987. [5] D. A. Jaffe and J. O. Smith III, "Extensions of the Karplus-Strong plucked-string algorithm,0 Computer Music J., vol. 7, no. 2, pp. 76-87, 1983. [6] B. Bank and V. VSlimSki, "Robust loss filter design for digital waveguide synthesis,OIEEE Signal Processing Letters, vol. 10, no. 1, pp. 18-20, 2003. [7] Matlab version 6.5, Mathworks Inc. 1984 -2002.