Grouped Feedback Delay Networks With Frequency-Dependent Coupling

Feedback Delay Networks are one of the most popular and efficient means of generating artificial reverberation. Recently, we proposed the Grouped Feedback Delay Network (GFDN), which couples multiple FDNs while maintaining system stability. The GFDN can be used to model reverberation in coupled spaces that exhibit multi-stage decay. The block feedback matrix determines the inter- and intra-group coupling. In this article, we expand on the design of the block feedback matrix to include frequency-dependent coupling among the various FDN groups. We show how paraunitary feedback matrices can be designed to emulate diffraction at the aperture connecting rooms. Several methods for the construction of nearly paraunitary matrices are investigated. The proposed method supports the efficient rendering of virtual acoustics for complex room topologies in games and XR applications.


I. INTRODUCTION
D ELAY network reverberators, such as the Feedback Delay Network (FDN), are often preferred for flexible, real-time simulations because of their computational efficiency. Feedback delay networks are efficient IIR structures for synthesizing room impulse responses (RIRs). RIRs consist of a set of sparse early reflections which increase in density over time, building toward late reverberation where the impulse density is high and statistically Gaussian. Feedback delay networks are composed of delay lines in parallel, which are connected through a feedback matrix (or mixing matrix), which is unitary to conserve system energy [1]. Schlecht and Habets expanded the set of matrices that preserve energy in FDNs to unilossless matrices [2]. The input signal traverses through the delay lines and the mixing matrix, building echo density over time. Jot proposed adding decay filters to the delay lines to yield a desired frequency dependent T 60 [3], [4]. Since then, FDNs have become one of the most popular structures for synthesizing reverberation due to the relative efficiency of the approach. Recent research on FDNs has focused on mixing matrix design to increase echo density [5], modal analysis [6], [7], timevarying FDNs [8], allpass FDNs for colorless reverberation [9] and scattering FDNs for dense reverberation [10]. Grouping of delay lines to control direction-dependent energy decay, known as the Directional Feedback Delay Network, was proposed by Alary et al. in [11].
Coupled spaces with distinctly different reverberation times exhibit multi-stage decay [12]. Coupled spaces are found in apartments, offices, concert halls, and churches. Eyring developed a model for energy decay in coupled rooms and rooms with distinctly different absorbing surfaces [13]. Xiang's work [12], [14], [15] has shed further light on the analysis and statistical modeling of energy decay in coupled spaces. In [16], a roundrobin study was conducted comparing coupled volume simulations in various geometrical and wave-based room acoustics software. The simulation of coupled room acoustics with delay network architectures is less investigated.
We proposed the Grouped Feedback Delay Network (GFDN) [17] to connect multiple FDNs whilst maintaining system stability. In a GFDN, groups of delay lines have the same target T 60 response associated with them, compared to traditional FDNs, in which all delay lines share the same decay characteristics. The interaction among the different delay line groups is controlled by a block mixing matrix. The diagonal sub-matrices of the block mixing matrix control the mixing in the individual rooms or among walls made of different materials, whereas the off-diagonal sub-matrices control the coupling among the rooms or wall materials. The GFDN was used to model coupled rooms and a room made of different materials where multiple concurrent decay times are observed. The discussion on coupled room modeling was further elaborated in [18] by comparing the GFDN's behavior to measurements. Similarly, a coupled volume extension to the scattering delay network was proposed recently in [19].
Schlecht has proposed adding delays [5] and filters [10] to the feedback matrix in an FDN. This has the advantage of increasing the echo density of the RIR more quickly, leading to dense, rich-sounding reverberation. Inspired by the idea, we incorporate frequency-dependent feedback matrices into the GFDN. Paraunitary matrices constitute a subset of lossless matrices. If the frequency-dependent feedback matrix is lossless, the overall T 60 is preserved despite introducing filtering effects [10]. This article describes how filter feedback matrices (FFMs) can be incorporated into the GFDN to preserve stability and emulate wave effects, such as diffraction at the aperture connecting coupled rooms. r Give applications of frequency-dependent feedback matrices by modeling diffraction in coupled rooms. Although this article focuses on coupled room modeling, the paraunitary matrix design principles discussed here have a broader application area. They can be used to model any frequency-dependent, lossless energy exchange between coupled systems, such as bridge-string coupling in musical instruments [20] or frequency-dependent scattering from objects [21].
The rest of the article is organized as follows. In Section II-A, the GFDN architecture is introduced. In Section III, some theorems related to paraunitary matrix construction are used to design a 2-tap FIR paraunitary feedback matrix with controllable frequency response. The rules of constructing a paraunitary matrix for coupling 2 FDNs are derived in Section III-C, followed by a general derivation for coupling K FDNs in Section III-D.
Here, a discussion on methods for solving the FIR paraunitary approximation problem is presented. The design of filter feedback matrices to model diffraction at the aperture connecting coupled rooms is discussed in Section IV. Asymmetric coupling and accurate scaling for level correction are discussed in Section IV-B. Section V presents GFDN evaluations modeling 2 and 3 coupled rooms. Section VI concludes the article.

