1 Introduction

Underwater acoustic communication is challenging due to the large complexity of the underwater acoustic channel and serious noise interference [1,2,3]. Spread spectrum communication technology has good anti-jamming performance and can ensure reliable communication in complex channels, so it is often used in underwater acoustic communication [4,5,6]. The underwater acoustic channel is a typical coherent multipath channel. Signals arriving along different paths have different propagation delays and directions of arrival (DOA). So, the received signal has spatiotemporal characteristics [7, 8]. The received signal has delay spread and angle spread. Coherent superposition of multipath signals can lead to severe inter-symbol interference (ISI) in the received signal [9, 10]. A space-time processor based on array reception can filter signals arriving along each path separately, and it can enhance the communication system’s reliability by considering the diversity of spatiotemporal clusters. The core problem of the space-time processor is DOA and source number estimation [11, 12]. The spread spectrum signal is wideband and usually works at a low signal-to-noise ratio (SNR). The conventional wideband signal DOA estimation algorithms and source estimation algorithms are usually developed based on narrowband signal estimation methods. Since the array is affected by frequency, the signal subspaces of different frequencies are different, and an array steering matrix cannot be constructed using a single frequency. As a result, many DOA and source number estimation methods of narrowband signals cannot be directly applied to wideband signal source number estimation [13,14,15].

For wideband signals, the common DOA estimation method usually uses a fast Fourier transform (FFT) to transform the received data into snapshot data of each subcarrier and then process narrowband signals at the sub-frequency points within the signal bandwidth [16,17,18,19,20]. However, this method’s estimation accuracy is determined by the division density of the signal frequency points and the structure of the focusing matrix, and the calculation amount is large.

At present, most of the research and experiments use uniform linear arrays as receivers. The receiving end and the sending end of the communication are far apart. In this case, the influence of the relative motion of the two communicating parties on the signal arrival angle can be ignored [21, 22]. If the distance between the two communicating parties is close and there is a fast relative motion, the signal arrival angle will change significantly in each communication. In this case, it is necessary to continuously perform DOA estimation in the communication and obtain the change of the signal’s DOA in real time for subsequent processing. If the conventional wideband signal DOA estimation method is applied, the computational complexity will be increased, and the real-time performance will be decreased.

Thus, to estimate the DOA of the signal quickly, the characteristics of anti-multipath, anti-interference, and good self-coherence of the spread spectrum sequence are used in the case of spread spectrum communication [23]. Based on this, a MUSIC algorithm is developed. It utilizes the good self-coherence of the baseband spread spectrum signal to accurately estimate the DOA under a low SNR. There is no need to transform the signal into the frequency domain sub-band, which effectively reduces the algorithm’s complexity. However, since the algorithm uses the received signal’s covariance matrix for DOA estimation, the number of sources needs to be predicted in advance [24,25,26,27].

In the case of low SNR, a commonly used method is to use the steering matrix estimated by the blind beamforming technology for estimating the covariance matrix [28,29,30]. Although this method can perform better signal source number estimation under a low SNR, it can only be applied to process narrowband signals. For wideband signals, the number of sources is estimated by adopting the Incoherent Signal Subspace Method (ISM). Then, the Eigenvalue Grads Method (EGM) based on the statistical characteristics of noise eigenvalues is proposed for signal source number estimation [31]. However, it only uses the sub-band information without comprehensive utilization of all bandwidth information, the estimation performance is limited, and the resolution is low. Furthermore, the Coherent Signal Subspace Method (CSM) and its improved method are developed, e.g., Two-side Correlation Transformation (TCT), Signal Subspace Transformation (SST), etc. However, in these methods, the array output needs to be focus-transformed first, and the estimation performance of the method is affected by the focusing matrix. In addition, when these methods are applied to broadband signal source number estimation, there are special requirements for the construction method of the focusing matrix. That is, most of them need to use the method based on angle pre-estimation to construct the focusing matrix, which generates a large computation overload. Also, the focus matrix involves a complicated transformation, which often introduces transformation errors. Besides, this type of algorithm can obtain good estimation results only under the condition of large snapshots [17, 18, 32, 33]. Therefore, this paper proposes a method that performs singular value decomposition (SVD) of the delay structure information of the array element to estimate the number of sources of the spread spectrum sequence under a low SNR.

The DOA estimation experiment of the spread spectrum sequence is conducted through simulation. Experiments show that our method can quickly track the change in the signal’s angle of arrival, and the method based on the SVD of the delay structure information of the array element performs well in signal source number estimation under a low SNR.

2 Method

