# Optimal and Adaptive Virtual Unidirectional Sound Source in

Authors: Dariusz Bismor
Journal: Advances in Acoustics and Vibration
 (2008)

---

## Abstract

One of the problems concerned with active noise control is the existence of acoustical feedback between the control value (“active” loudspeaker output) and the reference signal. Various experiments show that such feedback can seriously decrease effects of attenuation or even make the whole ANC system unstable. This paper presents a detailed analysis of one of possible approaches allowing to deal
with acoustical feedback, namely, virtual unidirectional sound source. With this method,two
loudspeakers are used together with control algorithm assuring that  the combined behaviour of the pair makes virtual propagation of sound only in one direction. Two different designs are presented
for the  application of active noise control in an acoustic duct: analytical (leading to fixed controller) and adaptive. The algorithm effectiveness in simulations and real experiments for both solutions is showed, discussed, and compared.


---

## Body

## 1. Introduction

Active noise control is mainly concerned with low-frequency sounds that cannot be suppressed in passive manner at reasonable
cost. Other prerequisite for using ANC is that it does not obstruct media flow
in ducts, for example, air conditioning ducts. Unfortunately, almost all duct
applications have acoustic feedback; a phenomenon that can heavily decrease
or even destroy the results of otherwise properly set-up control system
[1, 2].There are several methods for avoiding acoustic
feedback in duct applications [2–6]. The best one that avoids
acoustic feedback completely is to use nonacoustic reference sensor.
Unfortunately, this usually limits the application to cancel the tonal sounds
and does not allow for the  reduction of broadband noise. Another approach to
acoustic feedback is to internally process signal from the reference microphone
in a way allowing to compensate for this effect [2, 4–6]. This technique, usually called “feedback
neutralisation,” uses a model of acoustic feedback path to filter the
controller output and then subtract it from the reference signal. In the
simplest case, the model is obtained offline and does not account for time variation
of the feedback path. Moreover, the method requires a model of the acoustic
feedback path and therefore is computationally more expensive.The approach to acoustic feedback cancellation
presented in this paper is called virtual unidirectional sound source (VUSS).
The method, in authors opinion, has some advantages over well-known two
loudspeaker active noise control system in duct applications. First, it clearly
separates the subsystem responsible for neutralisation of undesired feedback
from active noise control part. Second, its adaptive version can also be
applied in case of varying acoustic path responses. Last, virtual
unidirectional sound source can also be applied in more general active control
of sound applications, not only in active noise control. The paper attempts to
explain the theory behind VUSS as well as to give some results showing its
excellent performance in case of ANC.The paper is organised as follows. Section2
summarises physical basics necessary to understand principles of VUSS operation
and its advantage over one-loudspeaker applications. In Section
3, the
idea of VUSS is explained together with  the  original detailed analysis of optimal
(fixed) design (Section 3.1) and adaptive approach
(Section 3.2).
Finally, Section 4 gives the results of simulations and real experiments of
feedforward active noise control system using virtual unidirectional sound
source.


## 2. Acoustical Phenomena in Ducts

The simplest
solution to acoustic feedback cancellation problem in case where omnidirectional
microphone must be used as the reference sensor is called feedback
neutralisation [2, 4–6]. With this approach, only one
loudspeaker is used, but the signal obtained from reference microphone is
processed in the way allowing for compensation of its influence. If no
compensation is performed, severe problems can result. In the simplest case of
attenuation of one tone, the pressure distribution upstream secondary source
[3] takes shape of
standing wave, presented with solid line on Figure 1. If, for particular
frequency of this tone, the reference microphone is placed in standing wave
nodal point, it will give measurement equal to zero regardless of waves' actual
magnitude.Figure 1
Acoustic pressure distribution in duct: one
loudspeaker case (solid) and two loudspeaker case (dashed).Another situation occurs when two loudspeakers are
used, as presented in Figure2. Assume first loudspeaker is placed in the duct at
position x=0 and the second at position x=l.
Denote pp+(x) as the complex pressure [3] at location x produced by primary source emitting tonal
sound:(1)pp+(x)=Ae−ikx, where A is the complex amplitude of the sound wave and
where acoustic wave number k is equal to k=ω/c0,
with ω being wave angular frequency and c0 being the sound speed in air.Figure 2
Acoustic duct with two active loudspeakers.Assume the two loudspeakers on Figure2 produce sound
waves propagating upstream the duct, marked as pa1−(x) and pa2−(x),
and sound waves propagating downstream the duct, marked as pa1+(x) and pa2+(x).
Assume also that  the sound wave produced by the primary source is propagating
downstream the duct, as marked on the figure.Assume that  the complex pressure of sound wave generated by
first “active” loudspeaker is given by(2)pa1+(x)=Be−ikxforx>0,pa1−(x)=Beikxforx<0,  and the complex
pressure of sound wave produced by the second “active” loudspeaker is given
by(3)pa2+(x)=Ce−ik(x−l)forx>l,pa2−(x)=Ceik(x−l)forx<l,  where B and C are the complex amplitudes of the sound waves.In the situation described above it is possible to set
another condition on the system regardless of primary sound source attenuation
requirement. In the case of unidirectional sound source, it may be the
requirement of self-cancellation of the secondary sound waves propagating
upstream the duct:(4)pa1−(x)+pa2−(x)=0forx<0.In this case, the secondary sources have to be driven
in the way providing [3](5)B=−Ce−ikl.In the same time to achieve the
main ANC goal, which is to cancel the primary sound downstream the second
loudspeaker completely, it is necessary to set the complex
amplitude:(6)C=−A2isin(kl). The acoustic pressure level distribution for two
loudspeakers case has been shown with dashed line on Figure 1. The
figure shows
that the sound pressure has been perfectly cancelled at coordinates x>l while the absolute value of the pressure level
has not been affected by the secondary sources at negative coordinates x.
The zone between the secondary loudspeakers (0<x<l) is a transient zone.Although two-loudspeaker system offers significant
improvement over one loudspeaker system, it has its drawbacks too.
Specifically, the above control law says that for some frequencies for whichsin(kl) is close to zero, both the secondary sources
will have to produce sound waves with amplitude very large, compared with the
amplitude A of the primary waveform. Consider, for
example, the case where the loudspeakers are separated by the distance l=0.3 m. In this case the sound waves of frequency
570 Hz and all its multiplies will result in sin(kl) equal to, or close to, zero, depending on the
actual speed of sound, and thus on temperature, humidity, and so on. This fact will
certainly influence the performance of ANC system trying to attenuate such
sound waves (see, e.g., the drop in attenuation just below 600 Hz on
Figure 10).


## 3. Virtual Unidirectional Sound Source

It is well
known that introducing the second “active” loudspeaker allows to deal with
many one-loudspeaker systems difficulties [4]. Moreover, the analysis outlined in
Section 2  proves that  it allows also for (hypothetically) perfect cancellation of the
acoustic feedback. One of possible systems using two secondary loudspeakers is
called virtual unidirectional sound source.The idea of virtual unidirectional sound source (VUSS)
is to use digital signal processing algorithm to drivetwo loudspeakers in such way that the
sound produced by them propagates only downstream the duct [7]. In fact, although the sound generated by each
loudspeaker propagates in both directions, the processing algorithm tries to
assure that the sound waves propagating upstream the duct are actively cancelled by
themselves while those propagating downstream the duct are amplified. The
additional advantage of this approach is that it is sometimes possible to
equalise the amplitude spectrum of the secondary path transfer function
(transfer function between signals t(i) and e(i) on Figure 3) that plays very important role in
active control of sound [4].Figure 3
Block diagram of virtual unidirectional sound source.The block diagram of VUSS has been shown in Figure3.
The system contains two “active” loudspeakers, L2a and L2b, as well as two
microphones: M1, generally called reference microphone, and M2 called error
microphone (the names of the microphones correspond to the role they play in
active noise control). The signal to be processed by VUSS is the set point
value (or antinoise signal) t(i) produced by ANC algorithm. It is filtered by
two filters with transfer functions, W1(z) and W2(z), and sent to the loudspeakers, L2a and L2b,
respectively. The error signals are obtained by comparing values acquired from
the microphones M1 and M2 with their required values (set points). As the goal
of VUSS system is to cancel the influence of the set point value t(i) at the microphone M1 point, the signal x(i) is compared with zero (e1(i)=x(i)−0=x(i)). Similarly, as the desired signal at the
microphone M2 point is the delayed set point value, the second error signal is
calculated as e2(i)=e(i)−to(i).
The need for delaying the set point value t(i) comes from condition on causality of the
system: no filtration can compensate for the delay time introduced by the
acoustic path between the loudspeakers and the M2 microphone. The value of the
delay ensuring causality of the algorithm should be chosen before running the
algorithm. The error signals e1(i) and e2(i) are used by adaptive algorithm to update
parameters of the transfer functions W1(z) and W2(z) (in case of optimal design they are used only
for observing the performance of the system).


### 3.1. Optimal and Suboptimal Filter Designs

In this
subsection, an attempt to derive analytical optimal and suboptimal filter
formulae will be made. The need for suboptimal formula will be explained and
feasibility of both the solutions will be discussed.


#### 3.1.1. Formulation Using Transfer Functions

Suppose that the
whole system presented in Figure3 is linear and the signal to be cancelled is
deterministic, or random wide-sense
stationary [8].
Suppose that  W1 and W2 are infinite impulse response (IIR) filters.
Figure 3 shows that the control values u1(i) and u2(i) can be expressed as a filtration of set point
value t(i) by the W1 and W2 filters. Using the notation of difference
equation, this filtration can be expressed as(7)u1(i)=W1(z)⋅t(i)=w1,0t(i)+w1,1t(i−1)+⋯=∑n=0∞w1,nt(i−n),u2(i)=W2(z)⋅t(i)=w2,0t(i)+w2,1t(i−1)+⋯=∑n=0∞w2,nt(i−n),where(8)W1(z)=w1,0+w1,1z−1+w1,2z−2+⋯,W2(z)=w2,0+w2,1z−1+w2,2z−2+⋯.In the above equation, w1,n and w2,n are nth elements of the W1 and W2 filter impulse response functions, and are
also called filter coefficients. Using this notation, the z−1 is treated as discrete shift operator so that z−nt(i)=t(i−n).From Figure3 it also follows that the reference signal x(i) due to the secondary sources measured by the
M1 microphone can be expressed with the following difference
equation:(9)x(i)=S11(z)⋅u1(i)+S12(z)⋅u2(i),where S11(z) and S12(z) are impulse response functions of
electroacoustical paths between the reference signal x(i) and the u1(i),u2(i) control values, including amplifiers, filters,
loudspeakers, and others.Similarly, the signale(i) measured by the microphone M2 due to the
secondary sources can be given by(10)e(i)=S21(z)⋅u1(i)+S22(z)⋅u2(i).Considering the above and the
definitions of error signals given in the previous subsection, the error
signals can be expressed as(11)e1(i)=S11⋅u1(i)+S12⋅u2(i),e2(i)=S21⋅u1(i)+S22⋅u2(i)−to(i),where the dependence of Snn transfer functions on z has been omitted to simplify the notation.The goal of the tuning algorithm is to drive error
signals to zero:(12)e1(i)=0,e2(i)=0.Substituting  (7) and
(11) into (12),
recognising that to(i)=z−Δt(i) gives for all t(i)≠0(13)S11⋅W1(z)+S12⋅W2(z)=0,S21⋅W1(z)+S22⋅W2(z)=z−Δ,where Δ samples delay time was introduced to assure
causality of the system.Solving the above equation set gives the transfer
functions of optimal filters:(14)W1opt(z)=−z−ΔS12S11⋅S22−S12⋅S21,W2opt(z)=z−ΔS11S11⋅S22−S12⋅S21.
Unfortunately, (14) cannot be used in practise due
to non-minimumphase nature of transfer functions Snn leading to instability of VUSS filters. Moreover, even when all the transfer functions Snn are minimumphase, the subtraction in
denominator of (14) can still result in instability of the whole filter.
This is referred to in the literature as “unconstrained controller”
[3].The suboptimal design of virtual unidirectional sound
source can be performed after omitting the second condition from (12) set
in mathematical derivations. In consequence, it means we do not expect the
transfer function from set point value t(i) to microphone M2 signal e(i) to be pure delay, but we allow it to take more
complex form. Furthermore, it means that the amplitude spectrum of this
transfer function will not be flat over the frequency range under
consideration. However, we still request cancellation of the secondary sources
influence on the reference signal x(i).With control goal stated above, there are two
equivalent solutions that will be called suboptimal, namely,(15)W1sub(z)=−S12(z),W2sub(z)=S11(z),or(16)W1sub(z)=S12(z),W2sub(z)=−S11(z).Moreover, assuming the Snn transfer functions are of FIR type, the
suboptimal solutions are also FIR type.It is, however, important to notice that in case of
suboptimal solution the secondary path transfer function (transfer function
between signalst(i) and e(i) on Figure 3) is given by(17)S(z)=−S12⋅S21+S11⋅S22,orS(z)=S12⋅S21−S11⋅S22.(Again, the dependence of Snn transfer functions on z has been omitted to clarify the notation.)
This shows that in case of no demand for equalisation of the secondary path
transfer function amplitude spectrum it can assume arbitrary shape, potentially
“worse” than in case of single secondary source. This effect must be taken
into account when estimating secondary path models required by active noise
control algorithm (e.g., by using longer filters).


#### 3.1.2. Formulation Using Reference Signal Matrix

The following
subsection presents result obtained using the approach presented in [3, 9], which is based on impulse
response of appropriate electroacoustic paths and on the  use of filtered-reference
signals. Using this approach, the lth error signal (l={1,2} in this case) can be expressed as a sum of
desired signals and contributions from m secondary sources (with m=2 in this case) as:(18)el(i)=−dl(i)+∑m=12∑j=0J−1slm(j)um(i−j),where slm(j) is the jth coefficient of impulse response of
electroacoustic path between the  lth sensor and the  mth actuator and dl(i) is desired value of the  lth sensor signal. The length of impulse
response can be arbitrary number J guaranteeing desired accuracy.In the same manner, the signals driving actuators can
be expressed as(19)um(i)=∑n=0N−1wm,nt(i−n),where wm,n is the nth parameter of the Nth order wm FIR filter impulse response.Substituting (19) into (18)
yields(20)el(i)=−dl(i)+∑m=12∑j=0J−1∑n=0N−1slm(j)wm,nt(i−j−n). To derive a matrix formulation for the above equation
it is necessary to reorganise the order of filtration of the set point value t(i).
To do this, it is necessary to assume that the filters slm and wm,n are time invariant. If the filtered reference
signal is defined as(21)rlm(i)=∑j=0J−1slm(j)t(i−j),the error signals can be written
as(22)el(i)=−dl(i)+∑m=12∑n=0N−1wm,nrlm(i−n).After defining(23)w(i)=[w1,iw2,i],rl(i)=[rl1(i)rl2(i)].Equation (22) can be expressed
as(24)el(i)=−dl(i)+∑n=0N−1wT(n)rl(i−n).Finally, the vector of two error
signals can be defined as(25)e(i)=[e1(i)e2(i)]Tand the vector of two desired
signals can be defined as(26)d(i)=[d1(i)d2(i)]T,leading to the conclusive
formula(27)e(i)=−d(i)+R(i)W,where(28)R(i)=[r1T(i)r1T(i−1)⋯r1T(i−N+1)r2T(i)r2T(i−1)⋯r2T(i−N+1)]is called matrix of filtered
reference signals, and(29)W=[wT(0)wT(1)…wT(N−1)]Tis the vector containing all the
coefficients of both W1 and W2 filters (the filters are time independent).It is interesting to notice that in case of virtual
unidirectional sound source the vector of desired signals can be instantiated
as(30)d(i)=[0t(i−Δ)].Following the methodology
presented in [9] the
optimal filter solution can be expressed as(31)Wopt=[E{RT(i)R(i)}]−1E{RT(i)d(i)},where E{⋅} denotes expectation operator.The matrix to be inverted is of dimensionN×N,
where N is the filter length. Fortunately, it can be
proved that  the matrix is block Toeplitz matrix [3], so effective iterative methods can be applied for
inversion. The expectation operator tells that statistical properties of
filtered set point value t(i) will be taken into account. The expression
under the right-hand side expectation operator, in the above equation, takes
particularly the simple form of(32)RT(i)d(i)=[r2T(i)r2T(i−1)⋯r2T(i−N+1)]Tt(i−Δ).


#### 3.1.3. Formulation in Frequency Domain

