DOA and Polarization Estimation for Non-Circular Signals in 3-D Millimeter Wave Polarized Massive MIMO Systems

1 DO A and Polarization Estimation for Non-Circular Signals in 3-D Millimete r W a v e Polarized Massi ve MIMO Systems Liangtian W an, Member , I EEE, Kaihui Liu, Y ing-Cha ng Liang F ellow , IEEE, and T on g Zhu . Abstract —In this paper , an algorithm of multiple signal classiﬁ- cation (MUSIC) is proposed fo r two-dimensional (2-D) direction- of-arriv al (DO A) and polarization estimation of non-circular signal in three-dimensional (3 -D) millimeter wa ve polarized large- scale/massiv e multipl e-input-multi ple-output (MIMO) systems. The traditional MUSIC-based a lgorithms can estimate either the DOA and p olarization for circular signal or the DOA fo r n on-circular signal by usin g sp ectrum search. By contrast, in the proposed algorithm only the DO A estimation needs spectrum search, and the polarization estimation has a closed- fo rm expression. First, a n ov el dimension-reduced M USIC (DR- MUSIC) is propose d f or DO A and polarization estimation of circular signal with l ow computational complexity . Next, based on the qu aternion th eory , a novel algorithm named q uaternion non-circular M USIC ( QNC- M USIC) is proposed for parameter estimation of n on-circular signal with high esti mation accuracy . Then based on the DO A estimation result u sing QNC-MU S IC, the polarization estimation of non-circular signal i s acqui red b y using the closed-form expression of th e polarization estimation in DR-MUSIC. In addition, t h e computational complexity analysis shows that compared with the conv entional DOA and p olarization estimation algorithms, our proposed QNC-MUSIC and DR- MUSIC h a ve much lower computational complexity , especially when the source number is large. The stochastic Cram ´ er -Rao Bound (CRB) for the estimation of th e 2-D DO A and polariza- tion parameters of the non-circular signals is deriv ed as well. Finally , numerical examples are p ro vided to demonstrate that the proposed algorithms can improv e the parameter estimation perfo rmance wh en th e large-scale/massi ve MIMO systems are employed. Index T erms —Direction-of-arriv al (DO A) and p olarization es- timation, polarized large-scale/ma ssiv e multiple-input-mul tiple- output ( L S-MIMO/massive M IMO), three-dimensional (3-D) mil- limeter wav e communication, circular and non-circular signals, quaternion. I . I N T RO D U C T I O N Driv en by the p roliferation of more s ophisticated smart- phone and social media, the wireless mobile trafﬁc will L. W an is with the K ey Laboratory for Ubiquitous Network and Service Softwa re of Liaoning Province, School of Software, Dalian Uni versity of T echnol ogy , Dalian 116620, China (e-mail: wanliangt ian@dlut.edu.cn). K. Liu is with National Key Laborato ry of Scienc e and T echnology on Communicat ions, Univ ersity of Electronic S cience and T echnology of China, Chengdu 611731, China (e-mail: kaihuil@std.uestc . edu.cn). Y .-C. Liang is with National Key Laboratory of Science and T echnology on Communications, Univ ersity of E lectro nic Science and T echnolog y of China, Chengdu 611731, China , and also with School of Elect rical and Information Engineering, Uni versity of Sydney , NSW 2006, Australia (e-mail: liangyc @ieee.org). T . Zhu is with Ti anjin Institute of Computing T echnol ogy , 300000, Ti anjin, China (e-mail: zhutonghe u@gmail.com). continue to grow at an expo nential p a c e, thus the capacity of cellular d ata n e tworks needs to in crease in o rders of magnitu de [ 1 ]. Millimeter w av e commu nication is a very pr omising approa c h for meetin g this challenge because of two reasons. First, there is a huge amount of a vailable spectrum ranging from 30 GHz to 300 GHz, which is mu ch mo re than those used by the existing wireless commun ication system s. Secon d , due to the small carrier wavelength of millimeter w av e , a large numb er of anten na elements can be arran ged in highly directive steerable arrays [ 2 ]. Thus the energy efﬁciency can be improved d ramatically when more en ergy is con c entrated in a particular direction. Furtherm ore, the highly d ir ectiv e steerable array m entioned above, are also kn own as massiv e mu ltiple- input-mu ltiple-outp u t (MIMO), can achieve extremely high spectrum ef ﬁciency , thus is the key enabling tech nology fo r gigabit-p er-second data transmission in millime te r wav e com- munication [ 3 ]. W ith the spatial free d om of fered by large antenna arrays, abundant mo bile terminals are expected to occupy the same set of tim e and frequen cy resources with negligible interfer e n ce [ 4 ], [ 5 ]. Extensive research has been conducted fo r massive MIMO systems, such as inter ference mitigation [ 6 ], multiuser be a m- forming [ 7 ] and join t spatial division and multiplexing [ 8 ]. Howe ver, in all the research works mentio ned ab ove, knowl- edge of channe l corr elations at the base stations (BSs) is required [ 5 ]. In order to model chan nel corr elations, the geometric sto c hastic chan nel model is widely used [ 6 ], [ 8 ]– [ 11 ], wh erein d irection-of -arriv als ( D OAs) of sign a l pa th s are crucial m odel p arameters. Thu s, accur ate DO A e stima tio n fo r dominan t sig nal paths is a prerequ isite for chann el cor r elation acquisitions in millimeter wave commun ic a tio n. In wireless com m unication systems, modulated signals based on binary ph ase sh ift ke ying (BPSK) and ampli- tude modu la tio n (AM) have been wide ly u sed . Dif ferent from quad r ature amp litude modulatio n (QAM) and quadr a- ture phase sh ift keying (QPSK) signals, the afo remention ed modulated signals are non-circ u lar in th e sense th at their uncon ju gated c ovariance m atrices are not equal to zero. This distinctiv e cha racteristic can be utilized to improve the par am- eter estimation perform ance. A. Rela ted W ork For DO A e stimation o f non-cir cular signals, the so -called non-cir cular m u ltiple signal c la ssiﬁcation (NC-MUSIC) algo- rithm was p roposed in [ 1 2 ], [ 13 ]. T he secon d-orde r asym ptot- 2 ically minimum variance (AMV) algorithms were pro p osed in [ 14 ], in which a closed-form expression o f the lower bound on the asymptotic covariance of estimatio ns giv en by arbitrary second-o rder algorithms was ev alu ated. Howe ver , the computatio nal c o mplexity of these algo rithms is tremendous because of the mu ltid imensional n onlinear optimizatio n [ 14 ]. T o reduce the co m putational comp lexity , the ro ot-NC-MUSIC algorithm was pr oposed in [ 15 ]. Based on a determin ant-based method, the DOAs of non - circular an d circular signals were simultaneou sly estimated in [ 1 6 ]. Howev e r, when two types of signals are too clo se, the DO A estimation pe rforman ce may degrade. By exploiting the non-circu larity , the DO As of no n- circular an d circular sign als were separately estima ted in [ 17 ]. By exploiting th e stro nger orthog onality in the biquater nion domain, the biquatern ion cumu la n t-MUSIC has b e en prop osed for DOA estimation [ 18 ]. Recently , the sparse r epresentation based m ethod has be e n propo sed in [ 19 ] with high e stimation accuracy and resolution. However , the com p utational complex- ity of this meth od is much larger than that of th e subspace- based meth ods. The one-d imensional (1 -D) NC Stan dard e stimation of sig- nal p arameters v ia rotational inv ariance techniq ues (ESPRIT) and two-dimensional (2 -D) NC Unitary ESPRIT h av e been propo sed in [ 20 ] and [ 2 1 ], r espectiv ely , for DO A estima- tion of non- c ircular signals. Recently , R multid im ensional ESPRIT -typ e algorith ms have been p roposed in [ 22 ], and the perturb ation analy sis of tensor -ESPRIT -ty p e alg orithms have been presented as well. Based on this, R multidimen sional ESPRIT -typ e algorithm s h ave b een applied to estima te strictly second-o rder non-circu lar sources [ 23 ], which is r egarded as an extension of methods in [ 20 ] and [ 21 ]. I n addition , in [ 23 ], the performance of th ese ESPRIT -type alg orithms has been analyzed as well. Recen tly , two ESPRIT -based algorith m s, termed CNC Standard ESPRIT an d CNC Unitary ESPRIT , have been d evised in [ 24 ] un der coexistence of circu lar and strictly non-circular signals b ased on NC-ESPRIT metho ds [ 25 ]. Th ey yield closed -form estimation w ith low compu ta- tional co mplexity . Howe ver , these methods [ 20 ]–[ 24 ] can not be used for polarization parameter estimatio ns. For DOA an d polarization estimation , th e E SPRIT algorithm has bee n used in the polariza tio n sensitiv e array . The polar- ization character istic of the signal and the relative inv ar iance between th e ortho gonal d ipole and the mag netic output have been explo ited f o r the p a r ameter estimation with unifo rm linear array (ULA) [ 26 ]. The root-M U SI C algorith m has been propo sed based on the diversely polariz e d c haracteristic [ 27 ]. In the case when the array m anifold is p artly known, the fourth - order statistics-based m ethod has been presented for joint para meter estimation [ 28 ]. Then 2 q th- order, q ≥ 2 , MUSIC methods h a ve been applied to array s having diversely polarized antennas for diversely polarized sources [ 29 ]. The parallel factor (P ARAF A C) an alysis (low or der ten so r) has been used for estimating DO A and polarizatio n par a meters, and the con ventional complex m atrix model is replaced b y the low or der tensor model [ 30 ]. The orthog onality amon g propag ation direction of electromagn etic wav e, electric ﬁeld and m a gnetic ﬁeld is reﬂected proﬁtab ly [ 31 ]. Based on the effective apertu re distribution function , an extension of root-MUSIC algo rithm was p roposed f o r DOA and po lar- ization estimation with arbitrary array co nﬁguration s [ 32 ]. The sparse r epresentation based me th od has b een p roposed in [ 33 ] for ULA by solv in g a weig hted grou p la sso problem in secon d-ord e r statistics d o main. Howe ver, the comp utational complexity of this m ethod is m uch larger than that of the subspace-b a sed metho ds. The q uaternion -MUSIC algo rithm has been proposed in [ 34 ], and a comparison between long vector ortho gonality and quaternion vector orthogonality is also per formed. Th e biquater n ion matr ix d iagonalization has been used f or DO A estimation b ased on vecto r-antenn as [ 35 ], [ 36 ]. However , these algorithm s assume th at the signal is circular, the e stima tio n perform ance of non-circu lar signal cannot b e improved any fu rther . T o the best of our knowledge, no con tributions have dealt yet with DO A and polarization fo r non-cir cular sign als. B. Movitation Recently , there has b een a grad ual d emand for the use of polarized an tenna systems, especially for 5G mo bile commu- nication systems [ 37 ], [ 38 ]. Th is is because of the fact that, for th e design of space- limited wireless devices, the antenn a polarization is a crucial resource to be exploited. T he degree- of-free dom and multip lexing could be incr eased by exploiting the antenn a polarization . In addition , a massiv e MIMO system equippe d with electroma g netic vetor sensors (EMVSs) cou ld generally form a uniform rectan g ular arr ay (URA). It sh o uld be noted that th is URA cou ld estimate not o nly the DOA of the incident signal, but also its po larization. The BS could use po - larization parame ters to disting uish d ifferent mobile terminals, since tho se p a rameters sh o uld contain un ique id entiﬁcation of mobile term inals. In a secure millimeter wave com m unication, the po larization parameter s can be used for encrypting the classiﬁed informatio n, only the po larized massi ve MIMO sys- tems could decode this encryptio n in formation . There should be other applications that the polar ization p arameters could b e used fo r . Thus, th e polarizatio n parameters’ estimation using polarized massi ve MIMO systems is a meaningful resear c h ﬁeld in millimeter wave commun ication as we ll. In this pap er , we a d opt the pola r ized massive MIMO sys- tems to estimate th e 2-D DO A and polarization of m ultiple no-circu lar source s, since ac curate DO A an d polarization estimations are particu larly cr itica l for chan nel cor relation acquisitions, as well a s for the m o bile terminal identiﬁcatio n and the infor mation security mentioned above in millimeter wa ve co mmunicatio n. C. Contribution In this paper, a MUSIC-b a sed algorithm is proposed for 2-D DO A an d p olarization estimation of multiple no-circu lar sources in polarized massiv e MIMO sy stem s em ploying very large URAs. The circu lar signal model contain ing DO A and polarization par ameter are constructed fo r po larized massive MIMO system s, and th e n on-circu la r signal model is con- structed based on quaternio n theory . The partial deri vativ e of the spectr um function is utilized to reduce th e dimension of par a meter search. The DOA p arameter is estimated at ﬁrst, a n d the po larization par ameter is estimated based on the 3 results of the DO A par ameter . T o be mo re spec iﬁc, the main contributions of th is paper are listed as follows. 1) For circ u lar signals, a dimension- reduced MUSIC (DR- MUSIC) algor ith m is pro posed for DOA an d p olarization estimation. Comp ared with classical lon g-vector MUSIC (L V - MUSIC) and q uaternion dimension-r e duced MUSIC (QDR- MUSIC) alg orithm [ 39 ], th e com p utational complexity of DR- MUSIC algorithm is further reduced, since the polarization estimation of DR-MUSIC ha s a closed-form expression. 2) For non- circular signals, an improved DO A estimation al- gorithm is p r oposed based o n th e URA equipp ed with EMVSs. Compared with the QDR-MUSIC algorithm, the estimation accuracy is f urther improved. Th is is bec ause a novel received data mo del is c o nstructed b a sed on quatern ion theory , and the u nconjug ated cov ar iance matrix o f n on-circu la r has been used in the pr oposed q uaternion n on-circu lar MUSIC (QNC- MUSIC) a lg orithm to im prove the DO A estimation accu racy . 3) By combin ing the QNC-MUSIC and DR-MUSIC algo - rithms, the po larization estimation can be achieved fo r non- circular sig n als. Based o n the result for the DO A estimation using QNC-MUSI C algorithm, the p olarization estimation o f non-cir cular signal can be acq uired by usin g the closed- form expre ssion of the p olarization estimation of DR-M USI C algorithm . 4) T h e com putational comp lexity of the L V -MUSIC, DR- MUSIC, QDR-MUSIC and QNC-MUSIC are analyzed. Com- pared with L V -MUSIC an d QDR-MUSIC, the com putational complexity o f DR-MUSIC and QNC-MUSIC is m uch lo wer . This advantage is particularly attractive in the massive MIMO systems, sinc e the p otentially proh ibitiv e comp utational com- plexity is o ne of the major challen ges faced by ma ssi ve MIMO systems. 5) The stochastic Cram ´ er-Rao Boun d (CRB) fo r the esti- mation of the 2- D DO A and polarization parameters of the non-cir cular signals is deriv e d, where a s the known CRB is only valid for th e estimation o f the 2-D DOA and polarizatio n parameters of the cir c ular signals. D. Organization of th e P aper This paper is organized as follows. The prob lem form ula- tions are given in Sectio n I I . T h e b asic concept and prop erty of quaternio n ar e given in Section III . The propo sed DR-MUSIC algorithm for circular sign al is p resented in Section IV . The propo sed QNC-MUSI C algorithm for non- circular sign al is presented in Section V . The comp utational com plexity analysis is given in Section VI. The stochastic CRB is der i ved in Section VI I. Th e simulation results are shown and an alyzed in Sectio n VIII. Th e co nclusions are d rawn in Section IX. E. Notatio n In th is paper, the operator ( · ) † , ( · ) ∗ , ( · ) T , ( · ) H and E {·} are co mplex matr ix pseudo - in verse, con ju gate, transp ose, con- jugate tr anspose an d expectatio n , respectively; the operator ( · ) # , ( · ) ⋄ , ( · ) ‡ and E are conjugate, transpo se, con jugate transpose a n d expectation for quaternion matrix , respectively . The boldface upperc ase letters and boldface lowercase letters denote m a trices a n d column vectors, r espectiv ely . Th e symbo l diag { z 1 , z 2 } stand s for a diagon al matrix whose diago nal entries are z 1 and z 2 , respecti vely . I M and J M stand for the M × M identity m atrix and the M × M matrix o f ones, respectively . |·| and k·k stand fo r the m o dule operator and the absolute value opera to r , respectiv ely . k·k F denotes Frobeniu s norm. ar g( · ) is the phase operato r of complex n u mbers, in radian. Symbols ⊙ and ⊗ stand for the Hadamard matrix produ ct an d the Kro n ecker product, respectively . ⊥ denotes the orth o-comp lement of a projec to r matrix. I I . P RO B L E M F O R M U L A T I O N As shown in Fig. 1 , we consider a 3-D millimeter wa ve polarized massive MI MO system with EMVSs ar ranged in a URA fo rm at the BS. Ther e are totally M = M x M y EMVSs, where M x and M y are the number s of antennas in the x - direction an d the y -directio n, respe cti vely . Obviously , the URA would b e degenerated to th e con ventional ULA when M x or M y are equal to 1 . For the m o bile terminals, each is eq uipped with one EMVS. The u plink sig n als o f the L mob ile termin a ls are non- circular signals such as BPSK mod ulated sign als going through L channels, and each has a corr esponding azimu th ang le θ l and ele vation angle ϕ l for the l th mobile terminal, which satisfy 0 ≤ θ l < π and 0 ≤ ϕ l < π/ 2 . The L channels are unco rrelated with each o ther . It shou ld be noted that the ran ges of th e DOAs are the localization ranges of the array , which means th at so u rces out of these ran g es c a n not b e localized b y the array . The steer ing vector a ( θ l , ϕ l ) ∈ C M × 1 is the response of th e a r ray co rrespond ing to th e azim uth an d elev ation DOAs of θ l and ϕ l . W ith respect to the EMVS at the origin o f the axes, the m th element of a ( θ l , ϕ l ) is deﬁned as [ a ( θ l , ϕ l )] m = exp( iu sin ϕ l [( m x − 1) cos θ l + ( m y − 1) sin θ l ]) , (1) where m = ( m y − 1) M x + m x , m x = 1 , 2 , . . . , M x , m y = 1 , 2 , . . . , M y , u = 2 π d/λ , d is the d istance between two adjacent EMVSs, λ is the wav elength. It can be seen that [ a ( θ l , ϕ l )] m correspo n ds to the r esponse of the ( m x , m y ) th EMVS in the co ordinate system shown in Fig. 1 . An EMVS equipp ed with two dipole antennas offers a good trade-o ff between perform ance and overall system de- velopment cost for the polarized massiv e MIMO system s equippe d with a large numb er o f EMVSs, thus we con sider the case that an EMVS e q uipped w ith two dipole antennas, which are arranged in the x -direction and the y -directio n , respectively , m easures the ho rizontal and vertical co mponen ts of th e electron ic ﬁeld . For the l th channe l, the comp onents of the electric ﬁeld rece ived on a n EMVS can be deﬁned as [ 26 ] ξ l ( θ l , ϕ l , γ l , η l ) =  ξ 1 l ( θ l , ϕ l , γ l , η l ) ξ 2 l ( θ l , ϕ l , γ l , η l )  =  cos θ l cos ϕ l − sin θ l sin θ l cos ϕ l cos θ l   sin γ l exp( iη l ) cos γ l  , (2) l = 1 , 2 , . . . , L , where ξ 1 l ( θ l , ϕ l , γ l , η l ) and ξ 2 l ( θ l , ϕ l , γ l , η l ) stand for the h orizontal and vertical comp o nents of the elec- tronic ﬁeld rece ived by an EMVS in the x -direction and the y -direction , respectively , while the 0 ≤ γ l < π/ 2 and 0 ≤ η l < 2 π are the rang es of the p olarization angle and p hase 4 x y z l T l M x M y M Fig. 1. Array geometry of the URA considered. The direction of the incident path is projecte d onto the array plane. T he azimuth angle , θ l , is deﬁned as the angle from the x-axis to the project ed line, and the ele vatio n DOA, ϕ l , is deﬁned as the angle from the z-axis to the incident path. The ranges of the two parameters are 0 ≤ θ l < π and 0 ≤ ϕ l < π / 2 . difference, respectively . Thu s the time-dom a in sig n als rec e i ved by th e m th E M VS eq uipped with two d ipole an tennas can be expressed as x 1 m ( t ) = L X l =1 [ a ( θ l , ϕ l )] m ξ 1 l ( θ l , ϕ l , γ l , η l ) s l ( t ) + n 1 m ( t ) , x 2 m ( t ) = L X l =1 [ a ( θ l , ϕ l )] m ξ 2 l ( θ l , ϕ l , γ l , η l ) s l ( t ) + n 2 m ( t ) , (3) m = 1 , 2 , . . . , M , where s l ( t ) is the comp lex envelope of the receiv ed signal, n 1 m ( t ) and n 2 m ( t ) are the additi ve white Gaussian no ise ( A WGN) of the m th E MVS consisting of two dipole antennas. It is to be no ted h ere that we will replace ξ 1 l ( θ, ϕ, γ , η ) , ξ 2 l ( θ, ϕ, γ , η ) and [ a ( θ l , ϕ l )] m with ξ 1 l , ξ 2 l and [ a l ] m , respectively , in th e following for notational con venienc e . Let x 1 ( t ) = [ x 11 , x 12 , . . . , x 1 M ] T , x 2 ( t ) = [ x 21 , x 22 , . . . , x 2 M ] T , then the data m odel th a t the array r eceiv ed can be expressed a s x ( t ) =  x 1 ( t ) x 2 ( t )  =  AV 1 AV 2  S ( t ) +  n 1 ( t ) n 2 ( t )  = AS ( t ) + N ( t ) , (4) where A = [ a 1 , a 2 , . . . , a L ] ∈ C M × L is the arr a y manifo ld matrix, S ( t ) = [ s 1 ( t ) , s 2 ( t ) , . . . , s L ( t )] T ∈ C L × 1 is the signal vector . V k ( t ) = diag { ξ k 1 , ξ k 2 , . . . , ξ kL } ∈ C L × L ( k = 1 , 2) is a diagona l matrix constructed by the compone nts of the electric ﬁeld. n k ( t ) = [ n k 1 ( t ) , n k 2 ( t ) , . . . , n kM ( t )] T ∈ C M × 1 ( k = 1 , 2) are the n oise vectors co m posed of temporally and spatially indepen d ent and identically d istributed (i.i.d.) circularly symmetric zero -mean Gaussian ran dom variables, who se covariance matrix is E  n k ( t ) n H k ( t )  = σ 2 n I M . It sh ould be noted that the un conjugated covariance matrix of n k ( t ) is E  n k ( t ) n T k ( t )  = 0 . Then the covariance matrix R x of x ( t ) in ( 4 ) can b e wr itten as R x = E  x ( t ) x ( t ) H  =  AV 1 AV 2  R s ( t )  AV 1 AV 2  H + σ 2 n I 2 M , (5) where R s ( t ) = E[ S ( t ) S H ( t )] is the complex covariance matrix of signal S ( t ) . Collecting N snapshots t 1 , t 2 , . . . , t N , the received d ata can be written in matrix form as X = [ x ( t 1 ) , x ( t 2 ) , . . . , x ( t N )] . (6) The DO A and p olarization estimation prob lem is stated as follows. Given the received data X , we ne ed to esti- mate the DO A an d p o larization p arameters ( θ l , ϕ l , η l , γ l ) , l = 1 , 2 , . . . , L . I I I . B A S I C C O N C E P T A N D P RO P E RT Y O F Q U A T E R N I O N Hamiltion’ s quatern io ns H is a no ntrivial generalization of comp lex numbers C . In general, qua te r nions are a four dimensiona l hyper complex num bers system, and they are an extension of com plex nu mbers to four-dimensio n al (4-D) space. Basics about quaternio n s an d their pr operties can b e found in [ 41 ] and [ 42 ]. Several basic deﬁnition s and prop erties which would be used in this pa per are introduced as follows. Deﬁnition 1 [ 34 ]: A quaternion is d escribed b y four c om- ponen ts (one real and three imaginar ies). It can be expr essed in its Cartesian f orm as α = α 0 + α 1 j 1 + α 2 j 2 + α 3 j 3 =( α 0 + α 1 j 1 ) + ( α 2 + α 3 j 1 ) j 2 = c 1 + c 2 j 2 , α 0 , α 1 , α 2 , α 3 ∈ R . (7) where j 1 j 2 = j 3 , j 2 j 1 = − j 3 , c 1 , c 2 ∈ C ar e co mplex number s. The modu lus of a qua ter nion α can be calculated as | α | = p | c 1 | 2 + | c 2 | 2 . The set of qu aternions forms a noncomm utativ e norm ed division a lg ebra, which means that gi ven two quaternions α and κ , we have ακ 6 = κα . Deﬁnition 2 [ 34 ]: I n this paper, the c onjugate of a qu ater- nion α , n oted α # , is g iv en b y α # = α 0 − α 1 j 1 − α 2 j 2 − α 3 j 3 . Deﬁnition 3 [ 41 ]: The m odulus of a qu aternion α is | α | , which can be expre ssed as | α | = √ αα # = p α 2 0 + α 2 1 + α 2 2 + α 2 3 . Property 1 [ 41 ]: For two quaternions α, κ ∈ H and a complex num ber c ∈ C with its im a ginary pa r t j 1 , noted c = a + b j 1 ∈ C , we have ( ακ ) # = κ # α # , cj 2 = j 2 c ∗ . Th us the c onjugate of α can be written as α # = ( c 1 + c 2 j 2 ) # = c ∗ 1 − j 2 c ∗ 2 = c ∗ 1 − c 2 j 2 . Remark 1 : It sho uld be n oted that th e imagin ary pa r t i used in section II is essentially identical with the imaginary part j 1 , j 2 , j 3 used in section III. In this paper, i is equ iv alent to j 1 . Deﬁnition 4 [ 34 ]: Deﬁne the ( p, q ) th entry of a P × Q matrix B as a qu aternion [ B ] p,q ∈ H , then B ∈ H P × Q is called the quatern ion ma trix. Property 2 [ 42 ]: F o r tw o qu aternion sq u are matrices B , C ∈ H P × P , α ∈ H a nd a com plex matrix D ∈ C P × Q with its imagin ary part j 1 , we ha ve ( α C ) ‡ = C ‡ α # , ( B C ) ‡ = C ‡ B ‡ , D j 2 = j 2 D ∗ . It shou ld be noted that, in general, f or the tr anspose o f two qua te r nion matrices B , C ∈ H P × P , we have ( B C ) ⋄ 6 = C ⋄ B ⋄ . Deﬁnition 5 [ 41 ]: Given a quater nion m atrix B ∈ H P × Q , The Cay ley-Dickson notatio n can be written as B = B 1 + j 2 B 2  B 1 , B 2 ∈ C P × Q  . T h en on e can d eﬁne th e complex 5 adjoint m atrix, denoted by B σ ∈ C 2 P × 2 Q , correspo nding to the q uaternion matrix , as follows B σ =  B 1 B ∗ 2 − B 2 B ∗ 1  . (8) Deﬁnition 6 [ 43 ]: Giv en a quaternion squar e ma tr ix B ∈ H P × P , there exists a qua ter nion λ ∈ H and a non-zer o vector b ∈ H P × 1 , which satisfy B b = λ b (Le f t) o r B b = b λ (Right), λ is the left or r ight eigenv alue o f the matrix B , and b is the eigenv ector o f the matrix B corresp onding to λ . Theorem 1 [ 43 ]: Gi ven a q uaternion square matrix B ∈ H P × P satisﬁes B = B ‡ , th en B is called self-con jugated matrix. Th e self-conjug ated matrix B h as equivalent rig ht and left eigenv alues. Moreover , they are both real n umbers, and they are the eigenv alues of th e co mplex representa tio n matr ix B σ as well. Pr oof: This pro of can b e found in [ 43 ]. For a quatern ion self-con jugated matrix B = B 1 + j 2 B 2 ∈ H P × P with B = B ‡ , we have B 1 = B H 1 and B 2 = − B T 2 . According to ( 8 ), it can be known that B σ is a 2 P × 2 P co m- plex Her m ite matr ix. Based on the eigenv alue de c o mposition (EVD) of the c omplex self-conjug ated matrix B σ , the EVD of the c orrespon ding q uaternion m atrix B can be acqu ired. The EVD o f the co mplex adjoint m atrix B σ can be expressed as B σ =  ˆ U 1 ˆ U ∗ 2 − ˆ U 2 ˆ U ∗ 1   ˆ Λ 0 0 ˆ Λ   ˆ U 1 ˆ U ∗ 2 − ˆ U 2 ˆ U ∗ 1  H , (9) where ˆ Λ = diag { λ 1 , λ 2 , · · · , λ P } , λ p , p = 1 , 2 , · · · , P are the eigenv alues of qu a ternion m atrix B . Th e n the EV D of a quaternio n self-con jugated matrix B = B 1 + j 2 B 2 ∈ H P × P can be expressed as B = ˆ U ˆ Λ ˆ U ‡ = ( ˆ U 1 + j 2 ˆ U 2 ) ˆ Λ ( ˆ U 1 + j 2 ˆ U 2 ) ‡ . (10) I V . D OA A N D P O L A R I Z AT I O N E S T I M A T I O N F O R C I R C U L A R S I G N A L S In this section, ﬁr st, a quatern ion dimension-re d uced MU- SIC (QDR-MUSIC) algorith m pro posed in [ 39 ] is introd uced; then a dimension -reduced MUSIC (DR-MUSIC) algo rithm is propo sed based on the partial deriv ati ve o f the spec tr um f u nc- tion [ 40 ]. Comp ared with QDR-MUSIC, th e computational complexity of our propo sed DR-MUSIC is further red uced. A. Qua ternion Dimension- Reduced MUS IC Algo rithm Based on the quater nion th e o ry , the sign al rece ived by an EMVS equipped with two dipole anten n as can b e co mbined. The received time-do main quaternion x m ( t ) ∈ H of the m th EMVS ca n be expressed as x m ( t ) = x 1 m ( t ) + j 2 x 2 m ( t ) = L X l =1 [ a l ] m ξ 1 l s l ( t ) + n 1 m ( t ) + j 2 L X l =1 [ a l ] m ξ 2 l s l ( t ) + n 2 m ( t ) ! ∆ = L X l =1 [ a l ] m ξ l s l ( t ) + n m ( t ) , (11) where ξ l = ξ 1 l + j 2 ξ 2 l ∈ H is a quaternio n , and th e q u aternion n m ( t ) = n 1 m ( t ) + j 2 n 2 m ( t ) ∈ H is the additive noise of the m th EMVS. It f ollows that x ( t ) = x 1 ( t ) + j 2 x 2 ( t ) = L X l =1 a l ξ l s l ( t ) + n ( t ) ∆ = Aξ s ( t ) + n ( t ) , (12) where ξ = diag { ξ 1 , ξ 2 , . . . , ξ L } ∈ H L × L is a quaternion matrix. When the self unc orrelated noise is indepen dent of th e signals, the cov ar iance matrix of x ( t ) that the array received can be expressed as R = E  x ( t ) x ‡ ( t )  = Aξ R s ( Aξ ) ‡ + E  n ( t ) n ‡ ( t )  = Aξ R s ( Aξ ) ‡ + 2 σ 2 n I M , (13) where R s is signal covariance matrix. E  n ( t ) n ‡ ( t )  = E [( n 1 ( t ) + j 2 n 2 ( t ))( n H 1 ( t ) − n H 2 ( t ) j 2 )] = 2 σ 2 n I M . In practical situations, the theoretical ar ray covariance ma- trices g i ven in ( 13 ) is u nav ailab le an d it can be estimated b y ˆ R = 1 N { X X ‡ } = 1 N ( X 1 + j 2 X 2 )( X 1 + j 2 X 2 ) ‡ = 1 N [ X 1 X H 1 + X ∗ 2 X T 2 + j 2 ( X 2 X H 1 − X ∗ 1 X T 2 )] ∆ = R 1 + j 2 R 2 , (14) where X ∈ H M × N , X 1 ∈ C M × N and X 2 ∈ C M × N are the snapshot data matrices con structed by x ( t ) , x 1 ( t ) and x 2 ( t ) , respectively , t = t 1 , . . . , t N . Giv en the qua ternion self-conju gated covariance matrix R ∈ H M × M , accordin g to th e quater nion for m deﬁned in ( 11 ) and the correspo nding EVD of its complex adjoin t matrix de ﬁned in ( 9 ) , th e EVD of complex Hermite adjoint matrix R σ is giv en by R σ =  R 1 R ∗ 2 − R 2 R ∗ 1  =  U 1 U ∗ 2 − U 2 U ∗ 1   Λ 0 0 Λ   U 1 U ∗ 2 − U 2 U ∗ 1  H . (15) The EVD of quatern ion cov ar iance matrix R can be ex- pressed as R = U Λ U ‡ = U S Λ S U ‡ S + 2 σ 2 n U N U ‡ N , (16) where U N = U N 1 + j 2 U N 2 ∈ H M × ( M − L ) and U N 1 , U N 2 ∈ C M × ( M − L ) . Similar to the prop erty o f comp lex matr ix MUSIC algorithm [ 44 ], the quater nion matrix U S and U N satisfy th e or th ogo- nality r elationship, i.e., the steering vector belong ing to Aξ is orthog onal to U N . The spectru m function of qu aternion dimension -reduced MUSIC (QDR-M USIC) is co nstructed as f QDR ( θ l , ϕ l ,γ l , η l ) =    ξ # l a H l U N    2 F =( ξ ∗ l 1 − ξ ∗ l 2 j 2 ) C ( θ l , ϕ l )( ξ ∗ l 1 − ξ ∗ l 2 j 2 ) ‡ (17) where C ( θ l , ϕ l ) = a H l ( U N 1 U H N 1 + U ∗ N 2 U T N 2 ) a l , l = 1 , 2 , · · · , L . In ge n eral, when 0 ≤ γ l < π / 2 , ( ξ ∗ l 1 − ξ ∗ l 2 j 2 ) 6 = 0 . 6 Thus C ( θ l , ϕ l ) = 0 means f Q ( θ l , ϕ l , γ l , η l ) = 0 . The spatial spectrum fun ction of DOA estimation is given b y f QDR ( θ l , ϕ l ) = a H l ( U N 1 U H N 1 + U ∗ N 2 U T N 2 ) a l , (18) l = 1 , 2 , · · · , L . In general, substitute the result of the DO A estimation in to ( 17 ), the polarizatio n estimation can be obtained by searching the sp e ctrum fun ction in ( 17 ). Howe ver, the pola r ization paramete r cann ot be o btained. The rea son is that after we take EVD of R , the polar ization inf o rmation is contained in the eigenv alues, not the eig en vector s, i.e. , the noise sub space U N does n ot con tain the polarization inf orma- tion. Th us the long vector (L V) metho d is r econsidered . In practical situations, the theoretical ar ray covariance ma- trices g i ven in ( 5 ) is unavailable and it can be estimated by ˆ R LV = 1 N { X X H } . (19) The EVD of covariance m atrix can be expr essed as R LV = ˆ U S ˆ Σ S ˆ U H S + ˆ U N ˆ Σ N ˆ U H N , (20) where ˆ U S and ˆ U N are the signal and noise subspace s, respectively , the diag onal matrices ˆ Σ S and ˆ Σ N respectively contain their corresp o nding eig en values. Since th e DO A param eter has already been estimated, the manifold matrix A ∈ C 2 M × L in ( 4 ) is merely the functio n of ( γ l , η l ) . T he spectr um func tio n of polarization estimation is given by f QDR ( γ l , η l ) = a H Ql ˆ U N ˆ U H N a Ql , l = 1 , 2 , · · · , L, (21) where a Ql is the steering vector belong ing to A . Thus the pseudo code of the QDR-MUSIC can be summarize d as Algorithm 1 . Algorithm 1 QDR-MUSIC fo r Circu lar Sig n al 1: Estimate ˆ R according to ( 14 ); 2: Constru ct th e complex a djoint matrix R σ of ˆ R acco rding to ( 9 ); 3: T ake the EVD of R σ accordin g to ( 15 ); 4: Constru ct the spectrum fun ction f QDR ( θ l , ϕ l ) accor d ing to ( 17 ); 5: Search ( 18 ) to obtain the spectrum extremum ( ˆ θ l , ˆ ϕ l ) , l = 1 , 2 , · · · , L ; 6: Estimate ˆ R LV accordin g to ( 19 ); 7: T ake the EVD of R LV accordin g to ( 20 ) to acq uire ˆ U N ; 8: Estimate ( ˆ γ l , ˆ η l ) , l = 1 , 2 , · · · , L by search ing ( 2 1 ). B. Dimen sion-Reduce d MUSIC A lgorithm 1) DO A Estimation: Based on th e princ iple of MUSIC algorithm , the subspace ˆ U N is ortho gonal to span  AV 1 AV 2  , which is th e space span ned by the array manifold matr ix, the n it ho lds that  AV 1 AV 2  H ˆ U N = V H 1  A AV 2 V − 1 1  H ˆ U N = 0 . (22) Partitioning ˆ U N into two block matr ices ˆ U N =  ˆ U N 1 ˆ U N 2  , wh ere ˆ U N 1 ∈ C M × (2 M − L ) and ˆ U N 2 ∈ C M × (2 M − L ) are noise subspaces, th en we hav e  a 1 a 2 · · · a L ξ 21 ξ 11 a 1 ξ 22 ξ 12 a 2 · · · ξ 2 L ξ 1 L a L  H  ˆ U N 1 ˆ U N 2  = 0 . (23) Deﬁne a complex numb er as follows ρ l exp ( iδ l ) = ξ 2 l ξ 1 l . (24) Based on ( 22 ), we c an construc t a spectru m function f DL ( θ, ϕ, γ , η ) as follows f DR ( θ l , ϕ l , γ l , η l ) =      a H l ρ l exp ( − iδ ) a H l   ˆ U N 1 ˆ U N 2      2 F = a l H ˆ U N 1 ˆ U H N 1 a l + ρ l exp ( − iδ ) a H l ˆ U N 2 ˆ U H N 1 a l + ρ l exp ( iδ ) a H l ˆ U N 1 ˆ U H N 2 a l + ρ 2 l a H l ˆ U N 2 ˆ U H N 2 a l . (25) By setting that the p artial derivati ve of f DR ( θ l , ϕ l , γ l , η l ) with respect to ( γ l , η l ) equal to zero, we can obtain from ( 25 ) that exp ( − iδ l ) a H l ˆ U N 2 ˆ U H N 1 a l = ex p ( i δ l ) a H l ˆ U N 1 ˆ U H N 2 a l = − ρ l a H l ˆ U N 2 ˆ U H N 2 a l , (26) the p roof is g i ven in the supplemental materials. Th e spatial spectrum f DR ( θ l , ϕ l , γ l , η l ) in ( 25 ) minimizes accordin g to ( γ l , η l ) when ( 2 6 ) takes the minus value. Substituting ( 26 ) into ( 2 5 ), we c a n simplify ( 25 ) as f DR ( θ l , ϕ l , ρ l ) = a H l ˆ U N 1 ˆ U H N 1 a l − ρ 2 l a H l ˆ U N 2 ˆ U H N 2 a l . (27) According to ( 26 ), exp ( iδ l ) a H ˆ U N 1 ˆ U H N 2 a l is a real num ber , we h a ve exp ( iδ ) a H ˆ U N 1 ˆ U H N 2 a l = r    a H ˆ U N 1 ˆ U H N 2 a l    2 F . (28) Then ρ l is given by ρ l = − r    a H l ˆ U N 1 ˆ U H N 2 a l    2 F a H l ˆ U N 2 ˆ U H N 2 a l . (29) Substituting ( 29 ) into ( 2 7 ), we can simplify the spectrum function ( 2 7 ) for the DO A estimation a s f DR ( θ l , ϕ l ) = a H l ˆ U N 1 ˆ U H N 1 a l −    a H l ˆ U N 1 ˆ U H N 2 a l    2 F a l H ˆ U N 2 ˆ U H N 2 a l . (30) Now the p a r ameter search dimension of ( 30 ) has been reduced to two d imensions. 2) P o larization Estimation: The polariza tio n estimation ( ˆ γ l , ˆ η l ) can be obtained when L DO A estimation s h av e be e n acquired . Based on ( 24 ), ( 28 ) and ( 29 ), we have c l = −    ˆ a H l ˆ U N 1 ˆ U H N 2 ˆ a l    2 F ˆ a H l ˆ U N 2 ˆ U H N 2 ˆ a l 1 ˆ a H l ˆ U N 1 ˆ U H N 2 ˆ a l = ξ 2 l ξ 1 l , (31) where the steerin g vecto r ˆ a l is constructed b y ( ˆ θ l , ˆ ϕ l ) which has alr eady been estimated from ( 30 ) . 7 Substituting ( 2 ) into ( 31 ), we have sin ˆ θ l cos ˆ ϕ l sin γ l exp ( iη l ) + co s ˆ θ l cos γ l cos ˆ θ l cos ˆ ϕ l sin γ l exp ( iη l ) − sin ˆ θ l cos γ l = c l , (32) which can be sim p liﬁed as tan( γ l )exp ( iη l ) = cos ˆ θ l + c l sin ˆ θ l − sin ˆ θ l cos ˆ ϕ l + c l cos ˆ θ l cos ˆ ϕ l . (33) Thus th e polar ization estimatio n of the l th signal is expr e ssed as ˆ γ l = arctan      cos ˆ θ l + c l sin ˆ θ l − sin ˆ θ l cos ˆ ϕ l + c l cos ˆ θ l cos ˆ ϕ l      ! , ˆ η l = arg cos ˆ θ l + c l sin ˆ θ l − sin ˆ θ l cos ˆ ϕ l + c l cos ˆ θ l cos ˆ ϕ l ! . (34) Based o n the m ethod men tio ned above, th e DO A an d polarization estimation are ach ie ved. Thu s the pseudo code of the DR-MUSIC c an be su mmarized as Algorithm 2 . Algorithm 2 DR-MUSIC for Circu lar Sig n al 1: Estimate ˆ R LV accordin g to ( 19 ); 2: T ake the EVD of R LV accordin g to ( 20 ); 3: Constru ct the spectrum fu nction f Q ( θ l , ϕ l ) ac c o rding to ( 30 ); 4: Search ( 3 0 ) to obtain the spec trum extrem um ( ˆ θ l , ˆ ϕ l ) , l = 1 , 2 , · · · , L . 5: Calculate c l accordin g to ( 31 ); 6: Estimate ( ˆ γ l , ˆ η l ) , l = 1 , 2 , · · · , L accordin g to ( 34 ). The QDR-MUSIC and the DR-MUSIC can b oth be used for cir cular and n on-circu lar signals. Howe ver, the property of n on-circu lar signal, which can be utilized to imp r ove the parameter estimation perfo rmance, has not b een u sed in the QDR-MUSIC an d the DR-M U SI C. V . D OA A N D P O L A R I Z A T I O N E S T I M A T I O N F O R N O N - C I R C U L A R S I G N A L S For a n on-circu lar signal s , it holds that [ 13 ] E [ s ( t ) s ( t )] = µ exp ( iβ )E [ s ( t ) s ∗ ( t )] , (35) in which β is th e non-cir cularity phase, µ is the non -circularity rate w ith µ = 1 fo r the maximal non-circ ularity rated signal and 0 < µ < 1 fo r the c o mmon non- circularity rated sign a l. For signal vector S ∈ C L × 1 consisting of L ind ependent compon ents, its unco n jugated covariance m atrix is gi ven b y R ′ s =E  S ( t ) S T ( t )  =diag { µ 1 exp ( iβ 1 )E [ s 1 ( t ) s 1 ∗ ( t )] , µ 2 exp ( iβ 2 )E [ s 2 ( t ) s 2 ∗ ( t )] , . . . , µ L exp ( iβ L )E [ s L ( t ) s L ∗ ( t )] } ∆ = P B R s , (36) where P = diag { µ 1 , µ 2 , . . . , µ L } is a r eal valued diag onal matrix consisting of the non-circularity rates of L signals. B = diag { exp ( i β 1 ) , e x p ( i β 2 ) , . . . , e x p ( iβ L ) } is a diagon a l matrix consisting of their no n-circular ity phases. F or the maximal n on-circu larity rated sig n als, we ha ve P = I L . A. DO A Estimation In order to utilize the inf o rmation contained in the un- conjuga te d covariance matr ix of th e no n-circular signals, two quaternio n vectors a re constru cted as y ( t ) = x 1 ( t ) + j 2 x 2 ( t ) , z ( t ) = x ∗ 1 ( t ) + j 2 x ∗ 2 ( t ) . (37) Substituting ( 4 ) into ( 37 ), we have y ( t ) = A ∗ ( V 1 + j 2 V 2 ) s ( t ) + ( n 1 ( t ) + j 2 n 2 ( t )) ∆ = A ∗ V y s ( t ) + n y ( t ) , (38) z ( t ) = A ( V ∗ 1 + j 2 V ∗ 2 ) s ∗ ( t ) + ( n ∗ 1 ( t ) + j 2 n ∗ 2 ( t )) ∆ = AV z s ∗ ( t ) + n z ( t ) . (39) where the d iagonal entries o f the q uaternion matrices V y and V z are v y l and v z l , l = 1 , 2 , · · · , L , respectively , w h ich can be expressed as v y l = ξ 1 l + j 2 ξ 2 l , v z l = ξ ∗ 1 l + j 2 ξ ∗ 2 l . (40) The qu aternion vector w ( t ) ∈ H 2 M × 1 is constructed as w ( t ) =  y ( t ) z ( t )  =  x 1 ( t ) x ∗ 1 ( t )  + j 2  x 2 ( t ) x ∗ 2 ( t )  = w 1 ( t ) + j 2 w 2 ( t ) . (41) The quaternion covariance m atrix R w ∈ H 2 M × 2 M of the extended vector w ( t ) is expressed as R w = E  w ( t ) w ‡ ( t )  =  E  y ( t ) y ‡ ( t )  E  y ( t ) z ‡ ( t )  E  z ( t ) y ‡ ( t )  E  z ( t ) z ‡ ( t )   . (42) According to ( 38 ) and ( 39 ), we hav e E  y ( t ) y ‡ ( t )  = A ∗ V y R S V ‡ y A T + E  n y ( t ) n y ‡ ( t )  = A ∗ V y R s V ‡ y A T + 2 σ 2 n I M , (43) E  y ( t ) z ‡ ( t )  = A ∗ V y R ′ S V ‡ z A H + E  n y ( t ) n z ‡ ( t )  = A ∗ V y R ′ s V ‡ z A H , (44) E  z ( t ) z ‡ ( t )  = AV z R S V ‡ z A H + E  n z ( t ) n z ‡ ( t )  = AV z R s V ‡ z A H + 2 σ 2 n I M . (45) Substituting ( 43 ), ( 44 ) and ( 45 ) into ( 42 ) an d u sing th e relationship R ′ s = P B R s , we can r ewrite R w as R w =  A ∗ V y R s V ‡ y A T A ∗ V y B R s V ‡ z A H AV z B ∗ R s V ‡ y A T AV z R s V ‡ z A H  + 2 σ 2 n I 2 M =  A ∗ V y AV z B ∗  R s  A ∗ V y AV z B ∗  ‡ + 2 σ 2 n I 2 M ∆ = ˆ AR s ˆ A ‡ + 2 σ 2 n I 2 M . (46) In practical situations, the theoretical ar ray covariance ma- trices given in ( 42 ) is unavailable an d it is usu a lly estimated by ˆ R w = 1 N { W W ‡ } = 1 N ( W 1 + j 2 W 2 )( W 1 + j 2 W 2 ) ‡ = 1 N [ W 1 W H 1 + W ∗ 2 W T 2 + j 2 ( W 2 W H 1 − W ∗ 1 W T 2 )] ∆ = R w 1 + j 2 R w 2 , (47) 8 where W ∈ H 2 M × N , W 1 ∈ C 2 M × N and W 2 ∈ C 2 M × N are the snapshot data matrices construc ted by w ( t ) , w 1 ( t ) and w 2 ( t ) , respectiv ely , t = t 1 , . . . , t N . Similar to ( 15 ), the EVD of co mplex Hermite adjoint matrix R σ w can be obtained as we ll, thus the EVD of R w can be expressed as R w = U S Σ S U ‡ S + U N Σ N U ‡ N . (48) Based on th e ortho gonality o f subsp a ce U S ∈ H 2 M × L and U N ∈ H 2 M × (2 M − L ) , it holds th a t ˆ A ‡ U N = 0 . (49) Let us partition the noise subspa c e into two block matrices U N =  U N 1 U N 2  . Then the sp e ctrum function of quaternio n non-cir cular MUSI C (QNC-MUSIC) is constructed as f QN C ( θ l ,ϕ l , γ l , η l ) =      a ∗ l v y l a l v z l exp ( − i β l )  ‡  U N 1 U N 2      2 F . (50) W e d e ﬁne three qu aternion n umbers Π 1 = a T U N 1 U ‡ N 1 a ∗ , Π 2 = a H U N 2 U ‡ N 2 a a nd Π 3 = a H U N 2 U ‡ N 1 a ∗ . Th en ( 50 ) can be written as f QN C ( θ l ,ϕ l , γ l , η l ) =  | ξ l 1 | 2 + | ξ l 2 | 2 + 2real { iξ l 1 ξ ∗ l 2 }  Π 1 +  | ξ l 1 | 2 + | ξ l 2 | 2 − 2real { iξ l 1 ξ ∗ l 2 }  Π 2 + exp ( iβ l )  ξ l 1 2 + ξ l 2 2  Π ∗ 3 + exp ( − iβ l )  | ξ ∗ l 1 | 2 + | ξ ∗ l 2 | 2  Π 3 . (51) Let th e partial deriv ative o f f QN C ( θ l , ϕ l , γ l , η l ) with r e sp ect to β l equal to zero, we can get ∂ f QN C ( θ l , ϕ l , γ l , η l ) ∂ β l = i exp ( iβ l )  ξ l 1 2 + ξ l 2 2  Π ∗ 3 − i exp ( − iβ l )  | ξ ∗ l 1 | 2 + | ξ ∗ l 2 | 2  Π 3 =0 . (52) The no n-circular ity phase after calculating ( 52 ) is giv en by exp ( iβ l ) = −  | ξ ∗ l 1 | 2 + | ξ ∗ l 2 | 2  Π 3 | ( | ξ l 1 | 2 + | ξ l 2 | 2 ) Π ∗ 3 | . (53) Substituting ( 53 ) into ( 51 ), we can obtain the simpliﬁed spectrum fu nction f or DO A and po larization estimation of non - circular sig n als as f QN C ( θ l , ϕ l , γ l , η l ) =  | ξ l 1 | 2 + | ξ l 2 | 2 + 2real { iξ l 1 ξ ∗ l 2 }  Π 1 +  | ξ l 1 | 2 + | ξ l 2 | 2 − 2real { iξ l 1 ξ ∗ l 2 }  Π 2 − 2    | ξ l 1 | 2 + | ξ l 2 | 2  Π ∗ 3   . (54) Howe ver, the DOAs of the in cident signals can be estimated by the propo sed a lg orithm, but the polarizatio ns of the inciden t signals ca n not be estimated by the prop osed algorithm . One p ossible re a son is that the co nstraint condition of the polarization p arameters is relatively weak. The DO A parameters ar e con tain ed not on ly in steering vector a l , b u t also in ξ 1 l and ξ 2 l , The polarization param e ters are merely contained in ξ 1 l and ξ 2 l . Thus the conditio n that the po - larization estimation with mu ltiple solutions may hap pen by searching ( 54 ). Anoth er p ossible reason is that the po la r ization parameters ar e presen ted as a quaternio n ( this can be fou nd in ( 40 )). After we take the EVD of R w , a majo r part of the polarization inf ormation is con tained in the eigenv alues, not the eigen vector s, i.e., th e noise sub space U N sparingly contains the polarization informa tio n. Thu s the p olarization parameter can not be estimated using ( 49 ). Algorithm 3 QNC-MUSI C for No n-Circular Signa ls 1: Estimate ˆ R w accordin g to ( 47 ); 2: Constru ct th e co mplex adjo int m atrix R σ w of ˆ R w accord- ing to ( 9 ); 3: T ake the EVD of R σ w accordin g to ( 48 ); 4: Constru ct th e spectrum function f QN C ( θ l , ϕ l , γ l , η l ) ac- cording to ( 54 ); 5: Search ( 5 4 ) to ob tain the sp ectrum extremum ( θ l , ϕ l ) , l = 1 , 2 , · · · , L . 6: Estimate ˆ R LV accordin g to ( 19 ); 7: T ake the EVD of R LV accordin g to ( 20 ); 8: Partition ˆ U into ˆ U N 1 and ˆ U N 2 ; 9: Calculate c l accordin g to ( 31 ); 10: Estimate ( γ l , η l ) , l = 1 , 2 , · · · , L accordin g to ( 34 ). B. P ola rization Estimation The n oise su b space co ntaining the po larization inf o rmation has been estimated fro m ( 2 0 ), thus the p olarization e stimation ( ˆ γ l , ˆ η l ) ca n be obtained as well w h en L DO As of non-c ircular signals h av e been acquired from ( 54 ). The steerin g vector a l is con structed by ( θ l , ϕ l ) . Based on ( 3 1 ), we hav e d l = −    a l H ˆ U N 1 ˆ U H N 2 a l    2 F a l H ˆ U N 2 ˆ U H N 2 a l a l H ˆ U N 1 ˆ U H N 2 a l . (55) Then the polarization estimation γ l , η l of the l th non-circu lar signal can be achieved as the step 5 and 6 in Alg o rithm 2 . Thus the DO A an d polarizatio n estima tio n for non-circu lar signals are achiev ed . Th e pseudo cod e o f the QNC-MUSIC can be summarize d as Algorithm 3 . V I . C O M P U TA T I O N A L C O M P L E X I T Y A N A LY S I S In this section, we an alyze the computational complexity of the QNC-MUSIC and th e DR-MUSIC and compare their computatio nal com plexity of th em to those of the classical L V - MUSIC and the QDR-MUSIC. The com putational comp lexity of th e MUSIC-ba sed algorithm concen trates on th e estimation of the covariance matrix, th e EVD of the covariance matrix, the spectru m search of th e DOA param e ter ( θ , ϕ ) and the spec- trum search of the polarization parame ter ( γ , η ) , r espectiv ely . W e comp are the fou r algorithms in the fo ur parts men tioned above. (1) The co mputation al complexity of the e stimation of th e covariance matrix For th e classical L V -MUSIC, the dimen sion of the received data m atrix is 2 M × N , thus 4 M 2 N ﬂop s are required for the 9 calculation of the covariance ma trix. For th e QDR-MUSIC and the QNC-MUSIC, th e comp lex adjoint m a tr ices ( 14 ) and ( 47 ) of the q uaternion covariance matrix R and R w require 4 M 2 N ﬂops to b e calculate d , respectively . For the DR-MUSIC, the dimension of the received data matr ix is iden tical with that of the L V -MUSIC, thus 4 M 2 N ﬂops are required for the calculation of the covariance m atrix as well. (2) The compu tational complexity of the EVD on the co- variance m atrix For the L V -MUSIC and the DR-MUSIC, the EVD of the covariance matrix are using the fast sub space decompo si- tion ( FSD) techniq ue [ 45 ], and their complexity is gi ven by 4 M 2 ( L + 2) ﬂo p s. For the QDR-MUSIC and the QNC-M USIC using the FSD techniqu e, the comp utational co mplexity o f the EVD of the complex adjoint matrices R σ ∈ C 2 M × 2 M and R σ w ∈ C 4 M × 2 M are 4 M 2 ( L + 2 ) and 16 M 2 ( L + 2) ﬂops, respectively . (3) The co mputationa l complexity of the spectrum sear ch of the D O A parameter ( θ, ϕ ) Let J 1 , J 2 , J 3 and J 4 stand fo r the number of spe c tral points for par ameter θ , ϕ , γ and η , respectively . For the L V - MUSIC, it has to compute k a H LV ( θ, ϕ, γ , η ) ˆ U N k 2 for each spectral point, the spectral search step costs J 1 J 2 J 3 J 4 (2 M + 1)(2 M − L ) ﬂops, where a LV ( θ, ϕ, γ , η ) ∈ C 2 M × 1 is the steering vector belon ging to the L V -MUSIC, and ˆ U N ∈ C 2 M × (2 M − L ) is the noise subspac e . For the QDR-MUSIC, it h as to co mpute a H ( θ, ϕ ) ( U N 1 U H N 1 + U ∗ N 2 U T N 2 ) a ( θ , ϕ ) for each spectr al point, where a ( θ , ϕ ) ∈ C M × 1 is the steering vector , a n d U N 1 , U N 2 ∈ C M × ( M − L ) are the noise subspaces, then th e sp ectral search step costs 2 J 1 J 2 ( M + 1)( M − L ) . For the QNC-MUSIC, it has to compute Π 1 = a T ( θ, ϕ ) U N 1 U ‡ N 1 a ∗ ( θ, ϕ ) , Π 2 = a H ( θ, ϕ ) U N 2 U ‡ N 2 a ( θ , ϕ ) and Π 3 = a H ( θ, ϕ ) U N 2 U ‡ N 1 a ∗ ( θ, ϕ ) for each spec- tral point, where a ( θ , ϕ ) ∈ C M × 1 is the steering vec- tor , U N 1 ∈ H M × (2 M − L ) and U N 2 ∈ H M × (2 M − L ) are the quaternion matrices. Since the intermediate re- sults o f the calcu lation o f Π 1 and Π 2 can be used for the calcu lation of Π 3 , th e spectral search step costs 4 J 1 J 2 ( M + 1)(2 M − L ) + 2 J 1 J 2 (2 M − L ) ﬂop s. For the DR-MUSIC, it has to co mpute a H ( θ, ϕ ) ˆ U N 1 ˆ U H N 1 a ( θ , ϕ ) , a H ( θ, ϕ ) ˆ U N 2 ˆ U H N 2 a ( θ , ϕ ) and a H ( θ, ϕ ) ˆ U N 1 ˆ U H N 2 a ( θ , ϕ ) for each spectr al point, where a ( θ, ϕ ) ∈ C M × 1 is the steering vector , ˆ U N 1 ∈ C M × (2 M − L ) and ˆ U N 2 ∈ C M × (2 M − L ) are complex m atrices. The spectral search step costs 2 J 1 J 2 ( M + 1)(2 M − L ) + J 1 J 2 (2 M − L ) ﬂops, which is similar as the calculation pro cess of the QNC-MUSIC. (4) The co mputationa l complexity of the spectrum sear ch of the p olarization pa rameter ( γ , η ) For the L V - MUSIC, the four dim ensional (4 -D) parameter search is u sed for jo in tly estimatin g θ , ϕ , γ and η , th us the total complexity f or the DO A and po la r ization estimation is J 1 J 2 J 3 J 4 (2 M + 1)(2 M − L ) ﬂo p s. For the QDR-MUSIC, ˆ U N ∈ C 2 M × (2 M − L ) in ( 20 ) is need ed to be calculated , thus the EVD of ( 20 ) costs 4 M 2 N + 4 M 2 ( L + 2) ﬂops. For each in cident signal, it has to compute k a H Ql ( γ , η ) ˆ U N k 2 for each spectral point, the spectral search step c osts J 3 J 4 (2 M + 1)(2 M − L ) ﬂops, wher e a Ql ∈ C 2 M × 1 is the steering vector of the l th signal. The total computational co mplexity of th e Q D R-M USIC is J 3 J 4 L (2 M + 1)(2 M − L ) + 4 M 2 N + 4 M 2 ( L + 2) . For the QNC-MU SI C, ˆ U N 1 ∈ C M × (2 M − L ) and ˆ U N 2 ∈ C M × (2 M − L ) in ( 5 5 ) are n eeded to be calcu lated, thus the EVD of ( 20 ) has to be taken . According to ( 3 4 ), it can be seen that an analytical solu tion exists fo r the p olarization parameter ( γ , η ) , the spectrum search is not needed. Thus the computationa l c omplexity f or the p olarization parameter ( γ , η ) estimation is 4 M 2 N + 4 M 2 ( L + 2) ﬂops. For the DR- MUSIC, the spectrum sear ch is no t need ed for the p olarization parameter ( γ , η ) estimation, which is similar as that of th e QNC-MUSIC. W e c a n summar ize the compu tational co mplexity of the four algorithm s as follows . The computationa l c omplexity of the L V -MUSIC, the QDR-MUSIC, QNC-MUSIC an d D R-M USIC are given by C LV = J 1 J 2 J 3 J 4 (2 M + 1)(2 M − L ) + 4 M 2 ( N + L + 2) ﬂop s, C QDR = 2 J 1 J 2 ( M + 1)( M − L ) + J 3 J 4 L (2 M + 1)(2 M − L ) + 8 M 2 ( N + L + 2) ﬂops, C QN C = 4 J 1 J 2 ( M + 1)(2 M − L ) + 2 J 1 J 2 (2 M − L ) + 2 0 M 2 ( L + 2) + 8 M 2 N ﬂo ps. and C DR = 2 J 1 J 2 (2 M + 1)(2 M − L ) + J 1 J 2 (2 M − L ) + 4 M 2 ( N + L + 2) ﬂops, resp ecti vely . The com putational com plexity of the L V - MUSIC ( C LV → 4 M 2 J 1 J 2 J 3 J 4 ) , QDR-MUSIC ( C QDR → 2 M 2 J 1 J 2 + 4 LM 2 J 3 J 4 ) , QNC-MUSIC ( C QN C → 8 M 2 J 1 J 2 ) and th e DR-MUSIC ( C DR → 8 M 2 J 1 J 2 ) , respectively , as ( M → ∞ ) . When the so u rce n u mber increases, the com putational com- plexity of the QDR-MUSIC would increase greatly . Obviously , in terms of implemen tation, the L V -MUSIC and the QDR- MUSIC algorithm s are signiﬁcan tly more complicated than that of th e QN C-M USIC an d the DR-MUSIC algorithm s in the con text o f the m assi ve MIMO systems. V I I . S T O C H A S T I C C R A M ´ E R - R AO B O U N D F O R N O N - C I R C L A R S O U R C E S The stochastic CR B [ 46 ]–[ 49 ] can b e a chiev ed asym ptot- ically (in the num ber of m e a surements) by the stochastic maximum likelihoo d (M L) method as the (asymp totic) covari- ance matrix of th e ML estimator . Th e parameter estimation can achieve a h igher estimation accu racy and can resolve up to twice as many sign al source s compa r ed to the traditional methods for a r bitrary signals via e xploiting the stru c ture of second-o rder non -circular (NC) signals [ 48 ], [ 49 ]. In the following, we will b e concer n ed with the sign a l model x ( t ) = A ( ϑ ) S ( t ) + N ( t ) , t = 1 , ..., N , (56) where A ( ϑ ) ∈ C 2 M × L is the steering matrix . The parameter ϑ T = [ θ T , ψ T , γ T , η T ] T ∈ R 4 L × 1 and θ , ψ , γ , η a r e all in R L . Non-Circlar Sources S ( t ) with R s = E { S ( t ) S ( t ) H } and R ′ s = E { S ( t ) S ( t ) T } , circular com plex noise N ( t ) with uncon ju gated covariance matrix E { N ( t ) N ( t ) T } = 0 . Th us we h a ve R x = AR s A H + σ 2 n I 2 M and R ′ x = AR ′ s A T . (57 ) If no a pr iori informa tio n is av ailable, ( R x , R ′ x ) is ge n er- ically parame te r ized by th e real unkn own parameter vector 10 Θ := [ ϑ T , ¯ α T , σ n ] T and ¯ α = h Re  [ R s ] ¯ i, ¯ j  , Im  [ R s ] ¯ i, ¯ j  , Re  [ R ′ s ] ¯ i, ¯ j  , Im  [ R ′ s ] ¯ i , ¯ j   1 ≤ ¯ j < ¯ i ≤ L ,   [ R s ] ¯ i , ¯ i  , Re  [ R ′ s ] ¯ i , ¯ i  , Im  [ R ′ s ] ¯ i , ¯ i   ¯ i =1 ,...,L i T ∈ R  L 2 + L ( L +1)  × 1 . (58) The prob ability d ensity function of x ( t ) is expre ssed as a function ˜ x ( t ) := h x ( t ) x ∗ ( t ) i ∈ C 4 M in the case of uniform white n oise. p ( ˜ x ( t )) = ( π ) − 2 M   R ˜ x   − 1 2 exp h − 1 2 ˜ x H R − 1 ˜ x ˜ x i , (59 ) where R ˜ x = E { ˜ x ˜ x H } = ˜ AR ˜ s ˜ A H + R ˜ n (60) with R ˜ s =  R s R ′ s R ′∗ s R ∗ s  ˜ A =  A 0 2 M × L 0 2 M × L A ∗  (61) and R ˜ n =  σ 2 n I 2 M 0 2 M × 2 M 0 2 M × 2 M σ 2 n I ∗ 2 M  = σ 2 n I 4 M . (62) W e note that the log-likelihood function associated with the PDF ( 59 ) can be classically written as L ( Θ ) = − N 2  ln  | R ˜ x |  + T r  R − 1 ˜ x R ˜ x ,N   (63) with R ˜ x ,N := 1 N · P N t =1 ˜ x ( t ) ˜ x H ( t ) , wh ere N is the snapsho t number . Due to the structures of R ˜ s and R ˜ n in R ˜ x , the ML estimation of Θ can be obtain ed in a sep arable for m. The ML estimation of ϑ is g ot by minimizin g with r espect to parame ter ϑ . F N ( ϑ ) = ln    ˜ A ˆ R ˜ s ,M L ˜ A H + ˆ σ 2 n,M L I 4 M    (64) where ˆ R ˜ s ,M L , an d ˆ σ 2 n,M L are gi ven b y ˆ R ˜ s ,M L =  ˜ A H ( ϑ ) ˜ A ( ϑ )  − 1 ˜ A H ( ϑ ) ·  R ˜ x ,N − ˆ σ 2 n,M L I 4 M  ˜ A ( ϑ )  ˜ A H ( ϑ ) ˜ A ( ϑ )  − 1 (65) and ˆ σ 2 n,M L = 1 2 M − L T r  Π ⊥ A ( ϑ ) R x ,N  (66) where Π A ( ϑ ) is the projection matrix A ( ϑ )  A H ( ϑ ) A ( ϑ )  − 1 A H ( ϑ ) . Follo wing the line of derivation giv en in [ 46 ], [ 48 ], the CRB of par ameter ϑ is given by C ϑ =  F ′′ ( ϑ )  − 1 ·  lim N →∞ E  F ′ N ( ϑ )  F ′ N ( ϑ )  T    F ′′ ( ϑ )  − 1 (67) where F ′ ( ϑ ) is the g r adient of F N ( ϑ ) , and F ′′ ( ϑ ) is the limit of the Hessian o f F N ( ϑ ) wh en N → ∞ . The deriv ation details are g i ven in the su pplemental materials, then th e CRB o f no n- circular sig n als C ( N C ) ϑ is given by C ( N C ) ϑ = σ 2 n 2 n Re h D H Π ⊥ A D ⊙  J 4 ⊗   R s A H , R ′ s A T  R − 1 ˜ x h AR s A ∗ R ′∗ s i T io − 1 ∈ R 4 L × 4 L , (68) where D := h ∂ A ∂ ϑ i ∈ C 2 M × 4 L (69) V I I I . S I M U L A T I O N R E S U LT S In this section, we show n u merical results to demonstrate the perform ance of the proposed algorithms. W e com pare the pro posed QNC-MUSIC, DR-MUSIC algorithms with the ESPRIT (On ly θ and ϕ is estimated by using the 2-D ESPRIT , the p olarization p arameter is based on spectr al search which is similar as the QNC alg orithm.) and Q-MUSIC [ 13 ], th e Q D R- MUSIC algorithm s introdu c ed in section IV . A [ 39 ], r e spec- ti vely . The L V - MUSIC is not considered in o ur simulations because of its prohibitive co mputation al co mplexity . For all simulations, we consider a URA equ ipped with EM VSs sh own in Fig. 1 . The horizonta l inter-EMVS sp a cing d 1 and vertical inter-EMVS sp acing d 2 are set to be λ/ 2 . The data symbols are BPSK signal (non-circular signal) with un it power . The signal-to-n oise r a tio (SNR) is deﬁned as 10lo g 10 (1 /σ 2 ) wh ere σ 2 is the no ise p ower . T he search step size 0 . 1 ◦ has been u sed for θ and ϕ ; 0 . 03 ◦ for γ and η . The nu mber of independ ent trials is 20 0. Th e metric o f root-m ean-square error ( RMSE) is ev aluated f or the estimations of various sour ce param e ters. The simulation parame te r s in the ﬁrst two simu lations as shown in Fig. 2 and Fig. 3 are given as follows. T h e number of the mo b ile termin als is L = 3 . The azimuth DOAs of thr ee mobile terminals are θ 1 = 20 ◦ , θ 2 = 75 ◦ , θ 3 = 115 ◦ , and the corr espondin g elev ation DO As are ϕ 1 = 1 0 ◦ , ϕ 2 = 15 ◦ , ϕ 3 = 20 ◦ ; the cor respondin g po larization angles are γ 1 = 40 ◦ , γ 2 = 7 0 ◦ , γ 3 = 2 5 ◦ , and th e co rrespond ing p hase differences are η 1 = 2 0 ◦ , η 2 = 40 ◦ , η 3 = 30 ◦ . In the ﬁrst test as shown in Fig. 2 , the average received SNR from each mob ile term inal is 0 dB, and the sn apshot numb e r is N = 20 0 . T h e numb e r of EMVSs of the BS in the x -direction and the y -direction satisfy M x = M y = √ M . The RMSEs of the estimated DO A a n d p olarization parameter versus the number of the BS an tennas M are d epicted in Fig. 2 . It c an be ob served that the RMSEs of these estimated para m eters of the three algorithms decrease as M increases. For the DO A estimation, the RMSEs of the DR-MUSIC are smaller than those of the Q-MUSIC and QDR-MUSIC. The reason is that the dimension of the covariance matrix ˆ R LV is larger than that of the covariance matrix ˆ R , the superio rity of the L V -MUSIC has been reserved in th e DR-MUSIC. Th e RMSEs of the QNC- MUSIC are smaller than those of the o ther algo rithms and approa c h the CRB closely . This is because that the prop erty of the no n-circular signal has b een used in the QNC-MUSIC to improve the accuracy of the DOA estimation. Howe ver, ESPRIT p erforms the w orst o f all, since the a r ray ap e r ture is not utilized completely . Th e R MSEs of the polarization parameter o f the ESPRIT , the Q-MUSI C and th e QDR-MUSIC are larger than those of th e DR-MUSIC and the QNC-M USIC 11         1XPEHURIWKH%6DQWHQQDV0         (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%         1XPEHURIWKH%6DQWHQQDV0       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%         1XPEHURIWKH%6DQWHQQDV0       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%         1XPEHURIWKH%6DQWHQQDV0       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5% 506(RI ș HVWLPDWHVGHJUHHV 506(RIĳHVWLPDWHVGHJUHHV 506(RIȖHVWLPDWHVGHJUHHV 506(RIȘHVWLPDWHVGHJUHHV D E F G Fig. 2. RMSEs versus the number of the BS antenna s M for the estimates of DOA and polarizatio n parameters when using dif ferent estimation algorithms, and the av erage recei ved SNR from each mobile terminal is 0 dB. (a), (b), (c) and (d) correspond to the estimation of azimuth, ele vat ion, polarizati on angle and phase diffe rence, respecti vely . with M = 2 5 , but the RMSEs of the polar ization para meter of the ESPRIT , the Q-MUSI C and the QDR-MUSIC are smaller than those of the two alg orithms a s M further increases. This is because that the po larization parameter e stima tion of the DR- MUSIC and the QNC-MUSIC has a c losed-form expression, and the polar iza tion par ameter estimation of the ESPRIT , the Q-MUSIC a nd the QDR-MUSIC based o n spectru m search exploiting the or thogon a lity between the signal an d noise subspaces, i.e., k a H Ql ( γ , η ) ˆ U N k 2 . For the three algorithms, the estimation accu racy of the po larization parameter depend s on n ot only the estimation accuracy of the DOA parame te r, but also the d imension o f the estimated noise subspace ˆ U N . The d imension of ˆ U N increases as M increases, and the orthog onality between the signal and noise subspace plays a m ore im portant role than the estimation accuracy of the DO A p arameter in the p o larization p arameter estimation as M increases. Thus th e RMSEs o f the p o larization p arameter a re smaller th an tho se of the o ther two algorith m s as M increases. Howe ver, as shown in the last section , the comp utational co m- plexity of the QDR-MUSIC is ab out LM 2 J 3 J 4 − 6 M 2 J 1 J 2 ﬂops larger than that of the DR-MUSIC and the QNC-MUSIC. This condition w ould be e ven worse as the source number increases. In the second test as shown in Fig. 3 , the numb er of EMVSs of the BS in the x -directio n and the y -d ir ection are M x = 8 and M y = 8 , re sp ecti vely , and hence M = 64 . The snap shot number is N = 200 . T he R MSEs of the e stimated DO A and po la r ization par ameter versus the average received SNR from each mob ile terminal are de picted in Fig. 3 . It can be observed th at the RMSEs of these estimated par ameters o f the thr ee algorithms decrease as SNR inc reases. For the DOA estimation, the RMSEs of the DR-MUSIC are smaller than those of th e ESPRIT , the Q-M U SI C and the QDR-MUSIC, and QNC-MUSIC o utperfor ms the o ther alg o rithms. The r eason is identical with the ﬁrst test. For the polarization p arameter estimation, the RMSEs of th e ESPRIT , th e Q-MUSIC and the QDR-MUSIC are smaller than those of the DR-MUSIC and the QNC-MUSIC when the SNR is less than 0 dB. This is cau sed by the fact th at the orth ogonality between        $YHUDJHGUHFHLYHG615G%       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        $YHUDJHGUHFHLYHG615G%         (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%         $YHUDJHGUHFHLYHG615G%         (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        $YHUDJHGUHFHLYHG615G%         (635,7 4086,& 4'5086,& '5086,& 41&086,& &5% 506(RI ș HVWLPDWHVGHJUHHV 506(RIĳHVWLPDWHVGHJUHHV 506(RIȖHVWLPDWHVGHJUHHV 506(RIȘHVWLPDWHVGHJUHHV D E F G Fig. 3. RMSE s versus the avera ge recei ved SNR from each mobile terminal for the estimates of DOA and polarizati on parameters when using differe nt estimati on algorithms, while the number of the BS antennas is M = 64 , and the snapshot number is N = 500 . (a), (b), (c) and (d) correspond to the estimation of azimuth, ele vation , polarizat ion angle and phase dif ference, respect i vely . the signal and noise subsp aces p lays a more importa n t r ole than th e estima tio n accuracy of the DO A parameter whe n th e SNR is less than 0 dB. The RMSEs of the ESPRIT , the Q- MUSIC and the QDR-MUSIC are larger than those of the DR-MUSIC an d th e QNC-MUSIC when the SNR is larger than 0 dB. T he reason is that the estimation accuracy of the DO A par ameter is very high at a high SNR, and the effect of the or thogon ality between the sign al and noise sub spaces is much weaker . These results d emonstrate that th e e stima tion accuracy of the polar ization par ameter deterior ates when the power of the received n oise is h igh. It should be noted tha t th e RMSEs of th e p o larization param eter beco me smaller as the number of the BS antennas inc reases. The r efore, the propo sed DR-MUSIC and QNC-M U SI C can poten tially achieve good perfor mance by employing a larger num ber of EMVSs. In other words, for th e p olarized massive MI MO systems th e transmitted power can be sign iﬁcantly reduced due to the application of an unp recedented large number o f EM VS at the BS. In the thir d te st as sh own in Fig. 4 , the nu mber of EMVSs of the BS in the x -directio n an d the y - d irection are M x = 8 and M y = 8 , respectively , and hence M = 64 . The average received SNR fro m each mo bile terminal is 5 dB. Th e RMSEs of the estimated DO A a n d p olarization parameter versus the snapshot number N re c e i ved at the BS are depicted in Fig. 4 . It can be observed that the RMSEs of th ese estimated parameters of the three a lg orithms decrease as the snap shot number increases. The curves sh own in Fig. 4 are ﬂatter than those shown in Fig. 3 , which means that f or the par ameter estimation performan c e, the ef fect of the snapshot nu mber is slightly less than that of the SNR. For the polarization parameter estimation , th e RMSEs of the ESPRIT , the Q- MUSIC and the QDR-MUSIC ar e smaller than tho se of the DR-MUSIC and the QNC-MUSIC when the snapshot n umber is less than 4 00; the RMSEs of the ESPRIT , the Q-MUSIC and the QDR-M U SI C are larger tha n those of the DR-MUSI C and the QNC-MUSIC when snapshot number is larger than 400. The effect of increasing the sn apshot nu mber is similar 12        6QDSVKRWQXPEHU       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        6QDSVKRWQXPEHU       (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        6QDSVKRWQXPEHU     (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        6QDSVKRWQXPEHU     (635,7 4086,& 4'5086,& '5086,& 41&086,& &5% 506(RI ș HVWLPDWHVGHJUHHV 506(RIĳHVWLPDWHVGHJUHHV 506(RIȖHVWLPDWHVGHJUHHV 506(RIȘHVWLPDWHVGHJUHHV D E F G Fig. 4. RMSEs versus the s napshot number N rec ei ved at the BS for the estimate s of DOA and polarizati on parameters when using dif ferent estimation algorit hm s, while the number of the BS antenn as is M = 64 , and the av erage recei ved SNR from each mobile terminal is 5 dB. (a), (b), (c) and (d) correspond to the estimation of azimuth, elev ation, polarizat ion angle and phase diffe rence, respecti vely . to the effect o f increasing the SNR, as co mpared with Fig. 3 . In th e fourth test as shown in Fig. 5 , some of th e par a m eters are changed for e valuating the perf ormance of the three algo- rithms with the incre ased numb er of mo bile term inals. A more realistic scenario considerin g the m ultipath pro pagation in 3D millimeter wa ve channels is used fo r si mulation. F or each mobile terminal, the p olarization par ameter is not chan g ed, but the DOA parameter s are different. Ther e are 3 do minant paths from each mobile terminal to the BS, and the sig nal among them a r e coheren t with each o ther . For the ﬁrst mobile terminal, the polarization angle is γ 1 = 40 ◦ , and the phase difference is η 1 = 20 ◦ ; the co r respondin g azimuth DO As are θ 11 = 20 ◦ , θ 12 = 30 ◦ , θ 13 = 50 ◦ , an d th e co rrespond ing elev ation DO As are ϕ 11 = 1 0 ◦ , ϕ 12 = 1 5 ◦ , ϕ 13 = 20 ◦ . For the second mob ile terminal, the p o larization ang le is γ 2 = 1 0 ◦ , and the phase d ifference is η 2 = 60 ◦ ; the correspon ding azimuth DOAs a r e θ 21 = 2 5 ◦ , θ 22 = 6 0 ◦ , θ 23 = 1 25 ◦ , and the corresp o nding elev ation DO As are ϕ 21 = 15 ◦ , ϕ 22 = 5 0 ◦ , ϕ 23 = 40 ◦ . For the third mobile terminal, the polarization angle is γ 3 = 25 ◦ , an d th e p h ase difference is η 3 = 30 ◦ ; the corr espondin g azimuth DOAs are θ 31 = 40 ◦ , θ 32 = 80 ◦ , θ 33 = 110 ◦ , and the correspond ing elev ation DOAs are ϕ 31 = 25 ◦ , ϕ 32 = 60 ◦ , ϕ 33 = 35 ◦ . For the fou rth mobile terminal, the polarization angle is γ 4 = 70 ◦ , and the phase difference is η 4 = 40 ◦ ; the co r respondin g azimuth DO As are θ 41 = 7 0 ◦ , θ 42 = 1 25 ◦ , θ 43 = 140 ◦ , and the c orrespon ding elev ation DOAs ar e ϕ 41 = 30 ◦ , ϕ 42 = 45 ◦ , ϕ 43 = 6 0 ◦ . For the ﬁfth mobile termin al, th e po larization angle is γ 5 = 80 ◦ , an d the pha se difference is η 5 = 50 ◦ ; the correspo nding a zimuth DO As are θ 51 = 75 ◦ , θ 52 = 13 0 ◦ , θ 53 = 150 ◦ , and the correspo n ding ele vation DO As are ϕ 51 = 20 ◦ , ϕ 52 = 5 5 ◦ , ϕ 53 = 7 5 ◦ . The numb er of EM VSs of the BS in the x -direction and the y -dire c tion are M x = 10 and M y = 10 , respectively , and hence M = 1 00 . The average r eceiv ed SNR fro m each mo bile terminal is 0 d B, and the snapsh o t number is N = 500 . The 2-D spatial smoo thing tech n ique has b een u sed for solving the coheren t signa ls, and the size o f sub array is 8 × 8 . The RMSEs          1XPEHURIPRELOHWHUPLQDOV     (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%          1XPEHURIPRELOHWHUPLQDOV     (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%        1XPEHURIPRELOHWHUPLQDOV     (635,7 4086,& 4'5086,& '5086,& &5%   1XPEHURIPRELOHWHUPLQDOV   (635,7 4086,& 4'5086,& '5086,& 41&086,& &5%   41&086,&          506(RI ș HVWLPDWHVGHJUHHV 506(RIĳHVWLPDWHVGHJUHHV 506(RIȖHVWLPDWHVGHJUHHV 506(RIȘHVWLPDWHVGHJUHHV D E F G Fig. 5. RMSEs versus the number of the mobile terminals for the estimates of DOA and polarizatio n parameters when using dif ferent estimation algorithms, while the number of the BS antenna s is M = 100 . T he ave rage recei ved SNR from each m obile terminal is 0 dB, and the snapshot number is N = 500 . (a), (b), (c) and (d) correspond to the estimation of azimuth, elev ation, polarization angle and phase diffe rence, respecti vely . of the estimated DO A a n d p olarization parameter versus the number o f the mobile ter m inals are d epicted in Fig. 5 . It can be observed that the RMSEs of the ﬁv e algorithms except ESPRIT incr ease slowly as the n umber of the mobile term inal increases. This is because a rough estimation of DO A pa- rameter and DO A, polariz a tion par ameter have been o btained for the DR-MUSIC, QNC-MUSIC and the Q-MUSIC, the QDR-MUSIC, respectively , based o n a large step size o f the spectrum sear ch. For the Q-MUSIC a nd th e QDR-MUSIC, the DO A and po larization param eter of the mobile ter m inals are only estimated by searching around the true v alu es. For the DR-MUSIC an d the QN C-M USIC, the DO A parameters of the mob ile terminals are on ly e stima te d by sear ching around the true values, while the polarizatio n para meter ha s a clo sed form. Th us the RMSEs o f th ese four algorith ms ch ange slowly . For the E SPRIT , the RMSEs increases gre a tly compare d with the other four algo rithms. Th is is mainly becau se tha t the array aperture is not utilized comp letely . Some m obile ter minals are too close in space, the ESPRIT cannot distinguish them exactly , thu s the RMSEs increases as the num ber of m o bile terminals incr eases. I X . C O N C L U S I O N In th is paper, we have pro posed a MUSIC-b a sed algo rithm, QNC-MUSIC, fo r the 2-D DO A and polarization estimation of the n on-circu lar sig n al in th e 3-D millimeter wa ve polarized massi ve MIMO s ystems. For th e DO A estimation of non- circular sig n als using the QNC-M USIC, the prop erty of non- circular signals is used to fu rther improve the DO A esti- mation accuracy . For the polarization estimation of th e non - circular signa l, the closed-for m expression of the DR-M USI C is ado pted based o n the DOA estimation result o f the QNC- MUSIC. Compa red with the trad itional L V -MUSIC, Q- MUSIC and th e QDR-MUSIC, the com putational complexity of the QNC-MUSIC and the DR-MUSIC are m uch lower with the help of th e der i vative of the spectru m functio n an d the closed- form expr e ssion of the polar ization estimation. Our an alysis and simulation show that the p erforma n ce o f the pr oposed QNC-MUSIC impr oves as the numb er of the B S an tennas 13 increases, and the DO A estimation accuracy is h igher than other alg orithms in massiv e MIM O systems in particular . In the future, the method of utilizing the pro p erty of non- circular sig n als to improve the estimation perfor mance o f the polarization p arameter shou ld b e stud ied. R E F E R E N C E S [1] M. R. Akdeniz, Y . Liu, M. K. Samimi, S. Sun, S. Rangan, T . S. Rappaport, and Elza E rkip, “ Millimeter wa ve channel modeling and cellular capacity e valuat ion, ” IEEE J . Sel. Areas Commun. , vol. 32, no. 6, pp. 1164–1179 , Jun. 2014. [2] Z. Marzi, D. Ramasamy , and U. Madhow , “Compressi ve channel estima- tion and tracking for large arrays in mm-W ave picocells, ” IEEE J. Sel. T opics Signal Proc ess. , vol. 10, no. 3, pp. 514–527, Apr . 2016. [3] R. Shaﬁn, L . Liu, C. Zhang, and Y . C. W u, “DoA estimati on and capac ity analysis for 3D mill imeter wav e massive -MIMO/FD-MIMO OFDM systems, ” IEEE Tr ans. W irel ess Commun. , vol. 15, no. 10, pp. 6963–6978, Oct. 2016. [4] E. G. L arsson, O. Edfors, F . T ufve sson, and T . L. Marzetta, “Massi ve MIMO for next generation wirele ss systems, ” IEEE Commun. Mag. , vol. 52, no. 2, pp. 186–195, Feb . 2014. [5] L. Cheng, Y . C. Wu, J. Zhang, and L . Liu, “Subspa ce ident iﬁcati on for DOA estimation in massiv e/full-dimension MIMO systems: bad data mitigat ion and auto matic source enumeration, ” IEEE T rans. Signal Pr ocess. , vol. 63, no. 22, pp. 5897–5909, Nov . 15, 2015. [6] H. Y in, D. Gesbert, M. Filippou, and Y . Liu, “ A coordi nated approach to channel estimation in larg e-scale m ultiple-ante nna systems, ” IEEE J. Sel. Areas Commun. , vol. 31, no. 2, pp. 264–273, Feb. 2013. [7] S. He, Y . Huang, H. W ang, S. Jin, and L. Y ang, “Leakage-a ware energy - ef ﬁcient beamformin g for het erogene ous multi cell multiuse r systems, ” IEEE J . Sel. A r eas Commun. , vol. 32, no. 6, pp. 1268–1281, Jun. 2014. [8] A. Adhikary , J. Nam, J.-Y . Ahn, and G. Caire, “Joint spatial div ision and multiple xing: T he large-sca le array regime, ” IEEE T rans. Inf. Theory , vol. 59, no. 10, pp. 6441–6463, Oct. 2013. [9] A. Alkhateeb, O. El A yach, G. Leus, and R. W . Heath, J r., “Channel estimati on and hybrid precoding for millimeter wav e cell ular systems, ” IEEE J . Sel. T opics Signal. Pr ocess. , vol. 8, no. 5, pp. 831–846, Oct. 2014. [10] R. W . Heath, Jr ., N. Gonz ´ alez-Pr elcic , S. Rangan, W . Roh, and A. M. Sayeed, “ An ov ervie w of signal processing techni ques for millimeter wave MIMO systems, ” IEEE J. Sel. T opics Signal Pr ocess. , vol. 10, no. 3, pp. 436–453, Oct. 2016. [11] J. G. Andrews, T . Bai, M. N. Kulka rni, A. Alkhateeb, A. K. Gupta, and R. W . Heath, Jr ., “Modeling and analyzing millimet er wa ve cellul ar systems, ” IEEE Tr ans. Commun. , vol. 65, no. 1, pp. 403C-430, Jan. 2017. [12] P . Gounon, C. Adnet and J. Galy , “Localisa tion angulai re de signaux non circula ires, ” T rait. Signal , vol. 15, no. 1, pp. 17–23, 1998. [13] H. Abeida and J. P . Delmas, “MUSIC-like estimat ion of direc tion of arri val for noncircu lar sources, ” IEEE Tr ans. Signal P r ocess. , vol. 54, no. 7, pp. 2678–2690 , Jul. 2006. [14] J. P . Delmas, “ Asymptotically m inimum varia nce second-order estima- tion for noncircul ar signals with applicat ion to DOA estimati on, ” IEEE T rans. Signal Pr ocess. , vol. 52, no. 5, pp. 1235–1241, May 2004. [15] P . Charge , Y . W ang and J. Saillard, “ A root-MUSIC algorithm for non circul ar sources, ” in P r oc. IEEE Int. Confer ence on A coustic s, Speech, and Signal Proc essing (ICASSP) , Salt Lake City , USA, May . 2001. [16] A. L iu, G. L iao, Q. Xu and C. Zeng, “ A circulari ty-based DOA estima- tion method under coexist ence of noncircul ar and circular signals, ” in Pr oc. IE EE Int. Confer ence on Acoustics, Speech , and Signal Pro cessing (ICASSP) ,Kyo to, Japan, Mar . 2012. [17] F . Gao, A. Nallanat han and Y . W ang, “Improv ed MUSIC under the coe xistenc e of both circula r and noncircular sources, ” IEEE T rans. Signal Pr ocess. , vol. 56, no. 7, pp. 3033–3038, Jul. 2008. [18] X. Gou, Z. Liu, and Y . Xu, “Biquaterni on cumulant-MUSIC for DOA estimati on of noncircul ar signals, ” Signal Pr ocess. , vol. 93, no. 4, pp. 874–881, Apr . 2013. [19] Z. M. Liu, Z. T . Huang, Y . Y . Zhou, and J. Liu, “Direction -of-arri val estimati on of noncirc ular signals via sparse representati on, ” IEEE T rans. Aer osp. E lect r on. Syst. , vol. 48, no. 3, pp. 2690–2698, J ul. 2012. [20] A. Zoubir , P . Charg´ e, and Y . W ang, “Non circula r sources localizati on with ESPRIT , ” in E ur . Conf . W irel ess T echn ol. (ECWT) , Munich, Germany , O ct. 2003. [21] M. Haardt and F . Roemer , “ Enhancemen ts of unitary ESPRIT for noncirc ular sourc es, ” in IEEE Int. Conf. A coust., Speec h, Signal Pr ocessing (ICASSP) , Montreal, QC, Canada, May 2004. [22] F . Roemer , M. Haardt, and G. Del Galdo, “ Analytic al performance assessment of multi-di mensional matrix- and tensor-b ased ESPRIT -type algorit hm s, ” IEE E T rans. Signa l P r ocess. , vol. 62, no. 10, pp. 2611–2625 , May 2014. [23] J. Steinwandt , F . Roemer , M. Haardt, and G. Del Galdo, “R dimensional ESPRIT -type algorith ms for strictly second-or der noncircular sources and their performance analysis, ” IEEE T rans. Signal Proce ss. , vol. 62, no. 18, pp. 4824–4838, Sep. 2014. [24] J. Steinwandt , F . Roemer , and M. Haardt, “ESPRIT -T ype Algorithms for a Recei ved Mixture of Circular and S trictly Non-Circu lar Signals, ” in Pr oc. IEEE Int. Conf. A coustic s , Speech and Sig. P r oc. (ICASSP 2015) , Brisbane , Australia , Apr . 2015. [25] F . Roemer and M. Haardt, “Efﬁcie nt 1-D and 2-D DOA estimation for non-circ ular sources with hexagonal shaped espar arrays, ” in Pr oc. IEEE Int. Confere nce on A coustic s, Speech, and Signal Proce s sing (ICASSP) , T oulouse, France, pp. 881–884, May 2006. [26] J. L i and R. T . C. J r, “ Angle and polarizatio n estimation using ESPRIT with a polarizati on sensiti ve array , ” IEEE T rans. Antennas Propa g. , vol. 39, no. 9, pp. 1376–1383, Sep. 1991. [27] K. T . W ong and M. D. Zolto wski, “Di versely polarized root-MUSIC for azimuth- ele vati on angle of arri v al estimati on, ” Dig. 1996 IEEE Antennas P r opagation Soc. Int. Symp. , pp. 1352–1355, Sep. 1996. [28] E. Gonen and J. M. Mendel, “ Applicati ons of cumulants to array processing. Part VI. Polarizati on and direction of arri v al estimation with minimally constraine d arrays, ” IEE E T rans. Signal Pro cess. , vol. 47, no. 9, pp. 2589–2592, Sep. 1999. [29] P . Chev alier , A. Ferr ´ eol, L. Albera, and Gw ´ ena ¨ el Birot, “Higher order direct ion ﬁnding from arrays with div ersely polarized antenna s : the PD- 2q-MUSIC algorit hms, ” IEEE T rans. Signal Pr ocess. , vol. 55, no. 11, pp. 5337–5350, Nov . 2007. [30] X. Liu and N. D. Sidiropo ulos, “Cram ´ er -Rao lower bounds for lo w - rank decomposition of multidimensiona l arrays, ” IEEE T rans. Signal Pr ocess. , vol. 49, no. 9, pp. 2074–2086, Sep. 2001. [31] X. Guo, S. Miron, D. Bri e, S. Zhu, and X. Liao, “ A CANDE- COMP/P ARAF AC perspec ti ve on uniqueness of DOA estimation using a vec tor s ensor array , ” IEEE T rans. Signal P r ocess. , vol. 59, no. 7, pp. 3475–3481, Jul. 2011. [32] M. Costa, A. Richter , and V . Koi vunen, “DoA and polarizat ion estima- tion for arbit rary array conﬁgura tions, ” IEE E T rans. Signal P r ocess. , vol. 60, no. 5, pp. 2330–2343, May 2012. [33] Y . Tia n, X. Sun, S. Zhao, “Sparse -reconst ruction -based directi on of arri val, polarisat ion and po wer esti m ation using a cross-dipole array , ” IET Radar Sonar Navig. , vol. 9, no. 6, pp. 727–731, Jul. 2015. [34] S. Miron, N. L . Bihan, and J. I. Mars, “Quaternion-M USIC for vector - sensor array processing, ” IEEE T rans. Signal Proce ss. , vol. 54, no. 4, pp. 1218–1229, Apr . 2006. [35] N. L. Bihan, S. Miron, and J. I. Mars, “MUSIC algorithm for vector - sensors array using biquat ernions, ” IEEE T rans. Signal Pr ocess. , vol. 55, no. 9, pp. 4523–4533, Sep. 2013. [36] X. Gong, Z. W . Liu, and Y . G. Xu, “Direction ﬁnding via biquaternion matrix diagonali zation with vec tor-sensors, ” Signal Pr ocess. , vol. 91, no. 4, pp. 821–831, Apr . 2011. [37] A. S. Y . Poon and D. N. C. T se, “Degree-of-fre edom gain from using polarime tric antenna elements, ” IEEE T rans. Inf. Theory , vol. 57, no. 9, pp. 5695–5709, Sep. 2011. [38] X. Su, D. Choi, X. Liu, and B Peng, “Cha nnel Model for Polarize d MIMO Systems With Power Radiation Patt ern Concern, ” IEEE Access , vol. 4, pp. 1061–1072, Mar . 2016. [39] J. Li and J. T ao, “The dimension reducti on quaternion MUS IC algorithm for jointly estimating DOA and polarizati on, ” J. Electr on. Inf. T echnol. , vol. 33, no. 1, pp. 106–111, Jan. 2011. [40] W . Si, T . Zhu, and M. Zhang, “Dimension-re ductio n MUSIC for jointly estimati ng DOA and polariza tion usin g plane polariz ed arrays, ” J . Commun. , vol. 35, no. 12, pp. 28–35, Dec. 2014. [41] N. L. Bihan and J . Mars, “Singular val ue decomposit ion of quaterni on matrice s: a new tool for ve ctor -s ensor signal processing, ” Signal Pr ocess. , vol. 84, no. 7, pp. 1177–1199, Jul. 2004. [42] J. P . W ard, Quaternions and Cayley Numbers, Algebra and Applica- tions . Norwell, MA: Kluwer , 1997. [43] F . Zhang, “Quatern ions and matrices of quaternions, ” Linear Algebr a Its A ppl. , vol. 251, pp. 21–57, Jan. 1997. [44] R. Schmidt, “Multiple emitte r location and signal parameter estimation, ” IEEE T rans. Antennas Propa g. , vol. 34, no. 3, pp. 276–280, Mar . 1986. [45] G. Xu and T . Kailath, “Fa s t subspace decomposition, ” IEE E T rans. Signal Pr ocess. , vol. 42, no. 3, pp. 539–551, Mar . 1994. 14 [46] P . Sto ica and A. Nehorai, “Performance study of co ndition al and uncondit ional direction- of-arri val estimation, ” IEEE T rans. Acoust., Speec h, Signal P r ocess. , vol. 38, no. 10, pp. 1783–1795, O ct. 1990. [47] P . Stoica, E. G. Larsson and A. B. Gershman, “The stochastic CRB for array proce s sing: a textboo k deriv ation, ” IEEE Signal Pro cess. Lett . , vol. 8, no. 5, pp. 148–150, May 2001. [48] H. Abeida and J. P . Delmas, “Stoch astic Cram ´ er -Rao Bound for noncir - cular signals w ith applicat ion to DOA estimatio n, ” IEEE T rans. Signal Pr ocess. , vol. 52, no. 11, pp. 3192–3199, Nov . 2004. [49] S. B. Hassen, F . Bellili, A. Samet and S. Affes, “DOA estimatio n of temporal ly and spatiall y corre lated narrowban d noncircular sources in spatial ly correla ted white noise, ” IE EE Tr ans. Signal Pr ocess. , vol. 59, no. 9, pp. 4108–4121, Sep. 2011.

DOA and Polarization Estimation for Non-Circular Signals in 3-D Millimeter Wave Polarized Massive MIMO Systems

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment