Spectrum Sensing Using Co-Prime Array Based Modulated Wideband Converter

Abstract

As known to us all, it is challenging to monitor wideband signals in frequency domain due to the restriction of hardware. Several practical sampling schemes, such as multicoset sampling and the modulated wideband converter (MWC), have been proposed. In this work, a co-prime array (CA) based modulated wideband converter (MWC) spectrum sensing method is suggested. Our proposed method has the same sampling principle as the MWC but has some advantages compared to MWC. Firstly, CA-based MWC is an array-based MWC system. Each sensor is usually corrupted by independent noise for an array system which can be used for noise averaging, while all channels in conventional MWC have the same receiving noise. Secondly, by incorporating the co-prime array, we can estimate the power spectrum of signal directly employing its second-order statistical properties. Moreover, the system minimal sampling rate can be reduced further because of the reduction of sampling channels. Simulation results show that our method has better performance than traditional methods.

1. Introduction

Nowadays, spectral resources traditionally allocated to licensed users by governmental organizations are becoming scant. Cognitive Radio (CR) is an emerging dynamic spectrum management technology which can make the best use of spectral resources in wireless communications [1,2]. How to monitor the spectrum reliably and fast is an essential goal for many scholars [3]. A number of spectrum sensing methods exist, such as filterbank spectrum sensing and multi-taper spectrum sensing, etc. [4]. However, all the above methods are based on Nyquist sampling. In cognitive radio, it is necessary to sense a wide band of spectrum, leading to prohibitively high Nyquist rates. So this will affect the cost and power assumption of high-speed analog-to-digital converters (ADCs). Therefore, it is desirable to design a sub-Nyquist sampling scheme which can effectively estimate the spectrum without loss of any information.

Recently, several sub-Nyquist sampling scheme for spectrum sensing are developed [5,6,7,8]. In [5], the multi-coset sampling (MS) was used, where m low-rate cosets were chosen out of L cosets of samples. These samples are obtained from time uniformly distributed samples taken at a same sampling in each channel. And there is a different time offset in the sampling origin of each channel. In [6], a synchronous multirate sampling (SMRS) scheme was developed. Unlike the MS method, samples in SMRS are obtained at m different sampling rates and the sampling of all channels starts simultaneously at a given time. The above multicoset or time-interleaved approach suffers from some practical issues. Firstly, the radio-frequency (RF) signal is sampled directly. As a result, the signal bandwidth can exceed the analog bandwidth of ADCs. Another practical issue stems from the time shift since it is difficult to maintain accurate delays or synchronization among the ADCs at such high rates. To solve this problem, an analog system, referred to as the modulated wideband converter (MWC) which is comprised of a bank of modulators and low-pass filters is adopted in [7,8,9]. However, conventional MWC systems mainly have two disadvantages. One is that each channel has the same noise, which can’t be used for noise averaging. The other is that we need to choose appropriate periodic functions in different channels. These problems are overcome by adopting a uniform linear array (ULA) based MWC system [10]. In [10], a ULA is used for spectrum sensing where a same periodic function is chosen in each channel. Nevertheless, when signals have high carrier frequencies, there may exist mutual coupling between sensor elements. Moreover, the number of transmissions detected can be no more than the number of physical elements.

In this paper, we propose a spectrum sensing method using an improved MWC system based on co-prime array. The radio frequency (RF) signal impinges on a co-prime array. The received signal in each sensor is multiplied by a same periodic function, low-pass filtered and sampled at a low rate. In the reconstruction process, an enhanced virtual ULA can be produced by vectorizing the data covariance matrix of the co-prime array. It can detect more transmissions than ULA-based MWC or can reduce the system sampling rate further when the number of transmissions is fixed. Besides, each channel in our method is corrupted by different noise, so it has an advantage that the noise can be averaged which will improve sensing performance.

Notations. ⊗ and ⊙ denote the Kronecker product and Khatri-rao product, respectively. vec(.) the vectorization operator that turns a matrix into a vector by stacking all columns on top of the another. ∘ denotes the element-wise product. The complex conjugate operator is ${(.)}^{*}$, the transpose operator is ${(.)}^T$, and the complex conjugate-transpose is ${(.)}^H$. $\angle (.)$ returns the phase of input. $\mathrm{eig}(.)$ returns the eigenvalue of input matrix. $\lceil .\rceil$ returns the nearest integer towards positive infinity.

2. Array Signal Model

Consider a co-prime linear array consisting of two uniform linear subarrays with $2M$ and N sensor elements where M and N are co-prime integers and $M<N$. One has the position set $S_1=\{Mnd,0\le n\le N-1\}$ and the other has position set $S_2=\{Nmd,0\le m\le 2M-1\}$ where $d={\lambda }_{\min }/2$. ${\lambda }_{\min }$ denotes the wavelength corresponding to the highest frequency of interest. Since these two subarrays share the first sensor, namely the reference sensor, there are $N+2M-1$ sensors totally. Such array configuration allows most adjacent elements to be spaced farther apart, which is attractive when it is necessary to reduce mutual coupling between sensor elements. Following each sensor includes an analog front-end composed of a mixer with the same periodic function $p\left(t\right)$, a low-pass filter and a sampler at rate $f_s$. The sensing system is depicted in Figure 1 and Figure 2. In practice, to ensure that all the mixing functions are the same without any distortions. We can first produce the satisfying analog mixing signal $p\left(t\right)$, then let it pass through a power divider to obtain several same versions of $p\left(t\right)$. Each version of the output of power divider can be used as the mixing function of each channel.

Assume that there are K uncorrelated transmissions ${\left\{s_i\left(t\right)\right\}}_{i=1}^K$ impinging on the array which have identical and known direction of arrival (DOA). Each transmission $s_i\left(t\right)$ with the bandwidth $B_i$ is modulated by a carrier frequency $f_i\in \mathbb{R}$. So the received RF signal denoted as $z\left(t\right)$ which is bandlimited to $\mathcal{F}=\left[-\frac{-f_{Nyq}}{2}\frac{f_{Nyq}}{2}\right]$ can be written as

$$\begin{array}{c}\hfill z\left(t\right)=\sum _{i=1}^K{s}_i\left(t\right)e^{j2\pi f_{i}t}\end{array}$$

(1)

We define $B_{\max }\triangleq {max}_{i=1}^K{B}_i$ and assume ${min}_{i\ne j}|f_i-f_j|>B_{\max }$. Then the Fourier transform of $z\left(t\right)$ can be written as

$$\begin{array}{c}\hfill Z\left(f\right)={\int }_{-\infty }^{\infty }z\left(t\right)e^{-j2\pi ft}\mathrm{d}t=\sum _{i=1}^K{S}_i(f-f_i)\end{array}$$

(2)

where $S_i\left(f\right)$ is the Fourier transform of $s_i\left(t\right)$.

3. MWC Based on Co-Prime Array

The co-prime array model is depicted in Section 2. The received signal in each sensor is multiplied by a periodic function $p\left(t\right)$ whose period is defined as $T_p=1/f_p$, low-pass filtered with a filter that has cut-off frequency $f_s/2$ and sampled at the low rate $f_s$. $T_s=1/f_s$ is the sampling interval. We define ${\mathcal{F}}_s\triangleq [-f_s/2,f_s/2]$ and ${\mathcal{F}}_p\triangleq [-f_p/2,f_p/2]$. Referring to [7], we know that the sampling rate of each channel $f_s$ must satisfy $f_s\ge f_p\ge B_{\max }$, through which all signal information can be reserved in the baseband ${\mathcal{F}}_s$ without any loss. Here, for simplicity of analysis, we choose $f_s=f_p$.

Consider the received signal $z_n\left(t\right)$ at the n-th sensor of the co-prime array

$$\begin{array}{cc}\hfill z_n\left(t\right)& =\sum _{i=1}^K{s}_i(t+{\tau }_n)e^{j2\pi f_i(t+{\tau }_n)}\hfill \\ & \approx \sum _{i=1}^K{s}_i\left(t\right)e^{j2\pi f_i(t+{\tau }_n)}\hfill \end{array}$$

(3)

where ${\tau }_n=d_{n}cos\left(\theta \right)/c$ is the time delay for signal arriving at the n-th sensor with respect to the reference sensor. $\theta$ is measured from the axis which is parallel to the linear array. We will show that $\theta$ should not be equal to ${90}^{\circ }$ because different time delays between different sensors are useful for the estimation of carrier frequencies. $d_n\in \{S_1\cup S_2\}$ is the position of the n-the sensor. c is the speed of wave propagation. The approximation in (3) stems from the narrowband assumption on $s_i\left(t\right)e^{j2\pi f_{i}t}$. We denote the Fourier transform of $z_n\left(t\right)$ as

$$\begin{array}{c}\hfill Z_n\left(f\right)=\sum _{i=1}^K{S}_i(f-f_i)e^{j2\pi f_i{\tau }_n}\end{array}$$

(4)

A typical $Z_n\left(f\right)$ is shown in Figure 3a. We only show the amplitude spectrum in the figure.

In each sensor channel, it works like one channel of traditional MWC system [7]. So the discrete-time Fourier transform (DTFT) of output $y_n\left[k\right]$ in the n-th channel can be written as

$$\begin{array}{c}\hfill Y_n\left(e^{j2\pi f{T}_s}\right)=\sum _{l=-L_0}^{L_0}{c}_l{Z}_n(f-l{f}_p),f\in {\mathcal{F}}_s\end{array}$$

(5)

where

$$\begin{array}{c}\hfill c_l=\frac{1}{T_p}{\int }_0^{T_p}p\left(t\right)e^{-j2\pi lt/T_p}\mathrm{d}t\end{array}$$

(6)

and $L_0$ is chosen as the smallest integer so that the sum contains all nonzero contributions of $Z_n\left(f\right)$ over ${\mathcal{F}}_s$. Here, we can choose $L_0=\lceil \frac{f_{\mathrm{Nyq}}}{2f_p}\rceil$. For the sake of concreteness, the mixing function $p\left(t\right)$ is chosen as a piecewise constant function that alternates between the levels $\pm 1$ randomly for each of P equal time intervals. Formally,

$$\begin{array}{c}\hfill p\left(t\right)={\alpha }_k,\frac{k{T}_p}{P}\le t\le \frac{(k+1)T_p}{P}\end{array}$$

(7)

where $0\le k\le P-1$. ${\alpha }_k\in \{+1,-1\}$. $p\left(t\right)=p(t+n{T}_p),n\in \mathbb{Z}$. Our only request for $p\left(t\right)$ is that its Fourier coefficients $c_l\ne 0,-L_0\le l\le L_0$.

Substituting (4) into (5), we have

$$\begin{array}{cc}\hfill Y_n\left(e^{j2\pi f{T}_s}\right)& =\sum _{l=-L_0}^{L_0}{c}_l\sum _{i=1}^K{S}_i(f-f_i-l{f}_p)e^{j2\pi f_i{\tau }_n}\hfill \\ & =\sum _{i=1}^K{\tilde{S}}_i\left(f\right)e^{j2\pi f_i{\tau }_n},f\in {\mathcal{F}}_s\hfill \end{array}$$

(8)

where ${\tilde{S}}_i\left(f\right)\triangleq {\sum }_{l=-L_0}^{L_0}{c}_l{S}_i(f-f_i-l{f}_p)$ is a cyclic shifted and scaled version of $S_i\left(f\right)$ in the interval ${\mathcal{F}}_s$. The whole processing flow in frequency domain is clearly shown in Figure 3.

We write (8) in matrix form as

$$\begin{array}{c}\hfill \mathbf{Y}\left(e^{j2\pi f{T}_s}\right)=\mathbf{A}\tilde{\mathbf{S}}\left(f\right),f\in {\mathcal{F}}_s\end{array}$$

(9)

where $\mathbf{Y}\left(e^{j2\pi f{T}_s}\right)={[Y_1\left(e^{j2\pi f{T}_s}\right),\dots ,Y_{N+2M-1}\left(e^{j2\pi f{T}_s}\right)]}^T$ and $\tilde{\mathbf{S}}\left(f\right)={[{\tilde{S}}_1\left(f\right),\dots ,{\tilde{S}}_K\left(f\right)]}^T$. The matrix $\mathbf{A}$ is defined by

$$\begin{array}{c}\hfill \mathbf{A}=\left(\begin{array}{ccc}{e}^{j2\pi f_1{\tau }_1}& \cdots & e^{j2\pi f_K{\tau }_1}\\ ⋮& \ddots & ⋮\\ e^{j2\pi f_1{\tau }_{N+2M-1}}& \cdots & e^{j2\pi f_K{\tau }_{N+2M-1}}\end{array}\right)\end{array}$$

(10)

From (9), we can easily get the discretized model under sampling rate $f_s$ in the time domain

$$\begin{array}{c}\hfill \mathbf{y}\left[k\right]=\mathbf{A}\tilde{\mathbf{s}}\left[k\right]\end{array}$$

(11)

4. Reconstruction Method

4.1. Carrier Frequency Recovery

By invoking (11), the temporal covariance matrix can be written as

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{y}}& =& \mathrm{E}\left\{\mathbf{y}\left[k\right]{\mathbf{y}}^H\left[k\right]\right\}\hfill \\ & =& \mathrm{E}\left\{\mathbf{A}\tilde{\mathbf{s}}\left[k\right]{\tilde{\mathbf{s}}}^H\left[k\right]{\mathbf{A}}^H\right\}\hfill \\ & =& \mathbf{A}\mathrm{E}\left\{\tilde{\mathbf{s}}\left[k\right]{\tilde{\mathbf{s}}}^H\left[k\right]\right\}{\mathbf{A}}^H\hfill \\ & =& \mathbf{A}{\mathbf{R}}_{\tilde{\mathbf{s}}}{\mathbf{A}}^H\hfill \end{array}$$

(12)

Because we assume that the baseband transmissions are uncorrelated with each other, ${\mathbf{R}}_{\tilde{\mathbf{s}}}=\mathrm{E}\left\{\tilde{\mathbf{s}}\left[k\right]{\tilde{\mathbf{s}}}^H\left[k\right]\right\}$ is a diagonal matrix. Then, by vectorizing the covariance matrix ${\mathbf{R}}_{\mathbf{y}}$, we can get

$$\begin{array}{ccc}\hfill {\mathbf{r}}_{\mathbf{y}}& =& \mathrm{vec}\left({\mathbf{R}}_{\mathbf{y}}\right)\hfill \\ & =& ({\mathbf{A}}^{*}\otimes \mathbf{A})\mathrm{vec}\left({\mathbf{R}}_{\tilde{\mathbf{s}}}\right)\hfill \\ & =& ({\mathbf{A}}^{*}\otimes \mathbf{A})\mathrm{diag}\left({\mathbf{R}}_{\tilde{\mathbf{s}}}\right)\hfill \\ & =& ({\mathbf{A}}^{*}\odot \mathbf{A}){\mathbf{r}}_{\tilde{\mathbf{s}}}\hfill \end{array}$$

(13)

where ${\mathbf{r}}_{\tilde{\mathbf{s}}}$ is a $K\times 1$ vector including the diagonal elements of ${\mathbf{R}}_{\tilde{\mathbf{s}}}$. Note that the above equations hold only when all transmissions are uncorrelated with each other. Here, we denote ${\mathbf{A}}_{\mathrm{ca}}\triangleq {\mathbf{A}}^{*}\odot \mathbf{A}$ as the manifold matrix of difference coarray. Obviously, the $(l,k)$-th element of ${\mathbf{A}}_{\mathrm{ca}}$ is given by

$$\begin{array}{ccc}\hfill {\left[{\mathbf{A}}_{\mathrm{ca}}\right]}_{l,k}& =& e^{-j2\pi f_k{d}_{n}cos\left(\theta \right)/c}·e^{j2\pi f_k{d}_{m}cos\left(\theta \right)/c}\hfill \\ & =& e^{j2\pi f_k(d_m-d_n)cos\left(\theta \right)/c},l=(n-1)(N+2M-1)+m\hfill \end{array}$$