Applying the Fourier transform to (7) and
(11)
and substituting the transformed (7) into
(11)
yields(33)E1(eiωTp)=S11(eiωTp)W1(eiωTp)T(eiωTp)+S12(eiωTp)W2(eiωTp)T(eiωTp),E2(eiωTp)=S21(eiωTp)W1(eiωTp)T(eiωTp)+S22(eiωTp)W2(eiωTp)T(eiωTp)−To(eiωTp)with z=eiωTp,
where Tp is a sampling period.The above equations can be expressed in more
convenient matrix form as(34)e(eiωTp)=−d(eiωTp)+S(eiωTp)W(eiωTp)T(eiωTp),where e(eiωTp),d(eiωTp), and W(eiωTp) are transformed vectors defined by (25),
(30), and
(29), respectively, and the matrix of transfer functions is defined
as(35)S(ejωTp)=[S11(eiωTp)S12(eiωTp)S21(eiωTp)S22(eiωTp)]. The unconstrained controller can be found by
minimising the cost function equal to expectation of the sum of squared errors,
independently at each frequency [3]:(36)J=E{eHe}=traceE{eeH}.Thus, the optimal controller is
equal to(37)Wopt(eiωTp)=S−1(eiωTp)Std(eiωTp)Stt−1(eiωTp),where matrices of power and
cross spectral densities are defined as(38)Std(eiωTp)=E{d(eiωTp)T*(eiωTp)},Stt(eiωTp)=E{T(eiωTp)T*(eiωTp)}.
Assume now that the set point value t(i) used during VUSS learning phase is a sequence
of white noise with expected value σ2.
The power spectral density matrix Stt is then equal to(39)Stt(eiωTp)=σ2,independent on ω,
and the cross spectral density matrix Std is equal to(40)Std(eiωTp)=[0e−iωTpΔσ2].After evaluating inverse of
transfer function block matrix S(eiωTp), the optimal controller obtained with white
noise excitation becomes(41)Wopt(eiωTp)=e−iωTpΔS11S22−S12S21[−S12S11],where the dependence of Snn transfer function on eiωTp has been omitted to clarify the notation.The above equation is consistent with the result
obtained in Section3.1.1. It is, however, still the unconstrained form of
the controller which is not feasible for practical implementation.


### 3.2. Adaptive Filter Solution

This
subsection presents an alternative approach to optimal and suboptimal filter
designs discussed above. Now we will try to develop an adaptive algorithm with
the goal defined by set of (12). We will assume that W1(z) and W2(z) are in form of finite impulse response filters
so the filter coefficients can be stored in a vector similar to this defined by
(29), but now with coefficients varying slowly (compared to timescales of
plant dynamics) in time.Although there are many algorithms that can be used to
tune filter coefficients [4], two-channel filtered-x LMS algorithm was chosen for
the following derivation and experiments. Its advantages are simplicity and
robustness, even if speed of convergence is not the best. Using two-channel
FXLMS algorithm, the update equation is given by(42)W(i+1)=W(i)+μRT(i)e(i),where μ is step size and R(i) and e(i) are defined by (28) and (25),
respectively.The matrix of reference signalsR(i) is generated by filtration of the set point
value by electroacoustic path transfer functions S11(z) to S22(z).
In practise, the latter can be only estimated yielding a matrix of estimated
plant responses R^(i).
The estimation must take place before running FXLMS algorithm and is usually
done by separate identification procedure. The same procedure can also be
turned on during active noise control phase after detecting substantial changes
in plant dynamics (e.g., changes in air temperature).To study convergence properties of the above algorithm,
it is necessary to substitute (27) into (42) with estimated plant
responses matrix, producing(43)W(i+1)=W(i)+μ[−R^T(i)d(i)+R^T(i)R(i)W(i)].If the adaptive algorithm is stable,
it will converge to the solution setting the expectation value of the term in
bracket to zero [3].
The steady-state vector of filter coefficients will therefore be equal
to(44)W∞={E[R^(i)R(i)]}−1E[R^(i)d(i)]. The above result for the adaptive algorithm steady
state is in coincidence with the optimal solution presented in
Section 3.1.2
(31) if and only if the estimated matrix of reference signals R^ is equal to the true matrix of reference
signal. This validates, however, the methodology.The precise consequences of nonperfect matching of
true plant responses contained in the matrix of reference signals on adaptation
process are still unknown. The strongest result has been shown by Wang and Ren
[10], the
sufficient condition for stability of adaptation. Nevertheless it is not a necessary
condition and it has been obtained with the small step size assumption.During this research it was useful to introduce an
additional weight parameterβ∈(0,1) to specify which of the goals defined in
(12) should have more bearing on adaptation process. For such case the
definition of error vector e(i) was modified as follows:(45)e′(i)=[(1−β)e1(i)βe2(i)]=[(1−β)00β]e(i)=Be(i).When β is close to zero, the algorithm is better
excited along the modes responsible for neutralisation of the acoustic feedback
effect, while when β is close to one the algorithm puts more stress
on equalising the secondary path transfer function.Introducing the weight parameterβ modifies the steady-state vector of filter
coefficients giving(46)W∞={E[R^(i)BR(i)]}−1E[R^(i)Bd(i)].The maximum step size parameter, μ,
assuring convergence of the whole algorithm is given by(47)0<μ<2Re(λmax)|λmax|2,where λmax denotes maximum eigenvalue of R^(i)BR(i) matrix rather than R^(i)R(i) as was in the  original solution [3]. The weight parameter β can therefore influence, to some degree, the
eigenvalues of the matrix under consideration.


## 3.1. Optimal and Suboptimal Filter Designs

In this
subsection, an attempt to derive analytical optimal and suboptimal filter
formulae will be made. The need for suboptimal formula will be explained and
feasibility of both the solutions will be discussed.


### 3.1.1. Formulation Using Transfer Functions

Suppose that the
whole system presented in Figure3 is linear and the signal to be cancelled is
deterministic, or random wide-sense
stationary [8].
Suppose that  W1 and W2 are infinite impulse response (IIR) filters.
Figure 3 shows that the control values u1(i) and u2(i) can be expressed as a filtration of set point
value t(i) by the W1 and W2 filters. Using the notation of difference
equation, this filtration can be expressed as(7)u1(i)=W1(z)⋅t(i)=w1,0t(i)+w1,1t(i−1)+⋯=∑n=0∞w1,nt(i−n),u2(i)=W2(z)⋅t(i)=w2,0t(i)+w2,1t(i−1)+⋯=∑n=0∞w2,nt(i−n),where(8)W1(z)=w1,0+w1,1z−1+w1,2z−2+⋯,W2(z)=w2,0+w2,1z−1+w2,2z−2+⋯.In the above equation, w1,n and w2,n are nth elements of the W1 and W2 filter impulse response functions, and are
also called filter coefficients. Using this notation, the z−1 is treated as discrete shift operator so that z−nt(i)=t(i−n).From Figure3 it also follows that the reference signal x(i) due to the secondary sources measured by the
M1 microphone can be expressed with the following difference
equation:(9)x(i)=S11(z)⋅u1(i)+S12(z)⋅u2(i),where S11(z) and S12(z) are impulse response functions of
electroacoustical paths between the reference signal x(i) and the u1(i),u2(i) control values, including amplifiers, filters,
loudspeakers, and others.Similarly, the signale(i) measured by the microphone M2 due to the
secondary sources can be given by(10)e(i)=S21(z)⋅u1(i)+S22(z)⋅u2(i).Considering the above and the
definitions of error signals given in the previous subsection, the error
signals can be expressed as(11)e1(i)=S11⋅u1(i)+S12⋅u2(i),e2(i)=S21⋅u1(i)+S22⋅u2(i)−to(i),where the dependence of Snn transfer functions on z has been omitted to simplify the notation.The goal of the tuning algorithm is to drive error
signals to zero:(12)e1(i)=0,e2(i)=0.Substituting  (7) and
(11) into (12),
recognising that to(i)=z−Δt(i) gives for all t(i)≠0(13)S11⋅W1(z)+S12⋅W2(z)=0,S21⋅W1(z)+S22⋅W2(z)=z−Δ,where Δ samples delay time was introduced to assure
causality of the system.Solving the above equation set gives the transfer
functions of optimal filters:(14)W1opt(z)=−z−ΔS12S11⋅S22−S12⋅S21,W2opt(z)=z−ΔS11S11⋅S22−S12⋅S21.
Unfortunately, (14) cannot be used in practise due
to non-minimumphase nature of transfer functions Snn leading to instability of VUSS filters. Moreover, even when all the transfer functions Snn are minimumphase, the subtraction in
denominator of (14) can still result in instability of the whole filter.
This is referred to in the literature as “unconstrained controller”
[3].The suboptimal design of virtual unidirectional sound
source can be performed after omitting the second condition from (12) set
in mathematical derivations. In consequence, it means we do not expect the
transfer function from set point value t(i) to microphone M2 signal e(i) to be pure delay, but we allow it to take more
complex form. Furthermore, it means that the amplitude spectrum of this
transfer function will not be flat over the frequency range under
consideration. However, we still request cancellation of the secondary sources
influence on the reference signal x(i).With control goal stated above, there are two
equivalent solutions that will be called suboptimal, namely,(15)W1sub(z)=−S12(z),W2sub(z)=S11(z),or(16)W1sub(z)=S12(z),W2sub(z)=−S11(z).Moreover, assuming the Snn transfer functions are of FIR type, the
suboptimal solutions are also FIR type.It is, however, important to notice that in case of
suboptimal solution the secondary path transfer function (transfer function
between signalst(i) and e(i) on Figure 3) is given by(17)S(z)=−S12⋅S21+S11⋅S22,orS(z)=S12⋅S21−S11⋅S22.(Again, the dependence of Snn transfer functions on z has been omitted to clarify the notation.)
This shows that in case of no demand for equalisation of the secondary path
transfer function amplitude spectrum it can assume arbitrary shape, potentially
“worse” than in case of single secondary source. This effect must be taken
into account when estimating secondary path models required by active noise
control algorithm (e.g., by using longer filters).


