Modal Decomposition of Feedback Delay Networks

TL;DR

This paper introduces an efficient numerical method for modal decomposition of feedback delay networks, enhancing sound synthesis and reverberation analysis by revealing their modal behavior.

Contribution

It presents a novel pole finding algorithm based on Ehrlich-Aberth iteration, significantly improving computational efficiency for FDN modal analysis.

Findings

The new method achieves up to three orders of magnitude faster computation.

Few modes dominate the impulse response energy in FDNs.

Modal analysis aids in optimizing artificial reverberation systems.

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 Ehrlich-Aberth 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.

Figures

Tables

Table 1. TABLE I: Probability 𝒞 HIST ​ ( κ ) subscript 𝒞 HIST 𝜅 \mathcal{C}_{\textrm{HIST}}(\kappa) of cluster numbers of mode frequencies

Cluster size $\kappa$	0	1	2	3	$\ge$ 4
Uniform Random	0.3690	0.3661	0.1854	0.0610	0.0186
Lossless 8-FDN	0.1694	0.6632	0.1653	0.0020	0.0001
Equidistributed	0.0000	1.0000	0.0000	0.0000	0.0000

Equations

$$h_i(n)=\text{absolutevalue}{\rho }_i\text{absolutevalue}{{\lambda }_i}^n{e}^{ı(n\angle {\lambda }_i+\angle {\rho }_i)},$$

$$h_{}(n)=\sum _{i=1}^{\mathfrak{N}}{h}_i(n)$$

$$y(n)={\mathbf{c}}^{\top }\mathbf{s}(n)+dx(n)$$

$$\mathbf{s}(n+\mathbf{m})=\mathbf{A}\mathbf{s}(n)+\mathbf{b}x(n),$$

$$H_{}(z)={\mathbf{c}}^{\top }{\left[{\mathbf{D}}_{\mathbf{m}}(z{)}^{-1}-\mathbf{A}\right]}^{-1}\mathbf{b}+d,$$

$$\mathfrak{N}\gg N.$$

$$H(z)=d+\sum _{i=1}^{\mathfrak{N}}\frac{{\rho }_i}{1-{\lambda }_i{z}^{-1}},$$

$$\mathrm{adj}(\mathbf{A})=\det \left(\mathbf{A}\right){\mathbf{A}}^{-1},$$

$$\mathbf{P}(z)={\mathbf{D}}_{\mathbf{m}}(z{)}^{-1}-\mathbf{A}.$$

$$H_{}(z)=\frac{q_{\mathbf{m},\mathbf{A},\mathbf{b},\mathbf{c},d}(z)}{p_{\mathbf{m},\mathbf{A}}(z)},$$

$$p_{\mathbf{m},\mathbf{A}}(z)=\det \left(\mathbf{P}(z)\right)$$

$$q_{\mathbf{m},\mathbf{A},\mathbf{b},\mathbf{c},d}(z)=d\det \left(\mathbf{P}(z)\right)+{\mathbf{c}}^{\top }\mathrm{adj}(\mathbf{P}(z))\mathbf{b}.$$

$$p_{\mathbf{m},\mathbf{A}}(z)=\det \left(z{\mathbf{I}}_{\mathfrak{N}}-\mathfrak{A}\right),$$

$$p_{\mathbf{m},\mathbf{A}}(z)=\sum _{i=0}^{\mathfrak{N}}{c}_i{z}^i,$$

$$\mathbf{P}(z)=\sum _{k=0}^K{\mathbf{P}}_k{z}^k,$$

$$p^{\prime }(z)=\text{derivative}zp(z)=$$

$$=$$

$${\lambda }_i^{(j+1)}={\lambda }_i^{(j)}-{\Delta }_i^{(j)}$$

$${\Delta }_i^{(j)}=\frac{\mathcal{N}\left({\lambda }_i^{(j)}\right)}{1-\mathcal{N}\left({\lambda }_i^{(j)}\right){\mathcal{D}}_i\left({\mathbf{\lambda }}^{(j)}\right)}.$$