2.1 Channel model

Generally, Underwater acoustic (UWA) channels are random channels that vary in time and space [7]. The received signal is a coherent superposition of the transmission signals by all the sound rays at the receiver. Due to the influence of different sound rays, the received signal exhibits spatiotemporal characteristics [34]. In general, this paper ignores the frequency characteristics of dielectric absorption and assumes that there is no dispersion along any path. Suppose there are P paths and the acoustic signal waveform along each propagation path remains unchanged. For the UWA channel, the impulse response can be expressed as:

$$\begin{aligned} \begin{aligned} h(t)=\sum \limits _{p=1}^{P}{{{A}_{p}}\delta (t-{{\tau }_{p}})} \end{aligned} \end{aligned}$$
(1)

where \({{A}_{p}}\) is the amplitude of the p-th path, and \({{\tau }_{p}}\) is the delay of the p-th path.

Another channel model takes the DOA of different paths at the receiving end as the classification standard. Then, for the coherent multipath channel, the impulse response is expressed as:

$$\begin{aligned} \begin{aligned} h(t)=\sum \limits _{j=1}^{J}{{{A}_{j}}\sum \limits _{l=1}^{{{L}_{\text {j}}}}{\delta (t-{{\tau }_{j,l}})}} \end{aligned} \end{aligned}$$
(2)

where J is the DOA number of all paths. \({{L}_{j}}\) is the number of paths in the j-th direction, \({{\tau }_{j,l}}\) is the delay of the corresponding path, and:

$$\begin{aligned} \begin{aligned} \sum _{j=1}^{L}L_{j}=P \end{aligned} \end{aligned}$$
(3)

The communication system is illustrated in Fig. 1. The receiver uses a uniform linear array (ULA) to receive UWA signals. Assuming that the ULA is located in the far-field, the UWA signal incident on the ULA is a parallel plane wave. In the subsequent analysis, the first channel model is adopted.

Fig. 1
figure 1

Schematic diagram of the multipath channel. In the communication system, the receiver uses a uniform linear array to receive signals

Full size image

2.2 Receive signal model

Assuming that the transmitted signal is a single-frequency signal spread by m-sequence, it can be represented as:

$$\begin{aligned} \begin{aligned} C\left( t \right) =c\left( t \right) \text {cos}\left( {{\omega }_{c}}t \right) \ \end{aligned} \end{aligned}$$
(4)

where c(t) is the i-th symbol of the m-sequence and \(c(t)=\sum \limits _{i=0}^{{{N}_{m}}-1}{{{c}_{i}}{{P}_{{{T}_{c}}}}(i-n{{T}_{c}})}\),\({{c}_{i}}\in (-1,1)\). \({{T}_{c}}\) is the symbol interval, \({{N}_{m}}\) is the period of the m-th sequence. \({{P}_{{{T}_{c}}}}\)is the symbol pulse shaping filter, and the root raised cosine filter is adopted. The roll-off coefficient is \(\beta\). Then, the bandwidth of the signal is \(B={(1+\beta )}/{{{T}_{c}}}\;\).

This paper takes the ULA as an example and assumes that the sound source is in the far-field. After passing through the channel, the spread spectrum signal on the \({{N}_{p}}\) path is incident on the M-element ULA [35], as demonstrated in Fig. 2.

Fig. 2
figure 2

A plane wave is an incident on a ULA. ULA: a uniform linear array

Full size image

The first array element is the reference array element. d represents the distance between the array elements. The signal received by the m-th array element is:

$$\begin{aligned} \begin{aligned} {{x}_{m}}(t)=\sum \limits _{p=1}^{{{N}_{p}}}{{{A}_{p}}C\left[ (1+a)t-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}) \right] }+{{n}_{m}}(t),\text {1}\le \text {m}\le \text {M}\ \end{aligned} \end{aligned}$$
(5)

where \({{A}_{p}}\) represents the gain of the p-th path. a is the Doppler compression factor. \({{\tau }_{p}}\) represents the time delay of the p-th path that reaches the reference array element. \({{\tau }_{m}}({{\theta }_{p}})\) represents the sound path difference of the m-th array element relative to the reference array element. \({{n}_{m}}(t)\) represents the noise of the m-th element.

2.3 MUSIC Algorithm Based on Spread Spectrum Sequence

