EURASIP Journal on Applied Signal Processing 2005:1, 61–68 c ○ 2005 Hindawi Publishing Corporation On Improvements of Cyclic MUSIC

Many man-made signals encountered in communications exhibit cyclostationarity. By exploiting cyclostationarity, cyclic MUSIC has been shown to be able to separate signals with different cycle frequencies, thus, to be able to perform signal selective direction of-arrival (DOA) estimation. However, as will be shown in this paper, the DOA estimation of cyclic MUSIC is actually biased. We show in this paper that by properly choosing the frequency for evaluating the steering vector, the bias of DOA estimation can be substantially reduced and the performance can be improved. Furthermore, we propose another algorithm exploiting cyclic conjugate correlation to further improve the performance of DOA estimation. Simulation results show the effectiveness of both of our methods.


INTRODUCTION
Among subspace-based DOA estimation methods, MUSIC [1,2] is relatively simple and effective, and is thus the most studied. To improve the performance of the conventional MUSIC, cyclic MUSIC [3] is shown to be effective to combat noise and interference by exploiting cyclostationary property possessed by most man-made signals in communications [4]. From then on, many papers have been developed to improve the performance of cyclic MUSIC, such as [5,6,7]. Reference [5] is asymptotically exact for either narrowband or wideband sources, but as stated in [5], as long as the crosscyclic correlations of the sources are not small enough, they will act consistently as small induced interferences with the same cycle frequency, especially when the source is narrowband. Thus, this method works better for the wideband case. By narrowband, we mean that the signal "fractional bandwidth," that is, the ratio of the signal bandwidth over its carrier frequency, is small. Reference [6] exploits both cyclostationarity and conjugate cyclostationarity to improve the ability of separating two closely impinging signals, and in certain condition, it can detect more signals than the number of antennas. Reference [7] makes use of spatial smoothing (SS) [8], and presents a scheme called Hankel approximation method (HAM). Then in conjunction with cyclic MUSIC, it solves the DOA estimation problem in the presence of coherent signals. Both [6,7] are based on the assumption that the signals are narrowband and the model x(t) = As(t) + n(t) holds exactly. Many papers refer to this model as a narrowband model; whereas, it actually holds exactly only for pure sinusoidal signals (see more details in Section 3). Thus in this paper, we will refer to this model as a sinusoidal model. Based on the fact that this model does not hold exactly for narrowband signals, we will show in this paper that cyclic MUSIC is actually biased and will provide appropriate ways to reduce the bias and improve the performance of cyclic MUSIC. Note that [6,7] did not provide any solution or improvement for this problem.
In this paper, we first propose an improved cyclic MU-SIC algorithm for narrowband signals. We assume that the sinusoidal model does not hold exactly and we calculate the cyclic correlation matrix of the exact received signal. Then by simplification and approximation, we can write the cyclic correlation matrix in a form, which is like that of the cyclic MUSIC, but the steering vector will be evaluated at f 0 + α/2 instead of f 0 as in cyclic MUSIC, where f 0 is the carrier frequency and α is the cycle frequency. Analysis and simulation results show that this method is effective in reducing the bias and improving the performance of the DOA estimation.
Next, by further studying the approximation we made in improved cyclic MUSIC, we notice that if we calculated the cyclic conjugate correlation matrix of the received signal, the error caused by approximation will be odd symmetric.
Then by transposing the matrix and averaging them, the error will be approximately canceled. We name this algorithm as improved conjugate cyclic MUSIC. Simulation results show that improved conjugate cyclic MUSIC performs even better than improved cyclic MUSIC. This paper is organized as follows. In Section 2, we briefly review the signal cyclostationarity and the existing cyclic MUSIC algorithm. Then the improved cyclic MUSIC and improved conjugate cyclic MUSIC algorithms are presented and discussed in Sections 3 and 4, respectively. Simulation results are shown in Section 5 and finally, we conclude this paper in Section 6.

Cyclostationarity
For a given signal s(t), its cyclic correlation and cyclic conjugate correlation are defined as [9] r α where (·) * denotes complex conjugate and · denotes infinite time average. Then, s(t) is said to be cyclostationary if r α ss (τ) or r α ss * (τ) is not zero at some time delay τ and cycle frequency α. Many man-made communication signals exhibit cyclostationarity due to modulation, periodic gating, and so forth. They usually have cycle frequency at twice the carrier frequency or multiples of the baud rate, or combinations of these. Moreover, some signals may have both nonzero cyclic correlation and nonzero cyclic conjugate correlation.
For a given vector x(t), we can calculate the cyclic correlation matrix and the cyclic conjugate correlation matrix which are defined as [3] where · H denotes conjugate transpose and · T denotes transpose.