### 3.1.2. Formulation Using Reference Signal Matrix

The following
subsection presents result obtained using the approach presented in [3, 9], which is based on impulse
response of appropriate electroacoustic paths and on the  use of filtered-reference
signals. Using this approach, the lth error signal (l={1,2} in this case) can be expressed as a sum of
desired signals and contributions from m secondary sources (with m=2 in this case) as:(18)el(i)=−dl(i)+∑m=12∑j=0J−1slm(j)um(i−j),where slm(j) is the jth coefficient of impulse response of
electroacoustic path between the  lth sensor and the  mth actuator and dl(i) is desired value of the  lth sensor signal. The length of impulse
response can be arbitrary number J guaranteeing desired accuracy.In the same manner, the signals driving actuators can
be expressed as(19)um(i)=∑n=0N−1wm,nt(i−n),where wm,n is the nth parameter of the Nth order wm FIR filter impulse response.Substituting (19) into (18)
yields(20)el(i)=−dl(i)+∑m=12∑j=0J−1∑n=0N−1slm(j)wm,nt(i−j−n). To derive a matrix formulation for the above equation
it is necessary to reorganise the order of filtration of the set point value t(i).
To do this, it is necessary to assume that the filters slm and wm,n are time invariant. If the filtered reference
signal is defined as(21)rlm(i)=∑j=0J−1slm(j)t(i−j),the error signals can be written
as(22)el(i)=−dl(i)+∑m=12∑n=0N−1wm,nrlm(i−n).After defining(23)w(i)=[w1,iw2,i],rl(i)=[rl1(i)rl2(i)].Equation (22) can be expressed
as(24)el(i)=−dl(i)+∑n=0N−1wT(n)rl(i−n).Finally, the vector of two error
signals can be defined as(25)e(i)=[e1(i)e2(i)]Tand the vector of two desired
signals can be defined as(26)d(i)=[d1(i)d2(i)]T,leading to the conclusive
formula(27)e(i)=−d(i)+R(i)W,where(28)R(i)=[r1T(i)r1T(i−1)⋯r1T(i−N+1)r2T(i)r2T(i−1)⋯r2T(i−N+1)]is called matrix of filtered
reference signals, and(29)W=[wT(0)wT(1)…wT(N−1)]Tis the vector containing all the
coefficients of both W1 and W2 filters (the filters are time independent).It is interesting to notice that in case of virtual
unidirectional sound source the vector of desired signals can be instantiated
as(30)d(i)=[0t(i−Δ)].Following the methodology
presented in [9] the
optimal filter solution can be expressed as(31)Wopt=[E{RT(i)R(i)}]−1E{RT(i)d(i)},where E{⋅} denotes expectation operator.The matrix to be inverted is of dimensionN×N,
where N is the filter length. Fortunately, it can be
proved that  the matrix is block Toeplitz matrix [3], so effective iterative methods can be applied for
inversion. The expectation operator tells that statistical properties of
filtered set point value t(i) will be taken into account. The expression
under the right-hand side expectation operator, in the above equation, takes
particularly the simple form of(32)RT(i)d(i)=[r2T(i)r2T(i−1)⋯r2T(i−N+1)]Tt(i−Δ).


### 3.1.3. Formulation in Frequency Domain

Applying the Fourier transform to (7) and
(11)
and substituting the transformed (7) into
(11)
yields(33)E1(eiωTp)=S11(eiωTp)W1(eiωTp)T(eiωTp)+S12(eiωTp)W2(eiωTp)T(eiωTp),E2(eiωTp)=S21(eiωTp)W1(eiωTp)T(eiωTp)+S22(eiωTp)W2(eiωTp)T(eiωTp)−To(eiωTp)with z=eiωTp,
where Tp is a sampling period.The above equations can be expressed in more
convenient matrix form as(34)e(eiωTp)=−d(eiωTp)+S(eiωTp)W(eiωTp)T(eiωTp),where e(eiωTp),d(eiωTp), and W(eiωTp) are transformed vectors defined by (25),
(30), and
(29), respectively, and the matrix of transfer functions is defined
as(35)S(ejωTp)=[S11(eiωTp)S12(eiωTp)S21(eiωTp)S22(eiωTp)]. The unconstrained controller can be found by
minimising the cost function equal to expectation of the sum of squared errors,
independently at each frequency [3]:(36)J=E{eHe}=traceE{eeH}.Thus, the optimal controller is
equal to(37)Wopt(eiωTp)=S−1(eiωTp)Std(eiωTp)Stt−1(eiωTp),where matrices of power and
cross spectral densities are defined as(38)Std(eiωTp)=E{d(eiωTp)T*(eiωTp)},Stt(eiωTp)=E{T(eiωTp)T*(eiωTp)}.
Assume now that the set point value t(i) used during VUSS learning phase is a sequence
of white noise with expected value σ2.
The power spectral density matrix Stt is then equal to(39)Stt(eiωTp)=σ2,independent on ω,
and the cross spectral density matrix Std is equal to(40)Std(eiωTp)=[0e−iωTpΔσ2].After evaluating inverse of
transfer function block matrix S(eiωTp), the optimal controller obtained with white
noise excitation becomes(41)Wopt(eiωTp)=e−iωTpΔS11S22−S12S21[−S12S11],where the dependence of Snn transfer function on eiωTp has been omitted to clarify the notation.The above equation is consistent with the result
obtained in Section3.1.1. It is, however, still the unconstrained form of
the controller which is not feasible for practical implementation.


## 3.1.1. Formulation Using Transfer Functions

Suppose that the
whole system presented in Figure3 is linear and the signal to be cancelled is
deterministic, or random wide-sense
stationary [8].
Suppose that  W1 and W2 are infinite impulse response (IIR) filters.
Figure 3 shows that the control values u1(i) and u2(i) can be expressed as a filtration of set point
value t(i) by the W1 and W2 filters. Using the notation of difference
equation, this filtration can be expressed as(7)u1(i)=W1(z)⋅t(i)=w1,0t(i)+w1,1t(i−1)+⋯=∑n=0∞w1,nt(i−n),u2(i)=W2(z)⋅t(i)=w2,0t(i)+w2,1t(i−1)+⋯=∑n=0∞w2,nt(i−n),where(8)W1(z)=w1,0+w1,1z−1+w1,2z−2+⋯,W2(z)=w2,0+w2,1z−1+w2,2z−2+⋯.In the above equation, w1,n and w2,n are nth elements of the W1 and W2 filter impulse response functions, and are
also called filter coefficients. Using this notation, the z−1 is treated as discrete shift operator so that z−nt(i)=t(i−n).From Figure3 it also follows that the reference signal x(i) due to the secondary sources measured by the
M1 microphone can be expressed with the following difference
equation:(9)x(i)=S11(z)⋅u1(i)+S12(z)⋅u2(i),where S11(z) and S12(z) are impulse response functions of
electroacoustical paths between the reference signal x(i) and the u1(i),u2(i) control values, including amplifiers, filters,
loudspeakers, and others.Similarly, the signale(i) measured by the microphone M2 due to the
secondary sources can be given by(10)e(i)=S21(z)⋅u1(i)+S22(z)⋅u2(i).Considering the above and the
definitions of error signals given in the previous subsection, the error
signals can be expressed as(11)e1(i)=S11⋅u1(i)+S12⋅u2(i),e2(i)=S21⋅u1(i)+S22⋅u2(i)−to(i),where the dependence of Snn transfer functions on z has been omitted to simplify the notation.The goal of the tuning algorithm is to drive error
signals to zero:(12)e1(i)=0,e2(i)=0.Substituting  (7) and
(11) into (12),
recognising that to(i)=z−Δt(i) gives for all t(i)≠0(13)S11⋅W1(z)+S12⋅W2(z)=0,S21⋅W1(z)+S22⋅W2(z)=z−Δ,where Δ samples delay time was introduced to assure
causality of the system.Solving the above equation set gives the transfer
functions of optimal filters:(14)W1opt(z)=−z−ΔS12S11⋅S22−S12⋅S21,W2opt(z)=z−ΔS11S11⋅S22−S12⋅S21.
Unfortunately, (14) cannot be used in practise due
to non-minimumphase nature of transfer functions Snn leading to instability of VUSS filters. Moreover, even when all the transfer functions Snn are minimumphase, the subtraction in
denominator of (14) can still result in instability of the whole filter.
This is referred to in the literature as “unconstrained controller”
[3].The suboptimal design of virtual unidirectional sound
source can be performed after omitting the second condition from (12) set
in mathematical derivations. In consequence, it means we do not expect the
transfer function from set point value t(i) to microphone M2 signal e(i) to be pure delay, but we allow it to take more
complex form. Furthermore, it means that the amplitude spectrum of this
transfer function will not be flat over the frequency range under
consideration. However, we still request cancellation of the secondary sources
influence on the reference signal x(i).With control goal stated above, there are two
equivalent solutions that will be called suboptimal, namely,(15)W1sub(z)=−S12(z),W2sub(z)=S11(z),or(16)W1sub(z)=S12(z),W2sub(z)=−S11(z).Moreover, assuming the Snn transfer functions are of FIR type, the
suboptimal solutions are also FIR type.It is, however, important to notice that in case of
suboptimal solution the secondary path transfer function (transfer function
between signalst(i) and e(i) on Figure 3) is given by(17)S(z)=−S12⋅S21+S11⋅S22,orS(z)=S12⋅S21−S11⋅S22.(Again, the dependence of Snn transfer functions on z has been omitted to clarify the notation.)
This shows that in case of no demand for equalisation of the secondary path
transfer function amplitude spectrum it can assume arbitrary shape, potentially
“worse” than in case of single secondary source. This effect must be taken
into account when estimating secondary path models required by active noise
control algorithm (e.g., by using longer filters).


