Direction-of-Departure and Direction-of-Arrival Estimation Algorithm Based on Compressive Sensing: Data Fitting

: In this paper, a compressive sensing-based data fitting direction-of-departure/ direction-of-arrival (DOD/DOA) estimation algorithm is proposed to apply the superior performance of compressive sensing method to the bistatic MIMO sonar systems. The algorithm proposed in this paper optimizes the output data via convex optimization-based sparse recovery, so that it is possible to estimate the DOD and the DOA for each target accurately. In order to minimize the amount of computation, the cost function with constraint condition is implemented in this paper. Furthermore, the constraint condition parameter of the cost function is analytically derived. Through various simulations, it is shown that the superior DOD and DOA estimation performance of the proposed algorithm and that the analytical derivation of the constraint condition parameter is useful for determination of regularization parameter.

The compressive sensing CS-based DOA estimation method is a technique for estimating DOAs by using the sparse property of a signal when incident signals are spatially sparse. Figure 1 shows the concept of the CS-based DOA estimation method.
The received signals used for DOA estimation methods are incident on the sensor array at specific angles. Therefore, the received signals can be said to be spatially sparse as the signals are only incident in specific angles, not in all the angles.
In [4], a method of extending single measurement or multiple measurements to beam space and performing sparse recovery using basis pursuit denoising (BPDN) has been presented. Furthermore, in order to solve the computational problem that occurs when using multiple measurements, a L1-SVD method that performs sparse recovery using singular vectors of multiple measurements has been proposed.
The convex optimization based sparse recovery used in [4] requires a large amount of computation. In [5], to decrease the computation of the sparse recovery, a method of approximating the l 0 norm In [6][7][8], compressive sensing-based covariance fitting methods have been studied. In [9], a compressive sensing-based covariance fitting algorithm using iterative operation sparse iterative covariance-based estimation (SPICE) has been presented. As SPICE does not perform sparse recovery using BPDN, it requires less computation that the existing covariance fitting algorithms [6][7][8].
In [10], a sparse spectrum fitting (SpSF) DOA estimation algorithm is introduced and its effective regularization has been studied. A weight used in the cost function of SpSF depends on the number of snapshots. In [11], an optimal weight of the cost function setting method according to the number of snapshots has been presented.
The phase error of the uniform linear sensor array has a great effect on the estimation performance of the DOA estimation algorithms. In [12], a new compressive sensing-based DOA estimation algorithm that repeatedly estimates the phase error and DOA has been proposed.
The conventional compressive sensing-based DOA estimation methods depend on the noise field. In [13], to make it possible to estimate DOA based on the compressive sensing method even in the unknown noise field, a new compressive sensing-based DOA estimation algorithm based on the generalized correlation decomposition algorithm has been presented.
The CS-based DOA estimation method is utilized for estimating DOAs by spectralizing the signal data vector itself using the property that the received signal is spatially sparse. In addition, the CS-based DOA estimation method estimates DOAs by optimizing a sparse output signal data with two constraints of sparsity reduction and sensor data model fitting. Therefore, the CS-based DOA estimation is more robust to sensor elements failure in comparison with the non-CS-based DOA algorithm. In addition, the performance of the CS-based DOA estimation algorithm does not degrade significantly with the decrease of sensor elements spacing and the decrease of the number of sensor elements in comparison with non-CS-based algorithms.
Target strength (TS) technology is currently being developed to avoid detection of the passive sonar system. Geometric modeling to minimize the target strength has been applied to the submarines and vessels. This technique significantly degrades the performance of the existing DOA estimation algorithms.
The bistatic multiple input multiple output (MIMO) sonar system is the solution to the problem of detecting the low TS submarines or vessels. This is a system in which two sensor arrays are divided and performed by a transmitter and a receiver, rather than performing both transmission and reception in a sensor array. When using the bistatic MIMO sonar system, the receiver estimates the direction-of-departure (DOD) and the DOA. The two sensor arrays can change roles and can detect a low TS target of the stealth technology by estimating the DOD and DOA.
In order to detect the low TS submarines or vessels, a number of DOA estimation algorithms implemented in bistatic MIMO sonar systems have been extensively studied. The conventional beamforming, the adaptive beamforming and MUSIC, which are representative DOA estimation algorithms, have been expanded to the DOD/DOA estimation algorithms [14][15][16][17][18][19].
In [14], a DOD/DOA estimation method has been proposed by applying ESPRIT to the bistatic MIMO radar system. In [15,16], MUSIC has been applied to the bistatic MIMO system and a method to estimate DOD and DOA separately has been proposed.
In [17], a joint DOD/DOA estimation method using ESPRIT and SVD of an array covariance matrix form transmit subarrays has been presented. The proposed scheme in [17] can estimate the DODs and DOAs of targets in closed-form and pair automatically.
In [18], the bistatic co-prime multiple input and multiple output arrays have been considered to estimate multiple far-field targets. Using the scheme presented in [18], many DOA estimation methods for ULA can be applied to the sparse bistatic MIMO array and with closed-form degree of freedom (DOF) presented in [18]; it can have an advantage in terms of array aperture.
In [19], to estimate the DOD and DOA and the Doppler shift of a moving target for a MIMO radar system, a reduced dimension and low complexity algorithm has been proposed.
In order to take advantage of the superior performance of the DOA estimation based on compressive sensing in the bistatic MIMO sonar system, this paper proposes a compressive sensing-based DOD/DOA estimation algorithm in the bistatic MIMO sonar system.
Among the compressive sensing-based algorithms, it is shown that the DOD/DOA estimation is possible in the case of orthogonal matching pursuit (OMP), which is a greedy algorithm [20]. However, in the case of greedy algorithms, there is a condition that the received signal is spare. That is, the number of targets must be known in advance. The proposed algorithm in this paper optimizes the output data as sparse by searching the solution through convex optimization, so that it is possible to estimate the DOD and the DOA for each target accurately. It has been shown that CS-based DOA estimation algorithms are superior to non-CS-based DOA estimation algorithms at the expense of much computational burden. In this paper, it is shown how the CS-based DOD/DOA estimation algorithm outperforms non-CS-based DOD/DOA estimation algorithms. The reduced-dimension cost functions based on fast Fourier transform (FFT) for moving targets was derived.