The MUSIC algorithm has been widely used because of its high resolution and precision. However, it also has several shortcomings [36]. For instance, the MUSIC algorithm suffers from a reduced resolution under a low SNR and few snapshots, and the estimation performance drops substantially. Meanwhile, the algorithm fails in a multipath environment because the signals are highly correlated [37]. These shortcomings affect the application of this algorithm in the communication field. However, using the MUSIC algorithm to estimate the modulated m-sequence does not require decoherence such as spatial smoothing, and the signal’s DOA can be estimated in the easiest and fastest way.

Before using the MUSIC algorithm, quadrature demodulation must first be performed to create a complex signal. In Fig. 3, the signal received by the m-th element is multiplied by sine and cosine signals and then processed by a low-pass filter. This process is called quadrature demodulation [38].

Fig. 3
figure 3

Quadrature demodulation. LPF: Low-pass filter

Full size image

The signal after quadrature demodulation can be represented as:

$$\begin{aligned} \begin{aligned} y_{mi}(t)&=2x_{m}(t)cos(\omega t)\\&\overset{LPF}{\mathop {=}}\,\sum \limits _{p=1}^{{{N}_{p}}}{{{A}_{p}}c\left[ (1+a)t-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}) \right] }\\&\cdot \cos \left[ a{{\omega }_{\text {c}}}t-{{\omega }_{\text {c}}}{{\tau }_{p}}-{{\omega }_{\text {c}}}{{\tau }_{m}}({{\theta }_{p}}) \right] +{{n}_{mI}}(t)\ \end{aligned} \end{aligned}$$
(6)
$$\begin{aligned} \begin{aligned} {{y}_{mQ}}(t)&=-2{{x}_{m}}(t)\sin ({{\omega }_{c}}t)\\&\overset{LPF}{\mathop {=}}\,\sum \limits _{p=1}^{{{N}_{p}}}{{{A}_{p}}c\left[ (1+a)t-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}) \right] }\\&\cdot \sin \left[ a{{\omega }_{\text {c}}}t-{{\omega }_{\text {c}}}{{\tau }_{p}}-{{\omega }_{\text {c}}}{{\tau }_{m}}({{\theta }_{p}}) \right] +{{n}_{mQ}}(t)\ \end{aligned} \end{aligned}$$
(7)

where \({{n}_{mIQ}}(t)={{n}_{mI}}(t)+j{{n}_{mQ}}(t)\), and \(j=\sqrt{-1}\). The output of quadrature demodulation is a complex number consisting of in-direction and quadrature components:

$$\begin{aligned} \begin{aligned} {{y}_{m}}(t)&={{y}_{mI}}(t)+j{{y}_{mQ}}(t)\\&=\sum \limits _{p=1}^{{{N}_{p}}}{{{A}_{p}}c\left[ (1+a)t-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}) \right] }\\&\cdot {{e}^{j{{\omega }_{\text {c}}}\left[ at-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}) \right] }}+{{n}_{mIQ}}(t)\ \end{aligned} \end{aligned}$$
(8)

Here, an approximation is made:

$$\begin{aligned} \begin{aligned} c(t+at-{{\tau }_{p}}-{{\tau }_{m}}({{\theta }_{p}}))\approx c(t+at-{{\tau }_{p}}) \end{aligned} \end{aligned}$$
(9)

Experiments indicate that this approximation influences DOA estimation performance. Still, the orientation can be estimated. After approximation, for the p-th path, the baseband response is:

$$\begin{aligned} \begin{aligned} {{\mathbf {Y}}_{p}}={{A}_{p}}c(t+{{a}_{p}}t-{{\tau }_{p}}){{e}^{-j{{\omega }_{c}}({{\tau }_{p}}-at)}}\mathbf {a}({{\theta }_{p}})+{{\mathbf {n}}_{IQ}}(t) \end{aligned} \end{aligned}$$
(10)

In Eq. (10), we have:

$$\begin{aligned} \begin{aligned} \mathbf {a}({{\theta }_{p}})={{\left[ \begin{matrix} {{e}^{-j{{\omega }_{c}}{{\tau }_{1}}({{\theta }_{p}})}} &{} {{e}^{-j{{\omega }_{c}}{{\tau }_{2}}({{\theta }_{p}})}} &{} \cdots &{} {{e}^{-j{{\omega }_{c}}{{\tau }_{M}}({{\theta }_{p}})}} \\ \end{matrix} \right] }^{T}}\ \end{aligned} \end{aligned}$$
(11)

Then, we get baseband responses for all paths:

$$\begin{aligned} \begin{aligned} \mathbf {Y}=\sum \limits _{p=1}^{{{N}_{p}}}{{{\mathbf {Y}}_{p}}}=\mathbf {AS}+\mathbf {n}\ \end{aligned} \end{aligned}$$
(12)

