# An Alternative Approach to the Balanced Model Truncation Algorithm for Acoustic Minimum-Phase Inverse Filters Order Reduction

Authors: Maamar Ahfir; Izzet Kale; Daoud Berkani
Journal: ISRN Signal Processing
 (2011)

---

## Abstract

We propose an alternative approach to the Balanced Model Truncation method (standard method). This approach reduces substantially the order of minimum-phase inverse filters for equalizing room acoustics. This method is based on a property of the filterz transform function, which modifies the corresponding FIR coefficients before the application of the standard technique to the modified FIR coefficients filter version. In the standard technique, the Hankel singular values plot is the chief guide for a user for the selection of a reduced filter order. Results for minimum-phase inverse filter corresponding to partial equalization of measured acoustic impulse response show the superiority of the proposed method over the standard technique, in terms of reduced filters order selection.


---

## Body

## 1. Introduction

Although partial equalization in room acoustics allows shorter length FIR digital inverse filters than that required in case of complete equalization (perfect dereverberation) [1–3], these are always of very high order (significantly 1000 coefficients or more for an audio-conferencing room), which makes their DSP implementation more complex and inefficient. Such FIR filters have a minimum-phase character, which makes them amenable for possible order reduction when converting them to IIR filters [4]. One very appealing method of FIR to IIR filter conversion is the Balanced Model Truncation technique presented in [4], which is called here the standard method. In this method, the only allowed way for the selection of a possible reduced filter order depends on the visual Hankel singular values plot, which is the chief guide for a user. The smallest they are the smallest is the filter order. In this paper, we propose an alternative approach to the standard method, which is based on a property of the FIR filter z transform function, where corresponding coefficients are modified before the application of the standard method. The goal is to make corresponding Hankel singular values much smaller to substantially reduce the filter order. The standard technique which converts a given single-input/single-output FIR filter to corresponding IIR approximations can be summarized as follows: convert the corresponding transfer function into a matrices system (state-space form), create a Hankel matrix by calculating the Markov parameters, calculate the Hankel singular values, plot and inspect the Hankel singular values and decide upon the reduced system order (equivalently the rank of the Hankel matrix), and finally convert the reduced order state-space form into a transfer function in order to compare the original and the reduced order system against various criteria, and a new order of approximation is attempted if necessary. The paper is organized as follows: Section 2 presents the FIR to IIR filter conversion technique based on Balanced Model Truncation (standard method). In Section 3, an alternative method based on a well-known property of the z transform of a digital filter is proposed. Results and comments are given in Section 4, and finally Section 5 concludes the paper.


## 2. Standard Method

### 2.1. FIR to IIR Filter Conversion Concept

Let be anN coefficients (ci) FIR digital filter where its transfer function in z domain is written as follows: (1)F(z)=c0+c1z-1+⋯+cnz-n,(2)F(z)=c0+F1(z).
This system with the order n=N-1 can be represented in discrete time domain by the following difference equations [4]:(3)y(k)=CX(k)+Du(k),X(k+1)=AX(k)+Bu(k),
where x(k) and y(k) represent, respectively, the input signal and the output signal of the system and u(k) represents the delta function.(4)X(k)=[x(k-1)x(k-2)⋮x(k-n)],A=[00⋯0010⋯00⋯00⋯10],B=[10⋮0],C=[c1c2⋯cn],D=c0.
Notice that F1(z) is representable by the A, B, and C matrices alone.Now applying thez transform to the system (3),                                                    (5)F1(z)=CBzI-A=z-1CBI-z-1A.
The inverse z transform of F1(z) gives the following infinite impulse response sequence (IIR):(6)Hk=CAk-1Bfork=1,2,3,…Hk represents the Markov parameters for the Hankel matrix which can be defined as follows:(7)H=[H1H2H3⋯H2H3H4⋯H3H4H5⋯⋮].
Now for the (A, B, C) system having a finite impulse response (FIR), the rows and columns of zeros are omitted and the following finite Hankel matrix is used: (8)H=[c1c2⋯cnc2c3⋯0⋮⋮⋮cn0⋯0].
Given such real symmetric matrix (8), the following properties are well known in the literature. (a)
H can be factorized as follows:
(9)H=V-1∑V,
whereV-1V=I, with∑ being a diagonal matrix containing the real eigenvaluesσi ofH,V an orthogonal matrix containing the corresponding eigenvectorsvi, andI a unit matrix,(b)
if the matrix∑ contains the eigenvalues σi of H, then ∑2 contains the eigenvalues σi2 of H2,(c)
the rank ofH is equal to n,(d)
the singular values ofH are the positive square roots of the eigenvalues of H2; then the absolute values of its eigenvalues are given by (10)
(10)|∑|=diag⁡(|σ1|,|σ2|,…,|σn|).