$$\mathcal{N}\left(z\right)=\frac{p(z)}{p^{\prime }(z)}=\frac{1}{\text{tr}\left(\mathbf{P}(z{)}^{-1}{\mathbf{P}}^{\prime }(z)\right)}$$

$${\mathcal{D}}_i\left({\mathbf{\lambda }}^{(j)}\right)=\sum _{l=1,l\ne i}^{\mathfrak{N}}\frac{1}{{\lambda }_i^{(j)}-{\lambda }_l^{(j)}}.$$

$${\Delta }_i^{(j)}=\frac{1}{\text{tr}\left(\mathbf{P}{\left({\lambda }_i^{(j)}\right)}^{-1}{\mathbf{P}}^{\prime }\left({\lambda }_i^{(j)}\right)\right)-{\mathcal{D}}_i\left({\mathbf{\lambda }}^{(j)}\right)}.$$

$$\text{absolutevalue}{\Delta }_i^{(j)}\le {\tau }_2\text{absolutevalue}{\lambda }_i^{(j)},$$

$$\sqrt[\min \mathbf{m}]{\min \sigma (\mathbf{A})}\le \text{absolutevalue}{\lambda }_i\le \sqrt[\max \mathbf{m}]{\max \sigma (\mathbf{A})}.$$

$${\mathbf{\lambda }}^{(0)}=\exp \left(ı2\pi \left[\frac{0}{\mathfrak{N}},\frac{1}{\mathfrak{N}},\dots ,\frac{\mathfrak{N}-1}{\mathfrak{N}}\right]\right).$$

$$\mathbf{A}={\mathbf{I}}_{\text{S}}=\left[\begin{array}{cccccc}0& 1& 0& \cdots & \cdots & 0\\ 0& 0& 1& \cdots & \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & 1& 0\\ 0& 0& 0& \cdots & 0& 1\\ 1& 0& 0& \cdots & 0& 0\end{array}\right]$$

$$p_{\mathbf{m},{\mathbf{I}}_{\text{S}}}(z)=z^{\mathfrak{N}}-1.$$

$$\text{absolutevalue}{\Delta }_i^{(j)}-{\tilde{\Delta }}_i^{(j)}\le {\tau }_3.$$

$${\mathcal{D}}_i\left({\mathbf{\lambda }}_{\text{far}}^{(0)}\right)=\frac{1}{{\lambda }_i^{(0)}}\frac{\mathfrak{N}-{\mathfrak{N}}_{\text{near}}-1}{2}.$$

$$\widetilde{{\mathcal{D}}_i}\left({\mathbf{\lambda }}^{(j)}\right)={\mathcal{D}}_i\left({\mathbf{\lambda }}_{\text{near}}^{(j)}\right)+{\mathcal{D}}_i\left({\mathbf{\lambda }}_{\text{far}}^{(0)}\right)$$

$${ϵ}_{\mathcal{D}}\ge \underset{i,j}{\max }\text{absolutevalue}\widetilde{{\mathcal{D}}_i}\left({\mathbf{\lambda }}^{(j)}\right)-{\mathcal{D}}_i\left({\mathbf{\lambda }}^{(j)}\right).$$

Full text

I Introduction