## 3.1.2. Formulation Using Reference Signal Matrix

The following
subsection presents result obtained using the approach presented in [3, 9], which is based on impulse
response of appropriate electroacoustic paths and on the  use of filtered-reference
signals. Using this approach, the lth error signal (l={1,2} in this case) can be expressed as a sum of
desired signals and contributions from m secondary sources (with m=2 in this case) as:(18)el(i)=−dl(i)+∑m=12∑j=0J−1slm(j)um(i−j),where slm(j) is the jth coefficient of impulse response of
electroacoustic path between the  lth sensor and the  mth actuator and dl(i) is desired value of the  lth sensor signal. The length of impulse
response can be arbitrary number J guaranteeing desired accuracy.In the same manner, the signals driving actuators can
be expressed as(19)um(i)=∑n=0N−1wm,nt(i−n),where wm,n is the nth parameter of the Nth order wm FIR filter impulse response.Substituting (19) into (18)
yields(20)el(i)=−dl(i)+∑m=12∑j=0J−1∑n=0N−1slm(j)wm,nt(i−j−n). To derive a matrix formulation for the above equation
it is necessary to reorganise the order of filtration of the set point value t(i).
To do this, it is necessary to assume that the filters slm and wm,n are time invariant. If the filtered reference
signal is defined as(21)rlm(i)=∑j=0J−1slm(j)t(i−j),the error signals can be written
as(22)el(i)=−dl(i)+∑m=12∑n=0N−1wm,nrlm(i−n).After defining(23)w(i)=[w1,iw2,i],rl(i)=[rl1(i)rl2(i)].Equation (22) can be expressed
as(24)el(i)=−dl(i)+∑n=0N−1wT(n)rl(i−n).Finally, the vector of two error
signals can be defined as(25)e(i)=[e1(i)e2(i)]Tand the vector of two desired
signals can be defined as(26)d(i)=[d1(i)d2(i)]T,leading to the conclusive
formula(27)e(i)=−d(i)+R(i)W,where(28)R(i)=[r1T(i)r1T(i−1)⋯r1T(i−N+1)r2T(i)r2T(i−1)⋯r2T(i−N+1)]is called matrix of filtered
reference signals, and(29)W=[wT(0)wT(1)…wT(N−1)]Tis the vector containing all the
coefficients of both W1 and W2 filters (the filters are time independent).It is interesting to notice that in case of virtual
unidirectional sound source the vector of desired signals can be instantiated
as(30)d(i)=[0t(i−Δ)].Following the methodology
presented in [9] the
optimal filter solution can be expressed as(31)Wopt=[E{RT(i)R(i)}]−1E{RT(i)d(i)},where E{⋅} denotes expectation operator.The matrix to be inverted is of dimensionN×N,
where N is the filter length. Fortunately, it can be
proved that  the matrix is block Toeplitz matrix [3], so effective iterative methods can be applied for
inversion. The expectation operator tells that statistical properties of
filtered set point value t(i) will be taken into account. The expression
under the right-hand side expectation operator, in the above equation, takes
particularly the simple form of(32)RT(i)d(i)=[r2T(i)r2T(i−1)⋯r2T(i−N+1)]Tt(i−Δ).


## 3.1.3. Formulation in Frequency Domain

Applying the Fourier transform to (7) and
(11)
and substituting the transformed (7) into
(11)
yields(33)E1(eiωTp)=S11(eiωTp)W1(eiωTp)T(eiωTp)+S12(eiωTp)W2(eiωTp)T(eiωTp),E2(eiωTp)=S21(eiωTp)W1(eiωTp)T(eiωTp)+S22(eiωTp)W2(eiωTp)T(eiωTp)−To(eiωTp)with z=eiωTp,
where Tp is a sampling period.The above equations can be expressed in more
convenient matrix form as(34)e(eiωTp)=−d(eiωTp)+S(eiωTp)W(eiωTp)T(eiωTp),where e(eiωTp),d(eiωTp), and W(eiωTp) are transformed vectors defined by (25),
(30), and
(29), respectively, and the matrix of transfer functions is defined
as(35)S(ejωTp)=[S11(eiωTp)S12(eiωTp)S21(eiωTp)S22(eiωTp)]. The unconstrained controller can be found by
minimising the cost function equal to expectation of the sum of squared errors,
independently at each frequency [3]:(36)J=E{eHe}=traceE{eeH}.Thus, the optimal controller is
equal to(37)Wopt(eiωTp)=S−1(eiωTp)Std(eiωTp)Stt−1(eiωTp),where matrices of power and
cross spectral densities are defined as(38)Std(eiωTp)=E{d(eiωTp)T*(eiωTp)},Stt(eiωTp)=E{T(eiωTp)T*(eiωTp)}.
Assume now that the set point value t(i) used during VUSS learning phase is a sequence
of white noise with expected value σ2.
The power spectral density matrix Stt is then equal to(39)Stt(eiωTp)=σ2,independent on ω,
and the cross spectral density matrix Std is equal to(40)Std(eiωTp)=[0e−iωTpΔσ2].After evaluating inverse of
transfer function block matrix S(eiωTp), the optimal controller obtained with white
noise excitation becomes(41)Wopt(eiωTp)=e−iωTpΔS11S22−S12S21[−S12S11],where the dependence of Snn transfer function on eiωTp has been omitted to clarify the notation.The above equation is consistent with the result
obtained in Section3.1.1. It is, however, still the unconstrained form of
the controller which is not feasible for practical implementation.


## 3.2. Adaptive Filter Solution

This
subsection presents an alternative approach to optimal and suboptimal filter
designs discussed above. Now we will try to develop an adaptive algorithm with
the goal defined by set of (12). We will assume that W1(z) and W2(z) are in form of finite impulse response filters
so the filter coefficients can be stored in a vector similar to this defined by
(29), but now with coefficients varying slowly (compared to timescales of
plant dynamics) in time.Although there are many algorithms that can be used to
tune filter coefficients [4], two-channel filtered-x LMS algorithm was chosen for
the following derivation and experiments. Its advantages are simplicity and
robustness, even if speed of convergence is not the best. Using two-channel
FXLMS algorithm, the update equation is given by(42)W(i+1)=W(i)+μRT(i)e(i),where μ is step size and R(i) and e(i) are defined by (28) and (25),
respectively.The matrix of reference signalsR(i) is generated by filtration of the set point
value by electroacoustic path transfer functions S11(z) to S22(z).
In practise, the latter can be only estimated yielding a matrix of estimated
plant responses R^(i).
The estimation must take place before running FXLMS algorithm and is usually
done by separate identification procedure. The same procedure can also be
turned on during active noise control phase after detecting substantial changes
in plant dynamics (e.g., changes in air temperature).To study convergence properties of the above algorithm,
it is necessary to substitute (27) into (42) with estimated plant
responses matrix, producing(43)W(i+1)=W(i)+μ[−R^T(i)d(i)+R^T(i)R(i)W(i)].If the adaptive algorithm is stable,
it will converge to the solution setting the expectation value of the term in
bracket to zero [3].
The steady-state vector of filter coefficients will therefore be equal
to(44)W∞={E[R^(i)R(i)]}−1E[R^(i)d(i)]. The above result for the adaptive algorithm steady
state is in coincidence with the optimal solution presented in
Section 3.1.2
(31) if and only if the estimated matrix of reference signals R^ is equal to the true matrix of reference
signal. This validates, however, the methodology.The precise consequences of nonperfect matching of
true plant responses contained in the matrix of reference signals on adaptation
process are still unknown. The strongest result has been shown by Wang and Ren
[10], the
sufficient condition for stability of adaptation. Nevertheless it is not a necessary
condition and it has been obtained with the small step size assumption.During this research it was useful to introduce an
additional weight parameterβ∈(0,1) to specify which of the goals defined in
(12) should have more bearing on adaptation process. For such case the
definition of error vector e(i) was modified as follows:(45)e′(i)=[(1−β)e1(i)βe2(i)]=[(1−β)00β]e(i)=Be(i).When β is close to zero, the algorithm is better
excited along the modes responsible for neutralisation of the acoustic feedback
effect, while when β is close to one the algorithm puts more stress
on equalising the secondary path transfer function.Introducing the weight parameterβ modifies the steady-state vector of filter
coefficients giving(46)W∞={E[R^(i)BR(i)]}−1E[R^(i)Bd(i)].The maximum step size parameter, μ,
assuring convergence of the whole algorithm is given by(47)0<μ<2Re(λmax)|λmax|2,where λmax denotes maximum eigenvalue of R^(i)BR(i) matrix rather than R^(i)R(i) as was in the  original solution [3]. The weight parameter β can therefore influence, to some degree, the
eigenvalues of the matrix under consideration.