In Eq. (12):

$$\begin{aligned} \begin{aligned} \mathbf {A}&=\left[ \begin{array}{llll} \mathbf {a}({{\theta }_{1}}) &{} \mathbf {a}({{\theta }_{2}}) &{} \cdots &{} \mathbf {a}({{\theta }_{{{N}_{p}}}}) \\ \end{array} \right] \ \end{aligned} \end{aligned}$$
(13)
$$\begin{aligned} \begin{aligned} \mathbf {S}&=\left[ \begin{matrix} {{A}_{1}}{{c}_{n}}(t+at-{{\tau }_{1}}){{e}^{-j{{\omega }_{c}}({{\tau }_{1}}-at)}} \\ \begin{matrix} {{A}_{2}}{{c}_{n}}(t+at-{{\tau }_{2}}){{e}^{-j{{\omega }_{c}}({{\tau }_{2}}-at)}} \\ \vdots \\ \end{matrix} \\ {{A}_{p}}{{c}_{n}}(t+at-{{\tau }_{{{N}_{p}}}}){{e}^{-j{{\omega }_{c}}({{\tau }_{{{N}_{p}}}}-at)}} \\ \end{matrix} \right] \ \end{aligned} \end{aligned}$$
(14)
$$\begin{aligned} \begin{aligned} \mathbf {n}&={{[\begin{matrix} {{n}_{IQ1}} &{} {{n}_{IQ2}} &{} \cdots &{} {{n}_{IQM}} \\ \end{matrix}]}^{T}}\ \end{aligned} \end{aligned}$$
(15)

The baseband signal’s covariance matrix is represented as:

$$\begin{aligned} \begin{aligned} {{\mathbf {R}}_{yy}}=\mathbf {Y}{{\mathbf {Y}}^{H}}=\mathbf {A}{{\mathbf {R}}_{\text {s}}}{{\mathbf {A}}^{H}}+{{\sigma }^{2}}\mathbf {I}\ \end{aligned} \end{aligned}$$
(16)

where H stands for the conjugate transpose; \({{\sigma }^{2}}\)stands for the noise power; \(\mathbf {I}\) stands for the identity matrix; \({{\mathbf {R}}_{\text {s}}}\) stands for the signal’s covariance matrix:

$$\begin{aligned} \begin{aligned} {{\mathbf {R}}_{s}}=\mathbf {S}{{\mathbf {S}}^{H}}\ \end{aligned} \end{aligned}$$
(17)

Since a uniform linear matrix is used,\(\mathbf {A}\) is a Vandermonde matrix. When the number of array elements is larger than that of signal sources and each path has different incident directions, the rank of \({{\mathbf {R}}_{yy}}\) is dependent on\({{\mathbf {R}}_{s}}\):

$$\begin{aligned} \begin{aligned} rank({{\mathbf {R}}_{yy}})=rank({{\mathbf {R}}_{s}})\ \end{aligned} \end{aligned}$$
(18)

When coherent signals are received by the array in different directions, a rank loss will be caused, and the signal eigenvectors will diverge to the noise subspace. In this case, the MUSIC algorithm is ineffective. Since the m-sequence has excellent autocorrelation characteristics, a rank-free matrix can be obtained. Then, the MUSIC algorithm can be applied as usual. Equation (19) indicates that each element in matrix \({{\mathbf {R}}_{s}}\) involves the spreading sequence correlation function \(R(\tau )\), and \(*\) stands for the conjugate operation.

$$\begin{aligned} \begin{aligned} {{\mathbf {R}}_{s}}&=E\left\{ \left[ \begin{matrix} {{A}_{1}}{{c}_{n}}(t+at-{{\tau }_{1}}){{e}^{-j{{\omega }_{c}}({{\tau }_{1}}-at)}} \\ {{A}_{2}}{{c}_{n}}(t+at-{{\tau }_{2}}){{e}^{-j{{\omega }_{c}}({{\tau }_{2}}-at)}} \\ \vdots \\ {{A}_{p}}{{c}_{n}}(t+at-{{\tau }_{{{N}_{p}}}}){{e}^{-j{{\omega }_{c}}({{\tau }_{{{N}_{p}}}}-at)}} \\ \end{matrix} \right] \right. \ \\&\cdot \left. \left[ \begin{matrix} {{A}_{1}}c_{n}^{*}(t+at-{{\tau }_{1}}){{e}^{j{{\omega }_{c}}({{\tau }_{1}}-at)}} &{} \cdots &{} {{A}_{p}}c_{n}^{*}(t+at-{{\tau }_{{{N}_{p}}}}){{e}^{j{{\omega }_{c}}({{\tau }_{{{N}_{p}}}}-at)}} \\ \end{matrix} \right] \right\} \ \\&=E\left\{ \left[ \begin{matrix} {{A}_{1}}^{2}R(0) &{} \cdots &{} {{A}_{1}}{{A}_{p}}R({{\tau }_{1}}-{{\tau }_{{{N}_{p}}}}){{e}^{-j{{\omega }_{c}}({{\tau }_{1}}-{{\tau }_{{{N}_{p}}}})}} \\ \vdots &{} \ddots &{} \vdots \\ {{A}_{p}}{{A}_{1}}R({{\tau }_{{{N}_{p}}}}-{{\tau }_{1}}){{e}^{j{{\omega }_{c}}({{\tau }_{1}}-{{\tau }_{{{N}_{p}}}})}} &{} \cdots &{} A_{p}^{2}R(0) \\ \end{matrix} \right] \right\} \ \end{aligned} \end{aligned}$$
(19)