A feedback delay network (FDN) consists of a set of delay lines with lengths $\bm{m}$ which are interconnected via a feedback matrix $\bm{A}$ (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 time-domain, e.g., with circular buffers resulting in a constant time complexity $\mathcal{O}(1)$ 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 one-pole resonators, so-called modes. The time-domain signal of such a resonator with pole $\lambda_{i}$ and residue $\rho_{i}$ is

[TABLE]

where $\angle$ indicates the argument of a complex number in radiant, $\absolutevalue{\cdot}$ is the magnitude, $\imath=\sqrt{-1}$ and $n$ indicates the discrete time index. Each individual resonating mode is governed by four parameters: mode frequency $\angle\lambda_{i}$, decay rate $\absolutevalue{\lambda_{i}}$, initial phase $\angle\rho_{i}$ and initial amplitude $\absolutevalue{\rho_{i}}$ (see Fig 1). With modal decomposition, we aim to uncover the specific parameters of each mode. The time-domain impulse response of the FDN

[TABLE]

is the sum of the complex modes $h_{i}(n)$, where $\mathfrak{{N}}$ 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 moving-average [14], Bayesian inference [15] and all-pole 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.

I-A Problem Statement

For a single input and single output, the time-domain recursion of an FDN with $N$ delay lines is given by

[TABLE]

where $n$ is the time index, $\cdot^{\mathstrut\scriptscriptstyle\top}$ denotes the transpose operation and $\bm{A}\in\mathbb{C}^{N\times N}$, $\bm{b},\bm{c},\bm{s}(n)\in\mathbb{C}^{N\times 1}$, $x(n),y(n),d\in\mathbb{C}$ [17]. The state vector is defined as $\bm{s}(n+\bm{m})=\left[s_{1}(n+m_{1}),\dots,s_{N}(n+m_{N})\right]$. We write $N$-FDN to denote an FDN of size $N$. The transfer function of an FDN is

[TABLE]

where $\bm{D}_{\bm{m}}(z)=\mathrm{diag}\left\lparen z^{-m_{1}},z^{-m_{2}},\dots,z^{-m_{N}}\right\rparen$. The system order is given by $\mathfrak{{N}}=\sum_{i=1}^{N}m_{i}$ [17]. For commonly used delays $\bm{m}$, the system order is much larger than the FDN size, i.e.,

[TABLE]

The modal decomposition of the FDN, i.e., the partial fraction decomposition (PFD) of the transfer function (4) is

[TABLE]

where $\rho_{i}$ is the residue of the pole $\lambda_{i}$. The time-domain 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).

I-B Direct Approach

We first review two standard methods for the modal decomposition [17, 18]. Let $\bm{A}$ be any invertible matrix, then

[TABLE]

where $\mathrm{adj}(\bm{A})$ is the adjugate of the matrix $\bm{A}$ [19]. In the following, we denote

[TABLE]

With (7) and (8), the transfer function (4) can be expressed as a rational polynomial

[TABLE]

where

[TABLE]

and

[TABLE]

For brevity, we occasionally omit the parameters and write $q(z)$ and $p(z)$. The FDN system poles $\lambda_{i}$, where $1\leq i\leq\mathfrak{{N}}$, are the roots of the generalized characteristic polynomial (GCP) $p_{\bm{m},\bm{A}}(z)$ in (10) such that they are fully characterized by the delay matrix $\bm{D}_{\bm{m}}(z)$ and the feedback matrix $\bm{A}$.

For a moment, let us assume that all delays are single time steps, i.e., $\bm{m}=\bm{1}$. The time-domain recursion in (3) reduces to the standard state-space description of a linear time-invariant (LTI) filter. The system poles $\lambda_{i}$ are the eigenvalues of the feedback matrix $\bm{A}$ such that the modal decomposition (6) is easily computed with standard methods. However, for longer delays $\bm{m}$ such that (5) holds, the modal decomposition becomes more involved.

The GCP $p_{\bm{m},\bm{A}}$ can be expressed in a linearized fashion

[TABLE]

where $\bm{I}_{\mathfrak{{N}}}$ is the identity matrix of size $\mathfrak{{N}}$ and ${\bm{\mathfrak{A}}}\in\mathbb{C}^{\mathfrak{{N}}\times\mathfrak{{N}}}$ such that the system poles are the eigenvalues of ${\bm{\mathfrak{A}}}$ [17]. Unfortunately, for large delays $\bm{m}$ this eigenvalue problem becomes quickly numerically intractable. Alternatively, the GCP can be expressed as a scalar polynomial

[TABLE]

where the coefficients $c_{i}$ are derived from the principal minors of $\bm{A}$ [18]. The system poles are the roots of the scalar polynomial. Again, the polynomial degree increases with longer delays $\bm{m}$ 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 $\mathfrak{{N}}$ and modest-sized $N$. 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 $p(z)$ and subsequently recover the residues $\rho_{i}$. We conclude this section with a generalization to additional filtering in the delay lines and feedback matrix.

II-A 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 I-B, we compute the system poles without expanding the problem. In fact, (10) is a polynomial eigenvalue problem of degree $K=\max\bm{m}$, i.e.,