## 4. Active Noise Control with VUSS

### 4.1. Active Noise Control Algorithms

Feedforward
active noise control system using virtual unidirectional sound source has been
shown in Figure4. This system uses a finite impulse response filter as control
filter with different adaptation algorithms. As in case of all feedforward
algorithms, only the reference value x(i) is used by the control filter F(z) to produce the set point value t(i).
The set point value is then processed by VUSS filters W1(z) and W2(z) to give two control values u1(i) and u2(i) that are amplified and sent to the
loudspeakers (via amplifiers and reconstruction filters not shown on the
figure). The adaptation algorithm on the other hand uses the reference signal x(i) as well as the error signal e(i) to tune the control filter coefficients.Figure 4
Active noise control system using virtual unidirectional sound source.The duct used in experiments described below was made
out of wood. It was 4 m long, with0.2×0.4 m rectangular section. One of the
duct ends was terminated with noise generating loudspeaker, while the other was
opened. The attenuating loudspeakers were of 0.16 m diameter and were situated
approximately in the middle of the duct. The distance between the middles of
the loudspeakers was equal to 0.3 m. The reference microphone was located 1.23 m from the middle of L2a loudspeaker, while the error microphone was separated
from L2b loudspeaker by the distance of 0.23 m.It should be emphasised that the goal of active noise
control algorithm is in contradiction to the goal of the virtual unidirectional
sound source adaptation algorithm [11]. The former tries to assure the compensation of sound
waves from both the primary and secondary sources at the microphone M2 point
without taking any notice of what happens at the microphone M1 point. The goal
of the latter is to compensate sound waves from the secondary sources only at
the microphone M1 point and try to assure the set point value appearing without
any alternation at the microphone M2 point. This leads to a conclusion that
both the algorithms should never operate at the same time. Indeed, the
experiments show that when both the adaptation algorithms are in operation the
whole system goes unstable. Therefore, the VUSS adaptation algorithm is usually
run on the beginning of experiments and switched off after adaptation is
completed. Next, only ANC algorithm is in use. Only when significant changes in
the environment are detected (e.g., temperature change above some level), the ANC
adaptation is temporarily frozen and the VUSS is retuned.FXLMS Algorithm
The first of
ANC filter adaptation algorithms tested was Filtered-x LMS algorithm [4, 9]. Assume that  the F(z) is an FIR filter with coefficients (see Figure 4)
in the  ith sample of time stored in a form of a vector f(i).
In that case, the FXLMS algorithm update equation is given by(48)f(i+1)=f(i)+μ⋅x′(i)⋅e(i),where μ is the step size and x′(i) is the vector of reference signal x(i) samples filtered by the secondary path
transfer function estimate, that is, an estimate of the transfer function
between set point value signal t(i) and error signal e(i).
If VUSS algorithm was working perfectly, the secondary
path transfer function would be equal to simple time delay. Unfortunately, the
study from  Section3.1 leads to a conclusion that such perfect situation is
impossible: the VUSS algorithm has to “pay more attention” to neutralising
acoustic feedback and secondary path equalisation is only its additional task.RLS Algorithm
The
second algorithm  of
ANC filter adaptation algorithms tested was recursive least squares (RLS)
algorithm (see, e.g., [9]). RLS algorithm uses very similar update equation in
the form(49)f(i+1)=f(i)+k(i)⋅e(i),where k(i) is the gain vector showing how much the value
of e(i) will modify different filter coefficients. The
gain vector is calculated as(50)k(i)=P(i)⋅x(i),where P(i) is a matrix updated in each step according to
the following equation:(51)P(i)=1λ(P(i−1)−P(i−1)x(i−1)xT(i−1)P(i−1)λ+xT(i−1)P(i−1)x(i−1)).
The RLS algorithm usually converges faster than the
LMS algorithm [12].
Because RLS algorithm needs to updateP matrix of size equal to the number of filter
coefficients in each adaptation step, it requires computational effort of N2,
whereas LMS algorithm requires computational effort of N only [13]. There is, however, a family of RLS algorithm
implementations called fast RLS algorithms that allows to omit computation of P matrix and use a selection of row vector
instead [13, 14]. This work used one of such
fast algorithms, called fast
transversal filter [15]. The algorithm was parametrised as follows. The
initial value of the minimum sum of backward a posteriori prediction-error
squares was equal to 1 and the exponential weighting factor was equal to
0.9999. The algorithm showed no stability issues, as it was expected on
floating-point arithmetic platform.
Both the FXLMS and RLS algorithms were implemented in
C programming language. Testing and debugging were performed using developed
simulation platform for Linux operating system. The platform allowed to emulate
DSP processor board behaviour and therefore the next step, moving the code
onto Texas Instruments TMS320C31 processor board, was purely automatic.
Another PC-computer program was used to tune various algorithm parameters and
to acquire the results. The sampling frequency of 2 kHz was chosen as frequency
band up to 1 kHz was of the author interest. It appeared that, due to hardware
limitations, appliable filter lengths were up to 70 with this sampling
frequency.


### 4.2. Testing Sounds

The set of signals chosen for experiments is presented
with dashed line on
Figures12–17. Testing signals N1 and N2
were generated offline as white noise filtered with high-order bandpass
filter. Bandpass was 200–300 Hz in case of N1 signal and 700–800 Hz in case of
signal N2. The former was chosen to show low-frequency attenuation capabilities
and the latter to show high frequency attenuation. The N3 is a signal acquired
in close vicinity of a food processor. It has the dominant frequency of about
270 Hz, but with substantial amount of broadband noise. The N4 is a signal
collected near a power transformer. It has all the harmonics of 50 Hz visible
with 250 Hz being the dominant. The N5 signal was recorded in water power plant
turbine proximity. It has the harmonics of 100 Hz distinguishable among
broadband noise. The last signal N6 was recorded in small bureau with an
electric propeller fan turned on. It is similar to N5, but the dominant
frequency is higher.All the signals were filed using 1000 samples. This
ensemble was repeated many times and (after amplification) sent to loudspeaker
L2 (see Figure4). The results of attenuation were measured with microphone M2
and calculated as(52)𝒯=10log(MSVofe(i)withoutattenuationMSVofe(i)withattenuation)[dB].


### 4.3. Simulation Results

The simulations
were performed using 150th order FIR filters models of the duct paths,
identified offline using methodology presented in [16]. The algorithm for
suboptimal VUSS filter solution described in Section 3.1 was tested first. As
expected, it allowed for very effective acoustical feedback cancellation but at
the cost of the secondary path being substantially degenerated. The amplitude
spectrum of transfer function between set point value t(i) and both the reference signal x(i) and error signal e(i) is presented on Figure 5.Figure 5
Amplitude spectrum of
transfer function between set point valuet(i) and x(i) (solid) and e(i) (dotted).After tuning phase, the system was engaged to active
noise control algorithms: FXLMS and Fast RLS. The ANC filter having 70
coefficients proved to be long enough to cancel out the noise in case of tonal
sounds, so this was the chosen value. However, in case of N1–N6 signals
attenuation was only (as for simulations of feedforward controller) 6–16 dB
(with FXLMS algorithm, which proved to be better), see
Table1. The examples
of spectra of the signals observed at the error microphone before and after
attenuation have been presented on Figures 6 and
7.Table 1
The
attenuation in simulations.
SignalAnalytical solution + FXLMSAnalytical solution + FXLMSAnalytical solution + RLSN19 dB15 dB16 dBN216 dB30 dB33 dBN310 dB15 dB16 dBN415 dB22 dB22 dBN59 dB17 dB18 dBN66 dB9 dB8 dBFigure 6
Noise N2
attenuation using analytical VUSS solution. —
Spectrum with ANC, - - -  spectrum
without ANC.Figure 7
Noise N5 attenuation using analytical VUSS
solution. — Spectrum with ANC, - - -  spectrum without ANC.In the following experiments, adaptive algorithm described
in Section3.2 was responsible for VUSS filters tuning. The value of the β parameter (see (45)) has been chosen as
0.5 to put the same effort into satisfying both the goals defined by (12).
The attenuation of tonal sounds was also infinite in case of adaptive VUSS
filters tuning. But the attenuation of testing signals was slightly better:
9–30 dB in case of FXLMS algorithm and 8–33 dB in case of RLS algorithm, see
Table 1. The examples of spectra of the signals observed at the error
microphone before and after attenuation have been shown on
Figures 8 and  9. The
difference in amplitude spectrum of the N2 and N5 signals without attenuation
visible on these figures and on Figure 6 and
7 comes from the fact that this is
the spectrum measured with the M2 microphone, not the spectrum of the signal
driving the L1 loudspeaker.Figure 8
Noise N2
attenuation using adaptive VUSS solution. —
Spectrum with ANC, - - - spectrum
without ANC.Figure 9
Noise N5 attenuation using adaptive VUSS solution. — Spectrum
with ANC, - - - spectrum without ANC.Figure 10
Attenuation (in dB) of tonal sounds for adaptive filter tuning. —
RLS algorithm, - - -  FXLMS
algorithm.Table1 summarises the results obtained during
simulations for broadband signals N1–N6. The simulation experiments proved
adaptive VUSS filters tuning superiority over analytical suboptimal filter
design.