(14)

where $e^{-j2\pi f_k{d}_{n}cos\left(\theta \right)/c}$ is the $(n,k)$-th element of ${\mathbf{A}}^{*}$ and $e^{j2\pi f_k{d}_{m}cos\left(\theta \right)/c}$ is the $(m,k)$-th element of $\mathbf{A}$. $d_m,d_n\in S_1\cup S_2$. So we can easily get that the aperture of the difference coarray extends from $-(2M-1)Nd$ to $(2M-1)Nd$. But this difference coarray is not filled, there are some holes in it. By referring to [11], we know that it has a contiguous set of elements from $-(MN+M-1)d$ to $(MN+M-1)d$, which acts like a filled virtual uniform linear array (ULA). To make it more clear, we first define the weight function at each element position of the virtual ULA.

Definition 1.

(The weight function, $w\left(l\right),-(MN+M-1)\le l\le MN+M-1$). Consider a co-prime array with its co-prime parameters M and N. Let $\mathbb{S}=S_1\cup S_2$ be the element position set of physical array and $\mathbb{D}=\left\{l\right|-(MN+M-1)\le l\le MN+M-1\}$ be the element position set of virtual ULA. The weight function is the number of pairs $(n,m)\in {\mathbb{S}}^2$ which have a difference $m-n=l$, defined as

$$\begin{array}{ccc}\hfill w\left(l\right)& =& \mathit{Cardinality}\mathit{of}\mathit{the}\mathit{set}\mathbb{M}\left(l\right)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \mathbb{M}\left(l\right)& =& \left\{(n,m)\in {\mathbb{S}}^2|m-n=l\right\}\hfill \end{array}$$

(16)

For example, if we choose $M=2,N=3$, then $\mathbb{S}=\{0,2,3,4,6,9\}$ and $\mathbb{D}=\left\{l\right|-7\le l\le 7\}$. The weight function $w\left(l\right)$ satisfies $w(-7)=1,w(-6)=2$, etc.

Based on (13), we denote ${\left[{\mathbf{r}}_{\mathbf{y}}\right]}_{(n,m)}$ as the row of ${\mathbf{r}}_{\mathbf{y}}$, which is produced by the n-th row of ${\mathbf{A}}^{*}$ and the m-th row of $\mathbf{A}$ via Khatri-Rao product. So by choosing the continuous lags from $-(MN+M-1)$ to $(MN+M-1)$ and taking the weight function into account, we can get a virtual ULA model

$$\begin{array}{c}\hfill {\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}={\mathbf{B}}_{\mathrm{ULA}}{\mathbf{r}}_{\tilde{\mathbf{s}}}\end{array}$$

(17)

where the i-th ($i=l+MN+M$) element of ${\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}$ is

$$\begin{array}{c}\hfill {\left[{\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}\right]}_i=\frac{1}{w\left(l\right)}\sum _{(n,m)\in \mathbb{M}\left(l\right)}{\left[{\mathbf{r}}_{\mathbf{y}}\right]}_{(n,m)},l\in \mathbb{D}\end{array}$$

(18)

and ${\mathbf{B}}_{\mathrm{ULA}}$, namely the array manifold of the virtual ULA, has the structure

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{e}^{j2\pi f_1{\tilde{\tau }}_{-MN-M+1}}& \cdots & e^{j2\pi f_K{\tilde{\tau }}_{-MN-M+1}}\\ ⋮& \ddots & ⋮\\ e^{j2\pi f_1{\tilde{\tau }}_{MN+M-1}}& \cdots & e^{j2\pi f_K{\tilde{\tau }}_{MN+M-1}}\end{array}\right)\end{array}$$

where ${\tilde{\tau }}_n=dncos\left(\theta \right)/c,-(MN+M-1)\le n\le MN+M-1$. Based on (17), we have the following sufficient condition for unique solution to ${\left\{f_i\right\}}_{i=1}^K$. Obviously, ${\mathbf{r}}_{\tilde{\mathbf{s}}}$ is equivalent to K coherent sources with only one snapshot.

Theorem 1.

Consider a co-prime array consisting of $N+2M-1$ sensor elements which can be transformed into a virtual filled ULA in (17). If

$$\begin{array}{c}\hfill MN+M>K\end{array}$$

(19)

then (17) has a unique solution of ${\left\{f_i\right\}}_{i=1}^K$.

Proof. 

Because ${\mathbf{r}}_{\tilde{\mathbf{s}}}$ is equivalent to K coherent sources, we have $\dim \left(\mathrm{span}\left({\mathbf{r}}_{\tilde{\mathbf{s}}}\right)\right)=1$. And the virtual array acts as a filled ULA with inter-element spacing d satisfies $d\le \frac{c}{2f_{\mathrm{Nyq}}}$ and the number of virtual elements is $2MN+2M-1$. So we refer the reader to [10] which deals with the physical ULA case. With these substitutions, the result follows from Theorem 1 in [10]. ☐

Next, we define

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{y}\mathrm{ULA}}& =& {\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}{\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}^H\hfill \\ & =& {\mathbf{B}}_{\mathrm{ULA}}{\mathbf{r}}_{\tilde{\mathbf{s}}}{\mathbf{r}}_{\tilde{\mathbf{s}}}^H{\mathbf{B}}_{\mathrm{ULA}}^H\hfill \\ & =& {\mathbf{B}}_{\mathrm{ULA}}{\mathbf{R}}_{\tilde{\mathbf{s}}\tilde{\mathbf{s}}}{\mathbf{B}}_{\mathrm{ULA}}^H\hfill \end{array}$$

(20)

Because

$$\begin{array}{c}\hfill \mathrm{rank}\left({\mathbf{R}}_{\mathbf{y}\mathrm{ULA}}\right)=\mathrm{rank}\left({\mathbf{R}}_{\tilde{\mathbf{s}}\tilde{\mathbf{s}}}\right)=1\end{array}$$

(21)

we need to implement a spatial smoothing step to enhance the rank of the covariance matrix. As analyzed above, the virtual ULA has the element position from $-(MN+M-1)d$ to $(MN+M-1)d$. Now, we divide this virtual array into $MN+M$ overlapping subarrays, each with $MN+M$ elements. The i-th subarray has sensors located at

$$\begin{array}{c}\hfill \left\{\right(-i+1+n)d,n=0,1,\dots ,MN+M-1\}\end{array}$$

(22)

which corresponds to the $(MN+M+1-i)$-th to $(2MN+2M-i)$-th rows of ${\mathbf{r}}_{\mathbf{y}\mathrm{ULA}}$. So we have

$$\begin{array}{c}\hfill {\mathbf{r}}_{\mathbf{y}\mathrm{ULA}i}={\mathbf{B}}_{\mathrm{ULA}i}{\mathbf{r}}_{\tilde{\mathbf{s}}}\end{array}$$

(23)

where ${\mathbf{B}}_{\mathrm{ULA}i}$ is a $(MN+M)\times K$ matrix consisting of the $(MN+M+1-i)$-th to $(2MN+2M-i)$-th rows of ${\mathbf{B}}_{\mathrm{ULA}}$ which has the structure

$$\begin{array}{c}\hfill {\mathbf{B}}_{\mathrm{ULA}i}=\left(\begin{array}{ccc}{e}^{j2\pi f_1(1-i)dcos\left(\theta \right)/c}& \cdots & e^{j2\pi f_K(1-i)dcos\left(\theta \right)/c}\\ ⋮& \ddots & ⋮\\ e^{j2\pi f_1(MN+M-i)dcos\left(\theta \right)/c}& \cdots & e^{j2\pi f_K(MN+M-i)dcos\left(\theta \right)/c}\end{array}\right)\end{array}$$

(24)

Obviously, from the above structure, we can get

