# Polarimetric Direction of Arrival Estimations Based on Adaptive Linear Time-Frequency Transforms

Authors: Shao Shuai; Liu Aijun; Wang Xiuhong; Yang Hongjuan
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

A spatially polarized time-frequency distribution (SPTFD) based on dual-polarized double-fed antenna arrays is adapted to deal with polarization-unstable signals. A linear time-frequency (TF) representation was used for an instantaneous frequency (IF) estimate, primarily due to its simplicity and immunity to cross-interference. Using a set of linear TF transformations using Gaussian windows and Fourier oscillation kernels, the IF estimated window widths of multiple unstable signals are obtained. This paper introduces a new method for estimating the direction of arrival (DOA) of polarized waves using adaptive linear time-frequency transforms. In this paper, a narrowband far-field point source on the receiving array is analyzed. The source signal is split into two orthogonally polarized components. The optimal window is determined by the first derivative of the IF; for this purpose, we take a simple algorithm for solving the derivative and optimize it. In developing TF-adaptive and fully automatic TF display technology, the first method is to use the time-adaptive window for minimizing the IF estimate mean square error (MSE) sum at each moment, while the second procedure is to adjust according to time and frequency and minimize estimate MSE sum at each position in the TF region. Due to its combination with signal polarization, the spatial time-frequency distribution (STFD) gains more freedom and thus perfects the phonon space estimation of noise and signal. On the SPTFD platform, polarized time-frequency multiple signal classification (PTF-MUSIC) is used for the estimation of signal direction of arrival, which outperforms conventional time-frequency MUSIC. Using the example of a synthesized signal, this method outperforms conventional techniques in DOA estimation.


---

## Body

## 1. Introduction

In many aspects of technology, we are studying changes in the spectrum. Descriptive information for such a signal can be obtained from a time-frequency representation (TFR), which projects the signal into the time-frequency region, exposing momentary changes in its spectral content [1]. In recent decades, various TFRs have been developed and used in a wide range of applications [2, 3]. These techniques have taken advantage of signal sparseness inside TF regions [4]. Time-frequency analysis is an effective method for DOA estimation of unstable signals. To improve the performance of DOA, [5–7] propose new DOA estimation methods.In TFR, the linear transform can be viewed as a window to the Fourier transform; prominent examples of this type are the short-time Fourier transform (STFT) and the S transform (ST) [8]. Cohen’s TF distribution is based on a generalized instantaneous autocorrelation function. During the analysis, these signals can generate unwanted error signals, which can be hidden, so it is necessary to have a suitable core [9]. The instantaneous frequency (IF) is the most fundamental concept in TF, reflecting the instantaneous change rate of the signal phase function, a simple and well-known instantaneous frequency estimation method by the maxima location of the TFR [10]. To accurately express the frequency law of the signal, TFR is generally required according to the design of the signal because no TFR can optimize all signals like all fixed graphs. Therefore, much research has been devoted to adaptive TFRs whose parameters can be changed with the signal without being disturbed by the user. In this case, energy concentration measurement (ECM) is widely used to determine the energy distribution over the signal components to automatically select good TFR criteria [11, 12]. A multiview approach to adaptive TFR is to build multiple TFRs and combine them according to predetermined criteria [13–15].The accuracy of IF is a primary measure of TF, so we can use it to measure the statistical performance of IF to get a more accurate TFR. In several works [16, 17], the accuracy of IF estimation for TFR was investigated, where the accuracy of the estimate depends on the width of the window and optimizing for this requires optimizing for unknown first derivatives of information. Because of this limitation, this paper uses the intersection of confidence intervals (ICI) to approximate the optimal width in [18]. In [19], a three-stage adaptive algorithm was proposed to reduce the computational complexity required to obtain the linear TFR improved in [20]. However, the optimal width depends on the statistical performance of the IF estimation. [21] utilized a new general formula that describes the error and variation of IF estimation in noisy environments, resulting in the optimal window width for linear TFR. Meanwhile, in [22, 23], polarization information is adopted, which further improves the estimation performance of DOA. Then, an adaptive method based on the combination of polarization and time-frequency domain information is proposed.In Materials and Methods, we will discuss the time-frequency distribution patterns of spatial polarization and two automation linear TFRs. The one utilizes a time-adaptive window for minimizing the MSE sum values. In another TFR, a TF-adaptive window is used to minimize the location of the estimated MSE in the TF region. The part also introduces the polarized time-frequency MUSIC. The method in this paper is evaluated in the Results and Discussion. Through numerical calculation, it is proved that the procedure is superior to general and complex algorithms.


## 2. Materials and Methods

### 2.1. Spatial Polarimetric Time-Frequency Distributions

