# Spatial-Frequency Estimation for OFDM System with Coprime Array

Authors: Gang Shao; Wen-qin Wu; Lei Yin; Chen-yang Ding; Xiao-yuan Zhang
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Orthogonal frequency division multiplexing (OFDM) is the key technique of the communication system. In this paper, coprime array is designed in OFDM system to obtain joint direction of arrival (DOA) and carrier frequency offset (CFO) estimates. The proposed system can achieve beamforming and adaptive signal processing, reduce the multipath fading, and attain angle diversity. Meanwhile, coprime array has the advantages of extended array aperture, increased degrees of freedom, and reduced mutual coupling effect, which is adaptable in OFDM system. Moreover, this paper proposes a low-complexity parameter algorithm with superior performance, which first exploits propagator method (PM) as the initialization and then parallel factor (PARAFAC) method is employed for the accurate DOA and CFO estimation. Simulation results verify the effectiveness of the OFDM coprime array system and the proposed low-complexity multiparameter algorithm.


---

## Body

## 1. Introduction

Array signal processing is to arrange several sensors into sensor array to collect spatial signals and output discrete observation data related to signal sources [1]. Compared with the wireless communication model with only a single receiving antenna, antenna arrays can provide spatial angle diversity and are widely used in radar, sonar, seismic detection, and radio astronomy [2, 3]. Direction of arrival (DOA) estimation is an important processing step in the signal processing of array antenna or sensor systems. The process is to estimate the spatial position of the signal of interest based on the signal received by the antenna/sensor array. Since Howells proposed an adaptive sidelobe canceller [4], a large number of scholars have developed various DOA methods in the past few decades, such as traditional spectrum estimation and subspace spectrum estimation. The dimensions of DOA estimation have also expanded from one-dimensional DOA estimation to multiparameter joint estimation, such as typical elevation/azimuth two-dimensional DOA estimation, while emerging technologies such as compressed sensing [5, 6] and deep neural network [7] have been introduced to achieve superresolution channel estimation and DOA estimate.With the development of array signal processing, more and more attention has been paid to the modification of array manifold. The idea of the coprime linear array (CLA) is introduced in [8–10], and an equivalent virtual array with a larger aperture can be generated through the difference of multiple subarray data [11–17]. The spacing of these subarrays meets the coprime relationship. The coprime array proposed by Vaidyanathan and Pal is composed of two subarrays with M1 and M2 array elements, respectively. The distances between adjacent elements in the two subarrays are M2λ/2 and M1λ/2, respectively. The whole coprime array can obtain OM1M2 degrees of freedom by generating the virtual sensors. Reference [11] introduces compressed sensing into coprime array to achieve DOA estimation. First, a random compressed sensing core compresses the received signal into low-dimensional measurement values, and then, the compressed measurement values are utilized to achieve high-resolution DOA estimation as well as maintain the enlarged array aperture produced by the relatively prime array. A search-free DOA estimation method for coprime array is proposed in [12], which avoids spectral peak search and thus reduces the computational complexity significantly. Reference [13] investigates the problem of joint DOD and DOA with bistatic MIMO radar using tensor decomposition. Reference [14] proposes an improved coprime array model. Compared with the general coprime array, the number of elements required can be greatly reduced under the goal of achieving the same virtual aperture and degree of freedom. The array manifolds can be used in subspace classes or DOA estimation algorithms based on compressed sensing. Reference [15] studied the impact of bit quantization technology widely used in massive MIMO on the DOA estimation performance of nested arrays and coprime arrays and pointed out that the sparse arrays can yield better DOA estimation performance by exploiting 1-bit quantization.Orthogonal frequency division multiplexing (OFDM) has superior ability to resist multipath fading and can support multiuser access, which represents an efficient technique for mobile communications. However, because even a small frequency offset estimation error can cause interference between OFDM subcarriers [18, 19], resulting in a substantial decline in the performance of the OFDM system, accurate carrier frequency offset (CFO) estimation is critical to the practical application of OFDM systems. Therefore, CFO estimation has attracted the attention of a large number of researchers. Generally speaking, CFO estimation algorithms can be divided into nonblind CFO estimation algorithm and blind CFO algorithm. Nonblind CFO estimation algorithms include pilot-based CFO estimation [20] and training sequence-based CFO estimation [21, 22] and CFO estimation methods based on cyclic prefix [23]. Blind CFO estimation algorithms include CFO estimation based on MUSIC [24], CFO estimation based on ESPRIT [25], blind CFO estimation algorithm based on PM [26], and Capon [27] blind CFO estimation algorithm, kurtosis-based and constant modulus-based CFO estimation [28, 29], and cycle statistics-based CFO estimation [30] algorithm. Compared with the methods in [21, 22, 24], blind CFO estimation methods have higher bandwidth efficiency. Some subspace-based blind estimation methods such as MUSIC and ESPRIT can achieve high-precision CFO estimation. However, they require a virtual carrier and will bring a certain loss of spectral efficiency. For non-Gaussian sources, reference [28] uses kurtosis to measure non-Gaussianity, while references [29, 30] use other properties of the source, such as constant mode and cyclic stability. In addition, the trilinear decomposition technique [31] is also widely used in the field of CFO estimation and can achieve accurate estimation.However, according to the current researches, there are few effective algorithms for joint estimation of DOA and CFO in OFDM system. Considering the fact that the extended and sparse coprime array avoids compact array structure, where severe mutual coupling effects of half-wavelength-based array are alleviated, and the OFDM model offers the possibility of the joint estimation of DOA and CFO, we combine the advantages of coprime array and OFDM systems to achieve high-precision DOA and CFO joint estimation in this paper. In terms of algorithms, we propose a low-complexity joint parameter algorithm based on the traditional PM and trilinear decomposition parallel factor algorithm. PM is employed to obtain the coarse estimates as initialization, which possesses low complexity since there is no need to perform eigenvalue decomposition on the covariance matrix or singular value decomposition of received data.Considering that the conventional trilinear decomposition has a large amount of calculation due to the large number of iterations, the proposed algorithm employs the estimates obtained by PM as the initial iteration value, which reduces the number of iterations of the trilinear decomposition and thereby reduces the computational complexity remarkably.To sum up, the contributions of this paper are concluded as follows:(1)
We model the problem of spatial-frequency estimation for OFDM system with parallel factor (PARAFAC)(2)
We employ coprime array to obtain joint CFO and DOA estimates. Compared with uniform linear array (ULA), coprime array has extended array aperture and reduced mutual coupling effect(3)
A low-complexity multiparameter estimation algorithm is proposed to reduce the computational burden, where PM is utilized as the coarse initialization and PARAFAC algorithm as the accuracy estimates. Benefiting from the initialization, the number of iterations can be greatly reduced and the computational complexity can be substantially reduced


### 1.1. Notations

We use lowercase (uppercase) bold character to denote vector (matrix).⋅T is transpose of a matrix or vector. ⋅−1 denotes matrix inverse and ⋅+ denotes matrix pseudoinverse, respectively. ⊙ represents the Khatri-Rao product. diag⋅ symbolizes a diagonal matrix that uses the elements of the matrix as its diagonal element. ⋅F denotes Frobenius norm and tr⋅ represents the trace of the matrix, which is the sum of the elements on the main diagonal of the matrix.


## 1.1. Notations

We use lowercase (uppercase) bold character to denote vector (matrix).⋅T is transpose of a matrix or vector. ⋅−1 denotes matrix inverse and ⋅+ denotes matrix pseudoinverse, respectively. ⊙ represents the Khatri-Rao product. diag⋅ symbolizes a diagonal matrix that uses the elements of the matrix as its diagonal element. ⋅F denotes Frobenius norm and tr⋅ represents the trace of the matrix, which is the sum of the elements on the main diagonal of the matrix.


## 2. Preliminaries

### 2.1. Data Model

Assuming that each antenna is affected by the same carrier frequency offset, an OFDM system with an array of antennas is established. Among theN subcarriers of OFDM, P subcarriers are used for data communication, and the other N−P are virtual carriers. The length of the system cyclic prefix (CP) is L sampling intervals, where L is greater than the maximum delay spread. The output signal has been inserted into the cyclic prefix before being transmitted through the multipath fading channel. Define the baseband signal received at the reference point as
(1)uit=∑l=1Lmhildt−l,where hil,l=1,⋯,Lm represents the fading of the l-th path and Lm is the maximum multipath delay.After removing the cyclic prefix, the output signal of thei-th antenna can be expressed as [32]
(2)xik=EFPdiaghiskej2πΔfk−1N+L,where sk=s1k,s2k,⋯,sPkT represents the k-th idle and Δf denotes the carrier frequency offset. E=diag1,e2πΔf,⋯,ej2πN−1Δf symbolizes the carrier frequency offset matrix. hi=Hi0,Hi0,⋯,HiP−1T denotes the frequency domain channel response vector of the i-th receiving antenna. Hin=∑l=1Lihile−j2πnl/N is the frequency domain channel response of the n-th subcarrier at the i-th receiving antenna, where Li denotes the number of multipath channels and hil represents the l-th channel fading between the i-th receiving antenna and the transmitting antenna. Taking the channel fading of element at the origin h0l as the reference, the channel fading of the i-th element is hil=h0le−j2πdisinθi/λ, where di is the distance between i-th element and the reference, θ denotes the direction of arrival angle, and λ represents the wavelength. Without loss of generality, assuming that the channel parameters remain unchanged during K idle intervals, the source matrix can be defined as S=s1,s2,⋯,sKT∈ℂK×P.DefineXi=xi1,xi2,⋯,xiK, and Xi can be written as
(3)Xi=AdiaghiBT=ADiHBT,i=1,2,⋯,I,where B=diag1,ej2πΔfN+L,⋯,ej2πΔfK−1N+LS and Di. means to extract i rows from a matrix and form a diagonal matrix. A is a Vandermonde matrix and can be specified by
(4)A=EFP=11⋯1ej2π0/N+Δfej2π1/N+Δf⋯ej2πP−1/N+Δf⋮⋮⋱⋮ej2πN−10/N+Δfej2πN−11/N+Δf⋯ej2πN−1P−1/N+Δf∈ℂN×P.Then, the output signals of allI antennas can be expressed as
(5)X=X1X2⋮XI=AD1HAD2H⋮ADIHBT=H⊙ABT.In the presence of noise, the received signal model is(6)X~i=ADiHBT+Wi,i=1,2,⋯,I,where Wi represents the received noise of the i-th antenna.


### 2.2. Coprime Array Structure

