<?xml version="1.0" encoding="utf-8" standalone="no"?>
<!DOCTYPE document PUBLIC "-//CNX//DTD CNXML 0.5 plus MathML//EN" "http://cnx.rice.edu/cnxml/0.5/DTD/cnxml_mathml.dtd">
<document xmlns="http://cnx.rice.edu/cnxml" xmlns:md="http://cnx.rice.edu/mdml/0.4" xmlns:m="http://www.w3.org/1998/Math/MathML" xmlns:bib="http://bibtexml.sf.net/" id="Module.2003-11-26.4600">
  <name>Delays and Effects</name>
  <metadata>
  <md:version>1.7</md:version>
  <md:created>2003/11/26 06:46:00 US/Central</md:created>
  <md:revised>2004-02-24T17:08:21Z</md:revised>
  <md:authorlist>
    <md:author id="drocchesso">
      <md:firstname>Davide</md:firstname>
      
      <md:surname>Rocchesso</md:surname>
      <md:email>Davide.Rocchesso@univr.it</md:email>
    </md:author>
  </md:authorlist>

  <md:maintainerlist>
    <md:maintainer id="drocchesso">
      <md:firstname>Davide</md:firstname>
      
      <md:surname>Rocchesso</md:surname>
      <md:email>Davide.Rocchesso@univr.it</md:email>
    </md:maintainer>
  </md:maintainerlist>
  
  <md:keywordlist>
    <md:keyword>delay lines</md:keyword>
    <md:keyword>comb filter</md:keyword>
  </md:keywordlist>

  <md:abstract>Delay lines and comb filters are the key elements of many digital audio effects. This module explains what they are and how they are typically implemented.</md:abstract>
</metadata>

  <content>
  <note>
  Extracted from 
    <link src="http://www.theassayer.org/cgi-bin/asbook.cgi?book=522"> Introduction to Sound Processing
    </link>
  </note>
    <para id="intro0">
      Most acoustic systems have some component where waves can propagate, such as a membrane, a string, or the air in an enclosure.
      If propagation in these media is ideal, i.e., free of losses, dispersion, and nonlinearities, it can be simulated by delay lines.
      A delay line is a linear time-invariant, single-input single-output system, whose output signal is a copy of the input signal delayed by 
      <m:math>
	<m:ci>τ</m:ci>
      </m:math>
      seconds. In continuous time, the frequency response of such system is
    </para>
    <equation id="eq1">
      <m:math display="block">
	<m:apply>
	  <m:eq/>
	  <m:apply>
	    <m:ci type="fn"><m:msub>
	      <m:mi>H</m:mi>
		<m:msub>
		  <m:mi>D</m:mi>
		  <m:mi>s</m:mi>
		</m:msub>
	      </m:msub></m:ci>
	    <m:apply>
	      <m:times/>
	      <m:imaginaryi/>
	      <m:ci>Ω</m:ci>
	    </m:apply>
	  </m:apply>
	  <m:apply>
	    <m:exp/>
	    <m:apply>
	      <m:minus/>
	      <m:apply>
		<m:times/>
		<m:imaginaryi/>
		<m:ci>Ω</m:ci>
		<m:ci>τ</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:apply>
	</m:apply>
      </m:math>
    </equation>
    <para id="intro1">
      <cnxn target="eq1"/> tells us that the magnitude response is
      unitary, and that the phase is linear with slope
      <m:math>
	<m:ci>τ</m:ci>
      </m:math>.
    </para>
		
    <section id="par31">
      <name>The Circular Buffer</name>
      <para id="par3100">
	A discrete-time realization of the <cnxn target="eq1"/> is given by a system that implements the transfer function
	<equation id="eq3">
	  <m:math display="block">
           <m:mrow>
	  <m:msub>
	  <m:mi>H</m:mi>
	  <m:mi>D</m:mi>
	</m:msub>
	  <m:mo stretchy="false">(</m:mo>
	  <m:mi>z</m:mi>
	  <m:mo stretchy="false">)</m:mo>
	  <m:mo>=</m:mo>
	  <m:msup>
	  <m:mi>z</m:mi>
	  <m:mrow>
	  <m:mo>−</m:mo>
	  <m:mi>τ</m:mi>
	  <m:msub>
	  <m:mi>F</m:mi>
	  <m:mi>s</m:mi>
	</m:msub>
	</m:mrow>
	</m:msup>
	  <m:mo>≜</m:mo>
	  <m:msup>
	  <m:mi>z</m:mi>
	  <m:mrow>
	  <m:mo>−</m:mo>
	  <m:mi>m</m:mi>
	</m:mrow>
	</m:msup>
	</m:mrow>