In Figure1, a electromagnetic wave incident into the array can be described as
(1)Et=Eθtθ+Eϕtϕ=Eθtcosθcosϕ−Eϕtsinϕx+Eθtcosθsinϕ+Eϕtcosϕy+Eθtsinθz,where ϕ and θ are the spherical azimuth and elevation unit vectors observed by the source. The x, y, and z are unit vectors along the x, y, and z directions, respectively. For generality and simplicity, we assumed that the signal is located in the x−y plane and the array is in the y−z plane. Accordingly θ=90° (θ=−z),
(2)Et=−Eϕtsinϕx+Eϕtcosϕy+Eθtz.Figure 1
Dual-polarized array.We expressst in terms of the amplitude of the source measured by the receiving reference sensor, which has a polarization angle γ∈0,π/2 and a polarization phase difference η∈−π,π. svt and sht are the source horizontal and vertical polarization part components, which can be described by spherical fields, Eθt and Eϕt, as
(3)Eθt=svt=stcosγ,Eϕt=sht=stsinγejη.A signal is linearly polarized ifη=0° or η=180°. Replacing (3) into (2) gets
(4)Et=st−cosγsinϕx+cosϕsinγejηy+cosγz.We assume thatN signals enter on the array of M dual-polarized antennas. The distance between sensors is d, and d is not more than λ/2 for ensuring that there is no array ambiguity. In another word, the expected direction of the incident signal is only one. When d is λ/2, there are high array sensor utilization rate and effect guarantee. The nth source vertical and horizontal components are
(5)snvt=sntcosγn≜cn1snt,snht=sntsinγnejηn≜cn2snt,where the parameters cn1=cosγn and cn2=sinγnejηn represent the vertical and horizontal polarization coefficients. The signal entering on the mth dual-polarized antenna is
(6)y¯mt=ymvt,ymhtT=∑n=1NanmvEn⋅z,anmhEn⋅yT=∑n=1Nanmvsnvt,anmhsnhtcosϕnT,where “⋅” denotes the dot product, En represents the nth source electric field vector, and anmv and anmh, respectively, denote the mth elements of the vertically and horizontally polarized array vectors, avϕn and ahϕn. We assume that both avϕ and ahϕ are known and normalized (avϕ2=ahϕ2=M) and the array has been calibrated. It should be pointed out that in the array calibration of the relevant area, the cosϕn term in the horizontally polarized array manifold can be introduced, so it will not be considered. Then, the above equation is simplified as
(7)y¯mt=anmvsnvt,anmhsnhtT=sntanmvanmhT⊙cn1cn2T≜sntanm⊙cn,where the vector cn=cn1,cn2T=cosγn,sinγnejηnT is the nth source polarization signature.For a dual-polarized sensor, we define the general short Fourier transform (GSFT) as(8)Dxit,f=∫−∞+∞wt−τ,σt,fxite−j2πfτdτ,where the Gaussian window is defined as
(9)wt,σt,f=12πσt,fexp−t22σ2t,f.Formulas (1)–(9) correspond to the single dual-polarization sensor case. With an M-sensor array, for each polarization i, i=vorh, the data vector is
(10)xit=x1it,x2it,⋯,xMitT=yit+nit=AiΦsit+nit.Using formula (10), the case of extending the polarization time-frequency distribution of a single sensor to multiple sensors is obtained. Instead of the scalar TFD of (8), we define the polarization STFD matrix of the data vector
(11)Dxit,f=∫−∞+∞wt−τ,σt,fxite−j2πfτdτ,which, for the no-noise case, can be represented as
(12)Dxit,f=AiΦDsit,fAiΦH.According to (10), the double polarization data vectors can be built as follows:
(13)xt=xvtxht=AvΦ00AhΦsvtsht+nvtnht=AvΦ00AhΦQvQhst+nvtnht=BΦQst+nt,where
(14)BΦ=AvΦ00AhΦ,represents block diagonal, and
(15)Q=QvQh,represents the signal polarization characteristic vector, where
(16)qv=cosγ1,⋯,cosγNT,Qv=diagqv,(17)qh=sinγ1ejη1,⋯,sinγNejηNT,Qh=diagqh.Correspondingly,(18)BΦQ=avϕ1cosγ1⋯avϕncosγNahϕ1sinγ1ejη1⋯ahϕnsinγNejηN=a~ϕ1⋯a~ϕN.The matrix (18) can be considered as the extended mixing matrix, within a~ϕn denoting the nth signal joint spatial polarization characteristic. For the nth signal, the extended spatial polarization characteristic vector is
(19)a~ϕn=avϕncosγnahϕnsinγnejηn.Obviously, the double polarization array can double the spatial dimension of the vector.The polarization, spatiality, and time-frequency characteristics can be combined with the source signal incident on the receiver array. The STFD of a dual-polarization data vectorxt can be as follows:
(20)Dxt,f=∫−∞+∞wt−τ,σt,fxte−j2πfτdτ.Dxt,f is called the SPTFD matrix. This matrix is a general method that can solve some common problems in array processing, as described in the following.When the noise effect is neglectful, the SPTFD matrix is connected with the signal TFD matrix by(21)Dxxt,f=BΦQDst,fQHBHΦ.Firstly, the time-frequency points are selected by the algorithm in the next section to construct the matrix. Then, eigenvectors are obtained by eigendecomposition of formula (21). Finally, the noise subspace Un constructed by eigenvectors is used for subsequent spatial spectrum estimation.


### 2.2. Three-Step Method for IF Estimation

Abdoush et al. [21] gave a general expression relating the error and variance of the IF estimate to the smoothness of the observation window and the IF. The correct IF estimation demands a data-based approach to choose at least reasonable window width. Therefore, the paper summarizes and extends the “low complexity” described in [19, 20]; the method includes three processes, as shown in Figure 2, which will be explained in this section. On this basis, this paper proposes a new linear TFR method based on estimation error and variance. The one uses a time-varying window for minimizing the MSE sum value estimated by IF at each moment, while in the latter TFR, the TF value is adjustable to minimize the MSE of each point in the TF region.Figure 2
Flowchart of adaptive TFRs based on linear TF transforms.For practical application, we will limit ourselves to a discrete version of GSFT, which we define as(22)Dxkt,kf=∑n=0N−1xlwkt−l,σkt,kfe−j2πkfl/N,where xkt represents a length Kt discrete time series corresponding to xktTs by sampling interval Ts and kf=−Kt/2,⋯,Kt/2−1 denoting frequency index.


#### 2.2.1. Preliminary TFR

For any TFR, the optimized formula (22) requires an understanding of the IF. So, in adaptation, the first step is to create a preliminary TFR (PTFR) that roughly estimates the IF orbit. Although this stage can tolerate minor estimation errors, however, more significant estimation errors are propagated to the after step, so the final IF estimation accuracy decreases. To this end, Pei and Huang in [20] adopted an STFT of optimizing window width according to the ECM. While it performed well, when we used different multicomponent signals to detect this ECM, we found that the approach was beneficial for a noisy element and detrimental for a rapidly changing part. That is, in TFR, only quasistable components are selected. In this paper, we design an STFT based on the standard deviation of conventional rates based on PTFR. The rate is defined as the ratio of the effectivity bandwidth Beff to the effectivity duration Teff of the signal being analyzed.
(23)ν0=BeffTeff=1N∑k=−Kt/2Kt/2−1k−k02Xk2∑kt=0Kt−1kt−kt02xkt21/2,where Xk represents the discrete Fourier transform of xkt:
(24)k0=∑k=−Kt/2Kt/2−1kXk2∑k=−Kt/2N/2−1Xk2,kt0=∑kt=0Kt−1ktxkt2∑kt=0Kt−1xkt2.If the normalized power spectral density (that is, the normalized unit) of a signal is regarded as a function of probability mass, the second central moment of the signalBeff2 is the second central moment of its frequency, Beff is its standard deviation, and it represents the bandwidth of the signal. We can think of Teff as the time width of a signal, as before. In this way, the conventional ratio ν0 can be understood as the linear chirp scanning the frequency range Beff throughout Teff. We use the regular rate to determine the overall constant change rate of the signal IF regularity, which can be used as a fixed observation window for the PTFR to set the standard deviation of the STFT according to σ=3/71/4/2πν0.


#### 2.2.2. IF Rate Estimation

In stage 2 of the adaptive algorithm, we get a rough estimate of IFR (IFR: IF rate, that is first derivative of IF). For this end, the IF is estimated by the main trajectory of PTFR. On this basis, the IFR is obtained by the IF first derivative. If the signal is multicomponent, each component must be found before differentiating each IF ridge. Here, the IF component separation is implemented according to the image processing algorithm described in [24]. It involves converting the TFR to a binary image and then performing a component connection process, in the case of TFR longer than a predetermined threshold, according to the adjacent connection criterion, extracting the connection peak from the TFR. This approach does not require prior knowledge of the relevant components but instead assumes that the IF ridge can be distinguishable in the TFR, which are many real-world signals, including most signals [25–27].For the derivation, [19] adopted the conventional difference operator. However, we know that this method is very sensitive to identifying discrete signals with noise. Therefore, the signal-to-noise ratio of the differential signal is usually lower than that of the original signal. Thus, in [20], the difference operator is replaced by principal component analysis, which executes eigenvalue decomposition for the 2 L estimations neighbouring f^Ikt to approximate fI′kt. In this article, we use a more straightforward way to find first-order differentials. First, the estimates f^Ikt are smoothed according to
(25)f^~Ikt=12L+1∑kt−Lkt+Lf^Ikt,kt=L,⋯,Kt−L,where L represents a positive integer. The above smoother represents a simple average moving filter used to improve the SNR of the differential signal [28]. After above smoothing the signal, the IFR is estimated using a traditional difference operator. As L increases, the degree of noise reduction increases, but the attenuation of the differential amplitude also becomes more remarkable. In addition, filtering the IF estimates according to (25) can help exclude possible outliers before derivation. Still, the accuracy of this method is limited because there is a certain relation between consecutive sampling errors in TFR [10]. Therefore, be careful when choosing L. Here, the smoothing width is 11.