Consider a pair of uniform linear arrays collocated at the origin, extended to the same direction withM1 and M2 sensor elements, respectively, where M1 and M2 are coprime integers. For the first linear subarray with M1 sensor elements, the interval between adjacent sensor elements is d1=M2λ/2, and for the second linear subarray with M2 sensor elements, the interval between adjacent sensor elements is d2=M1λ/2. The two subarrays share the first sensor element at the origin; hence, the total number of sensor element is I=M1+M2−1. Figure 1 depicts the typical array structure model of coprime array.Figure 1
Coprime array structure model.Compared with the traditional half-wavelength uniform linear array with spacing between adjacent elements of half wavelength, coprime array has extended array aperture and low mutual coupling effect. Parameter estimation performance is greatly improved due to the sparse array structure. The ambiguity problem can be solved according to the coprime characteristic.Proof.
Assume that there areK far-field signals impinging on the array with θk∈−90∘,90∘,k=1,2,⋯,K. Suppose there is an ambiguous angle θ^k except for θk, θk and θ^k can result in the same array manifold, which means the following relationship exists [33]:
(7)sinθk−sinθ^k=2PM2M2,sinθk−sinθ^k=2PM1M1,where M1=−M1−1,⋯,−1,1,⋯,M1−1 and M2=−M2−1,⋯,−1,1,⋯,M2−1. Then, there may have the equation that 2PM2/M2=2PM1/M1. In other words that is
(8)PM1M2=PM2M1.
SinceM1 and M2 are coprime numbers, no PM1,PM2 pair satisfies (8). Therefore, θ^k does not exist.☐


## 2.1. Data Model

Assuming that each antenna is affected by the same carrier frequency offset, an OFDM system with an array of antennas is established. Among theN subcarriers of OFDM, P subcarriers are used for data communication, and the other N−P are virtual carriers. The length of the system cyclic prefix (CP) is L sampling intervals, where L is greater than the maximum delay spread. The output signal has been inserted into the cyclic prefix before being transmitted through the multipath fading channel. Define the baseband signal received at the reference point as
(1)uit=∑l=1Lmhildt−l,where hil,l=1,⋯,Lm represents the fading of the l-th path and Lm is the maximum multipath delay.After removing the cyclic prefix, the output signal of thei-th antenna can be expressed as [32]
(2)xik=EFPdiaghiskej2πΔfk−1N+L,where sk=s1k,s2k,⋯,sPkT represents the k-th idle and Δf denotes the carrier frequency offset. E=diag1,e2πΔf,⋯,ej2πN−1Δf symbolizes the carrier frequency offset matrix. hi=Hi0,Hi0,⋯,HiP−1T denotes the frequency domain channel response vector of the i-th receiving antenna. Hin=∑l=1Lihile−j2πnl/N is the frequency domain channel response of the n-th subcarrier at the i-th receiving antenna, where Li denotes the number of multipath channels and hil represents the l-th channel fading between the i-th receiving antenna and the transmitting antenna. Taking the channel fading of element at the origin h0l as the reference, the channel fading of the i-th element is hil=h0le−j2πdisinθi/λ, where di is the distance between i-th element and the reference, θ denotes the direction of arrival angle, and λ represents the wavelength. Without loss of generality, assuming that the channel parameters remain unchanged during K idle intervals, the source matrix can be defined as S=s1,s2,⋯,sKT∈ℂK×P.DefineXi=xi1,xi2,⋯,xiK, and Xi can be written as
(3)Xi=AdiaghiBT=ADiHBT,i=1,2,⋯,I,where B=diag1,ej2πΔfN+L,⋯,ej2πΔfK−1N+LS and Di. means to extract i rows from a matrix and form a diagonal matrix. A is a Vandermonde matrix and can be specified by
(4)A=EFP=11⋯1ej2π0/N+Δfej2π1/N+Δf⋯ej2πP−1/N+Δf⋮⋮⋱⋮ej2πN−10/N+Δfej2πN−11/N+Δf⋯ej2πN−1P−1/N+Δf∈ℂN×P.Then, the output signals of allI antennas can be expressed as
(5)X=X1X2⋮XI=AD1HAD2H⋮ADIHBT=H⊙ABT.In the presence of noise, the received signal model is(6)X~i=ADiHBT+Wi,i=1,2,⋯,I,where Wi represents the received noise of the i-th antenna.


## 2.2. Coprime Array Structure

Consider a pair of uniform linear arrays collocated at the origin, extended to the same direction withM1 and M2 sensor elements, respectively, where M1 and M2 are coprime integers. For the first linear subarray with M1 sensor elements, the interval between adjacent sensor elements is d1=M2λ/2, and for the second linear subarray with M2 sensor elements, the interval between adjacent sensor elements is d2=M1λ/2. The two subarrays share the first sensor element at the origin; hence, the total number of sensor element is I=M1+M2−1. Figure 1 depicts the typical array structure model of coprime array.Figure 1
Coprime array structure model.Compared with the traditional half-wavelength uniform linear array with spacing between adjacent elements of half wavelength, coprime array has extended array aperture and low mutual coupling effect. Parameter estimation performance is greatly improved due to the sparse array structure. The ambiguity problem can be solved according to the coprime characteristic.Proof.
Assume that there areK far-field signals impinging on the array with θk∈−90∘,90∘,k=1,2,⋯,K. Suppose there is an ambiguous angle θ^k except for θk, θk and θ^k can result in the same array manifold, which means the following relationship exists [33]:
(7)sinθk−sinθ^k=2PM2M2,sinθk−sinθ^k=2PM1M1,where M1=−M1−1,⋯,−1,1,⋯,M1−1 and M2=−M2−1,⋯,−1,1,⋯,M2−1. Then, there may have the equation that 2PM2/M2=2PM1/M1. In other words that is
(8)PM1M2=PM2M1.
SinceM1 and M2 are coprime numbers, no PM1,PM2 pair satisfies (8). Therefore, θ^k does not exist.☐


## 3. Joint Parameter Estimation Algorithm

In this section, we first exploit PM to obtain the coarse initialization, and the PARAFAC is employed to get the accurate DOA and CFO estimates. Hence, the times of iterations can be reduced greatly.


### 3.1. PM

This section introduces the PM which is utilized as the initialization of CFO. In the case of noise, the received signal in (6) can be expressed as
(9)X~i=AS+Wi,i=1,2,⋯,I,where S=DiHBT,i=1,2,⋯,I and Wi denotes the noise received of i-th antenna. A is the direction matrix described in (4).It can be known from (4) that the direction matrix A is a Vandermonde matrix, and CFO estimation can be carried out by using the propagation method. The direction matrix A can be divided into
(10)A=A1A2,where A1∈CP×P is full rank, A2∈CN−P×P. There is a linear factor between the two matrices, denoted by
(11)A2=PcA1,where Pc is the propagator matrix.Define the matrixP∈CN×P as
(12)P=IPPc,where IP∈CP×P is an identity matrix. According to (10) and (12), the following expression can be obtained
(13)PA1=A1A2=A.The firstN−1 and last N−1 rows of P are represented by Pa∈CN−1×P and Pb∈CN−1×P, and the first N−1 and next N−1 rows of A are represented by Αa∈CN−1×P and Ab∈CN−1×P, respectively. According to Equation (13), we have
(14)PaPbA1=AaAb=AaAaΦr,where Φr=diag1,ej2πΔf+1/N,⋯,ej2πΔf+P−1/N∈CP×P. Subsequently, we can deduce the following relation:
(15)Pa+Pb=A1ΦrA1−1.DefineΨr=Pa+Pb∈CP×P. Since Ψr and Φr have the same eigenvalues, Φr can be obtained by eigen decomposition of Ψr, and then, the estimate of 2πΔf can be obtained:
(16)ej2πΔf=trΦr∑p=0P−1ej2πp/N.


### 3.2. PARAFAC Algorithm

#### 3.2.1. Trilinear Decomposition

One of the commonly used data detection methods for trilinear models is the trilinear alternating least squares (TALS) algorithm [11]. The basic ideas of the algorithm are as follows: (1) each residual matrix is updated based on least squares; (2) start to update another matrix; (3) repeat the above steps until the result converges. The TALS algorithm is described as follows. According to the signal model,
(17)X1X2⋮XI=AD1HAD2H⋮ADIHBT.The cost function in iterative decomposition is(18)minA,H,BX~1X~2⋮X~I−AD1HAD2H⋮ADIHBTF,where F represents the Frobenius norm. X~i,i=1,2,⋯,I are received signals with noise. Based on the least square criterion, the matrix B can be updated by
(19)B∧T=AD1H∧AD2H∧⋮ADIH∧+X~1X~2⋮X~I,where ⋅+ represents pseudoinverse in (19). H^ represents the estimate of the matrix H. Similarly, according to the second form of slicing in the parallel factor method Zm=BDmAHT,m=1,2,⋯,N, which is denoted as
(20)Z1Z2⋮ZN=BD1ABD2A⋮BDNAHT.Based on (20), we construct the least squares fitting as
(21)minF,H,SZ~1Z~2⋮Z~N−BD1ABD2A⋮BDNAHTF,where Z~m,m=1,2,⋯,N are the signal slices containing noise. The least squares of H can be updated to
(22)H∧T=B∧D1AB∧D2A⋮B∧DNA+Z~1Z~2⋮Z~N,where B^ represents the estimate of matrix B.The conditional least squares are used to update the matricesB and H until the algorithm converges, and the updated results are shown in Equations (19) and (22). The channel matrix H and source matrix B in frequency domain can be obtained according to the TALS algorithm. The TALS algorithm is optimal in the case of additive Gaussian noise [34]. The advantages of the algorithm include easy to implement, good convergence, and easy to extend to high order data. The main drawback is that the convergence process is sometimes slow [35]. In order to accelerate the convergence rate, the algorithm proposed in this paper obtains the estimation of matrix A according to (16).


#### 3.2.2. DOA Estimation

The channel is modeled as a finite impulse response (FIR) filter of lengthLm. The frequency domain channel vector Ht is processed through IFFT to get the time domain channel response vector, hi1,hi2,⋯,hiLm. According to the frequency domain channel matrix H, the time domain channel matrix Ht is shown as follows:
(23)Ht=h11h12⋯h1Lmh21h22⋯h2Lm⋮⋮⋱⋮hI1hI2⋯hILm,where hjl=hile−j2πj−idsinθl/λ.In order to get the direction matrixA, we need to normalize the matrix Ht. The direction vector of DOA θl is as follows.
(24)αθl=1,e−j2πdsinθl/λ,⋯,e−j2πI−1dsinθl/λT.According to (24), we can obtain
(25)g=−imaglnαθl=0,2πdsinθlλ,⋯,2πdI−1sinθlλT,where imag⋅ means to get its imaginary part of a complex number and ln⋅ represents the natural logarithm. Equation (25) should be arithmetic progression sequence. sinθl can be estimated by the least squares principle, and then, DOAs can be estimated according to sinθl.The frequency domain channel matrixH^ can be estimated according to the TALS algorithm. By performing IFFT processing of H^, the time domain channel response matrix H^t can be obtained. The direction matrix A^ can be obtained by normalizing the matrix H^t, which also solves the scale ambiguity problem. Suppose that the estimation result of the direction vector is α^θl (the l-th column of the estimated direction matrix A^), then according to (25), g^ can be obtained by processing α^θl. sinθl can be estimated according to the least squares principle.Least squares fitting result is(26)PC=g^,where
(27)P=11⋮102πdλ⋮I−12πdλ.The solution forC by least square is
(28)C=PTP−1PTg^.We defineD=PTP−1PT, (28) can be simplified as C=Dg^. Because the matrix D is a constant matrix in practical application, the computation of C can be greatly reduced.Finally, the DOA estimates are obtained by(29)θ^l=sin−1c1.


#### 3.2.3. CFO Estimation

According to (4), the p-th column of A is apΔf=1,ej2πp−1/N+Δf,⋯,ej2πN−1p−1/N+ΔfT, and from that, we define
(30)gCFO=imaglnapΔf.The specific expression ofgCFO is as follows:
(31)gCFO=0,2πp−1N+Δf,⋯,2πN−1p−1N+ΔfT=p−1N+Δfq,where q=0,2π,⋯,2πN−1T, CFO Δf can be estimated based on the LS criterion. The scale ambiguity problem can be solved by dividing each element of a^l (the l-th column of matrix A^) by a^l,1 (the first member of vector a^l). According to (31), g^CFO can be obtained. Then, we use the LS criterion to measure the CFO.Define LS fitting asQc=g^CFO, where Q=1N×1,q, c=c0,fpT, and fp is an estimate of p−1/N+Δf,p∈1,2,⋯,P. Subsequently, we can obtain the LS solution for c:
(32)c^0f^p=Q+g^CFO.Then, we can estimatef^p,p=1,2,⋯,P. Subsequently, the CFO estimation can be achieved by
(33)Δf^=1P∑p=1Pf^p−∑p=1PP−1N.


### 3.3. Performance Analysis

In this paper, a PM-PARAFAC algorithm is proposed for OFDM systems with array antenna, by which CFO and DOA can be estimated jointly. In this section, we provide the main steps of the proposed algorithm, complexity analysis, and the advantages of the algorithm.The detailed steps are shown as follows.Step 1.
According to the received signalX~i, the initial value of CFO estimation Δf^ini is obtained by PM.Step 2.
According to (16), the initial estimation of matrix A is obtained.Step 3.
The source matrixB and the frequency domain channel matrix H are initialized, respectively.Step 4.
According to (19), the source matrix B is updated by LS.Step 5.
According to (22), the frequency domain channel matrix H is updated by LS.Step 6.
Repeat Step4 to Step 5 until convergence.Step 7.
Estimate DOA and CFO estimates by frequency domain channel matrixH and matrixA, respectively.The calculation complexity of the proposed algorithm and the other three algorithms to be compared is shown in Table1, where N and P are the number of subcarriers and channels employed by the OFDM system and K denotes the number of snapshots.Table 1
Complexity calculation of each algorithm.
AlgorithmComplexityPMOKN2+PNP+N+3P2P−1+P3ESPRITOKN2+N3+3P2N−1+P3PARAFACO3P3+3PNIK+P2IN+IK+NK+I+N+KCOMFACO3P3+3PN′I′K′+P2I′N′+I′K′+N′K′+I′+N′+K′PM-PARAFACOKN2+N3+3P2N−1+P3+3P3+3PN″I″K″+P2I″N″+I″K″+N″K″+I″+N″+K″N′<N, I′<I, K′<K, N′′<N, I′′<I, and K′′<K.Finally, we summarize the merits of the proposed PM-PARAFAC algorithm as follows:(a)
The proposed algorithm initializes the CFO estimation by PM. It needs no eigen decomposition which takes a great amount of calculation burden(b)
Joint DOA and CFO estimates can be obtained by the PARAFAC algorithm. The estimation accuracy is relatively high and the initialization reduces the number of iterations(c)
The principle of coprime array keeps the ambiguity-free characteristic, which requires no extra pairing steps and avoids paring failure


## 3.1. PM

This section introduces the PM which is utilized as the initialization of CFO. In the case of noise, the received signal in (6) can be expressed as
(9)X~i=AS+Wi,i=1,2,⋯,I,where S=DiHBT,i=1,2,⋯,I and Wi denotes the noise received of i-th antenna. A is the direction matrix described in (4).It can be known from (4) that the direction matrix A is a Vandermonde matrix, and CFO estimation can be carried out by using the propagation method. The direction matrix A can be divided into
(10)A=A1A2,where A1∈CP×P is full rank, A2∈CN−P×P. There is a linear factor between the two matrices, denoted by
(11)A2=PcA1,where Pc is the propagator matrix.Define the matrixP∈CN×P as
(12)P=IPPc,where IP∈CP×P is an identity matrix. According to (10) and (12), the following expression can be obtained
(13)PA1=A1A2=A.The firstN−1 and last N−1 rows of P are represented by Pa∈CN−1×P and Pb∈CN−1×P, and the first N−1 and next N−1 rows of A are represented by Αa∈CN−1×P and Ab∈CN−1×P, respectively. According to Equation (13), we have
(14)PaPbA1=AaAb=AaAaΦr,where Φr=diag1,ej2πΔf+1/N,⋯,ej2πΔf+P−1/N∈CP×P. Subsequently, we can deduce the following relation:
(15)Pa+Pb=A1ΦrA1−1.DefineΨr=Pa+Pb∈CP×P. Since Ψr and Φr have the same eigenvalues, Φr can be obtained by eigen decomposition of Ψr, and then, the estimate of 2πΔf can be obtained:
(16)ej2πΔf=trΦr∑p=0P−1ej2πp/N.


## 3.2. PARAFAC Algorithm

### 3.2.1. Trilinear Decomposition