[TABLE]

where $\bm{P}_{k}\in\mathbb{C}^{N\times N}$ for $0\leq k\leq K$. For a proper matrix polynomial $\bm{P}(z)$, i.e., $\det\left\lparen\bm{P}_{K}\right\rparen\neq 0$, the number of roots is $KN$ [21]. For FDNs, however, $\bm{P}_{K}$ is singular such that the actual number of roots is lower, respectively, many roots are infinite. In fact, if $\det\left\lparen\bm{A}\right\rparen\neq 0$, the number of finite roots is $\mathfrak{{N}}$ which is also the degree of the scalar polynomial in (10) [18].

In the following, we use the derivative of the polynomial $p(z)=\det\left\lparen\bm{P}(z)\right\rparen$. According to Jacobi’s formula [22], we have

[TABLE]

where $\bm{P}^{\prime}(z)=\derivative{\bm{P}(z)}{z}$ and $\tr(\bm{X})$ denotes the trace of matrix $\bm{X}$. Stewart [23] showed that the adjugate of $\bm{A}$ can be well-conditioned even when $\bm{A}$ is ill-conditioned, and he shows how $\mathrm{adj}(\bm{A})$ can be computed in a numerically stable way from a rank revealing decomposition of $\bm{A}$ [22]. With the definitions of the polynomial eigenvalue problem introduced, we present the proposed root finding algorithm.

II-B Ehrlich-Aberth Method

The polynomial eigenvalue problem can be solved with the Ehrlich-Aberth 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 $\bm{\lambda}^{(0)}\in\mathbb{C}^{\mathfrak{{N}}}$ be a vector of initial estimates for the $\mathfrak{{N}}$ roots of the polynomial $p(z)$ and $\bm{\lambda}^{(j)}=[\lambda_{1}^{(j)},\lambda_{2}^{(j)},\dots,\lambda_{\mathfrak{{N}}}^{(j)}]$ be the $j$-th EAI iteration. The EAI provides the sequence of estimates

[TABLE]

with the EAI step being

[TABLE]

Using the identity in (15), the Newton correction term is

[TABLE]

and the deflation term is

[TABLE]

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

[TABLE]

The method, given here in the Jacobi version, is known to converge cubically for simple roots and linearly for multiple roots [21]. The Gauss-Seidel version of EAI [21], which updates the estimates as soon as they become available, may converge even slightly faster.

II-C Stopping Criteria

The system poles $\lambda_{i}$ are the roots of the polynomial $p(z)$ in (10), i.e., $\det\left\lparen\bm{P}(\lambda_{i})\right\rparen=0$. In other words, $\bm{P}(\lambda_{i})$ is a singular matrix for all system poles. Thus, the natural stopping criteria is the reciprocal of the condition number $\kappa(\bm{P}(z))$ being less than a prescribed tolerance $\tau_{1}$. 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

[TABLE]

where $\tau_{2}$ 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].

II-D 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 $\sigma(\bm{A})$ of the feedback matrix (see Appendix A)

[TABLE]

Equation (22) is a generalization on the relation of eigenvalues and singular values as given in the Weyl-Horn Theorem [26] in the case of unit delays $\bm{m}=\bm{1}$. 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 $\max\sigma(\bm{A})$ may be arbitrarily large while all system poles lie on the unit circle [9]. For large delays $\bm{m}$ however, (22) shows that the pole magnitudes tend to be close to the unit circle.

We can further derive from (22) that if $\max\sigma(\bm{A})\leq 1$ 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 $\bm{A}$ being unitary, i.e., $\bm{A}^{\mathstrut\scriptscriptstyle H}\bm{A}=\bm{I}$, all system poles lie on the unit circle regardless of the delays $\bm{m}$. 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 $\bm{\lambda}^{(0)}$ uniformly on the unit circle. More precisely, we chose the roots of unity

[TABLE]

It is worthwhile to note that $\bm{\lambda}^{(0)}$ is the solution of a particular FDN with a circular shift matrix

[TABLE]

such that the GCP is

[TABLE]