#### 2.2.3. Optimizing Window Width

After extracting the IF and IFR, they can be optimized for appropriate window widths. At this step, two ways adjust the TFR window width. In one case, the window width can be set in terms of time or frequency. For all signals, the primary purpose of this width optimization is to perform an optimal tradeoff at each moment or frequency. Therefore, when the window width is optimal, its cost function is the smallest. Since both adaptivity STFT (ASTFT) and ST can be done rapidly with FFT, the advantage of this method is that its complexity is slight. On the other hand, you can set a window to change at different times and frequencies. Within the TF region, width adjustment is allowed, and the second method requires more flexibility; however, sinceσkt,kf is related to kt and kf, in (22), GSFT cannot be implemented with FFT.We design an ASTFT that minimizes the IF estimation MSE sum at any moment. Because the TFRs generated by ASTFTs generally do not have crossterms, the TFRs can be approximated to the TFRs corresponding to the individual components of the overall signal. In this example, the formula can be generalized to the case of multicomponent signals with TF regions. That is, with ASTFT, an estimate of the IF with no offset can be obtained, and the total MSE can be obtained:(26)MSEtotkt,σkt=∑l=1Ln2πfIl′ktσ2kt2+15/232π5/2σ3ktC0,where Ln represents the signal component number at the moment instant ktTs and fIl′kt represents the lth single IF component IFR. The optimal window width is defined by minimizing MSEtotkt,σkt, which may be approximated by
(27)σoptkt≈argminσkt∑l=1Ln2πf∧Il′ktσ2kt2+15/2σ3kt.Note thatMSEtotkt,σkt represents the summation of all functions, each of which is completely convex and has an overall minimum value on σlkt, and there also is an overall minimum value for MSEtotkt,σkt. Therefore, it is contained in the following range:
(28)σoptkt∈minσ1kt,⋯,σLnkt,maxσ1kt,⋯,σLnkt.A good approximation ofσoptkt is obtained in two steps: first, σlkt is found corresponding to the estimated IFR minimum and maximum; thus, for width between the above two limits, the optimal candidates depend on the objective function form (27). The TFR we obtained is called optimal ASTFT (OASTFT).To reduce the sum of MSE estimated by IF, the width of the frequency domain correlation window in ST is optimized. It is required to obtain the 1st- and 2nd-order frequency derivative information of the window widths. Not knowing the nature of the differential equation, it has to be designed accordingly, making the method impractical; therefore, we cannot go on like this.The previous OASTFT algorithm has achieved a good balance among the signal components in the TF region, but it needs to be adjusted more to make it more adaptable. If the IF ridges can be separated within the TF region, then the window width of the GSFT can be adjusted individually. We can roughly estimate the IFs and IFRs, so the problem is how to defend the standard deviation of a window about the on-ridge and off-ridge time-frequency points. At IF, the first and second derivatives of the standard deviation and frequency are 0, and then, window width is optimal. Therefore, we define a window width intervalΔf determined by the following formula on each upper ridge point of the IF component according to
(29)σkt,kf=3/71/42πf^Il′kt,0≤kf−f^Ilkt≤Δf2,whereas outside the section, the window width is achieved by two-dimensional linear interpolation in [20]. The last interpolation ensures that the σkt,kf 1st- and 2nd-order frequency derivatives are close to 0 neighbouring f^Ilkt, so the width chosen is best for MSE. Also, the coarse IF estimate can be far from its actual value due to noise. Therefore, preferably, not only the width of (29) is specified for f^Ilkt but also for adjacent frequencies. In experiments, we found that several frequency bands are wide enough to obtain accurate IF estimates. The resulting transformation is expressed by optimal GSFT (OGSFT).The adaptation plan adopted consists of the three phases mentioned above. Its main work is to realize PTFR, extract and track IF components, thus obtain the corresponding IFR value, and realize adaptive TFR. The PTFR designed in this paper is STFT, and its computational complexity isON2logN. The chosen IF component extraction and tracking algorithm is presented in [24]; the computational cost for this analysis procedure is outside the scope of this work. The final TFR can be a time-adaptive ASTFT or a TF-adaptive GSFT. OASTFT can be done by FFT. That is, OASTFT needs to operate on the order of ON2logN. The TF-adaptive GSFT algorithm has high adaptability, but because the FFT program cannot complete the conversion, it will lead to ON3 level of computational complexity.


### 2.3. Polarimetric Time-Frequency Music and Spatial Polarimetric Correlations