Cyclic MUSIC
Consider a uniform linear array (ULA) of size N that receives I s signals from directions θ i , i = 1, . . . , I s . The incident waves are assumed to be narrowband plane waves from far-field sources. Without losing generality, we assume that among these signals, only the first I signals have cycle frequency α, which is of interest. Other signals either have different cycle frequencies or do not exhibit cyclostationarity. If we define where x n (t) is the signal received at the nth antenna for n = 1, . . . , N, then, for the narrowband case, x(t) can be approximated by the following form which will be referred to as the sinusoidal model: x(t) = As(t) + n(t). (5) In this model, the vector s(t) = [s 1 (t), . . . , s I (t)] T contains the first I signals that have cycle frequency α, that is, the signals-of-interest (SOI), the vector n(t) contains interfering sources and noise, and the matrix A = [a(θ 1 ), . . . , a(θ I )] contains the steering vectors defined as a θ i = 1, e j2π f0d sin θi/c , . . . , e j2π f0(N−1)d sin θi/c T (6) for i = 1, . . . , I, where d is the antenna spacing, c is the speed of propagation, and f 0 is the carrier frequency of the narrowband signals. Instead of calculating the correlation matrix of the received signal in the conventional MUSIC, cyclic MUSIC calculates the cyclic correlation matrix which is estimated by (2). Substituting (5) into (2), the cyclic correlation matrix can be approximated by As evaluating the cyclic correlation at α retains only those SOI, n(t) can be safely neglected. Then applying SVD to (7), we obtain where E s spans the signal subspace and E n spans the interference and noise subspace. Then, cyclic MUSIC finds the maxima of 1/ E H n a(θ) 2 or E H s a(θ) 2 in terms of θ as the DOA estimate.