The shift matrix thus combines the FDN delays into a single long delay line.

II-E Approximate Deflation

For a high system order $\mathfrak{{N}}$, the computational complexity of the deflation term (19) may become excessive. We propose an approximate deflation (AD) according to a maximum error tolerance $\tau_{3}$ for the resulting EAI step $\widetilde{\Delta}_{i}^{(j)}$ in (20), i.e.,

[TABLE]

The magnitude of the deflation term summands decreases with the pole distance $\lambda_{i}^{(j)}-\lambda_{l}^{(j)}$. The idea is then to divide the poles into a near and far pole sets, $\bm{\lambda}_{\textrm{near}}^{(j)}$ and $\bm{\lambda}_{\textrm{far}}^{(j)}$, respectively, and approximate the deflation of the less significant far poles $\bm{\lambda}_{\textrm{far}}^{(j)}$ by a default term, e.g., $\bm{\lambda}_{\textrm{far}}^{(0)}$. For symmetry, the number of near poles $\mathfrak{{N}}_{\textrm{near}}$ is assumed to be an even number.

It can be shown, that for equidistributed poles such as $\bm{\lambda}^{(0)}$ in (23), the far deflation is (see Appendix B)

[TABLE]

Thus, the total deflation may be approximated by

[TABLE]

if the far poles $\bm{\lambda}_{\textrm{far}}^{(j)}$ are sufficiently uniformly distributed. By sorting the system poles iterations along pole angles, we can find the near poles $\bm{\lambda}_{\textrm{near}}^{(j)}$. To establish the quality of this approximation, let us assume that there exists an upper bound $\epsilon_{\mathcal{D}}$ for the approximation error of the deflation term, i.e.,

[TABLE]

We can show that the error tolerance $\tau_{3}$ in (26) is satisfied if (see Appendix C)

[TABLE]

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 $\epsilon_{\mathcal{D}}$. As $\epsilon_{\mathcal{D}}$ depends on many factors such as matrix size $N$, system order $\mathfrak{{N}}$, feedback matrix $\bm{A}$ and number of near poles $\mathfrak{{N}}_{\textrm{near}}$, 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 $\bm{\lambda}^{(0)}$ are equidistributed, $\bm{\lambda}^{(1)}$ may be computed only from the estimated far deflation with $\mathfrak{{N}}_{\textrm{near}}=0$.

II-F Residues

Once we have found the system poles, the residues of the modal decomposition (6) are computed by

[TABLE]

where we assume that all poles are unique. Similar, but more intricate solutions exists for non-unique poles [28]. The undriven residue, i.e., the system response without excitation, is

[TABLE]

The undriven residue is a valuable intermediate step to analyze the mode initial amplitude independent from the input and output drives $q(\lambda_{i})$. Since, $\bm{P}(\lambda_{i})$ is a singular matrix, the derivative of the GCP $p^{\prime}(\lambda_{i})$ in (15) may only be computed by the adjugate formulation. Since $\det\left\lparen\bm{P}(\lambda_{i})\right\rparen=0$, the input-output drives in (11) are

[TABLE]

The difference between the driven and undriven residues may be expressed as a linear combination of the matrix entries of $\mathrm{adj}(\bm{P}(\lambda_{i}))$. Alternatively, the driven residues may also be computed by a least linear squares fit for the time-domain impulse response since the sum of complex resonators in (2) depends linearly on the residues [29].

II-G Polynomial Feedback and Delay Matrices

Although, the focus of this work is on frequency-independent feedback matrices $\bm{A}$, much of the development in Section II is applicable to general polynomial matrices. Therefore, it is easy to include further filtering such as a frequency-dependent feedback matrix $\bm{A}(z)$. There also exists a singular value decomposition for polynomial matrices $\bm{A}(z)$ [30]. Alternatively, the delay lines are often extended with an attenuation or allpass filter $\alpha_{i}(z)$, i.e.,

[TABLE]

It is important to note that additional filters may increase the number of system poles. Further, if $P(z)$ 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].

For a unitary feedback matrix $\bm{A}$ and for attenuation filters in (34), it is possible to improve the pole magnitude bounds (22) to

