The Predict and Invert Feedback Active Noise and Vibration Control Algorithm

This work proposes an algorithm for feedback ANC that does not require a prior secondary path model and usually remains stable after fast secondary path changes, as other algorithms proposed for feedforward ANC. This is achieved using a recursive least squares algorithm to model the secondary path and the primary noise with an autoregressive moving average model. The resulting model allows for predicting future values of the primary noise. Finally, the primary noise values predicted are filtered by a non-causal inverse of the secondary path model to generate the anti-noise signal. Simulation results attest to the validity of the algorithm in reducing narrowband noise.


Introduction
Active noise and vibration control (ANVC) or simply active noise control (ANC) [11,12,15,17,23,24,33,42] is a class of methods to reduce noise and vibration by generating an anti-noise signal with a phase opposite to the original noise so that both interfere destructively, reducing the noise level.ANC works best at low frequencies where passive techniques are less effective and can be used together with the last.It is finding applications in many fields like noise reduction in cars [40], homes [8], headphones [10,25,39], and even vibration of airplane tails [35].
ANC methods can be divided into two categories [23], feedforward ANC and feedback ANC.In feedforward ANC, a reference sensor captures a reference signal B Paulo A. C. Lopes paulo.lopes@tecnico.ulisboa.pt 1 Instituto de Engenharia de Sistemas e Computadores, Investigação e Desenvolvimento, Instituto Superior Técnico, Universidade de Lisboa, INESC-ID/IST/UL, Rua Alves Redol n. 9, 1000-029 Lisbon, Portugal early in the noise path.Then a controller uses information from this sensor to generate an anti-noise signal that is fed to an actuator to reduce the noise.The residual noise level is measured by an error sensor and used by the controller when generating the anti-noise.In feedback ANC, there is no reference sensor, and the anti-noise signal is generated using only the information from the error sensor.This work focuses on feedback ANC since in many applications may be difficult to obtain a good reference.Feedback ANC has been extensively applied to ANC headphones, for instance, [22,41,43,44] and many other systems [21,35].
Most feedforward and feedback ANC algorithms (like the filtered-x least mean squares (FxLMS) algorithm [32]) require a model of the secondary path.This model may be obtained offline before the ANC controller starts to reduce noise or online during the operation of the ANC controller.Online secondary path modeling allows the ANC system to cope with slow changes in the secondary path.However, only a few algorithms can cope with fast changes.Eriksson [13] first proposed online secondary path modeling using a small auxiliary white noise signal as input to an adaptive filter.After this, it was proposed to improve the modeling by using an additional adaptive filter to remove the primary noise from the secondary path model error signal [5,45].The auxiliary noise power can be chosen so that residual noise to auxiliary noise at the error sensor is constant [9] and it can also be turned off after convergence has been reached [2,26] and on after the secondary path changes.Finally, it can also be used to model the path from the anti-noise transducer to the reference sensor (feedback path) [1].
In feedforward ANC, the secondary path can also be modeled without the auxiliary noise by the simultaneous equations method [14] or the overall modeling algorithm (OMA) [23].However, the secondary path model of the OMA algorithm may be incorrect [6,23,27].Nevertheless, even with an incorrect secondary path model, the resulting overall model (primary and secondary path) can be used to determine the optimum controller as in the mirror-modified FxLMS (MMFxLMS) algorithm [27,28].
There is very little work on online secondary path modeling for feedback ANC.However, when using auxiliary noise, the same algorithms as in feedforward ANC can be used [18,43].Unfortunately, to the authors' knowledge, there is no work on online secondary path modeling for feedback ANC without auxiliary noise.Moreover, the anti-noise signal and the (unknown) disturbance signal are usually highly correlated, making it difficult to use the first to obtain an accurate secondary path model.This work proposes to use a variant of adaptive model predictive control (MPC) [30,38] to determine the anti-noise signal.In MPC, the plant output (the ANC error signal) future values are predicted using a model and set close to a predetermined signal (zero in ANC) by carefully selecting the present and future values of the plant input signal (anti-noise).The predictions consider the effect that present inputs have on future outputs.The ANC's primary noise signal is the MPC disturbance.However, unlike standard model predictive control applications, in ANC, the plant plus disturbance model is unknown and needs to be estimated, making this an adaptive model predictive control application [30,38] that has yet to receive much attention.Nevertheless, the model is linear, allowing for simpler solutions, as shown in this work.ANC also falls in the more general field of adaptive control [4] and stochastic adaptive control [7] that, in general, requires complex solutions with caution and probing (a dual controller) [20].Adaptive model predictive control applications have been presented in [34] and [31] for ANC.This work builds on these two papers but presents a different approach with lower computational complexity.
Regarding notation, vectors are boldface lowercase letters as u.Given a signal u(n), u(n) is the vector with its past samples u(n) = [u(n) . . .u(n − N + 1)] T where N is the vector length, which should be determined from context.The operator * stands for convolution.