$$\begin{array}{c}\hfill {\mathbf{B}}_{\mathrm{ULA}i}={\mathbf{B}}_{\mathrm{ULA}1}{\mathbf{\Phi }}^{i-1}\end{array}$$

(25)

where $\mathbf{\Phi }$ is a diagonal matrix with its diagonal elements as $\{e^{-j2\pi f_{1}dcos\left(\theta \right)/c},\dots ,e^{-j2\pi f_{K}dcos\left(\theta \right)/c}\}$. So, we rewrite (23) as

$$\begin{array}{c}\hfill {\mathbf{r}}_{\mathbf{y}\mathrm{ULA}i}={\mathbf{B}}_{\mathrm{ULA}1}{\mathbf{\Phi }}^{i-1}{\mathbf{r}}_{\tilde{\mathbf{s}}}\end{array}$$

(26)

Then, we can get the spatially smoothed matrix

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathrm{ss}}& =& \frac{1}{MN+M}\sum _{i=1}^{MN+M}{\mathbf{R}}_{\mathbf{y}\mathrm{ULA}i}\hfill \end{array}$$

(27)

where

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{y}\mathrm{ULA}i}& =& {\mathbf{r}}_{\mathbf{y}\mathrm{ULA}i}{\mathbf{r}}_{\mathbf{y}\mathrm{ULA}i}^H\hfill \\ & =& {\mathbf{B}}_{\mathrm{ULA}1}{\mathbf{\Phi }}^{i-1}{\mathbf{r}}_{\tilde{\mathbf{s}}}{\mathbf{r}}_{\tilde{\mathbf{s}}}^H{\left({\mathbf{\Phi }}^{i-1}\right)}^H{\mathbf{B}}_{\mathrm{ULA}1}^H\hfill \end{array}$$

(28)

The spatially smoothed matrix ${\mathbf{R}}_{\mathrm{ss}}$ can be used to estimate carrier frequencies by the following theorem.

Theorem 2.

Consider the spatially smoothed matrix ${\mathbf{R}}_{\mathrm{ss}}$ in (27) and define a $K\times K$ diagonal matrix ${\Lambda }_s$ with its diagonal elements as the covariances of K targets. Then, we have

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathrm{ss}}& =& {\tilde{\mathbf{R}}}^2\hfill \\ \hfill \tilde{\mathbf{R}}& =& \frac{1}{\sqrt{MN+M}}{\mathbf{B}}_{\mathrm{ULA}1}\Lambda {\mathbf{B}}_{\mathrm{ULA}1}^H\hfill \end{array}$$

(29)

Proof. 

The proof follows the same lines as Theorem 1 in [12], only substituting the values of ${\mathbf{B}}_{\mathrm{ULA}1}$ and ${\Lambda }_s$ in our paper. ☐

By decomposing $\tilde{\mathbf{R}}$ using the singular value decomposition, we have

$$\begin{array}{c}\hfill \tilde{\mathbf{R}}=\left[{\mathbf{U}}_1{\mathbf{U}}_2\right]\left[\begin{array}{cc}{\Lambda }_s& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]{\mathbf{V}}^H\end{array}$$

(30)

The columns of the matrix $\left[{\mathbf{U}}_1{\mathbf{U}}_2\right]$ are the left singular vectors of $\tilde{\mathbf{R}}$, where ${\mathbf{U}}_1$ contains the vectors corresponding to the first K singular values, ${\Lambda }_s$ is a $K\times K$ diagonal matrix with the K first singular values of $\tilde{\mathbf{R}}$, and $\mathbf{V}$ contains the right singular vector of $\tilde{\mathbf{R}}$. Based on (29) and (30), we know that there exists an invertible $K\times K$ matrix $\mathbf{T}$ such that

$$\begin{array}{c}\hfill {\mathbf{U}}_1={\mathbf{B}}_{\mathrm{ULA}1}\mathbf{T}\end{array}$$

(31)

Consider the first $MN+M-1$ rows of ${\mathbf{U}}_1$, we have

$$\begin{array}{c}\hfill {\mathbf{U}}_{11}={\mathbf{B}}_{\mathrm{ULA}11}\mathbf{T}\end{array}$$

(32)

Similarly, we can have the last $MN+M-1$ rows of ${\mathbf{U}}_1$

$$\begin{array}{c}\hfill {\mathbf{U}}_{12}={\mathbf{B}}_{\mathrm{ULA}12}\mathbf{T}\end{array}$$

(33)

where ${\mathbf{B}}_{\mathrm{ULA}11}$ is the virtual sub-array consisting of element positions $\{0,d,\dots ,(MN+M-2\left)d\right\}$ and ${\mathbf{B}}_{\mathrm{ULA}12}$ is the virtual sub-array consisting of elements positions $\{d,2d,\dots ,(MN+M-1\left)d\right\}$. So, we can get the relationship between ${\mathbf{B}}_{\mathrm{ULA}11}$ and ${\mathbf{B}}_{\mathrm{ULA}12}$ as

$$\begin{array}{c}\hfill {\mathbf{B}}_{\mathrm{ULA}11}={\mathbf{B}}_{\mathrm{ULA}12}\mathbf{\Phi }\end{array}$$

(34)

where $\mathbf{\Phi }$ is a diagonal matrix which is defined in (26). So we rewrite (32) as

$$\begin{array}{c}\hfill {\mathbf{U}}_{11}={\mathbf{B}}_{\mathrm{ULA}12}\mathbf{\Phi }\mathbf{T}\end{array}$$

(35)

Here, we use the least squares recovery

$$\begin{array}{c}\hfill \mathbf{\Psi }={\mathbf{U}}_{12}^{†}{\mathbf{U}}_{11}\end{array}$$

(36)

Then, we have

$$\begin{array}{c}\hfill {\widehat{f}}_i=\angle \left(\mathrm{eig}\left({\mathbf{\Psi }}_{i,i}\right)\right)\frac{c}{2\pi dcos\left(\theta \right)}\end{array}$$

(37)

where ${\mathbf{\Psi }}_{i,i}$ is the i-th diagonal element of $\mathbf{\Psi }$.

Remark 1.

It can be seen from (37) that θ can not be equal to ${90}^{\circ }$. And the performance of carrier frequency estimation is affected by θ. Because $cos\left(\theta \right)$ is the denominator term in (37), a small $cos\left(\theta \right)$ will amplify the error which is caused by the calculation of Ψ. Assuming $\theta \in [0,{90}^{\circ })\cup ({90}^{\circ },{180}^{\circ }]$, the closer to $0^{\circ }$ or ${180}^{\circ }$ the impinging direction θ is, the smaller the estimation error is. Conversely, the closer to ${90}^{\circ }$, the larger the error is. In practice, if we know that θ is approaching ${90}^{\circ }$, we can add an adjustable known time delay line after each sensor which is equivalent to rotating the array with a known angle. If we denote the man-made time delay as $d_{n}cos\left({\theta }_a\right)/c$, then the denominator term in (37) is modified as $cos\left(\theta \right)+cos\left({\theta }_a\right)$. In the following discussion, we consider the case that θ is close to $0^{\circ }$ or ${180}^{\circ }$ for simplicity.

4.2. Signal Power Spectrum Recovery

Once the carrier frequencies are recovered, the steering matrix defined in (17) can be constructed. So in this subsection, we will first consider the power spectrum recovery of $\tilde{\mathbf{s}}\left[k\right]$. After that, we will investigate how to recover the power spectrum of $\mathbf{s}\left[k\right]$ from $\tilde{\mathbf{s}}\left[k\right]$.

