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 filter z 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.

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 ๐‘ง 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 ๐‘ง 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 an ๐‘ coefficients (๐‘๐‘–) FIR digital filter where its transfer function in ๐‘ง domain is written as follows: ๐น(๐‘ง)=๐‘0+๐‘1๐‘งโˆ’1+โ‹ฏ+๐‘๐‘›๐‘งโˆ’๐‘›,(1)๐น(๐‘ง)=๐‘0+๐น1(๐‘ง).(2) This system with the order ๐‘›=๐‘โˆ’1 can be represented in discrete time domain by the following difference equations [4]:๐‘‹๐‘ฆ(๐‘˜)=๐ถ๐‘‹(๐‘˜)+๐ท๐‘ข(๐‘˜),(๐‘˜+1)=๐ด๐‘‹(๐‘˜)+๐ต๐‘ข(๐‘˜),(3) where ๐‘ฅ(๐‘˜) and ๐‘ฆ(๐‘˜) represent, respectively, the input signal and the output signal of the system and ๐‘ข(๐‘˜) represents the delta function.โŽกโŽขโŽขโŽขโŽฃโ‹ฎ๐‘ฅโŽคโŽฅโŽฅโŽฅโŽฆ,โŽกโŽขโŽขโŽขโŽฃโ‹ฏโŽคโŽฅโŽฅโŽฅโŽฆโŽกโŽขโŽขโŽขโŽฃ10โ‹ฎ0โŽคโŽฅโŽฅโŽฅโŽฆ,๎€บ๐‘๐‘‹(๐‘˜)=๐‘ฅ(๐‘˜โˆ’1)๐‘ฅ(๐‘˜โˆ’2)(๐‘˜โˆ’๐‘›)๐ด=00โ‹ฏ0010โ‹ฏ0000โ‹ฏ10,๐ต=๐ถ=1๐‘2โ‹ฏ๐‘๐‘›๎€ป,๐ท=๐‘0.(4) Notice that ๐น1(๐‘ง) is representable by the ๐ด, ๐ต, and ๐ถ matrices alone.

Now applying the ๐‘ง transform to the system (3), ๐น1(๐‘ง)=๐ถ๐ต=๐‘ง๐‘ง๐ผโˆ’๐ดโˆ’1๐ถ๐ต๐ผโˆ’๐‘งโˆ’1๐ด.(5) The inverse ๐‘ง transform of ๐น1(๐‘ง) gives the following infinite impulse response sequence (IIR):๐ป๐‘˜=๐ถ๐ด๐‘˜โˆ’1๐ตfor๐‘˜=1,2,3,โ€ฆ(6)๐ป๐‘˜ represents the Markov parameters for the Hankel matrix which can be defined as follows:โŽกโŽขโŽขโŽขโŽฃ๐ป๐ป=1๐ป2๐ป3โ‹ฏ๐ป2๐ป3๐ป4โ‹ฏ๐ป3๐ป4๐ป5โ‹ฏโ‹ฎโŽคโŽฅโŽฅโŽฅโŽฆ.(7) Now for the (๐ด, ๐ต, ๐ถ) system having a finite impulse response (FIR), the rows and columns of zeros are omitted and the following finite Hankel matrix is used: โŽกโŽขโŽขโŽขโŽฃ๐‘๐ป=1๐‘2โ‹ฏ๐‘๐‘›๐‘2๐‘3๐‘โ‹ฏ0โ‹ฎโ‹ฎโ‹ฎ๐‘›โŽคโŽฅโŽฅโŽฅโŽฆ0โ‹ฏ0.(8) Given such real symmetric matrix (8), the following properties are well known in the literature. (a)๐ป can be factorized as follows: ๐ป=๐‘‰โˆ’1๎“๐‘‰,(9) where๐‘‰โˆ’1๐‘‰=๐ผ, withโˆ‘ being a diagonal matrix containing the real eigenvalues๐œŽ๐‘– of๐ป,๐‘‰ an orthogonal matrix containing the corresponding eigenvectors๐‘ฃ๐‘–, and๐ผ a unit matrix, (b)if the matrix โˆ‘ contains the eigenvalues ๐œŽ๐‘– of ๐ป, then โˆ‘2 contains the eigenvalues ๐œŽ๐‘–2 of ๐ป2,(c)the rank of ๐ป is equal to ๐‘›,(d)the singular values of ๐ป are the positive square roots of the eigenvalues of ๐ป2; then the absolute values of its eigenvalues are given by (10) ||๎“||๎€ท||๐œŽ=diag1||,||๐œŽ2||||๐œŽ,โ€ฆ,๐‘›||๎€ธ.(10)

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 ๐ป is equal to ๐‘›, so the same degree (order) of ๐น1(๐‘ง) which represents the system (๐ด, ๐ต, ๐ถ) then, reducing the rank ๐‘› to an order ๐‘˜ by applying a partition to the system according to the following partitioning form (11), results in a truncated Hankel matrix with a rank ๐‘˜: ||๎“||=โŽกโŽขโŽขโŽขโŽขโŽขโŽขโŽฃ|||||๎“๐‘˜|||||00|||||๎“๐‘›โˆ’๐‘˜|||||โŽคโŽฅโŽฅโŽฅโŽฅโŽฅโŽฅโŽฆ๎€บ๐‘‰associatedwith๐‘‰=๐‘˜๐‘‰๐‘›โˆ’๐‘˜๎€ป,(11) where |โˆ‘๐‘˜|=diag(|๐œŽ1|,|๐œŽ2|,โ€ฆ,|๐œŽ๐‘˜|) with the associated ๐‘›ร—๐‘˜ rectangular matrix ๐‘‰๐‘˜ corresponding to the truncated system (๐ป๐‘˜) and |โˆ‘๐‘›โˆ’๐‘˜|=diag(|๐œŽ๐‘˜+1|,|๐œŽ2|,โ€ฆ,|๐œŽ๐‘›|) with the associated ๐‘›ร—๐‘›โˆ’๐‘˜ rectangular matrix ๐‘‰๐‘›โˆ’๐‘˜ corresponding to the rejected system (๐ป๐‘›โˆ’๐‘˜), which is associated with the smallest Hankel singular values.