### 4.4. Real Experiments

To perform the
real experiments using analytical filter design it was necessary to identify
the models of electroacoustic pathsS11(z)−S22(z)  (see Figure 3) first. The identification of FIR
models of these paths was performed using LMS algorithm prior to each
experiment. Moreover, in case of both the analytical and adaptive filter
designs it was necessary to identify the model of the secondary path before the
ANC. The identification of FIR model of this path was performed after VUSS was
tuned by online identification procedure. The procedure, however, was disabled
during the ANC.The first ANC laboratory experiments used tonal
sounds. The frequency range between 100 and 950 Hz has been checked with
resolution of 20 Hz. The VUSS filter design was performed using adaptive
method. The attenuation obtained during these experiments has been shown on
Figure10. The figure shows that in frequency range from 150 to 950 Hz the
attenuation was between 14 and 33 dB (average 25 dB) for FXLMS ANC algorithm
and between 15 and 44 dB (average 30 dB) for RLS ANC algorithm.The drop in attenuation above a frequency of about 420 Hz, especially distinct in case of RLS algorithm, can be explained as the ducts
first cut-on frequency appears at 425 Hz. The figure shows that the performances
of both FXLMS and RLS algorithms above first cut-on frequency are similar.
However, below a frequency of 400 Hz the performance of RLS algorithm is up to
14 dB better than in case of FXLMS algorithm.The next step was to compare the performance of ANC
system in case of analytical and adaptive filter design. During these
experiments frequency range from 100 to 900 Hz was tested with resolution of
100 Hz. The results are shown on Figure11. In all cases the efficiency of ANC
algorithm with analytical VUSS tuning was worse than efficiency of the same
algorithm with adaptive VUSS tuning. Therefore, only adaptive VUSS design was
considered for the following tests.Figure 11
Comparison of attenuation
of tonal sounds for analytical and adaptive VUSS filter tuning. —
FXLMS algorithm with analytical tuning, - - -  RLS algorithm with analytical tuning, ...
RLS algorithm with adaptive tuning.Figure 12
Noise N1
(broadband noise 200–300 Hz) attenuation. —
Spectrum with ANC, - - - spectrum
without ANC.Figure 13
Noise N2 (broadband noise
700–800 Hz) attenuation. —  Spectrum with
ANC, - - -  spectrum without ANC.Figure 14
Noise N3 (food
processor) attenuation. — Spectrum with
ANC, - - - spectrum without ANC.Figure 15
Noise N4 (power transformer)
attenuation. —  Spectrum with ANC, - - -  spectrum without ANC.Figure 16
Noise N5 (power
plant turbine) attenuation. — Spectrum
with ANC, - - - spectrum without ANC.Figure 17
Noise N6 (electric propeller fan)
attenuation. — Spectrum with ANC, - - - spectrum without ANC.Finally, the noise signals described in Section4.2
were used in real experiments. Although both FXLMS and RLS ANC algorithms were
checked for performance, the only spectra presented on Figure 12
through Figure 17 are
those obtained with RLS algorithm, as this algorithm performance was superior to
FXLMS in all experiments. The attenuation factors for both FXLMS and RLS
algorithms are presented in Table 2.Table 2
The
attenuation in real experiments.
SignalAdaptive solution + FXLMSAdaptive solution + RLSN112 dB14 dBN27 dB7 dBN34 dB8 dBN413 dB16 dBN53 dB8 dBN63 dB5 dB


## 4.1. Active Noise Control Algorithms

Feedforward
active noise control system using virtual unidirectional sound source has been
shown in Figure4. This system uses a finite impulse response filter as control
filter with different adaptation algorithms. As in case of all feedforward
algorithms, only the reference value x(i) is used by the control filter F(z) to produce the set point value t(i).
The set point value is then processed by VUSS filters W1(z) and W2(z) to give two control values u1(i) and u2(i) that are amplified and sent to the
loudspeakers (via amplifiers and reconstruction filters not shown on the
figure). The adaptation algorithm on the other hand uses the reference signal x(i) as well as the error signal e(i) to tune the control filter coefficients.Figure 4
Active noise control system using virtual unidirectional sound source.The duct used in experiments described below was made
out of wood. It was 4 m long, with0.2×0.4 m rectangular section. One of the
duct ends was terminated with noise generating loudspeaker, while the other was
opened. The attenuating loudspeakers were of 0.16 m diameter and were situated
approximately in the middle of the duct. The distance between the middles of
the loudspeakers was equal to 0.3 m. The reference microphone was located 1.23 m from the middle of L2a loudspeaker, while the error microphone was separated
from L2b loudspeaker by the distance of 0.23 m.It should be emphasised that the goal of active noise
control algorithm is in contradiction to the goal of the virtual unidirectional
sound source adaptation algorithm [11]. The former tries to assure the compensation of sound
waves from both the primary and secondary sources at the microphone M2 point
without taking any notice of what happens at the microphone M1 point. The goal
of the latter is to compensate sound waves from the secondary sources only at
the microphone M1 point and try to assure the set point value appearing without
any alternation at the microphone M2 point. This leads to a conclusion that
both the algorithms should never operate at the same time. Indeed, the
experiments show that when both the adaptation algorithms are in operation the
whole system goes unstable. Therefore, the VUSS adaptation algorithm is usually
run on the beginning of experiments and switched off after adaptation is
completed. Next, only ANC algorithm is in use. Only when significant changes in
the environment are detected (e.g., temperature change above some level), the ANC
adaptation is temporarily frozen and the VUSS is retuned.FXLMS Algorithm
The first of
ANC filter adaptation algorithms tested was Filtered-x LMS algorithm [4, 9]. Assume that  the F(z) is an FIR filter with coefficients (see Figure 4)
in the  ith sample of time stored in a form of a vector f(i).
In that case, the FXLMS algorithm update equation is given by(48)f(i+1)=f(i)+μ⋅x′(i)⋅e(i),where μ is the step size and x′(i) is the vector of reference signal x(i) samples filtered by the secondary path
transfer function estimate, that is, an estimate of the transfer function
between set point value signal t(i) and error signal e(i).
If VUSS algorithm was working perfectly, the secondary
path transfer function would be equal to simple time delay. Unfortunately, the
study from  Section3.1 leads to a conclusion that such perfect situation is
impossible: the VUSS algorithm has to “pay more attention” to neutralising
acoustic feedback and secondary path equalisation is only its additional task.RLS Algorithm
The
second algorithm  of
ANC filter adaptation algorithms tested was recursive least squares (RLS)
algorithm (see, e.g., [9]). RLS algorithm uses very similar update equation in
the form(49)f(i+1)=f(i)+k(i)⋅e(i),where k(i) is the gain vector showing how much the value
of e(i) will modify different filter coefficients. The
gain vector is calculated as(50)k(i)=P(i)⋅x(i),where P(i) is a matrix updated in each step according to
the following equation:(51)P(i)=1λ(P(i−1)−P(i−1)x(i−1)xT(i−1)P(i−1)λ+xT(i−1)P(i−1)x(i−1)).
The RLS algorithm usually converges faster than the
LMS algorithm [12].
Because RLS algorithm needs to updateP matrix of size equal to the number of filter
coefficients in each adaptation step, it requires computational effort of N2,
whereas LMS algorithm requires computational effort of N only [13]. There is, however, a family of RLS algorithm
implementations called fast RLS algorithms that allows to omit computation of P matrix and use a selection of row vector
instead [13, 14]. This work used one of such
fast algorithms, called fast
transversal filter [15]. The algorithm was parametrised as follows. The
initial value of the minimum sum of backward a posteriori prediction-error
squares was equal to 1 and the exponential weighting factor was equal to
0.9999. The algorithm showed no stability issues, as it was expected on
floating-point arithmetic platform.
Both the FXLMS and RLS algorithms were implemented in
C programming language. Testing and debugging were performed using developed
simulation platform for Linux operating system. The platform allowed to emulate
DSP processor board behaviour and therefore the next step, moving the code
onto Texas Instruments TMS320C31 processor board, was purely automatic.
Another PC-computer program was used to tune various algorithm parameters and
to acquire the results. The sampling frequency of 2 kHz was chosen as frequency
band up to 1 kHz was of the author interest. It appeared that, due to hardware
limitations, appliable filter lengths were up to 70 with this sampling
frequency.