By invoking (9), we consider the signal model in the frequency domain. Define the autocorrelation matrix of $\mathbf{Y}\left(e^{j2\pi f{T}_s}\right)$ as ${\mathbf{R}}_{\mathbf{Y}}\left(e^{j2\pi f{T}_s}\right)\triangleq \mathrm{E}\left\{\mathbf{Y}\left(e^{j2\pi f{T}_s}\right){\mathbf{Y}}^H\left(e^{j2\pi f{T}_s}\right)\right\}$. Similarly, define ${\mathbf{R}}_{\tilde{\mathbf{S}}}\left(f\right)\triangleq \mathrm{E}\left\{\tilde{\mathbf{S}}\left(f\right){\tilde{\mathbf{S}}}^H\left(f\right)\right\}$ and ${\mathbf{R}}_{\mathbf{S}}\left(f\right)\triangleq \mathrm{E}\left\{\mathbf{S}\left(f\right){\mathbf{S}}^H\left(f\right)\right\}$ for $f\in {\mathcal{F}}_s$. Then, we have

$$\begin{array}{c}\hfill {\mathbf{R}}_{\mathbf{Y}}\left(e^{j2\pi f{T}_s}\right)=\mathbf{A}{\mathbf{R}}_{\tilde{\mathbf{S}}}\left(f\right){\mathbf{A}}^H,f\in {\mathcal{F}}_s\end{array}$$

(38)

Due to the assumption that all transmissions are uncorrelated with each other, so ${\mathbf{R}}_{\tilde{\mathbf{S}}}\left(f\right),f\in {\mathcal{F}}_s$ and ${\mathbf{R}}_{\mathbf{S}}\left(f\right),f\in {\mathcal{F}}_s$ are both diagonal matrixes. Then, similar to the processing steps in (13) and (17), by vectorization, removing the redundancies and choosing the continuous lags, we can get the virtual array model in the frequency domain

$$\begin{array}{c}\hfill {\mathbf{r}}_{\mathbf{Y}\mathrm{ULA}}\left(f\right)={\mathbf{B}}_{\mathrm{ULA}}{\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f\right),f\in {\mathcal{F}}_s\end{array}$$

(39)

where ${\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f\right)$ is a $K\times 1$ vector which contains the diagonal elements of ${\mathbf{R}}_{\tilde{\mathbf{S}}}\left(f\right)$. Similarly, we denote ${\mathbf{r}}_{\mathbf{S}}\left(f\right)$ as a $K\times 1$ vector which contains the diagonal elements of ${\mathbf{R}}_{\mathbf{S}}\left(f\right)$. From (17), ${\mathbf{B}}_{\mathrm{ULA}}$ is a Vandermonde matrix, it has full column rank if and only if $(2MN+2M-1)\ge K$. Referring to Theorem 1, if the sufficient condition (19) is satisfied, ${\mathbf{B}}_{\mathrm{ULA}}$ will have full column rank. Then we can obtain the power spectrum of $\tilde{\mathbf{S}}\left(f\right)$ by inverting the steering matrix,

$$\begin{array}{c}\hfill {\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f\right)={\mathbf{B}}_{\mathrm{ULA}}^{†}{\mathbf{r}}_{\mathbf{Y}\mathrm{ULA}}\left(f\right),f\in {\mathcal{F}}_s\end{array}$$

(40)

As analyzed in the third section, $\tilde{\mathbf{S}}\left(f\right)$ is a cyclic and shifted version of $\mathbf{S}\left(f\right)$.

Consider the i-th transmission $S_i\left(f^{\prime }\right),f^{\prime }\in {\mathcal{F}}_p$. It holds that

$$\begin{array}{c}\hfill {\tilde{S}}_i\left(f^{\prime }\right)=c_{l_a}{S}_i(f^{\prime }-f_i-l_a{f}_p),\end{array}$$

(41)

where $l_a$ is known as

$$\begin{array}{c}\hfill l_a=⌊\frac{f_i+f^{\prime }+f_p/2}{f_p}⌋\end{array}$$

(42)

Then we have the relationship between the power spectrum of $S_i\left(f^{\prime }\right)$ and ${\tilde{S}}_i\left(f^{\prime }\right)$,

$$\begin{array}{ccc}\hfill {\left[{\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f^{\prime }\right)\right]}_i& \triangleq & \mathrm{E}\left\{{\tilde{S}}_i\left(f^{\prime }\right){\tilde{S}}_i^{*}\left(f^{\prime }\right)\right\}\hfill \\ & =& \mathrm{E}\left\{c_{l_a}{S}_i(f^{\prime }-f_i-l_a{f}_p)S_i^{*}(f^{\prime }-f_i-l_a{f}_p)c_{l_a}^{*}\right\}\hfill \\ & =& |c_{l_a}{|}^2\mathrm{E}\left\{S_i(f^{\prime }-f_i-l_a{f}_p)S_i^{*}(f^{\prime }-f_i-l_a{f}_p)\right\}\hfill \\ & =& |c_{l_a}{|}^2{\left[{\mathbf{r}}_{\mathbf{S}}(f^{\prime }-f_i-l_a{f}_p)\right]}_i\hfill \end{array}$$

(43)

where ${\left[{\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f^{\prime }\right)\right]}_i$ is the i-th element of the $K\times 1$ vector ${\mathbf{r}}_{\tilde{\mathbf{S}}}\left(f^{\prime }\right)$ and ${\left[{\mathbf{r}}_{\mathbf{S}}(f^{\prime }-f_i-l_a{f}_p)\right]}_i$ is the the i-th element of the $K\times 1$ vector ${\mathbf{r}}_{\mathbf{S}}(f^{\prime }-f_i-l_a{f}_p)$. After a change of variables,

$$\begin{array}{c}\hfill {\left[{\mathbf{r}}_{\mathbf{S}}\left(f^{\prime }\right)\right]}_i=\frac{1}{|c_{l_a}{|}^2}{\left[{\mathbf{r}}_{\tilde{\mathbf{S}}}(f^{\prime }-f_i-l_a{f}_p)\right]}_i,f^{\prime }\in {\mathcal{F}}_p\end{array}$$

(44)

Observing (43) and (44), the equality in (44) holds if and only if $c_{l_a}\ne 0$.

4.3. Comparison with Previous MWC Systems

By referring to conventional MWC [7] and ULA-based MWC [10], we can have the following conclusions. Firstly, we compare our proposed CA-based MWC with ULA-based MWC and conventional MWC. Our method processes the signal in the co-array domain, while the latter two methods process signal in physical sensor (channel) domain. That means, if we fix the number of physical sensors or channels as $N+2M-1$, then our proposed CA-based method can produce a virtual ULA which has $MN+M$ elements. It is much larger than that of ULA-based MWC and conventional MWC which can increase the system’s robustness to noise. Another difference is that our method can directly recover the power spectrum of impinging signal, while the latter methods must first recover the signal itself after which the power spectrum is calculated. Here, we need to point out a disadvantage as shown in (13) that the impinging signal for our method must be uncorrelated with each other. Secondly, we compare CA-based MWC, ULA-based MWC with conventional MWC. In our proposed CA-based MWC and ULA-based MWC, carrier frequencies are first estimated, then the baseband transmissions are estimated. For conventional MWC, there’s no need to estimate carrier frequencies, all RF signals are estimated directly. In addition, each channel of CA-based MWC and ULA-based MWC is corrupted by independent noise, while each channel of conventional MWC is corrupted by the same noise. Lastly, we compare our method with ULA-based MWC. Besides a difference about the number of sensors, another difference is that our proposed system is a sparse array system while ULA-based MWC is a filled array system. As we all know, the closer the sensors are, the more correlated their samples are, which can affect the performance. The differences among these three methods are shown clearly in

Table 1. Comparisons of CA-based MWC, ULA-based MWC and MWC.

	CA-based MWC	ULA-based MWC	MWC
Periodic function	Same for all sensors	Same for all sensors	Different for all sensors
Number of channels	$N+2M-1$	$N+2M-1$	$N+2M-1$
Noise of each channel	Independent	Independent	Same
Processing domain	Coarray	Physical	Physical
Number of virtual elements	$MN+M$	×	×
Reconstruction	Power spectrum	Original signal	Original signal

where × denotes “Not exist”.

5. Robustness Analysis for Imperfect $p\left(t\right)$ Among Different Channels

Because the mixing function $p\left(t\right)$ in each channel is produced by a power divider, there may be amplitude and phase imbalances of $p\left(t\right)$ among different channels. In this section, we will analyze this issue.

