A Stable IIR Filter Design Approach for High-Order Active Noise Control Applications

: In commercial non-adaptive active noise control (ANC) applications, an IIR ﬁlter structure is often used to reduce real-time computations. On the contrary, an FIR ﬁlter structure is usually preferred in the ﬁlter design phase because the FIR ﬁlter design formulation can be convex and is simple to solve. To combine the beneﬁts of both FIR and IIR ﬁlter structures, one common approach in ANC applications is to use an IIR ﬁlter structure to ﬁt a pre-designed FIR ﬁlter. However, to ensure stability, most of the common IIR ﬁlter ﬁtting approaches involve the computation and relocation of poles which can be difﬁcult for high-order cases. In this current work, a stable IIR ﬁlter design approach that does not need the computation and relocation of poles is improved to be applicable in ANC applications. The results demonstrate that the proposed method can achieve better ﬁtting accuracy and steady-state noise control performance in high-order non-adaptive applications when the pre-designed noise control FIR ﬁlter is ﬁtted. Besides ﬁtting the noise control ﬁlter, the proposed method can also be used to ﬁt the secondary path and acoustic feedback path to reduce the required real-time computations if adaptive controllers are applied.


Introduction
Active noise control (ANC) techniques have been implemented in a wide range of applications in the past decades, such as headrests [1], automobiles [2], motorcycle helmets [3], air-conditioners [4,5], ventilation windows [6], etc.In most of these practical applications, finite impulse response (FIR) filters were used because of their simplicity and stability in real-time implementation [7].Besides the real-time implementation phase, the use of a FIR filter structure can be more beneficial in the design phase because the filter response will be a linear function of the FIR filter coefficients to be designed and it simplifies the filter design problem formulation.Moreover, various practical constraints shall be considered when designing the ANC filter coefficients in many commercial ANC applications, e.g., headphones [8], vehicles [2], and air-conditioners [4,5].The constrained filter design problem can be formulated to be a convex problem if an FIR filter structure is used in the control filter and the global optimal filter coefficients can be obtained [2,4].
However, compared with the infinite impulse response (IIR) filter, the FIR filter usually requires a longer filter length, thus, a higher real-time computational power.In many commercial ANC applications, a high sampling rate may be chosen for multiple reasons [9], e.g., to reduce the electronic delay to improve the noise control performance [10], to achieve a wider frequency band [11], to be compatible with the sampling rate in audio applications, etc.The use of a high sampling rate in ANC systems significantly increases the required order of control filters if an FIR filter structure is used [9,11].Thus, an IIR filter structure is preferred in high-sampling cases to reduce the real-time computations at the sacrifice of simplicity in the filter design phase.
One common approach in ANC applications to combine the benefits of both FIR and IIR filter structures is to use an IIR filter structure to fit pre-designed FIR filter coefficients (or the desired frequency responses).Thus, the desired time-domain filter coefficients or frequency-domain responses can be designed first using the FIR filter structure to simplify computation in the filter design phase.Then, an IIR filter is fitted from the desired timedomain coefficients or the frequency-domain responses.The most challenging part of using this approach is that the designed IIR filter should be stable.Many traditional IIR filter fitting methods cannot ensure stability directly, such as the prony's method [12], stmcb method [13,14], the least-square approximation in the frequency domain using equation error method [15], or the damped Gauss-Newton method [16].To ensure stability, most of the methods compute the poles and map the unstable poles to the inside of the unit circle in the z-transform domain at each iteration.However, for high-order cases, numerically computing the poles (i.e., computing the roots of a high-order polynomial) can be difficult [17,18], which makes many traditional stable IIR filter approximation approaches not suitable to be used in high sampling rate ANC applications.Thus, a stable IIR approximation approach for high sampling rate ANC applications should be developed.
In this paper, to fit a stable IIR filter for non-adaptive ANC applications, the IIR filter approximation method proposed by Brandenstein and Unbehauen [19,20] is adopted.Improvements are proposed based on their method to make it applicable in ANC applications.A non-adaptive ANC system installed on a commercial earphone prototype is used to investigate the performance of the proposed method.The result shows that the proposed improved IIR filter approximation method can achieve better fitting accuracy and noise control performance compared with the traditional IIR fitting approach.Besides fitting the noise control filter, the proposed method can also be applied to the fitting of secondary paths and acoustic feedback paths to further reduce the real-time computations if adaptive controllers are applied.

