Modal Decomposition of Feedback Delay Networks
Abstract
Feedback delay networks (FDNs) belong to a general class of recursive filters which are widely used in sound synthesis and physical modeling applications. We present a numerical technique to compute the modal decomposition of the FDN transfer function. The proposed pole finding algorithm is based on the EhrlichAberth iteration for matrix polynomials and has improved computational performance of up to three orders of magnitude compared to a scalar polynomial root finder. We demonstrate how explicit knowledge of the FDN’s modal behavior facilitates analysis and improvements for artificial reverberation. The statistical distribution of mode frequency and residue magnitudes demonstrate that relatively few modes contribute a large portion of impulse response energy.
¿\SplitArgument1;m \setargsaux#1 \NewDocumentCommand\setargsauxmm \IfNoValueTF#2#1 #1 \delimsize— #2
I Introduction
A feedback delay network (FDN) consists of a set of delay lines with lengths which are interconnected via a feedback matrix (see Fig. 1). FDNs arise in many physical modeling applications where geometrically distributed components are approximated by time delays, e.g., strings [1], plates and membranes [2], springs [3], and air volume [4, 5]. The interest in FDNs is fueled by the highly efficient implementation of delays in the timedomain, e.g., with circular buffers resulting in a constant time complexity independent of its length. Therefore, the computational complexity of the FDN scales with the number of delay lines and not the system order. FDNs are a popular choice for artificial reverberation applications particularly because of the favorable relation between FDN size and system order [6, 7, 8, 9].
In this work, we present a modal decomposition technique for FDNs. The modal decomposition of a system is an equivalent representation as the sum of complex onepole resonators, socalled modes. The timedomain signal of such a resonator with pole and residue is
(1) 
where indicates the argument of a complex number in radiant, is the magnitude, and indicates the discrete time index. Each individual resonating mode is governed by four parameters: mode frequency , decay rate , initial phase and initial amplitude (see Fig 1). With modal decomposition, we aim to uncover the specific parameters of each mode. The timedomain impulse response of the FDN
(2) 
is the sum of the complex modes , where is the system order. In sound synthesis applications for instance, the human auditory system can recognize the spectral quality composed of the individual modes and this representation is therefore termed additive or modal synthesis [10]. Modal analysis of recursive systems is applied in various system modeling applications, ranging from acoustics and digital filter design to mechanical modeling [11]. A particularly challenging application for modal decomposition is room acoustics, where even medium room sizes exhibits millions of modes [12]. Only for simple room geometries, an analytic expression for the system poles and residues can be stated [13]. System poles may also be recovered from the impulse response by various techniques such as an autoregressive movingaverage [14], Bayesian inference [15] and allpole modeling [16]. Whereas these techniques may be able to successfully compute partial solutions or compute the solution for specific configurations, the computation of the entire set of modes is in general challenging. In the following, we give the precise problem statement of this work.
Ia Problem Statement
For a single input and single output, the timedomain recursion of an FDN with delay lines is given by
(3)  
where is the time index, denotes the transpose operation and , , [17]. The state vector is defined as . We write FDN to denote an FDN of size . The transfer function of an FDN is
(4) 
where . The system order is given by [17]. For commonly used delays , the system order is much larger than the FDN size, i.e.,
(5) 
The modal decomposition of the FDN, i.e., the partial fraction decomposition (PFD) of the transfer function (4) is
(6) 
where is the residue of the pole . The timedomain representation of the sum in (6) is given in (2) as the sum of complex resonators. The objective of this work is to present an efficient numerical method to compute the modal decomposition (6) from the transfer function (4).
IB Direct Approach
We first review two standard methods for the modal decomposition [17, 18]. Let be any invertible matrix, then
(7) 
where is the adjugate of the matrix [19]. In the following, we denote
(8) 
With (7) and (8), the transfer function (4) can be expressed as a rational polynomial
(9) 
where
(10) 
and
(11) 
For brevity, we occasionally omit the parameters and write and . The FDN system poles , where , are the roots of the generalized characteristic polynomial (GCP) in (10) such that they are fully characterized by the delay matrix and the feedback matrix .
For a moment, let us assume that all delays are single time steps, i.e., . The timedomain recursion in (3) reduces to the standard statespace description of a linear timeinvariant (LTI) filter. The system poles are the eigenvalues of the feedback matrix such that the modal decomposition (6) is easily computed with standard methods. However, for longer delays such that (5) holds, the modal decomposition becomes more involved.
The GCP can be expressed in a linearized fashion
(12) 
where is the identity matrix of size and such that the system poles are the eigenvalues of [17]. Unfortunately, for large delays this eigenvalue problem becomes quickly numerically intractable. Alternatively, the GCP can be expressed as a scalar polynomial
(13) 
where the coefficients are derived from the principal minors of [18]. The system poles are the roots of the scalar polynomial. Again, the polynomial degree increases with longer delays and finding the roots of the polynomial becomes numerically intractable [20].
In the remainder of this paper, we present a numerically stable and computationally efficient method to compute the modal decomposition for large system order and modestsized . In Section II, we derive the fundamental algorithm based on a polynomial matrix formulation. In Section III, we evaluate the performance of the proposed algorithm. In Section IV, we apply modal decomposition to analyze the effects of attenuation filters and to study the statistical distributions of mode frequencies and residue magnitudes.
Ii Numerical Modal Decomposition
In the following, we present a root finding algorithm for the GCP and subsequently recover the residues . We conclude this section with a generalization to additional filtering in the delay lines and feedback matrix.
Iia Polynomial Matrix Formulation
It is a common heuristic in numerical computation that the inherent problem structure shall be preserved as much as possible throughout all computation steps to improve numerical performance. In contrast to Section IB, we compute the system poles without expanding the problem. In fact, (10) is a polynomial eigenvalue problem of degree , i.e.,
(14) 
where for . For a proper matrix polynomial , i.e., , the number of roots is [21]. For FDNs, however, is singular such that the actual number of roots is lower, respectively, many roots are infinite. In fact, if , the number of finite roots is which is also the degree of the scalar polynomial in (10) [18].
In the following, we use the derivative of the polynomial . According to Jacobi’s formula [22], we have
(15)  
where and denotes the trace of matrix . Stewart [23] showed that the adjugate of can be wellconditioned even when is illconditioned, and he shows how can be computed in a numerically stable way from a rank revealing decomposition of [22]. With the definitions of the polynomial eigenvalue problem introduced, we present the proposed root finding algorithm.
IiB EhrlichAberth Method
The polynomial eigenvalue problem can be solved with the EhrlichAberth Iteration (EAI) method, i.e., a combination of Newton method and a deflation term which prevents that two eigenvalues converge to the same solution [21]. Let be a vector of initial estimates for the roots of the polynomial and be the th EAI iteration. The EAI provides the sequence of estimates
(16) 
with the EAI step being
(17) 
Using the identity in (15), the Newton correction term is
(18) 
and the deflation term is
(19) 
The deflation term may be interpreted as a penality term if two eigenvalues approach each other too closely and guarantees that the all eigenvalues reached are unique. Using (18) we can expand (17) to
(20) 
The method, given here in the Jacobi version, is known to converge cubically for simple roots and linearly for multiple roots [21]. The GaussSeidel version of EAI [21], which updates the estimates as soon as they become available, may converge even slightly faster.
IiC Stopping Criteria
The system poles are the roots of the polynomial in (10), i.e., . In other words, is a singular matrix for all system poles. Thus, the natural stopping criteria is the reciprocal of the condition number being less than a prescribed tolerance . This stopping condition is also computationally favorable as the condition number can be estimated highly efficiently [22]. However, for multiple eigenvalues this stopping condition may result in a premature halt [21].
An alternative stopping condition says that the computed correction is too tiny and would not change the significant digits of the current estimate
(21) 
where is a small positive tolerance threshold [21]. In practice, good global convergence properties are observed; a theoretical analysis of global convergence, though, is still missing and constitutes an open problem. There is empirical evidence that the number of Newton iterations heavily depends on the choice of the initial estimates [21].
IiD Initialization
Aberth [24] proposed to choose initial estimates placed along a circle centered at the origin of sufficiently large radius so that it contains all the roots. In case the magnitude of the roots vary largely, multiple circles with suitable radii may be chosen instead [25]. With Rouché’s theorem, we can derive upper and lower bounds on the pole magnitudes for the FDN depending on the singular values of the feedback matrix (see Appendix A)
(22) 
Equation (22) is a generalization on the relation of eigenvalues and singular values as given in the WeylHorn Theorem [26] in the case of unit delays . The bound is tight for a diagonal feedback matrix where the minimum and maximum delays coincide with the minimum and maximum diagonal element, respectively. However, the bound may be arbitrarily loose. For instance, the maximum singular value of a triangular matrix may be arbitrarily large while all system poles lie on the unit circle [9]. For large delays however, (22) shows that the pole magnitudes tend to be close to the unit circle.
We can further derive from (22) that if then all poles lie on the closed unit disk which is equivalent to the FDN being marginally stable [27]. In particular, if all singular values are 1, which is equivalent to being unitary, i.e., , all system poles lie on the unit circle regardless of the delays . Such an FDN is called lossless, and represents an important special case [9].
In this work, we are interested in lossless and stable FDNs for their practical relevance. In combination with (5), such FDNs have all poles in the unit disk, but close to the unit circle. Thus, we place the initial estimates uniformly on the unit circle. More precisely, we chose the roots of unity
(23) 
It is worthwhile to note that is the solution of a particular FDN with a circular shift matrix
(24) 
such that the GCP is
(25) 
The shift matrix thus combines the FDN delays into a single long delay line.
IiE Approximate Deflation
For a high system order , the computational complexity of the deflation term (19) may become excessive. We propose an approximate deflation (AD) according to a maximum error tolerance for the resulting EAI step in (20), i.e.,
(26) 
The magnitude of the deflation term summands decreases with the pole distance . The idea is then to divide the poles into a near and far pole sets, and , respectively, and approximate the deflation of the less significant far poles by a default term, e.g., . For symmetry, the number of near poles is assumed to be an even number.
It can be shown, that for equidistributed poles such as in (23), the far deflation is (see Appendix B)
(27) 
Thus, the total deflation may be approximated by
(28) 
if the far poles are sufficiently uniformly distributed. By sorting the system poles iterations along pole angles, we can find the near poles . To establish the quality of this approximation, let us assume that there exists an upper bound for the approximation error of the deflation term, i.e.,
(29) 
We can show that the error tolerance in (26) is satisfied if (see Appendix C)
(30) 
In other words, if the deflation approximation is sufficiently far from the inverse Newton term, the deflation error becomes negligible. On the contrary, if the deflation term is close to the inverse Newton step, the EAI step error can be large even for small deviations in the deflation term. In case the error tolerance in (30) is not satisfied, we compute the exact deflation instead.
To implement the proposed approximate deflation, we need a priori knowledge of the approximation error bound . As depends on many factors such as matrix size , system order , feedback matrix and number of near poles , in this work, it is determined experimentally from random FDNs. The performance of the approximate deflation method may be quantified by the number of exact deflations versus approximate deflations. As the initial estimates are equidistributed, may be computed only from the estimated far deflation with .
IiF Residues
Once we have found the system poles, the residues of the modal decomposition (6) are computed by
(31) 
where we assume that all poles are unique. Similar, but more intricate solutions exists for nonunique poles [28]. The undriven residue, i.e., the system response without excitation, is
(32) 
The undriven residue is a valuable intermediate step to analyze the mode initial amplitude independent from the input and output drives . Since, is a singular matrix, the derivative of the GCP in (15) may only be computed by the adjugate formulation. Since , the inputoutput drives in (11) are
(33) 
The difference between the driven and undriven residues may be expressed as a linear combination of the matrix entries of . Alternatively, the driven residues may also be computed by a least linear squares fit for the timedomain impulse response since the sum of complex resonators in (2) depends linearly on the residues [29].
IiG Polynomial Feedback and Delay Matrices
Although, the focus of this work is on frequencyindependent feedback matrices , much of the development in Section II is applicable to general polynomial matrices. Therefore, it is easy to include further filtering such as a frequencydependent feedback matrix . There also exists a singular value decomposition for polynomial matrices [30]. Alternatively, the delay lines are often extended with an attenuation or allpass filter , i.e.,
(34) 
It is important to note that additional filters may increase the number of system poles. Further, if is a rational polynomial, in other words consists of IIR filters, then the transfer function in (4) is no longer proper, i.e., the polynomial degree of the nominator is larger than the polynomial degree of the denominator [31]. Nonetheless, improper partial fraction decomposition can be solved with a delayed parallel form by separating the FIR and IIR part of the transfer function [29].
Iii Modal Synthesis and Evaluation
The following evaluation uses realvalued FDN parameters such that the system poles appear in complex conjugate pairs.
Iiia Modal Synthesis and Accuracy
A numerically accurate way to verify the modal decomposition is to synthesize each mode in timedomain as expressed in (1) and compare the sum of all modes with the impulse response computed by the timedomain recursion in (3). The concept of modal synthesis and verification is depicted in Fig. 1. The error is given by the maximum difference^{1}^{1}1The maximum error is chosen as it is an upper bound for the root mean squared error (RMSE) and as such a strict error measure. between the two impulse responses, i.e.,
(36) 
In this work, we use double precision floating point arithmetic and the modal decomposition is regarded successful if the maximum error . The EAI is numerically stable as the matrix inversion in (20) is only necessary if the matrix is sufficiently nonsingular due to the first stopping criteria. For large delays , the evaluation of may become extremely large or small if is too far away from the unit circle. At the same time, the poles tend to be close to the unit circle for large delays due to the bounds given in (22). As a practical intervention, the pole location is clipped to the magnitude bounds if the EAI step causes the pole location to exceed the bounds.
IiiB Numerical Evaluation
For the FDN, a single EAI step in (20) can be evaluated in : an evaluation of and is merely an evaluation of the delay matrix in ; a numerical matrix inversion can be performed in ; and the deflation term is evaluated in . Thus, a full iteration from can be evaluated in . This compares favorably with the bound of a matrixbased algorithm applied to the linearization in (12). For a high number of system poles , the complexity of computing the deflation term in (19) becomes the dominating part. The complexity of the approximate deflation in (28) is similar asymptotically, however in practice, the computational complexity is reduced significantly.
Fig. 2 shows the average number of full iterations depending on the matrix size , total delays between 50 and samples and a random orthogonal feedback matrix. It can be seen that the number of iterations is largely dependent on the parity of the matrix size and . This illustrates how the initialization (23) influences the performance of the EAI. Overall, about 4 to 5 iterations per root may be expected for the EAI to converge.
Fig. 3 depcits a comparison of measured computation time with the MATLAB^{2}^{2}2Matlab is a registered trademark of The MathWorks Inc. All computations were performed with Matlab R2016b on a desktop machine with an Intel Core i7 @ 3,40 GHz and 32 GB of RAM. functions eig and roots solving the direct problems (12) and (13), respectively. The total number of delays were distributed randomly among eight delay lines and the feedback matrix was a random orthogonal matrix. All methods gave the correct answer with the required accuracy. For the approximate deflation, the number of near poles was set to . The maximum deflation error was determined a priori by probing an independent set of random FDNs of similar configuration. The EAI step tolerance was set to .
The EAI implementation utilizes only standard MATLAB functions and no Coptimization which explains relatively poor performance for small system order . For high system order such as , the standard EAI and EAI with AD outperform the MATLAB’s eig function by a factor of more than 300 and 1300, respectively. Further, the memory requirements of the EAI are only linear in and cubic in such that it is possible to perform modal decomposition up to . Whereas the memory requirements for eig become prohibitive for . For , more the of the computation time of the standard EAI was spent on the deflation term. For the EAI with AD, the number of exact iterations never exceeded 1% of the total number of iterations proofing the chosen heuristic parameters effective. The EAI with AD performs similar for small delays but outperforms the standard EAI by a factor of 100 for large delays. Each EAI step is independent and only requires synchronization at every full iteration step such that the overall performance of the EAI might further be improved by parallelization.
Iv Analysis of Feedback Delay Networks
We study two applications of modal decomposition in artificial reverberation: Firstly, we study the effect of attenuation filters on the poles and residues of an FDN. Secondly, we study the statistical distribution of poles and residues of random lossless FDNs.
Iva Attenutation
Attenuation filters in FDNs, as they are typically applied in artificial reverberation, aim to control the frequencydependent reverberation time [7, 32]. As expressed in (34), all delays are extended with absorption filters and the feedback matrix is orthogonal or more generally unilossless [9]. We study three types of attenuation: homogeneous, nearhomogeneous and inhomogeneous attenuation.
IvA1 Homogeneous Attenuation
The attenuation filters are called homogeneous if there exists an attenuationpersample such that
(37) 
The attenuated delay lines can be expressed as plain delay lines with a mapped argument, i.e.,
(38) 
Consequently, the system poles with attenuation can be related to the system poles without attenuation by
(39) 
If we assume that the attenuation filters have a purely real frequency response^{3}^{3}3Although such a frequency response is not realizable with a digital filter in general, but useful for the theoretical analysis., i.e., then the mode frequencies are unaltered by the attenuation, i.e., . For a unilossless , all unattenuated system poles lie on the unit circle such that
(40) 
and the attenuated FDN is stable if . For homogeneous attenuation, the magnitude bounds in (35) are tight.
IvA2 Nearhomogeneous Attenuation
Typically, the attenuation filters are implemented with relatively low order, such as onepole filters [7], however higher order filters were proposed as well [33, 32]. The attenuation filters are designed to match the magnitude response
(41) 
where the attenuationpersample is derived from a target reverberation time
(42) 
where is the sampling frequency and is the time in seconds for the energy decay curve of the impulse response at frequency to decay by 60 dB [34]. For illustration, we compute the onepole filter according to [7] given the target reverberation time at DC and Nyquist frequency . Figure 4 depicts the resulting modal decomposition for an 8FDN with an orthogonal feedback matrix and a target reverberation time seconds and seconds. The system pole magnitudes are modified according to the target reverberation time. However, the attenuation varies especially in the transition band due to errors in the magnitude response caused by the limited filter order. The magnitude of the residues are depicted in Fig. 4. For nearhomogeneous attenuation, it can be observed that the residues with and without attenuation are rather similar. Although the attenuation filters are not completely homogeneous, their phase component is small compared to the phase of the delays such that the overall behavior is well approximated by (40). This suggests that studies on residues of lossless systems may translate well to results for moderately lossy systems.
IvA3 Inhomogeneous Attenuation
While the homogeneous attenuation has perceptually desirable properties in artificial reverberation, more physically oriented FDN designs such as scattering delay networks [35] and radiance transfer [36] employ attenuation filters which are unrelated to the delay lengths but related to the boundary materials of the simulated space. Figure 5 depicts the modal decay rate of the same 8FDN as in Fig. 4 with different attenuation filters. Instead of the delay proportional design in (41), all onepole filters have the same target frequency response corresponding to an average delay length. As a consequence, the decay time of the neighboring modes are largely different, while the overall shape still follows the target reverberation time.
IvB Statistical Distribution of Poles and Residues
We present a set of statistical analyses of lossless FDNs which rely on the proposed largescale numerical computation of the modal decomposition and are difficult to derive by analytic methods. The statistical analysis answers a longstanding question in artificial reverberation design [37]: Why do some FDNs have an unpleasant metallic ringing despite a sufficiently high modal density? While ideal late reverberation has been characterized as Gaussian white noise [38], the metallic ringing is caused by excessive energy at few frequencies. In terms of modal decomposition, metallic ringing may be caused by either clustering of multiple poles at the ringing frequencies or largely varying energy of neighboring modes. We study the following two questions:

