Table of Contents
ISRN Signal Processing
Volume 2011, Article ID 971051, 6 pages
http://dx.doi.org/10.5402/2011/971051
Research Article

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

1Department of Informatics, University of Laghouat, BP 37G, Laghouat 03000, Algeria
2Applied DSP and VLSI Research Group, Department of Electronic Systems, University of Westminster, 115 New Cavendish Street, London W1W 6UW, UK
3Department of Electronics, ENPA, BP 182, Algiers 16000, Algeria

Received 3 March 2011; Accepted 21 March 2011

Academic Editors: Y. H. Ha and F. Perez-Cruz

Copyright © 2011 Maamar Ahfir et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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) [13], 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.𝑥,100,𝑐𝑋(𝑘)=𝑥(𝑘1)𝑥(𝑘2)(𝑘𝑛)𝐴=000010000010,𝐵=𝐶=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𝑛00.(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.

971051.fig.001
Figure 1: Acoustic impulse response measured in an audio-conferencing room.
971051.fig.002
Figure 2: Hankel singular values plot of both methods for the FIR sequence corresponding to Figure 6(a).
fig3
Figure 3: IIR filter spectral reproduction obtained by different methods for 𝑘=6 as compared to that of the original FIR prototype of Figure 6(a), (a) Magnitude and (b) Phase.

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).

971051.fig.004
Figure 4: Hankel singular values plot of the standard method for the transformed IIR (delivered by the alternative method for 𝑘=6) to FIR sequence.
fig5
Figure 5: IIR filter spectral reproduction obtained by the standard method for 𝑘=6 as compared to that of the transformed IIR (delivered by the alternative method for 𝑘=6) to FIR sequence—(a) Magnitude, (b) Phase, (c) Magnitude error, and (d) Phase error.
fig6
Figure 6: 𝑧-plane plot of the delivered IIR filter by the standard method for the selected order 𝑘=6 in (b) as compared to that of the original FIR prototype in (a).

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.

References

  1. B. D. Radlovic and R. A. Kennedy, “Nonminimum-phase equalization and its subjective importance in room acoustics,” IEEE Transactions on Speech and Audio Processing, vol. 8, no. 6, pp. 728–737, 2000. View at Publisher · View at Google Scholar · View at Scopus
  2. B. D. Radlovic and R. A. Kennedy, “Iterative cepstrum-based approach for speech dereverberation,” in Proceedings of the ISSPA, Brisbane, Australia, August 1999.
  3. M. Ahfir, I. Kale, and B. Daoud, “Room equalization based on iterative simple complex smoothing of acoustic impulse responses,” in Proceedings of the IEEE ICASSP, vol. 5, pp. 109–112, Toulouse, France, May 2006.
  4. B. Beliczynski, I. Kale, and G. D. Cain, “Approximation of FIR by IIR digital filters: an algorithm based on balanced model reduction,” IEEE Transactions on Signal Processing, vol. 4, no. 3, pp. 532–542, 1992. View at Google Scholar · View at Scopus
  5. V. N. Faddeeva, Computational Method of Linear Algebra, Dover, New York, NY, USA, 1959.
  6. M. Bellanger, Traitement Numerique du Signal, Dunod, Paris, France, 1998.
  7. M. Kunt, Traitement Numerique des Signaux, Dunod, Paris, France, 1981.
  8. J. P. Jullien, A. Gilloire, and A. Saliou, “Mesure des réponses impulsionnelles en acoustique,” Note technique NT/LAA/TSS/181, 1984. View at Google Scholar
  9. I. Kale, G. D. Cain, and R. C. S. Morling, “Minimum-phase filter design from linear-phase startpoint via balanced model truncation,” Electronics Letters, vol. 31, no. 20, pp. 1728–1729, 1995. View at Publisher · View at Google Scholar · View at Scopus