The Influence of Random Element Displacement on DOA Estimates Obtained with (Khatri–Rao-)Root-MUSIC

Although a wide range of direction of arrival (DOA) estimation algorithms has been described for a diverse range of array configurations, no specific stochastic analysis framework has been established to assess the probability density function of the error on DOA estimates due to random errors in the array geometry. Therefore, we propose a stochastic collocation method that relies on a generalized polynomial chaos expansion to connect the statistical distribution of random position errors to the resulting distribution of the DOA estimates. We apply this technique to the conventional root-MUSIC and the Khatri-Rao-root-MUSIC methods. According to Monte-Carlo simulations, this novel approach yields a speedup by a factor of more than 100 in terms of CPU-time for a one-dimensional case and by a factor of 56 for a two-dimensional case.


Introduction
DOA estimation is a major application of the sensor array, since there are many real-world problems where an accurate estimate of the source direction is essential: for example, in radar, electroencephalogram, sonar and microphone array systems. The number of sources that can be resolved by an N -element uniform linear array using traditional subspace-based methods, such as MUSIC, is N − 1. Recently, more advanced methods have been presented for DOA estimation, for example, the Khatri-Rao (KR) product approach. By assuming quasi-stationary sources, this concept can identify up to 2N − 1 sources using an N -element uniform linear array, without computing higher-order statistics [1].
Irrespective of the DOA estimation algorithm applied, deviations between the actual steering vector and the presumed steering vector are unavoidable in most applications that rely on an array of sensors. Errors during fabrication, uncalibrated arrays and arrays undergoing deformations are among the potential factors that can contribute to such errors. Moreover, the arrays may be composed of sensor elements with directional radiation patterns [2][3][4]. In most papers on DOA estimation, the authors assume that there is no steering vector error or that steering vector deviations are deterministic and constant in time. This has resulted in the development of algorithms that estimate such errors and, subsequently, remove them by a calibration procedure [5][6][7][8][9]. In many cases, however, these errors are random and, therefore, statistical in nature. The motivation of our work is that most calibration methods are computationally expensive and/or time consuming. Moreover, to estimate the modeling errors in the array, often, certain assumptions have to be made, additional experiments have to be carried out and/or specific training sequences must be transmitted. Therefore, these techniques are difficult to apply in real-time applications or in systems where the random errors change quickly. Hence, it is interesting to take into account the effect of uncertainty in the placement of the sensor elements by means of a stochastic framework. Hence, there is a need for stochastic simulation tools that provide the statistical data to quantify the effect of the random displacement of these sensor elements on the DOA estimates.
Conventionally, one may use Monte-Carlo simulations to quantify the effect of position errors of one of the sensor elements to characterize the resulting distribution of the estimated DOAs. However, Monte-Carlo simulations require a large amount of realizations to accurately capture the statistics of the random process. Hence, the method is time consuming.
In this paper, we introduce the stochastic collocation method (SCM) [10] as a more effective way to rapidly model uncertainty in the estimated DOAs, due to variations in the elements of the steering vector. Thanks to this method, we can reduce the CPU time by a factor of more than 100 for a one-dimensional displacement and by a factor of 56 for a two-dimensional displacement of one of the sensor elements.
In this paper, we apply the nominal steering vector of a uniform linear sensor array to estimate the DOAs of a signal impinging on a sensor array where one of the sensors has a random position error. We establish a stochastic framework to find the probability density function (pdf) of the DOA-estimates. To the best of the authors' knowledge, such a framework has not been established in the open literature yet. In [11], a first-order sensitivity analysis is carried out to study the effect of modeling errors on the conventional MUSIC algorithm. This sensitivity analysis results in a sensitivity parameter and failure threshold, but does not provide the complete statistical distribution of the modeling errors. In [12], the authors examine a resolution threshold of the MUSIC algorithm for situations in which the array response is perturbed from its assumed value. This could be a boundary approximation for our framework. However, they assume that the arrays are illuminated with two equal power emitters.
In order to illustrate and validate the SCM method, we derive the statistical distribution on the DOA estimates for the root-MUSIC algorithm and for the KR-root-MUSIC algorithm. The latter algorithm applies the root-MUSIC method to the noise subspace matrix of the KR algorithm, as described in [1]. For the (KR)-root-MUSIC algorithm, we derive the cumulative distribution function (cdf), for one-and two-dimensional displacements of one of the sensor elements.
In earlier work, it has been proven that underdetermined DOA estimation is possible when the source signals are non-Gaussian stationary. In [13], by the use of fourth-order cumulants, a virtual array with an increase in the degrees of freedom is achieved. In [14], a new array geometry, which is capable of increasing the degrees of freedom of linear arrays, is proposed. The structure is obtained by systematically nesting two or more uniform linear arrays. It can provide N 2 degrees of freedom using only N physical sensors and the second-order statistics of the data. Blind source separation, by use of quasi-stationarity, has also received attention. One technique is based on the least squares fitting (LSF) criterion, using parallel factor analysis (PARAFAC) [15]. One can prove that quasi-stationarity enables the identification of sources when the number of sensors is lower than the number of sources. However, the LSF is based on a multi-dimensional non-linear optimization problem. In [16], the authors have used a method called focusing Khatri-Rao subspace (FKR) for wideband array processing. They calculate a focusing matrix using a rotational signal-subspace [17], and then, the covariance matrices of different frequencies are transformed and combined by means of the KR-product. Recently, in [18], the authors used sparse covariance fitting in order to estimate underdetermined DOA estimation, and in [19], a sparse representation of the array covariance vector is used to obtain the DOAs. In [20,21], the Khatri-Rao approach is also extended to a uniform circular array. In [22], the Khatri-Rao approach is used for 2D DOA-estimation. In [23], the authors make use of a maximum likelihood method to obtain the underestimated DOAs for a multiple-input multiple-output radar.
Notations: We denote matrices and vectors by boldfaced capital letters and lower-case letters, respectively. Superscript H denotes the transpose conjugate, whereas superscript T denotes the transpose without the conjugate. For a given vector x ∈ C, ||x|| denotes its Euclidean norm. For a given matrix A ∈ C M ×N , the i-th column is denoted by a i . The notation vec(A) stands for vectorization; i.e., if The remainder of the paper is organized as follows. In Section 2, we present the signal model and the KR-MUSIC algorithm [1]. Section 3 describes how the steering vector and the manifold array change if there is a displacement in one of the sensor elements. A brief description of the SCM, which is applied to model the uncertainty due to statistical variations of the input parameters, is given in Section 4. The proposed method is validated in Section 5, and conclusions are drawn in the last section.