Algorithm
In order to analyze the existing cyclic MUSIC algorithm and develop our improved cyclic MUSIC algorithm, we again consider a ULA of size N that receives I signals with the cycle frequency α, and other interference signals. If we choose the first antenna as the reference element, and assume the signal induced on this element due to the ith narrowband source being written in a complex form as where a i (t) is the amplitude and f 0 is the carrier frequency, then the signal induced on the nth element due to the ith narrowband source will be where ∆ i = d sin θ i /c. Adding all the sources together, the signal induced on the nth antenna will be where n n (t) contains interference signals and noise.
On the other hand, the nth entry of (5) can be written as Comparing (11) with (12), we notice that the following approximation is made: Thus, we conclude that (5) does not hold exactly for narrowband signals unless the signals are pure sinusoidal, that is, a i (t) is constant. In our analysis for improved cyclic MUSIC, instead of approximating the signals to fit the sinusoidal model, we start with calculating the cyclic correlation of x(t) of the exact model. Thus, the elements of x(t) are given in (11), not in (12). Then, the (p, n)th element of R α xx (τ) can be written as where the shift property of cyclic correlation is applied, that is, if y(t) = x(t + T), then r α yy (τ) = r α xx (τ)e j2παT , n n (t) is neglected, and the cross terms are also neglected as we assume that the signals are mutually cyclically uncorrelated. Now we look at the term r α To validate the approximation in (15), we need to assume As we know, for a given signal a(t), to make the approximation a(t ± τ d ) ≈ a(t) appropriate, the time difference τ d must be less than the coherence time τ c of the signal. We know τ c ∝ 1/B w , where B w is the bandwidth of the signal a(t), so in order to measure the accuracy of the approximation, we define ρ = τ d B w . The smaller the ρ is, the more accurate the approximation will be. Then, in improved cyclic MUSIC, . . , I, and p, n = 1, . . . , N.
As we know, ∆ i = d sin θ i /c, it is obvious that ρ i reaches its maximum when |p − n| = N − 1 and θ i = ±90 • . If d is half of the signal wavelength, then the maximum value of ρ i will be (N − 1)B w /(4 f 0 ). As the signal is narrowband, B w / f 0 will be small, the number of antennas N could not be too large, so we can say that ρ i is small, thus, the approximation in (16) could be made. Compared with (13), the maximum value of ρ i for the approximation in (16) is only half of that of the cyclic MUSIC in (13). Thus, the new approximation (16) is more accurate than (13). Now substituting (15) into (14), we obtain the (p, n)th element of R α xx (τ) as In order to write R α xx (τ) more compactly, we define where M is a diagonal matrix with its (i, i)th element to be r α aiai (τ)e j2π f0τ , a( f , θ) is the steering vector evaluated at frequency f and impinging direction θ, and A( f ) is constructed with the steering vectors for all I signals. Now putting all the elements of (17) together and using (18), (19), and (20), we obtain To see the difference between our expression of cyclic correlation matrix (21), and the one for cyclic MUSIC (7), we study the middle factor of (7), R α ss (τ). With the assumption that the sources are mutually cyclically uncorrelated, we could obtain that R α ss (τ) is actually equal to M by simple deduction. Thus, when the bandwidth of the signal becomes narrower, that is, α becomes smaller, (21) would approach (7). However, as long as the signal is not pure sinusoidal, α will not be 0, and (21) will be different from (7). The difference is a consequence of the factor e j2π(α/2)(p+n−2)∆i in (14), which is lost in cyclic MUSIC due to the approximation made in (13) before calculating the cyclic correlation.
If A( f 0 + α/2) and A( f 0 − α/2) are full rank, which is true with the assumption that the array manifold is unambiguous, the steering vector a( f 0 + α/2, θ i ) which is connected with the direction of the ith signal will be in the signal subspace or orthogonal to the noise subspace of R α xx (τ).
Then applying SVD to (21) to estimate the signal subspace and the noise subspace the same way as in cyclic MUSIC, we obtain (8). We then search θ and find the maxima of 1/ E H n a( f 0 + α/2, θ) 2 or E H s a( f 0 + α/2, θ) 2 in terms of θ as the DOA estimate. We notice that here, the steering vector is evaluated at f 0 + α/2, not at f 0 as in cyclic MUSIC. Alternatively, as the steering vector a( f 0 − α/2, θ i ) is lying in the subspace spanned by V s or orthogonal to the subspace spanned by V n , we can find the maxima of 1/ V H n a( f 0 − α/2, θ) 2 or V H s a( f 0 − α/2, θ) 2 in terms of θ as the DOA estimate. This alternative way provides DOA estimates that are almost identical as the first way ( f 0 + α/2). Therefore, this method will not be considered in the discussion hereafter. The improved cyclic MUSIC algorithm could be summarized as follows: (1) estimate R α xx (τ); (2) apply SVD to R α xx (τ) to estimate E s or E n ; (3) find the maxima of 1/ E H n a( f 0 +α/2, θ) 2 or E H s a( f 0 + α/2, θ) 2 in terms of θ as DOA estimates.
In cyclic MUSIC, on the other hand, x n (t) is approximated by thus, the cyclic correlation of the approximated x p (t) and x n (t) in cyclic MUSIC is r α xpxn (τ) ≈ e j2π f0(p−1)∆ r α aa (τ)e j2π f0τ e − j2π f0(n−1)∆ .
Comparing (23) and (25) with (22), we notice that both improved cyclic MUSIC and cyclic MUSIC make the approximation r α aa (τ + (p − n)∆) ≈ r α aa (τ), but cyclic MUSIC also approximates f 0 ± α/2 by f 0 . This is due to the loss of phase information by making approximation (24) before calculating the cyclic correlation. Therefore, an error of α/2 in the frequency for evaluating the steering vector has been resulted in cyclic MUSIC.
From the summary of our improved cyclic MUSIC algorithm, we notice that the first two steps are exactly the same as those of cyclic MUSIC, thus, the estimated signal subspace E s or noise subspace E n is the same for both algorithms. Therefore, if we define a steering vector a(ω) = [1, e jω , . . . , e j(N−1)ω ] T , where ω = 2π f d sin θ/c, then the estimated ω, by finding the maxima of 1/ E H n a(ω) 2 or E H s a(ω) 2 , is the same too. But since ω is a function of f and θ, the error in f will cause estimation error in θ. In our deduction above, we know that an error of α/2 in the frequency for evaluating the steering vector has resulted in cyclic MU-SIC, so we can conclude that the estimated DOA for cyclic MUSIC is biased.
Specifically, if we define γ = α/ f 0 , the actual DOA being θ 0 , the estimated DOAs for cyclic MUSIC and improved cyclic MUSIC being θ 1 0 and θ 2 0 , respectively, and the error of the estimated DOAs for cyclic MUSIC and improved cyclic MUSIC being δ 1 0 and δ 2 0 , respectively, then since estimated ω is equal, we obtain and then, If δ is small, we can write sin(θ 0 + δ) ≈ sin θ 0 + δ cos θ 0 . So simplifying (27) using this approximation, we obtain then the mean error and the standard deviation (STD) of the estimated DOAs for cyclic MUSIC and improved cyclic MU-SIC will be related as std where a denotes the mean of a, and std(a) denotes the STD of a. In improved cyclic MUSIC, we approximate the (p, n)th element of the cyclic correlation matrix R α xx (τ) by (23). Comparing (22) and (23), we see that this approximation is fairly accurate. Therefore, MUSIC should be able to estimate ω = 2π( f 0 + α/2)d sin θ/c with only a small error. Note that this statement may also be true for cyclic MUSIC which may also estimate the same ω accurately. However, in improved cyclic MUSIC, θ is obtained from ω = 2π( f 0 + α/2)d sin θ/c which is identical to the left exponent of (22). In other words, δ 2 0 should be very small. Thus from (29), δ 1 0 ≈ (γ/2) tan θ 0 , which depends on γ and θ 0 . Since α is chosen to be equal to the baud rate (or multiples of the baud rate), γ is related to the bandwidth of the signal. Although for narrowband signals, γ is very small, when the direction of the impinging signal θ 0 is large, (γ/2) tan θ 0 could be of several degrees, so the bias of DOA estimation for cyclic MUSIC is obvious. On the other hand as γ is small for narrowband signals, the difference of STD for cyclic MUSIC and improved cyclic MUSIC is not obvious. The simulation results shown in Section 5 will illustrate our conclusions here.

