## Loop Filter Identification

In §6.7 we discussed damping filters for vibrating string models, and in §6.9 we discussed dispersion filters. For vibrating strings which are well described by a linear, time-invariant (LTI) partial differential equation, damping and dispersion filtering are the only deviations possible from the ideal string string discussed in §6.1.

The ideal damping filter is ``zero phase'' (or linear phase)
[449],^{7.10}while the ideal dispersion filter is ``allpass'' (as described in
§6.9.1). Since
every desired frequency response can be decomposed into a zero-phase
frequency-response in series with an allpass frequency-response, we
may design a single loop filter whose amplitude response gives the
desired damping as a function of frequency, and whose phase response
gives the desired dispersion vs. frequency. The next subsection
summarizes some methods based on this approach. The following two
subsections discuss methods for the design of damping and dispersion
filters separately.

### General Loop-Filter Design

For general loop-filter design in vibrating string models (as well as in woodwind and brass instrument bore models), we wish to minimize [428, pp. 182-184]:

- Errors in partial overtone
*decay times* - Errors in partial overtone
*tuning*

There are numerous methods for designing the string loop filter
based on measurements of real string behavior. In
[428], a variety of methods for *system identification*
[288] were explored for this purpose, including
``periodic linear prediction'' in which a linear combination of a
small group of samples is used to predict a sample one period away
from the midpoint of the group. An approach based on a genetic
algorithm is described in [378]; in that work,
the error measure used with the genetic algorithm is based on
properties of human perception of short-time spectra, as is now
standard practice in digital audio coding [62].
Overviews of other approaches appear in [29] and
[508].

Below is an outline of a simple and effective method used (ca. 1995) to design loop filters for some of the Sondius sound examples:

- Estimate the fundamental frequency (see §6.11.4 below)
- Set a Hamming FFT-window length to approximately four periods
- Compute the short-time Fourier transform (STFT)
- Perform a sinusoidal modeling analysis [466] to
- -
- detect peaks in each spectral frame, and
- -
- connect peaks through time to form
*amplitude envelopes*

- Fit an exponential to each amplitude envelope
- Prepare the desired frequency-response, sampled at the harmonic
frequencies of the delay-line loop
*without*the loop filter. At each harmonic frequency,- -
- the nearest-partial decay-rate gives the desired loop-filter gains,
- -
- the nearest-partial peak-frequency give the desired loop-filter phase delay.

- Use a phase-sensitive filter-design method such as
`invfreqz`in matlab to design the desired loop filter from its frequency-response samples (further details below).

Physically, amplitude envelopes are expected to decay
*exponentially*, although coupling phenomena may obscure the
overall exponential trend. On a dB scale, exponential decay is a
*straight line*. Therefore, a simple method of estimating the
exponential decay time-constant for each overtone frequency is to fit
a straight line to its amplitude envelope and use the slope of the
fitted line to compute the decay time-constant. For example, the
matlab function `polyfit` can be used for this purpose (where
the requested polynomial order is set to 1). Since ``beating'' is
typical in the amplitude envelopes, a refinement is to replace the raw
amplitude envelope by a piecewise linear envelope that connects the
upper local maxima in the raw amplitude envelope. The estimated
decay-rate for each overtone determines a sample of the loop-filter
amplitude response at the overtone frequency. Similarly, the overtone
frequency determines a sample of the loop-filter phase response.

Taken together, the measured overtone decay rates and tunings
determine samples of the complex frequency response of the desired
loop filter. The matlab function `invfreqz`^{7.11}
can be used to convert these complex
samples into recursive filter coefficients (see §8.6.4 for a related example).
A weighting function inversely proportional to frequency is recommended. Additionally,
Steiglitz-McBride iterations can improve the results
[287], [428, pp. 101-103]. Matlab's version of
`invfreqz` has an iteration-count argument for specifying the number of
Steiglitz-McBride iterations. The maximum filter-gain versus
frequency should be computed, and the filter should be renormalized,
if necessary, to ensure that its gain does not exceed 1 at any
frequency; one good setting is that which matches the overall decay rate
of the original recording.

### Damping Filter Design

When dispersion can be neglected (as it typically can in many cases,
such as for guitar strings), we may use *linear-phase filter
design* methods, such as provided by the functions `remez`,
`firls`, and `fir2` in matlab.
Since strings are usually very lightly damped, such linear-phase
filter designs give high quality at very low orders. Another approach
is to fit a low-order IIR filter, such as by applying
`invfreqz` to a *minimum-phase* version of the desired
amplitude response, and then subtract the phase-response of the
resulting filter from the desired phase-response used in a subsequent
allpass design. (Alternatively, the phase response of the loop filter
can simply be neglected, as long as tuning is unaffected. If tuning
is affected, the tuning allpass can be adjusted to compensate.)

### Dispersion Filter Design

A pure *dispersion filter* is an ideal allpass filter. That is,
it has a gain of 1 at all frequencies and only delays a signal in a
frequency-dependent manner. The need for such filtering in
*piano string models* is discussed in §9.4.1.

There is a good amount of literature on the topic of *allpass
filter design*. Generally, they fall into the categories of
*optimized parametric*, *closed-form parametric*, and
*nonparametric* methods. Optimized parametric methods can
produce allpass filters with optimal group-delay characteristics in
some sense [272,271]. Closed-form parametric
methods provide coefficient formulas as a function of a desired
parameter such as ``inharmonicity'' [368].
Nonparametric methods are generally based on measured signals and/or
spectra, and while they are suboptimal, they can be used to design
very large-order allpass filters, and the errors can usually be made
arbitrarily small by increasing the order
[551,369,42,41,1],
[428, pp. 60,172]. In music applications, it is usually the
case that the ``optimality'' criterion is unknown because it depends
on aspects of sound perception (see, for example,
[211,384]). As a result, perceptually
weighted nonparametric methods can often outperform optimal parametric
methods in terms of cost/performance [2].