A. Grouped Feedback Delay Network
A standard feedback delay network consists of N delay lines of length τ i samples each, i = 1, 2, . . . , N, with its associated decay filter, g i (z), z = exp(jω), connected through an N × N feedback matrix, M . For a frequency-dependent T 60 (z) in seconds, the decay filter gains are related to the delay line length as [3] where F s stands for the sampling frequency. In the grouped feedback delay network (GFDN), we use multiple sets of delay lines each with a different T 60 (z). In Fig. 1, a GFDN with two sets of delay lines is shown. Out of the total of N delay lines, N 1 delay lines have a decay response, T 60 1 (z), and N 2 delay lines have a decay response, T 60 2 (z), such that N 1 + N 2 = N . The two groups of decay filter gains, g 1 (z) and g 2 (z) are calculated according to the different T 60 (z)s, see (1). The mixing matrix M is now an N × N block matrix made of 4 submatrices, M ij ∈ R N i ×N j , i, j = 1, 2, where block diagonal matrices M 11 , M 22 control intra-group coupling and the block off-diagonal matrices M 12 , M 21 control inter-group coupling. With where I is the identity matrix and 0 is a matrix of zeros.

B. Coupled Mixing Matrix Design
The mixing matrix, or the feedback matrix, is an important property of the FDN, which determines the amount of coupling between various delay lines and system stability. This matrix controls the rate at which the echo density increases.
In the GFDN, we can choose different, independent unitary mixing matrices for each delay line group (the diagonal submatrices M 11 and M 22 ). The off-diagonal submatrices (M 12 and M 21 ) control how strongly coupled the groups are to each other. This gives us independent control over the intra-and inter-group mixing characteristics. To preserve energy in the system, the coupled mixing matrix needs to be unilossless [2]. Examples of such lossless matrices are orthonormal or unitary matrices [1], [2].
For two coupled rooms, each with individual mixing matrices M 2 1 and M 2 2 , and a coupling coefficient of α, the coupled mixing matrix can be characterized as [17], Here, φ = arctan(α) is the coupling angle. φ = 0 gives zero coupling, and φ = π 4 gives maximum coupling between the two FDNs. Note that for any unitary matrix, M , its square, M 2 , is also unitary, since its eigenvalues are the squares of unity.
In [18], the derivation of the orthonormal coupled mixing matrix was extended to an arbitrary number of coupled rooms. It was shown that a coupled feedback matrix, parameterized by a coupling coefficient matrix and individual unitary feedback matrices for each room, is unitary if the coupling coefficient matrix is unitary by design.