[TABLE]

where all vector operations are element-wise (see Appendix A).

III Modal Synthesis and Evaluation

The following evaluation uses real-valued FDN parameters such that the system poles appear in complex conjugate pairs.

III-A Modal Synthesis and Accuracy

A numerically accurate way to verify the modal decomposition is to synthesize each mode $h_{i}(n)$ in time-domain as expressed in (1) and compare the sum of all modes with the impulse response $h_{\textrm{}}(n)$ computed by the time-domain recursion in (3). The concept of modal synthesis and verification is depicted in Fig. 1. The error is given by the maximum difference111The 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.,

[TABLE]

In this work, we use double precision floating point arithmetic and the modal decomposition is regarded successful if the maximum error $\epsilon<10^{-10}$. The EAI is numerically stable as the matrix inversion in (20) is only necessary if the matrix is sufficiently non-singular due to the first stopping criteria. For large delays $\bm{m}$, the evaluation of $\bm{D}_{\bm{m}}(z)$ may become extremely large or small if $z$ 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 $\bm{m}$ 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.

III-B Numerical Evaluation

For the FDN, a single EAI step in (20) can be evaluated in $\mathcal{O}(\mathfrak{{N}}+N^{3})$: an evaluation of $\bm{P}(z)$ and $\bm{P}^{\prime}(z)$ is merely an evaluation of the delay matrix $\bm{D}_{\bm{m}}(z)$ in $\mathcal{O}(N)$; a numerical matrix inversion can be performed in $\mathcal{O}(N^{3})$; and the deflation term is evaluated in $\mathcal{O}(\mathfrak{{N}})$. Thus, a full iteration from $\bm{\lambda}^{(j)}\rightarrow\bm{\lambda}^{(j+1)}$ can be evaluated in $\mathcal{O}(\mathfrak{{N}}^{2}+\mathfrak{{N}}N^{3})$. This compares favorably with the bound $\mathcal{O}(\mathfrak{{N}}^{3})$ of a matrix-based algorithm applied to the linearization in (12). For a high number of system poles $\mathfrak{{N}}\gg N^{3}$, 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 $N$, total delays $\mathfrak{{N}}$ between 50 and $10^{4}$ 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 $\det\left\lparen\bm{A}\right\rparen$. 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 MATLAB222Matlab 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 $\mathfrak{{N}}$ were distributed randomly among eight delay lines and the feedback matrix $\bm{A}\in\mathbb{R}^{8\times 8}$ was a random orthogonal matrix. All methods gave the correct answer with the required accuracy. For the approximate deflation, the number of near poles $\mathfrak{{N}}_{\textrm{near}}$ was set to $\mathfrak{{N}}/100$. The maximum deflation error $\epsilon_{\mathcal{D}}=10^{3}$ was determined a priori by probing an independent set of random FDNs of similar configuration. The EAI step tolerance $\tau_{3}$ was set to $10^{-3}$.

The EAI implementation utilizes only standard MATLAB functions and no C-optimization which explains relatively poor performance for small system order $\mathfrak{{N}}<10^{3}$. For high system order $\mathfrak{{N}}$ such as $~{}5\cdot 10^{4}$, 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 $\mathfrak{{N}}$ and cubic in $N$ such that it is possible to perform modal decomposition up to $\mathfrak{{N}}=10^{6}$. Whereas the memory requirements for eig become prohibitive for $\mathfrak{{N}}>5\cdot 10^{4}$. For $\mathfrak{{N}}>10^{5}$, more the $95\%$ 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.

IV-A Attenutation

Attenuation filters in FDNs, as they are typically applied in artificial reverberation, aim to control the frequency-dependent reverberation time [7, 32]. As expressed in (34), all delays are extended with absorption filters $\bm{\alpha}(z)$ and the feedback matrix $\bm{A}$ is orthogonal or more generally unilossless [9]. We study three types of attenuation: homogeneous, near-homogeneous and inhomogeneous attenuation.

IV-A1 Homogeneous Attenuation

The attenuation filters $\bm{\alpha}(z)$ are called homogeneous if there exists an attenuation-per-sample $\Gamma(z)$ such that

