Design of causal preconditioning filters for adaptive filtering algorithms in real-time multi-channel active noise control applications

Yiming Wang 1

Beijing Ancsonic Technology Co.Ltd 621-1, floor 6, No. 15, Xinxi road, Haidian District, Beijing, 100080, China

Yongjie Zhuang Ray W. Herrick Laboratories, School of Mechanical Engineering, Purdue University 177 S. Russell Street, West Lafayette, IN, 47907-2099, USA

Yangfan Liu Ray W. Herrick Laboratories, School of Mechanical Engineering, Purdue University 177 S. Russell Street, West Lafayette, IN, 47907-2099, USA

ABSTRACT Least-mean-square (LMS) algorithm has been widely used in the area of active noise control (ANC). One practical concern of LMS-based algorithms is its slow convergence rate in multi-channel broadband noise control applications. To improve the convergence speed of traditional LMS algorithm, preconditioning filters were added to the LMS system in previous studies to remove the correlation between the reference signals and decouple the plant responses. However, the preconditioning filters implemented previously are usually obtained through singular value decomposition of the cross spectral matrix of reference signals and the plant response matrix, which does not lead to causal preconditioning filters and, thus, it can only be applied in spatial audio applications when delay is not an important concern, but cannot be implemented in real-time active noise control applications. In the current work, a method is proposed to obtain a casual multi-channel preconditioning filter for real-time active noise control applications through a numerically robust algorithm to perform a spectral factor decomposition to the reference signal cross spectral matrix and a minimum-phase and all-pass decomposition to the multi-channel secondary path. Simulation results show that, by applying the proposed preconditioning filter, the convergence speed of the LMS algorithm can be significantly improved.

1. INTRODUCTION

Least-mean-square (LMS) algorithm is one of the most commonly applied adaptive filtering algorithms. Due to its simplicity and e ff ectiveness, it has been widely used in fields of active noise control (ANC), cross talk cancellation, adaptive equalization, speech analysis, etc [1, 2]. One disadvantage of LMS algorithm for a multiple-input-multiple-output (MIMO) system is its

1 wangyiming.acoustics@gmail.com

inter.noise scorn ROSE DODD wa

slow convergence rate. Several studies have been conducted previously regarding the convergence analysis of MIMO LMS algorithms [3, 4], most of which attribute the main reasons of the slow convergence in ANC applications to two factors: the correlation between di ff erent reference signals and the coupling of di ff erent input and output channel pairs of the secondary plant, i.e., the response to the control sources at error microphones. To improve the convergence rate, Douglas et al. has developed a self-whitening algorithm to preprocess the reference signals with an adaptive filter before feeding them into LMS algorithm [5]. Later, Bai and Elliott have introduced a method to calculate preconditioning filters for the reference signals and the plant responses using singular value decomposition (SVD) [6], which removes the correlation among the reference signals as well as the coupling in the secondary path responses. A limitation of this method is that the preconditioning filters obtained from SVD are not causal filters so that delays must be chosen carefully and applied to those filters to guarantee the causality of the preconditioning filters. This causality issue makes the preconditioning method only applicable in applications where system delay does not a ff ect filter performance, such as spatial audio applications, however, it may become a critical problem for real-time active noise control applications where the noise control performance is quite sensitive to filter time delays. To prevent the introduction of additional delays to the signal path, an alternative preconditioning method is to use spectral factorization techniques to find the causal preconditioning filters. The decomposition process is straightforward for single-input-and-single-output (SISO) systems. For Multi-input Multi-output (MIMO) system, the factorization is significantly more di ffi cult. Youla was among the first to give the decomposition method for a MIMO system by taking advantages of the properties of Smith canonic form of a polynomial matrix [7]. However, Youla’s method is numerically cumbersome and di ffi cult to implement when the order of the matrix is high. Recently, several new methods have also been developed to reduce the calculation time of the spectral factorization [8, 9]. Among these methods, the group of methods based on solving a Riccati equation summarized by Kailath [10] has unique advantages over other methods due to the convenience in formulation and the development of algorithms in recent years to solve a discrete-time Riccati equation (DARE) [11]. The method can be extended to MIMO case from SISO case conveniently due to the state space representation of the algorithm thanks to the earlier work of Kalman [12]. In the method proposed in the current work, causal preconditioning filters for the reference signals and for the secondary path plant response in a MIMO real-time active noise control system are designed with a spectral factorization technique based on solving a discrete-time Riccati equation. Then the causal preconditioning filters are converted back to the time domain so that they can be included into the filtered reference LMS (FxLMS) algorithm to accelerate the algorithm’s convergence for adaptation. This paper includes three sections. The proposed causal preconditioning filter design approach based on polynomial matrix spectral factorization is discussed in Sec. 2. Numerical simulations using the proposed preconditioning techniques are given in Sec. 3 to show the performance of the proposed preconditioning filters in accelerating the convergence of FxLMS algorithms and improving the noise reduction performance in a real-time MIMO active noise control system. The conclusion of the current study is given in Sec. 4.