The truncated Hankel matrix with a rank ๐‘˜ is equivalent to a truncated system (๐ด๐‘˜, ๐ต๐‘˜, ๐ถ๐‘˜) whose transfer function is given by๐น1,๐‘˜๐ถ(๐‘ง)=๐‘˜๐ต๐‘˜๐‘ง๐ผ๐‘˜โˆ’๐ด๐‘˜,(12) where ๐ด๐‘˜=๐‘‰๐‘‡๐ต(2โˆถ๐‘›,1โˆถ๐‘˜)๐‘‰(1โˆถ๐‘›โˆ’1,1โˆถ๐‘˜),๐‘˜=๐‘‰(1,1โˆถ๐‘˜)๐‘‡,๐ถ๐‘˜=๐ถ๐‘‰(1โˆถ๐‘›,1โˆถ๐‘˜),(13) with ๐‘‰(๐‘–โˆถ๐‘—,๐‘˜โˆถ๐‘š) being an extraction of the matrix ๐‘‰โ€™s row from ๐‘– to ๐‘— and its columns from ๐‘˜ to ๐‘š.

Now having ๐น1,๐‘˜(๐‘ง), the transfer function of the truncated system (๐ด๐‘˜, ๐ต๐‘˜, ๐ถ๐‘˜, ๐ท), with ๐ท=๐‘0, can be written as follows:๐น๐‘˜(๐‘ง)=๐‘0+๐น1,๐‘˜(๐‘ง),(14) and the corresponding 2๐‘˜+1 implementable IIR filter coefficients ๐‘0, (๐‘Ž๐‘–,๐‘๐‘–), ๐‘–=1,๐‘˜, determined after the conversion from the system state-space representation into the transfer function form in ๐‘ง domain are given by [5]๐น๐‘˜๐‘(๐‘ง)=0+๐‘1๐‘งโˆ’1+๐‘2๐‘งโˆ’2+โ‹ฏ+๐‘๐‘˜๐‘งโˆ’๐‘˜1+๐‘Ž1๐‘งโˆ’1+๐‘Ž2๐‘งโˆ’2+โ‹ฏ+๐‘Ž๐‘˜๐‘งโˆ’๐‘˜,(15) with ๐‘˜โ€‰โ€‰being the reduced filter order.

The original system whose transfer function is ๐น(๐‘ง) and the reduced order system whose transfer function is ๐น๐‘˜(๐‘ง) 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].

Given ๐‘ coefficients ๐‘๐‘– of an FIR filter(1)create the Hankel matrix ๐ป, (8),(2)decompose the matrix ๐ป to obtain ๐‘‰ and โˆ‘, (9),(3)display the Hankel singular values, (10), and choose a desired order of approximation,(4)calculate the matrices ๐ด๐‘˜, ๐ต๐‘˜, and ๐ถ๐‘˜, (13),(5)convert the system state-space representation (๐ด๐‘˜, ๐ต๐‘˜, ๐ถ๐‘˜, ๐ท), where ๐ท=๐‘0, 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 the ๐‘ง transform of a filter, being an FIR or IIR, in order to modify the original filter coefficients, that is [6, 7],๐‘˜๐‘”(๐‘˜)โŸบโˆ’๐‘ง๐‘‘๐บ(๐‘ง)๐‘‘๐‘ง,(16) with ๐‘”(๐‘˜) being the sequence of the inverse ๐‘ง transform of ๐บ(๐‘ง).