### 2.2. Hankel Singular Values and System Order Reduction

The primary principle used for a system order reduction lies in the elimination of a subsystem associated with the smallest Hankel singular values [4]. Since the rank of H is equal to n, so the same degree (order) of F1(z) which represents the system (A, B, C) then,  reducing the rank n to an order k by applying a partition to the system according to the following partitioning form (11), results in a truncated Hankel matrix with a rank k: (11)|∑|=[|∑k|00|∑n-k|]associatedwithV=[VkVn-k],
where |∑k|=diag(|σ1|,|σ2|,…,|σk|) with the associated n×k rectangular matrix Vk corresponding to the truncated system (Hk) and |∑n-k|=diag(|σk+1|,|σ2|,…,|σn|) with the associated n×n-k rectangular matrix Vn-k corresponding to the rejected system (Hn-k), which is associated with the smallest Hankel singular values.The truncated Hankel matrix with a rankk is equivalent to a truncated system (Ak, Bk, Ck) whose transfer function is given by(12)F1,k(z)=CkBkzIk-Ak,
where        (13)Ak=VT(2:n,1:k)V(1:n-1,1:k),Bk=V(1,1:k)T,Ck=CV(1:n,1:k),
with V(i:j,k:m) being an extraction of the matrix V’s row from i to j and its columns from k to m.Now havingF1,k(z), the transfer function of the truncated system (Ak, Bk, Ck, D), with D=c0, can be written as follows:(14)Fk(z)=c0+F1,k(z),
and the corresponding 2k+1 implementable IIR filter coefficients b0, (ai,bi), i=1,k, determined after the conversion from the system state-space representation into the transfer function form in z domain are given by [5](15)Fk(z)=b0+b1z-1+b2z-2+⋯+bkz-k1+a1z-1+a2z-2+⋯+akz-k,
with k  being the reduced filter order.The original system whose transfer function isF(z) and the reduced order system whose transfer function is Fk(z) can then be compared against various criteria, particularly the spectrum magnitude and phase. The plot of Hankel singular values plays the key role in the selection of the filter order reduction.The standard algorithm can be implemented as follows [4].GivenN coefficients ci of an FIR filter(1)
create the Hankel matrixH, (8),(2)
decompose the matrixH to obtain V and ∑, (9),(3)
display the Hankel singular values, (10), and choose a desired order of approximation,(4)
calculate the matricesAk, Bk, and Ck, (13),(5)
convert the system state-space representation (Ak, Bk, Ck, D), where D=c0, into the transfer function form (15) to compare it with that of the original FIR filter (1).


## 2.1. FIR to IIR Filter Conversion Concept

Let be anN coefficients (ci) FIR digital filter where its transfer function in z domain is written as follows: (1)F(z)=c0+c1z-1+⋯+cnz-n,(2)F(z)=c0+F1(z).
This system with the order n=N-1 can be represented in discrete time domain by the following difference equations [4]:(3)y(k)=CX(k)+Du(k),X(k+1)=AX(k)+Bu(k),
where x(k) and y(k) represent, respectively, the input signal and the output signal of the system and u(k) represents the delta function.(4)X(k)=[x(k-1)x(k-2)⋮x(k-n)],A=[00⋯0010⋯00⋯00⋯10],B=[10⋮0],C=[c1c2⋯cn],D=c0.
Notice that F1(z) is representable by the A, B, and C matrices alone.Now applying thez transform to the system (3),                                                    (5)F1(z)=CBzI-A=z-1CBI-z-1A.
The inverse z transform of F1(z) gives the following infinite impulse response sequence (IIR):(6)Hk=CAk-1Bfork=1,2,3,…Hk represents the Markov parameters for the Hankel matrix which can be defined as follows:(7)H=[H1H2H3⋯H2H3H4⋯H3H4H5⋯⋮].
Now for the (A, B, C) system having a finite impulse response (FIR), the rows and columns of zeros are omitted and the following finite Hankel matrix is used: (8)H=[c1c2⋯cnc2c3⋯0⋮⋮⋮cn0⋯0].
Given such real symmetric matrix (8), the following properties are well known in the literature. (a)
H can be factorized as follows:
(9)H=V-1∑V,
whereV-1V=I, with∑ being a diagonal matrix containing the real eigenvaluesσi ofH,V an orthogonal matrix containing the corresponding eigenvectorsvi, andI a unit matrix,(b)
if the matrix∑ contains the eigenvalues σi of H, then ∑2 contains the eigenvalues σi2 of H2,(c)
the rank ofH is equal to n,(d)
the singular values ofH are the positive square roots of the eigenvalues of H2; then the absolute values of its eigenvalues are given by (10)
(10)|∑|=diag⁡(|σ1|,|σ2|,…,|σn|).


