Analysis and Simulation of Frost's Beamformer

RADIOENGINEERING, VOL. 12, NO. 2, JUNE 2003 1 Analysis and Simulation of Frost’s Beamformer Miroslav ˇ STRUPL, Pavel SOVKA Dept. of Circuit Theory, Czech Technical University, Technick´ a 2, 166 27 Praha 6, Czech Republic [email protected], [email protected] Abstract. Sensor arrays are often used for a signal separation from noises using the information about the direction of arrival. The aim of this paper is to analyze Frost’s beamformer with respect to the speech preprocessing for the hearing impaired people. The frequency response of the system including the background noise attenuation are derived as functions of the direction of arrival. The derivation supposes a uniform linear array of sensors and plane waves. It is shown that the number of possible configurations can be decreased by using some symmetries. The impact of the used algorithm constraint on the frequency response and subsequently on the directional noise suppression is derived analytically. Keywords Adaptive arrays, Frost’s beamformer, speech enhancement, noise reduction. 1. Introduction This article deals with the application of Frost’s beamformer for the background noise reduction. Frost’s beamformer belongs to the broad field of beamforming systems, see. e.g. [7], [8], [9]. The application of this beamforming system discussed in this paper is the speech preprocessing for hearing impaired people. Other possible applications of this system can be speech processing for conferencing and a mobile telephony (direct sequence CDMA) [1], where an array of sensors helps to discriminate between signal and noises sources. For example, a hearing impaired person using a small microphone array mounted on glasses is able to discriminate between different speakers. In these applications the broadband rather than narrowband input signals should be assumed. Thus the frequency response and SNRE (Signal to Noise Ratio Enhancement) evaluated over the whole frequency band have to be used. The following restrictions about an environment, a geometry of array and input signals are made. •In order to obtain the frequency response easily the environment should be homogeneous, isotropic, lossless, nondispersive with respect to the wave propagation (phase) velocity. •The array of sensors are assumed to be a linear with uniformly spaced sensors. Sensors are assumed to have omnidirectional characteristics. •Both the directional signals and noises are assumed to propagate with the plane wave front. This assumption is not too restrictive because spherical waves generated in a far field from the array can be approximated with plane waves. •Noises are supposed to be either directional (forming the plane wave - see preceding item) or uncorrelated between sensors. •Frequency response is derived under assumption of linearity and time invariance - LTI (no coefficients update) of the filter. The paper is organized as follows. In section 2 Frost’s beamformer is described including a linear constraint. In section 3 its frequency response for a plane wave input is derived. In section 4 two symmetries of configurations are proved. In section 5 properties of the frequency response caused by the constraint are discussed. In section 6 the influence of a noise uncorrelated between sensors and samples is shown. Finally, the conclusion is given. 2. Frost’s Beamformer Frost’s beamformer [2] (see Fig. 1a) consists of an array with Ksensors, where each sensor is followed by a transversal filter with Jweights. The number of weights is equal for all transversal filters. The sum of the filter outputs is the beamformer output. Weights are updated by Frost’s constrained least mean square (CLMS) algorithm which minimizes the mean square error of the output signal while satisfying a constraint [2]. In order the input signal s(t)(arriving to all sensors with the same delay) to be passed without any distortion, the impulse response of the whole system must be equal to the unit impulse. This impulse response represents the constraint for the weights of all filters. The whole system can be replaced by one transversal FIR filter for the signal s(t). The replacement is shown in Fig. 1b, where f1, f2, . . . , fJis the impulse response for the signal. Constraint equations (see Fig. 1b) can be written also in 2 M. ˇ STRUPL, P. SOVKA, ANALYSIS AND SIMULATION OF FROST’S BEAMFORMER matrix form W     1 1 . . . 1      =     f1 f2 . . . fJ      ,(1) where Wstands for weight matrix with real elements W=     w1w2. . . wK wK+1 wK+2 . . . w2K . . .. . ..... . . wJK−K+1 wJK−K+2 . . . wJK      .(2) To discuss the Frost’s beamformer behavior in details, let us define some terms needed. The digitized input noisy signals xi[n], i = 1,2, . . . , JK are formed by components of both clean signal s(t)and noise n(t). The vector ˜ x[n] represents noisy signals on taps, the vector wconsists of weights value, the vector Frepresents the constrained impulse response and the matrix Cwill be used in constraint formulation ˜ xT[n] = [ x1[n]x2[n]. . . xJK [n]], wT= [ w1w2. . . wJK ], FT= [ f1f2. . . fJ], C= [ c1c2. . . cJ].(3) Elements ciare the column vectors of length JK with (i−1)Kzeros followed by Kones and (J−i)Kzeros cT i= [ 0 0 . . . 0 | {z } (i−1)Kzeros 1 1 . . . 1 | {z } Kones 0 0 . . . 0 | {z } (J−i)Kzeros ]. (4) Now the problem of finding the optimum weight vector for a stationary signal wopt (Wiener solution) can be formulated. The weight vector minimizing E[y2[n]] = wTE[˜ x[n]˜ x[n]T]w=wTRxxwand satisfying the constraint CTw=Fhave to be found. Rxx stands for the autocorrelation matrix. In [2] the method of Lagrange multipliers was used to obtain the Wiener solution wopt =R−1 xx C(CTR−1 xx C)−1F(5) and the adaptive CLMS algorithm w[0] = f, w[n+ 1] = P(w[n]−µy[n]˜ x[n]) + f.(6) The vector fand the projection matrix Pare defined as f=C(CTC)−1F, P=E−C(CTC)−1CT.(7) Positive scalar µis a step-size parameter. The choice of µ is the tradeof between the convergence time and the missadjustment of weights from Wiener solution. An easily computable upper bound for µis given by µ < 2/(3E[˜ xT˜ x]). a xK+1[n] xK+2[n] x2K[n]xK[n] x2[n] x1[n]x1(t) x2(t) xK(t) s(t) xJK−K+1[n] xJK−K+2[n] xJK [n] y[n] wK w2 w1 wK+2 wK+1 wJK−K+2 wJK−K+1 w2KwJK n(t) b s[n]s(t) y[n] fJ=PJK i=JK−K+1 wi f1=PK i=1 wif2=P2K i=K+1 wi Fig. 1. a Frost’s beamformer structure, bFrost’s beamformer from s(t)view - constraint formulation The convergence behavior and the choice of µis deeply discussed in [2]. The alternative form of equation (6) for the implementation is wi[n+ 1] = wi[n]−µy[n]xi[n] −1 K (bi Kc+1)K X j=bi KcK+1 (wj[n]−µy[n]xj[n]) + fbi Kc+1 K.(8) 3. Frequency Response Analysis Now the frequency response for the LTI case (fixed weights, fixed geometry, ...) and plane wave fronts will be derived. Since there are digital filters behind sensors, it will be more convenient for the analysis to replace the transfer paths from the source to sensors by a discrete time model. Thus the variable twill be replaced by the nTs, where Ts is the sampling period used for filters behind sensors. The main reason is to understand the constraint influence on the system behavior. Now the coordinate system used for the sources description will be given (see Fig. 2). Note that the origin is placed in the center of the array and y-axis goes through all the sensors of the array, which are spaced uniformly with the distance dand are numbered from 0to K−1. The sensors are assumed to have omnidirectional characteristics. Due to the rotation symmetry of the array (about y-axis) every threedimensional scenario with arbitrary plane wave sources can be described in x∈<0,+∞),y∈(−∞,+∞)half-plane (signals at sensors will be the same). Thus without the loose of generality only wave fronts generated by sources lying in this half-plane and having normal vectors in xy-plane can be considered. Let us assume one source generating a plane wave with a waveform u(t). Because the normal vector of RADIOENGINEERING, VOL. 12, NO. 2, JUNE 2003 3 a wave front at the source lies in xy-plane, the whole wave front is projected onto a single line (depicted by the thick line labeled by u(t)). Thus for the source description it is sufficient to give the position rof the nearest point of this wave front to the origin (note that the vector ris perpendicular to the wave front). It seems to be the most appropriate to express rin polar coordinates by giving angle ϕ∈<−π 2,π 2> and distance r∈ R+. The distance rcan be measured using delay between the source and the center of the array R=r Tsv. The direction of arrival can be also measured using delay C(in samples) between adjacent sensors. C relates to ϕaccording to C=dsin ϕ Tsv,(9) where ddenotes the distance between adjacent sensors, vthe speed of sound. Due to the assumed environment (homogeneous, isotropic, lossless, nondispersive) u(t)reaches the sensor mwith the delay Ts[R+C(m−K−1 2)] without any distortion. Thus the frequency response (measured from the source to a sensor m) is ˆ Pm(ω) = e−jωTs[R+C(m−K−1 2)], where ωstands for the angular frequency. To obtain corresponding frequency response Pmof the discrete time model ˆ Pmhas to be restricted to the range (−ωs 2,ωs 2)(antialiasing filters must be placed behind sensors). Using the normalized frequency Θ = ωTsfor ω∈(−ωs 2,ωs 2)and the periodic extension of the frequency interval we obtain ∀k∈ Z,∀Θ∈<−π, π >: Pm(Θ + 2kπ) = e−jΘ[R+C(m−K−1 2)].(10) It is convenient to use the normalized frequency when describing a digital system. This approach assures the independence of the system description on a sample period Ts (or on a sample frequency fs=1 Ts). The frequency response for the whole system (see Fig. 1a) can be derived using frequency responses of transversal filters and (10) ∀k∈ Z,∀Θ∈<−π, π >:H(Θ + 2kπ) = = K−1 X m=0 e−jΘ[R+C(m−K−1 2)]wm+wK+me−jΘ+. . . +wK(J−1)+me−j(J−1)Θ= =e−jΘRejΘC(K−1 2) ·1e−jΘ. . . e−j(J−1)Θ W     1 e−jCΘ . . . e−j(K−1)CΘ      , (11) where the matrix Wis defined in (2). Now the power gain AP will be derived for a system with the frequency response given by (11). The input waveform u[n]of the system is assumed to be the white stationary process with zero mean and variance σu. Let us define AP as the ratio of expectations of sensor K 2 s.0 0 ϕ x sensor K−1 r d ∼R y ∼ |C| u(t) Fig. 2. Configuration used to derive frequency response (depicted for even number of sensors) squared output and input AP =E[v2[n]] E[u2[n]] = =1 σ2 u ∞ X k=−∞ ∞ X l=−∞ h[k]h[l]E[u[n−k]u[n−l]] = = ∞ X k=−∞ h2[k] = 1 2π π Z −π |H(Θ)|2dΘ, (12) where h[n]and v[n]are the impulse response of the system and the system output respectively (Parceval’s theorem was used in the last step). The H(Θ) should be written in more suitable form than (11) to proceed with the integration ∀Θ∈<−π, π >:H(Θ) = JK X i=1 wie−jl(i)Θ.(13) The l(i)is a new function expressing the delay connected with the weight wi1 l(i) = C(i−1) −Kbi−1 Kc+bi−1 Kc.(14) We can continue with the integration of equation (12) AP =1 4πZπ −πH(Θ)H(Θ)∗+ (H(Θ)H(Θ)∗)∗= =1 2π JK X r=1 JK X s=1 wrwsZπ −π cos([l(r)−l(s)]Θ)dΘ = = JK X r=1 JK X s=1 wrws sin([l(r)−l(s)]π) [l(r)−l(s)]π. (15) 4. Configurations The verification of the system behavior by the simulation cannot be carried out for all the possible configurations. 1An arbitrary sensor placement could be expressed by changing l(i)definition 4 M. ˇ STRUPL, P. SOVKA, ANALYSIS AND SIMULATION OF FROST’S BEAMFORMER Thus the number of configurations have to be decreased. This can be done by using some types of equivalences. The most convenient choice seems to be the equivalence with regard to generated output, because it seems to be invaluable to simulate set of configurations giving the same output. Of course, the precise definition depends on what we are interested in. Such the equivalence is moreover conditioned by the equality of input signals usually. This means, the system is considered to be a black box with inputs and one output. Needed definitions will be introduced now. Definition 1 The configuration is an ordered triplet of input signals used, transfer functions from sources to sensors (a function of an environment, a placement of sources and sensors and generated wave fronts), the filter behind sensors (the filter is described by the filter order and weight values, filtration and weight update equations). Definition 2 Configurations Kand K0are equivalent K≡ K0if and only if Kand K0having the same input signals generate the same output ∀n∈ Z+:y[n] = y0[n]. 4.1 Equivalent Configurations with Different Distance between Sensors The frequency response from the plane wave source to the chosen sensor (10) depends only on the delay Cand neither on the angle ϕnor on the distance between adjacent sensors d. Consider the configuration Kwith the adjacent sensor distance d,splane wave sources with directions of arrival given by angles ϕ1, ϕ2, . . . ϕsleaving delays C1, C2, . . . Cs. Let us construct the configuration K0from Kby changing dto d0and choosing ϕ0 1, ϕ0 2, . . . ϕ0 sunder the condition C0 1=C1, C0 2=C2, . . . C0 s=Cs. Thus the frequency response (10) will be the same. Therefore the filter behind sensors gets the same input and generates the same output and consequently K≡K0according Def. 2. But how do we chose ϕ0 iin order to Ciremain the same? Equation (9) must hold for such ϕi:Ci=d0sin ϕ0 i Tsv. Since ϕiis restricted to <−π 2,π 2>, preceding relation between Ciand ϕ0 ican be rewritten in the form ϕ0 i= arcsinvCiTs d0i= 1,2, . . . s. (16) 4.2 Symmetry of Filters behind Sensors Consider two configurations K, K0which are identical except for directions of arrival of plane waves ϕ0 1=−ϕ1, ϕ0 2=−ϕ2,. . . ϕ0 s=−ϕs. Thus from (9) C0 1=−C1, C0 2=−C2,. . . C0 s=−Cs. Then equation (10) must hold for spectra of waveforms on sensors ∀k∈ Z,∀Θ∈<−π, π >:XK−1−m(Θ + 2kπ) = = s X i=1 e−jΘCi(K−1−m)ejΘCiK−1 2e−jΘRiUi(Θ) = s X i=1 e−jΘ(−Ci)mejΘ(−Ci)K−1 2e−jΘRiUi(Θ) =X0 m(Θ + 2kπ). (17) It can be concluded, that waveforms on sensors are the same except for the permutation x0 0[n] = xK−1[n], x0 1[n] = xK−2[n], . . . x0 K−1[n] = x0[n]. (18) If a noise uncorrelated between sensors and samples is present, its components on sensors should satisfy (18) in order to waveforms on sensors satisfy (18). If there is the pure summation behind sensors (delay and sum beamformer), the permutation takes no effect due to the commutativity of addition. Thus generated outputs are the same in both configurations and K≡K0from Def. 2. The case of the filter driven by Frost’s CLMS algorithm will be discussed now. It can be stated after looking at the filter structure Fig. 1aand considering the property (18), that ˜ x0[n]is the permutated version of ˜ x[n](at most after J−1iterations). The permutation occurs within each column. The permutation vector πand projection πperforming that permutation are defined for easier manipulation π= [π1, . . . πJK ]T= [K, K−1, . . . 1,2K, 2K−1, . . . K+1 . . . JK, JK −1, . . . JK −K+ 1]T, π([x1x2. . . xJK ]T) = [xπ1xπ2. . . xπJK ]T.(19) Let us assume that Kand K0initial conditions fulfill the described property ˜ x0[0] = π(˜ x[0]),w0[0] = π(w[0]).(20) The usual initial conditions (zero xi[0] values and weight setup f) satisfy (20), because these initial conditions have the same values within each column. The observation made about ˜ x0[n]and ˜ x[n]components now holds for every iteration ∀n∈ Z+:˜ x0[n] = π(˜ x[n]).(21) In order to show K≡K0, the following conditions remain to be proved ∀n∈ Z+:y0[n] = y[n],w0[n] = π(w[n]).(22) The induction scheme will be followed here. The desired equation for the weight vector is included in (20) for n= 0. Using (20) following equation can be written for outputs y0[0] = w0[0]T˜ x0[0] = π(w[0]T)π(˜ x[0]) = y[0].(23) The last equality follows from the fact that if components of both vectors are equally permutated, then scalar multiplication leaves the same result due to the commutativity of addition. Now suppose that the property holds for some n≥0 y0[n] = y[n], w0 i[n] = wπi[n].(24) RADIOENGINEERING, VOL. 12, NO. 2, JUNE 2003 5 Notice that any column number (from 0to J−1) for a weight wicould be expressed by bi Kc. The updating equation (8) for a weight w0 ican be written w0 i[n+ 1] = w0 i[n]−µy0[n]x0 i[n] −1 K (bi Kc+1)K X j=bi KcK+1 (w0 j[n]−µy0[n]x0 j[n]) + fbi Kc+1 K.(25) Using the substitution from (24) and (21) into (25) leads to (notice that bi Kc=bπi Kc) w0 i[n+ 1] = wπi[n]−µy[n]xπi[n] −1 K (bπi Kc+1)K X j=bπi KcK+1 (wπj[n]−µy[n]xπj[n]) + fbπi Kc+1 K. (26) Because the permutation occurs only within one column and the summation goes through the whole column, summing non permutated version leaves the same result w0 i[n+ 1] = wπi[n]−µy[n]xπi[n] −1 K (bπi Kc+1)K X j=bπi KcK+1 (wj[n]−µy[n]xj[n]) + fbπi Kc+1 K =wπi[n+ 1]. (27) The last equality was obtained using (8). This can be rewritten using vector notation w0[n+ 1] = π(w[n+ 1]).(28) Using (21),(28) and following the same reasoning as in obtaining (23) the equation for y0[n+ 1] can be written y0[n+ 1] = w0[n+ 1]T˜ x0[n+ 1] = =π(w[n+ 1]T)π(˜ x[n+ 1]) = y[n+ 1].(29) This completes the induction step. Equations (22) are proved and through Def. 2 K≡K0. The consequences of introduced properties are discussed now for the uniform linear array, plan waves, the filter behind sensors driven by CLMS and the signal coming from ϕ= 0◦. If one directional noise source is present (noise uncorrelated between sensors and samples may be present too), then it is sufficient to simulate configurations with noise directions ϕfrom 0◦to 90◦, because other configurations are equivalent with them due to 4.2 (see sect. Symmetry of Filters behind Sensors). Obtained results can be further generalized to different distances between sensors according to 4.1 (see sect. Equivalent Configurations with Different Distance between Sensors). If there are an arbitrary number of noise sources, then it is sufficient to simulate system for directions ϕbetween −90◦and 90◦. Obtained results can be generalized by using 4.1 similarly as in previous case. Non zero distance r(proportional to delay R) causes only superposition of the linear trend to the phase of frequency response (see (11)). Thus it seems sometimes to be uninteresting to consider different r. 5. Constraint Influence Now some properties of the frequency response for plane wave inputs (11) will be discussed. If the elements of the rightmost vector (11) became ones for some Θ, then the resulting expression simplifies greatly. But for which normalized frequency Θthis simplification can be made? This is satisfied for Θ = 2π Cl l ∈ Z,|l| ≤ bC 2c.(30) Now substitution from (1) to (11) for Θgiven by (30) can be made H(2π Cl+ 2kπ) = =ej2π ClC(K−1 2)e−j2π ClR 1e−j2π Cl. . . e−j2π Cl(J−1)  ·f1f2. . . fJT =ej2π Cl(C(K−1 2)−R)J−1 X i=0 fi+1e−j2π Cli l∈ Z,|l| ≤ bC 2c. (31) If fiis the unit impulse, then this expression further simplifies to H(2π Cl+ 2kπ) = ej2π Cl(C(K−1 2)−R)l∈ Z,|l| ≤ bC 2c. (32) It can be concluded that due to the imposed constraint the frequency response is fixed at some normalized frequencies. These normalized frequencies are determined by the direction C(30). Fig. 3 illustrates this for some directions Cwhen the unit impulse constraint is assumed (black dots stands for fixed points). Using (30) it can be concluded about normalized frequencies fixed by the constraint (assuming one half period of frequency response): 1) The zero normalized frequency is always included. 2) Fixed normalized frequencies are uniformly spaced with spaces 2π C. Thus as Cincreases the spaces between fixed normalized frequencies decreases. 3) The number of fixed normalized frequencies increases with increasing Ceven. The values of frequency response at these normalized frequencies are determined by (31). Following can be noted: 4) The value of frequency response at zero normalized frequency (l= 0) is the same as the value of frequency response at zero normalized frequency for the signal s(t)coming through the path with impulse response fi. 5) The unit impulse constraint causes all magnitudes of fixed normalized frequencies to be one (32). If the signal should be passed without any distortion then the unit impulse constraint is the most convenient choice. Assume that the plane wave of the white stationary noise is arriving from the direction C. Then the output power 6 M. ˇ STRUPL, P. SOVKA, ANALYSIS AND SIMULATION OF FROST’S BEAMFORMER minimization could occur only at bands excluding fixed normalized frequencies. This leads to sharp local maxima at these normalized frequencies (see Fig. 3, dashed curves). Thus the power gain for the white stationary noise is almost determined by the number of these maxima. Therefore the dependence of the power gain on the direction Cis similar to the stair function where the rapid increasing occurs at C even (see item 3)). The existence of local maxima has also the impact on SNRE (Signal to Noise Ratio Enhancement) which therefore falls with every Ceven (see Fig. 4). It can be concluded that the delay range should be restricted to C between 0 and 2 samples. In the case of an input colour noise the situation is more complicated. The behavior of the system depends, on whether the noise spectrum hits or miss the normalized frequencies fixed by the constraint. In general, the delay range should be restricted similarly as in the previous case, because there is only the zero normalized frequency to hit (see Tab. 1 for examples). The car noise is the worst example. Most power of its spectra is spread around the zero normalized frequency, which is fixed by the unit impulse constraint to one (see 5)), thus most of the noise passes through the system. It motivates us to abandon the unit impulse constraint and turn rather to the high-pass constraint. However it excludes the zero normalized frequency from the signal spectrum and causes some distortion, it fixes the frequency response at the zero normalized frequency to zero value for directional noises according 4). If some distortion is present in the signal s(t)(due to the small length of transversal filters) it could be equalized by a subsequent filter. But this equalization must be done carefully, because this subsequent filtration affects the noise too. The proposed frequency response of Frost’s beamformer for the signal, the equalization filter and the resulting frequency response, respectively, are shown in Fig. 9. Results for this constraint are summarized also in Tab. 1. The choice of the frequency response for the signal assuming Griffiths-Jim beamformer [3] is discussed also in [4], leading to the same result. The examples in Tab. 1 are obtained by simulating following configuration (The same array and filter setup was used to obtain curves in Fig. 4). The array consists of three uniformly spaced sensors with the adjacent sensor distance d= 76mm ( the sound speed v= 330ms−1). The noise arrives from the direction ϕ= 30◦, which corresponds to delay C= 0.92. The adaptive filter driven by Frost’s CLMS algorithm used behind sensors has the convergence constant µ= 0.001,63 weights and the sampling frequency fs= 8000Hz. A speech (words: ”0,1,0,1”) is used as the input signal. Three types of input noise are used (see Tab. 1). Input waveforms are scaled to have the unit average power. Therefore the global SNR is 0dB. The steady state SSNRE is computed from last ”0,1” segments by time averaging from one realization. The length of used waveforms is 2.47 s. The case of the white stationary noise and the unit impulse constraint is further illustrated by showing the signal and the noise vs. time in Fig. 5a,SNRE vs. time in Fig. 5b, and spectrograms Fig. 6. The magnitude frequency response and the power gain for the white stationary noise vs. the direction of arrival are shown in Fig. 7 and Fig. 8. 1 0 Θ π C= 0.9|H(Θ)| 0 C= 2 π |H(Θ)| 1 Θ 2π 3 0 C= 3 π |H(Θ)| 1 Θ Fig. 3. Impact of constraint on frequency response - examples for different C, unit impulse constraint is assumed, black dots denotes fixed points CLMS initial setup Wiener solution Fig. 4. SNRE in dependence on direction Cfor unit impulse constraint and zero mean white stationary unit variance stochastic processes used for both signal and noise noise white st. fan car 0 Θ 0.1 π0 Θ 0.1 π0 Θ 0.1 π unit imp. |H|1 0 Θ |H| π 1 0 Θ |H| π 1 0 Θ |H| π s.s. SNRE 12.36 dB 12.10 dB 6.57 dB constraint type high-pass |H|1 0 Θ |H| π 1 0 Θ |H| π 1 0 Θ |H| π s.s. SNRE 15.18 dB 12.82 dB 22.22 dB Tab. 1 Reached frequency responses and s.s. SNRE for some directional noises arriving from direction (delay) C= 0.9[samples] and constraint type. RADIOENGINEERING, VOL. 12, NO. 2, JUNE 2003 7 a b Fig. 5. a Signals vs. time (desired speech - light gray, white stationary noise - dark gray), bSNRE vs. time (SNRE in each segment is shown) a b c Fig. 6. Spectrograms: adesired speech,bmixture on sensor, cfilter output Fig. 7. Magnitude frequency response on direction for weights sampled at time instant 1.75 s. Darker color corresponds to lower magnitude. Directions of desired signal and noise are labeled sig1 and sig2 respectively and marked with vertical lines. Fig. 8. Power gain on direction for weights sampled at time instant 1.75 s. Directions of desired signal and noise are labeled sig1 and sig2 respectively and marked with rays. Fig. 9. Frequency response for: desired signal (constraint) - dotted, equalization filter - dashed, result - solid 8 M. ˇ STRUPL, P. SOVKA, ANALYSIS AND SIMULATION OF FROST’S BEAMFORMER They are computed from weights using (11) and (15) sampled at the time instant 1.75 s (after the convergence is reached). The decreasing of the system gain close to the noise direction is evident in both cases. 6. Noise Uncorrelated between Sensors and Samples This section deals with the influence of the noise uncorrelated between sensors and the samples on the system. Noise moments are assumed to be the same across both samples (time stationarity) and sensors (spatial stationarity). The following trend can be observed. As the this type of noise dominates against the directional one, the weights approach to the initial ones. This could be clarified by finding the Wiener solution. The autocorrelation matrix takes the form Rxx =σ2 oE, where σ2 odenotes the variance of the this uncorrelated noise. Using (7) the Wiener solution (5) takes the following form wopt = (σ2 oE)−1CCTσ2 oE−1C−1F= =C(CTC)−1F=f.(33) Note that fstands for the initial weight vector. Thus it can be concluded that the adaptation process takes no effect and the filter remains in the initial state. The noise power at the output is then E[y2[n]] = wT optRxxwopt =fTσ2 oEf = σ2 oFT[(CTC)−1]TF=σ2 oFT[(KE)−1]TF=σ2 o K J X i=1 f2 i. (34) While the unit impulse constraint gives the power gain for the signal APs= 1, the power gain for the noise is APn=E[y2[n]] σ2 o =1 K(35) and SNRE SNRE =APs APn =K. (36) Thus for the unit impulse constraint, SNRE equals to the number of sensors used. It is the same as for the delay and sum beamformer. It is not so surprising, because the initial weights ffor the unit impulse constraint are non zero only in one column and their values are equal. In order to compare the SNRE (for the dominating uncorrelated noise between sensors and samples) with results summarized in Tab. 1 (for the directional noise), the equation (36) must be evaluated for the same number of sensors K= 3. This evaluation yields 4.8dB. 7. Estimating SNRE The estimation of SNRE computed directly from waveforms require many realizations to be sufficiently smooth. In the case of limited small number of signal realizations the batch or recursive estimates can be used with the assumption of ergodicity. But what one can do when only one realization is at hand (e.g. speech), and the use of mentioned estimates are rather devastating to give readable results? In this paragraph one alternative method is described which is derived under following assumptions: 1) The system under consideration is Frost’s beamformer driven by CLMS algorithm with the unit impulse constraint so the power gain for the signal is APs[n] = 1. 2) Input noises are mutually independent white stationary noises. 3) The current effective window of input noises (see below) is independent on current weights representing the instantaneous impulse responses. Of course non-LTI system has not any impulse response in general. But something like the impulse response can be assigned to each iteration. The current output sample is then computed by multiplying input waveform with current impulse response. Assumption 3) means that the effective window of the input waveform which goes into the multiplication is independent on the current impulse response. This is called the independence assumption and is used for example in [5] for LMS MSE derivation. Reasoning behind this is as follows. When the step size parameter µis very small then current weights are almost determined by past input samples and the current effective window has negligible impact. When the input samples are mutually independent then the current effective window can be thought to be independent of current weights (or current impulse response). Let us assume Ndirectional noises n11[n],...,n1N[n] at inputs and noise uncorrelated between sensors and samples n1o[n]. Let standard deviations of these noises are σ11, ...,σ1N,σ1o, impulse responses h1n[k], ...,hNn[k],hon[k](for iteration n) and power gains AP1[n], ...,APN[n],APo[n](see (15), (35) ) for iteration n. Let us define the sum of noises at the input n1[n]and the sum of the noises at the output n2[n]. n1[n] = N X i=1 n1i[n] + n1o[n](37) n2[n] = N X i=1 ∞ X ki=−∞ hin[ki]n1i[n−ki]+ ∞ X ko=−∞ hon[ko]n1o[n−ko](38) Let APn[n] = E[n2 2[n]] E[n2 1[n]] to be the power gain for noises in iteration n. And finally let us define the current SNRE [n] SNRE [n] = APs[n] APn[n]=E[n2 1[n]] E[n2 2[n]] =PN i=1 σ2 1i+σ2 1o E[n2 2[n]] . (39) RADIOENGINEERING, VOL. 12, NO. 2, JUNE 2003 9 During the last two steps in (39) the assumptions 1),2) and equation (37) were used. When assumptions 3),2) and equation (38) are used for the denominator of (39) it can be written E[n2 2[n]] = E( N X i=1 ∞ X ki=−∞ hin[ki]n1i[n−ki]+ ∞ X ko=−∞ hon[ko]n1o[n−ko])2= = N X i=1 E∞ X ki=−∞ h2 in[ki]σ2 i+E∞ X ko=−∞ h2 on[ko]σ2 1o = N X i=1 EAPi[n]σ2 i+EAPo[n]σ2 1o. (40) The substitution (40) to (39) gives SNRE [n] = PN i=1 σ2 1i+σ2 1o PN i=1 EAPi[n]σ2 i+EAPo[n]σ2 1o . (41) Because weights are smoothed by CLMS itself and second order moments of noises can be computed by time averaging due to their ergodicity, this estimate is much more smoother then previous mentioned batch and recursive estimates. The described approach was used to obtain time dependence of SNRE in Fig. 5. 8. Conclusion Uniform linear arrays and plane waves propagating in homogeneous, isotropic lossless nondispersive environment were assumed throughout this paper. Frequency response (11) and white stationary process power gain (15) were derived as functions of the direction of arrival. Some symmetries were shown using the configuration concept. It was found (4.2) that for the desired signal, one directional noise (noise uncorrelated between sensors and samples may be also present) and Frost’s CLMS algorithm behind sensors, it is sufficient to consider noise directions ϕbetween 0◦and 90◦. For greater number of directional noises it is sufficient to consider noise directions ϕbetween −90◦and 90◦. Obtained results can be further extended to different adjacent sensor distances using (4.1). The impact of constraints for the directional noise suppression was described. Discussion of this issue leads to the restriction of possible delays Cbetween adjacent sensors to the delay range between 0and 2samples. This can be done by restricting possible direction of arrival ϕor by changing adjacent sensor distance daccording to (9). The choice of the constraint was also discussed. The high-pass constraint seems to fit most applications. Bad performance for noise uncorrelated between sensors and samples was clarified by finding the Wiener solution giving the initial weights. This means that when the noise uncorrelated between sensors and samples dominates then the filter remains in its initial state. Thus adaptation process takes no effect. The problem of estimating SNRE was considered and the particular solution was given (41). 9. Acknowledgment This research was supported by grant No. MSM:210000012 Transdisciplinary research in biomedical engineering and grant GACR 102/02/0124 Voice Technologies for Support of Information Society. References [1] WERNER, S., APOLINARIO, J., Laakso, T. Multipleantenna CDMA mobile reception using constrained normalized adaptive algorithms. In Proceedings SBT/IEEE International Telecommunications Symposium ITS’98, Sao Paulo (Brazil), 1998, pp. 354-358. [2] FROST, O. L. An algorithm for linearly constrained adaptive array processing. In Proceedings of IEEE, 1972, vol. 60, no. 8, p. 926-934. [3] GRIFFITHS, L. J., JIM, C. W. An alternative approach to linearly constrained adaptive beamforming. IEEE Transactions on Anntenas and Propagation, 1982, vol. AP-30, p. 27-34. [4] JOHO, M., MOSCHYTZ, G. S. On the design of the target-signal filter in adaptive beamforming. In Proceedings of IEEE International Symposium on Circuits and Systems, 1998, vol. 5, p. 166-169. [5] FAHRANG-BOROUJENY, B. Adaptive Filters, Theory and Applications. Wiley and Sons, 1996. [6] ˇ STRUPL, M. Analysis and simulation of Frost beamformer algorithm. Master’s thesis, FEE CTU in Prague, 2002. [7] WIDROW, B., DUVALL, K. M., GOOCH, R. P., NEUMAN, W. C. Signal cancelation phenomena in adaptive antenas: causes and cures. IEEE Transactions on Antenas and Propagation, 1982, vol. AP-30, no. 3, p. 469478. [8] GANNOT, S.,BRUSHTEIN, D., WEINSTEIN, E. Beamforming for multichannel speech enhancement. In Proceedings of IWAENC’99, 1999, p. 96-99. [9] LASENBY, J., FITZGERALD, W. J. A Bayesian approach to high resolution beamforming. In IEE Proceedings-F, Radar and Signal Processing, 1991, vol. 138, no. 6.