Here, we denote the amplitude and phase imbalances of the n-th ($n=1,\dots ,N+2M-1$) channel as ${\delta }_n$ and ${\phi }_n$, respectively. We assume that the amplitude imbalances ${\left\{{\delta }_n\right\}}_{n=1}^{N+2M-1}$ are i. i. d random variables with Gaussian distribution $\mathcal{N}(0,{\sigma }_{\delta }^2)$ and the phase imbalances ${\left\{{\phi }_n\right\}}_{n=1}^{N+2M-1}$ are i. i. d random variables distributed uniformly in $[-\rho /2,\rho /2]$. Moreover, ${\left\{{\phi }_n\right\}}_{n=1}^{N+2M-1}$ and ${\left\{{\delta }_n\right\}}_{n=1}^{N+2M-1}$ are independent with each other. So, for the n-th channel, the imperfect mixing function $p_n\left(t\right)$ can be modeled as

$$\begin{array}{c}\hfill p_n\left(t\right)=(1+{\delta }_n)e^{-j{\phi }_n}p\left(t\right)\end{array}$$

(45)

where $p\left(t\right)$ is the ideal mixing function. We denote the Fourier coefficients of $p_n\left(t\right)$ as ${\tilde{c}}_{n,l}$. Then, we have

$$\begin{array}{ccc}\hfill {\tilde{c}}_{n,l}& =& \frac{1}{T_p}{\int }_0^{T_p}{p}_n\left(t\right)e^{-j2\pi lt/T_p}\mathrm{d}t\hfill \\ & =& \frac{1}{T_p}{\int }_0^{T_p}(1+{\delta }_n)e^{-j{\phi }_n}p\left(t\right)e^{-j2\pi lt/T_p}\mathrm{d}t\hfill \\ & =& (1+{\delta }_n)e^{-j{\phi }_n}{c}_l\hfill \end{array}$$

(46)

where $c_l$ is the ideal Fourier coefficients of $p\left(t\right)$. Based on (8), we can obtain the perturbed model

$$\begin{array}{ccc}\hfill {\widehat{Y}}_n\left(e^{j2\pi f{T}_s}\right)& =& \sum _{l=-L_0}^{L_0}(1+{\delta }_n)e^{-j{\phi }_n}{c}_l\sum _{i=1}^K{S}_i(f-f_i-l{f}_p)e^{j2\pi f_i{\tau }_n}\hfill \\ & =& \sum _{i=1}^K{\tilde{S}}_i\left(f\right)e^{j2\pi f_i{\tau }_n}(1+{\delta }_n)e^{-j{\phi }_n},f\in {\mathcal{F}}_s\hfill \end{array}$$

(47)

We write (47) in matrix form as

$$\begin{array}{c}\hfill \widehat{\mathbf{Y}}\left(e^{j2\pi f{T}_s}\right)=\sum _{i=1}^K({\mathbf{a}}_i\circ \mathbf{p}){\tilde{S}}_i\left(f\right),f\in {\mathcal{F}}_s\end{array}$$

(48)

where ${\mathbf{a}}_i$ is the i-th column of $\mathbf{A}$ defined in (10) and $\mathbf{p}={[(1+{\delta }_1)e^{-j{\phi }_1},\dots ,(1+{\delta }_{N+2M-1})e^{-j{\phi }_{N+2M-1}}]}^T$. Then, we can easily get the temporal perturbed model under sampling rate $f_s$ as

$$\begin{array}{c}\hfill \widehat{\mathbf{y}}\left[k\right]=\sum _{i=1}^K({\mathbf{a}}_i\circ \mathbf{p}){\tilde{s}}_i\left[k\right]\end{array}$$

(49)

Now we consider the covariance matrix ${\widehat{\mathbf{R}}}_{\mathbf{y}}=\mathrm{E}\left\{\widehat{\mathbf{y}}\left[k\right]{\widehat{\mathbf{y}}}^H\left[k\right]\right\}$. The structure of ${\widehat{\mathbf{R}}}_{\mathbf{y}}$ is given by the following corollary.

Theorem 3.

In the CA-based MWC system, if we consider the amplitude and phase imbalances for mixing function $p\left(t\right)$ among different channels, then the perturbed covariance matrix is given by

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{R}}}_{\mathbf{y}}& =& \sum _{i=1}^K{\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i\left({\mathbf{a}}_i{\mathbf{a}}_i^H\right)\circ \mathbf{E}\hfill \\ & =& {\mathbf{R}}_{\mathbf{y}}\circ \mathbf{E}\hfill \end{array}$$

(50)

where ${\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i$ is the i-th element of ${\mathbf{r}}_{\tilde{\mathbf{s}}}$ defined in (13). $\mathbf{E}$ is a matrix with $1+{\sigma }_{\delta }^2$ on its diagonal and ${\mathrm{sinc}}^2(\rho /2)$ elsewhere.

Proof. 

The perturbed covariance matrix can be obtained by

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{R}}}_{\mathbf{y}}& =& \mathrm{E}\left\{\widehat{\mathbf{y}}\left[k\right]{\widehat{\mathbf{y}}}^H\left[k\right]\right\}\hfill \\ & =& \mathrm{E}\left\{\sum _{i=1}^K({\mathbf{a}}_i\circ \mathbf{p})({\mathbf{a}}_i^H\circ {\mathbf{p}}^H){\left|{\tilde{s}}_i\left[k\right]\right|}^2\right\}\hfill \\ & =& \mathrm{E}\left\{\sum _{i=1}^K\left({\mathbf{a}}_i{\mathbf{a}}_i^H\right)\circ \left(\mathbf{p}{\mathbf{p}}^H\right){\left|{\tilde{s}}_i\left[k\right]\right|}^2\right\}\hfill \\ & =& \sum _{i=1}^K\left({\mathbf{a}}_i{\mathbf{a}}_i^H\right)\circ \mathrm{E}\left\{\left(\mathbf{p}{\mathbf{p}}^H\right)\right\}{\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i\hfill \end{array}$$

(51)

We denote $\mathbf{P}=\mathbf{p}{\mathbf{p}}^H$, our task is to find the expectation of $\mathbf{P}$. We first investigate the distribution of diagonal elements of $\mathbf{P}$. We have

$$\begin{array}{ccc}\hfill \mathrm{E}\left\{{\left[\mathbf{P}\right]}_{n,n}\right\}& =& \mathrm{E}\left\{(1+{\delta }_n)e^{-j{\phi }_n}{(1+{\delta }_n)}^{*}{e}^{j{\phi }_n}\right\}\hfill \\ & =& \mathrm{E}\left\{{(1+{\delta }_n)}^2\right\}\hfill \\ & =& \mathrm{E}\{1+{\delta }_n^2+2{\delta }_n\}\hfill \\ & =& 1+{\sigma }_{\delta }^2,n=1,\dots ,N+2M-1\hfill \end{array}$$

(52)

For the off-diagonal elements of $\mathbf{P}$, we have

$$\begin{array}{ccc}\hfill \mathrm{E}\left\{{\left[\mathbf{P}\right]}_{n,m}\right\}& =& \mathrm{E}\left\{(1+{\delta }_n)e^{-j{\phi }_n}{(1+{\delta }_m)}^{*}{e}^{j{\phi }_m}\right\}\hfill \\ & =& \mathrm{E}\left\{(1+{\delta }_n){(1+{\delta }_m)}^{*}{e}^{j({\phi }_m-{\phi }_n)}\right\}\hfill \\ & =& \mathrm{E}\left\{e^{j({\phi }_m-{\phi }_n)}\right\}\hfill \\ & =& \mathrm{E}\left\{e^{j{\beta }_{mn}}\right\},n,m=1,\dots ,N+2M-1,m\ne n\hfill \end{array}$$

(53)

where ${\beta }_{mn}\triangleq {\phi }_m-{\phi }_n$, namely the difference of two independent random variables with uniform distribution in $[-\rho /2,\rho /2]$. So we can easily obtain that the probability density function (pdf) of ${\beta }_{mn}$ is a triangular function in the interval $[-\rho ,\rho ]$:

$$f_{{\beta }_{mn}}\left(\beta \right)=\{\begin{array}{cc}\hfill \frac{\rho +\beta }{{\rho }^2},& -\rho \le \beta <0\\ \hfill \frac{\rho -\beta }{{\rho }^2},& 0\le \beta \le \rho \end{array}$$

(54)

By integration, the expectation of ${\left[\mathbf{P}\right]}_{n,m}$ is calculated as

$$\begin{array}{ccc}\hfill \mathrm{E}\left\{{\left[\mathbf{P}\right]}_{n,m}\right\}& =& \mathrm{E}\left\{e^{j{\beta }_{mn}}\right\}\hfill \\ & =& {\int }_{-\rho }^{\rho }{e}^{j\beta }{f}_{{\beta }_{mn}}\left(\beta \right)\mathrm{d}\beta \hfill \\ & =& {\left(\frac{sin(\rho /2)}{\rho /2}\right)}^2\hfill \\ & =& {\mathrm{sinc}}^2(\rho /2)\hfill \end{array}$$

(55)

In (50), $\mathbf{E}=\mathrm{E}\left\{\mathbf{P}\right\}$ which completes the proof. ☐

Using Theorem 3, we can know that the structured noise caused by the amplitude and phase imbalances of $p\left(t\right)$ can be written as

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}}& =& \sum _{i=1}^K{\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i\left({\mathbf{a}}_i{\mathbf{a}}_i^H\right)\circ (\mathbf{1}-\mathbf{E})\hfill \\ & =& {\mathbf{R}}_{\mathbf{y}}\circ (\mathbf{1}-\mathbf{E})\hfill \end{array}$$

(56)

where $\mathbf{1}$ is an all-ones matrix with the size $(N+2M-1)\times (N+2M-1)$. Obviously, the structured noise acts like the additive colored correlated noise which can degrade the final estimation performance.

Corollary 1.

Let $\rho \ll 1$, then the deviation of perturbed covariance matrix from the ideal covariance matrix is given by

$$\begin{array}{c}\hfill {∥{\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}}∥}_F={\parallel {\mathbf{r}}_{\tilde{\mathbf{s}}}\parallel }_2\sqrt{(N+2M-1)[{\sigma }_{\delta }^4+(N+2M-2){\rho }^4/144]}\end{array}$$