## 2.2. Hankel Singular Values and System Order Reduction

The primary principle used for a system order reduction lies in the elimination of a subsystem associated with the smallest Hankel singular values [4]. Since the rank of H is equal to n, so the same degree (order) of F1(z) which represents the system (A, B, C) then,  reducing the rank n to an order k by applying a partition to the system according to the following partitioning form (11), results in a truncated Hankel matrix with a rank k: (11)|∑|=[|∑k|00|∑n-k|]associatedwithV=[VkVn-k],
where |∑k|=diag(|σ1|,|σ2|,…,|σk|) with the associated n×k rectangular matrix Vk corresponding to the truncated system (Hk) and |∑n-k|=diag(|σk+1|,|σ2|,…,|σn|) with the associated n×n-k rectangular matrix Vn-k corresponding to the rejected system (Hn-k), which is associated with the smallest Hankel singular values.The truncated Hankel matrix with a rankk is equivalent to a truncated system (Ak, Bk, Ck) whose transfer function is given by(12)F1,k(z)=CkBkzIk-Ak,
where        (13)Ak=VT(2:n,1:k)V(1:n-1,1:k),Bk=V(1,1:k)T,Ck=CV(1:n,1:k),
with V(i:j,k:m) being an extraction of the matrix V’s row from i to j and its columns from k to m.Now havingF1,k(z), the transfer function of the truncated system (Ak, Bk, Ck, D), with D=c0, can be written as follows:(14)Fk(z)=c0+F1,k(z),
and the corresponding 2k+1 implementable IIR filter coefficients b0, (ai,bi), i=1,k, determined after the conversion from the system state-space representation into the transfer function form in z domain are given by [5](15)Fk(z)=b0+b1z-1+b2z-2+⋯+bkz-k1+a1z-1+a2z-2+⋯+akz-k,
with k  being the reduced filter order.The original system whose transfer function isF(z) and the reduced order system whose transfer function is Fk(z) can then be compared against various criteria, particularly the spectrum magnitude and phase. The plot of Hankel singular values plays the key role in the selection of the filter order reduction.The standard algorithm can be implemented as follows [4].GivenN coefficients ci of an FIR filter(1)
create the Hankel matrixH, (8),(2)
decompose the matrixH to obtain V and ∑, (9),(3)
display the Hankel singular values, (10), and choose a desired order of approximation,(4)
calculate the matricesAk, Bk, and Ck, (13),(5)
convert the system state-space representation (Ak, Bk, Ck, D), where D=c0, into the transfer function form (15) to compare it with that of the original FIR filter (1).


## 3. Alternative Method

In this method, we use a well-known property of thez transform of a filter, being an FIR or IIR, in order to modify the original filter coefficients, that is [6, 7],(16)kg(k)⟺-zdG(z)dz,
with g(k) being the sequence of the inverse z transform of G(z).The transfer functionF(z) of an FIR filter can be written as follows: (17)F(z)=c0+c1z-1+⋯+cnz-n=c0+F1(z)=c0-zddz(c1z-1+c22z-2+c33z-3+⋯+cnnz-n),F(z)=c0-zddzG(z),
with  (18)G(z)=c1z-1+c22z-2+c33z-3+⋯+cnnz-n,
or (19)G(z)=∫z-1F1(z)dz.Alternatively, the standard algorithm stated in Section2 can be applied to the transfer function G(z), ((18), (19)), which represents the z transform function of the modified coefficients of the original FIR filter, equivalently the average of the delayed original sequence ci, i=1,n, in z domain. The corresponding reduced order transfer function, which is the output of the algorithm, will be Gk(z), with k being the reduced order which will be selected from the Hankel singular values plot. The reduced order system whose transfer function is Fk(z) is given by             (20)Fk(z)=c0-zddzGk(z).
In time domain, this is equivalent to(21)h(k)=c0u(k)+kg(k),
with g(k) being the inverse z transform of Gk(z).The transfer function of the corresponding implementable IIR filter is now given by(22)Fk(z)=c0-zddz(b0+b1z-1+b2z-2+⋯+bkz-k1+a1z-1+a2z-2+⋯+akz-k),
and the implementation complexity is 2(2k+1)-1 coefficients, because of the derivative function.