A. Lossless Coupling Matrix
Let us consider the case of K coupled rooms with a total number of N delay lines in the GFDN. With frequency-dependent coupling coefficients, a polynomial mixing matrix of degree p, F (z) ∈ R N ×N ×(p+1) , is defined as, In the above equation, M i is any unitary matrix that determines how diffuse the ith room is, and M 2 i is its square that is also unitary. Note that M i and M j need to have the same dimensions, i.e., each group in the GFDN should have the same number of delay lines. φ ij (z) is an FIR filter of degree p that determines the coupling between rooms i and j. We The feedback matrix is energy-preserving if it is paraunitary [5]. A paraunitary matrix must satisfy the following criteria where ω ∈ [−π, π] denotes frequency in rad/s. The Hermitian operator F (z) H denotes the transposition of the elements of F (z) and the complex conjugation of their coefficients.
Theorem 1 states that the feedback matrix is guaranteed to be paraunitary if the matrix of coupling filters is paraunitary. The proof of Theorem 1 is included in Appendix A. This gives us a paradigm for designing lossless frequency-dependent feedback matrices for the GFDN by simply designing a paraunitary coupling matrix.
Paraunitary coupling matrices ensure that the desired T 60 is maintained while introducing frequency-dependent effects. Let us recall that the maximum deviation in the FDN mode magnitude is given by [6], where Λ represents the FDN modes, σ min (F (z)) and σ max (F (z)) are the minimum and maximum singular values of the feedback matrix F (z), and τ min and τ max are the shortest and longest delay line lengths. In the undamped case, the largest and smallest singular values of F (z) correspond to the largest and smallest singular values of Φ(z). In the damped case, the mode amplitude decays to −60 dB after T 60 F s samples.
where g(e jω ) is the frequency response of the absorption filter in the delay lines. The deviation from the desired reverberation time, T 60 des (e jω) is, γ where γ = ln(0.001)/F s . Since paraunitary matrices have singular values with unit magnitude, the achieved reverberation time is equal to the desired reverberation time, T 60 (e jω ) = T 60 des (e jω ). This expression is generally useful for finding the deviation in the reverberation time from its desired value based on the feedback matrix.

B. Simplest Paraunitary Matrix Construction
First, we report some theorems related to the construction of paraunitary matrices [22] in Appendix B. Using these theorems, a first-order (two-tap) paraunitary FIR coupling matrix is designed for coupling two FDNs. Let us recall the 2 × 2 orthonormal coupling matrix, characterized by the coupling The following paraunitary matrix can be constructed from this orthogonal set of idempotents, such that Since the filter coefficients are complex, the real and imaginary parts of the incoming and outgoing signals from the feedback matrix must be propagated. This increases the computational complexity of the GFDN. However, only the real part of the signal is used when taking the final output, y(n) in Fig. 1.
The magnitude response of the paraunitary matrix is explicitly derived in (11). The parameter β controls the position of the coupling filter notches. The magnitude response and the filter coefficients in the complex plane for two different values of φ and increasing values of β ∈ [0, π 2 ] are shown in Fig. 2. The diagonal filters have an almost flat magnitude response. When β = 0, the response is a straight line with a 6 dB/octave slope. The parameter φ controls the radii of the filter coefficients and, in turn, the average magnitude of the off-diagonal filters. This provides a method for a very efficient and flexible design of a parameterized FIR coupling matrix that is paraunitary and requires only two taps. The diagonal filter response is flat, and the off-diagonal response is a notch filter whose cut-off frequency can be controlled by the parameter β. It can be implemented in real-time and has been incorporated into a GFDN plugin [23]. Higher-order causal filters for more than two coupled FDNs can be designed by following the rules of Theorem 3 (see Appendix B).

C. Coupling Matrix for 2 FDNs With Desired Coupling Filter
While the previous section discussed the design of a coupling matrix for a specific low-pass filter, the general principles of designing a 2 × 2 paraunitary coupling matrix with any desired coupling filter, α(z), is discussed here. The coupling filter determines the energy exchange between the two FDNs. First, the coupling filters need to be placed on the off-diagonal elements of the coupling matrix. Then, the diagonal elements must be designed so that the resulting matrix is paraunitary. Let us denote the diagonal elements as β(z), and its in-place time flipped version as, β (z) = z −L β(z −1 ). We assume that all the filter coefficients are real. The coupling matrix and its para-hermitian version are, . (12) Now, the product of these matrices, which should be identity, is given by where we have omitted the dependence on z for brevity. We desire Φ(z) to be paraunitary, i.e., and equally is also real-valued and nonnegative. Therefore, the Fejér-Riesz lemma (see Appendix C) applies, and a filter β(z) can be derived from β(z)β(z −1 ). The filter β(z) is not unique. Note that regardless of the choice of the phase of the filter, the matrix in (13) is identity. This is because, for a real filter with an antisymmetric phase, the phase cancels out in the convolution of the filter with its time-flipped version. The coupling matrix thus constructed from α(z) and β(z) is paraunitary if |α(z)| ≤ 1 and |β(z)| ≤ 1.