Signal Model and Summary of the KR-MUSIC Algorithm
In this section, we give a short summary of the KR-MUSIC algorithm (see [1]), preceded by the signal model used.

Signal Model
Consider K sources impinging on a uniform linear sensor array, consisting of N omnidirectional sensors and neglecting mutual coupling.
For a signal s k (t) emitted by the k-th source, let s(t) = [s 1 (t), . . . , s K (t)] T denote the K × 1 source signal vector.
Let a(θ) = 1, e j2π d λ sin(θ) , . . . , e j2π d λ (N −1) sin(θ) T be the N × 1 steering vector corresponding to the direction θ, with d the spacing between the N sensors and λ the signal wavelength. In this paper, we fix the element spacing to d/λ = 1/2. Our aim is to study the effect of random variations in the array geometry on algorithms estimating the directions of arrival θ k corresponding to the K source signals, in the presence of measurement noise. Therefore, A = [a(θ 1 ) a(θ 2 ) . . . a(θ K )] describes the N × K array manifold matrix and v(t) = [v 1 (t), . . . , v N (t)] T ∈ C N ×1 represents the spatial and temporal white noise. For the observed signal x n (t) of the n-th sensor, let x(t) = [x 1 (t), . . . , x N (t)] T denote the N × 1 receiver signal vector. This vector equals: The DOA-estimation algorithms under study are the root-MUSC and KR-root-MUSIC algorithms. These algorithms rely on the following assumptions: • The noise v(t) is white Gaussian and zero-mean wide-sense stationary with covariance matrix C E{||v(t)|| 2 }. It is statistically independent of the source signals.
• All source DOAs are distinct: θ k = θ l for k = l.
• The source signals s k (t) are mutually uncorrelated and have a zero mean.
• Each source signal is wide-sense quasi-stationary with frame length L, such that: with d m,k the frame-dependent and sensor-dependent average normalized signal power. This assumption means that the second-order statistics of the source signals are time varying, but that they remain static over a short period of time. Under this quasi-stationarity assumption, we define a local covariance matrix: In the remainder of the paper, we apply the DOA estimation algorithms to synthetic signals generated by the procedure proposed in [1] to construct source signals s k (t) satisfying the above constraints, with sequence length T = 25, 600 and an allowable range of the frame periods [L low , L upp ] = [300, 700]. Algorithm 1 is applied to generate the synthetic signals.
Set T cur = 0. T cur is the start of a new frame of the signal; while T cur < T do Randomly generate L f following a uniform distribution on [L low , L upp ]; Randomly generate σ s following a uniform distribution on [0,1]; for t = T cur to T cur + L f − 1 do Randomly generate s k (t) = s R (t) + js I (t), where s R (t) and s I (t) are independent and identically Laplacian distributed with zero mean and variance σ 2 s /2; end T cur := T cur + L f ; end Algorithm 1: The algorithm to generate synthetic signals.
Observe that the frame intervals of the different impinging sources are not fixed and not synchronized. We apply the KR subspace methods by choosing a fixed frame period of L = 512. The number of frames is set to M = T /512 = 50.