2. CAUSAL MULTICHANNEL PRECONDITIONING FILTER DESIGN

2.1. Filtered-reference LMS algorithm A typical filtered-reference LMS algorithm for an ANC application is illustrated in Figure 1. The ANC filter W ( z ) is a M × K dimension system which takes K real-valued reference signals x ( n ) to generate M control signals y ( n ) as the input to the secondary path G ( z ). The K discrete-time reference signals can be described by the vector

x ( n ) = [ x 1 ( n ) ... x K ( n )] T , (1)

Figure 1: Filtered-reference LMS algorithm.

Figure 2: Active noise control system using preconditioned Filtered-reference LMS algorithm.

where n is the time index. Similarly, the M control signals y ( n ) can be represented by the vector

L disturbance M sina Mxk LxM contrat fires SOMO! plane (M2) F re error e(n) G@) LxM plant response model

y ( n ) = [ y 1 ( n ) ... y M ( n )] T . (2)

The plant response matrix G ( z ) is a L × M dimension system which takes M input control signals to generate L output ’anti-noise’ signals. G ( z ) includes the response of the data converters, the analogue anti-aliasing, reconstruction filters, processing delays, and the acoustic propagation. The L disturbance signals to be suppressed by the ANC system can be expressed by

d = [ d 1 ( n ) ... d L ( n )] T . (3)

The L error signals are the same as the disturbance signals when the ANC system is not activated. The error signal vector can be represented as

e = [ e 1 ( n ) ... e L ( n )] T . (4)

The adaptive algorithm of the filtered-reference LMS system is then expressed as

L X

w mk ( n + 1) = w mk ( n ) − µ

l = 1 [ e l ( n ) r lmk ( n )] , (5)

where w mk is the impulse response of the ANC filter W ( z ). r lmk is calculated from discrete convolution of the modeled plant impulse responses and the reference signals, which is expressed as

J − 1 X

r lmk ( n ) =

j = 0 ˆ g lm j ( n ) x k ( n − j ) , (6)

where ˆ g lm j is the impulse response of the modeled secondary path, ˆ G ( z ) with the m -th control signal, the l -th output channel at the j -th time sample.

a(n) zo) FQ)

2.2. Causal preconditioning filters for MIMO FxLMS algorithms If causal preconditioning filters can be obtained from the information of the reference signals and secondary path without introducing significant additional delays, a preconditioning approach can be proposed to accelerate the convergence of FxLMS algorithms for real-time ANC applications. A FxLMS system with causal preconditioning filters is presented in Figure 2. In the proposed preconditioning algorithm, matrix F ( z ) is a minimum phase system derived from S xx ( z ), which is the cross spectral density matrix of the reference signals. The operation to obtain minimum phase F ( z ) from S xx ( z ) is usually called spectral factorization, and it can be expressed as

S xx ( z ) = F ( z ) F T (1 / z ) . (7)

According to the definition of the cross spectral density matrix, S xx ( z ) can also be expressed as

S xx ( z ) = E [ x ( z ) x T (1 / z )] , (8)

where x ( z ) is the z transforms of references signals x ( n ) and E is the expectation operator. Based on the spectral factorization operation, the reference signals x ( n ) can be interpreted as signals to be generated by passing uncorrelated white noise signals with unit signal power v ( n ) through a minimum phase shaping filter F ( z ), i.e., x ( z ) = F ( z ) v ( z ) . (9)

It is noted that the cross spectral density matrix of the above mentioned white noise signals v ( n ) is S vv ( z ) = I K , where I K is an identity matrix of size K × K . Since F ( z ) is causal and minimum phase, its inverse is also causal and minimum phase. Thus, the reference signals can be whitened by the F − 1 ( z ) filter as

v ( z ) = F − 1 ( z ) x ( z ) , (10)