Signal Modeling of Bistatic MIMO Sonar System
In this section, the bistatic MIMO sonar signal model [14] used in this paper is presented. The transmit and receive arrays are uniform linear array (ULA). The transmitter array has M elements and receiver array has N elements. λ denotes the wavelength. Let d t and d r be the inter-element spacing at the receiver and transmitter, respectively. The number of targets is P. Assume that the DOD and the DOA for the p-th target be expressed as θ tp and θ rp , respectively. The array vector for the n-th sensor of receiver can be written as a n θ rp = exp j 2π The array vector for the m-th sensor of transmitter can be written as Using (1), the receiver steering matrix A r can be written as Using (2), the transmitter steering matrix A t can be written as In the case of P targets, α is defined as where α 1 , · · · , α P are the target strengths of each target. A coded pulse of the transmitter can be expressed as a matrix S: where s i is the transmitted signal from i-th sensor of the transmitter and s j is the transmitted signal from j-th sensor of the transmitter. These vectors are designed as orthogonal to each other to distinguish signals of each of the transmitters after matched filtering. The number of snapshots is L: where N is the matrix of noise. The noise is complex-valued and the real part and the imaginary part of the noise are independent Gaussian random variables with zero mean. The variance of the real part is denoted by σ 2 /2, which is equal to the variance of the imaginary part, σ 2 /2. The received signal matrix can be written as After matched filtering the signal, the output of the filter can be expressed as By transforming the matched filtered signal to an independent sufficient statistic vector, the signal data becomes applicable to the DOD/DOA estimation methods. In case the number of targets is one, the independent sufficient statistic vector can be expressed as The signal covariance matrix R s can be written as The array vector reflecting the bistatic MIMO sonar system can be expressed as row(·) denotes the operator that stacks the rows of a matrix in a column vector. The noise vector passed through the matched filter can be written as When the number of targets is P (more than two), the array manifold of the bistatic MIMO sonar system can be written as When Q transmit pulses are used, the target strengths matrix H can be expressed as In this case, the matrix of the noise can be written as Only stationary targets are considered in this paper. Therefore the target strengths are constant in a pulse, and the target strengths are independent from pulse to pulse, rather than scan to scan.
When the number of targets is P (more than two) and multiple transmit pulses are used, the independent sufficient statistic matrix Y η can be expressed as When a single transmit pulse is used, the independent sufficient statistic vector y η can be expressed as