Summary of the KR-MUSIC Algorithm
For a given receiver signal vector {x(t)} T −1 t=0 , with frame length L, where L divides T , the KR-MUSIC algorithm proceeds as follows [1]: • Compute the local variance estimateŝ • Denote an orthogonal complement projector P ⊥ , withŶ the data matrix formed by vectorize the covariance matrix. This projection eliminates the unknown noise covariance. • As in [1], perform a dimension reduction on Y obtainingỸ.
• To apply MUSIC toỸ, perform a singular value decomposition (SVD) on: where U ∈ C (2N −1)×K and V ∈ C M ×K are the left and right singular matrices, respectively, and Σ ∈ R K×K is the diagonal singular values matrix with the singular values being arranged in descending order. • By extracting the noise subspace matrix we compute the DOA spectrum: and b(θ) = b r (e j2π d λ sin(θ) ). The K largest local maxima of the DOA spectrum P KR−MUSIC (θ) are the DOA estimates.
Instead of searching the peaks of the MUSIC spectrum, we can also apply the root-MUSIC method [24] to ||U n H W|| 2 . In the remainder of the paper, we call the resulting method the KR-root-MUSIC method.

Displacement of One of the Sensor Elements
In this section, we assume that one of the sensor elements is displaced due to an unintentional movement of an element. Figure 1 depicts an array of N sensor elements, where the first sensor element is placed at the origin of the xy-plane and the other elements along the x-axis of this plane. The l-th sensor element is displaced by a vector r h with respect to its original position.  Denote r h = || r h || and α = arg( r h ). The new coordinates of the displaced l-th sensor are given by r = [(l − 1)d + r h cos(α), r h sin(α)]. Define β = arg( r) and r = || r||, such that 2π λ u · r = 2π λ r sin(θ + β). u is the unit vector with the same direction of the position vector of the source, which is located in the array's far field. The l-th row of the array manifold matrix then becomes e j2π r λ sin(θ 1 +β) , e j2π r λ sin(θ 2 +β) , . . . , e j2π r λ sin(θ K +β) .

Stochastic Collocation Method
In order to investigate the influence of a stochastic displacement of one of the sensor elements, we need stochastic simulation tools to characterize the statistical distribution of the DOAs due to the distribution of the random displacement vectors.
A Monte-Carlo simulation for such a task requires a large amount of realizations to accurately capture the statistics of the random process. Therefore, we apply the SCM proposed in [10] as a more effective way to rapidly model geometrical uncertainty, due to random variations of the array element positions.