The transfer function ๐น(๐‘ง) of an FIR filter can be written as follows: ๐น(๐‘ง)=๐‘0+๐‘1๐‘งโˆ’1+โ‹ฏ+๐‘๐‘›๐‘งโˆ’๐‘›=๐‘0+๐น1(๐‘ง)=๐‘0๐‘‘โˆ’๐‘ง๎‚€๐‘๐‘‘๐‘ง1๐‘งโˆ’1+๐‘22๐‘งโˆ’2+๐‘33๐‘งโˆ’3๐‘+โ‹ฏ+๐‘›๐‘›๐‘งโˆ’๐‘›๎‚,๐น(๐‘ง)=๐‘0๐‘‘โˆ’๐‘ง๐‘‘๐‘ง๐บ(๐‘ง),(17) with ๐บ(๐‘ง)=๐‘1๐‘งโˆ’1+๐‘22๐‘งโˆ’2+๐‘33๐‘งโˆ’3๐‘+โ‹ฏ+๐‘›๐‘›๐‘งโˆ’๐‘›,(18) or ๎€œ๐‘ง๐บ(๐‘ง)=โˆ’1๐น1(๐‘ง)๐‘‘๐‘ง.(19)

Alternatively, the standard algorithm stated in Section 2 can be applied to the transfer function ๐บ(๐‘ง), ((18), (19)), which represents the ๐‘ง transform function of the modified coefficients of the original FIR filter, equivalently the average of the delayed original sequence ๐‘๐‘–, ๐‘–=1,๐‘›, in ๐‘ง domain. The corresponding reduced order transfer function, which is the output of the algorithm, will be ๐บ๐‘˜(๐‘ง), with ๐‘˜ being the reduced order which will be selected from the Hankel singular values plot. The reduced order system whose transfer function is ๐น๐‘˜(๐‘ง) is given by ๐น๐‘˜(๐‘ง)=๐‘0๐‘‘โˆ’๐‘ง๐บ๐‘‘๐‘ง๐‘˜(๐‘ง).(20) In time domain, this is equivalent toโ„Ž(๐‘˜)=๐‘0๐‘ข(๐‘˜)+๐‘˜๐‘”(๐‘˜),(21) with ๐‘”(๐‘˜) being the inverse ๐‘ง transform of ๐บ๐‘˜(๐‘ง).

The transfer function of the corresponding implementable IIR filter is now given by ๐น๐‘˜(๐‘ง)=๐‘0๐‘‘โˆ’๐‘ง๎‚ต๐‘๐‘‘๐‘ง0+๐‘1๐‘งโˆ’1+๐‘2๐‘งโˆ’2+โ‹ฏ+๐‘๐‘˜๐‘งโˆ’๐‘˜1+๐‘Ž1๐‘งโˆ’1+๐‘Ž2๐‘งโˆ’2+โ‹ฏ+๐‘Ž๐‘˜๐‘งโˆ’๐‘˜๎‚ถ,(22) and the implementation complexity is 2(2๐‘˜+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 (Figure 1) [8]. The corresponding minimum-phase inverse impulse response (FIR prototype) was obtained partially by the Homomorphic method for an optimum ๐ฟ=16. This particular choice ofโ€‰ ๐ฟ โ€‰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 ๐‘ง-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 ๐‘˜=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.

Both methods were also compared to the direct average of the delayed original FIR sequence (๐‘๐‘–, ๐‘–=1,๐‘›) in ๐‘ง domain, making sense of noise reduction, that is ๐‘0+๐บ(๐‘ง) 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(2๐‘˜+1)โˆ’1 coefficients are needed rather than (2๐‘˜+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, ๐‘˜=6. In Figure 5, we can see a perfect spectral reproduction obtained by the new IIR filter for the reduced order ๐‘˜=6, as compared to that of the transformed IIR to FIR prototype, but with less implementation complexity, that is, (2๐‘˜+1) coefficients rather than 2(2๐‘˜+1)โˆ’1 coefficients for the alternative method. This means that the relation (15) is used for DSP implementation rather than (22).

The corresponding magnitude and phase error are both less than โˆ’70โ€‰dB in the entire frequency band (Figures 5(c) and 5(d)). Figure 6(b) shows the ๐‘ง-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 filter ๐‘ง 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.

Acknowledgment

The authors would like to thank the reviewers for their constructive comments.