The speed of the ship is usually within 10 m/s. The speed of sound in water is 1500 m/s. So the Doppler factor ranges from \(-\)0.014 to 0.014. In Fig. 4, for a large Doppler compression factor, there is still good autocorrelation in the spread spectrum sequence [39]. So, the MUSIC algorithm is not affected by Doppler. Meanwhile, the correlation output value decreases quickly with the increase of \(|\tau |\). Ideally, it will be approximated as a diagonal matrix. Based on this, the rank deficient problem can be solved, and the MUSIC algorithm can efficiently estimate DOA.

Fig. 4
figure 4

Correlation output of m-sequence at different doppler factors. The a is doppler factor (x axis: The sampling point of the signal in the time domain; y axis: Correlative output of the signal)

Full size image

2.4 Estimated number of sources

Based on the spread spectrum sequence’s coherence, the MUSIC algorithm can quickly estimate the signal’s DOA. However, because the spread spectrum sequence usually works at a low SNR, conventional source estimation methods cannot achieve accurate source number estimation. Aiming at this problem, this paper proposes to perform source estimation based on the Hankel matrix’s SVD of element delay structure information. The general form of the Hankel matrix [40] is as follows:

$$\begin{aligned} \begin{aligned} {{H}_{n}}=[{{h}_{i+j-1}}]_{i,j=1}^{n}=\left[ \begin{matrix} {{h}_{1}} &{} {{h}_{2}} &{} \cdots &{} {{h}_{n}} \\ {{h}_{2}} &{} {{h}_{3}} &{} \cdots &{} {{h}_{n+1}} \\ \vdots &{} \vdots &{} \cdots &{} \vdots \\ {{h}_{n}} &{} {{h}_{n+1}} &{} \cdots &{} {{h}_{2n+1}} \\ \end{matrix} \right] \ \end{aligned} \end{aligned}$$
(20)

where \([{{h}_{i+j-1}}]_{i,j=1}^{n}\in C\).

By using the ULA time delay fixed value feature and the cross-correlation function of the specific array element signal (reference signal) and other array element signals, the Hankel matrix is reconstructed through the delay difference operation. The reconstructed Hankel matrix is in the following form:

$$\begin{aligned} \begin{aligned} {{H}_{i+k+q,j}}={{\left[ \begin{matrix} {{h}_{i,j}} &{} {{h}_{i+1,j}} &{} \cdots &{} {{h}_{i+q-1,j}} \\ {{h}_{i+1,j}} &{} {{h}_{i+2,j}} &{} \cdots &{} {{h}_{i+q,j}} \\ \vdots &{} \vdots &{} \cdots &{} \vdots \\ {{h}_{i+k-1,j}} &{} {{h}_{i+k,j}} &{} \cdots &{} {{h}_{i+(k-1)+(q-1),j}} \\ \end{matrix} \right] }_{k\times n}}\ \end{aligned} \end{aligned}$$
(21)

Among them, i, k, and q, respectively, represent the starting sequence number of matrix elements, the number of matrix rows, and the number of matrix columns. The subscript j of \({{H}_{i+k+q,j}}\) stands for the serial number of the reference signal, and the expression of \({{h}_{i,j}}\) is:

$$\begin{aligned} \begin{aligned} {{h}_{i,j}}=E\{{{x}_{i}}(t)x_{j}^{H}(t)\}-{{\sigma }^{2}}{{\delta }_{ij}},j=1,2,\cdots ,M \end{aligned} \end{aligned}$$
(22)