One-Dimensional Stochastic Input
Assume that an input random variable X (in our case, e.g., the displacement in the xor y-direction) is given. We want to determine the variation of the output Y = f (X) (in this paper, the DOA estimates), due to the statistical variation of X. The random variable X follows the cumulative distribution P X and probability density function dP X in the sample space Ω. To determine the statistics of Y , we rely on the Askey scheme [25] to approximate the transformation Y = f (X) by a polynomial expansion of order P .
with φ X k (X) the expansion polynomials and y X k the weights belonging to the k-th expansion polynomial.
An optimal expansion is obtained when the set of expansion polynomials forms a complete orthogonal basis in Ω with orthogonality relation: In this case, the Cameron-Martin [26] convergence theorem ensures exponential convergence to the function Y = f (X) for P → ∞. In the remainder of the paper, we assume that the statistical variation of X follows a Gaussian distribution in the sample space Ω. Relying on the Askey scheme, we know that the probabilistic Hermite polynomials H k (X) provide an optimal expansion for the normal Gaussian distribution [25].
To determine the unknown expansion coefficients y X k , we apply Galerkin weighting to Equation (7) [27]: We can approximate this integral by a Q-point Gauss-Hermite quadrature rule, being [28]: where the quadrature points x i are given by the Q zeros of φ X Q (x) in Ω and with w i the corresponding weights. In order to evaluate Equation (10), we must determine Y = f (X) for Q realizations of the random variable X, corresponding to the quadrature points. When utilizing a Gaussian quadrature, Equation (10) exactly integrates all polynomials of degree equal or less than 2Q − 1. Given an expansion order P , the highest order coefficient evaluation of Equation (9) can be assumed to involve integrands of at least polynomial order 2P , being φ X (x) of order P and Y (x) modeled to order P , such that a minimal Gaussian quadrature order of P +1 will be required to obtain good accuracy in these coefficients. Given the exponential convergence of the expansion (7), using the SCM theory, we need to apply the KR-root-MUSIC algorithm in the Q = P + 1 quadrature points. The Monte-Carlo analysis may then be applied to the computationally cheap expansion (7) instead of to the KR-root-MUSIC. This results in huge savings in CPU-time, as the Monte-Carlo method is only slowly converging with respect to the number of samples taken, being inversely proportional to the square root of the number of samples [29].

Two-Dimensional Stochastic Input
We assume that the variation of the output Y = f (X 1 , X 2 ) depends on two independent stochastic variables X 1 and X 2 , respectively, in the sample space Ω 1 and Ω 2 , respectively (e.g., a displacement in the xand y-direction).
We first extend the GPoC expansion (7) to two dimensions. For Y = f (X 1 , X 2 ), we obtain: which is a polynomial expansion of degree P 1 in X 1 and P 2 in X 2 . The weight coefficients are given by: or if we approximate Equation (12) by a tensor product quadrature rule: with the quadrature points u k and u l given by the Q 1 zeros of φ X 1 Q 1 (x 1 ) and the Q 2 zeros of φ X 2 Q 2 (x 2 ), respectively. Equation (13) requires Q 1 · Q 2 function evaluations. If we only need a small number of quadrature points (whenever we can model Y as a polynomial of low degree), this is a very effective numerical tool. However, the number of function evaluations grows quadratically with the number of quadrature points Q (in the case that Q = Q 1 = Q 2 ).
In the case that X 1 and X 2 are independent Gaussian distributed variables with zero mean and dispersions σ 1 and σ 2 , respectively, we can rewrite Equation (12) as: with, We approximate this integral by using cubature formulas for the plane with weight function e −r 2 ( [30,31]): with (x 1,k , x 2,k ) the quadrature points and w k the corresponding weights. For these formulas, the number of function evaluations is limited to Q. We call d e the degree of the cubature formula, meaning that Equation (16) is exact for polynomials of degree d e .
The cubature formulas that we have applied in this paper are those from [32] with Q = 44, d e = 15, those from [33] with Q = 99, d e = 21 and those from [33] with Q = 172, d e = 31.
Given the exponential convergence of the polynomial expansions, for the two-dimensional expansion, we have to apply the KR-root-MUSIC algorithm only Q 1 · Q 2 times, when making use of the tensor rule Equation (13). Whereas, if we use the cubature formulas following Equation (16), we need to apply the KR-root-MUSIC algorithm Q times. Once we have estimated the DOAs in the quadrature points, the whole distribution can be obtained by applying the Monte-Carlo method to the computationally cheap polynomial expansion instead of to the DOA estimation algorithms.