where F − 1 ( z ) is the inverse system of F ( z ). The prewhittened reference signals v ( z ) can thus be used as the new reference signals for the FxLMS algorithm. The convergence rate can be improved with these new reference signals since they are uncorrelated. This improvement can be explained by the reduction of the eigenvalue spread of the auto-correlation matrix of the reference signals [4]. The M × M matrix G min ( z ) and the L × M matrix G all ( z ) are minimum phase and all-pass components of the plant responses G ( z ) respectively. Their relation with the plant responses G ( z ) can be expressed as G ( z ) = G all ( z ) G min ( z ) . (11)

Here the all-pass component satisfies the property that G T all ( z ) G all ( z ) = I M . In the proposed FxLMS algorithm (or proposed LMS in short), the filter coe ffi cients can be updated according to the expression in Equation 5. The only di ff erence is, in contrast to the original reference signals filtered by G ( z ), that a new sets of filtered reference signals should be used in the proposed LMS algorithm, which is calculated from equation

J − 1 X

r lmk ( n ) =

j = 0 g all , lm j v k ( n − j ) , (12)

where g all , lm j is the lmj -th entry of the all-pass component in the time domain, and v k is the k -th signal of the whitened reference signal. The decomposition process in Equation 7 can be accomplished through solving a Riccati equation of the form Σ = C Σ C H − ( C Σ h H − ¯ N )( h Σ h H − R 0 ) − 1 ( C Σ h H − ¯ N ) H (13)

, where each term in Equation 13 can be calculated using the cross spectral density matrix S xx ( z ) according to [13,14]. Similarly, the minimum / all-pass decomposition of the plant response matrix G , can be performed by applying the spectral factorization to a para-conjugate matrix formed by

G H ( z ) G ( z ) = G H min ( z ) G min ( z ) . (14)

It is obvious that the spectral factorization of G H ( z ) G ( z ) results in the minimum phase component of the plant response matrix, G min ( z ). Then the all-pass component can be derived as

G all = G ( z ) G − 1 min ( z ) . (15)

3. NUMERICAL SIMULATION

In this section, two simulations are demonstrated. In the first simulation, simple white noise and pink noise source models and free field propagation models of the secondary path are used. The all-pass / minimum phase decomposition of the plant response of the secondary path can be obtained through observations. The purpose of the first simulation is to verify the accuracy of the proposed spectral factorization approach with analytical solution of the decomposed all-pass filter and minimum phase filter. In the second simulation, to prove that the decomposition method can be applied to more general types of plant response functions, a realistic plant response is used which was constructed through experimental measurement.

3.1. Simulation study on systems with a simple plant response A numerical simulation is conducted to validate the key steps in the proposed preconditioning filter design method for multichannel LMS algorithms and to demonstrate the associated convergence rate improvement brought by the preconditioning process. A white noise sequence and an uncorrelated pink noise sequence are used as the two independent sources v ( n ) = [ v 1 ( n ) , v 2 ( n )] T . Two sources are then mixed through a mixing filter M , which is defined as

 0 . 75 0 . 3

 . (16)

M =

0 . 3 1

The resulted two mixed signals x 1 ( n ) and x 2 ( n ) are then used as the reference sources for the preconditioning LMS system. The sampling rate, f s , of the simulation is 1 kHz. Two error microphones are symmetrically arranged. The same transfer functions for the secondary path are used as those used in the work of Bai and Elliot [6], which can be described as

 , (17)

 z − N 1 l 1 l 2 z − N 2

G ( z ) =

l 1 l 2 z − N 2 z − N 1

where N 1 is the nearest integer value of l 1 f s / c 0 and N 2 is the nearest integer value of l 2 f s / c 0 . l 1 = 2 m and l 2 = 3 m are the distances from the first secondary source to the first and the second error microphone respectively. c 0 = 340 m / s is the speed of sound. Since the secondary sources and error microphones are symmetrically arranged, the G ( z ) matrix is symmetric. Additionally, the disturbance signals d 1 ( n ) and d 2 ( n ) are calculated using the same transfer functions in Equation 17, but with di ff erent distances: l 1 = 3 m and l 2 = 4 m . The polynomial form of S xx matrix is then identified using the factorization techniques. The identified S xx ( z ) matrix is subsequently decomposed using the spectral factorization algorithm by solving the corresponding discrete Riccati equation. The decomposition matrix F ( z ) can then be used