(57)

where ${\parallel .\parallel }_F$ is the Frobenius norm.

Proof. 

Obviously, the n-th diagonal element of ${\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}}$ have the following equation

$$\begin{array}{ccc}\hfill {\left|{({\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}})}_{n,n}\right|}^2& =& \sum _{i=1}^K{\left|{\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i\right|}^2{\sigma }_{\delta }^4\hfill \\ & =& {\sigma }_{\delta }^4{\parallel {\mathbf{r}}_{\tilde{\mathbf{s}}}\parallel }_2^2\hfill \end{array}$$

(58)

For the $(n,m)$-th off-diagonal element of ${\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}}$, we have

$$\begin{array}{ccc}\hfill {\left|{({\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}})}_{n,m}\right|}^2& =& \sum _{i=1}^K{\left|{\left[{\mathbf{r}}_{\tilde{\mathbf{s}}}\right]}_i\right|}^2{\left(1-{\mathrm{sinc}}^2(\rho /2)\right)}^2\hfill \\ & =& {\left(1-{\mathrm{sinc}}^2(\rho /2)\right)}^2{\parallel {\mathbf{r}}_{\tilde{\mathbf{s}}}\parallel }_2^2\hfill \end{array}$$

(59)

Then, we can easily get

$$\begin{array}{ccc}\hfill {∥{\mathbf{R}}_{\mathbf{y}}-{\widehat{\mathbf{R}}}_{\mathbf{y}}∥}_F& =& \parallel {\mathbf{r}}_{\tilde{\mathbf{s}}}{\parallel }_2\sqrt{(N+2M-1){\sigma }_{\delta }^4+(N+2M-2)(N+2M-1){\left(1-{\mathrm{sinc}}^2(\rho /2)\right)}^2}\hfill \\ & =& \parallel {\mathbf{r}}_{\tilde{\mathbf{s}}}{\parallel }_2\sqrt{(N+2M-1)[{\sigma }_{\delta }^4+(N+2M-2){\left(1-{\mathrm{sinc}}^2(\rho /2)\right)}^2]}\hfill \end{array}$$

(60)

Due to the assumption that $\rho \ll 1$, so $1-{\mathrm{sinc}}^2(\rho /2)$ in (60) can be approximated as ${\rho }^2/12$. Then we can obtain (57). ☐

6. Choice of Co-Prime Parameters

In Theorem 1, we know how to choose co-prime parameters M and N to satisfy the unique recovery condition. In application, the number of total sensors is usually fixed, namely $N+2M-1=Q$. As analyzed in Section 3, a virtual array whose elements are given by the difference coarray can be produced by vectorizing the data covariance matrix. Although the coarray of co-prime array has some missing elements or ’holes’, we can only employ that part of the coarray which has the continuous elements from $-(MN+M-1)$ to $MN+M-1$. So it acts like a filled ULA which has $2MN+2M-1$ elements. If we can increase the aperture of the virtual ULA, then the number of transmissions which can be detected also increases. So we have the following optimization

$$\begin{array}{c}\hfill \underset{M,N}{max}2MN+2M-1\\ \hfill \mathrm{s}.\mathrm{t}.2M+N-1=Q,M<N\\ \hfill \gcd (M,N)=1,M>0,N>0\end{array}$$

(61)

where gcd(.) is a function which returns the greatest common divisor. According to Cauchy inequality [13], we have the optimal solutions $M^{*}$ and $N^{*}$ satisfying $2M^{*}=N^{*}+1=(Q+2)/2$, namely $M^{*}=(Q+2)/4$ and $N^{*}=Q/2$ without considering the constrictions. Taking that M and N are co-prime integers into account, we can search the satisfying optimal solution around $(M^{*},N^{*})$.

7. Numerical Results

In this section, we assume $f_{\mathrm{Nyq}}=10\mathrm{GHz}$, $\theta =0^{\circ }$ and $B_{\max }=50\mathrm{MHz}$. In traditional MWC system, because the signal in each sampling channel comes from the same sensor, we can assume that all sampling channels are corrupted by the same additive Gaussian white noise. However, in our proposed array-based MWC system, we can assume that each sampling channel has uncorrelated Gaussian white noise because we use different sensors to receive the signal in different sampling channels. For simplicity, we choose $f_s=f_p=1.3B_{\max }$ in all simulations. The Matlab codes of conventional MWC system can be referred to [14].

7.1. Detection Performance

In this experiment, we will examine the detection performance of our proposed method. Here, we set $N=3$ and $M=2$. So the first subarray has the sensor position $[0,3,6,9]d$, and the second subarray has the sensor position $[0,2,4]d$, with d taken as half of the wavelength corresponding to the Nyquist rate. The first sensor of these two subarrays are co-located. So the total number of sensors in our co-prime array is $N+2M-1=6$. Here, we compare our method with ULA-based MWC system. In the latter case, the number of physical sensors is also set as 6, which is the same as that of co-prime array. The Signal-to-noise rate (SNR) is set as 10 dB. Firstly, we assume that there are 3 transmissions with the carrier frequencies as $[-1,1,1.5]$ $\mathrm{GHz}$. In Figure 4, it can be seen that 3 transmissions can be detected clearly both in traditional ULA-based MWC system and our proposed CA-based MWC system. In addition, in Figure 4 (bottom), one transmission is a little higher than 1.5 GHz. This is caused by the randomness of $p\left(t\right)$. For a specific run of CA-based MWC, it is possible that some elements in ${\left\{c_l\right\}}_{l=-L_0}^{L_0}$ are close to zero, which will degrade the estimating performance. But the probability is small. Then, we increase the number of transmissions to $K=7$ with the carrier frequencies as $[-4,-2,-0.5,1,1.5,3,4]$ $\mathrm{GHz}$. We can see from Figure 5 that our proposed method can detect all the transmissions successfully while traditional ULA-based MWC fails.