Results
In all simulation examples below, the signal-to-noise ratio SNR (in dB) is defined as: In order to validate the efficiency of the SCM method, we will provide several numerical examples. We first consider the case where the displacement vector r h is a constant vector. Next, we focus on the more general case where the xand/or the y-coordinates of the displacement vector are random variables.
In Section 5.1, we investigate the influence on the DOA estimates obtained by root-MUSIC and KR-root-MUSIC for a displacement of one sensor element, in absence of noise. In Section 5.2, we use the SCM method to approximate the curves obtained in Section 5.1. In Section 5.3, we analyze the cdf for a displacement of the sensor element in one direction, for noisy data signals. In Section 5.4, we derive the GPoC expansion by means of the SCM method for a two-dimensional shift.

One Shifted Sensor Element: DOA Estimation in Absence of Noise
Consider the case of K = 3 sources and N = 4 sensor elements. The true DOAs are assumed to be [−18 • , 5 • , 25 • ]. We apply the synthetic signal generation procedure of Section 2.1 to generate a random source signal and use this signal vector throughout the simulations. We furthermore assume that no noise has been added to the impinging signal. We assume that one of the sensor elements has moved along the x-axis. For e rel x , being the relative displacement (with respect to the distance d between two sensor elements) of the sensor element, r h corresponds to coordinates [e rel x d, 0]. In Figures 2 and 3, we plot the estimated DOAs, as a function of e rel x for a relative displacement of each of the sensor elements for the root-MUSIC algorithm and the MUSIC algorithm. As a reference, the figures also show the true DOAs as dashed lines. We observe that even without any noise, the MUSIC algorithm is sensitive to a shift of one of the sensor elements and that if that shift is too large, in some cases, we do not find peaks in the MUSIC spectrum, such as in Figure 3 for e rel x > 0.2. On the same figure, we observe that the root-MUSIC algorithm [24] estimates all DOAs, even in these cases where the conventional MUSIC algorithm does not find enough peaks. Furthermore, we observe that the DOAs estimated by the root-MUSIC algorithm and those estimated by the MUSIC algorithm are almost the same. The figures for a displacement of the first and fourth array element and those for a displacement of the second and third array element will almost be symmetric. In Figures 4 and 5, we have plotted the estimated DOAs, as a function of e rel x , for a random relative displacement of the first and second sensor element for the KR-MUSIC and KR-root-MUSIC algorithm. In contrast to MUSIC, the KR-MUSIC algorithm is more robust to shifts in one of the sensor elements. The Khatri-Rao transformation expands the dimension of the vector on which we perform the SVD and the extracted noise subspace matrix. This increases the probability of finding K peaks. We also observe that in the KR-MUSIC case, the estimated curves are almost linear as a function of the shifts. Furthermore, we observe that DOAs estimated by using the KR-MUSIC algorithm are almost the same as those estimated by the KR-root-MUSIC algorithm.