According to the above formula, when \(i=j\), the element \({{h}_{i,j}}\) contains the noise component; otherwise, there is only the signal component in the element \({{h}_{i,j}}\) because the signal and the noise are not related. Therefore, the Hankel matrix can be constructed as a matrix with only signal components. Through SVD of the Hankel matrix, we have:

$$\begin{aligned} \begin{aligned} {{H}_{i+k+q,j}}={{U}_{H}}{{\sum }_{H}}V_{H}^{\text {H}}\ \end{aligned} \end{aligned}$$
(23)

where \({{U}_{H}}=[{{u}_{H1}},{{u}_{H2}},\cdots ,{{u}_{Hk}}]\)is a \(k\times k\) left singular vector matrix. \({{V}_{H}}=[{{u}_{H1}},{{u}_{H2}},\cdots ,{{u}_{Hn}}]\) is a \(q\times q\) right singular vector matrix. \({{\sum }_{H}}\)is a \(k\times q\) singular value matrix, and it satisfies:

$$\begin{aligned} \begin{aligned} {{\sum }_{H}}=diag(\underbrace{{{\lambda }_{1}},{{\lambda }_{2}},\cdots ,{{\lambda }_{{N}_{p}}}}_{{N}_{p}}\underbrace{0,\cdots ,0}_{\min (k-{{N}_{p}},q-{{N}_{p}})}) \end{aligned} \end{aligned}$$
(24)

where the first \({{N}_{p}}\) singular values correspond to the signal, and the last \(\min (k-{{N}_{p}},q-{{N}_{p}})\) singular values correspond to the noise.

According to the analysis of the Hankel matrix, we can get that the Hankel matrix has the following characteristics:

  1. (1)

    It is different from having only one spatial covariance matrix of the received data and one eigenvalue decomposition operation for source number estimation. For enough array elements M, a large number of matrices of different reference signals and different dimensions \({{H}_{i+k+q,j}}\) can be constructed, which improves the variety of judgment methods and greatly enhances data utilization. Based on this, the ability to estimate the signal source number accurately is increased.

  2. (2)

    In the Hankel matrix \({{H}_{i+k+q,j}}\), the number of nonzero singular values equals the matrix’s rank and the signal source number \({{N}_{p}}\). In other words, the signal subspace covered by the left singular vector corresponding to the nonzero singular value equals that covered by the steering vector, and it is orthogonal to the noise subspace.

  3. (3)

    In practical applications, the noise variance \({{\sigma }^{2}}\) is unknown. Due to the flexibility of the construction method of the Hankel matrix, the Hankel matrix can be constructed by selecting elements \({{h}_{i,j}}\) that do not contain \(i=j\). In this approach, the influence of noise is avoided.

According to the singular values obtained after the decomposition of the Hankel matrix, the signal source number is estimated by heuristic decision criteria. That is, when the first C(i) satisfies \(C(i)\le 0\), there are \({\widehat{q}}=i-1\) signal sources. C(i) is defined as:

$$\begin{aligned} \begin{aligned} C(i)=\widehat{{{\lambda }_{i}}}-\frac{1}{\min (k,q)}\sum \limits _{i=1}^{\min (k,q)}{\widehat{{{\lambda }_{i}}}} \end{aligned} \end{aligned}$$
(25)

The algorithm is described as follows:

Step 1 Perform quadrature demodulation on the signals received by the M array elements to obtain a complex baseband signal.

Step 2 Construct a Hankel matrix according to the M complex baseband signals, and different Hankel matrices can be constructed by selecting different values.

Step 3 Conduct SVD on the obtained Hankel matrix in step 2 to obtain singular values.

Step 4 Determine the signal source number \({{N}_{p}}\) through heuristic criteria.

Step 5: Obtain the array’s covariance matrix based on the M complex baseband signals following Eq. (16).

Step 6: Conduct eigenvalue decomposition on the covariance matrix obtained according to Eq. (21).

$$\begin{aligned} \begin{aligned} {{\mathbf {R}}_{yy}}=\mathbf {U\Sigma }{{\mathbf {U}}^{H}} \end{aligned} \end{aligned}$$
(26)

where \(\mathbf {\Sigma }=diag\left\{ {{\lambda }_{1}},{{\lambda }_{2}},\cdots ,{{\lambda }_{M}} \right\}\) stands for the eigenvalue matrix, and \(\mathbf {U}\) stands for the eigenvector matrix.