## 4.2. Testing Sounds

The set of signals chosen for experiments is presented
with dashed line on
Figures12–17. Testing signals N1 and N2
were generated offline as white noise filtered with high-order bandpass
filter. Bandpass was 200–300 Hz in case of N1 signal and 700–800 Hz in case of
signal N2. The former was chosen to show low-frequency attenuation capabilities
and the latter to show high frequency attenuation. The N3 is a signal acquired
in close vicinity of a food processor. It has the dominant frequency of about
270 Hz, but with substantial amount of broadband noise. The N4 is a signal
collected near a power transformer. It has all the harmonics of 50 Hz visible
with 250 Hz being the dominant. The N5 signal was recorded in water power plant
turbine proximity. It has the harmonics of 100 Hz distinguishable among
broadband noise. The last signal N6 was recorded in small bureau with an
electric propeller fan turned on. It is similar to N5, but the dominant
frequency is higher.All the signals were filed using 1000 samples. This
ensemble was repeated many times and (after amplification) sent to loudspeaker
L2 (see Figure4). The results of attenuation were measured with microphone M2
and calculated as(52)𝒯=10log(MSVofe(i)withoutattenuationMSVofe(i)withattenuation)[dB].


## 4.3. Simulation Results

The simulations
were performed using 150th order FIR filters models of the duct paths,
identified offline using methodology presented in [16]. The algorithm for
suboptimal VUSS filter solution described in Section 3.1 was tested first. As
expected, it allowed for very effective acoustical feedback cancellation but at
the cost of the secondary path being substantially degenerated. The amplitude
spectrum of transfer function between set point value t(i) and both the reference signal x(i) and error signal e(i) is presented on Figure 5.Figure 5
Amplitude spectrum of
transfer function between set point valuet(i) and x(i) (solid) and e(i) (dotted).After tuning phase, the system was engaged to active
noise control algorithms: FXLMS and Fast RLS. The ANC filter having 70
coefficients proved to be long enough to cancel out the noise in case of tonal
sounds, so this was the chosen value. However, in case of N1–N6 signals
attenuation was only (as for simulations of feedforward controller) 6–16 dB
(with FXLMS algorithm, which proved to be better), see
Table1. The examples
of spectra of the signals observed at the error microphone before and after
attenuation have been presented on Figures 6 and
7.Table 1
The
attenuation in simulations.
SignalAnalytical solution + FXLMSAnalytical solution + FXLMSAnalytical solution + RLSN19 dB15 dB16 dBN216 dB30 dB33 dBN310 dB15 dB16 dBN415 dB22 dB22 dBN59 dB17 dB18 dBN66 dB9 dB8 dBFigure 6
Noise N2
attenuation using analytical VUSS solution. —
Spectrum with ANC, - - -  spectrum
without ANC.Figure 7
Noise N5 attenuation using analytical VUSS
solution. — Spectrum with ANC, - - -  spectrum without ANC.In the following experiments, adaptive algorithm described
in Section3.2 was responsible for VUSS filters tuning. The value of the β parameter (see (45)) has been chosen as
0.5 to put the same effort into satisfying both the goals defined by (12).
The attenuation of tonal sounds was also infinite in case of adaptive VUSS
filters tuning. But the attenuation of testing signals was slightly better:
9–30 dB in case of FXLMS algorithm and 8–33 dB in case of RLS algorithm, see
Table 1. The examples of spectra of the signals observed at the error
microphone before and after attenuation have been shown on
Figures 8 and  9. The
difference in amplitude spectrum of the N2 and N5 signals without attenuation
visible on these figures and on Figure 6 and
7 comes from the fact that this is
the spectrum measured with the M2 microphone, not the spectrum of the signal
driving the L1 loudspeaker.Figure 8
Noise N2
attenuation using adaptive VUSS solution. —
Spectrum with ANC, - - - spectrum
without ANC.Figure 9
Noise N5 attenuation using adaptive VUSS solution. — Spectrum
with ANC, - - - spectrum without ANC.Figure 10
Attenuation (in dB) of tonal sounds for adaptive filter tuning. —
RLS algorithm, - - -  FXLMS
algorithm.Table1 summarises the results obtained during
simulations for broadband signals N1–N6. The simulation experiments proved
adaptive VUSS filters tuning superiority over analytical suboptimal filter
design.


## 4.4. Real Experiments

To perform the
real experiments using analytical filter design it was necessary to identify
the models of electroacoustic pathsS11(z)−S22(z)  (see Figure 3) first. The identification of FIR
models of these paths was performed using LMS algorithm prior to each
experiment. Moreover, in case of both the analytical and adaptive filter
designs it was necessary to identify the model of the secondary path before the
ANC. The identification of FIR model of this path was performed after VUSS was
tuned by online identification procedure. The procedure, however, was disabled
during the ANC.The first ANC laboratory experiments used tonal
sounds. The frequency range between 100 and 950 Hz has been checked with
resolution of 20 Hz. The VUSS filter design was performed using adaptive
method. The attenuation obtained during these experiments has been shown on
Figure10. The figure shows that in frequency range from 150 to 950 Hz the
attenuation was between 14 and 33 dB (average 25 dB) for FXLMS ANC algorithm
and between 15 and 44 dB (average 30 dB) for RLS ANC algorithm.The drop in attenuation above a frequency of about 420 Hz, especially distinct in case of RLS algorithm, can be explained as the ducts
first cut-on frequency appears at 425 Hz. The figure shows that the performances
of both FXLMS and RLS algorithms above first cut-on frequency are similar.
However, below a frequency of 400 Hz the performance of RLS algorithm is up to
14 dB better than in case of FXLMS algorithm.The next step was to compare the performance of ANC
system in case of analytical and adaptive filter design. During these
experiments frequency range from 100 to 900 Hz was tested with resolution of
100 Hz. The results are shown on Figure11. In all cases the efficiency of ANC
algorithm with analytical VUSS tuning was worse than efficiency of the same
algorithm with adaptive VUSS tuning. Therefore, only adaptive VUSS design was
considered for the following tests.Figure 11
Comparison of attenuation
of tonal sounds for analytical and adaptive VUSS filter tuning. —
FXLMS algorithm with analytical tuning, - - -  RLS algorithm with analytical tuning, ...
RLS algorithm with adaptive tuning.Figure 12
Noise N1
(broadband noise 200–300 Hz) attenuation. —
Spectrum with ANC, - - - spectrum
without ANC.Figure 13
Noise N2 (broadband noise
700–800 Hz) attenuation. —  Spectrum with
ANC, - - -  spectrum without ANC.Figure 14
Noise N3 (food
processor) attenuation. — Spectrum with
ANC, - - - spectrum without ANC.Figure 15
Noise N4 (power transformer)
attenuation. —  Spectrum with ANC, - - -  spectrum without ANC.Figure 16
Noise N5 (power
plant turbine) attenuation. — Spectrum
with ANC, - - - spectrum without ANC.Figure 17
Noise N6 (electric propeller fan)
attenuation. — Spectrum with ANC, - - - spectrum without ANC.Finally, the noise signals described in Section4.2
were used in real experiments. Although both FXLMS and RLS ANC algorithms were
checked for performance, the only spectra presented on Figure 12
through Figure 17 are
those obtained with RLS algorithm, as this algorithm performance was superior to
FXLMS in all experiments. The attenuation factors for both FXLMS and RLS
algorithms are presented in Table 2.Table 2
The
attenuation in real experiments.
SignalAdaptive solution + FXLMSAdaptive solution + RLSN112 dB14 dBN27 dB7 dBN34 dB8 dBN413 dB16 dBN53 dB8 dBN63 dB5 dB


## 5. Conclusions

The idea of
virtual unidirectional sound source presented in this paper is based on
theoretical study of wave propagation in duct. VUSS itself is a special case of
two-reference, two-output system with a detailed analysis presented in Section3.
Its application to active noise control in an acoustic duct proved to be
effective, resulting in 20–40 dB attenuation of tonal sounds and 5–16 dB
attenuation of complex signals with broadband noise. The best results were
obtained with adaptive VUSS design together with RLS active noise control
algorithm. These results are comparable with similar results reported by other
authors. The results, however, do not show one of the strengths of VUSS, the
ability to adapt to time varying feedback paths.



---

*Source: 647318-2008-06-22.xml*