Proposed Algorithm
In feedback ANC, the residual noise signal e(n) is given by the sum of the primary noise signal d(n) (disturbance) with the anti-noise signal u(n) after passing through the secondary path (or plant) with impulse response s(n) (1) In the proposed algorithm, d(n) is modeled by an order N d autoregressive (AR) model: The signal d(n) is taken to have small power making d(n) predictable.This happens in narrowband systems, for instance.The secondary path is modeled by an autoregressive exogenous input (ARX) model so that the anti-noise signal y assuming the model coefficients are constant or slowly varying.Taking the z-transform [36] results in and taking the inverse z-transform result in where a = a s * a d and b = b s * a d .
Estimates for a and b, â, and b can be obtained by minimizing the residual noise estimation error square sum ξ These estimates can be obtained using the RLS algorithm [19,31] This model is the least-squares (LS) one-step-ahead predictor of e(n) given u(n), and it can be used to predict future values of e(n).Replacing the LS with the least mean squares (LMS) predictor (since they are equivalent for large M and stationary signals) and since the last is given by the expected value [3] results in the following.Using measures up to time n, the predicted residual noise at time n + 1 is The predicted residual noise at time n + 2 is etc. So, ê(n +i) given u(n) can be obtained simply by iterating the difference equation of the estimated model.Namely, where Note, however, that the model formed by âx and b is the LS predictor of e(n) and may not be the best predictor of e(n + i) even for i = 1.It can be tuned to a set of u(n) components that change in u(n + i), especially after fast changes.Regardless, this work will use this predictor in the remaining text. Let and Then, since (9) describes a linear IIR filter [36], then where i > 0, ŝ(i) is the impulse response of the IIR filter, and the signal ê0 (n + i) can be obtained from iterating the difference equation with ê(n + i) = e(n + i) for i ≤ 0. Finally, let Since the goal is to drive ê(n + i, n) to e 0 (n + i), keeping past values and driving future values to zero, then u 1 (n) should be where ŝ−1 (n) is the impulse response of the inverse of the IIR filter in (9).Let ŝ(z) = b(z)/ â(z) and ŝ−1 (z) = 1/ŝ(z) be the z-transform of ŝ(n) and ŝ−1 (n), then ŝ−1 (n) can be obtained using the inverse z-transform of ŝ−1 (z).To obtain a finite energy signal (and a stable filter), the region of convergence (ROC) of the z-transform should be selected to include the unit circle resulting in general in a non-causal signal (and filter).This (non-causal) filter can be implemented since ê0 (n) is known for past and future.Finally, u(n + 1) is set to u 1 (n + 1) and used as input to the plant at time n + 1.
In the proposed algorithm, the non-causal filter ŝ−1 (z) = â(z)/ b(z) is implemented by numerically calculating the roots of b(z), z i .Then forming bst (z) by the set of roots with |z i | < α with α slightly greater than one (inside the unit circle or close), and but (z) by the set of roots with |z i | ≥ α (outside the unit circle).Choosing α > 1 assures that canceling poles and zeros at the unit circle (that form oscillators to predict sinusoidal signals) stay in the same filtering operation.Then ê0 (n + i, n) − e 0 (n + i) is filtered by 1/ but (z) by flipping (time reversing, x(−n)) all the signals to make the filter stable; and then by â(z)/b st (z).All the signals extend from n − L + 1 to n + L (length 2L) with L large enough so that the transients of the filters are negligible at time n + 1.The source code for the proposed algorithm is available in [29].