What is the distribution of the mode frequencies?

What is the distribution of residue magnitudes?
In the analyses, we rely on Monte Carlo simulations of randomly generated lossless FDNs.
IvB1 Mode Frequency Distribution
The nearequidistribution of mode frequencies has been conjectured before [39] and the authors have given an analytical bound on the equidistribution based on Hayman’s theorem [18]. The cluster number
(43) 
is a measure on how equally distributed the mode frequencies are. Here, denotes the cardinality of a set. The higher the cluster number, the more poles cluster around the frequency . In contrast, a mode gap occurs if , i.e., no mode lies in this frequency interval. For perfectly equidistributed poles for all . We evaluate the distribution of mode frequencies by computing the histogram of cluster numbers
(44) 
where is the dirac function, is the integer cluster size, and is the number of observations. For large enough , the histogram converges towards the probability of cluster numbers. The random 8FDNs have delays between 50 and 1000 samples and an orthogonal feedback matrix. The probabilities are averaged over 100 random instances each.
In Table I, the probability of cluster numbers for randomly generated FDNs are compared to cluster numbers of a pseudouniform random number generator with equal sampling size. The discrepancy of the cluster number from an equidistribution is relatively low for the FDN modes compared to the random number generator. In fact for FDNs, it is very rare to find an interval of width with more than two modes. In stark contrast, acoustic mode density of physical spaces increase quadratically with frequency [12].
Cluster size  0  1  2  3  4 