IMPROVED CONJUGATE CYCLIC MUSIC
In improved cyclic MUSIC, if we do not make the approximation as in (16), (21) should be rewritten as where C is the error term. Clearly, if C can be forced to zero, the performance of DOA estimation can be further improved. We now show how we achieve this objective in our improved conjugate cyclic MUSIC algorithm. First, we calculate the cyclic conjugate correlation of x(t) to obtain R α xx * (τ). Here, we evaluate the cyclic conjugate correlation at cycle frequency α and choose α = α + 2 f 0 , where α is the cycle frequency for the cyclic correlation and f 0 is the carrier frequency. Then, the (p, n)th element of R α xx * (τ) is where the shift property of cyclic conjugate correlation is applied, that is, if y(t) = x(t + T), then r α yy * (τ) = r α xx * (τ)e j2παT , and as in improved cyclic MUSIC, n n (t) and the cross terms are neglected due to cyclic uncorrelation between different sources. Now we define the derivative of r α aia * i (τ) as d i (τ), and approximate r α We also define A( f ) as in (20), and Then, R α xx * (τ) can be written as where the (p, n)th element of C i , C i (p, n), can be approximated by It is obvious that C i is odd symmetric, that is, C i (p, n) ≈ −C i (n, p). Now transposing R α xx * (τ), then averaging R α xx * (τ) and R α xx * (τ) T , we obtain Now applying SVD to (37), we can perform DOA estimation the same way as in improved cyclic MUSIC. The improved conjugate cyclic MUSIC algorithm could be summarized as follows: apply SVD to R to estimate E s or E n ; (4) find the maxima of 1/ E H n a( f 0 +α/2, θ) 2 or E H s a( f 0 + α/2, θ) 2 in terms of θ as DOA estimates.
Note that the steering vector is evaluated at the same frequency as in improved cyclic MUSIC, that is, f 0 + α/2, so the bias of the estimated DOA for improved conjugate cyclic MUSIC is also much reduced from cyclic MUSIC. Moreover, as the approximation in (37) is more accurate than that of (21), the performance of improved conjugate cyclic MUSIC is potentially better than that of improved cyclic MUSIC. The simulation results shown in the next section will illustrate our conclusions here.

SIMULATION RESULTS
Four sets of simulations are performed to compare the performance of cyclic MUSIC, improved cyclic MUSIC, and improved conjugate cyclic MUSIC.