Figure 3: Sound pressure levels at microphone 1. From top to bottom: (a) SPL of the disturbance signal; (b) original LMS algorithm; (c) LMS with G − 1 min and G all ; (d) LMS with F − 1 ; (e) LMS with F − 1 , G − 1 min and G all .

to find its inverse F − 1 ( z ), which, due to its minimum phase nature, is a causal prewhitening filter for the reference signals. The reference signals can then be fed into the prewhitening filter F − 1 ( z ), which produces a new pair of white noise signals ¯ v 1 and ¯ v 2 . The power level of the filtered reference signals become almost a constant for both S ¯ v ¯ v (1 , 1) and S ¯ v ¯ v (2 , 2). At the same time the S ¯ v ¯ v (1 , 2) term is more than 10 dB smaller than the diagonal terms, which suggests the whitening process is successful. It also suggests that although the whitening filter can be calculated with spectral factorization algorithm, the original pre-mixed signals v 1 ( n ) and v 2 ( n ) cannot be recovered through the approach. The next step is to remove the coupling in the secondary path by applying a minimum phase / all- pass decomposition. Using the techniques introduced in Section 2, the minimum phase component G min can be obtained first through a spectral factorization process. Subsequently, the all pass component G all ( z ) can be calculated with Equation 15. The impulse response of G all ( z ), which is modeled as a FIR filter is the same as the analytical solution obtained via observation [3]:

 z − N 1 0

 . (18)

G all =

0 z − N 1

In order to demonstrate the improvements in ANC performance as well as convergence rate brought by the proposed preconditioning filters, the sound pressure level at the first error microphone is plotted in Figure 3. The curves a, b, c, d, and e denote respectively the SPLs (sound pressure levels) of the disturbance signal, the error signal of the traditional LMS algorithm, the error signal of the LMS algorithm preconditioned with G − 1 min and G all , the error signal of the LMS algorithm preconditioned with F − 1 , and the error signal of the LMS algorithm preconditioned with G − 1 min , G all and F − 1 . These levels are averaged over 30 samples for better readability. The results at the second microphone is similar and omitted. The convergence coe ffi cient µ in Equation 5 is set to one third of the lowest value that results in instability. It can be observed from Figure 3 that the original LMS algorithm has the worst performance, which has an around 10dB reduction of noise after 4000 samples. After introducing the prewhitening filter F − 1 , the reduction increases to 25 dB after 4000 samples. If G all is introduced without the prewhitening filter F − 1 , the reduction can reach approximately 14 dB, which is slightly better than the original FxLMS algorithm. If both F − 1 and G all are used, the noise reduction reaches 30dB after only 2000 samples. In the case, the introduction of the prewhitening filter F − 1 is more e ff ective than the introduction of decoupling filters G all and G − 1 min .

‘SPL at the error micphone (dB) Ait Cc atl 2000 4000 6000 ‘2000 +0000 ‘sample number

Figure 4: Sound pressure levels at microphone 1. From top to bottom: (a) SPL of the disturbance signal; (b) original LMS algorithm (indistinguishable to the original disturbance signal); (c) LMS with G all and G − 1 min ; (d) LMS with F − 1 , G − 1 min and G all .