## 4. Results

Both techniques were applied to a partial equalization FIR digital filter for an acoustic impulse response, representing the context of an audio-conferencing room (Figure1) [8]. The corresponding minimum-phase inverse impulse response (FIR prototype) was obtained partially by the Homomorphic method for an optimum L=16. This particular choice of  L  is for avoiding the undesirable effects of the circular deconvolution of the inverse impulse response with the original room impulse response, for some predefined threshold level of magnitude equalization, and which provides some desired sound quality [2]. Figure 6(a) shows its z-plane plot. The plot of Hankel singular values of both methods depicted in Figure 2 shows clearly that the alternative method allows the choice of an IIR filter order which is much more reduced than that allowed when using the standard method. When comparing the spectrum magnitude and phase obtained by both methods for the same reduced IIR filter order (with a selected order k=6 and a root mean square error of some 0.04) to that of the original FIR prototype, we can see in Figure 3 a better spectral reproduction obtained by the alternative method than that obtained by the standard method, that is, the smoothed form against the affected form (peaks), which is the result of the truncating effect caused by the selection of an underestimated IIR filter order, as displayed by the corresponding Hankel singular values plot in Figure 2.Figure 1
Acoustic impulse response measured in an audio-conferencing room.Figure 2
Hankel singular values plot of both methods for the FIR sequence corresponding to Figure6(a).IIR filter spectral reproduction obtained by different methods fork=6 as compared to that of the original FIR prototype of Figure 6(a), (a) Magnitude and (b) Phase.
(a)(b)Both methods were also compared to the direct average of the delayed original FIR sequence (ci, i=1,n) in z domain, making sense of noise reduction, that is c0+G(z) for the full order. We can see in Figure 3 that despite the fact that the spectral average has completely removed the noise from the original FIR spectrum, which resulted in better smoothed spectral form, the proposed alternative method shows its superiority in tracking nearly closely the form of the original spectrum of the FIR prototype, for both magnitude and phase.To overcome the complexity implementation problem when using the alternative method, (2(2k+1)-1 coefficients are needed rather than (2k+1) coefficients in the standard method); the standard method has been applied to the transformed version of the reduced order IIR impulse response to an FIR filter (infinite impulse response truncated to a finite number of samples) [9]. The Hankel singular values displayed in Figure 4 allow choosing a new reduced IIR filter order being the same as that delivered by the alternative method, that is, k=6. In Figure 5, we can see a perfect spectral reproduction obtained by the new IIR filter for the reduced order k=6, as compared to that of the transformed IIR to FIR prototype, but with less implementation complexity, that is, (2k+1) coefficients rather than 2(2k+1)-1 coefficients for the alternative method. This means that the relation (15) is used for DSP implementation rather than (22).Figure 4
Hankel singular values plot of the standard method for the transformed IIR (delivered by the alternative method fork=6) to FIR sequence.IIR filter spectral reproduction obtained by the standard method fork=6 as compared to that of the transformed IIR (delivered by the alternative method for k=6) to FIR sequence—(a) Magnitude, (b) Phase, (c) Magnitude error, and (d) Phase error.
(a)(b)(c)(d)z-plane plot of the delivered IIR filter by the standard method for the selected order k=6 in (b) as compared to that of the original FIR prototype in (a).
(a)(b)The corresponding magnitude and phase error are both less than −70 dB in the entire frequency band (Figures5(c) and 5(d)). Figure 6(b) shows the z-plane plot for this approximation IIR filter, with only 6 zeros and 6 poles (13 IIR coefficients) rather than more than 1000 zeros (1000 FIR coefficients) for the original FIR prototype (Figure 6(a)).


## 5. Conclusion

An alternative approach to the Balanced Model Truncation algorithm (standard method) to substantially reduce minimum-phase inverse filters order for room acoustics equalization is proposed. This method is based on a property of the FIR filterz transform function, which modifies the original coefficients of a given FIR sequence, by taking the average of its corresponding transfer function before the conversion to its IIR coefficients approximations. Results obtained by the alternative method show better specifications reproduction of the FIR prototype, with a very low IIR filter order, as compared to that obtained by the standard method.



---

*Source: 971051-2011-04-13.xml*