<!--	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:ci type="fn"><m:msub>
		    <m:mi>H</m:mi>
		    <m:mi>D</m:mi>
		  </m:msub></m:ci>
		<m:ci>z</m:ci>
	      </m:apply>
	      <m:apply>
		<m:mo>≜</m:mo>
		<m:apply>
		  <m:power/>
		  <m:ci>z</m:ci>
		  <m:apply>
		    <m:minus/>
		    <m:apply>
		      <m:times/>
		      <m:ci>τ</m:ci>
		      <m:ci><m:msub>
			  <m:mi>F</m:mi>
			  <m:mi>s</m:mi>
			</m:msub></m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
		<m:apply>
		  <m:power/>
		  <m:ci>z</m:ci>
		  <m:apply>
		    <m:minus/>
		    <m:ci>m</m:ci>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:apply> -->
	  </m:math>
	</equation>
	where
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	is the number of samples of delay. When the delay
	<m:math>
	  <m:ci>τ</m:ci>
	</m:math>
	is an integral multiple of the sampling quantum, 
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	is an integer number and it is straightforward
	to implement the <cnxn target="eq3"/> by means of a memory buffer.
      </para>
      <para id="par3102">
	In fact, an
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>-samples delay line can be implemented by means of a circular buffer, that is a set of
	<m:math>
	  <m:ci>M</m:ci>
	</m:math>
	contiguous memory cells accessed by a write pointer
	<m:math>
	  <m:mrow>
	    <m:mtext>IN</m:mtext>
	  </m:mrow>
	</m:math>
	and a read pointer
	<m:math>
	  <m:mrow>
	    <m:mtext>OUT</m:mtext>
	  </m:mrow>
	</m:math>, such that
	<equation id="eq5">
	  <m:math display="block">
	    <m:mtext>IN</m:mtext>
	    <m:mo>=</m:mo>
	    <m:mo stretchy="false">(</m:mo>
	    <m:mtext>OUT</m:mtext>
	    <m:mo>+</m:mo>
	    <m:mi>m</m:mi>
	    <m:mo stretchy="false">)</m:mo>
	    <m:mo>%</m:mo>
	    <m:mi>M</m:mi>
	  </m:math>
	</equation>
      </para>
      <para id="par3103">
	where the symbol
	<m:math>
	  <m:ci>ρ</m:ci>
	</m:math>
	is used for the quotient modulo
	<m:math>
	  <m:ci>M</m:ci>
	</m:math>. At each sampling instant, the input is written in the location pointed by
	<m:math>
	  <m:mrow>
	    <m:mtext>IN</m:mtext>
	  </m:mrow>
	</m:math>, the output is taken
	from the location pointed by
	<m:math>
	  <m:mrow>
	    <m:mtext>OUT</m:mtext>
	  </m:mrow>
	</m:math>, and the two pointers are updated with
	<equation id="eq6">
	  <m:math display="block">
	    <m:mtable columnalign="left">
	      <m:mtr>
		<m:mtd>
		  <m:mtext>IN</m:mtext>
		  <m:mo>=</m:mo>
		  <m:mo stretchy="false">(</m:mo>
		  <m:mtext>IN</m:mtext>
		  <m:mo>+</m:mo>
		  <m:mn>1</m:mn>
		  <m:mo stretchy="false">)</m:mo>
		  <m:mo>%</m:mo>
		  <m:mi>M</m:mi>
		</m:mtd>
	      </m:mtr>
	      <m:mtr>
		<m:mtd>
		  <m:mtext>OUT</m:mtext>
		  <m:mo>=</m:mo>
		  <m:mo stretchy="false">(</m:mo>
		  <m:mtext>OUT</m:mtext>
		  <m:mo>+</m:mo>
		  <m:mn>1</m:mn>
		  <m:mo stretchy="false">)</m:mo>
		  <m:mo>%</m:mo>
		  <m:mi>M</m:mi>
		</m:mtd>
	      </m:mtr>
	    </m:mtable>
	  </m:math>
	</equation>
      </para>
      <para id="par3104">
	In words, the pointers are incremented respecting the circularity of the buffer.
	In some architectures dedicated to sound processing, memory organization
	is optimized for wavetable synthesis, where a stored waveform is read with
	variable increments of the reading pointer. In these architectures, a quantity of
	<m:math>
	  <m:apply>
	    <m:power/>
	    <m:cn>2</m:cn>
	    <m:ci>r</m:ci>
	  </m:apply>
	</m:math>
	memory locations is available, and from these
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>M</m:ci>
	    <m:apply>
	      <m:power/>
	      <m:cn>2</m:cn>
	      <m:ci>s</m:ci>
	    </m:apply>
	  </m:apply>
	</m:math>
	locations (with
	<m:math>
	  <m:apply>
	    <m:lt/>
	    <m:ci>s</m:ci>
	    <m:ci>r</m:ci>
	  </m:apply>
	</m:math>) are uniformly chosen among the
	<m:math>
	  <m:apply>
	    <m:power/>
	    <m:cn>2</m:cn>
	    <m:ci>r</m:ci>
	  </m:apply>
	</m:math>
	available cells. In this case the locations of the circular
	buffer are not contiguous, and the update of the pointers is
	done with the operations
	<equation id="eq7">
	  <m:math display="block">
	    <m:mtable columnalign="left">
	      <m:mtr>
		<m:mtd>
		  <m:mtext>IN</m:mtext>
		  <m:mo>=</m:mo>
		  <m:mo stretchy="false">(</m:mo>
		  <m:mtext>IN</m:mtext>
		  <m:mo>+</m:mo>
		  <m:msup>
		    <m:mn>2</m:mn>
		    <m:mrow>
		      <m:mi>r</m:mi>
		      <m:mo>−</m:mo>
		      <m:mi>s</m:mi>
		    </m:mrow>
		  </m:msup>
		  <m:mo stretchy="false">)</m:mo>
		  <m:mo>%</m:mo>
		  <m:msup>
		    <m:mn>2</m:mn>
		    <m:mi>r</m:mi>
		  </m:msup>
		</m:mtd>
	      </m:mtr>
	      <m:mtr>
		<m:mtd>
		  <m:mtext>OUT</m:mtext>
		  <m:mo>=</m:mo>
		  <m:mo stretchy="false">(</m:mo>
		  <m:mtext>OUT</m:mtext>
		  <m:mo>+</m:mo>
		  <m:msup>
		    <m:mn>2</m:mn>
		    <m:mrow>
		      <m:mi>r</m:mi>
		      <m:mo>−</m:mo>
		      <m:mi>s</m:mi>
		    </m:mrow>
		  </m:msup>
		  <m:mo stretchy="false">)</m:mo>
		  <m:mo>%</m:mo>
		  <m:msup>
		    <m:mn>2</m:mn>
		    <m:mi>r</m:mi>
		  </m:msup>
		</m:mtd>
	      </m:mtr>
	    </m:mtable>
	  </m:math>
	</equation>
      </para>
      <para id="par3105">
	In practice, since the addresses are 
	<m:math>
	  <m:ci>r</m:ci>
	</m:math>-bit long, there is no need to compute
	the modulo explicitly. It is sufficient to do the sum neglecting any possible
	overflow.
      </para>
      <para id="par3106">
	Of course, the <cnxn target="eq5"/> is also replaced by
	<equation id="eq9">
	  <m:math display="block">
	    <m:mrow>
	      <m:mtext>IN</m:mtext>
	      <m:mo>=</m:mo>
	      <m:mo stretchy="false">(</m:mo>
	      <m:mtext>OUT</m:mtext>
	      <m:mo>+</m:mo>
	      <m:mi>m</m:mi>
	      <m:msup>
		<m:mn>2</m:mn>
		<m:mrow>
		  <m:mi>r</m:mi>
		  <m:mo>−</m:mo>
		  <m:mi>s</m:mi>
		</m:mrow>
	      </m:msup>
	      <m:mo stretchy="false">)</m:mo>
	      <m:mo>%</m:mo>
	      <m:msup>
		<m:mn>2</m:mn>
		<m:mi>r</m:mi>
	      </m:msup>
	    </m:mrow>
	  </m:math>
	</equation>
      </para>
    </section>
    <section id="par32">
      <name>Fractional-Length Delay Lines</name>
      <para id="par3200">
	It might be thought that, choosing a sufficiently high sampling rate, it is
	always possible to use delay lines having an integer number of samples. Actually,
	there are some good reasons that lead us to state that this is not the case in
	sound synthesis and processing.
	In sound synthesis, the models have to be carefully tuned without resorting
	to very high sample rates. In particular, it is easy to verify that using integer length
	delays in physical models we get errors in fundamental frequencies that
	go well beyond the just noticeable difference in pitch.
	<note id="n1" type="note">To figure this out,
	  the reader can consider an
	  <m:math><m:ci>m</m:ci></m:math>-sample
	  delay line in a feedback loop. It
	  gives a harmonic series of partials
	  whose fundamental is
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:ci><m:msub> 
		  <m:mi>f</m:mi>
		  <m:mn>0</m:mn>
		</m:msub></m:ci>
	      <m:apply>
		<m:divide/>
		<m:ci><m:msub>
		    <m:mi>F</m:mi>
		    <m:mi>s</m:mi>
		  </m:msub></m:ci>
		<m:ci>m</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	  (see <cnxn target="par34"/>). The set of integer
	  delay lengths that give the best approximation to a tempered scale can be found and the curve of
	  fundamental frequency errors can be drawn.</note>
      </para>
      <para id="par3201">
	For instance, for a pressure wave propagating in air at normal temperature conditions,
	the spatial discretization given by the sampling rate
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci><m:msub>
		<m:mi>F</m:mi>
		<m:mi>s</m:mi>
	      </m:msub></m:ci>
	    <m:apply>
	      <m:times/>
	      <m:cn>44100</m:cn>
	      <m:ci>Hz</m:ci>
	    </m:apply>
	  </m:apply>
	</m:math>
	gives intervals of
	<m:math>
	  <m:apply>
	    <m:times/>
	    <m:cn>0.0075</m:cn>
	    <m:ci>m</m:ci>
	  </m:apply>
	</m:math>, a distance that can produce well-perceivable pitch
	differences in a wind instrument.
	Another reason for using fractional delays is that we often want to vary
	the delay lengths continuously, in order to reproduce effects such as glissando
	or vibrato. The adoption of integer-length delays would produce annoying discontinuities.
	The most widely used techniques for implementing fractional delays are
	interpolation by FIR filters or by allpass filters. These two techniques are, in
	some sense, complementary. The choice of one of the two has to be made according
	to the peculiarities of the system to be simulated or of the architecture
	chosen for the implementation. In any case, a delay of length
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	is obtained by means of a delay line whose length is equal to the integer part of
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>, cascaded with a block capable to approximate a constant phase delay equal to the fractional part of
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>. We recall that the phase delay at a given frequency
	<m:math>
	  <m:ci>ω</m:ci>
	</m:math>
	is the delay in time samples experienced by the sinusoidal component at frequency
	<m:math>
	  <m:ci>ω</m:ci>
	</m:math>. For instance, consider a linear filtering block enclosed in a feedback loop
	(see <cnxn target="par34"/>): the frequency of the
	<m:math>
	  <m:ci>k</m:ci>
	</m:math>-th resonance
	<m:math>
	  <m:ci><m:msub>
	      <m:mi>f</m:mi>
	      <m:mi>k</m:mi>
	    </m:msub></m:ci>
	</m:math>
	of the whole feedback system is found at the points where the phase response equates the multiples of
	<m:math>
	  <m:apply>
	    <m:times/>
	    <m:cn>2</m:cn>
	    <m:pi/>
	  </m:apply>
	</m:math>. At these frequencies, the components reappear in phase every round trip
	in the loop, thus reinforcing their amplitude at the output. The phase delay at frequency
	<m:math>
	  <m:ci><m:msub>
	      <m:mi>f</m:mi>
	      <m:mi>k</m:mi>
	    </m:msub></m:ci>
	</m:math>
	is therefore the effective delay length at that frequency, that is
	the length of an ideal (linear phase) delay line that gives the same
	<m:math>
	  <m:ci>k</m:ci>
	</m:math>-th resonance. <cnxn target="fig1"/> shows a phase curve and its crossings with multiples of
	<m:math>
	  <m:apply>
	    <m:times/>
	    <m:cn>2</m:cn>
	    <m:pi/>
	  </m:apply>
	</m:math>
	giving a distribution of resonances.
      </para>
      <section id="par321">
	<name>FIR Interpolation Filters</name>
	<para id="par3210">
	  The easiest and most intuitive way to obtain a variable-length delay is to
	  linearly interpolate the output of the line with the content of its preceding cell
	  in the memory buffer. This corresponds to using the first-order FIR filter
	  <equation id="eq10">
	    <m:math display="block">
	      <m:apply>
		<m:eq/>
		<m:apply>
		  <m:ci type="fn"><m:msub>
		      <m:mi>H</m:mi>
		      <m:mi>l</m:mi>
		    </m:msub></m:ci>
		  <m:ci>z</m:ci>
		</m:apply>
		<m:apply>
		  <m:plus/>
		  <m:ci><m:msub>
		      <m:mi>c</m:mi>
		      <m:mn>0</m:mn>
		    </m:msub></m:ci>
		  <m:apply>
		    <m:times/>
		    <m:ci><m:msub>
			<m:mi>c</m:mi>
			<m:mn>1</m:mn>
		      </m:msub></m:ci>
		    <m:apply>
		      <m:inverse/>
		      <m:ci>z</m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	  <figure id="fig1">
            <media type="image/gif" src="fase.gif"/>
	    <caption>Graphical construction to find the series of resonances produced by
	      a linear block in a feedback loop. The slope of the dashed lines indicates the
	      phase delay at each resonance frequency.</caption>
	  </figure>
	</para>
	<para id="par3211">
	  Given a certain phase delay
	  <equation id="eq11">
	    <m:math display="block">
	      <m:apply>
		<m:eq/>
		<m:ci><m:msub>
		    <m:mi>τ</m:mi>
		    <m:msub>
		      <m:mi>ph</m:mi>
		      <m:mn>0</m:mn>
		    </m:msub>
		  </m:msub></m:ci>
		<m:apply>
		  <m:minus/>
		  <m:apply>
		    <m:times/>
		    <m:apply>
		      <m:divide/>
		      <m:cn>1</m:cn>
		      <m:ci><m:msub>
			  <m:mi>ω</m:mi>
			  <m:mn>0</m:mn>
			</m:msub></m:ci>
		    </m:apply>
		    <m:apply>
		      <m:arctan/>
		      <m:apply>
			<m:divide/>
			<m:apply>
			  <m:minus/>
			  <m:apply>
			    <m:times/>
			    <m:ci><m:msub>
				<m:mi>c</m:mi>
				<m:mn>1</m:mn>
			      </m:msub></m:ci>
			    <m:apply>
			      <m:sin/>
			      <m:ci><m:msub>
				  <m:mi>ω</m:mi>
				  <m:mn>0</m:mn>
				</m:msub></m:ci>
			    </m:apply>
			  </m:apply>
			</m:apply>
			<m:apply>
			  <m:plus/>
			  <m:ci><m:msub>
			      <m:mi>c</m:mi>
			      <m:mn>0</m:mn>
			    </m:msub></m:ci>
			  <m:apply>
			    <m:times/>
			    <m:ci><m:msub>
				<m:mi>c</m:mi>
				<m:mn>1</m:mn>
			      </m:msub></m:ci>
			    <m:apply>
			      <m:cos/>
			      <m:ci><m:msub>
				  <m:mi>ω</m:mi>
				  <m:mn>0</m:mn>
				</m:msub></m:ci>
			    </m:apply>
			  </m:apply>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	  that has to be obtained at a given frequency
	  <m:math>
	    <m:ci><m:msub>
		<m:mi>ω</m:mi>
		<m:mn>0</m:mn>
	      </m:msub></m:ci>
	  </m:math>, the following formulas give the coefficient values:
	  <equation id="eq12">
	    <m:math display="block">
	      <m:mtable columnalign="left">
		<m:mtr>
		  <m:mtd>
		    <m:apply>
		      <m:eq/>
		      <m:apply>
			<m:plus/>
			<m:ci><m:msub>
			    <m:mi>c</m:mi>
			    <m:mn>0</m:mn>
			  </m:msub></m:ci>
			<m:ci><m:msub>
			    <m:mi>c</m:mi>
			    <m:mn>1</m:mn>
			  </m:msub></m:ci>
		      </m:apply>
		      <m:cn>1</m:cn>
		    </m:apply>
		  </m:mtd>
		</m:mtr>
		<m:mtr>
		  <m:mtd>
		    <m:apply>
		      <m:eq/>
		      <m:ci><m:msub>
			  <m:mi>c</m:mi>
			  <m:mn>1</m:mn>
			</m:msub></m:ci>
		      <m:apply>
			<m:approx/>
			<m:apply>
			  <m:divide/>
			  <m:cn>1</m:cn>
			  <m:apply>
			    <m:minus/>
			    <m:apply>
			      <m:plus/>
			      <m:cn>1</m:cn>
			      <m:apply>
				<m:divide/>
				<m:apply>
				  <m:sin/>
				  <m:ci><m:msub>
				      <m:mi>ω</m:mi>
				      <m:mn>0</m:mn>
				    </m:msub></m:ci>
				</m:apply>
				<m:apply>
				  <m:tan/>
				  <m:apply>
				    <m:times/>
				    <m:ci><m:msub>
					<m:mi>τ</m:mi>
					<m:msub>
					  <m:mi>ph</m:mi>
					  <m:mn>0</m:mn>
					</m:msub>
				      </m:msub></m:ci>
				    <m:ci><m:msub>
					<m:mi>ω</m:mi>
					<m:mn>0</m:mn>
				      </m:msub></m:ci>
				  </m:apply>
				</m:apply>
			      </m:apply>
			    </m:apply>
			    <m:apply>
			      <m:cos/>
			      <m:ci><m:msub>
				  <m:mi>ω</m:mi>
				  <m:mn>0</m:mn>
				</m:msub></m:ci>
			    </m:apply>
			  </m:apply>
			</m:apply>
			<m:ci><m:msub>
			    <m:mi>τ</m:mi>
			    <m:msub>
				<m:mi>ph</m:mi>
				<m:mn>0</m:mn>
			      </m:msub>
			  </m:msub></m:ci>
		      </m:apply>
		    </m:apply>
		  </m:mtd>
		</m:mtr>
	      </m:mtable>
	    </m:math>
	  </equation>
	  where the approximation is valid in the low-frequency range. The first part of <cnxn target="eq12"/>
	  is needed in order to normalize the low-frequency response to one. In the special case that
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:ci><m:msub>
		  <m:mi>c</m:mi>
		  <m:mn>0</m:mn>
		</m:msub></m:ci>
	      <m:ci><m:msub>
		  <m:mi>c</m:mi>
		  <m:mn>1</m:mn>
		</m:msub></m:ci>
	      <m:cn type="rational">1<m:sep/>2</m:cn>
	    </m:apply>
	  </m:math>
	  (averaging filter) the phase is linear and the delay
	  is of half a sample. Unfortunately, the magnitude response of this interpolator
	  is lowpass with a zero at the Nyquist frequency. <cnxn target="fig2"/>
	  shows the magnitude, phase, and phase delay responses for several first-order linear interpolators.  We
	  can see that the phase is linear in most of the audio range, but the magnitude
	  varies from the allpass to the lowpass with a zero at the Nyquist rate. When
	  the interpolator is inserted within a feedback loop, its lowpass behavior can be
	  treated as an additional frequency-dependent loss, which should be somewhat
	  taken into account.
	  <figure id="fig2">
	    <media type="image/gif" src="linintresp.gif"/>
	    <caption>Magnitude, phase, and phase delay responses of a linear interpolation filter
	      <m:math>
		<m:apply>
		  <m:plus/>
		  <m:apply>
		    <m:minus/>
		    <m:cn>1</m:cn>
		    <m:ci>α</m:ci>
		  </m:apply>
		  <m:apply>
		    <m:times/>
		    <m:mi>α</m:mi>
		    <m:apply>
		      <m:inverse/>
		      <m:ci>z</m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:math> for 
	      <m:math>
		<m:apply>
		  <m:in/>
		  <m:apply>
		    <m:eq/>
		    <m:ci>α</m:ci>
		    <m:ci>k</m:ci>
		  </m:apply>
		  <m:set>
		    <m:cn>0</m:cn>
		    <m:ci>…</m:ci>
		    <m:cn>16</m:cn>
		  </m:set>
		</m:apply>
	      </m:math>
	    </caption>
	  </figure>
	</para>
	<para id="par3214">
	  Interpolation filters can be of order higher than the first. We can do quadratic,
	  cubic, or other polynomial interpolations. In general, the problem of
	  designing an interpolator can be turned into the design of an
	  <m:math>
	    <m:ci>l</m:ci>
	  </m:math>-th order FIR filter approximating a constant-magnitude and linear-phase frequency response. Several
	  criteria can be adopted to drive the approximation problem. One approach is 
	  to impose that the first
	  <m:math>
	    <m:ci>L</m:ci>
	  </m:math>
	  derivatives of the error function will be zero at zero
	  frequency. In this way we obtain maximally-flat filters whose coefficients are
	  the same used in Lagrange interpolation as it is taught in numerical analysis
	  courses. For a thorough treatment of interpolation filters we suggest reading the article
	  <cite src="#bib1"/>. Here we only point out that using high orders allows to keep
	  the magnitude response close to unity and a phase response close to linear in
	  a wide frequency band. Of course, this is paid in terms of computational complexity.
	  In special architectures, where the access to delay lines is governed by
	  <cnxn target="eq7"/> and
	  <cnxn target="eq9"/>, the linear interpolation is implemented very efficiently by using the
	  <m:math>
	    <m:apply>
	      <m:minus/>
	      <m:ci>r</m:ci>
	      <m:ci>s</m:ci>
	    </m:apply>
	  </m:math>
	  bits that are not used to access the
	  <m:math>
	    <m:apply>
	      <m:power/>
		<m:cn>2</m:cn>
		<m:ci>s</m:ci>
	    </m:apply>
	  </m:math>
	  samples delay line. In fact, if the address is computed using
	  <m:math>
	    <m:ci>r</m:ci>
	  </m:math>
	  bits, the
	  <m:math>
	    <m:apply>
	      <m:minus/>
	      <m:ci>r</m:ci>
	      <m:ci>s</m:ci>
	    </m:apply>
	  </m:math>
	  least significant bits represent the fractional part of the delay or, equivalenty, the coefficient
	  <m:math>
	    <m:ci><m:msub>
		<m:mi>c</m:mi>
		<m:mn>1</m:mn>
	      </m:msub></m:ci>
	  </m:math>
	  of the interpolator. Therefore, it is sufficient to access two consecutive delay cells and keep the values	
	  <m:math>
	    <m:ci><m:msub>
		<m:mi>c</m:mi>
		<m:mn>0</m:mn>
	      </m:msub></m:ci>
	  </m:math>
	  and
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:ci><m:msub>
		  <m:mi>c</m:mi>
		  <m:mn>1</m:mn>
		</m:msub></m:ci>
	      <m:apply>
		<m:minus/>
		<m:cn>1</m:cn>
		<m:ci><m:msub>
		    <m:mi>c</m:mi>
		    <m:mn>0</m:mn>
		  </m:msub></m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	  in two registers. The implementation of a glissando with these architectures is immediate and free from complications.
	</para>
      </section>

      <section id="par322">
	<name>Allpass Interpolation Filters</name>
	<para id="par3220">
	  Another widely used technique to obtain the fractional part of a desired
	  delay length makes use of unit-magnitude IIR filters, i.e., allpass filters. Since
	  the magnitude of these filters is constant there is no frequency-dependent attenuation,
	  a property that can never be ensured by FIR filters. The simplest allpass
	  filter has order one, and it has the following transfer function:
	  <equation id="eq15">
	    <m:math display="block">
	      <m:apply>
		<m:eq/>
		<m:apply>
		  <m:ci type="fn"><m:msub>
		      <m:mi>H</m:mi>
		      <m:mi>a</m:mi>
		    </m:msub></m:ci>
		  <m:ci>z</m:ci>
		</m:apply>
		<m:apply>
		  <m:divide/>
		  <m:apply>
		    <m:plus/>
		    <m:ci>c</m:ci>
		    <m:apply>
		      <m:inverse/>
		      <m:ci>z</m:ci>
		    </m:apply>
		  </m:apply>
		  <m:apply>
		    <m:plus/>
		    <m:cn>1</m:cn>
		    <m:apply>
		      <m:times/>
		      <m:ci>c</m:ci>
		      <m:apply>
			<m:inverse/>
			<m:ci>z</m:ci>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	</para>
	<para id="par3221">
	  In order to make sure that the filter is stable, the coefficient
	  <m:math>
	    <m:ci>c</m:ci>
	  </m:math>
	  has to stay within the unit circle. Moreover, if we stick with real coefficients,
	  <m:math>
	    <m:ci>c</m:ci>
	  </m:math>
	  belongs to the real axis. The phase delay given by the filter <cnxn target="eq15"/>
	  is shown in <cnxn target="fig3"/> for several values of the coefficient
	  <m:math>
	    <m:ci>c</m:ci>
	  </m:math>.  It is clear that the phase delay is not as flat as in the case of the FIR interpolator, depicted in
	  <cnxn target="fig2"/>.
	</para>
	<figure id="fig3">
	  <media type="image/gif" src="apresp.gif"/>
	  <caption>Phase response and phase delay of a first-order allpass filter for the values of the coefficient
	    <m:math>
	      <m:apply>
		<m:eq/>
		<m:ci>c</m:ci>
		<m:apply>
		  <m:divide/>
		  <m:apply>
		    <m:times/>
		    <m:cn>1.998</m:cn>
		    <m:ci>k</m:ci>
		  </m:apply>
		  <m:apply>
		    <m:minus/>
		    <m:cn>17</m:cn>
		    <m:cn>0.999</m:cn>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>, 
	    <m:math>
	      <m:apply>
		<m:in/>
		<m:ci>k</m:ci>
		<m:set>
		  <m:cn>0</m:cn>
		  <m:ci>…</m:ci>
		  <m:cn>16</m:cn>
		</m:set>
	      </m:apply>
	    </m:math>.
	  </caption>
	</figure>
	<para id="par3222">
	  It is easy to verify that, at frequencies close to dc, the phase response of
	  <cnxn target="eq15"/> takes the approximate form
	  <equation id="eq17">
	    <m:math display="block">
	      <m:apply>
		<m:approx/>
		<m:apply>
		  <m:arg/>
		  <m:apply>
		    <m:ci type="fn">H</m:ci>
		    <m:ci>ω</m:ci>
		  </m:apply>
		</m:apply>
		<m:apply>
		  <m:approx/>
		  <m:apply>
		    <m:plus/>
		    <m:apply>
		      <m:minus/>
		      <m:apply>
			<m:divide/>
			<m:apply>
			  <m:sin/>
			  <m:ci>ω</m:ci>
			</m:apply>
			<m:apply>
			  <m:plus/>
			  <m:ci>c</m:ci>
			  <m:apply>
			    <m:cos/>
			    <m:ci>ω</m:ci>
			  </m:apply>
			</m:apply>
		      </m:apply>
		    </m:apply>
		    <m:apply>
		      <m:divide/>
		      <m:apply>
			<m:times/>
			<m:ci>c</m:ci>
			<m:apply>
			  <m:sin/>
			  <m:ci>ω</m:ci>
			</m:apply>
		      </m:apply>
		      <m:apply>
			<m:plus/>
			<m:cn>1</m:cn>
			<m:apply>
			  <m:times/>
			  <m:ci>c</m:ci>
			  <m:apply>
			    <m:cos/>
			    <m:ci>ω</m:ci>
			  </m:apply>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		  <m:apply>
		    <m:minus/>
		    <m:apply>
		      <m:times/>
		      <m:ci>ω</m:ci>
		      <m:apply>
			<m:divide/>
			<m:apply>
			  <m:minus/>
			  <m:cn>1</m:cn>
			  <m:ci>c</m:ci>
			</m:apply>
			<m:apply>
			  <m:plus/>
			  <m:cn>1</m:cn>
			  <m:ci>c</m:ci>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	  where the first approximation is obtained by replacing the argument of the
	  <m:math>
	    <m:apply>
	      <m:arctan/>
	    </m:apply>
	  </m:math>
	  with the function value and the second approximation, valid in an even
	  smaller neighborhood, is obtained by approximating
	  <m:math>
	    <m:apply>
	      <m:sin/>
	      <m:ci>x</m:ci>
	    </m:apply>
	  </m:math>
	  with
	  <m:math>
	    <m:ci>x</m:ci>
	  </m:math>
	  and
	  <m:math>
	    <m:apply>
	      <m:cos/>
	      <m:ci>x</m:ci>
	    </m:apply>
	  </m:math>
	  with
	  <m:math>
	    <m:cn>1</m:cn>
	  </m:math>. The phase and group delay around dc are
	  <equation id="eq18">
	    <m:math display="block">
	      <m:apply>
		<m:approx/>
		<m:apply>
		  <m:ci type="fn"><m:msub>
		      <m:mi>τ</m:mi>
		      <m:mi>ph</m:mi>
		    </m:msub></m:ci>
		  <m:ci>ω</m:ci>
		</m:apply>
		<m:apply>
		  <m:approx/>
		  <m:apply>
		    <m:ci type="fn"><m:msub>
			<m:mi>τ</m:mi>
			<m:mi>gr</m:mi>
		      </m:msub></m:ci>
		    <m:ci>ω</m:ci>
		  </m:apply>
		  <m:apply>
		    <m:divide/>
		    <m:apply>
		      <m:minus/>
		      <m:cn>1</m:cn>
		      <m:ci>c</m:ci>
		    </m:apply>
		    <m:apply>
		      <m:plus/>
		      <m:cn>1</m:cn>
		      <m:ci>c</m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	</para>
	<para id="par3224">
	  Therefore, the filter coefficient
	  <m:math>
	    <m:ci>c</m:ci>
	  </m:math>
	  can be easily determined from the desired low-frequency delay as
	  <equation id="eq19">
	    <m:math display="block">
	      <m:apply>
		<m:eq/>
		<m:ci>c</m:ci>
		<m:apply>
		  <m:divide/>
		  <m:apply>
		    <m:minus/>
		    <m:cn>1</m:cn>
		    <m:apply>
		      <m:ci type="fn"><m:msub>
			  <m:mi>τ</m:mi>
			  <m:mi>ph</m:mi>
			</m:msub></m:ci>
		      <m:cn>0</m:cn>
		    </m:apply>
		  </m:apply>
		  <m:apply>
		    <m:plus/>
		    <m:cn>1</m:cn>
		    <m:apply>
		      <m:ci type="fn"><m:msub>
			  <m:mi>τ</m:mi>
			  <m:mi>ph</m:mi>
			</m:msub></m:ci>
		      <m:cn>0</m:cn>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>
	  </equation>
	</para>

	<para id="par3225">
	  <cnxn target="fig3"/> shows that the delay of the allpass filter is approximately constant
	  only in a narrow frequency range. We can reasonably assume that such range,
	  for positive values of
	  <m:math>
	    <m:mi>c</m:mi>
	  </m:math>
	  smaller than one, extends from
	  <m:math>
	    <m:cn>0</m:cn>
	  </m:math> to 
	  <m:math>
	    <m:apply>
	      <m:divide/>
	      <m:ci><m:msub>
		  <m:mi>F</m:mi>
		  <m:mi>s</m:mi>
		</m:msub></m:ci>
	      <m:cn>5</m:cn>
	    </m:apply>
	  </m:math>. With
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:ci><m:msub>
		  <m:mi>F</m:mi>
		  <m:mi>s</m:mi>
		</m:msub></m:ci>
	      <m:apply>
		<m:times/>
		<m:cn>50</m:cn>
		<m:ci>kHz</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	  we see that at
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:divide/>
		<m:ci><m:msub>
		    <m:mi>F</m:mi>
		    <m:mi>s</m:mi>
		  </m:msub></m:ci>
		<m:cn>5</m:cn>
	      </m:apply>
	      <m:apply>
		<m:times/>
		<m:cn>10</m:cn>
		<m:ci>kHz</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	  we have an error of about
	  <m:math>
	    <m:cn>0.05</m:cn>
	  </m:math>
	  samples. In a note at that frequency produced by a feedback delay line, such an
	  error produces a pitch deviation smaller than
	  <m:math>
	    <m:apply>
	      <m:times/>
	      <m:cn>1</m:cn>
	      <m:ci>%</m:ci>
	    </m:apply>
	  </m:math>.  
	  For lower fundamental frequencies, such as those found in actual musical instruments, the error is smaller
	  than the just noticeable difference measured with slow pitch modulations. If the first-order filter represents an elegant and efficient solution to the
	  problem of tuning a delay line, it has also the relevant side effect of detuning
	  the upper partials, due to the marked phase nonlinearity. Such detuning can be
	  tolerated in most cases, but has to be taken into account in some other contexts.
	  If a phase response closer to linear is needed, we can use higher-order allpass filters
	  <cite src="#bib1"/>.
	  In some cases, especially in sound synthesis by physical modeling,
	  a specific inharmonic distribution of resonances has to be approximated. This
	  can be obtained by designing allpass filters that approximate a given phase
	  response along the whole frequency axis. In these cases the problem of tuning
	  is superseded by the more difficult problem of accurate partial positioning
	  <cite src="#bib2"/>.
	  With allpass interpolators it is more complicated to handle continuous delay
	  length variations, since the recursive structure of the filter does not show an obvious
	  way of transferring memory cells from and to the delay line, as it was in
	  the case of the FIR interpolator, which is constructed on the delay line by a
	  certain number of taps. Indeed, the glissando can be implemented with the allpass
	  filter by adding a new cell to the delay line whenever the filter coefficient
	  becomes one and, at the same time, zeroing out the filter state variable and the
	  coefficient. What is really more complicated with allpass filters is to handle
	  sudden variations of the delay length, as they are found, for instance, when a
	  finger hole is opened in a wind instrument. In this case, the recursive nature of
	  allpass filters causes annoying transients in the output signal. Ad hoc structures
	  have been devised to cancel these transients
	  <cite src="#bib1"/>.
	</para>
      </section>
    </section>

    <section id="par33">
      <name>The Non-Recursive Comb Filter</name>
      <para id="par3300">
	Sounds, propagating in the air, come into contact with surfaces and objects
	of various kinds and this interaction produces physical phenomena such as reflection, refraction, and diffraction. A simple and very important phenomenon
	is the reflection of sound about a planar surface. Due to a reflection such as
	this, a listener receives two delayed copies of the same signal. If the delay is
	larger than about a hundred milliseconds, the second copy is perceived as a
	distinguished echo, while if the delay is smaller than about ten milliseconds,
	the effect of a single reflection is perceived as a spectral coloration.
	A simple model of single reflection can be constructed starting from the
	basic blocks described in this and in the preceding chapters. It is constructed as an
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>-samples delay line, with the incidental fractional part of
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	obtained by FIR interpolation or allpass filtering, cascaded with an attenuation coefficient
	<m:math>
	  <m:ci>g</m:ci>
	</m:math>, possibly replaced by a filter if a frequency-dependent absorption has to be
	simulated. The output of this lossy delay line is summed to the direct signal.
	Let us analyze the structure in the case that
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	is integer and						
	<m:math>
	  <m:ci>g</m:ci>
	</m:math>
	is a positive constant not exceeding
	<m:math>
	  <m:cn>1</m:cn>
	</m:math>. The difference equation is expressed as
	<equation id="eq20">
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:ci type="fn">y</m:ci>
		<m:ci>n</m:ci>
	      </m:apply>
	      <m:apply>
		<m:plus/>
		<m:apply>
		  <m:ci type="fn">x</m:ci>
		  <m:ci>n</m:ci>
		</m:apply>
		<m:apply>
		  <m:times/>
		  <m:ci>g</m:ci>
		  <m:apply>
		    <m:ci type="fn">x</m:ci>
		    <m:apply>
		      <m:minus/>
		      <m:ci>n</m:ci>
		      <m:ci>m</m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
	and, therefore, the transfer function is
	<equation id="eq21">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:ci type="fn">H</m:ci>
		<m:ci>z</m:ci>
	      </m:apply>
	      <m:apply>
		<m:plus/>
		<m:cn>1</m:cn>
		<m:apply>
		  <m:times/>
		  <m:ci>g</m:ci>
		  <m:apply>
		    <m:power/>
		    <m:ci>z</m:ci>
		    <m:apply>
		      <m:minus/>
		      <m:ci>m</m:ci>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3302">
	In the case that
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>, it is easy to see by using the De Moivre formula
	that the frequency response of the comb filter has the following magnitude and group delay:
	<equation id="eq22">
	  <m:math display="block">
	    <m:mtable columnalign="left">
	      <m:mtr>
		<m:mtd>
		  <m:apply>
		    <m:eq/>
		    <m:apply>
		      <m:abs/>
		      <m:apply>
			<m:ci type="fn">H</m:ci>
			<m:ci>ω</m:ci>
		      </m:apply>
		    </m:apply>
		    <m:apply>
		      <m:root/>
		      <m:apply>
			<m:times/>
			<m:cn>2</m:cn>
			<m:apply>
			  <m:plus/>
			  <m:cn>1</m:cn>
			  <m:apply>
			    <m:cos/>
			    <m:apply>
			      <m:times/>
			      <m:ci>ω</m:ci>
			      <m:ci>m</m:ci>
			    </m:apply>
			  </m:apply>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:mtd>
	      </m:mtr>
	      <m:mtr>
		<m:mtd>
		  <m:apply>
		    <m:eq/>
		    <m:apply>
		      <m:ci type="fn"><m:msub>
			  <m:ci>τ</m:ci>
			  <m:mrow>
			    <m:mi>gr</m:mi>
			    <m:mo>,</m:mo>
			    <m:mi>H</m:mi>
			  </m:mrow>
			</m:msub></m:ci>
		      <m:ci>ω</m:ci>
		    </m:apply>
		    <m:apply>
		      <m:divide/>
		      <m:ci>m</m:ci>
		      <m:cn>2</m:cn>
		    </m:apply>
		  </m:apply>
		</m:mtd>
	      </m:mtr>
	    </m:mtable>
	  </m:math>
	</equation>
	and it is straightforward to verify that the frequency band ranging from dc to the Nyquist rate comprises
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	zeros (antiresonances), equally spaced by
	<m:math>
	  <m:apply>
	    <m:times/>
	    <m:apply>
	      <m:divide/>
	      <m:ci><m:msub>
		  <m:mi>F</m:mi>
		  <m:mi>s</m:mi>
		</m:msub></m:ci>
	      <m:ci>m</m:ci>
	    </m:apply>
	    <m:ci>Hz</m:ci>
	  </m:apply>
	</m:math>
	The phase response is piecewise linear with discontinuities of
	<m:math>
	  <m:pi/>
	</m:math>
	at the odd multiples of
	<m:math>
	  <m:apply>
	    <m:divide/>
	    <m:ci><m:msub>
		<m:mi>F</m:mi>
		<m:mi>s</m:mi>
	      </m:msub></m:ci>
	    <m:apply>
	      <m:times/>
	      <m:cn>2</m:cn>
	      <m:ci>m</m:ci>
	    </m:apply>
	  </m:apply>
	</m:math>. If
	<m:math>
	  <m:apply>
	    <m:lt/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>, it is easy to see that the amplitude of the resonances is
	<equation id="eq23">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:ci>P</m:ci>
	      <m:apply>
	      <m:plus/>
		<m:cn>1</m:cn>
		<m:ci>g</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
	while the amplitude of the points of minimum (halfway between contiguous resonances) is
	<equation id="eq24">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:ci>V</m:ci>
	      <m:apply>
		<m:minus/>
		<m:cn>1</m:cn>
		<m:ci>g</m:ci>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3305">
	An important parameter of this filtering structure, called <term>non-recursive comb filter</term> (or FIR comb), is the peak-to-valley ratio
	<equation id="eq25">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:divide/>
		<m:ci>P</m:ci>
		<m:ci>V</m:ci>
	      </m:apply>
	      <m:apply>
		<m:divide/>
		<m:apply>
		  <m:plus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
		<m:apply>
		  <m:minus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3306">
	<cnxn target="fig4"/>
	shows the response of a non-recursive comb filter having length
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>m</m:ci>
	    <m:cn>11</m:cn>
	  </m:apply>
	</m:math>
	samples and a reflection attenuation
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>g</m:ci>
	    <m:cn>0.9</m:cn>
	  </m:apply>
	</m:math>.
	The shape of the frequency response justifies the name comb given to the filter.
      <figure id="fig4">
	<media type="image/gif" src="combfirresp.gif"/>
	<caption>Magnitude of the frequency response of the comb FIR filter having coefficient
	    <m:math>
	      <m:apply>
		<m:eq/>
		<m:ci>g</m:ci>
		<m:cn>0.9</m:cn>
	      </m:apply>
	    </m:math>
	    and delay length
	    <m:math>
	      <m:apply>
		<m:eq/>
		<m:ci>m</m:ci>
		<m:cn>11</m:cn>
	      </m:apply>
	    </m:math>
	  </caption>
	</figure>
      </para>
      <para id="par3307">
	The zeros of the comb filter are evenly distributed along the unit circle at the			
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>-th roots of
	<m:math>
	  <m:apply>
	    <m:minus/>
	    <m:ci>g</m:ci>
	  </m:apply>
	</m:math>
	as shown in <cnxn target="fig5"/>.
	<figure id="fig5">
	  <media type="image/gif" src="pzcombfir.gif"/>
	  <caption>Zeros and poles of an FIR comb filter</caption>
	</figure>
      </para>
    </section>
    <section id="par34">
      <name>The Recursive Comb Filter</name>
      <para id="par3400">
	A simple model of one-dimensional resonator can be constructed using the
	basic blocks presented in this and in the preceding chapters. It is composed by an
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>-samples delay line, with the incidental fractional part of
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	obtained by FIR interpolation or allpass filtering, in feedback loop with an attenuation coefficient
	<m:math>
	  <m:ci>g</m:ci>
	</m:math>
	possibly replaced by a filter in order to give different decay times
	at different frequencies. Let us analyze the whole filtering structure in the case that
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	is integer and
	<m:math>
	  <m:ci>g</m:ci>
	</m:math>
	is a positive constant not exceeding
	<m:math>
	  <m:cn>1</m:cn>
	</m:math>.
	The difference equation is expressed as
	<equation id="eq26">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:ci type="fn">y</m:ci>
		<m:ci>n</m:ci>
	      </m:apply>
	      <m:apply>
		<m:plus/>
		<m:apply>
		  <m:times/>
		  <m:ci>x</m:ci>
		  <m:apply>
		    <m:minus/>
		    <m:ci>n</m:ci>
		    <m:ci>m</m:ci>
		  </m:apply>
		</m:apply>
		<m:apply>
		  <m:times/>
		  <m:ci>g</m:ci>
		  <m:ci>y</m:ci>
		  <m:apply>
		    <m:minus/>
		    <m:ci>n</m:ci>
		    <m:ci>m</m:ci>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
	and the transfer function is
	<equation id="eq27">
	  <m:math>
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:ci type="fn">H</m:ci>
		<m:ci>z</m:ci>
	      </m:apply>
	      <m:apply>
		<m:divide/>
		<m:apply>
		  <m:power/>
		  <m:ci>z</m:ci>
		  <m:apply>
		    <m:minus/>
		    <m:ci>m</m:ci>
		  </m:apply>
		</m:apply>
		<m:apply>
		  <m:minus/>
		  <m:cn>1</m:cn>
		  <m:apply>
		    <m:times/>
		    <m:ci>g</m:ci>
		    <m:apply>
		      <m:power/>
		      <m:ci>z</m:ci>
		      <m:apply>
			<m:minus/>
			<m:ci>m</m:ci>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3402">
	Whenever
	<m:math>
	  <m:apply>
	    <m:lt/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>, the stability is ensured. In the case that
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>, the frequency response of the filter has the following magnitude and group delay:
	<equation id="eq28">
	  <m:math display="block">
	    <m:mtable columnalign="left">
	      <m:mtr>
		<m:mtd>
		  <m:apply>
		    <m:eq/>
		    <m:apply>
		      <m:abs/>
		      <m:apply>
			<m:ci type="fn">H</m:ci>
			<m:ci>ω</m:ci>
		      </m:apply>
		    </m:apply>
		    <m:apply>
		      <m:divide/>
		      <m:cn>1</m:cn>
		      <m:apply>
			<m:times/>
			<m:cn>2</m:cn>
			<m:apply>
			  <m:sin/>
			  <m:apply>
			    <m:divide/>
			    <m:apply>
			      <m:times/>
			      <m:ci>ω</m:ci>
			      <m:ci>m</m:ci>
			    </m:apply>
			    <m:cn>2</m:cn>
			  </m:apply>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:mtd>
	      </m:mtr>
	      <m:mtr>
		<m:mtd>
		  <m:apply>
		    <m:eq/>
		    <m:apply>
		      <m:ci type="fn"><m:msub>
			  <m:mi>τ</m:mi>
			  <m:mrow>
			    <m:mi>gr</m:mi>
			    <m:mo>,</m:mo>
			    <m:mi>H</m:mi>
			  </m:mrow>
			</m:msub></m:ci>
		      <m:ci>ω</m:ci>
		    </m:apply>
		    <m:apply>
		      <m:divide/>
		      <m:ci>m</m:ci>
		      <m:cn>2</m:cn>
		    </m:apply>
		  </m:apply>
		</m:mtd>
	      </m:mtr>
	    </m:mtable>
	  </m:math>
	</equation>
	and it is easy to verify that the frequency band ranging from dc to the Nyquist rate comprises
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>
	vertical asymptotes (resonances), equally spaced by
	<m:math>
	  <m:apply>
	    <m:times/>
	    <m:apply>
	      <m:divide/>
	      <m:ci><m:msub>
		  <m:mi>F</m:mi>
		  <m:mi>s</m:mi>
		</m:msub></m:ci>
	      <m:ci>m</m:ci>
	    </m:apply>
	    <m:ci>Hz</m:ci>
	  </m:apply>
	</m:math>. If
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>
	the filter is at the limit of stability, and this is the only case when the phase response is piecewise linear, starting with the value
	<m:math>
	  <m:apply>
	    <m:minus/>
	    <m:apply>
	      <m:divide/>
	      <m:pi/>
	      <m:cn>2</m:cn>
	    </m:apply>
	  </m:apply>
	</m:math>
	at dc, with discontinuities of
	<m:math>
	  <m:pi/>
	</m:math>
	at the even multiples of
	<m:math>
	  <m:apply>
	    <m:divide/>
	    <m:ci><m:msub>
		<m:mi>F</m:mi>
		<m:mi>s</m:mi>
	      </m:msub></m:ci>
	    <m:apply>
	      <m:times/>
	      <m:cn>2</m:cn>
	      <m:ci>m</m:ci>
	    </m:apply>
	  </m:apply>
	</m:math>. If
	<m:math>
	  <m:apply>
	    <m:lt/>
	    <m:ci>g</m:ci>
	    <m:cn>1</m:cn>
	  </m:apply>
	</m:math>, it is easy to verify that the amplitude of the resonances is
	<equation id="eq29">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:ci>P</m:ci>
	      <m:apply>
		<m:divide/>
		<m:cn>1</m:cn>
		<m:apply>
		  <m:minus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
	while the amplitude of the points of minimum (halfway between contiguous resonances) is
	<equation id="eq30">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:ci>V</m:ci>
	      <m:apply>
		<m:divide/>
		<m:cn>1</m:cn>
		<m:apply>
		  <m:plus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3405">
	An important parameter of this filtering structure, called <term>recursive comb filter</term> (or IIR comb), is the peak-to-valley ratio
	<equation id="eq31">
	  <m:math display="block">
	    <m:apply>
	      <m:eq/>
	      <m:apply>
		<m:divide/>
		<m:ci>P</m:ci>
		<m:ci>V</m:ci>
	      </m:apply>
	      <m:apply>
		<m:divide/>
		<m:apply>
		  <m:plus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
		<m:apply>
		  <m:minus/>
		  <m:cn>1</m:cn>
		  <m:ci>g</m:ci>
		</m:apply>
	      </m:apply>
	    </m:apply>
	  </m:math>
	</equation>
      </para>
      <para id="par3406">
	<cnxn target="fig6"/>
	shows the frequency response of a recursive comb filter having a delay line of
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>m</m:ci>
	    <m:cn>11</m:cn>
	  </m:apply>
	</m:math>
	samples and feedback attenuation
	<m:math>
	  <m:apply>
	    <m:eq/>
	    <m:ci>g</m:ci>
	    <m:cn>0.9</m:cn>
	  </m:apply>
	</m:math>
	The shape of the magnitude response justifies the name comb given to the filter.
	<figure id="fig6">
	  <media type="image/gif" src="combresp.gif"/>
	  <caption>Magnitude and phase delay response of the recursive comb filter having coefficient
	    <m:math>
	      <m:apply>
		<m:eq/>
		<m:ci>g</m:ci>
		<m:cn>0.9</m:cn>
	      </m:apply>
	    </m:math>
	    and delay length
	    <m:math>
	      <m:apply>
		<m:eq/>
		<m:ci>m</m:ci>
		<m:cn>11</m:cn>
	      </m:apply>
	    </m:math>
	  </caption>
	</figure>
      </para>
      <para id="par3407">
	The poles of the comb filter are evenly distributed along the unit circle at the
	<m:math>
	  <m:ci>m</m:ci>
	</m:math>-th roots of 
	<m:math>
	  <m:ci>g</m:ci>
	</m:math>, as shown in <cnxn target="fig7"/>.
	<figure id="fig7">
	  <media type="image/gif" src="pzcombiir.gif"/>
	  <caption>Zeros and poles of an IIR comb filter</caption>
	</figure>
      </para>
      <para id="par3408">
	In sound synthesis by physical modeling, a recursive comb filter can be interpreted
	as a simple model of lossy one-dimensional resonator, like a string,
	or a tube. This model can be used to simulate several instruments whose resonator
	is not persistently excited. In fact, if the input is a short burst of filtered
	noise, we obtain the basic structure of the plucked string synthesis algorithm
	due to Karplus and Strong <cite src="#bib3"/>.
      </para>

      <section id="par341">
	<name>The Comb-Allpass Filter</name>
	<para id="par3410">
	  The filter given by the difference  <cnxn target="eq20"/> 
	  has a frequency response characterized by evenly-distributed resonances.
	</para>
	<para id="par3411">
	  With a slight modification of its structure, such filter can be made allpass. In other words, the magnitude response
	  of the filter can be made flat even though the impulse response remains almost the same
	  <cnxn target="eq20"/>.
	  The modification is just a direct path connecting the input
	  of the delay line to the filter output, as it is depicted in <cnxn target="fig8"/>.
	  <figure id="fig8">
	    <media type="image/gif" src="combap.gif"/>
	    <caption>Allpass comb filter</caption>
	  </figure>
	</para>
	<para id="par3412">It is easy to see that the transfer function of the filter of
	  <cnxn target="fig8"/>, called the allpass comb filter can be written as
	  <equation id="eq33">
	    <m:math display="block">
	      <m:apply>
		<m:eq/>
		<m:apply>
		  <m:ci type="fn">H</m:ci>
		  <m:ci>z</m:ci>
		</m:apply>
		<m:apply>
		  <m:divide/>
		  <m:apply>
		    <m:plus/>
		    <m:apply>
		      <m:minus/>
		      <m:ci>g</m:ci>
		    </m:apply>
		    <m:apply>
		      <m:power/>
		      <m:ci>z</m:ci>
		      <m:apply>
			<m:minus/>
			<m:ci>m</m:ci>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		  <m:apply>
		    <m:minus/>
		    <m:cn>1</m:cn>
		    <m:apply>
		      <m:times/>
		      <m:ci>g</m:ci>
		      <m:apply>
			<m:power/>
			<m:ci>z</m:ci>
			<m:apply>
			  <m:minus/>
			  <m:ci>m</m:ci>
			</m:apply>
		      </m:apply>
		    </m:apply>
		  </m:apply>
		</m:apply>
	      </m:apply>
	    </m:math>  
	  </equation>
	  which has the structure of an allpass filter. It is interesting to note that the
	  direct path introduces a nonzero sample at the time instant zero in the impulse
	  response. All the following samples are just a scaled version of those of the
	  impulse response of the comb filter, with a scaling factor equal to
	  <m:math>
	    <m:apply>
	      <m:minus/>
	      <m:cn>1</m:cn>
	      <m:apply>
		<m:power/>
		<m:ci>g</m:ci>
		<m:cn>2</m:cn>
	      </m:apply>
	    </m:apply>
	  </m:math>.
	  The time properties, such as the time decay, are substantially unvaried. The allpass
	  comb filter does not introduce any coloration in stationary signals. On the other
	  hand, its effect is evident on signals exhibiting rapid transients, and for these
	  signals we can not state that the filter is transparent.
	</para>
      </section>
    </section>

    <section id="par35">
      <name>Sound Effects Based on Delay Lines</name>
      <para id="par350">Many of the effects commonly used in electroacoustic music are obtained
	by composition of time-varying delay lines, i.e., by lines whose length is modulated
	by slowly-varying signals. In order to avoid discontinuities in the signals,
	it is necessary to interpolate the delay lines in some way. The interpolation
	by means of allpass filters is applicable only for very slow modulations or
	for narrow-width modulations, since sudden changes in the state of allpass filters
	give rise to transients that can be perceived as signal distortions <cite src="#bib4"/>. On
	the other hand, linear (or, more generally, polynomial) interpolation introduces
	frequency-dependent losses whose magnitude is dependent on the fractional
	length of the delay line. As the delay length is varied, these variable losses give
	an amplitude distortion due to amplitude modulation of the various frequency
	components. Coupled to amplitude modulation, there is also phase modulation
	due to phase nonlinearity of the interpolator, in both cases of FIR and IIR interpolation.
	The terminology used for audio effects is not consistent, as terms such as
	flanger, chorus, and phaser are often associated with a large variety of effects,
	that can be quite different from each other. A flanger is usually defined as an
	FIR comb filter whose delay length is sinusoidally modulated between a minimum
	and a maximum value. This has the effect of expanding and contracting
	the harmonic series of notches of the frequency response. The name flanger
	derives from the old practice, used long ago in the analog recording studios,
	to alternatively slow down the speed of two tape recorders or two turntables
	playing the same music track by pressing a finger on the flanges.
	The name phaser is most often reserved for structures similar to the comb
	FIR filter, with the difference that the notches are not harmonically distributed.
	Orfanidis <cite src="#bib5"/> proposes to use, instead of the delay line, a bunch of parametric
	notch filters. Each notch is controllable in its frequency position and width. Smith
	<cite src="#bib6"/>, instead, proposes to use
	a large allpass filter instead of the delay line. If this allpass filter is obtained as
	a cascade of second-order allpass sections, it becomes possible to control and
	modulate the position of any single pole couple, which represent all the single
	notches of the overall response. A common feature of flangers and phasers is
	the relatively large distance between the notches. Vice versa, if the notches are
	very dense, the term chorus is preferred. Orfanidis <cite src="#bib5"/>, suggests to implement
	a chorus as a parallel of FIR comb filters, where the delay lengths are randomly
	modulated around values that are slightly different from each other. This should
	simulate the deviations in time and height that are found in performances of a
	choir singing in unison. Vice versa, Dattorro <cite src="#bib4"/> says that a chorus can be obtained
	by the same structure used for the flanger, with a difference that the delay
	lengths have to be set to larger values than for the flanger. In this way, the
	notches are made more dense. For the flanger the suggested nominal delay is
	1msec and for the chorus it is 5msec. If the objective is to recreate the effect of
	a choir singing in unison, the fact of having many notches in the spectrum is
	generally disliked. Dattorro <cite src="#bib4"/> proposes a partial solution that makes use of a
	recursive allpass filter, where the delay line is read by two pointers, one is kept
	fixed and produces the feedback signal, the other is varied to pick up the signal
	that is fed directly to the output. In this way, when both the pointers are at the
	nominal position, the structure does not introduce any coloration for stationary
	signals. A final remark is reserved to the spatialization of these comb-based effects.
	In general, flanging, phasing, and chorusing effects can be obtained from two
	different time-varying allpass chains, whose outputs feed different loudspeakers.
	In this case, sums and subtractions between signals at the different frequencies
	happen “on air” in a way dependent from position. Therefore, the spatial
	sensation is largely due to the different spectral coloration found in different
	points of the listening area.
      </para>
    </section>
  </content>

  <bib:file>
    <bib:entry id="bib1">
      <bib:article>
	<bib:author>Timo I. Laakso, Vesa V\"alim\"aki, Matti Karjalainen, and Unto K. Laine</bib:author>
	<bib:title>Splitting  the Unit Delay - Tools for Fractional Delay Filter Design</bib:title>
	<bib:journal>IEEE Signal Processing  Magazine</bib:journal>
	<bib:year>1999</bib:year>
	<bib:month>Sep</bib:month>
      </bib:article>
    </bib:entry>
    <bib:entry id="bib2">
      <bib:article>
	<bib:author>Davide Rocchesso and Francesco Scalcon</bib:author>
	<bib:title>Bandwidth of perceived inharmonicity for physical modeling of dispersive strings</bib:title>
	<bib:journal>IEEE Transactions on Speech and Audio Processing</bib:journal>
	<bib:year>1996</bib:year>
	<bib:month>Jan</bib:month>
      </bib:article>
    </bib:entry>
    <bib:entry id="bib3">
      <bib:article>
	<bib:author>K. Karplus and A. Strong</bib:author>
	<bib:title>Digital Synthesis of Plucked String and Drum Timbres</bib:title>
	<bib:journal>Computer Music J.</bib:journal>
	<bib:year>1983</bib:year>
      </bib:article>
    </bib:entry>
    <bib:entry id="bib4">
      <bib:article>
	<bib:author>Jon Dattorro</bib:author>
	<bib:title>Effect design - part 2: Delay-line modulation and chorus</bib:title>
	<bib:journal>J. Audio Eng. Soc.</bib:journal>
	<bib:year>1997</bib:year>
	<bib:month>Oct</bib:month>
      </bib:article>
    </bib:entry>
    <bib:entry id="bib5">
      <bib:book>
	<bib:author>S. J. Orfanidis</bib:author>
	<bib:title>Introduction to Signal Processing</bib:title>
	<bib:publisher>Prentice Hall</bib:publisher>
	<bib:year>1996</bib:year>
      </bib:book>
    </bib:entry>
    <bib:entry id="bib6">
      <bib:article>
	<bib:author>Julius O. Smith</bib:author>
	<bib:title>An allpass approach to digital phasing and flanging</bib:title>
	<bib:journal>In Proc. International Computer Music Conference - Also available as Rep. STAN-M-21, CCRMA, Stanford University</bib:journal>
	<bib:year>1984</bib:year>
      </bib:article>
    </bib:entry>
  </bib:file>
</document>
