Partial Equalization of Non-Minimum-Phase Impulse Responses

We propose a modiﬁed version of the standard homomorphic method to design a minimum-phase inverse ﬁlter for non-mini-mum-phase impulse responses equalization. In the proposed approach some of the dominant poles of the ﬁlter transfer function are replaced by new ones before carrying out the inverse DFT. This method is useful when partial magnitude equalization is intended. Results for an impulse response measured in the car interior show that by using the modiﬁed version we can control the sound quality more precisely than when using the standard method.


INTRODUCTION
In sound-reproduction systems an equalization filter is often used to modify the frequency spectrum of the original source before feeding it to the loudspeaker. The purpose is to make the impulse response of the equalized sound-reproduction chain as close as possible to the desired one [1]. In principle the direct inversion of mixed-phase (or non-minimumphase) measured impulse responses of the systems is not possible since it leads to unstable equalization filter realizations. Since any mixed-phase impulse response can be represented mathematically by the convolution of a minimum-phase sequence and a maximum-phase (or all-pass) sequence [2], it is possible to derive and implement an approximate and stable inverse filter for such systems [3]. This is because a causal and stable sequence can invert the minimum-phase component of any mixed-phase sequence and an infinite acausal (anticipatory) and stable sequence can similarly invert the maximum-phase component of such sequences [3]. For the reason of the implementation complexity of such combined equalization filters as it will be discussed in Section 2, the work presented in this paper focuses on the equalization of the minimum-phase component of the system and its partial equalization importance. One method to design such a minimum-phase equalization filter is the homomorphic one based on the measured impulse response of the system. This method known as standard used for the case of single-point equalization is described in Section 2. In Section 3, a modified version of the standard homomorphic method is proposed. It takes into account that the listener is able to detect gradual response variations of less than 0.5 dB [4,5] and hence is able to control the sound quality more accurately. Section 4 shows the magnitude equalization performance results for an impulse response measured in a car interior using both objective and subjective measurements.

STANDARD HOMOMORPHIC METHOD
A non-minimum-phase discrete impulse response, h(n), of a system can be described as [2] where ⊗ denotes the discrete convolution. This can be shown in the frequency domain as where h mp (n) is a minimum-phase sequence, such that its DFT, H mp (k), satisfies the relation where N is the length of h(n) and h ap (n) is an all-pass sequence of |H ap (k)| = 1, for k = 0, 1, . . . , N − 1.
The convolution operation of h mp (n) and h ap (n) can be expressed as the algebraic addition of their corresponding complex cepstra h mp (n) and h ap (n) by the homomorphic transformation [6]. This leads to a decomposition of a non-minimum-phase impulse response into its minimumphase and all-pass components. The standard homomorphic method algorithm is outlined as follows [4,7,8].
(3) Compute the real part of the complex cepstrum of h(n), for n = 0, 1, . . . , N − 1. (4) Compute the corresponding real cepstrum of the minimum-phase h mp (n), where L is a positive real parameter [8].
(7) Compute the equalized response, H eq (k), where G mp (k) represents the inverse of H mp (k), In the time domain, this is equivalent to a deconvolution, with g mp (n), being the inverse DFT of G mp (k).
In the case of L = 1, the algorithm corresponds to a magnitude equalization. If a sufficiently large number N is used for DFT computation, the effect of magnitude distortion caused by the system can be perfectly removed in practice by convolving h(n) with the inverse minimum-phase impulse response g mp (n) [3,7]. The effect of phase distortion can also be solved by convolving the all-pass sequence, h ap (n), (obtained after deconvolution of h(n) with g mp (n)) with its time reversed version, h ap (−n), [4,9]. As a result, implementation of such combined equalization (complete equalization) requires very long FIR filters. But this is not always required in practice. For this reason, the equalization of the all-pass component (phase equalization) will not be considered in this work.
In the case of L > 1, the algorithm corresponds to a partial magnitude equalization. This requires a shorter FIR filter to keep the phase distortion below the threshold of audibility [4].
Sometimes, the frequency response of the system, H(k), and hence its minimum-phase part, H mp (k), can be represented by a low number of isolated dominant zeros. In such a case, increasing the parameter L by a significant value during the control process may shorten the length of the equalization filter too, resulting in an unsatisfactory equalization performance. This is because an increase in L results in a decrease of all the radii of the complex poles of G mp (k) together according to the relation derived from (3), (7), and (9) (see the appendix), This means that the complex poles of G mp (k) appear to be pushed together towards the origin of the unit circle.
In the next section, we propose an alternative approach in which, instead of pushing all the poles of G mp (k), we push the most dominant of them selectively and slightly towards the origin of the unit circle by decreasing the corresponding high values of the Q factors (values of the steady-state resonances). This allows controlling the magnitude equalization performance more precisely-especially applicable in practice. The reason is that the listener is able to detect gradual response variations of less than 0.5 dB [4,5]. Furthermore, the proposed technique is advantageous when the parameter L cannot be calculated theoretically, for example, for the case of the direct inverse filtering (no cepstral analysis, i.e., no steps 2 to 6) of a small reverberant room where the dominant poles can be identified even if they are closely spaced [10,11].