One-Dimensional GPoC Expansion
In this section, we start from the same assumptions as in Section 5.1, concerning the impinging signal and the sensor array.
Let f (e rel x ) be the estimation of the second DOA (true DOAs = [−18 • , 5 • , 25 • ]) as a function of the relative displacement of the second sensor element along the x-axis e rel x . As plotted in the top figure of Figure 6, the relative displacement e rel x is assumed to be Gaussian distributed with dispersion σ pl = 0.08. As outlined in Section 4, we therefore apply Hermite polynomials in the GPoC expansion of f (e rel x ). We use the SCM method (with Q quadrature points and a polynomial expansion of order P = Q − 1) to determine the expansion coefficients for the GPoC expansion. In the bottom figure of Figure 6, we have plotted the function f (e rel x ) and the GPoC expansion of f (e rel x ). We remark that, for e rel x between −0.16 and 0.16, all plotted curves coincide within 1%. We also observe that, if we make use of a higher expansion order for the polynomials, the GPoC expansions better approximate the curve representing f (e rel x ). For a polynomial expansion of order 16, the GPoC curve coincides with the f (e rel x ) curve within 8% in the observed interval.  e rel In Figure 7, we apply the SCM method with P = Q − 1 for e rel x Gaussian distributed with dispersion σ pl = 0.2. We observe that, in the 95% confidence interval [−0.4, 0.4] for e rel x , the GPoC based on the Hermite polynomials of order 12 or higher provides an accurate approximation for f (e rel x ). However, beyond the 95% confidence interval, the Hermite polynomials rapidly oscillate, and they do not yield a good approximation for the f (e rel x ) curve. In Figure 8, we have made the same assumptions as in Figure 7, but we have estimated the DOA using the KR-root-MUSIC algorithm. The curves coincide, even when restricting the GPoC to only P = 4.