Step 7: Sort the eigenvalues from large to small. Use the signal source number obtained in step 4 and form the signal subspace \({{\mathbf {U}}_{s}}\) by taking the eigenvectors corresponding to the first N eigenvalues. Then, form the noise subspace \({{\mathbf {U}}_{n}}\) by using the eigenvectors corresponding to the remaining smaller eigenvalues.

Step 8: Calculate the spatial spectrum function according to Eq. (27). Obtain the DOA estimate value of the signal by finding the peak value.

$$\begin{aligned} \begin{aligned} {{P}_{MUSIC}}=\frac{1}{{{a}^{H}}(\theta ){{\mathbf {U}}_{N}}\mathbf {U}_{N}^{H}a(\theta )} \end{aligned} \end{aligned}$$
(27)

3 Results and discussion

This section validates the effectiveness of our method and presents its performance. In the simulation experiments, the receiving array is composed of a linear array of 10 sensors. The distance between the sensors is half of the wavelength. There are two sources, which are incidents on the ULA from 10 and 20. The test SNR is − 30  dB to 10 dB, and the step size is 1 dB. In this paper, the performance of the proposed source estimation algorithm and the CSM method is compared under different SNRs. The parameter settings of the spread spectrum communication system are shown in Table. 1

Table 1 The parameter settings of the communication
Full size table

First, set \({{x}_{1}}\) as the reference signal selected in the simulation (The signal received by the first array element), and the Hankel matrix used in the simulation is:

$$\begin{aligned} H={{\left[ \begin{matrix} E\left\{ {{x}_{2}}(t)x_{1}^{H}(t) \right\} &{} E\left\{ {{x}_{3}}(t)x_{1}^{H}(t) \right\} &{} \cdots &{} E\left\{ {{x}_{6}}(t)x_{1}^{H}(t) \right\} \\ E\left\{ {{x}_{3}}(t)x_{1}^{H}(t) \right\} &{} E\left\{ {{x}_{4}}(t)x_{1}^{H}(t) \right\} &{} \cdots &{} E\left\{ {{x}_{7}}(t)x_{1}^{H}(t) \right\} \\ \vdots &{} \vdots &{} \cdots &{} \vdots \\ E\left\{ {{x}_{6}}(t)x_{1}^{H}(t) \right\} &{} E\left\{ {{x}_{7}}(t)x_{1}^{H}(t) \right\} &{} \cdots &{} E\left\{ {{x}_{10}}(t)x_{1}^{H}(t) \right\} \\ \end{matrix} \right] }_{5\times 5}}\ \end{aligned}$$
(28)
Fig. 5
figure 5

The detection probability of the proposed method and the CSM method under different SNRs. The figure compares the performance of the proposed source estimation algorithm and the CSM method (x axis: Signal-to-noise ratio form \(-30\) dB ~10 dB; y axis: Probability of successful detection). CSM: coherent signal subspace method. Proposed: MUSIC algorithm based on spread spectrum sequence

Full size image

Figure 5 indicates that the performance of the CSM method drops rapidly when the SNR is lower than 5 dB; our method achieves a high detection performance when the SNR is above − 20 dB, and a detection probability of more than 50% is guaranteed when the SNR is below − 20 dB.

Next, signals received by different arrays are selected as reference signals. Compare the performance of the proposed algorithm when the reference signals are different. The simulation results are shown in Fig. 6. As shown in Fig. 6, the choice of reference signal has no significant effect on performance.

Fig. 6
figure 6

Detection success rate of different reference signals (x axis: Signal-to-noise ratio form − 30 dB ~10 dB; y axis: Probability of successful detection)

Full size image

Next, the detection performance of our method is investigated under a different number of snapshots. Because our method is based on spread spectrum sequences, the detection ability of m-sequences of different orders under the same SNR is tested. The simulation results are presented in Fig. 7.

Fig. 7
figure 7

The detection probability of the number of sources under different snapshot (x axis: The Order of Spreading sequence; y axis: Probability of successful detection)

Full size image

As shown in Fig. 7, attributed to the excellent autocorrelation property of the spread spectrum sequence, the m-sequences of different orders have almost the same probability of successful detection under the same SNR. Meanwhile, the detection performance does not decrease with the number of data snapshots and declines.

Next, the probability that our method overestimates and underestimates the number of sources is investigated under a low SNR. The probability of overestimating (underestimating) is that the estimated source number is larger (smaller) than the actual one.

Fig. 8
figure 8

Overestimation (underestimation) probability at different SNR (x axis: Signal-to-noise ratio form − 30 dB ~10 dB; y axis: Probability of Overestimation (underestimation)). Overestimation: The estimated number of sources is more than the actual number of sources. Underestimation: The estimated number of sources is less than the actual number of sources