For improving the signal spatial resolution with well-defined time-frequency characteristics, time-frequency multiple signal classification (TF-MUSIC) has recently been proposed [29]. This method provides an essential summary for applying MUSIC to polarization arrays, which searches for the minimum of the array vector derived from the SPTFD matrix defined in the combined space and polarization domains. The advantages of the algorithm mainly include the following two points: the algorithm integrates polarization information and time-frequency information to increase the utilization rate of multidimensional information. The algorithm adaptively recognizes the time-frequency points of autoterms and avoids the trouble of manual selection.Consider the following spatial signature matrix:(30)Fϕ=1Mavϕ00ahϕ.Becauseaiϕ2=M, FHϕFϕ represents the identity matrix. For searching from the spatial and polarization joint domains, the spatial polarization vector is defined as
(31)fϕ,c=FϕcFϕc=Fϕc,where the vector c=c1c2T represents an unknown unit norm vector of polarization coefficients. In (31), Fϕc=cHFHϕFϕc1/2=cHc1/2=1.The PTF-MUSIC spectrum is provided from the following function:(32)Pϕ=mincfHϕ,cUnUnHfϕ,c−1=minccHFHϕUnUnHFϕc−1,where Un represents the noise subspace from the SPTFD matrix of (20) utilized the selected time-frequency points. For DOA estimation, different STFD or SPTFD matrices can be combined using TF average and joint block diagonalization techniques [29]. Picking these points from regions of high energy concentration associated with global or local sources can improve the signal-to-noise ratio and make the MUSIC algorithm more noise-robust than its traditional MUSIC counterpart [30]. In (32), by acquiring the smallest eigenvalue for the matrix FHϕUnUnHFϕ, this method can perform simple eigendecomposition on the 2×2 matrix, thus avoiding performing a large number of operations in the polarization region. In this way, the spectrum of PTF-MUSIC can be described as
(33)Pϕ=λmin−1FHϕUnUnHFϕ,where λmin⋅ represents the smallest eigenvalue solution. The DOA value of the source is the highest in the spectrum. Corresponding to each angle ϕn of arrival of N signals, n=1,2,⋯,N, each polarized signal parameter can be estimated as follows:
(34)c^ϕk=vminFHϕkUnUnHFϕk,where vmin⋅ represents the eigenvector corresponding to the minimality eigenvalue λmin⋅.The spatial resolution capability of the array has a lot to do with the propagation characteristics of the signal [30]. This depends on the manifold vector normalized inner product of each array. For potential problems involving spatial and polarization dimensions, the extended array manifold a~ϕ is used to define the joint spatial polarization correlation coefficient between the source nl and nk, that is,
(35)βnl,nk=1Ma~Hϕnka~ϕnl=cnk1∗cnl1βnl,nkv+cnk2∗cnl2βnl,nkh=1Mcnk1∗cnl1avϕnkHavϕnl+cnk2∗cnl2ahϕnkHahϕnl,where βnl,nki=1/MaiϕnkHaiϕnl represents the spatial relation coefficient.An interesting case phenomenon occurs with the identical array manifolds between vertically and horizontally polarization, i.e.,avϕ=ahϕ. For the condition, βnl,nkv=βnl,nkh, and the spatial polarization joint relation coefficient is the individual spatial polarimetric correlations product:
(36)βnl,nk=βnl,nkvρnl,nk,within
(37)ρnl,nk=cnkHcnl=cosγnlcosγnkejηnl−ηnk+sinγnlsinγnk,denoting the polarization correlation coefficient. Particularly, for linear polarization, ηnl=ηnk=0, and (37) reduces to
(38)ρnl,nk=cosγnl−γnk.Sinceρnl,nk≤1, the equality holds only when the polarization states of the two light sources are the same; compared with the single spatial correlation coefficient, the spatial polarization correlation coefficient is small. Due to the polarization diversity dropping correlation value, the introduced method ρnl,nk was adopted, which was converted into enhanced signal source identification. In this way, it is difficult to resolve the two sources with a single monopolar spatial array manifold avϕ or ahϕ with an extended spatial polarized array manifold, defined by a~ϕ, which can be easily separated. This improvement is more pronounced when the signal spatial relation coefficient is large, but the individual polarization correlations are low.


## 2.1. Spatial Polarimetric Time-Frequency Distributions

In Figure1, a electromagnetic wave incident into the array can be described as
(1)Et=Eθtθ+Eϕtϕ=Eθtcosθcosϕ−Eϕtsinϕx+Eθtcosθsinϕ+Eϕtcosϕy+Eθtsinθz,where ϕ and θ are the spherical azimuth and elevation unit vectors observed by the source. The x, y, and z are unit vectors along the x, y, and z directions, respectively. For generality and simplicity, we assumed that the signal is located in the x−y plane and the array is in the y−z plane. Accordingly θ=90° (θ=−z),
(2)Et=−Eϕtsinϕx+Eϕtcosϕy+Eθtz.Figure 1
Dual-polarized array.We expressst in terms of the amplitude of the source measured by the receiving reference sensor, which has a polarization angle γ∈0,π/2 and a polarization phase difference η∈−π,π. svt and sht are the source horizontal and vertical polarization part components, which can be described by spherical fields, Eθt and Eϕt, as
(3)Eθt=svt=stcosγ,Eϕt=sht=stsinγejη.A signal is linearly polarized ifη=0° or η=180°. Replacing (3) into (2) gets
(4)Et=st−cosγsinϕx+cosϕsinγejηy+cosγz.We assume thatN signals enter on the array of M dual-polarized antennas. The distance between sensors is d, and d is not more than λ/2 for ensuring that there is no array ambiguity. In another word, the expected direction of the incident signal is only one. When d is λ/2, there are high array sensor utilization rate and effect guarantee. The nth source vertical and horizontal components are
(5)snvt=sntcosγn≜cn1snt,snht=sntsinγnejηn≜cn2snt,where the parameters cn1=cosγn and cn2=sinγnejηn represent the vertical and horizontal polarization coefficients. The signal entering on the mth dual-polarized antenna is
(6)y¯mt=ymvt,ymhtT=∑n=1NanmvEn⋅z,anmhEn⋅yT=∑n=1Nanmvsnvt,anmhsnhtcosϕnT,where “⋅” denotes the dot product, En represents the nth source electric field vector, and anmv and anmh, respectively, denote the mth elements of the vertically and horizontally polarized array vectors, avϕn and ahϕn. We assume that both avϕ and ahϕ are known and normalized (avϕ2=ahϕ2=M) and the array has been calibrated. It should be pointed out that in the array calibration of the relevant area, the cosϕn term in the horizontally polarized array manifold can be introduced, so it will not be considered. Then, the above equation is simplified as
(7)y¯mt=anmvsnvt,anmhsnhtT=sntanmvanmhT⊙cn1cn2T≜sntanm⊙cn,where the vector cn=cn1,cn2T=cosγn,sinγnejηnT is the nth source polarization signature.For a dual-polarized sensor, we define the general short Fourier transform (GSFT) as(8)Dxit,f=∫−∞+∞wt−τ,σt,fxite−j2πfτdτ,where the Gaussian window is defined as
(9)wt,σt,f=12πσt,fexp−t22σ2t,f.Formulas (1)–(9) correspond to the single dual-polarization sensor case. With an M-sensor array, for each polarization i, i=vorh, the data vector is
(10)xit=x1it,x2it,⋯,xMitT=yit+nit=AiΦsit+nit.Using formula (10), the case of extending the polarization time-frequency distribution of a single sensor to multiple sensors is obtained. Instead of the scalar TFD of (8), we define the polarization STFD matrix of the data vector
(11)Dxit,f=∫−∞+∞wt−τ,σt,fxite−j2πfτdτ,which, for the no-noise case, can be represented as
(12)Dxit,f=AiΦDsit,fAiΦH.According to (10), the double polarization data vectors can be built as follows:
(13)xt=xvtxht=AvΦ00AhΦsvtsht+nvtnht=AvΦ00AhΦQvQhst+nvtnht=BΦQst+nt,where
(14)BΦ=AvΦ00AhΦ,represents block diagonal, and
(15)Q=QvQh,represents the signal polarization characteristic vector, where
(16)qv=cosγ1,⋯,cosγNT,Qv=diagqv,(17)qh=sinγ1ejη1,⋯,sinγNejηNT,Qh=diagqh.Correspondingly,(18)BΦQ=avϕ1cosγ1⋯avϕncosγNahϕ1sinγ1ejη1⋯ahϕnsinγNejηN=a~ϕ1⋯a~ϕN.The matrix (18) can be considered as the extended mixing matrix, within a~ϕn denoting the nth signal joint spatial polarization characteristic. For the nth signal, the extended spatial polarization characteristic vector is
(19)a~ϕn=avϕncosγnahϕnsinγnejηn.Obviously, the double polarization array can double the spatial dimension of the vector.The polarization, spatiality, and time-frequency characteristics can be combined with the source signal incident on the receiver array. The STFD of a dual-polarization data vectorxt can be as follows:
(20)Dxt,f=∫−∞+∞wt−τ,σt,fxte−j2πfτdτ.Dxt,f is called the SPTFD matrix. This matrix is a general method that can solve some common problems in array processing, as described in the following.When the noise effect is neglectful, the SPTFD matrix is connected with the signal TFD matrix by(21)Dxxt,f=BΦQDst,fQHBHΦ.Firstly, the time-frequency points are selected by the algorithm in the next section to construct the matrix. Then, eigenvectors are obtained by eigendecomposition of formula (21). Finally, the noise subspace Un constructed by eigenvectors is used for subsequent spatial spectrum estimation.