Bistatic Data Fitting DOD/DOA Estimation Algorithm
As in (18), in case of using one pulse, the signal data used for the DOD/DOA estimation can be obtained in the form of a vector in bistatic MIMO sonar system even if the transmitter sent a coded pulse into multiple snapshots. When only one pulse is used, the cost functions of the compressive sensing-based data fitting DOD/DOA estimation algorithm in bistatic MIMO sonar system can be obtained by applying the compressive sensing-based data fitting DOA estimation algorithm using a single snapshot [4].
The cost function of compressive sensing-based data fitting DOD/DOA estimation algorithm for determination of h without restrict condition can be defined as In the case of (19), the weight parameter w controls the trade-off between the sparsity of the estimated data h 1 and the difference between y η and Kh according to a given value.
By repeating the optimization process through the cost function of (19), we find the sparse signal covariance vector among countless solutions.
In the case of (19), as the convex optimization needs to be performed many times while changing the weight value w, it requires a large amount of operation time.
The cost function of compressive sensing-based data fitting DOD/DOA estimation algorithm for determination of h with constraint condition can be defined as Minimizing the number of non-zero elements of the signal vector Parameter of the constraint condition In the case of (20), it is possible to obtain a sparse output value through a single optimization. To use this algorithm, the parameter of the constraint condition β 2 needs to be determined.
In [4], the parameter of the constraint condition of the compressive sensing-based data fitting DOA estimation algorithm has been derived by using the fact that the real part and the imaginary part of the noise are independent Gaussian random variable with zero mean.
On the other hand, in the case of the bistatic MIMO sonar system, as the matched filter is used, the characteristic of the complex Gaussian noise can not be utilized as in [4].
To achieve the parameter of the constraint condition β 2 of (20), the mean of v 2 2 and standard deviation of v 2 2 are used in this paper.
The parameter of the constraint condition β 2 can be written as In the case of the constraint condition parameter β 2 , it should be able to accommodate most random variables of v 2 2 . Usually, the mean of v 2 2 + 6· and standard deviation of v 2 2 cover most of the its distribution, so the parameter of the constraint condition β 2 has been set as in (21).
The mean of v 2 2 and standard deviation of v 2 2 are derived in the Appendix A.

Numerical Results
This section shows numerical results for the two cases: The first case shows the accuracy of the derived parameter of the constraint condition. To show the superior performance of the DOD/DOA estimation of the proposed scheme, the second case illustrates the spectrum of the proposed scheme and the spectrum of conventional DOD/DOA estimation method. The simulation structure based on the bistatic sonar system is illustrated in Figure 2.

Checking the Analytically Derived Parameter of the Constraint Condition
The parameter of the constraint condition is dependent on the number of transmitters, the number of receivers, and the number of snapshots. Therefore, in order to show the accuracy of the analytically derived parameter of the constraint condition, based on various cases of the three parameters, the results of the analytic derivation have been compared with results of Monte Carlo simulation. Simulations have been made 10,000 times for the Monte Carlo simulation. Figures 3-6 show that the mean of v 2 2 , standard derivation of v 2 2 , and the parameter of the constraint condition β 2 derived in this paper with respect to the number of snapshots when the number of elements of the transmitter and receiver are variously set. In order to show the validity of each result of the analytic derivation, the results of each analytic derivation are compared with the results obtained by the Monte Carlo simulation. Figures 3-6 show that the analytically derived β 2 , which is the parameter of the constraint condition, and β 2 obtained from the Monte Carlo simulation are the same for all cases. In Section 4.2, using (20) and (21), the proposed scheme has been implemented by using MATLAB to show the performance of the DOD and DOA estimation.

Performance of the Bistatic Data Fitting DOD/DOA Estimation Algorithm
Figures 7-10 illustrate the spectrums of the proposed scheme and the bistatic conventional beamforming method when multiple targets (more than three) exist. Figure 7 shows the DOD/DOA estimation result of the proposed scheme based on the simulation parameters in Table 1. The number of targets is set to 4, and the SNR of each target is set differently. The figure shows that the proposed scheme accurately estimates the DODs and DOAs of the four targets.
The spectrum of the bistatic conventional beamforming algorithm based on the simulation parameters on Table 1 is illustrated in Figure 8. In the case of the beam forming method, as the main lobe and side lobe exist, the broad spectrum is usually formed. In Figure 8, it can be seen that the broad main lobe is formed around target 4, which has the highest SNR. However, due to side lobes, the DODs and DOAs of the other three targets cannot be estimated. Figure 9 illustrates the spectrum of the proposed scheme based on the simulation parameters in Table 2. The number of targets is four, and the SNR of each target is set differently. In Figure 9, it can be seen that the proposed scheme accurately estimates the DODs and DOAs of the four targets which have different SNRs to each other.