[TABLE]

The attenuated delay lines can be expressed as plain delay lines with a mapped argument, i.e.,

[TABLE]

Consequently, the system poles with attenuation $\lambda^{\Gamma}_{i}$ can be related to the system poles $\lambda_{i}$ without attenuation by

[TABLE]

If we assume that the attenuation filters have a purely real frequency response333Although such a frequency response is not realizable with a digital filter in general, but useful for the theoretical analysis., i.e., $\Gamma(e^{\imath\omega})\in\mathbb{R}$ then the mode frequencies are unaltered by the attenuation, i.e., $\angle\lambda_{i}=\angle\lambda^{\Gamma}_{i}$. For a unilossless $\bm{A}$, all unattenuated system poles $\lambda_{i}$ lie on the unit circle such that

[TABLE]

and the attenuated FDN is stable if $\absolutevalue{\Gamma(e^{\imath\omega})}<1$. For homogeneous attenuation, the magnitude bounds in (35) are tight.

IV-A2 Near-homogeneous Attenuation

Typically, the attenuation filters are implemented with relatively low order, such as one-pole filters [7], however higher order filters were proposed as well [33, 32]. The attenuation filters are designed to match the magnitude response

[TABLE]

where the attenuation-per-sample is derived from a target reverberation time

[TABLE]

where $f_{s}$ is the sampling frequency and $T_{60}(\omega)$ is the time in seconds for the energy decay curve of the impulse response at frequency $\omega$ to decay by 60 dB [34]. For illustration, we compute the one-pole filter according to [7] given the target reverberation time at DC $T_{60}(0)$ and Nyquist frequency $T_{60}(\pi)$. Figure 4 depicts the resulting modal decomposition for an 8-FDN with an orthogonal feedback matrix and a target reverberation time $T_{60}(0)=2$ seconds and $T_{60}(\pi)=0.4$ 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. 4b. For near-homogeneous 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 $z^{-\bm{m}}$ 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.

IV-A3 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 8-FDN as in Fig. 4 with different attenuation filters. Instead of the delay proportional design in (41), all one-pole 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.

IV-B Statistical Distribution of Poles and Residues

We present a set of statistical analyses of lossless FDNs which rely on the proposed large-scale numerical computation of the modal decomposition and are difficult to derive by analytic methods. The statistical analysis answers a long-standing 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? 2. 2)

What is the distribution of residue magnitudes?

In the analyses, we rely on Monte Carlo simulations of randomly generated lossless FDNs.

IV-B1 Mode Frequency Distribution

The near-equidistribution 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

[TABLE]

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 $\omega$. In contrast, a mode gap occurs if $\mathcal{C}(\omega)=0$, i.e., no mode lies in this frequency interval. For perfectly equidistributed poles $\mathcal{C}(\omega)=1$ for all $\omega$. We evaluate the distribution of mode frequencies by computing the histogram of cluster numbers

[TABLE]

where $\delta(\cdot)$ is the dirac function, $\kappa$ is the integer cluster size, and $\mathcal{L}$ is the number of observations. For large enough $\mathcal{L}$, the histogram converges towards the probability of cluster numbers. The random 8-FDNs 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 pseudo-uniform 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 $2\pi/\mathfrak{{N}}$ with more than two modes. In stark contrast, acoustic mode density of physical spaces increase quadratically with frequency [12].

IV-B2 Residue Magnitude Distribution

In Section II-F, 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 input-output drives for a random 8-FDN. The input-output drives are comprised of all individual input-output combinations, i.e., $\mathrm{adj}(\bm{P}(\lambda_{i}))$ in (33). The total residues $\rho(\lambda_{i})$ result from unit input and output gains, i.e., $\bm{b}=\bm{1}$ and $\bm{c}=\bm{1}$, or in other words, $\rho(\lambda_{i})=\rho^{\mathrm{u}}(\lambda_{i})(\bm{1}^{\mathstrut\scriptscriptstyle\top}\mathrm{adj}(\bm{P}(\lambda_{i}))\bm{1})$. The magnitude distributions of the inverse undriven residues $1/\rho^{\mathrm{u}}(\lambda_{i})$, the total residues $\rho(\lambda_{i})$ and input-output drives $\mathrm{adj}(\bm{P}(\lambda_{i}))$ all resemble log-Rayleigh distributions [40]. However, just by altering the feedback matrix $\bm{A}$, 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.