## 2.2. Three-Step Method for IF Estimation

Abdoush et al. [21] gave a general expression relating the error and variance of the IF estimate to the smoothness of the observation window and the IF. The correct IF estimation demands a data-based approach to choose at least reasonable window width. Therefore, the paper summarizes and extends the “low complexity” described in [19, 20]; the method includes three processes, as shown in Figure 2, which will be explained in this section. On this basis, this paper proposes a new linear TFR method based on estimation error and variance. The one uses a time-varying window for minimizing the MSE sum value estimated by IF at each moment, while in the latter TFR, the TF value is adjustable to minimize the MSE of each point in the TF region.Figure 2
Flowchart of adaptive TFRs based on linear TF transforms.For practical application, we will limit ourselves to a discrete version of GSFT, which we define as(22)Dxkt,kf=∑n=0N−1xlwkt−l,σkt,kfe−j2πkfl/N,where xkt represents a length Kt discrete time series corresponding to xktTs by sampling interval Ts and kf=−Kt/2,⋯,Kt/2−1 denoting frequency index.


### 2.2.1. Preliminary TFR

For any TFR, the optimized formula (22) requires an understanding of the IF. So, in adaptation, the first step is to create a preliminary TFR (PTFR) that roughly estimates the IF orbit. Although this stage can tolerate minor estimation errors, however, more significant estimation errors are propagated to the after step, so the final IF estimation accuracy decreases. To this end, Pei and Huang in [20] adopted an STFT of optimizing window width according to the ECM. While it performed well, when we used different multicomponent signals to detect this ECM, we found that the approach was beneficial for a noisy element and detrimental for a rapidly changing part. That is, in TFR, only quasistable components are selected. In this paper, we design an STFT based on the standard deviation of conventional rates based on PTFR. The rate is defined as the ratio of the effectivity bandwidth Beff to the effectivity duration Teff of the signal being analyzed.
(23)ν0=BeffTeff=1N∑k=−Kt/2Kt/2−1k−k02Xk2∑kt=0Kt−1kt−kt02xkt21/2,where Xk represents the discrete Fourier transform of xkt:
(24)k0=∑k=−Kt/2Kt/2−1kXk2∑k=−Kt/2N/2−1Xk2,kt0=∑kt=0Kt−1ktxkt2∑kt=0Kt−1xkt2.If the normalized power spectral density (that is, the normalized unit) of a signal is regarded as a function of probability mass, the second central moment of the signalBeff2 is the second central moment of its frequency, Beff is its standard deviation, and it represents the bandwidth of the signal. We can think of Teff as the time width of a signal, as before. In this way, the conventional ratio ν0 can be understood as the linear chirp scanning the frequency range Beff throughout Teff. We use the regular rate to determine the overall constant change rate of the signal IF regularity, which can be used as a fixed observation window for the PTFR to set the standard deviation of the STFT according to σ=3/71/4/2πν0.


### 2.2.2. IF Rate Estimation

In stage 2 of the adaptive algorithm, we get a rough estimate of IFR (IFR: IF rate, that is first derivative of IF). For this end, the IF is estimated by the main trajectory of PTFR. On this basis, the IFR is obtained by the IF first derivative. If the signal is multicomponent, each component must be found before differentiating each IF ridge. Here, the IF component separation is implemented according to the image processing algorithm described in [24]. It involves converting the TFR to a binary image and then performing a component connection process, in the case of TFR longer than a predetermined threshold, according to the adjacent connection criterion, extracting the connection peak from the TFR. This approach does not require prior knowledge of the relevant components but instead assumes that the IF ridge can be distinguishable in the TFR, which are many real-world signals, including most signals [25–27].For the derivation, [19] adopted the conventional difference operator. However, we know that this method is very sensitive to identifying discrete signals with noise. Therefore, the signal-to-noise ratio of the differential signal is usually lower than that of the original signal. Thus, in [20], the difference operator is replaced by principal component analysis, which executes eigenvalue decomposition for the 2 L estimations neighbouring f^Ikt to approximate fI′kt. In this article, we use a more straightforward way to find first-order differentials. First, the estimates f^Ikt are smoothed according to
(25)f^~Ikt=12L+1∑kt−Lkt+Lf^Ikt,kt=L,⋯,Kt−L,where L represents a positive integer. The above smoother represents a simple average moving filter used to improve the SNR of the differential signal [28]. After above smoothing the signal, the IFR is estimated using a traditional difference operator. As L increases, the degree of noise reduction increases, but the attenuation of the differential amplitude also becomes more remarkable. In addition, filtering the IF estimates according to (25) can help exclude possible outliers before derivation. Still, the accuracy of this method is limited because there is a certain relation between consecutive sampling errors in TFR [10]. Therefore, be careful when choosing L. Here, the smoothing width is 11.


### 2.2.3. Optimizing Window Width