One of the commonly used data detection methods for trilinear models is the trilinear alternating least squares (TALS) algorithm [11]. The basic ideas of the algorithm are as follows: (1) each residual matrix is updated based on least squares; (2) start to update another matrix; (3) repeat the above steps until the result converges. The TALS algorithm is described as follows. According to the signal model,
(17)X1X2⋮XI=AD1HAD2H⋮ADIHBT.The cost function in iterative decomposition is(18)minA,H,BX~1X~2⋮X~I−AD1HAD2H⋮ADIHBTF,where F represents the Frobenius norm. X~i,i=1,2,⋯,I are received signals with noise. Based on the least square criterion, the matrix B can be updated by
(19)B∧T=AD1H∧AD2H∧⋮ADIH∧+X~1X~2⋮X~I,where ⋅+ represents pseudoinverse in (19). H^ represents the estimate of the matrix H. Similarly, according to the second form of slicing in the parallel factor method Zm=BDmAHT,m=1,2,⋯,N, which is denoted as
(20)Z1Z2⋮ZN=BD1ABD2A⋮BDNAHT.Based on (20), we construct the least squares fitting as
(21)minF,H,SZ~1Z~2⋮Z~N−BD1ABD2A⋮BDNAHTF,where Z~m,m=1,2,⋯,N are the signal slices containing noise. The least squares of H can be updated to
(22)H∧T=B∧D1AB∧D2A⋮B∧DNA+Z~1Z~2⋮Z~N,where B^ represents the estimate of matrix B.The conditional least squares are used to update the matricesB and H until the algorithm converges, and the updated results are shown in Equations (19) and (22). The channel matrix H and source matrix B in frequency domain can be obtained according to the TALS algorithm. The TALS algorithm is optimal in the case of additive Gaussian noise [34]. The advantages of the algorithm include easy to implement, good convergence, and easy to extend to high order data. The main drawback is that the convergence process is sometimes slow [35]. In order to accelerate the convergence rate, the algorithm proposed in this paper obtains the estimation of matrix A according to (16).


### 3.2.2. DOA Estimation

The channel is modeled as a finite impulse response (FIR) filter of lengthLm. The frequency domain channel vector Ht is processed through IFFT to get the time domain channel response vector, hi1,hi2,⋯,hiLm. According to the frequency domain channel matrix H, the time domain channel matrix Ht is shown as follows:
(23)Ht=h11h12⋯h1Lmh21h22⋯h2Lm⋮⋮⋱⋮hI1hI2⋯hILm,where hjl=hile−j2πj−idsinθl/λ.In order to get the direction matrixA, we need to normalize the matrix Ht. The direction vector of DOA θl is as follows.
(24)αθl=1,e−j2πdsinθl/λ,⋯,e−j2πI−1dsinθl/λT.According to (24), we can obtain
(25)g=−imaglnαθl=0,2πdsinθlλ,⋯,2πdI−1sinθlλT,where imag⋅ means to get its imaginary part of a complex number and ln⋅ represents the natural logarithm. Equation (25) should be arithmetic progression sequence. sinθl can be estimated by the least squares principle, and then, DOAs can be estimated according to sinθl.The frequency domain channel matrixH^ can be estimated according to the TALS algorithm. By performing IFFT processing of H^, the time domain channel response matrix H^t can be obtained. The direction matrix A^ can be obtained by normalizing the matrix H^t, which also solves the scale ambiguity problem. Suppose that the estimation result of the direction vector is α^θl (the l-th column of the estimated direction matrix A^), then according to (25), g^ can be obtained by processing α^θl. sinθl can be estimated according to the least squares principle.Least squares fitting result is(26)PC=g^,where
(27)P=11⋮102πdλ⋮I−12πdλ.The solution forC by least square is
(28)C=PTP−1PTg^.We defineD=PTP−1PT, (28) can be simplified as C=Dg^. Because the matrix D is a constant matrix in practical application, the computation of C can be greatly reduced.Finally, the DOA estimates are obtained by(29)θ^l=sin−1c1.


### 3.2.3. CFO Estimation

According to (4), the p-th column of A is apΔf=1,ej2πp−1/N+Δf,⋯,ej2πN−1p−1/N+ΔfT, and from that, we define
(30)gCFO=imaglnapΔf.The specific expression ofgCFO is as follows:
(31)gCFO=0,2πp−1N+Δf,⋯,2πN−1p−1N+ΔfT=p−1N+Δfq,where q=0,2π,⋯,2πN−1T, CFO Δf can be estimated based on the LS criterion. The scale ambiguity problem can be solved by dividing each element of a^l (the l-th column of matrix A^) by a^l,1 (the first member of vector a^l). According to (31), g^CFO can be obtained. Then, we use the LS criterion to measure the CFO.Define LS fitting asQc=g^CFO, where Q=1N×1,q, c=c0,fpT, and fp is an estimate of p−1/N+Δf,p∈1,2,⋯,P. Subsequently, we can obtain the LS solution for c:
(32)c^0f^p=Q+g^CFO.Then, we can estimatef^p,p=1,2,⋯,P. Subsequently, the CFO estimation can be achieved by
(33)Δf^=1P∑p=1Pf^p−∑p=1PP−1N.


## 3.2.1. Trilinear Decomposition