MODIFIED VERSION
A replacing method of some dominant poles of the inverse minimum-phase function G mp (k) is described in this section. They are identified using the standard method (L = 1) Ahfir Maamar et al. 3 and then replaced before doing the inverse DFT in order to calculate the corresponding discrete time sequence g mp (n) representing the impulse response of the new equalization filter.
The z transform function of a complex pole pair is expressed as [10,12] where |a| is the pole radius in the z plane and θ = 2π( f p / f e ) is its phase angle with f e being the sampling frequency and f p the frequency of the complex pole.
Taking the inverse z transform of H p (z) the corresponding impulse response is [10,12] where u(n) is a unit step function. The transfer function of the selective filter for a complex pole pair is where the transfer function H p (z) contains a new complex pole pair at the same frequency of the old pair but at a desired smaller radius, | a|. This technique allows us to decrease selectively the Q factors values of a low order of isolated pole pairs in the frequency response G mp (k). The new inverse minimum-phase function becomes and its discrete version where P is the number of identified and replaced dominant pole pairs from G mp (k), and H (P) s (k), p = 1, . . . , P, are the sampled frequency responses of selective filters equal to the number of replaced pole pairs, P.
This function is then inverted using the inverse DFT in order to obtain its discrete time domain equivalent g mp (n), shorter than g mp (n) calculated by standard method for L = 1 and longer than g mp (n) obtained for L > 1.
One method to identify frequencies of the isolated poles is to iteratively search for the increased magnitude response level of G mp (k) caused by poles (peaks) residing within the frequency range of interest (in our case below 4 kHz). In each iteration a maximum magnitude level G mp ( f p ) corresponding to the highest pole frequency f p is found. This technique was found robust even in the case of very closely spaced poles [10,11]. After determining the frequency f p of the highest pole, the corresponding pole radius must be determined based on the Q factor value according to the following relation, since our work here is restricted to a low order of isolated poles [10,[13][14][15], The replacing method means that the dominant poles of G mp (k) are identified one by one and then replaced iteratively by new ones, where each corresponds to a desired Q factor, Q = 1/(1 − | a|), starting from the most dominant one.
The implementation algorithm of the proposed modified method (useful for partial magnitude equalization) is as follows.
(1) Compute the steps 1 to 6 as in the standard method for using G mp (k) = G mp (k).
In the time domain, this is equivalent to a deconvolution In the next section we present the performance evaluation of the magnitude equalization performed by the proposed version as compared to that from standard method, using both objective measures based on an error criterion and subjective tests of speech quality.

RESULTS
In order to assess the performance of our algorithms, we used a frequency domain error criterion, which estimates the standard deviation of the magnitude response from a constant level [4]. The error criterion Δ(dB) is given as follows: where 10 log 10 H eq (k) .
Two examples of non-minimum-phase impulse responses were used to compare both algorithms. The first one was synthetic, used just to enlighten the replacing approach of a low order of isolated poles, and the second one used real measurements taken in a car interior. We also introduced in the proposed version a real parameter l (l > 1), in order to selectively decrease the highest Q factors of dominant poles. This means that the new replacement poles correspond to desired Q factors, Q = Q/l.

Synthetic impulse response
We first considered a simple synthetic impulse response. This was obtained by successive convolution of six known zeros sequences somewhat isolated in the complex z plane ( Figure 1) and with at least one placed outside the unit circle, which make the impulse response a non-minimum-phase one. These are defined as ( f e = 8 kHz):  These two (P = 2) poles are selectively replaced by two new poles of smaller radius corresponding to Q 1 /l and Q 2 /l factors, respectively, with a significant value of l, (l = 2) ( Figure 1). In Figure 2, even with some error in the estimation of poles, we still can observe a decrease in the Q factors depending on the position of the new poles. This corresponds to a reduction of g mp (n) length when compared to g mp (n). When using the standard method and considering the same significant value of L (L = l = 2), we can see in Figure 2 that all the poles have been pushed together towards the origin of the unit circle too, resulting also in the reduction of the g mp (n) length, but this is considerably shorter than that of g mp (n).
The evaluation of the objective error criterion for this example is not considered because of its little practical interest.

Practical impulse response
A real impulse response was measured for the car interior at a sampling frequency of 8 kHz. A record of 1024 samples was zero padded up to N = 2048. This impulse response is shown in Figure 3. In Figure 4 an unstable direct inverse impulse response is shown, demonstrating its non-minimumphase character. Figure 5 shows the inverse minimum-phase frequency response G mp (k). It was calculated using the standard method (L = 1). The most dominant pole can be clearly seen there. The search was limited to a single pole, that is, P = 1, such that |a 1 | = 0.9993 at f 1 = 70.38 Hz that caused the inverse Ahfir Maamar et al. minimum-phase impulse response g mp (n) to be of a very long duration (Figure 7).
When using the standard version for a significant value of L, (L = 2), ( Figure 6) we observed that all the poles of G mp (k) were pushed together towards the origin of the unit circle too, resulting in an inverse minimum-phase impulse response g mp (n) (Figure 7) to be reduced in time too.
When using a modified version, in order to gradually reduce the length of the inverse minimum-phase impulse response g mp (n), only the most dominant pole needed to be replaced by a new pole with smaller radius. This pole corresponded to a Q/l factor with the same value of l, (l = L = 2), but at the same frequency. In Figure 5 we can see the inverse minimum-phase frequency response of G mp (k), where only the most dominant pole appears to be pushed towards the origin, with l = 2. Figure 7 also shows the corresponding inverse minimum-phase impulse response of g mp (n). Interestingly, its duration is not reduced here too. This may correspond to a desired magnitude equalization (Figure 8), if the system impulse response was minimum phase (no phase distortion effects). This is because the magnitude spectrum of the second case (modified method, l = 2, Δ(dB) = 0.7) is flatter than that of the first case (standard method, L = 2, Δ(dB) = 2.4). This means less magnitude distortion of the system.

Performance testing
The experiment was performed by developing models in Matlab and Simulink and carrying out listening tests using headphones. A reproduced speech signal of few seconds in duration was generated by filtering a clean speech (male and female measured in anechoic chamber) by the measured impulse response of the car interior ( Figure 3). In order to avoid undesirable convolution effects, we considered a sufficient large number N = 8192 for DFT computations. The reproduced speech signal was then filtered using equalizing filters calculated by the standard method with L = 1 and L = 2  Figure 7) and the modified version (P = 1), respectively. For the latter case, the inverse impulse responses corresponded to each error criterion, function of the parameter l such as that of Figure 7 with l = 2, for example. Test signals were played to ten untrained listeners with normal hearing at a comfortable listening level. The qualitative assessment of the test signals was based on subjective judgment of three listening sessions per each recording scheduled on six consecutive days.
The first signal was always chosen to be clean speech, while the reproduced unequalized and partially equalized speech signals were played in random order. The reproduced speech signals corresponded to the objective error criteria of 5 dB (unequalized signal for l = 0), 0 dB (magnitude equalized signal for l = 1), and 0.3 ≤ Δ(dB) ≤ 3 (partially equalized   speech. The final result was calculated as a mean of the individual listening results (18 each) for each of 10 subjects. The results are shown in Figure 9 as a function of parameter l, ranging from l = 0 (unequalized signal) to l = 5 (partially equalized signal).
The results confirmed those reported in [4] with higher accuracy. The highest score corresponds to the optimal quality of speech. That means no perception of phase distortion (like a bell chime sounded at the background, when l < 3), no echo and less magnitude distortion caused by the system.
The results also show the sensitivity of the listener's ear to small gradual response variations (a variation of less than 7 0.5 dB of objective error corresponds to a significant variation of subjective score of the sound quality); although the participants in the experiment were nonexpert listeners.

CONCLUSIONS
In this paper a modified version of the standard homomorphic method for minimum-phase inverse filter design for non-minimum-phase impulse responses equalization is presented. This version is useful in cases of partial magnitude equalization, where the dominant zeros density of the system is not very high. Although it is used in this work as an additional optimizing tool for the psychoacoustic quality measurement of speech, this alternative approach is advantageous in case of the direct inverse filtering (minimumphase system) when perfect equalization of a small reverberant room is not desired.

APPENDIX
Proof for the relation (13).
The real part of the complex cepstrum of h(n) in (6)

ACKNOWLEDGMENT
The authors would like to thank the reviewers for their helpful comments and suggestions.