Full size image

According to Fig. 8, under an extremely low SNR, there is a high probability that the signal source number will be overestimated. This is because at an extremely low SNR after the Hankel matrix is processed by SVD, there is not a large difference between the singular values corresponding to the signal and the noise. Using heuristic criteria can lead to overestimation. With the gradual increase in the SNR, due to the selection of elements that do not contain \(i=j\) in the construction of the Hankel matrix, the singular value regarding the noise will decrease rapidly. However, the singular value regarding the signal will not rapidly increase, thus resulting in underestimation. In the communication system, overestimation for communications has the same effect of adding noises to the signal. Underestimation corresponds to a portion of the weak signal energy has not been utilized. Therefore, overestimation will lead to a decrease in the performance of the system, while underestimation has less impact on the performance of the system.

Next, the DOA estimation performance is analyzed by using spreading sequences. First, based on the spreading sequence the MUSIC algorithm’s estimation performance is analyzed under a low SNR. Assuming that there is only one signal source. The incident angle is -10, and the SNR is -25 ~10 dB. Assuming two types of signals, the spread spectrum signal and the single-frequency signal based on m-sequence, both are under the same SNR. The length of the m sequence is 63. The bandwidth of the spread spectrum signal is 2.5 kHz. The distance between the array elements is 18 cm, and the sound speed is 1500 m/s. The comparison of the average error of the two types of signals is shown in Fig. 9. The spatial spectral scan interval is 0.5.It can be found that although the spread spectrum signal is wideband, the estimated effect is only slightly worse than that of the single-frequency signal, which is completely within the allowable range of the project.

Fig. 9
figure 9

Comparison of average errors of different signal MUSIC algorithms with different SNR (x axis: Signal-to-noise ratio form − 25 dB ~10 dB; y axis: mean absolute error). MAE: mean absolute error

Full size image

Assuming that the signal is incident from -10, 40, and 10. The SNR is 0 dB and − 15 dB, and other conditions remain unchanged. Also, the MUSIC algorithm does not use decoherence algorithms such as spatial smoothing. The method proposed in this paper is compared with the improved ISM method based on spatial smoothing algorithm. As shown in Fig. 10, the single-frequency signal used MUSIC algorithms can no longer be distinguished. The improved ISM algorithm can distinguish the DOA of the signal in the case of 0 dB. When the SNR is reduced to − 15 dB, the improved ISM algorithm can not distinguish. Meanwhile, because of the excellent correlation of m-sequence, the MUSIC algorithm based on m-sequence can still be well distinguished.

Fig. 10
figure 10

The proposed method compared with the improved ISM method based on spatial smoothing algorithm. a is the comparison of the two methods under the 0 dB SNR. b is the comparison of the two methods under the − 15 dB SNR

Full size image

Through simulation experiments, it can be verified that the MUSIC algorithm based on the spread spectrum sequence can perform DOA estimation of the signal under a low SNR. Meanwhile, based on the good autocorrelation of the spread spectrum sequence, the algorithm can achieve automatic decoherence without using a smoothing algorithm and converting the signal to a sub-band for processing. In view of the problem that the MUSIC algorithm cannot estimate the number of sources, the Hankel matrix is constructed by using the spread spectrum sequence baseband complex signal. Besides, through SVD, the signal source number can be more accurately estimated under a low SNR. Moreover, in the situation where the two communication parties are relatively close and there is mutual movement, the signal source number and the DOA of the signal can be estimated in real time and quickly.

4 Conclusion

This paper proposes a spread spectrum sequence-based method for DOA estimation under a low SNR. For two communication parties that are not far apart, the relative motion greatly affects the arrival angle of the signal, and real-time DOA estimation is required. The conventional wideband coherent signal DOA estimation algorithm suffers from a large amount of computation. Based on the good correlation of the spread spectrum sequence and by applying the MUSIC algorithm to DOA estimation of the spread spectrum signal, our method can well estimate the DOA of the signal under a low SNR. Meanwhile, there is no need to convert the signal to the individual sub-bands, which effectively reduces the calculation overhead. Under a low SNR, the signal source number cannot be accurately estimated for the MUSIC algorithm. To address this issue, the Hankel matrix is constructed by using the spreading sequence, and the signal source number is obtained by the singular value. Besides, constructing the Hankel matrix by delay difference will reduce the matrix dimension compared with the signal’s covariance matrix. This helps to estimate the signal source number under a low SNR by sacrificing the maximum number of sources that the system can resolve.