One of the commonly used data detection methods for trilinear models is the trilinear alternating least squares (TALS) algorithm [11]. The basic ideas of the algorithm are as follows: (1) each residual matrix is updated based on least squares; (2) start to update another matrix; (3) repeat the above steps until the result converges. The TALS algorithm is described as follows. According to the signal model,
(17)X1X2⋮XI=AD1HAD2H⋮ADIHBT.The cost function in iterative decomposition is(18)minA,H,BX~1X~2⋮X~I−AD1HAD2H⋮ADIHBTF,where F represents the Frobenius norm. X~i,i=1,2,⋯,I are received signals with noise. Based on the least square criterion, the matrix B can be updated by
(19)B∧T=AD1H∧AD2H∧⋮ADIH∧+X~1X~2⋮X~I,where ⋅+ represents pseudoinverse in (19). H^ represents the estimate of the matrix H. Similarly, according to the second form of slicing in the parallel factor method Zm=BDmAHT,m=1,2,⋯,N, which is denoted as
(20)Z1Z2⋮ZN=BD1ABD2A⋮BDNAHT.Based on (20), we construct the least squares fitting as
(21)minF,H,SZ~1Z~2⋮Z~N−BD1ABD2A⋮BDNAHTF,where Z~m,m=1,2,⋯,N are the signal slices containing noise. The least squares of H can be updated to
(22)H∧T=B∧D1AB∧D2A⋮B∧DNA+Z~1Z~2⋮Z~N,where B^ represents the estimate of matrix B.The conditional least squares are used to update the matricesB and H until the algorithm converges, and the updated results are shown in Equations (19) and (22). The channel matrix H and source matrix B in frequency domain can be obtained according to the TALS algorithm. The TALS algorithm is optimal in the case of additive Gaussian noise [34]. The advantages of the algorithm include easy to implement, good convergence, and easy to extend to high order data. The main drawback is that the convergence process is sometimes slow [35]. In order to accelerate the convergence rate, the algorithm proposed in this paper obtains the estimation of matrix A according to (16).


## 3.2.2. DOA Estimation

The channel is modeled as a finite impulse response (FIR) filter of lengthLm. The frequency domain channel vector Ht is processed through IFFT to get the time domain channel response vector, hi1,hi2,⋯,hiLm. According to the frequency domain channel matrix H, the time domain channel matrix Ht is shown as follows:
(23)Ht=h11h12⋯h1Lmh21h22⋯h2Lm⋮⋮⋱⋮hI1hI2⋯hILm,where hjl=hile−j2πj−idsinθl/λ.In order to get the direction matrixA, we need to normalize the matrix Ht. The direction vector of DOA θl is as follows.
(24)αθl=1,e−j2πdsinθl/λ,⋯,e−j2πI−1dsinθl/λT.According to (24), we can obtain
(25)g=−imaglnαθl=0,2πdsinθlλ,⋯,2πdI−1sinθlλT,where imag⋅ means to get its imaginary part of a complex number and ln⋅ represents the natural logarithm. Equation (25) should be arithmetic progression sequence. sinθl can be estimated by the least squares principle, and then, DOAs can be estimated according to sinθl.The frequency domain channel matrixH^ can be estimated according to the TALS algorithm. By performing IFFT processing of H^, the time domain channel response matrix H^t can be obtained. The direction matrix A^ can be obtained by normalizing the matrix H^t, which also solves the scale ambiguity problem. Suppose that the estimation result of the direction vector is α^θl (the l-th column of the estimated direction matrix A^), then according to (25), g^ can be obtained by processing α^θl. sinθl can be estimated according to the least squares principle.Least squares fitting result is(26)PC=g^,where
(27)P=11⋮102πdλ⋮I−12πdλ.The solution forC by least square is
(28)C=PTP−1PTg^.We defineD=PTP−1PT, (28) can be simplified as C=Dg^. Because the matrix D is a constant matrix in practical application, the computation of C can be greatly reduced.Finally, the DOA estimates are obtained by(29)θ^l=sin−1c1.


## 3.2.3. CFO Estimation

According to (4), the p-th column of A is apΔf=1,ej2πp−1/N+Δf,⋯,ej2πN−1p−1/N+ΔfT, and from that, we define
(30)gCFO=imaglnapΔf.The specific expression ofgCFO is as follows:
(31)gCFO=0,2πp−1N+Δf,⋯,2πN−1p−1N+ΔfT=p−1N+Δfq,where q=0,2π,⋯,2πN−1T, CFO Δf can be estimated based on the LS criterion. The scale ambiguity problem can be solved by dividing each element of a^l (the l-th column of matrix A^) by a^l,1 (the first member of vector a^l). According to (31), g^CFO can be obtained. Then, we use the LS criterion to measure the CFO.Define LS fitting asQc=g^CFO, where Q=1N×1,q, c=c0,fpT, and fp is an estimate of p−1/N+Δf,p∈1,2,⋯,P. Subsequently, we can obtain the LS solution for c:
(32)c^0f^p=Q+g^CFO.Then, we can estimatef^p,p=1,2,⋯,P. Subsequently, the CFO estimation can be achieved by
(33)Δf^=1P∑p=1Pf^p−∑p=1PP−1N.


## 3.3. Performance Analysis