7.2. Sensing Accuracy

In the second experiment, we test the spectrum sensing accuracy of our proposed CA-based MWC system, compared with ULA-based MWC [10], traditional MWC [6] and non-compressive technique [9]. Because the number of uncorrelated transmissions which ULA-based MWC system can detect is no more than the number of physical sensors [9], we assume there are 3 transmissions with the carrier frequencies $[-1,1,1.5]$ $\mathrm{GHz}$. The SNR is ranged from 0 dB to 20 dB with step size 2 dB. The other simulation parameters are the same as those in the first experiment, e.g. $N=3,M=2$. So the total number of sensors used is $N+2M-1=6$. The number of sensors in ULA-based MWC system is also chosen as 6. In addition, the non-compressive method comes from [9] where we make the number of sampling channels equal to the decimator factor, namely $⌈\frac{f_{\mathrm{Nyq}}}{f_{\mathrm{sys}}}⌉=154$. Here, we define the spectrum sensing accuracy as

$$\begin{array}{c}\hfill \epsilon =\frac{1}{J}\sum _{j=1}^J\sqrt{\frac{1}{K}\sum _{i=1}^K\left[{\left({\widehat{f}}_{i,j}^l-f_i^l\right)}^2+{\left({\widehat{f}}_{i,j}^u-f_i^u\right)}^2\right]}\end{array}$$

(62)

where ${\widehat{f}}_{i,j}^l$ and ${\widehat{f}}_{i,j}^u$ are the estimated lower boundary and upper boundary of the i-th transmission in the j-th Monte Carlo simulation, respectively. $f_i^l$ and $f_i^u$ are the true lower boundary and upper boundary of i-th transmission. $f_i$ is the true carrier frequency of the i-th transmission. Here, 200 Monte Carlo simulations are used.

It can be seen that in Figure 6 non-compressive method has the best performance because it uses all Nyquist samples in signal processing. Array-based MWC methods, including CA-based MWC and ULA-based MWC, outperform traditional MWC method. This is on one hand due to the noise averaging in array-based system. On the other hand, array-based MWC methods have a two-step processing procedure, namely estimating carrier frequencies before recovering baseband transmissions, while conventional MWC method recovers the RF signal directly. Figure 6 also shows that our method has better performance than ULA-based MWC system. This is because in our method, a virtual ULA which has much larger aperture than ULA-based MWC system can be produced. A larger amount of sensors can increase the robustness to noise and allows it to handle a greater amount of sources. It can enhance CR performance.

Figure 6 is the result in the case of $\theta =0^{\circ }$. As analyzed in Section 4.1, the choice of $\theta$ can affect the performance of carrier frequency estimation, which will affect the sensing accuracy eventually. To investigate the sensing accuracy under different $\theta$ where $\theta \in [0^{\circ },{90}^{\circ })\cup ({90}^{\circ },{180}^{\circ }]$, we fix SNR = 10 dB. For the interval $[0^{\circ },{90}^{\circ })$, we choose $\theta$ from $0^{\circ }$ to ${80}^{\circ }$ with step size ${10}^{\circ }$. For the interval $({90}^{\circ },{180}^{\circ }]$, we choose $\theta$ from ${100}^{\circ }$ to ${180}^{\circ }$ with step size ${10}^{\circ }$. The other parameters are the same as those of Figure 6. It can be seen from Figure 7 that the closer to $0^{\circ }$ or ${180}^{\circ }$ the impinging direction $\theta$ is, the smaller the estimation error is. Conversely, the closer to ${90}^{\circ }$, the larger the error is. And our proposed CA-based MWC method has a better performance than ULA-based MWC. In addition, we can see that the direction-of-arrival has no effect on conventional MWC method. This is because conventional MWC is a one-sensor, multichannel system, it does not use the time delays among different sensors which are caused by impinging angles. So conventional MWC method outperforms our proposed CA-based MWC when direction-of-arrival is close to ${90}^{\circ }$.

7.3. Parameter Choice Demonstration

In the third experiment, we demonstrate the co-prime parameter choice criterion in the fifth section. Unlike the above two experiments, we assume there are 22 physical sensors, namely $Q=2M+N-1=22$ where $M<N$. So we can list all the co-prime pairs $(M,N)$ which satisfy the conditions, that’s $(1,21),(2,19),(3,17),(4,15),(5,13),(6,11),(7,9)$. For each co-prime pair, the virtual ULA aperture is $2MN+2M-1$. Figure 8 shows the virtual ULA apertures for different $(M,N)$. It can be seen that the virtual ULA apertures vary in different co-prime parameters and the largest virtual ULA aperture appears when $(M,N)=(6,11)$. This result coincides with our co-prime parameter choice criterion which says the largest virtual ULA aperture can be obtained when $(M,N)=(\frac{Q+2}{4},\frac{Q}{2})=(6,11)$.

Here, we also investigate the sensing performance under different co-prime parameters $(M,N)$ for a fixed number of physical sensors $Q=22$. Because the co-prime parameter $(1,21)$ has no sense in practice, we only compare the sensing accuracy under the co-prime parameters $(2,19),(3,17),(4,15),(5,13),(7,9),(6,11)$. The other simulation parameters are the same as that of Section 7.2. It can been from Figure 9 that different co-prime parameters have different sensing performance. The larger the virtual ULA aperture is, the better the performance is.

7.4. Minimal System Sampling Rate Comparison

In the last experiment, we compare our proposed CA-based MWC with ULA-based MWC on minimal system sampling rate. As analyzed above, the system sampling rate $f_{\mathrm{sys}}$ is determined by the number of sampling channels Q and the low sampling rate $f_s$ of each channel, namely $f_{\mathrm{sys}}=Q\times f_s$. As shown in Section 3, our only requirement for $f_s$ is $f_s\ge f_p\ge B_{\max }$. Due to the existence of edge effect for signals, we make $f_s$ ($f_p$) a little larger than $B_{\max }$, that’s $f_s=f_p=1.3B_{\max }=65$ $\mathrm{MHz}$. Obviously, $f_s$ is the same in both ULA-based MWC system and CA-based MWC system. The difference between these two systems is the required number of sampling channels Q for a fixed number of targets K. In the ULA-based MWC system, referring to [10], the minimal number of sampling channels required for K uncorrelated targets is $K+1$. So the corresponding minimal system sampling rate is $f_{\mathrm{sys}}=(K+1)f_s$. In our CA-based MWC system, based on Theorem 1, we can obtain the minimal number of sampling channels by solving the following optimization for a fixed K,

$$\begin{array}{c}\underset{M,N}{min}2M+N-1\\ \mathrm{s}.\mathrm{t}.MN+M>K,\gcd (M,N)=1,\\ N>M>0\end{array}$$

(63)

Denote the optimal solution as $(M^{*},N^{*})$, then the minimal system sampling rate is $f_{\mathrm{sys}}=Q\times f_s=(2M^{*}+N^{*}-1)f_s$. Here, we make the number of targets from 10 to 100 with step size 10. It can be seen from Figure 10 that the system sampling rate of our proposed system is reduced largely compared to the ULA-based MWC system, especially when K is large.

8. Conclusions

In this paper, we used a co-prime array system for spectrum sensing of multiband signals. Each channel of our proposed CA-based MWC system is equivalent to one channel of traditional MWC system. The performance of our system outperforms traditional MWC system because the noise among all channels in our method is uncorrelated with each other, the noise can be averaged to increase the SNR. In addition, due to the incorporation of co-prime arrays, a virtual ULA which has a much larger aperture can be produced. Such virtual ULA can be used to detect more targets than traditional ULA-based MWC. Simulation results demonstrate the performance of our methods. In this paper, we mainly focus on the demonstration of feasibility of our proposed method which is based on computer experiment. Our next work is to design the whole hardware device of CA-based MWC system.