After extracting the IF and IFR, they can be optimized for appropriate window widths. At this step, two ways adjust the TFR window width. In one case, the window width can be set in terms of time or frequency. For all signals, the primary purpose of this width optimization is to perform an optimal tradeoff at each moment or frequency. Therefore, when the window width is optimal, its cost function is the smallest. Since both adaptivity STFT (ASTFT) and ST can be done rapidly with FFT, the advantage of this method is that its complexity is slight. On the other hand, you can set a window to change at different times and frequencies. Within the TF region, width adjustment is allowed, and the second method requires more flexibility; however, sinceσkt,kf is related to kt and kf, in (22), GSFT cannot be implemented with FFT.We design an ASTFT that minimizes the IF estimation MSE sum at any moment. Because the TFRs generated by ASTFTs generally do not have crossterms, the TFRs can be approximated to the TFRs corresponding to the individual components of the overall signal. In this example, the formula can be generalized to the case of multicomponent signals with TF regions. That is, with ASTFT, an estimate of the IF with no offset can be obtained, and the total MSE can be obtained:(26)MSEtotkt,σkt=∑l=1Ln2πfIl′ktσ2kt2+15/232π5/2σ3ktC0,where Ln represents the signal component number at the moment instant ktTs and fIl′kt represents the lth single IF component IFR. The optimal window width is defined by minimizing MSEtotkt,σkt, which may be approximated by
(27)σoptkt≈argminσkt∑l=1Ln2πf∧Il′ktσ2kt2+15/2σ3kt.Note thatMSEtotkt,σkt represents the summation of all functions, each of which is completely convex and has an overall minimum value on σlkt, and there also is an overall minimum value for MSEtotkt,σkt. Therefore, it is contained in the following range:
(28)σoptkt∈minσ1kt,⋯,σLnkt,maxσ1kt,⋯,σLnkt.A good approximation ofσoptkt is obtained in two steps: first, σlkt is found corresponding to the estimated IFR minimum and maximum; thus, for width between the above two limits, the optimal candidates depend on the objective function form (27). The TFR we obtained is called optimal ASTFT (OASTFT).To reduce the sum of MSE estimated by IF, the width of the frequency domain correlation window in ST is optimized. It is required to obtain the 1st- and 2nd-order frequency derivative information of the window widths. Not knowing the nature of the differential equation, it has to be designed accordingly, making the method impractical; therefore, we cannot go on like this.The previous OASTFT algorithm has achieved a good balance among the signal components in the TF region, but it needs to be adjusted more to make it more adaptable. If the IF ridges can be separated within the TF region, then the window width of the GSFT can be adjusted individually. We can roughly estimate the IFs and IFRs, so the problem is how to defend the standard deviation of a window about the on-ridge and off-ridge time-frequency points. At IF, the first and second derivatives of the standard deviation and frequency are 0, and then, window width is optimal. Therefore, we define a window width intervalΔf determined by the following formula on each upper ridge point of the IF component according to
(29)σkt,kf=3/71/42πf^Il′kt,0≤kf−f^Ilkt≤Δf2,whereas outside the section, the window width is achieved by two-dimensional linear interpolation in [20]. The last interpolation ensures that the σkt,kf 1st- and 2nd-order frequency derivatives are close to 0 neighbouring f^Ilkt, so the width chosen is best for MSE. Also, the coarse IF estimate can be far from its actual value due to noise. Therefore, preferably, not only the width of (29) is specified for f^Ilkt but also for adjacent frequencies. In experiments, we found that several frequency bands are wide enough to obtain accurate IF estimates. The resulting transformation is expressed by optimal GSFT (OGSFT).The adaptation plan adopted consists of the three phases mentioned above. Its main work is to realize PTFR, extract and track IF components, thus obtain the corresponding IFR value, and realize adaptive TFR. The PTFR designed in this paper is STFT, and its computational complexity isON2logN. The chosen IF component extraction and tracking algorithm is presented in [24]; the computational cost for this analysis procedure is outside the scope of this work. The final TFR can be a time-adaptive ASTFT or a TF-adaptive GSFT. OASTFT can be done by FFT. That is, OASTFT needs to operate on the order of ON2logN. The TF-adaptive GSFT algorithm has high adaptability, but because the FFT program cannot complete the conversion, it will lead to ON3 level of computational complexity.


## 2.2.1. Preliminary TFR

For any TFR, the optimized formula (22) requires an understanding of the IF. So, in adaptation, the first step is to create a preliminary TFR (PTFR) that roughly estimates the IF orbit. Although this stage can tolerate minor estimation errors, however, more significant estimation errors are propagated to the after step, so the final IF estimation accuracy decreases. To this end, Pei and Huang in [20] adopted an STFT of optimizing window width according to the ECM. While it performed well, when we used different multicomponent signals to detect this ECM, we found that the approach was beneficial for a noisy element and detrimental for a rapidly changing part. That is, in TFR, only quasistable components are selected. In this paper, we design an STFT based on the standard deviation of conventional rates based on PTFR. The rate is defined as the ratio of the effectivity bandwidth Beff to the effectivity duration Teff of the signal being analyzed.
(23)ν0=BeffTeff=1N∑k=−Kt/2Kt/2−1k−k02Xk2∑kt=0Kt−1kt−kt02xkt21/2,where Xk represents the discrete Fourier transform of xkt:
(24)k0=∑k=−Kt/2Kt/2−1kXk2∑k=−Kt/2N/2−1Xk2,kt0=∑kt=0Kt−1ktxkt2∑kt=0Kt−1xkt2.If the normalized power spectral density (that is, the normalized unit) of a signal is regarded as a function of probability mass, the second central moment of the signalBeff2 is the second central moment of its frequency, Beff is its standard deviation, and it represents the bandwidth of the signal. We can think of Teff as the time width of a signal, as before. In this way, the conventional ratio ν0 can be understood as the linear chirp scanning the frequency range Beff throughout Teff. We use the regular rate to determine the overall constant change rate of the signal IF regularity, which can be used as a fixed observation window for the PTFR to set the standard deviation of the STFT according to σ=3/71/4/2πν0.


## 2.2.2. IF Rate Estimation

In stage 2 of the adaptive algorithm, we get a rough estimate of IFR (IFR: IF rate, that is first derivative of IF). For this end, the IF is estimated by the main trajectory of PTFR. On this basis, the IFR is obtained by the IF first derivative. If the signal is multicomponent, each component must be found before differentiating each IF ridge. Here, the IF component separation is implemented according to the image processing algorithm described in [24]. It involves converting the TFR to a binary image and then performing a component connection process, in the case of TFR longer than a predetermined threshold, according to the adjacent connection criterion, extracting the connection peak from the TFR. This approach does not require prior knowledge of the relevant components but instead assumes that the IF ridge can be distinguishable in the TFR, which are many real-world signals, including most signals [25–27].For the derivation, [19] adopted the conventional difference operator. However, we know that this method is very sensitive to identifying discrete signals with noise. Therefore, the signal-to-noise ratio of the differential signal is usually lower than that of the original signal. Thus, in [20], the difference operator is replaced by principal component analysis, which executes eigenvalue decomposition for the 2 L estimations neighbouring f^Ikt to approximate fI′kt. In this article, we use a more straightforward way to find first-order differentials. First, the estimates f^Ikt are smoothed according to
(25)f^~Ikt=12L+1∑kt−Lkt+Lf^Ikt,kt=L,⋯,Kt−L,where L represents a positive integer. The above smoother represents a simple average moving filter used to improve the SNR of the differential signal [28]. After above smoothing the signal, the IFR is estimated using a traditional difference operator. As L increases, the degree of noise reduction increases, but the attenuation of the differential amplitude also becomes more remarkable. In addition, filtering the IF estimates according to (25) can help exclude possible outliers before derivation. Still, the accuracy of this method is limited because there is a certain relation between consecutive sampling errors in TFR [10]. Therefore, be careful when choosing L. Here, the smoothing width is 11.


## 2.2.3. Optimizing Window Width