3.2. Simulation study on systems with a realistic G ( z ) In the previous simulation case study, a simple secondary path frequency response function based on free field propagation method is used. The simulations suggests that the minimum phase / all-pass decomposition using spectral factorization for the simple frequency response function is successful. The improvement by introducing a secondary path decomposition can be observed, but is not particular obvious, which could possibly due to the weak coupling in the secondary path. In the following simulation, a realistic plant response obtained from experimental measurement in an anechoic chamber is used. In this case, G ( z ) has stronger inter-channel coupling than the previously used propagation matrix. The reference signals and the disturbance singles used in this simulation are the same as those used in the first simulation. The filter G ( z ) is then decomposed into the minimum phase component G min ( z ) and all-pass component G all ( z ) by the spectral factorization method described in Section 2. All the decoupling filters are causal and no additional delays are needed to be added to the impulse responses. The final comparison of the received noise signals obtained by di ff erent versions of FxLMS algorithms at a representative error microphone is plotted in Figure 4, where the convergence coe ffi cient µ is set to half of the lowest value that results in instability. The curves a, b, c, and d denote respectively the SPLs of the disturbance signal, the error signal of the traditional LMS algorithm, the error signal of the LMS algorithm preconditioned with G − 1 min and G all , and the error signal of the LMS algorithm preconditioned with G − 1 min , G all and F − 1 . In this simulation, the noise cannot be e ff ectively reduced by the original FxLMS algorithm due to the strong coupling in both the reference signals and in the plant response. By introducing the decoupling filters for the plant G − 1 min ( z ) and G all ( z ), the reduction could reach 10 dB after 4000 samples, which suggests a much faster convergence than the traditional FxLMS algorithm. No significant improvement was found in this simulation case study by solely incorporating prewhitening filter F − 1 ( z ) and thus this result is not plotted in the figure. By introducing preconditioning filters for both the reference signals and the plant response (curve d of Figure 4), the convergence speed is improved and an around 5dB additional noise reduction performance was achieved when compared with curve c. It should also be noted that, although the decoupling filters are e ff ective in increasing the speed of the convergence, the required filter length increases with the complexity of the plant response G ( z ), thus a trade-o ff between the convergence speed and the length of the preconditioning filters must be considered in actual engineering applications. In the two simulation cases described in this work, the time spent to calculate the spectral factorization of S xx ( z ) and G H ( z ) G ( z ) on a personal computer using Ryzen 3600x CPU is around 4 seconds.

Time (s) “4000 ‘6000 ‘2000 10000 2000 (gp) auoydonu 0149 ou9 12 1S ‘sample number

4. CONCLUSION

A prewhitening / preconditioning technique is proposed based on spectral factorization to improve the convergence speed and the performance of the FxLMS algorithm for real-time multi-channel active noise control applications. The causal nature of these derived prewhitening filters allows them to be directly implemented in a real-time active noise control application. By applying the prewhitening filter and the decoupling filter to the system, the convergence speed of the LMS algorithm and the associated active noise control performance are shown to be significantly improved.

5. ACKNOWLEDGMENTS

The authors thank Beijing Anesonic Technology Co. Ltd for providing financial support for the present work.

REFERENCES

[1] John C Burgess. Active adaptive sound control in a duct: A computer simulation. The Journal of the Acoustical Society of America , 70(3):715–726, 1981. [2] Stephen Elliott, IANM Stothers, and Philip Nelson. A multiple error lms algorithm and its application to the active control of sound and vibration. IEEE Transactions on Acoustics, Speech, and Signal Processing , 35(10):1423–1434, 1987. [3] Stephen J Elliott. Optimal controllers and adaptive controllers for multichannel feedforward control of stochastic disturbances. IEEE Transactions on signal Processing , 48(4):1053–1060, 2000. [4] Stephen Elliott. Signal processing for active control . Elsevier, 2000. [5] Scott C Douglas, Andrzej Cichocki, and Shun-ichi Amari. Self-whitening algorithms for adaptive equalization and deconvolution. IEEE Transactions on Signal processing , 47(4):1161– 1165, 1999. [6] MR Bai and SJ Elliott. Preconditioning multichannel adaptive filtering algorithms using evd- and svd-based signal prewhitening and system decoupling. Journal of sound and vibration , 270(4-5):639–655, 2004. [7] D Youla. On the factorization of rational matrices. IRE Transactions on Information Theory , 7(3):172–189, 1961. [8] Gigla Janashia, Edem Lagvilava, and Lasha Ephremidze. A new method of matrix spectral factorization. IEEE Transactions on information theory , 57(4):2318–2326, 2011. [9] Zeliang Wang, John G McWhirter, and Stephan Weiss. Multichannel spectral factorization algorithm using polynomial matrix eigenvalue decomposition. In 2015 49th Asilomar conference on signals, systems and computers , pages 1714–1718. IEEE, 2015. [10] Thomas Kailath, Ali H Sayed, and Babak Hassibi. Linear estimation . Number BOOK. Prentice Hall, 2000. [11] Hisham Abou-Kandil, Gerhard Freiling, Vlad Ionescu, and Gerhard Jank. Matrix Riccati equations in control and systems theory . Birkhäuser, 2012. [12] Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. 1960. [13] Ali H Sayed and Thomas Kailath. A survey of spectral factorization methods. Numerical linear algebra with applications , 8(6-7):467–496, 2001. [14] Vladimir Kucera. Factorization of rational spectral matrices: a survey of methods. In International Conference on Control 1991. Control’91 , pages 1074–1078. IET, 1991.