IV-B3 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 high-energy 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 online444www.audiolabs-erlangen.de/resources/2018-IEEE-Modal.

The residue distribution may be optimized in two steps: Firstly, optimization of the undriven residues by choosing delays $\bm{m}$ and feedback matrix $\bm{A}$. Secondly, optimization of the total residue by choosing the input and output gains, $\bm{b}$ and $\bm{c}$. While the first step is a non-linear 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 Ehrlich-Aberth 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 $~{}5\cdot 10^{4}$, 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, near-homogeneous 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 log-Rayleigh 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 $\absolutevalue{\bm{\lambda}}$ of an FDN. The bounds are based on the generalization of Rouché’s theorem to matrix polynomials.

Theorem 1** (see [41]).**

Let $\bm{S}(z)$ and $\bm{Q}(z)$ be matrix polynomials and let $r$ be a positive real number . If $\bm{S}(z)^{\mathstrut\scriptscriptstyle H}\bm{S}(z)-\bm{Q}(z)^{\mathstrut\scriptscriptstyle H}\bm{Q}(z)$ is positive definite for $|z|=r$, then the polynomials $\det\left\lparen\bm{S}(z)\right\rparen$ and $\det\left\lparen\bm{S}(z)+\bm{Q}(z)\right\rparen$ have the same number of roots of modulus less than $r$.

An immediate consequence of the above theorem applied to the polynomial $\bm{P}(z)$ of (8) with $\bm{S}(z)=-\bm{A}$ and $\bm{Q}(z)=\bm{D}_{\bm{m}}(z)^{-1}$ is [21]: If

[TABLE]

where $\bm{X}\succ\bm{Y}$ means that $\bm{X}-\bm{Y}$ is positive definite, then $\bm{P}(z)$ has no eigenvalues in the open disk with center [math] and radius $r$. The criterium in (45) is equivalent to

[TABLE]

which in turn is equivalent to [42]

[TABLE]

where $\rho(\bm{X})$ denotes the spectral radius of a matrix $\bm{X}$. Using properties of the spectral norm [42] we can give an upper bound on this expression by

[TABLE]

Thus, the criterium (47) is satisfied if

[TABLE]

Therefore with (45), the pole magnitude lower bound may be given as

[TABLE]

Analogously, applying the same arguments to the reversed matrix polynomial $z^{\max\bm{m}}\bm{P}(z^{-1})$ yields an upper bound

[TABLE]

For additional attenuation filters $\bm{\alpha}(z)$ as in (34) and a unitary feedback matrix $\bm{A}$, these bounds can be tightened further. Rouché’s criterion (45) gives the relation

[TABLE]

Thus, the lower bound of the pole magnitude is

[TABLE]

The corresponding upper bound may be derived similar to (51):

[TABLE]

These bounds are tight for a diagonal matrix $\bm{A}$.

Appendix B Far Deflation Estimation

We are given the equidistributed poles $\bm{\lambda}^{(0)}$ as defined in (23) and an even number of near poles $\mathfrak{{N}}_{\textrm{near}}$. We compute the far deflation for pole $\lambda_{i}^{(j)}$. First we state a useful identity. For any real $x$,

[TABLE]

The total deflation is

[TABLE]

Similarly, as each conjugate pair of poles contribute equally to the deflation, the far deflation is

[TABLE]

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 $\widetilde{\mathcal{D}_{i}}$, which satisfy (30), we obtain

[TABLE]

We show that (57) satisfies EAI step error tolerance (26). Because $\epsilon_{\mathcal{D}}\geq 0$, it is

[TABLE]

Further, as $\absolutevalue{\mathcal{D}_{i}-\widetilde{\mathcal{D}_{i}}}-\epsilon_{\mathcal{D}}\leq 0$,

[TABLE]

Eventually, we can show that

[TABLE]