Methods
In this section, the commonly used traditional stable IIR filter fitting approach is first reviewed.Then, the method proposed by Brandenstein and Unbehauen is reviewed and the proposed improvements based on their method are explained.Only FIR and IIR filters with real-valued time-domain coefficients are considered in this article.
Let the transfer function of the FIR filter that needs to be approximated be The transfer function of a causal and stable IIR filter is expressed as where, the numerator polynomial P(z) is and the denominator polynomial Q(z) is The objective of fitting an FIR filter using an IIR filter structure is to minimize the approximation error with respect to p µ (µ = 0, 1, . . ., M) and q ν (ν = 1, 2, . . ., N), where W(e jΩ ) denotes the weighting function and W(e −jΩ ) = W(e jΩ ) if the fitted FIR filter has only real coefficients.

Review of Traditional Stable IIR Filter Fitting Method
One commonly used traditional stable IIR filter fitting method is reviewed here.For simplicity, this approach is referred to as the "traditional method" in this article.
For this traditional method, instead of the objective function value E used in Equation ( 5), a discrete form of fitting error is formulated for numerical computation: where Ω i is the i-th sampled frequency, and N f is the total number of considered discretized frequency points.If the total number of frequencies is chosen to be sufficiently large, then minimizing E s in Equation ( 6) can lead to approximately the same optimal filter coefficients compared with minimizing E in Equation (5).To find the optimal filter coefficients that minimize E s in Equation ( 6), some iterative search algorithms can be implemented, e.g., the damped Gauss-Newton method [16].The initial point for the iterative search algorithms can be chosen using the optimal coefficients solved by minimizing the equation error formulation [15]: W(e jΩ i ) F(e jΩ i )Q(e jΩ i ) − P(e jΩ i ) Better fitting accuracy may be achieved by changing the limit of the integral in Equation ( 5), or the lower and upper bounds in summation in Equations ( 6) and (7) to the desired noise control range instead of the whole frequency range.However, this may lead to inaccurate fitting outside the desired frequency band and cause noise amplifications.Thus, in this presented work, the fitting procedure is still targeting the whole frequency range but a weighting function W is used to improve the fitting accuracy in the desired frequency range.
To ensure stability, at each iteration, the roots of the polynomial Q(z), i.e., the poles, should be computed and checked.If some roots z i are outside the unit circle in the z-transform domain (i.e., |z i | > 1), they should be stabilized by replacing z i with 1/z * i where z * i is the complex conjugate of z i .This traditional method is used in many IIR filter fitting toolboxes such as the invfreqz() function in Matlab R2022b.
The challenge of using this traditional method for high-order IIR filter cases is that computing the roots of a high-order polynomial can be difficult and there is no guarantee that the traditional method can converge when the pole stabilizing process is applied.