Simulation 1
This simulation is tested on a narrowband BPSK signal with raised-cosine pulse shape of a roll-off factor 0.7. The carrier frequency of the signal f 0 is 20 MHz and the symbol rate f b , which is also the cycle frequency, is 0.5 MHz. The bandwidthto-carrier frequency ratio of this signal is around 5%. Signalto-noise ratio (SNR) is assumed to be 10 dB. Three antennas equally separated at a distance of half wavelength of the signal are assumed to receive this signal from a DOA of θ 0 = 60 • . A number of 3200 snapshots are used in the computations. We run the program fifty times and the resulted spatial spectra are shown in Figure 1. The mean of the estimated DOA is 61.26 • for cyclic MUSIC, 60.03 • for improved cyclic MUSIC, and 60.01 • for improved conjugate cyclic MUSIC. If we substitute θ 0 = 60 • , δ 2 = 0.03 • , and γ = α/ f 0 = 0.5 MHz/20 MHz = 0.025 into (29), we obtain δ 1 = 1.2709 • . We can see that the simulated mean error of the estimated DOA for cyclic MUSIC is almost the same as the theoretical value. We also notice that the peaks of the spatial spectra resulted from improved conjugate cyclic MUSIC is much sharper than that resulted from cyclic MU-SIC or improved cyclic MUSIC. Simulation 4 will further elaborate this.

Simulation 2
From the performance analysis and (29), we can see that the mean error of the estimated DOA for cyclic MUSIC varies with the impinging direction of the signal, yet by choosing the correct frequency for evaluating the steering vector, improved cyclic MUSIC and improved conjugate cyclic MUSIC should not have this problem. To see this effect, we vary θ 0 from −60 • to 60 • , run the program for 200 times, and obtain the result as shown in Figure 2. We notice that the mean error of the estimated DOA for cyclic MUSIC is proportional to tan θ 0 , yet the mean error for improved cyclic MUSIC or improved conjugate cyclic MUSIC are approximately zero with all values of θ 0 .

Simulation 3
From the performance analysis and (29), we see that the mean error of the estimated DOA for cyclic MUSIC also varies with cycle frequency α, yet by choosing the correct frequency for evaluating the steering vector, improved cyclic MUSIC and improved conjugate cyclic MUSIC should not have this problem. To see this effect, we vary f b from 0.1 MHz to 0.5 MHz. Here we also choose α to be f b . We know that the raised-cosine shaping pulse with a roll-off factor a for the BPSK signal with a symbol rate f b , has a bandwidth (1+a) f b . If we define the bandwidth-to-carrier frequency ratio as ξ, then with the changing of f b , ξ also changes, but the value of ξ will not exceed 5% in this simulation. So the BPSK signals with f b from 0.1 MHz to 0.5 MHz used here are still narrowband. We run the program 200 times and the results are shown in Figures 3 and 4. Figure 3 is the mean error of the estimated DOA versus γ defined in Section 3.2 and expressed in percentage here. We notice that the mean error of the estimated DOA for cyclic MUSIC is changing proportionally with γ, while the mean error for improved cyclic MUSIC or improved conjugate cyclic MUSIC is almost zero with all   values of γ. Figure 4 is the STD of the estimated DOA versus γ. We cannot see much difference between cyclic MUSIC and improved cyclic MUSIC, which is what we predicted in Section 3.2. But the STD of the estimated DOA for improved conjugate cyclic MUSIC is much smaller than that for cyclic MUSIC or improved cyclic MUSIC.

Simulation 4
This simulation is tested on two closely impinging narrowband BPSK signals with raised-cosine pulse shape of a roll-off factor 0.7. The carrier frequency of the signal f 0 is 20 MHz and the symbol rate f b , which is also the cycle frequency,  is 0.5 MHz for both signals. SNR is assumed to be 15 dB. Six antennas equally separated at a distance of half the wavelength of the signals are assumed to receive these two signals from DOAs of 35 • and 40 • . A number of 3200 snapshots are used in the computations. The result is shown in Figure 5. We notice that improved conjugate cyclic MUSIC separates these two DOAs successfully, but cyclic MUSIC or improved cyclic MUSIC fails to separate them.

CONCLUSION
This paper has presented an improved cyclic MUSIC algorithm by choosing a more appropriate frequency for evaluating the steering vector used in DOA estimation such that the bias of the estimated DOA has been reduced. The improved conjugate cyclic MUSIC further improves the performance of DOA estimation by approximately cancelling the error part in the cyclic conjugate correlation matrix. The effectiveness of the proposed methods have been demonstrated by the simulation results.