In historical order, some of the allpass filter-design methods are as follows: A modification of the method in [551] was suggested for designing allpass filters having a phase delay corresponding to the delay profile needed for a stiff string simulation [428, pp. 60,172]. The method of [551] was streamlined in [369]. In [77], piano strings were modeled using finite-difference techniques. An update on this approach appears in [45]. In [340], high quality stiff-string sounds were demonstrated using high-order allpass filters in a digital waveguide model. In [384], this work was extended by applying a least-squares allpass-design method [272] and a spectral Bark-warping technique [459] to the problem of calibrating an allpass filter of arbitrary order to recorded piano strings. They were able to correctly tune the first several tens of partials for any natural piano string with a total allpass order of 20 or less. Additionally, minimization of the norm [271] has been used to calibrate a series of allpass-filter sections [42,41], and a dynamically tunable method, based on Thiran's closed-form, maximally flat group-delay allpass filter design formulas (§4.3), was proposed in [368]. An improved closed-form solution appears in [1] based on an elementary method for the robust design of very high-order allpass filters.

### Fundamental Frequency Estimation

As mentioned in §6.11.2 above, it is advisable to
estimate the fundamental frequency of vibration (often called
``*F0*'') in order that the partial overtones are well resolved
while maintaining maximum time resolution for estimating the decay
time-constant.

Below is a summary of the F0 estimation method used in calibrating loop filters with good results [471]:

- Take an FFT of the middle third of a recorded plucked string tone.
- Find the frequencies and amplitudes of the largest peaks, where is chosen so that the retained peaks all have a reasonable signal-to-noise ratio.
- Form a histogram of peak spacing
- The pitch estimate is defined as the most common spacing in the histogram.

#### Approximate Maximum Likelihood F0 Estimation

In applications for which the fundamental frequency F0 must be measured very accurately in a periodic signal, the estimate obtained by the above algorithm can be refined using a gradient search which matches a so-called ``harmonic comb'' to the magnitude spectrum of an interpolated FFT :

Note that freely vibrating strings are not exactly periodic due to exponenential decay, coupling effects, and stiffness (which stretches harmonics into quasiharmonic overtones, as explained in §6.9). However, non-stiff strings can often be analyzed as having approximately harmonic spectra ( periodic time waveform) over a limited time frame.

Since string spectra typically exhibit harmonically spaced
*nulls* associated
with the excitation and/or observation points, as well as from other
phenomena such as recording multipath and/or reverberation, it is
advisable to restrict to a range that does not include any
spectral nulls (or simply omit index when
is
too close to a spectral null),
since even one spectral null can push the product of
peak amplitudes to a very small value. As a practical matter, it is
important to inspect the magnitude spectra of the data manually to
ensure that a robust row of peaks is being matched by the harmonic
comb. For example, a display of the frame magnitude spectrum overlaid
with vertical lines at the optimized harmonic-comb frequencies yields
an effective picture of the F0 estimate in which typical problems
(such as octave errors) are readily seen.

#### References on F0 Estimation

An often-cited book on classical methods for pitch detection,
particularly for voice, is that by Hess [192]. The harmonic
comb method can be considered an approximate maximum-likelihood pitch
estimator, and more accurate maximum-likelihood methods have been
worked out
[114,547,376,377].
More recently, Klapuri has been developing some promising methods for
multiple pitch estimation
[254,253,252].^{7.12}A comparison of real-time pitch-tracking algorithms applied to guitar
is given in [260], with consideration of
latency (time delay).

#### Extension to Stiff Strings

An advantage of the harmonic-comb method, as well as other frequency-domain maximum-likelihood pitch-estimation methods, is that it is easily extended to accommodate stiff strings. For this, the stretch-factor in the spectral-peak center-frequencies can be estimated--the so-called coefficient of inharmonicity, and then the harmonic-comb (or other maximum-likelihood spectral-matching template) can be stretched by the same amount, so that when set to the correct pitch, the template matches the data spectrum more accurately than if harmonicity is assumed.

### EDR-Based Loop-Filter Design

This section discusses use of the Energy Decay Relief (EDR) (§3.2.2) to measure the decay times of the partial overtones in a recorded vibrating string.

First we derive what to expect in the case of a simplified string model along the lines discussed in §6.7 above. Assume we have the synthesis model of Fig.6.12, where now the loss factor is replaced by the digital filter that we wish to design. Let denote the contents of the delay line as a vector at time , with denoting the initial contents of the delay line.

For simplicity, we define the EDR based on a length DFT of the delay-line
vector
, and use a rectangular window with a ``hop size'' of samples,
*i.e.*,

Applying the definition of the EDR (§3.2.2) to this short-time spectrum gives

We therefore have the following recursion for successive EDR
time-slices:^{7.13}

This analysis can be generalized to a time-varying model in which the
loop filter is allowed to change once per ``period''
.^{7.14}

An online laboratory exercise covering the practical details of measuring overtone decay-times and designing a corresponding loop filter is given in [280].

**Next Section:**

String Coupling Effects

**Previous Section:**

The Externally Excited String