After extracting the IF and IFR, they can be optimized for appropriate window widths. At this step, two ways adjust the TFR window width. In one case, the window width can be set in terms of time or frequency. For all signals, the primary purpose of this width optimization is to perform an optimal tradeoff at each moment or frequency. Therefore, when the window width is optimal, its cost function is the smallest. Since both adaptivity STFT (ASTFT) and ST can be done rapidly with FFT, the advantage of this method is that its complexity is slight. On the other hand, you can set a window to change at different times and frequencies. Within the TF region, width adjustment is allowed, and the second method requires more flexibility; however, sinceσkt,kf is related to kt and kf, in (22), GSFT cannot be implemented with FFT.We design an ASTFT that minimizes the IF estimation MSE sum at any moment. Because the TFRs generated by ASTFTs generally do not have crossterms, the TFRs can be approximated to the TFRs corresponding to the individual components of the overall signal. In this example, the formula can be generalized to the case of multicomponent signals with TF regions. That is, with ASTFT, an estimate of the IF with no offset can be obtained, and the total MSE can be obtained:(26)MSEtotkt,σkt=∑l=1Ln2πfIl′ktσ2kt2+15/232π5/2σ3ktC0,where Ln represents the signal component number at the moment instant ktTs and fIl′kt represents the lth single IF component IFR. The optimal window width is defined by minimizing MSEtotkt,σkt, which may be approximated by
(27)σoptkt≈argminσkt∑l=1Ln2πf∧Il′ktσ2kt2+15/2σ3kt.Note thatMSEtotkt,σkt represents the summation of all functions, each of which is completely convex and has an overall minimum value on σlkt, and there also is an overall minimum value for MSEtotkt,σkt. Therefore, it is contained in the following range:
(28)σoptkt∈minσ1kt,⋯,σLnkt,maxσ1kt,⋯,σLnkt.A good approximation ofσoptkt is obtained in two steps: first, σlkt is found corresponding to the estimated IFR minimum and maximum; thus, for width between the above two limits, the optimal candidates depend on the objective function form (27). The TFR we obtained is called optimal ASTFT (OASTFT).To reduce the sum of MSE estimated by IF, the width of the frequency domain correlation window in ST is optimized. It is required to obtain the 1st- and 2nd-order frequency derivative information of the window widths. Not knowing the nature of the differential equation, it has to be designed accordingly, making the method impractical; therefore, we cannot go on like this.The previous OASTFT algorithm has achieved a good balance among the signal components in the TF region, but it needs to be adjusted more to make it more adaptable. If the IF ridges can be separated within the TF region, then the window width of the GSFT can be adjusted individually. We can roughly estimate the IFs and IFRs, so the problem is how to defend the standard deviation of a window about the on-ridge and off-ridge time-frequency points. At IF, the first and second derivatives of the standard deviation and frequency are 0, and then, window width is optimal. Therefore, we define a window width intervalΔf determined by the following formula on each upper ridge point of the IF component according to
(29)σkt,kf=3/71/42πf^Il′kt,0≤kf−f^Ilkt≤Δf2,whereas outside the section, the window width is achieved by two-dimensional linear interpolation in [20]. The last interpolation ensures that the σkt,kf 1st- and 2nd-order frequency derivatives are close to 0 neighbouring f^Ilkt, so the width chosen is best for MSE. Also, the coarse IF estimate can be far from its actual value due to noise. Therefore, preferably, not only the width of (29) is specified for f^Ilkt but also for adjacent frequencies. In experiments, we found that several frequency bands are wide enough to obtain accurate IF estimates. The resulting transformation is expressed by optimal GSFT (OGSFT).The adaptation plan adopted consists of the three phases mentioned above. Its main work is to realize PTFR, extract and track IF components, thus obtain the corresponding IFR value, and realize adaptive TFR. The PTFR designed in this paper is STFT, and its computational complexity isON2logN. The chosen IF component extraction and tracking algorithm is presented in [24]; the computational cost for this analysis procedure is outside the scope of this work. The final TFR can be a time-adaptive ASTFT or a TF-adaptive GSFT. OASTFT can be done by FFT. That is, OASTFT needs to operate on the order of ON2logN. The TF-adaptive GSFT algorithm has high adaptability, but because the FFT program cannot complete the conversion, it will lead to ON3 level of computational complexity.


## 2.3. Polarimetric Time-Frequency Music and Spatial Polarimetric Correlations

For improving the signal spatial resolution with well-defined time-frequency characteristics, time-frequency multiple signal classification (TF-MUSIC) has recently been proposed [29]. This method provides an essential summary for applying MUSIC to polarization arrays, which searches for the minimum of the array vector derived from the SPTFD matrix defined in the combined space and polarization domains. The advantages of the algorithm mainly include the following two points: the algorithm integrates polarization information and time-frequency information to increase the utilization rate of multidimensional information. The algorithm adaptively recognizes the time-frequency points of autoterms and avoids the trouble of manual selection.Consider the following spatial signature matrix:(30)Fϕ=1Mavϕ00ahϕ.Becauseaiϕ2=M, FHϕFϕ represents the identity matrix. For searching from the spatial and polarization joint domains, the spatial polarization vector is defined as
(31)fϕ,c=FϕcFϕc=Fϕc,where the vector c=c1c2T represents an unknown unit norm vector of polarization coefficients. In (31), Fϕc=cHFHϕFϕc1/2=cHc1/2=1.The PTF-MUSIC spectrum is provided from the following function:(32)Pϕ=mincfHϕ,cUnUnHfϕ,c−1=minccHFHϕUnUnHFϕc−1,where Un represents the noise subspace from the SPTFD matrix of (20) utilized the selected time-frequency points. For DOA estimation, different STFD or SPTFD matrices can be combined using TF average and joint block diagonalization techniques [29]. Picking these points from regions of high energy concentration associated with global or local sources can improve the signal-to-noise ratio and make the MUSIC algorithm more noise-robust than its traditional MUSIC counterpart [30]. In (32), by acquiring the smallest eigenvalue for the matrix FHϕUnUnHFϕ, this method can perform simple eigendecomposition on the 2×2 matrix, thus avoiding performing a large number of operations in the polarization region. In this way, the spectrum of PTF-MUSIC can be described as
(33)Pϕ=λmin−1FHϕUnUnHFϕ,where λmin⋅ represents the smallest eigenvalue solution. The DOA value of the source is the highest in the spectrum. Corresponding to each angle ϕn of arrival of N signals, n=1,2,⋯,N, each polarized signal parameter can be estimated as follows:
(34)c^ϕk=vminFHϕkUnUnHFϕk,where vmin⋅ represents the eigenvector corresponding to the minimality eigenvalue λmin⋅.The spatial resolution capability of the array has a lot to do with the propagation characteristics of the signal [30]. This depends on the manifold vector normalized inner product of each array. For potential problems involving spatial and polarization dimensions, the extended array manifold a~ϕ is used to define the joint spatial polarization correlation coefficient between the source nl and nk, that is,
(35)βnl,nk=1Ma~Hϕnka~ϕnl=cnk1∗cnl1βnl,nkv+cnk2∗cnl2βnl,nkh=1Mcnk1∗cnl1avϕnkHavϕnl+cnk2∗cnl2ahϕnkHahϕnl,where βnl,nki=1/MaiϕnkHaiϕnl represents the spatial relation coefficient.An interesting case phenomenon occurs with the identical array manifolds between vertically and horizontally polarization, i.e.,avϕ=ahϕ. For the condition, βnl,nkv=βnl,nkh, and the spatial polarization joint relation coefficient is the individual spatial polarimetric correlations product:
(36)βnl,nk=βnl,nkvρnl,nk,within
(37)ρnl,nk=cnkHcnl=cosγnlcosγnkejηnl−ηnk+sinγnlsinγnk,denoting the polarization correlation coefficient. Particularly, for linear polarization, ηnl=ηnk=0, and (37) reduces to
(38)ρnl,nk=cosγnl−γnk.Sinceρnl,nk≤1, the equality holds only when the polarization states of the two light sources are the same; compared with the single spatial correlation coefficient, the spatial polarization correlation coefficient is small. Due to the polarization diversity dropping correlation value, the introduced method ρnl,nk was adopted, which was converted into enhanced signal source identification. In this way, it is difficult to resolve the two sources with a single monopolar spatial array manifold avϕ or ahϕ with an extended spatial polarized array manifold, defined by a~ϕ, which can be easily separated. This improvement is more pronounced when the signal spatial relation coefficient is large, but the individual polarization correlations are low.