Review of a Weighted Least-Square Approximation Method
A least-square method for approximating an FIR filter using an IIR filter structure that guarantees the fitted IIR filter is stable without mapping unstable poles inside the unit circle was proposed by Brandenstein and Unbehauen [19] and then extended to its weighted version [20].Since no computation and mapping of poles are required, it is applicable even in cases where the IIR filter length is large [19,20].Another advantage of their proposed method is that the optimal denominator coefficients can be fitted first.The optimal numerator coefficients can then be computed using the optimal denominator coefficients.For simplicity, their original proposed method is referred to as the "BU's method" in this article.Their method is briefly summarized as follows (See Brandenstein and Unbehauen's work for more details [19,20]).
A user-defined weighting function W(e jΩ ) is usually real and positive (i.e., having zero phases).However, the phase of the weighting function W(e jΩ ) in Equation ( 6) can be arbitrary as long as the magnitude satisfies the required weighting shape.Brandenstein and Unbehauen proved that if a maximum phase weighting function is used, it can make their proposed method likely to converge to the optimal solution.[19,20] Thus, assume that there is a maximum phase function G(z) with order K satisfying |G(e jΩ )| 2 = W(e jΩ ).Then, a function X(z) can be defined as where G * (z) denotes the polynomial which is obtained from G(z) by replacing each coefficient with its conjugate complex counterpart, and F * (z) is defined in a similar way.By using X(z), starting from an initial guess Q 0 (z) = 1, the optimal coefficients in the denominator of the IIR filter Q(z) can be solved iteratively [19,20].More specifically, suppose at the k-th iteration, X (k) (z) is defined as where In practice, Equation ( 9) means the X (k) can be obtained by filtering the coefficients of X(z) using an IIR filter with 1 in the numerator and Q (k−1) (z) in the denominator.Then, Q (k) (z) can be obtained by solving the set of equations in a least-square sense, i.e., q where and x (k) (n) is the sequence obtained by the inverse z-transform of X (k) (z).The solved q (k) consists of the coefficients of Q (k) (z), The approximation error norm at k-th iteration is defined as where u (k) (n) is the sequence obtained by the inverse z-transform of U (k) (z) and After the optimal denominator coefficients, Q(z) is obtained, and the optimal numerator coefficients P(z) can be obtained by solving the M + 1 linear equations for λ = 0, 1, . . ., M, where |Q(e jΩ )| 2 e jnΩ dΩ (17) 2.3.Proposed Improvements in Weighted Least-Square Approximation Method for ANC Applications In this section, some improvements are made to the BU's method such that it is more suitable for practical ANC applications.For simplicity, these improvements on BU's method are referred to as the "proposed method" in this article.
When using the iteration process described in the previous section, the approximation error norm E (k) does not always converge to the minimum value, especially for some maximum phase FIR filters with high orders (larger than 50) [19].Thus, instead of using the Q (k) (z) solved at the last iteration, Brandenstein and Unbehauen suggested trying a sufficiently large number of iterations and then choosing the one with the lowest approximation error norm [19].In ANC applications, the order of the FIR filter is usually larger than 50.This phenomenon is more likely to happen.Since this phenomenon occurs when the matrix (k) is near singular at the current iteration, then a gradient descent is applied to continue solving Q(z).More specifically, the goal of solving Equation ( 10) is to minimize Thus, ˆ g, an estimation of the gradient of the function J( q (k) ) can be obtained by using The approximated q (k) can then be obtained by where α > 0 is the step length such that Q (k) (z) is stable and can have a smaller approximation error norm compared with Q (k) (z − 1).Repeat this process until some prescribed convergence criteria are satisfied, e.g., step length, norm of gradient, etc.Note that this proposed gradient method is implemented once in each step but if it is implemented sufficiently many times in each step and well-conditioned, it will likely approach the solution obtained using BU's original method.The ability to use prescribed convergence criteria in the proposed gradient method can be helpful to prevent the issues caused by the large condition number.
The selection of maximum phase weighting function G(z) also needs to be considered.In ANC applications, the required weighting function is usually in the form of |W(e −jΩ )|.For simplicity, a linear-phase FIR filter with real coefficients can be easily designed to have the magnitude |W(e −jΩ )|, such that W(e −jΩ ) = W(e jΩ ) is satisfied.Because of the linear phase, the coefficients of this designed filter will be symmetric such that a spectral factorization can be applied to the coefficient sequence of this FIR filter to obtain a minimum phase decomposition.After flipping the coefficients of decomposed minimum phase filter, the maximum phase counterpart G(z) can be obtained.
Since the signal used in the ANC controller is usually in a discrete format, the computation of c(n) and b(n) in Equations ( 17) and ( 18) is done by FFT and IFFT after a sufficiently large zero padding.
In summary, the whole proposed process of fitting an FIR filter using an IIR filter structure is presented in Algorithms 1-3.Algorithm 2: Finding the optimal IIR filter denominator coefficients the IIR filter numerator order M and denominator order N, the weighting filter coefficients G = g 0 g 1 . . .g K T in Algorithm 1.
← Filtering sequence X using IIR filter with numerator 1 and denominator coefficients Q (k−1) , keep the first L + N − M points.Construct A (k) , b (k) using Equations ( 11) and (12).
E temp ← Filtering X using IIR filter with numerator coefficients flip( Q temp ) and denominator coefficients Q temp , STOP ← 0, and break the loop.else Reduce α E (k−1) ← Filtering X using IIR filter with numerator coefficients flip( Q (k) ) and denominator coefficients Q (k) , compute the 2-norm of first L + N − M points.if Stoping criterion satisfied then Break the loop Algorithm 3: Finding the optimal IIR filter numerator coefficients the order of IIR numerator M, the weighting function G = g 0 g 1 . . .g K T in Algorithm 1, IIR denominator Q = q 0 q 1 . . .q N T in Algorithm 2.
Output : fitted IIR numerator P = p 0 p 1 . . .p M T S W ← square of the magnitude of FFT of sufficiently zero-padded sequence G S Q ← FFT of sufficiently zero-padded sequence Q S F ← FFT of sufficiently zero-padded sequence F c ← IFFT of S W pointwise divided by the square of magnitude of S Q , truncate to M + 1 points d ← IFFT of S W pointwise multiplied by S F and then pointwise divided by the complex conjugate of S Q , truncate to M + 1 points C p ← construct a square real symmetric Toeplitz matrix using c T as the first row P ← solve C p P = d, Toeplitz properties can be applied

Experimental Setup Description
To investigate the fitting performance and its effect on the noise control performance of the proposed method, an earphone is used to collect the system response data.For this earphone, there is one microphone used as the reference sensor and one speaker used as the noise control speaker.The error microphone used in this experiment that measured the noise control performance is in an artificial ear shown in Figure 1.The artificial ear used is Hangzhou Aihua AWA6162 (IEC711).The broadband noise control performance is considered in this section.Because of the limited dimension of the earphone, the delay in the signal path introduced by the electronic devices and anti-aliasing/reconstruction filters should be sufficiently small to achieve broadband noise control performance.Thus, the sampling rate of the data acquisition system and the noise control system was set to 48 kHz.Although the sampling rate is high (48 kHz), the desired noise control band is set from 50 Hz to 3 kHz for practical considera-tions.When measuring and computing the system responses, two million sampling points (about 42 s) are acquired and a hamming window of 48,000 points (around 1 s) is used for averaging with 50% overlap.
When the acoustic feedback path effect is canceled by an estimate of the acoustic feedback path, the block diagram of the ANC system installed in the tested earphone is shown in Figure 2. x denotes the noise signals from the primary noise source at the reference microphone location.White noise played by the speaker is used as the noise signal at the noise source in this experimental setup.y denotes the output of the noise control filter and W. G e denotes the secondary path which is the acoustical response of the secondary source at the error sensor position.The frequency responses of the measured secondary path are shown in Figure 3. d e denotes the disturbance signals and e denotes the sum of primary and secondary sound fields at error sensor locations, whose total power is attenuated by the ANC system.

Comparison of Fitting Accuracy and Noise Control Performance
Because the desired noise control band is from 50 Hz to 3 kHz and the sampling rate is 48 kHz, which is much higher than the desired noise control band, some constraints, e.g., the filter responses amplitude constraints or the disturbance enhancement constraints should be applied at the higher frequency band to prevent the amplification of the highfrequency noises outside the desired noise control frequency band.Thus, the time-domain coefficients of the noise control filter W were designed firstly using an FIR filter structure with constrained optimal filter design methods [4,5].Since the FIR filter is not the filter implemented in real-time, a sufficiently large order should be chosen.Thus, the FIR filter length is increased until the noise control performance cannot be improved any further.The nominal noise control performance using the designed FIR filter is shown in Figures 4 and 5. From Figure 4, the designed FIR filter can achieve satisfactory noise control performance in the desired noise control frequency band and does not amplify noises outside the desired noise control frequency band.Figure 5 demonstrates that a large filter length (i.e., 7200) is needed because, compared with using 2400 and 4800 coefficients, using 7200 coefficients (equivalent to 0.15 s in the time-domain filter response length) can still improve the noise control performance, especially around the lower frequency band (3-4 dB better than 2400 coefficients case and around 1 dB better than 4800 coefficient cases).Since the improvement from using 4800 coefficients to using 7200 coefficients are not significant, 7200 coefficients are considered sufficient for this ANC application.The impulse response and frequency response of the designed FIR filter using 7200 coefficients are shown in Figure 6.The impulse response and the frequency response of designed optimal control filter using an FIR filter structure using 7200 time-domain coefficients.
Although using 7200 coefficients can achieve satisfactory noise control performance at the desired frequency band without amplifying noises outside the frequency band, the real-time computational requirement is significant if a 7200-point FIR filter is implemented at a 48 kHz sampling rate.Thus, an IIR filter structure should be used to fit the responses of this 7200-point FIR filter.When using the IIR filter structure to fit the pre-designed FIR filter in Figure 6, the order of the numerator polynomial is set to be equal to the order of the denominator polynomial.Because the desired noise control (50 Hz to 3 kHz) is relatively narrow compared with the full band, 100 times higher weighting is applied to frequencies below 3 kHz to ensure an accurate fitting in the desired noise control band.When using the proposed method, a condition number higher than 10 12 is considered large and the estimated gradient is used.Since the condition number is used as a threshold and the computation of coefficients does not need the exact value of the condition number, efficient estimates of condition numbers will be sufficient, e.g., the rcond() function in Matlab [21].
When using the traditional method, if the order of the IIR filter is chosen to be higher than 64, the algorithm cannot converge to a stable IIR filter.However, the proposed method can successfully fit the frequency responses using IIR filters with higher orders (e.g., 128 or 256).The comparison between pre-designed FIR filter frequency responses and the fitted IIR filter frequency responses using different methods and filter orders are shown in Figure 7a.From Figure 7a, all methods and filter orders can achieve satisfactory overall fitting accuracy for the whole frequency band, especially at higher frequencies.A zoomed-in frequency band from 100 Hz to 800 Hz is shown in Figure 7b to have a clearer comparison of the fitting accuracy.It demonstrates that during this relatively lower frequency range, the IIR filter fitted using the proposed method has a much better fitting accuracy compared with that of using the traditional method because the proposed method enables a higher-order stable IIR filter to be fitted.The comparison of nominal noise control performance using the traditional method and the proposed method is shown in Figure 8 for the full frequency band and the desired noise control band.The zoomed-in frequency band from 100 Hz to 800 Hz is also shown in Figure 9.This demonstrates that the proposed method can achieve a better noise control performance compared with the traditional method, especially for certain frequency bands, e.g., 3-5 dB better from 400 Hz to 600 Hz.

Discussion and Conclusions
In this article, a stable IIR fitting approach in high-order ANC applications is proposed by improving Brandenstein and Unbehauen's method [19,20] such that the stability of the designed IIR filter can be guaranteed without implementing root-finding.A gradient-based searching method is proposed in this article in addition to the original solving approach in BU's method such that convergence can be guaranteed, which is more applicable in practical applications.Practical approaches to obtain the required maximum phase weighting function and compute the optimal IIR filter numerator coefficients are also proposed in this article.An ANC system installed on an earphone prototype is used to investigate the performance of the proposed method.The result shows that the proposed improved IIR filter approximation method can achieve better fitting accuracy and noise control performance compared with the traditional IIR fitting approach.
Besides the fitting of the noise control filter, the proposed method can also be used for fitting secondary paths and acoustic feedback paths to further reduce real-time computations.Moreover, the proposed method can be used for fitting the secondary paths and acoustic feedback paths in both non-adaptive and adaptive applications.
Since no root relocation is required during the fitting process, the proposed method is more efficient and reliable and it may be extended in the real-time continuous redesign process where a noise control filter is designed by solving the optimal filter coefficients continuously [22].

Algorithm 1 :
Finding an appropriate weighting function Input : desired weighting curve |W(e jΩ )| in the frequency domain; the order of weighting filter K Output : the weighting coefficients G = g 0 g 1 . . .g K T • Designed a linear phase filter coefficients W z = w 0 w 1 . . .w 2K T that has frequency magnitude |W(e jΩ )| (E.g., using Parks-McClellan or least-square optimal FIR filter design).• Apply spectral factorization to W z to obtain a minimum phase decomposition G min = g0 g1 . . .gK T .• The maximum phase weighting filter coefficients are G = flip( G min ) = gK gK−1 . . .g0 T .

Figure 1 .
Figure 1.A picture of the experimental setup showing an earphone mounted on an artificial ear for experimental data collection.

Figure 2 .
Figure 2. The block diagram of an ANC system installed in the earphone when the acoustic feedback path effect is canceled by an estimate of the acoustic feedback path.

Figure 3 .
Figure 3.The frequency responses of measured secondary path G e .The reference value for dB scale is 1.

Figure 4 .Figure 5 .
Figure 4.The nominal noise control performance of designed noise control filter using an FIR filter structure with different filter lengths.

Figure 6 .
Figure6.The impulse response and the frequency response of designed optimal control filter using an FIR filter structure using 7200 time-domain coefficients.
(a) The full frequency band from 0 Hz to 24 kHz.(b) Zoomed-in frequency band from 100 Hz to 800 Hz.

Figure 7 .
Figure 7.The comparison of pre-designed FIR filter frequency responses and the fitted IIR filter frequency responses using different methods and filter orders.
(a) The full frequency band from 0 Hz to 24 kHz.(b) Desired noise control frequency band from 0 Hz to 3 kHz.

Figure 8 .
Figure 8.The comparison of nominal noise control performance using the traditional method and the proposed method.

Figure 9 .
Figure9.Zoomed-in frequency band from 100 Hz to 800 Hz for the comparison of nominal noise control performance using the traditional method and the proposed method.