Parameter Value
The  Table 2. Simulation parameters of multiple targets (case 2).

Parameter Value
The  Figure 10 shows the spectrum of the bistatic conventional beamforming algorithm based on the simulation parameters on Table 2. In Figure 10, it can be seen that broad main lobe is formed around target 1, which has the highest SNR. However, due to side lobes, the DODs and DOAs of the other three targets cannot be estimated.  The spectrum of the proposed scheme is formed into a line spectrum due to sparse recovery. Therefore, the proposed scheme has superior resolution and can accurately estimate all four targets' DODs and DOAs. In addition, the spectrums of the proposed scheme in Figures 7 and 9 prove that the parameter of the constraint condition derived in this paper are appropriate.
In order to clearly compare the performance of the proposed scheme and the conventional method, the simulations were conducted under poor conditions in which the SNRs were 0 and −5 dB and the locations of the two targets were very close as tabulated in Table 3.
The spectrum of the proposed DOD/DOA estimation algorithm when the SNR is 0 dB is illustrated in Figure 11. The figure shows that the proposed scheme accurately estimates the DODs and DOAs of the two adjacent targets. The difference of the noise floor and the peak values of the two targets is greater than 200 dB. Figure 12 shows the spectrum of the conventional bistatic beamforming D0D/DOA estimation algorithm for SNR = 0 dB. Comparing Figures 11 and 12, Figure 12 shows that many spurious peaks occur in this case. Moreover, unlike the proposed scheme, the conventional method shows poor resolutions that cannot resolve two adjacent targets.
In Figure 13, the spectrum of the proposed DOD/DOA estimation algorithm for SNR = −5 dB is shown. Although the SNR in Figure 13 is worse than that in Figure 11, the spectrum of Figure  13 shows that the proposed scheme accurately estimates the DODs and DOAs of the two adjacent targets. Furthermore, the difference of the noise floor and the peak values of the two targets is greater than 200 dB. Table 3. Simulation parameters for two adjacent targets.

Parameter Value
The    Figure 14 illustrates the spectrum of the bistatic conventional beamforming method for SNR = −5 dB. The result of Figure 14 shows that the conventional method can not estimate the DODs and DOAs of the two adjacent targets. Figure 15 shows the operation times of the proposed scheme using (20); the cost function, which uses the parameter (β 2 ) of the constraint condition derived in this paper; and the scheme using the cost function without constraint condition in (19). The operation times of the two ways depend on the number of total elements of the transmitter and receiver. The simulation parameters of Figure 15 are in Table 4.
It can be confirmed using Figure 15 that the operation times increase as the number of receivers increases both in the proposed scheme using (20) and the other scheme using (19).
In the case of the scheme using (19), as the convex optimization needs to be performed many times while changing the weight value w, it requires a large amount of computation. On the other hand, in the case of the proposed scheme using the parameter (β 2 ) of the constraint condition derived in this paper, as the convex optimization is performed only once, the amount of computation can be minimized. Figure 15 illustrates that using the proposed scheme in (20) is much less computationally intensive than the other scheme using the cost function in (19), which justifies why the proposed scheme using the parameter (β 2 ) of the constraint condition should be employed for DOD/DOA estimation.  Eq. 19 Figure 15. The operation times Comparison between the proposed scheme using the parameter (β 2 ) of the constraint condition and the other scheme without using the constraint condition.