Computational Complexity
The computational complexity of the proposed algorithm is approximately (for large N and L) 32N 2 for the system identification component, 4/3N 3 for the calculation of the b(z) polynomial roots by calculating the eigenvalues of the companion matrix [16,37] and 16L N for filtering and prediction.This results in a total of 32N 2 + 4/3N 3 + 16L N .As a comparison, Mohseni's algorithm [31] has a computational complexity of L h N 2 + (8/3 + 6)L 3  h , where L h is the horizon length.Using the same parameters as in the simulation section, f s = 2000 Hz, N = 16, L = 64, L h = 64 results in 60 MFLOPS for the proposed algorithm and 4577 MFLOPS for Mohseni's algorithm.
Note that it is still possible to reduce the computational complexity by noting that ê0 (z, n)−e 0 (z) = g(z)/ â(z) where g(z) is an order N polynomial related to the initial conditions for ê0 (n + i, n) and that u 1 (n + i) = − b−1 (n) * g(z) but this will not be further discussed in this work.

Simulation Results
This section presents a set of simulation results.In the simulations, the sampling frequency is f s = 2000 Hz.The plant is formed by an order N x = 6 IIR filter.The plant parameters a s and b s coefficients are initially set Gaussian random values with zero mean and variance 1/N 2 x and one, respectively, with a s first coefficient equal to one.The parameters change during the simulation according to random walk model (a s (n where the state noise r a and r b is formed by a set independent white noise signals with variance q r (n).The plant changes fast from 2 to 3 s with q r = 1/ f s and slowly for the remaining time with q r (n) = 1/(400 f s ).The filter poles and zeros are kept at a distance greater than β = 0.1 of the unit circle at any time n.Three sinusoidal signals with frequencies of 200, 400, and 600 Hz and amplitudes of 0.5, 1.3, and 0.3 plus white background noise with variance q v = 0.01 are added at the plant's output.The charts plot the result for 100 simulations of several algorithms.The anti-noise signal saturates at 10 and −10.Before ANC is started, at 0.5 s, the input of the plant is set to a unit-power white noise signal, and the model identification component of the algorithm is left running.The percentile charts show the percentile 10%, 25%, 50%, 75%, 90%, and 99% residual noise level at each time, given some information about its probability density function (PDF).The algorithms parameters are presented in Table 1.These are based on the values presented in the corresponding papers and tuned by trial and error.In Carini's algorithm, the step size μ was scaled by the factor μ SC to improve stability.
Figure 1 shows the residual noise level power versus time of three typical simulations of the proposed algorithm.As can be seen, the residual noise quickly settles to the minimum value after ANC is started, rises while the plant changes fast, and settles to the minimum again as desired.However, there are some narrow spikes and even divergence (not shown).Figure 2 shows the percentile chart for the 100 simulations.As can be seen, performance is not always as good as in Fig. 1.The residual noise power takes large values in less than 10% of the simulations.Narrow spikes contribute to a rise in the percentile curves, but they are not visible on the plots.The very high percentile 99% curve is primarily due to divergence in some simulations.Nevertheless, the residual noise takes low values (at each time instant) in more than 90% of the simulations.This is confirmed by the histogram of the residual noise power at the simulation end, as shown in Fig. 3.It shows that the final residual noise power in 10 out Curves for percentile 10%, 25%, 50% (solid), 75%, 90%, and 99% (gray) of 100 simulations is considerably larger than the minimum value (with 6 divergence cases), while in the remaining 90, it is close to the optimum.
Figures 4 and 5 show percentile plots of the residual noise power versus time of Carini's [9] and Lopes' [26] algorithms simulation results when used in feedback ANC.The reference signal equals the estimated primary noise signal [23].This signal equals the error signal minus the anti-noise signal filtered by the plant These algorithms have very low computational complexity.They deal well with slow plant changes and are intended for feedforward (not feedback) ANC systems.So, performance is not good and considerably worse than the proposed algorithm, as can be seen by comparing (for instance) the percentile 90% curve.

Conclusion
This work presents an algorithm for narrowband feedback ANC systems that improves the performance of state-of-the-art algorithms when dealing with fast secondary path changes while keeping a moderate computational complexity.However, performance is as good as desired in every case.This algorithm uses an estimated model to predict the primary noise and the same model to estimate a non-causal secondary path inverse that filters the predicted noise to calculate the anti-noise signal.
. The vectors b and âx have length N + 1 and N resulting in a order N model.For small d(n) (as in narrowband disturbances), persistent excitation, and N large enough, the minimization results in â = a and b = b.The resulting model has canceling poles and zeros at the zeros of a d (z) and uncontrollable states that are only excited by d(n) in the original model, which may be described by a zero input filter in the new model.

Fig. 1
Fig. 1 Three typical simulations residual noise power versus time