One Shifted Element: Statistical Distribution of DOA Estimates
In this section, we again apply the synthetic signal generation procedure of Section 2.1 to generate a source signal. We construct an arbitrary complex Gaussian distributed noise vector v(t) = v R (t) + v I (t). v R (t) and v I (t) are both Gaussian distributed with zero mean and dispersion σ = 0.3 providing an SNR = 7.78 dB. Throughout the simulations, we use this source signal vector and noise vector. We also assume that the first sensor element is displaced along the x-axis and that the relative displacement e rel x is Gaussian distributed with zero mean and dispersion σ pl = 0.12. To obtain the cdfs, we use two different methods: Monte-Carlo (MC) simulation and SCM. We carry out an MC simulation (where we vary e rel  We also apply the SCM theory for different values of the quadrature points Q and a polynomial expansion of Y of order P = Q−1. We perform a two-sample Kolmogorov-Smirnov test (KS test) [34], to compare the samples obtained by the MC method with those obtained by the SCM method. The two-sample KS test verifies the null hypothesis that the two samples come from the same distribution. The test statistic is defined as: with sup x (f (x)) the supremum of the function f (x). F 1,n (x) is the empirical cdf obtained with the MC and F 2,n (x) the empirical cdf obtained with SCM. n = n = 10, 000 are the number of samples used for the MC and SCM method. The null hypothesis is rejected at level α if D n,n > c(α) n+n nn , where, for example, α = 0.05 corresponds to c(α) = 1.36. In our case, the critical value is D n,n = 0.0192. In Tables 1 and 2, we show the values of the test statistic D n,n , for different values of the number of quadrature points Q and orders of polynomial expansion P = Q − 1. We consider the root-MUSIC algorithm and the KR-root-MUSIC algorithm, respectively. The two-sample KS test accepts, with a confidence of 95%, that, for the root-MUSIC algorithm and for P = Q − 1 = 4, the SCM samples and the MC samples have the same distribution. We also show the time necessary to obtain the cdf using the SCM method. Note that the simulation time to obtain the samples for the SCM algorithm is a factor more than 100 smaller than the simulation time to calculate the samples obtained by the MC method.
For the KR-root-MUSIC algorithm, both empirical distributions are the same for P = Q − 1 = 4. In Figures 9-11, we have plotted the cdfs both obtained by the SCM method (with Q = 9 quadrature points and a polynomial expansion of order P = 8) and the MC method. We observe that, even for a small number of quadrature points, we see almost no difference between the curves obtained by the SCM theory and the curves obtained by MC simulation, as indicated by the true sample KS test. Moreover, for the first and third estimated DOA, the root-MUSIC algorithm estimates the DOA very accurately. For the estimation of the second DOA, however, the KR-root-MUSIC algorithm performs better.

Two-Dimensional Displacement of One Sensor Element
In this section, we assume that the first sensor element is displaced by r h = [e rel x d, e rel y d] (Figure 1). The relative displacements in both the x-direction e rel x and y-direction e rel y are independently Gaussian distributed, with zero mean and dispersion σ pl . This means that the relative length of the displacement vector is Rayleigh distributed with mean σ pl π 2 and variance 4−π 2 σ 2 pl and that the argument of the displacement vector is uniformly distributed in [−π, π].  For the same impinging signal as in Section 5.3, we have plotted the cdfs of the different estimated DOAs, for σ pl = 0.12. We have carried out an MC simulation with 10,000 realizations. The simulation time is about 128 s for the root-MUSIC algorithm and 157 s for the KR-root-MUSIC algorithm. We also apply the two-dimensional GPoC expansion (11) with degree P in both variables to obtain the cdfs for the DOAs. We calculate the weight coefficients using Q 1 = Q 2 = P − 1 quadrature points in the xand y-directions using Equation (13). Tables 3 and 4 show the KS test to compare the samples obtained by the SCM method and the MC method. Furthermore, we compare this with the two-dimensional SCM theory that relies on the cubature formulas with Q = 44, Q = 99 and Q = 172, respectively, and with polynomial expansions of order P = 7, P = 10 and P = 15, respectively. The measured value of the KS test, which compares the samples of the MC method with those obtained with the cubature formulas, is also presented in Tables 3 and 4. Furthermore, the simulation times are shown in the last column of the tables. We observe that now, for the root-MUSIC algorithm, we need 14 quadrature points in each direction to satisfy the null hypothesis of the KS test (α = 0.05). Using cubature formulas, only for the estimates of the third DOA, 172 quadrature points are insufficient to satisfy the null hypothesis. For the KR-root-MUSIC algorithm, four quadrature points in each direction are sufficient. If we use cubature formulas, Q = 44 quadrature points are sufficient to satisfy the null hypothesis. We also observe that, using the SCM theory, the simulation time is reduced by a factor of about 20, and using the cubature formulas, we can reduce the simulation time by a factor of 56.
In Figures 12-14, we have plotted the empirical cdf for the MC simulation and for the SCM method P = Q − 1 = 13. We observe that both curves almost coincide.  , such that 68% of the estimates lie in an interval of length 2S around the true DOA. For a Gaussian distribution, the value S will correspond to the dispersion. In Figures 15 and 16, we have plotted the value S as a function of σ pl , for a two-dimensional displacement of the first sensor element. Again, we observe that the SCM method provides a good approximation for the estimates, both for the root-MUSIC algorithm as for the KR-root-MUSIC algorithm. Observe that for values σ pl less than 0.4, the root-MUSIC yields better estimates for the first and third DOA than the KR-root-MUSIC algorithm. However, the KR-root-MUSIC estimates the second DOA better. We furthermore observe that, for the KR-root-MUSIC algorithm, the S-values are straight lines as a function of σ pl .   Observe that, with the same array configuration as in Figure 16, we estimate twice as many DOAs. Looking at the S = Q 84 −Q 16 2 , we observe that the KR-root-MUSIC algorithm finds the first estimate forθ 1 = −65 • with a large error, even for small values of σ pl . The estimation error for the remainder of the DOAs starts to increase dramatically for values of σ pl larger than 0.2.

Conclusions
In this paper, we have investigated the effect of a displacement of one of the sensor elements. We have applied the SCM formalism to efficiently predict the cdfs of the estimated DOAs, in the case that the relative displacement is Gaussian distributed. Comparison of the SCM method to the Monte-Carlo method demonstrates that SCM is an accurate method to predict the cdfs of the DOAs. For the one-dimensional case, the SCM method decreases the simulation time by a factor of more than 100. For the two-dimensional case and using cubature formulas, we can reduce the simulation time by a factor of 56. In the future, we will extend the applied method to array configurations in which more sensor elements are displaced. Moreover, the SCM will be applied to other DOA estimation algorithms and for other array configurations. The flexibility and non-intrusive character of the proposed method only requires that DOA-estimates need to be computed in the quadrature points, by use of the appropriate algorithm.