Uniform Random  0.3690  0.3661  0.1854  0.0610  0.0186 
Lossless 8FDN  0.1694  0.6632  0.1653  0.0020  0.0001 
Equidistributed  0.0000  1.0000  0.0000  0.0000  0.0000 
IvB2 Residue Magnitude Distribution
In Section IIF, we have presented the computation of the mode residues for a given set of system poles. Figure 6 depicts the magnitude histogram of the total and undriven residues as well as the inputoutput drives for a random 8FDN. The inputoutput drives are comprised of all individual inputoutput combinations, i.e., in (33). The total residues result from unit input and output gains, i.e., and , or in other words, . The magnitude distributions of the inverse undriven residues , the total residues and inputoutput drives all resemble logRayleigh distributions [40]. However, just by altering the feedback matrix , it is possible to encounter various other distributions of the residue magnitude. Figure 7, depicts the residue magnitude distribution of four selected orthogonal feedback matrices.
IvB3 Discussion
For randomly generated FDNs, the mode frequencies are nearly equidistributed such that every frequency band has energy contributions from a similar number of modes. On the other hand, the high dynamic range of the residue magnitudes suggest that a small number of poles contribute a large portion of the impulse response energy. In the context of artificial reverberation, the highenergy modes dominate the frequency spectrum such that the audible modal density is considerably lower than theoretic modal density, i.e., the number of modes per frequency. For illustration, we have synthesized audio examples from the four instances depcited in Fig. 7 and provided them online^{4}^{4}4www.audiolabserlangen.de/resources/2018IEEEModal.
The residue distribution may be optimized in two steps: Firstly, optimization of the undriven residues by choosing delays and feedback matrix . Secondly, optimization of the total residue by choosing the input and output gains, and . While the first step is a nonlinear process which requires further research, the second step may be readily solved by linear least square fitting.
V Conclusion
We presented a numerically efficient technique for modal decomposition of the FDN. Standard methods such as eigenvalue decomposition of the linearized system and polynomial root finding methods applied to the characteristic polynomial require significant computational resources when the system order is large. The proposed method applies the EhrlichAberth Iteration to the polynomial matrix formulation of the FDN. Further we proposed, an efficient approximate deflation technique based on the estimation of far poles. For high system order such as , the standard EAI and approximate EAI outperform the MATLAB’s eig function by a factor of more than 300 and 1300, respectively. The approximate EAI was able to give reliable results up to a system order of 1 million. The modal decomposition was applied to FDNs in the context of artificial reverberation. Three types of attenuation were studied: homogeneous, nearhomogeneous and inhomogeneous. The potential for explicit analysis of the pole and residues was demonstrated for attenuation filter design. Statistical analysis showed that for randomly generated FDNs, the mode frequencies are nearly equidistributed and the residue magnitudes follow a logRayleigh distribution. This analysis suggests that relatively few modes are contributing a large portion of the late reverberation energy.
Appendix A Lower bound of Pole Magnitude
We present lower and upper bounds on the pole magnitudes of an FDN. The bounds are based on the generalization of Rouché’s theorem to matrix polynomials.
Theorem 1 (see [41]).
Let and be matrix polynomials and let be a positive real number . If is positive definite for , then the polynomials and have the same number of roots of modulus less than .
An immediate consequence of the above theorem applied to the polynomial of (8) with and is [21]: If
(45) 
where means that is positive definite, then has no eigenvalues in the open disk with center and radius . The criterium in (45) is equivalent to
(46) 
which in turn is equivalent to [42]
(47) 
where denotes the spectral radius of a matrix . Using properties of the spectral norm [42] we can give an upper bound on this expression by
(48)  
Thus, the criterium (47) is satisfied if
(49) 
Therefore with (45), the pole magnitude lower bound may be given as
(50) 
Analogously, applying the same arguments to the reversed matrix polynomial yields an upper bound
(51) 
For additional attenuation filters as in (34) and a unitary feedback matrix , these bounds can be tightened further. Rouché’s criterion (45) gives the relation
(52) 
Thus, the lower bound of the pole magnitude is
(53) 
The corresponding upper bound may be derived similar to (51):
(54) 
These bounds are tight for a diagonal matrix .
Appendix B Far Deflation Estimation
We are given the equidistributed poles as defined in (23) and an even number of near poles . We compute the far deflation for pole . First we state a useful identity. For any real ,
(55) 
The total deflation is
Similarly, as each conjugate pair of poles contribute equally to the deflation, the far deflation is
(56) 
Appendix C Deflation Error
We show that inequality (30) is sufficient for inequality (26). For the sake of brevity, we omit the pole arguments in the following. Given the deflation approximation , which satisfy (30), we obtain
(57) 
We show that (57) satisfies EAI step error tolerance (26). Because , it is
(58) 
Further, as ,
Eventually, we can show that
References
 [1] K. Karplus and A. Strong, “Digital synthesis of pluckedstring and drum timbres,” Comput. Music J., vol. 7, no. 2, pp. 43–55, 1983.
 [2] J. O. Smith III, “Physical modeling using digital waveguides,” Comput. Music J., vol. 16, no. 4, pp. 74–91, 1992.
 [3] J. D. Parker, “Efficient dispersion generation structures for spring reverb emulation,” EURASIP J. Adv. Signal Process., vol. 2011, no. 1, pp. 1–8, 2011.
 [4] L. Savioja and U. P. Svensson, “Overview of geometrical room acoustic modeling techniques,” J. Acoust. Soc. Amer., vol. 138, no. 2, pp. 708–730, Aug. 2015.
 [5] V. Välimäki, J. D. Parker, L. Savioja, J. O. Smith III, and J. S. Abel, “Fifty years of artificial reverberation,” IEEE/ACM Trans. Audio, Speech, Lang. Proc., vol. 20, no. 5, pp. 1421–1448, Jul. 2012.
 [6] M. A. Gerzon, “Synthetic stereo reverberation: Part One,” Studio Sound, vol. 13, pp. 632–635, 1971.
 [7] J. M. Jot and A. Chaigne, “Digital delay networks for designing artificial reverberators,” in Proc. Audio Eng. Soc. Conv., Paris, France, Feb. 1991, pp. 1–12.
 [8] S. J. Schlecht and E. A. P. Habets, “Feedback delay networks: Echo density and mixing time,” IEEE/ACM Trans. Audio, Speech, Lang. Proc., vol. 25, no. 2, pp. 374–383, 2017.
 [9] ——, “On lossless feedback delay networks,” IEEE Trans. Signal Process., vol. 65, no. 6, pp. 1554–1564, Mar. 2017.
 [10] J.M. Adrien, “The missing link: modal synthesis.” MIT Press, May 1991, pp. 269–298.
 [11] J. He and Z.F. Fu, Modal Analysis. Oxford, UK: ButterworthHeinemann, 2001.
 [12] H. Kuttruff, Room Acoustics, Fifth Edition. CRC Press, Jun. 2009.
 [13] Y. Naka, A. A. Oberai, and B. G. ShinnCunningham, “Acoustic eigenvalues of rectangular rooms with arbitrary wall impedances using the interval Newtongeneralized bisection method,” J. Acoust. Soc. Amer., vol. 118, no. 6, pp. 3662–3671, Dec. 2005.
 [14] M. Karjalainen, P. A. A. Esquef, P. Antsalo, A. Makivirta, and V. Välimäki, “AR/ARMA analysis and modeling of modes in resonant and reverberant systems,” in Proc. Audio Eng. Soc. Conv., Munich, Germany, May 2002, pp. 1–16.
 [15] D. Beaton and N. Xiang, “Room acoustic modal analysis using Bayesian inference,” J. Acoust. Soc. Amer., vol. 141, no. 6, pp. 4480–4493, Jun. 2017.
 [16] Y. Haneda, S. Makino, and Y. Kaneda, “Common acoustical pole and zero modeling of room transfer functions,” IEEE Trans. Speech, Audio Process., vol. 2, no. 2, pp. 320–328, Apr. 1994.
 [17] D. Rocchesso and J. O. Smith III, “Circulant and elliptic feedback delay networks for artificial reverberation,” IEEE Trans. Speech, Audio Process., vol. 5, no. 1, pp. 51–63, 1997.
 [18] S. J. Schlecht and E. A. P. Habets, “Timevarying feedback matrices in feedback delay networks and their application in artificial reverberation,” J. Acoust. Soc. Amer., vol. 138, no. 3, pp. 1389–1398, Sep. 2015.
 [19] G. H. Golub and C. F. Van Loan, Matrix Computations. Baltimore : Johns Hopkins University Press, 1996.
 [20] J. M. McNamee, Numerical Methods for Roots of Polynomials, ser. Part I. New York, USA: Elsevier, Aug. 2007.
 [21] D. A. Bini and V. Noferini, “Solving polynomial eigenvalue problems by means of the Ehrlich–Aberth method,” Linear Algebra Appl., vol. 439, no. 4, pp. 1130–1149, Aug. 2013.
 [22] N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. Society for Industrial and Applied Mathematics, May 2012.
 [23] G. W. Stewart, “On the adjugate matrix,” Linear Algebra Appl., vol. 283, no. 13, pp. 151–164, Nov. 1998.
 [24] O. Aberth, “Iteration Methods for Finding all Zeros of a Polynomial Simultaneously ,” Mathematics of Computation, vol. 27, no. 122, pp. 339–344, 1973.
 [25] D. A. Bini, “Numerical computation of polynomial zeros by means of Aberth’s method,” Numerical Algorithms, vol. 13, no. 2, pp. 179–200, Feb. 1996.
 [26] A. Horn, “On the eigenvalues of a matrix with prescribed singular values,” Proc. Amer. Math. Soc., vol. 5, no. 1, pp. 4–7, 1954.
 [27] J. G. Proakis and D. G. Manolakis, Digital Signal Processing, 4th ed., Upper Saddle River, New Jersey, 2007, vol. 2.
 [28] Y. Ma, J. Yu, and Y. Wang, “Efficient recursive methods for partial fraction expansion of general rational functions,” Journal of Applied Mathematics, vol. 2014, pp. 1–18, Oct. 2014.
 [29] B. Bank, “Converting infinite impulse response filters to parallel form [Tips & Tricks],” IEEE Signal Process. Mag., vol. 35, no. 3, pp. 124–130, May 2018.
 [30] J. A. Foster, J. G. McWhirter, M. R. Davies, and J. A. Chambers, “An algorithm for calculating the QR and singular value decompositions of polynomial matrices,” IEEE Trans. Signal Process., vol. 58, no. 3, pp. 1263–1274, Mar. 2010.
 [31] B. Shahrrava, “Closedform impulse responses of linear timeinvariant systems: A unifying approach [Lecture Notes],” IEEE Signal Process. Mag., vol. 35, no. 4, pp. 126–132, Jul. 2018.
 [32] S. J. Schlecht and E. A. P. Habets, “Accurate reverberation time control in feedback delay networks,” in Proc. Int. Conf. Digital Audio Effects (DAFx), Edinburgh, UK, Aug. 2017, pp. 337–344.
 [33] J. M. Jot, “Proportional parametric equalizers  Application to digital reverberation and environmental audio Processing,” in Proc. Audio Eng. Soc. Conv., New York, NY, USA, Oct. 2015, pp. 1–8.
 [34] M. R. Schroeder, “New method of measuring reverberation time,” J. Acoust. Soc. Amer., vol. 37, no. 3, pp. 409–412, 1965.
 [35] E. De Sena, H. Hacıhabiboğlu, Z. Cvetkovic, and J. O. Smith III, “Efficient synthesis of room acoustics via scattering delay networks,” IEEE/ACM Trans. Audio, Speech, Lang. Proc., vol. 23, no. 9, pp. 1478–1492, 2015.
 [36] H. Bai, G. Richard, and L. Daudet, “Late reverberation synthesis: From radiance transfer to feedback delay networks,” IEEE/ACM Trans. Audio, Speech, Lang. Proc., vol. 23, no. 12, pp. 2260–2271, 2015.
 [37] M. Karjalainen and H. Jarvelainen, “More about this reverberation science: Perceptually good late reverberation,” in Proc. Audio Eng. Soc. Conv., New York, NY, USA, Nov. 2001, pp. 1–8.
 [38] J. A. Moorer, “About this reverberation business,” Comput. Music J., vol. 3, no. 2, pp. 13–17, Jun. 1979.
 [39] J. O. Smith III, Physical Audio Signal Processing, ser. For Virtual Musical Instruments And Audio Effects. W3K Publishing, 2010.
 [40] B. Rivet, L. Girin, and C. Jutten, “LogRayleigh distribution: A simple and efficient statistical representation of logspectral coefficients.” IEEE/ACM Trans. Audio, Speech, Lang. Proc., vol. 15, no. 3, pp. 796–802, 2007.
 [41] Y. Monden and S. Arimoto, “Generalized Rouche’s theorem and its application to multivariate autoregressions,” IEEE Trans. Acoust., Speech, Signal Process., vol. 28, no. 6, pp. 733–738, Dec. 1980.
 [42] R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. Cambridge University Press, 2013.