Conclusions
In this paper, a compressive sensing-based data fitting DOD/DOA estimation algorithm has been proposed to apply the superior performance of the compressive sensing method to the bistatic MIMO sonar systems. The proposed algorithm in this paper optimizes the output data as sparse by searching the solution through convex optimization, so that it is possible to estimate the DOD and the DOA for each target accurately.
Through simulation, the DOD/DOA estimation performance of two adjacent targets and multiple targets (more than three) has been compared with the bistatic conventional beamforming algorithm, and in all cases, the estimation performance of the proposed algorithm has been confirmed to be superior to that of the bistatic conventional beamforming algorithm.
In order to minimize the amount of computation, the cost function with constraint condition has been implemented in this paper. Furthermore, the constraint condition parameter of the cost function has been analytically derived. The validity of the parameter (β 2 ) of the constraint condition analytically derived in this paper has been shown by comparison with the value obtained from Monte Carlo simulation. The operation times of the cost function with the constraint condition and the cost function without the constraint condition have been compared. Through the numerical result, it has been confirmed that the proposed scheme in this paper is computationally efficient.
Through various simulations, it has been shown that the proposed scheme is superior to the previous DOD/DOA estimation algorithms and that the value of the regularization parameter can be analytically derived. The k-th row and l-th column element of the v can be written as E v 2 2 can be expressed as E v k,l 2 can be expressed as 1) s l (1) * + · · · + n k (L) s l (L) * n k (1) * s l (1) + · · · + n k (L) * s l (L) E v k,l 2 can be expanded to a term of sum of conjugate multiplication with oneself, A, and a term of sum of conjugate multiplication with each other, B.
The mean of A can be obtained as follows. At first, E n k (i) s l (i) * n k (i) * s l (i) can be written as E n k (i) s l (i) * n k (i) * s l (i) Using the fact that the coded pulse matrix S is based on Hadamard code and that every L cases of conjugate multiplication with oneself has the same value, the mean of A can be achieved as The mean of B can be obtained as follows. At first, E n k (i) s l (i) * n k (i ) * s l (i ) can be derived as i and i denote the different indices of the snapshots to each other. Using the fact that every L 2 − L cases of conjugate multiplication each other has the same value, the mean of B can be written as Using (A8) and (A10), E v 2 2 can be derived as can be written as To evaluate E v 2 2 2 , (A12) can be classified into two parts as follows, (A13) can be expressed as Derivation of E A 2 (A17) is an example of cases of the squared term.
E n k (l) s k (l) * n k (l) * s k (l) n k (l) s k (l) * n k (l) * s k (l) = |s k (l)| 2 2 · E |n k (l)| 2 2 = |s k (l)| 2 2 · E n Re k (l) 2 + n Im k (l) 2 2 = |s k (l)| 2 2 · 2σ 4 (A17) By using S is based on Hadamard code, (A17) can be simplified to An example of cases of the cross term of multiplication can be written as E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k l * s k l = E |n k (l)| 2 |s k (l)| 2 n k l 2 s k l 2 = |s k (l)| 2 s k l 2 E |n k (l)| 2 n k l Using |s k (l)| 2 = 1 and |s k (l )| 2 = 1, (A19) can be simplified to E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k l * s k l = |s k (l)| 2 s k l 2 σ 4 = σ 4 (A20) The number of squared terms is L and the number of cross terms of multiplication is L 2 − L. Therefore, E A 2 can be written as