D. Coupling Matrix for K FDNs With Desired Coupling Filters
In this section, we discuss the design of a paraunitary FIR coupling matrix to model room-to-room coupling for each pair of coupled rooms in a K coupled room system. An initial, non-paraunitary matrix is constructed starting from desired coupling filters. It is preconditioned to be as close to paraunitary as possible. Then, its closest paraunitary counterpart is found. This ensures that the coupling matrix is lossless while maintaining the desired frequency response of the off-diagonal elements (that represent the coupling filters).
The first part of the solution is to design an initial polynomial matrix based on chosen coupling filters. For K coupled rooms, the design of an initial, non-paraunitary coupling matrix, A(z) ∈ R K×K×P , involves placing the desired coupling filters in the off-diagonal elements. Note that this matrix is symmetric due to acoustic reciprocity. Here, P is the degree of the polynomial matrix. Let us denote the Discrete Fourier Transform (DFT) of this matrix at M discrete frequency bins (M = 2P ) as A(e jω m ) ∈ C K×K , m = 0, . . . , M − 1.
We have some freedom in choosing the diagonal elements and the sign of the matrix elements. We want to precondition this matrix as close to paraunitary as possible. A necessary condition of unitary matrices is that the squared modulus of each row and column sum to one. The diagonal entries are designed so that the total energy of each row of the matrix is matched to one in the frequency domain: The only other degree of freedom is the choice of the sign of each matrix element. For a symmetric K × K matrix, a total of 2 K(K−1) sign combinations of ±1 are possible. We represent this with an K × K sign matrix, S, S ij ∈ {−1, 1}. We check for all possible sign combinations and choose the one that minimizes the distance of the initial matrix from a unitary matrix at each frequency bin, This is equivalent to finding a matrix that has the minimum Frobenius error from its Procrustes solution at every frequency bin (see Appendix D).
If K is moderately large, then checking for all possible sign combinations is computationally expensive. In [24], it has been shown that a globally optimal solution exists with a positive first row and column; this reduces the search space to 2 (K−1)(K−2) sign matrices. Furthermore, every pair of sign row vectors (or sign column vectors) allows orthogonality if neither all signs are equal nor opposite. These are known as sign potentially orthogonal matrices [24]. This further reduces the possible S choices.
The second part of the solution is to find a paraunitary matrix that is closest to this preconditioned matrix. A solution based on iterative and greedy optimization has been proposed by Tkacenko et al. [25]. They minimize a weighted mean-squared error given by the Frobenius norm of the error between the magnitude response of a desired polynomial matrix, A(e jω ), and a paraunitary matrix with frequency response Φ(e jω ). The optimization problem is, The authors exploit the complete parameterization of causal FIR paraunitary systems in terms of Householder-like degreeone building blocks [26], where a paraunitary matrix can be constructed from the cascade of degree-one Householder-like building blocks, and a unitary matrix. The Householder building blocks and unitary matrix are solved for iteratively with a greedy approach.
Vouras et al. [27] simplify the problem further by replacing the integration over all frequencies to be a summation over the discrete frequency bins, where ω m = 2πk M , and the weighting function is chosen to be unity at the desired frequency bins. This reduces to a linear combination of individual Procrustes problems, which can be solved separately.
The standard Procrustes problem and its solution are, where the subscript m is the shorthand for e jω m , A m = U m Σ m V H m , i.e., U m and V m are full-rank, square matrices that are the left and right singular vectors of A m [28] and Φ m is the optimal solution at frequency bin. The matrix thus designed, Φ(e jω m ), can be inverted using the Inverse Discrete Fourier Transform (IDFT) to find its time-domain counterpart. This matrix is unitary only at discrete points on the unit circle, and losslessness is not guaranteed in-between frequency bins. However, if M is large, then the sampling is fine enough to approximate paraunitarity for many practical applications [27].