In this paper, a PM-PARAFAC algorithm is proposed for OFDM systems with array antenna, by which CFO and DOA can be estimated jointly. In this section, we provide the main steps of the proposed algorithm, complexity analysis, and the advantages of the algorithm.The detailed steps are shown as follows.Step 1.
According to the received signalX~i, the initial value of CFO estimation Δf^ini is obtained by PM.Step 2.
According to (16), the initial estimation of matrix A is obtained.Step 3.
The source matrixB and the frequency domain channel matrix H are initialized, respectively.Step 4.
According to (19), the source matrix B is updated by LS.Step 5.
According to (22), the frequency domain channel matrix H is updated by LS.Step 6.
Repeat Step4 to Step 5 until convergence.Step 7.
Estimate DOA and CFO estimates by frequency domain channel matrixH and matrixA, respectively.The calculation complexity of the proposed algorithm and the other three algorithms to be compared is shown in Table1, where N and P are the number of subcarriers and channels employed by the OFDM system and K denotes the number of snapshots.Table 1
Complexity calculation of each algorithm.
AlgorithmComplexityPMOKN2+PNP+N+3P2P−1+P3ESPRITOKN2+N3+3P2N−1+P3PARAFACO3P3+3PNIK+P2IN+IK+NK+I+N+KCOMFACO3P3+3PN′I′K′+P2I′N′+I′K′+N′K′+I′+N′+K′PM-PARAFACOKN2+N3+3P2N−1+P3+3P3+3PN″I″K″+P2I″N″+I″K″+N″K″+I″+N″+K″N′<N, I′<I, K′<K, N′′<N, I′′<I, and K′′<K.Finally, we summarize the merits of the proposed PM-PARAFAC algorithm as follows:(a)
The proposed algorithm initializes the CFO estimation by PM. It needs no eigen decomposition which takes a great amount of calculation burden(b)
Joint DOA and CFO estimates can be obtained by the PARAFAC algorithm. The estimation accuracy is relatively high and the initialization reduces the number of iterations(c)
The principle of coprime array keeps the ambiguity-free characteristic, which requires no extra pairing steps and avoids paring failure


## 4. Simulation Results

We perform some simulations to confirm the effectiveness of the proposed PM-PARAFAC algorithm in the section. First, we compare the proposed algorithm with PM and ESPRIT algorithms to highlight the prominent performance of the proposed algorithm. Then, we demonstrate the effectiveness of CFO estimation by comparing the proposed algorithm with COMFAC and PARAFAC algorithms. Finally, we compare the DOA estimation performance of the proposed algorithm with uniform linear array and coprime array. In the following simulations, root mean square error (RMSE) is employed to evaluate the performance, which is defined as(34)RMSEa=1K∑k=1K1I∑i=1Iθ∧k,i−θk2,RMSEΔf=1K∑k=1K1I∑i=1IΔf∧k,i−Δf∧k2,where θ^k,i and Δf^k,i are the estimated values of θk and Δfk during the i‐th simulation and I is the number of independent simulations.


### 4.1. CFO Performance Estimation with Different Algorithms

We simulate the RMSE of CFO with the proposed CLA and the proposed algorithm in Figure2, where θ,Δf=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. Compared with the PM, the proposed algorithm improves the estimation accuracy by the cascading PARAFAC algorithm with the initial estimation of CFO.Figure 2
CFO performance comparison with different algorithms.


### 4.2. DOA Performance Estimation with Different Algorithms

In Figure3, we verify the effectiveness of DOA estimation with the proposed algorithm, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. As depicted in Figure 3, the proposed PM-PARAFAC algorithm outperforms the other three algorithms. Neither PARAFAC nor COMFAC alone can reach the performance of this algorithm because it does not narrow the iteration range through initialization.Figure 3
DOA performance comparison with different algorithms.


### 4.3. CFO Performance Estimation with Different Arrays

Figure4 demonstrates the CFO performance of the proposed algorithm with uniform linear array and coprime array, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. Coprime array shows better performance in RMSE due to its extended array aperture and large degrees of freedom.Figure 4
CFO performance comparison with different arrays.


### 4.4. DOA Performance Estimation with Different Arrays

Figure5 depicts the RMSE of DOA estimation of the proposed algorithm with uniform linear array and coprime array, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. It is indicated explicitly that the coprime array outperforms the uniform linear array in RMSE.Figure 5
DOA performance comparison with different arrays.


## 4.1. CFO Performance Estimation with Different Algorithms

We simulate the RMSE of CFO with the proposed CLA and the proposed algorithm in Figure2, where θ,Δf=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. Compared with the PM, the proposed algorithm improves the estimation accuracy by the cascading PARAFAC algorithm with the initial estimation of CFO.Figure 2
CFO performance comparison with different algorithms.


## 4.2. DOA Performance Estimation with Different Algorithms

In Figure3, we verify the effectiveness of DOA estimation with the proposed algorithm, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. As depicted in Figure 3, the proposed PM-PARAFAC algorithm outperforms the other three algorithms. Neither PARAFAC nor COMFAC alone can reach the performance of this algorithm because it does not narrow the iteration range through initialization.Figure 3
DOA performance comparison with different algorithms.


## 4.3. CFO Performance Estimation with Different Arrays

Figure4 demonstrates the CFO performance of the proposed algorithm with uniform linear array and coprime array, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. Coprime array shows better performance in RMSE due to its extended array aperture and large degrees of freedom.Figure 4
CFO performance comparison with different arrays.


## 4.4. DOA Performance Estimation with Different Arrays

Figure5 depicts the RMSE of DOA estimation of the proposed algorithm with uniform linear array and coprime array, where θ1,Δf1=10∘,0.3∘, θ2,Δf2=20∘,0.3∘, M1=8, M2=7, and SNR=0,20dB. It is indicated explicitly that the coprime array outperforms the uniform linear array in RMSE.Figure 5
DOA performance comparison with different arrays.


## 5. Conclusions

In this paper, we apply OFDM system to the coprime array to obtain the joint estimation of DOA and CFO. The coprime array structure extends the array aperture, reduces mutual coupling effect, and has the ability of achieving ambiguity-free estimates. Meanwhile, the proposed PM-PARAFAC algorithm uses the CFO initialization estimation of the PM algorithm to obtain high-precision estimates while easing the computational burden. Both theoretical analysis and simulation results substantiate the superiority of this algorithm.



---

*Source: 9998342-2021-05-27.xml*