## 3. Results and Discussion

This paper focuses on the improvement of the adaptive time-frequency analysis algorithm, so the polarization parameters are assumed to be known prior parameters in the simulation process. The results of the three-step adaptive time-frequency analysis are analyzed by simulation. Simulation conditions are as follows: the signal component 1 normalized frequency is from 0.15 to 0.1; the signal component 2 normalized frequency is from 0.2 to 0.4; the number of sampling points is 1024; the SNR is 0 dB. The simulation results are shown in Figure3. The black asterisk line in the figure represents the truth value; the red box line represents the estimated result. It can be seen that the algorithm can accurately estimate the instantaneous frequency of the signal. It is possible to select time-frequency points to construct the spatially polarized time-frequency distribution matrix. Then, the spatial spectrum estimation algorithm is simulated and analyzed under the following conditions: the normalized frequency of the two signal components is the same as before; the incidence angles are 20 degrees and 40 degrees, respectively. The number of snapshots is 512; the SNR is 0 dB. Simulation results of spatial spectrum estimation of polarization time-frequency MUSIC algorithm based on adaptive time-frequency analysis algorithm are shown in Figure 4. The blue triangle curve is the spatial spectrum estimation result of signal component 1. The green circle curve is the result of spatial spectrum estimation of signal component 2. It can be seen from the marked DOA estimation results of the figure that the algorithm can accurately distinguish and estimate the angle of the incident signal.Figure 3
The results of time-frequency analysis.Figure 4
The spatial spectrum estimation.Finally, the performance of the proposed algorithm is simulated and analyzed. The relationship between RMSE and SNR and snapshot number is analyzed, and the time-frequency MUSIC algorithm based on traditional time-frequency analysis is compared. The simulation analysis conditions of RMSE and SNR are as follows: the normalized frequency and incident angle of the two signals are the same as before; the number of snapshots is fixed at 512. The SNR increased from -11 dB to 15 dB with a step size of 1 dB. The simulation results are shown in Figure5. In the figure, the black star curve is the estimation result of the proposed algorithm, and the red circle curve is the result of the traditional algorithm. The results show that the RMSE of both methods decreases in the case of SNR increasing, and the proposed algorithm is significantly better than the traditional algorithm. The effect is more evident at low SNR. At about 0 dB, there is still about 0.2 degree of improvement. The simulation analysis conditions of RMSE and snapshot number are as follows: the normalized frequency and incident angle of the two signals are the same as before; the SNR is fixed at 0 dB; the number of snapshots increases gradually from 64 to 1024, taking the whole power of 2. The simulation results are shown in Figure 6. In the figure, the black star curve is the estimation result of the proposed algorithm, and the red circle curve is the result of the traditional algorithm. The results show that the RMSE of both methods decreases in the case of snapshot number increasing, and the proposed algorithm is obviously better than the traditional algorithm. When the number of snapshots is low, due to the small number of samples, the results of the two algorithms are relatively close, but the proposed algorithm is still better than the traditional algorithm. As the number of snapshots increases, the RMSE of the proposed algorithm decreases rapidly compared with the traditional algorithm. In practice, considering the actual finite length received data, the spatially polarized time-frequency distribution matrix is estimated according to the received data obtained by many snapshot number. As the snapshot number increases, the matrix estimation is more accurate. After the eigendecomposition, the noise subspace estimation is more accurate, thus improving DOA estimation performance.Figure 5
The curve between RMSE and SNR.Figure 6
The curve between RMSE and snapshot number.Additionally, the relationships between the estimation success probability and SNR and snapshot number are analyzed, and the time-frequency MUSIC algorithm based on traditional time-frequency analysis is compared. The simulation analysis conditions of the estimation success probability and SNR are as follows: the normalized frequency and incident angle of the two signals are the same as before; the number of snapshots is fixed at 512. The SNR increased from -11 dB to 5 dB with a step size of 1 dB. When the estimation error is less than 1.5 degrees, the estimation is judged to be successful. The simulation results are shown in Figure7. In the figure, the black star curve is the estimation result of the proposed algorithm, and the red circle curve is the result of the traditional algorithm. The results show that the estimation success probability of both methods increases with SNR increasing, and the proposed algorithm is significantly better than the traditional algorithm. The improvement of the estimation success probability is more evident at low SNR. The simulation analysis conditions of the estimation success probability and the snapshot number are as follows: the normalized frequency and incident angle of the two signals are the same as before; the SNR is fixed at -2 dB; the number of snapshots increases gradually from 64 to 1024, taking the whole power of 2. When the estimation error is less than 1.5 degrees, the estimation is judged to be successful. The simulation results are shown in Figure 8. In the figure, the black star curve is the estimation result of the proposed algorithm, and the red circle curve is the result of the traditional algorithm. The results show that the estimation success probability of both methods increases with snapshot number increasing, and the proposed algorithm is obviously better than the traditional algorithm. When the number of snapshots is low, due to the small number of samples, the results of the two algorithms are general, but the proposed algorithm is still better than the traditional algorithm. As the number of snapshots increases, the estimation success probability of the proposed algorithm increases rapidly compared with the traditional algorithm. From many aspects, the performance of the proposed algorithm is obviously due to the traditional algorithm. The polarization spatial time-frequency DOA estimation algorithm combines the three-step adaptive time-frequency analysis algorithm to achieve good results.Figure 7
The curve between the estimation success probability and SNR.Figure 8
The curve between the estimation success probability and snapshot number.


## 4. Conclusions

This paper introduces a new method for estimating the direction of arrival of polarized waves using time-frequency adaptive linear time-frequency transforms. A linear TFR was used as an IF estimate, primarily due to its simplicity and immunity to cross-interference. The optimal window relies on the IF 1st-order derivative. Therefore, a simple algorithm is used for deriving the derivative and optimizing it accordingly. Due to its combination with signal polarization, STFD gains more freedom and thus promotes the estimation of the signal and noise phonon space. On the SPTFD platform, PTF-MUSIC is used for the estimation of signal DOA, which outperforms conventional time-frequency MUSIC. Through the analysis of the simulation results in the previous section, it is proved that the proposed algorithm is superior to the traditional algorithm in many aspects. Compared with the traditional algorithm, the improvement is more significant.



---

*Source: 3463438-2022-08-08.xml*