IV. NEARLY PARAUNITARY ROOM COUPLING MATRIX TO MODEL DIFFRACTION
Now that we know how to design paraunitary matrices from a (preconditioned) desired matrix, we investigate how to design the desired matrix for a practical case of interest -modeling diffraction between coupled rooms. In the case of coupled spaces connected through an aperture, the diffraction of sound waves around the aperture should be modeled. Neither of the delay network coupled volume simulators proposed in the literature so far -GFDNs and coupled-volume SDNs -have considered diffraction. The net perceptual effect of diffraction is to dampen some frequencies; frequencies are damped depending on the ratio between the aperture size and the wavelength. Since diffraction is a wave phenomenon that exhibits frequency-dependent effects, it can be modeled in the coupling matrix with filters.
Morse and Ingard were the first to use transmission coefficients to account for wave effects at the coupling apertures [29]. We rely on the diffraction model proposed by Miles in [30]. The variations in the power transmission coefficients of the coupling apertures can affect variations in the shapes of decay curves, as shown in [31].

A. Transmission Model
The energy transmission coefficient, ζ, through an aperture connecting two rooms is a function of the angle of incidence with respect to the normal, η, aperture size, a and the wave number, k [31]. When the wavelength is comparable to the size of the aperture, diffraction (bending of waves around the aperture) can be observed. The model for a circular aperture is [30]: where σ is the transmission cross-section, k is the wave number (= 2πf c ), f is the frequency of the wave in Hz, c is the speed of sound in air.
We design an FIR filter from (21) to mimic the frequencydependent diffraction effects at the aperture between two rooms. With increasing mixing in an FDN, reflections from all directions Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. will be incident on the aperture. Assuming an isotropic diffuse field, we can integrate the transmission coefficient over all angles of incidence and elevation angles to find only a frequencydependent value [31], This integral can be solved numerically for discrete frequencies. The time domain FIR filter coefficients, ζ avg (n), can be found by taking the IDFT of ζ avg (k).
This transmission coefficient is only valid for small values of ka. At high frequencies, the transmission coefficients of all apertures approach the geometrical limit η = 1 [31]. The transition occurs around ka ≈ 5 [31]. Therefore, we enforce η avg (k) = 1 at ka = 5. To smooth out the discontinuity at ka = 5, an FIR filter of order L = 4aF s c with the magnitude response of η avg (k) can be designed using the Parks McClellan (PM) algorithm [32].
The time-domain filters and their magnitude responses for F s = 44.1 kHz and varying aperture sizes are shown in Fig. 3. The frequency response has a band-stop shape. The band edges are determined by the aperture size. Lower frequency transmission coefficients for larger aperture sizes have values greater than one but are attenuated for small apertures. For a slit, the transmission coefficient values at low frequencies can approach infinity; for other shapes, the values are limited [31]. Nevertheless, this model can give unusually large transmission coefficient values for larger aperture sizes at lower frequencies, leading to instabilities. In fact, (21) is valid only for a small value of ka [30]. Hence, the values of ζ avg (k) should be saturated to a maximum of 0 dB to avoid amplification in the feedback path.

B. Asymmetric Transmission
In coupled room systems, if the connecting aperture size is much smaller compared to the room surface area, then most of the RIR's energy will not be incident on the aperture; hence, only a fraction of the energy will pass through to the coupled room. A scaling factor must be applied to the diffraction filter to mimic this energy exchange. This is motivated by the fact that one of the most prominent perceptual effects of closing a door between two rooms is the fall in sound pressure level.
The paraunitary matrices designed so far are symmetric, i.e., the coupling filter between rooms i and j, is the same as that between rooms j and i. Theory by Cremer and Müller [33] suggests that the coupling factors between the two rooms are not the same and depend on the relative absorption area in each room. For 2 coupled rooms, let S 1 and S 2 be the total surface areas of two rooms with absorption coefficients α 1 , α 2 and a be the aperture area connecting the two rooms. The coupling factors are defined as [34] Here, k 1 indicates the coupling factor for energy going from room 1 to room 2, and k 2 is the coupling factor for energy going from room 2 to room 1. The mean coupling factor is defined For simplicity, we denote the denominator (total absorbing area excluding the aperture, plus the aperture that does not absorb any sound) for each room as S i , so the coupling factor becomes k i = a/S i . It has been shown in [2], that premultiplication of the feedback matrix with a diagonal matrix, and post-multiplication with its inverse preserves the unilossless property, i.e., conjugation with a non-singular diagonal matrix preserves energy. Therefore, for a diagonal matrix, E, and a paraunitary matrix, Φ(z),Φ(z) = E Φ(z) E −1 is unilossless. We use this property to design the appropriate diagonal matrix, E, to introduce asymmetric coupling.
For connecting N rooms, let us assume an N × N diagonal matrix, E = diag(e 1 , . . . , e N ). Therefore, E −1 = diag (1/e 1 , . . . , 1/e N ). The (i, j) th element ofΦ(z) is given byΦ Now, if the diffraction filter of length L ij between rooms i and j is denoted by FIR(L ij ), we want, Note that the FIR filters between rooms i and j are identical to those between rooms j and i because the filter only depends on the connecting aperture size, which is the same between rooms i and j. On multiplying and dividingΦ ij (z) andΦ ji (z), we get, From this, the appropriate scaling factor for each off-diagonal element in the coupling matrix is determined by the ratio of the aperture area and the mean coupling factor. The diffraction filters are scaled such that If we fix any element of the diagonal matrix to be unity, say e 1 = 1, then the calculation of the other elements becomes trivial, e i = S 1 /S i . Therefore, the diagonal matrix is constructed as, This ensures that, (E Φ(z) E −1 ) ij = a S i FIR(L ij ). Hence, the appropriate asymmetric coupling coefficient is maintained while maintaining the unilossless property of the feedback matrix.

V. EVALUATION
In this section, we present two examples of the GFDN with filter feedback matrices. In the first example, we investigate the case of two coupled rooms and model diffraction for different aperture sizes connecting the rooms. In the second example, we explore the scenario of three rooms all coupled with each other. We investigate the design of an appropriate paraunitary coupling matrix that maintains the desired T 60 of the three rooms while introducing filtering effects in the RIR 1 .

A. Two-Room Coupling
We designed a 16 delay line GFDN, with 8 delay lines each representing the smaller and larger room respectively, with the source placed in the smaller room and the listener placed in the bigger room. The relative source and listener locations are determined by the b 1  respectively. M 1 and M 2 are scaled Hadamard matrices for maximum diffusion. The absorption coefficients are calculated according to Sabine's equation [35] using the maximum T 60 values and room dimensions.
Assuming a circular aperture connecting the rooms, the aperture radius between the two rooms is increased in steps. This, along with the room area and absorption coefficients, is used to calculate the coupling coefficients according to (23). The frequency response of the Filter Feedback Matrix's (FFM) coupling matrix is shown in Fig. 4 for increasing aperture sizes. Apart from the obvious band-stop characteristic, the lower frequencies are amplified as the aperture size increases.
For the Scalar Feedback Matrix (SFM), the coupling matrix with asymmetric coupling coefficients is The two-stage decay plots of the coupled RIRs as the aperture size increases are shown in Fig. 5. All RIRs have been generated with the FDNToolbox [36]. As the aperture increases, the early decay generally becomes shorter for the RIRs. This is expected, since the rate of energy exchange between the coupled rooms increases with larger aperture size. In [31], Summers observes that for a larger total coupling area, the turning point, which is the point at which the two slopes intersect, shifts earlier, but the difference in the two slopes reduces. A larger coupling area allows faster energy exchange between the two rooms. The T 60 s approach each other, and the coupled rooms essentially act as a single space with single-slope decay when the aperture is large enough. This behavior can be observed in both the energy decay curves in Fig. 5. The T 60 values of the two slopes and their differences are reported in Table I. The difference in the T 60 s generally reduces with increasing aperture size for both kinds of matrices. The significant difference is in the late-stage T 60 , which is longer for the FFM. The RIRs with the SFM and FFM have different spectral characteristics depending on the aperture size.