Derivation of E [2AB]
Each term of E [2AB] can be classified into three cases: Case 1 = E n k (l) s k (l) * n k (l) * s k (l) n k (l) s k (l) * n k l * s k l Case 2 = E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k (l) * s k (l) Case 3 = E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k l * s k l , l = l = l (A22) By using E n k (l ) * = 0, case 1 can be written as E n k (l) s k (l) * n k (l) * s k (l) n k (l) s k (l) * n k l * s k l = E |n k (l)| 2 |s k (l)| 2 n k (l) s k (l) * n k l * s k l = |s k (l)| 2 s k (l) * s k l E |n k (l)| 2 n k (l) n k l * = |s k (l)| 2 s k (l) * s k l E |n k (l)| 2 n k (l) E n k l * = |s k (l)| 2 s k (l) * s k l E |n k (l)| 2 n k (l) · 0 = 0 (A23) By using E [n k (l )] = 0, case 2 can be written as E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k (l) * s k (l) = E |n k (l)| 2 |s k (l)| 2 n k l s k l * n k (l) * s k (l) = |s k (l)| 2 s k l * s k (l) E |n k (l)| 2 n k l n k (l) * = |s k (l)| 2 s k l * s k (l) E |n k (l)| 2 n k (l) * E n k l = |s k (l)| 2 s k l * s k (l) E |n k (l)| 2 n k (l) * · 0 = 0 (A24) By using E [n k (l )] = 0 and E [n k (l )] * = 0, case 3 can be written as E n k (l) s k (l) * n k (l) * s k (l) n k l s k l * n k l * s k l = E |n k (l)| 2 |s k (l)| 2 n k l s k l * n k l * s k l = |s k (l)| 2 s k (l) * s k l E |n k (l)| 2 n k l n k l * = |s k (l)| 2 s k (l) * s k l E |n k (l)| 2 E n k l E n k l * is identically zero for all cases: Derivation of E B 2 There are two cases in E B 2 . The example of the first case can be written as E n k (l) s k (l) * n k l * s k l n k (l) s k (l) * n k l * s k l = s k (l) * s k l 2 E n k (l) 2 n k l * 2 = s k (l) * s k l 2 E n Re k (l) + jn Im k (l) The second case can be subdivided into five cases: Case I E n k (l) s k (l) * n k l * s k l n k l s k l * n k (l) * s k (l) Case II E n k (l) s k (l) * n k l * s k l n k (l) s k (l) * n k l * s k l Case III E n k (l) s k (l) * n k l * s k l n k l s k l * n k (l) * s k (l) Case IV E n k (l) s k (l) * n k l * s k l n k l s k l * n k l * s k l Case V E n k (l) s k (l) * n k l * s k l n k l s k l * n k l * s k l (A28) By using |s k (l)| 2 = 1 and |s k (l )| 2 = 1, case I can be derived as E n k (l) s k (l) * n k l * s k l n k l s k l * n k (l) * s k (l) As the elements of the noise matrix are independent and identically distributed (IID) random variables, the results of cases II-V in (A28) turn to zero.
The number of terms of E B 2 is L 2 − L . Therefore, E B 2 can be written as From (A21), (A26), and (A30), (A15) is simplified to After expanding (A14), each term of (A14) can be written as For k = k It can be shown that E v k,l 2 v k ,l 2 is constant as long as k is equal to k irrespective of whether l is equal to l or not. For example, for evaluating E |v 1,1 | 2 |v 1,2 | 2 + |v 1,3 | 2 + · · · + |v N,M | 2 , we can use E |v 1,1 | 2 |v 1,1 | 2 = E |v 1,1 | 2 |v 1,2 | 2 = · · · = E |v 1,1 | 2 |v 1,M | 2 . (A33) The number of these terms is M (M − 1) for k (k = 1, · · · , N)-th sensor . By using (A21) and (A26), E [AA ], E [AB ], E [BA ] of E v k,l 2 v k ,l 2 for k = k , l = l and k = k , l = l can be simplified to In the case of E [BB ], terms of coded pulse's elements are not 1. It can be shown as E n k (l) s k (l) * n k l * s k l n k l s k l * n k (l) * s k (l) = s k (l) * s k l s k l * s k (l) E n k (l) n k l * n k l n k (l) * = s k (l) * s k l s k l * s k (l) E |n k (l)| 2 n k l 2 = s k (l) * s k l s k l * s k (l) σ 4 , s k (l) * s k l s k l * s k (l) = 1. (A36) l takes L − 1 values of 1, · · · , l − 1, l + 1, · · · , L. When l is fixed and the number of snapshots, L, is even, the sum of L − 1 cases of s k (l) * s k (l ) s k (l ) * s k (l) is −1. For k = k For k = k , (A38), an example case of the remaining E [AA ] of (A32) can be written as E |n k (l)| 2 |s k (l)| 2 |n k (l)| 2 |s k (l)| 2 = |s k (l)| 2 |s k (l)| 2 E |n k (l)| 2 |n k (l)| 2 = E |n k (l)| 2 E |n k (l)| 2 = σ 2 /2 + σ 2 /2 σ 2 /2 + σ 2 /2 = σ 4 . (A39) The number of terms is L 2 . The E [AA ] for k = k can be written as Expectation of a term in E [AB ] and E [BA ] For k = k can be written as E |n k (l)| 2 |s k (l)| 2 n k (l) s k (l) * n k l * s k l = |s 1 (1)| 2 s k (l) * s k l E |n k (l)| 2 n k (l) n k l * = |s 1 (1)| 2 s k (l) * s k l E |n 1 (1)| 2 E [n k (l)] E n k l = 0. (A41) Expected value of a term in E [BB ] of (A32) For k = k can be written as E n k (l) s k (l) * n k l * s k l n k (l) s k (l) * n k l * s k l = s k (l) * s k l s k (l) * s k l E n k (l) n k l * n k (l) n k l * = s k (l) * s k l s k (l) * s k l E [n k (l)] E n k l * E [n k (l)] E n k l * = 0. (A42) The sum of terms of (A14) for k = k can be given as