B. Multi-Room Coupling
We simulate an example of a 3 coupled room system, in which all rooms are interconnected, with the source in the dryest room and the microphone in the largest and most reverberant room. The dimensions of the rooms are 3 × 1.8 × 2 m 3 , 4 × 3.2 × 3.8 m 3 and 6 × 5.5 × 4.5 m 3 respectively. Each room is modeled by 4 delay lines in the GFDN and has M 1 , M 2 , M 3 equal to scaled Hadamard matrices.
We designed a coupling matrix with an aperture radius of 0.22 m between rooms 1 and 2, 0.35 m between rooms 2 and 3, and 0.5 m between rooms 1 and 3. The diffraction filters (with appropriate scaling based on aperture and room area) are Fig. 8. GFDN example for three coupled rooms. Fig. 8(a) -desired T 60 filters for individual rooms. Fig. 8(b) -T 60 error range for the three rooms due to per-frequency bin Procrustes unitarization. Fig. 8(c) -Spectrograms of the three rooms and coupled space.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
placed on the off-diagonal elements of the polynomial matrix. To determine the diagonal elements of the matrix (self-coupling filters), the energy of each column is matched to one. A paraunitary matrix requires some negative elements -there are 2 6 possible sign combinations. We choose the most optimal sign combination described in Section III-D. Starting from this matrix, the Procrustes problem is solved at every frequency bin to give us a nearly paraunitary matrix with the desired magnitude response. The initial polynomial coupling matrix designed is shown in Fig. 6(a). The matrix given by the per-frequency bin Procrustes solution is shown in Fig. 6(b).
We also tried Tkacenko's method with the degree of the desired FIR PU interpolant as the maximum length of filters in the initial polynomial matrix. Even for a large number of iterations, it failed to give us close fits the desired frequency response, and the singular values did not converge to one. This is likely because no FIR PU interpolant exists for the given matrix.
The magnitude responses are shown in Fig. 7. The perfrequency Procrustes solution matches the desired magnitude response exactly except at positions (1, 2) and (2, 1), where the desired response is significantly attenuated. In Fig. 6(c) and (d), we show the magnitude of the singular values per-frequency bin of the off-diagonal matrix elements with 2x zero-padding. These magnitudes should be unity for a paraunitary matrix. The per-frequency Procrustes solution has its singular value magnitudes within 10 −2 dB.
The three rooms' desired T 60 filters are shown in Fig. 8(a). Room 1 has the shortest T 60 , and Room 3 has the longest. The maximum error in the GFDN's T 60 by using a nearly paraunitary coupling matrix, calculated according to (8), is shown in Fig. 8(b) by solid lines. The dashed lines show the the Just Noticeable Difference (JND) in reverberation time, 5% [37]. These errors are perceptually negligible, and in the worst case, touch the JND but never exceed it. The spectrograms of the individual room RIRs and the resulting coupled room RIR when the source is placed in Room 1, and the receiver is placed in Room 3 are shown in Fig. 8(c). The net decay time is affected by all three rooms. The RIR observed in Room 3 has a noticeable band-stop effect, but the lower frequencies are also attenuated because of the small aperture sizes.

VI. CONCLUSION
Grouped Feedback Delay Networks consist of multiple FDNs connected in parallel, each with its own absorption filter and a block mixing matrix that controls inter-and intra-group coupling. In this article, we have extended the design of the GFDN to include filter feedback matrices (FFM). We have shown that for the resulting mixing matrix to be energy-preserving, the matrix of coupling filters can be designed to be paraunitary. Methods of designing low-order, causal FIR FFMs that are ideal for real-time implementation have been discussed. A closed-form coupling matrix with an FIR coupling filter for coupling 2 FDNs has been derived. Some methods for designing the closest paraunitary (or nearly paraunitary) matrix to a given polynomial matrix have been discussed for coupling K FDNs. We have used FFMs to simulate diffraction in coupled rooms based on a mid-frequency model by Miles. The inclusion of asymmetric coupling coefficients for gain adjustment has been proposed. Numerical examples have been included to exhibit the effects of including filters in the feedback matrix and how they compare against scalar feedback matrices.
Including filters in the feedback matrix allows a wide range of design choices. GFDNs with the proposed filter feedback matrices can be used to model late reverberation for complex coupled spaces in video games and XR applications. The theory related to the design of paraunitary coupling matrices extends beyond FDNs. It is useful for modeling any energy-preserving interaction in coupled systems, such as coupling between an instrument's bridge and its strings, since losses can be easily added to an already stable and lossless system. Further studies should look into diffraction filter design for coupled rooms since models in the literature are only valid for a small frequency range or have stability issues for larger aperture sizes. Future work should also include validating the proposed methods with perceptual and numerical evaluations, against room impulse responses obtained with a physically accurate method that can model diffraction, such as the Finite Difference Time Domain (FDTD) method.
where · * denotes complex conjugation and