Nonlinear Blind Compensation for Array Signal Processing Application

Abstract

Recently, nonlinear blind compensation technique has attracted growing attention in array signal processing application. However, due to the nonlinear distortion stemming from array receiver which consists of multi-channel radio frequency (RF) front-ends, it is too difficult to estimate the parameters of array signal accurately. A novel nonlinear blind compensation algorithm aims at the nonlinearity mitigation of array receiver and its spurious-free dynamic range (SFDR) improvement, which will be more precise to estimate the parameters of target signals such as their two-dimensional directions of arrival (2-D DOAs). Herein, the suggested method is designed as follows: the nonlinear model parameters of any channel of RF front-end are extracted to synchronously compensate the nonlinear distortion of the entire receiver. Furthermore, a verification experiment on the array signal from a uniform circular array (UCA) is adopted to testify the validity of our approach. The real-world experimental results show that the SFDR of the receiver is enhanced, leading to a significant improvement of the 2-D DOAs estimation performance for weak target signals. And these results demonstrate that our nonlinear blind compensation algorithm is effective to estimate the parameters of weak array signal in concomitance with strong jammers.

1. Introduction

During the past few decades, array signal processing (ASP) has been a world-wide hot spot of research [1,2], with discoveries as diverse in application as wireless communications, multiple-input multiple-output (MIMO) radar and sonar [3,4,5,6]. ASP is concerned with the problem of extracting high dimensional information of interest which is received from an array of spatially distributed sensors and transmitted to a multi-channel array receiver. There are normally three units in an array system: antenna arrays, array receiver and array signal processing unit [7]. Moreover, an array receiver, comprising of multi-channel radio frequency (RF) front-ends, contains many nonlinear circuits, bringing about the vast majority of nonlinear distortion components in an array signal processing system.

Currently, much research regarding array signal processing focuses on the development of parameter estimation methods [8,9,10,11]. Amongst them, two-dimensional direction-of-arrival (2-D DOA) estimation is the hottest topic and most of the published literature focuses on parametric algorithms to deal with it [9,10,11]. These methods are proven to perform well in source direction finding under ideal settings, which ignore the negative effect of the nonlinear distortion components of an array signal processing system on source localization accuracy. However, due to the nonlinearity of an array receiver, the signal-to-interference-plus-noise ratio (SINR) of several target signals is reduced. Such imperfect methods lead to sharply deteriorated performance under a condition of low SINR. Therefore, further research is required to better address the parameters estimation problem of an array signal in such a scenario as an electromagnetic environment in which an array signal processing system works.

Some effective methods [12,13] for such a scenario adopt high-order statistics to improve the direction-finding performance of array signal which, however, increase the computational complexity enormously compared with the methods such as the multiple signal classification (MUSIC) algorithm applying second-order statistics. To resolve the difficulty without increasing the computation burdens of the direction finding, the authors enhance the spurious-free dynamic range (SFDR) of an array receiver, which leads to significant SINR improvements of weak target signals. Therefore, to promote the development of array signal processing, an array receiver with high SFDR is a vital undertaking. Fewer researches for blindly mitigating the nonlinear distortion of multi-channel array receiving systems have been conducted, to the best of the authors’ knowledge.

The nonlinearity mitigation technique for a one-channel receiving system has been reported in some papers [14,15,16,17,18,19,20,21,22,23]. Several preliminary ideas of Reference [14] offer some insight for this research. This paper proposes a nonlinear blind compensation method to eliminate the distortion product of band-limited signal by suppressing the out-of-band distorted signal. However, this method is only applicable to compensating memoryless and monotonic nonlinearity under the premise of the known input signal bandwidth. A digital post-compensation technique is adopted to mitigate the nonlinear distortion of the multicarrier receiving system [15]. Moreover, for linearizing a direct-conversion receiver (DCR), a full-band adaptive interference cancellation (FB-AIC) method has been proposed [16]. This method utilizes the out-band distortions to identify and extract the nonlinear model coefficients. Furthermore, in Reference [17], some effective works have been done to degrade the third-order intermodulation distortion of DCR severely. Recently, the paper [18] developed a blind nonlinear post-distortion algorithm with follow-up blind finite impulse response (FIR) equalizers to compensate for the distortion product in the concurrent dual-band receiver, the memory effects and the non-ideal wireless channel response. Although References [15,16,17,18,19] mainly focus on cancelling the third-order inter-modulation distortion, they fail to be applied to a wideband nonlinear system. The current authors’ initial research [20] proposed that minimizing the short-time energy of the nonlinearity in the frequency domain is regarded as the blind identification criterion. However, the identification strategy might have difficultly to coping with the situation where multiple unknown signals with different power levels exist simultaneously, due to the non-ideal frequency response characteristic of FIR. Moreover, in References [18,20,21,22,23], the nonlinear behavior model parameters are identified and extracted in the frequency domain, resulting in heavy computational costs.

To solve these identification problems, separating large signals (strong array signals) and small signals (nonlinear distortion components and weak array signals) in time domain according to the power level is appropriate. Large signals produce most of the nonlinear components and are usually of far higher power level compared with small signals. Thereupon, the instantaneous power residual error between the small signals and nonlinear model of large signals can be considered as an objective function of nonlinear model kernels. Consequently, a blind identification criterion is designed to adaptively resize nonlinear model coefficients and an improved weighed iterative method (IWI) is presented to minimize the instantaneous power residual error based on the method in research [24].

In References [25,26,27,28,29], some blind signal separation algorithms have been investigated; for instance, the singular value decomposition (SVD)-based and the eigenvalue decomposition (EVD)-based methods. The SVD-based algorithm is the standard due to its good performance and applicability, as confirmed in References [26,28,29]. As learned from the References [26,28,29], the performance for large signal extraction is mainly determined by an adequate choice of the number of the largest singular values, $r$. However, it is hard to precisely determine $r$ on account of the unknown array signal. A smaller $r$ will bring about a part of the large signals, while a bigger one will draw some unwanted nonlinearity or weak array signals into the large signals. Both the conditions will result in incorrect extracted outcomes. Furthermore, the engineering implementation of a SVD-based algorithm needs massive computational loads [29].

A novel nonlinear blind compensation algorithm is suggested to cope with array signal processing. Under the existing semiconductor process levels and device performance status, this study focuses on significantly improving the SFDR of an array receiver and the parameters estimation accuracy of the weak target signals, (for example, their 2-D DOAs). To the best of the authors’ knowledge, the MUSIC algorithm is one of the most widely-used high-resolution 2-D DOA estimation approaches. Therefore, the authors consider it to be an efficient azimuth-elevation angles estimation method for an array signal.

Additionally, the authors propose a novel blind signal separation algorithm, a spectrum reduction algorithm based on time-frequency conversion (SRA-TFC) algorithms wherein large signals and small signals are extracted based on the power level in the frequency domain. To follow, their time domain waveforms are obtained by the Inverse Discrete Fourier transform (IDFT) technique. Comparing with the traditional SVD-based method, the proposed method is only necessary to set the power threshold rather than the number of the largest singular values of constant change. Therefore, the method is highly suitable to blindly separate signals in the electromagnetic environment.

Moreover, the identification criterion should be applied to each channel for extracting its identification coefficients if each channel is regarded to be an independent and irrelevant nonlinear system. The disadvantages are as follows: first, the system hardware resource is increased significantly. Most importantly, the inconsistent compensation performances of different channels, especially the phase disturbance of the post-compensation array signal, are caused by the iterative computation error of the proposed identification criterion. However, the input array signal and the linear/nonlinear transfer functions of each channel are similar. A more suitable method is to blindly identify and extract the nonlinear model kernel coefficients of any one of the channels and then to employ the extracted kernel coefficients to synchronously mitigate the nonlinearity of all channels. Consequently, the authors propose an effective digital compensation method for an array signal processing system.

The arrangement of this paper is as follows: Section 2 describes the nonlinear blind compensation processing architecture for an array system and the theoretical demonstration of the proposed algorithm. It is followed up in Section 3 with the nonlinear blind mitigation architecture. The experimental results are given in Section 4 to demonstrate the performance improvement of the proposed algorithm. Finally, a short conclusion is drawn in Section 5.

2. System Model and Problem Analysis

2.1. ASP Architecture Based on Nonlinear Blind Compensation

Figure 1 shows a structural diagram of the nonlinear blind compensation processing for an array system, which is composed of an antenna arrays unit, array receiver and blind compensation unit. Where, $x_i(t)$, $y_i(n)$ are the analog input, digital output in the $i$th channel RF front-end, respectively; $z_i(n)$ is the $i$th channel post-compensation signal; $f_s$ means the sampling rate of the synchronous clock; $w_c$, $B_t$ stands for the center frequency and bandwidth of the target signal.

Additionally, in the sub-figure across the top of Figure 1, it contains three power spectrum schematic diagrams of different signals which are input signals, output signals and completely compensated signals of an array receiver, separately. All red regions in these diagrams mean the power spectrum of a large array signal, while all purple regions are that of weak target signals. Moreover, the blue region shows the power spectrum of the nonlinear distortion stemming from the array signal through the array receiver.

The antenna arrays unit exists at the head of the system, where it is to capture the array signal. Then, the array data is transmitted to the array receiver. Each channel RF front-end of the array receiver has some nonlinear analog circuits and an ADC (analog to digital converter). The analog module usually consists of a bandpass filter (BPF) with high suppression ratio, a low-noise amplifier (LNA), a variable-gain amplifier (VGA), a low-pass filter (LPF) and more. The full-band array signal is achieved by filtering the input signal with a BPF in each channel of the array receiver. It is then amplified to the required level by a first-stage LNA and a first-stage VGA. Additionally, the amplified signal without out-band clutter obtained by an LPF is transmitted to the synchronization ADC module in which all ADCs’ sampling clocks are from the same crystal oscillator. These procedures are applied to synchronously gather and output the array signal. The blind compensation unit is designed to mitigate or even eliminate the distortion product supplied from the multi-channel array processing system with the aim to improve the parametric estimation accuracy of the array signal.

2.2. Theoretical Analysis of the Proposed Method

This section offers the feasibility analysis of the novel nonlinear blind compensation method which adaptively extracts the nonlinear model kernel coefficients of any one of the channels to mitigate the nonlinearity of all channels synchronously.

The bandwidth of the multi-channel nonlinear systems, tens of MHz, is much bigger than that of the target signals, several KHz. Thereby a same memoryless nonlinear model representing for every single channel of the nonlinear systems could validate the presented nonlinear compensation method [30]. Additionally, it is well known that the $p$th order orthogonal polynomials with completely uncorrelated terms contributes to only the $p$th order nonlinear distortion term. Thus, the orthogonal polynomials modeling for nonlinear systems can offer an intuitive means of compensation improvement analysis to every nonlinear term. The authors adopt the Chebyshev polynomials, one of the most frequently used orthogonal polynomials.

It is generally accepted that the second-order and third-order nonlinear terms are the key factors to exacerbate the SFDR performance of the nonlinear system [20]. Therefore, a Chebyshev nonlinear model with third-order non-linear terms could be modeled to approximate the multi-channel nonlinear systems. Therefore, the authors consider the $i$th channel RF front-end a discrete-time third–order Chebyshev orthogonal nonlinear system:

$$\begin{array}{ccc}\hfill y_i& =& f_i(x_i)\hfill \\ & =& a_i{T}_0(x_i)+b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)\hfill \end{array}$$

(1)

where $f_i(•)$ represents the nonlinear transformation function of the $i$th channel RF front-end; $T_p(•)$ denotes $p$th order Chebyshev node; $x_i$ and $y_i$ are the simplified forms of $x_i(t)$ and $y_i(n)$, respectively; $a_i$ means the coefficient of the direct-current (DC) component; $b_i$ represents the linear term parameter, $c_i$ and $d_i$ are the coefficients of the second-order and third-order nonlinear terms.

The first three Chebyshev orthonormal polynomials are shown as:

$$\begin{array}{l}{T}_0(x)=1\\ T_1(x)=x\\ T_2(x)=2x^2-1\\ T_3(x)=4x^3-3x\end{array}$$

(2)

substituting Equation (2) into Equation (1) gives:

$$\begin{array}{ccc}\hfill y_i& =& a_i+b_i{x}_i+c_i(2x_i^2-1)+d_i(4x_i^3-3x_i)\hfill \\ & =& (a_i-c_i)+(b_i-3d_i)x_i+2c_i{x}_i^2+4d_i{x}_i^3\hfill \end{array}$$

(3)

Equation (3) has its DC component assumed to be zero. Furthermore, the coefficient of linear term ($b_i$) is much bigger than that of nonlinear terms ($c_i$, $d_i$) in a weak nonlinear system. Thus, Equation (3) can be rewritten:

$$\begin{array}{ccc}\hfill y_i& =& (b_i-3d_i)x_i+2c_i{x}_i^2+4d_i{x}_i^3\hfill \\ & \approx & b_i{x}_i+2c_i{x}_i^2+4d_i{x}_i^3(a_i=c_i)\hfill \end{array}$$

(4)

The model parameters of any one of the channels are extracted, in this paper the first channel for an example. Via the blind compensation for the first channel nonlinear system, its inverse behavior model of the nonlinear system is obtained approximately [20]. Therefore, according to the knowledge of Reference [31], the inverse nonlinear behavior model can be defined as:

$$\begin{array}{ccc}\hfill z_i& =& g_1(y_i)\hfill \\ & =& {\rho }_0{T}_0(y_i)+{\rho }_1{T}_1(y_i)+{\rho }_2{T}_2(y_i)+{\rho }_3{T}_3(y_i)\hfill \end{array}$$

(5)

where, $z_i(n)$ is reduced to $z_i$; $g_1(•)$ is the inverse of $f_1(•)$ (the nonlinear transformation function of the first channel) and its coefficient vector is represented by $\mathit{g}={[{\rho }_0,{\rho }_1,{\rho }_2,{\rho }_3]}^{\mathrm{T}}$; $\mathrm{T}$ denotes transposition operation.

Corresponding to Equation (4), the vector of power series coefficients of the first channel $\mathbf{F}$ is obtained:

$$\mathbf{F}={[0,b_1,2c_1,4d_1]}^{\mathrm{T}}$$

(6)

The p-fold convolution vectors are calculated by:

$${\mathbf{F}}^{(p)}(p=0,1,2,3)$$

and the matrix $\mathbf{H}$ assembled by juxtaposing these column vectors as shown below:

$$\mathbf{H}=[{\mathbf{F}}^{(0)},{\mathbf{F}}^{(1)},{\mathbf{F}}^{(2)},{\mathbf{F}}^{(3)}]$$

(7)

Normally, for inversion, the dimensions of $\mathbf{H}$ is truncated to a 4 × 4 matrix, denoted as:

$${\mathbf{H}}_{4\times 4}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& b_1& 0& 0\\ 0& 2c_1& {(b_1)}^2& 0\\ 0& 4d_1& 4b_1{c}_1& {(b_1)}^3\end{array}\right]$$

(8)

Furthermore, Chebyshev series-power series conversion matrix $\mathbf{T}$ is truncated as the following:

$$\mathbf{T}=\left[\begin{array}{cccc}1& 0& -1& 0\\ 0& 1& 0& -3\\ 0& 0& 2& 0\\ 0& 0& 0& 4\end{array}\right]$$

(9)

$\mathit{g}$ is gained by choosing from the second column of matrix $\mathsf{\Gamma }=\mathbf{T}·({\mathbf{H}}^{-1})·({\mathbf{T}}^{-1})$ as described in Reference [32]:

$$\begin{array}{l}{\rho }_0=\frac{c_1}{b_1^3}\\ {\rho }_1=\frac{b_1^4-3b_1{d}_1+6c_1^2}{b_1^5}\\ {\rho }_2=-\frac{c_1}{b_1^3}\\ {\rho }_3=\frac{2c_1^2-b_1{d}_1}{b_1^5}\end{array}$$

(10)

Substituting Equation (10) into Equation (5) gives:

$$z_i=\frac{c_1}{b_1^3}{T}_0(y_i)+\frac{b_1^4-3b_1{d}_1+6c_1^2}{b_1^5}{T}_1(y_i)-\frac{c_1}{b_1^3}{T}_2(y_i)+\frac{2c_1^2-b_1{d}_1}{b_1^5}{T}_3(y_i)$$

(11)

According to the knowledge of Chebyshev polynomials:

$$\begin{array}{l}{T}_0(x)=1\\ T_1(x)=x\\ T_n(x)\times T_m(x)=\frac{T_{n+m}(x)+T_{n-m}(x)}{2}n\ge m\end{array}$$

(12)

Therefore, the result is:

$$\begin{array}{ccc}\hfill T_0(y_i)& =& 1\hfill \\ \hfill T_1(y_i)& =& y_i=c_i+b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)\hfill \\ \hfill T_2(y_i)& =& 2{(y_i)}^2-1\hfill \\ & =& 2{[(c_i+b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)]}^2-1\hfill \\ & =& 2[c_i^2+b_i^2{(T_1(x_i))}^2+c_i^2{(T_2(x_i))}^2+d_i^2{(T_3(x_i))}^2\hfill \\ & & +2b_i{c}_i{T}_1(x_i)+2c_i^2{T}_2(x_i)+2c_i{d}_i{T}_3(x_i)+2b_i{c}_i{T}_1(x_i)T_2(x_i)\hfill \\ & & +2b_i{d}_i{T}_1(x_i)T_3(x_i)+2c_i{d}_i{T}_2(x_i)T_3(x_i)]-1\hfill \\ \hfill T_3(y_i)& =& 4{(y_i)}^3-3y_i=y_i[2(2{(y_i)}^2-1)-1)=T_1(y_i)[2T_2(y_i)-1]\hfill \end{array}$$

(13)

This process is just concerned about linear term and lower nonlinear terms (the second-order and third-order ones).

Since:

$$\left|b_i\right|>>\left|d_i\right|,\left|b_i\right|>>\left|c_i\right|,$$

Thereby, it can be obtained:

$$T_1(y_i)\approx b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)$$

(14)

$$\begin{array}{ccc}\hfill T_2(y_i)& \approx & b_i^2{T}_2(x_i)+4b_i{c}_i{T}_1(x_i)+4c_i^2{T}_2(x_i)+4c_i{d}_i{T}_3(x_i)\hfill \\ & & +2b_i{c}_i(T_1(x_i)+T_3(x_i))+2b_i{d}_i{T}_2(x_i)+2c_i{d}_i{T}_1(x_i)\hfill \\ & =& (6b_i{c}_i+2c_i{d}_i)T_1(x_i)+(b_i^2+4c_i^2+2b_i{d}_i)T_2(x_i)+(4c_i{d}_i+2b_i{c}_i)T_3(x_i)\hfill \\ & \approx & 6b_i{c}_i{T}_1(x_i)+b_i^2{T}_2(x_i)+2b_i{c}_i{T}_3(x_i)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill T_3(y_i)& \approx & [c_i+b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)][(12b_i{c}_i+4c_i{d}_i)T_1(x_i)\hfill \\ & & +(2b_i^2+8c_i^2+4b_i{d}_i)T_2(x_i)+(8c_i{d}_i+4b_i{c}_i)T_3(x_i)\hfill \\ & & +2b_i{d}_i{T}_4(x_i)+2c_i{d}_i{T}_5(x_i)+6c_i^2+2b_i^2+2d_i^2-3]\hfill \\ & \approx & [c_i+b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i)](12b_i{c}_i{T}_1(x_i)+2b_i^2{T}_2(x_i)+4b_i{c}_i{T}_3(x_i)\hfill \\ & & +2b_i{d}_i{T}_4(x_i)+2c_i{d}_i{T}_5(x_i)+2b_i^2-3)\hfill \\ & \approx & 12b_i{c}_i^2{T}_1(x_i)+2b_i^2{c}_i{T}_2(x_i)+4b_i{c}_i^2{T}_3(x_i)+6b_i^2{c}_i{T}_2(x_i)+b_i^3{T}_1(x_i)+b_i^3{T}_3(x_i)\hfill \\ & & +2b_i^2{c}_i{T}_2(x_i)+b_i^2{d}_i{T}_3(x_i)+(2b_i^2-3)b_i{T}_1(x_i)+6b_i{c}_i^2{T}_1(x_i)+6b_i{c}_i^2{T}_3(x_i)\hfill \\ & & +2b_i{c}_i^2{T}_1(x_i)+b_i{c}_i{d}_i{T}_2(x_i)+c_i^2{d}_i{T}_3(x_i)+(2b_i^2-3)c_i{T}_2(x_i)+6b_i{c}_i{d}_i{T}_2(x_i)\hfill \\ & & +b_i^2{d}_i{T}_1(x_i)+b_i{d}_i^2{T}_1(x_i)+c_i{d}_i^2{T}_2(x_i)+(2b_i^2-3)d_i{T}_3(x_i)\hfill \\ & \approx & (2b_i^3-3b_i)T_1(x_i)+(12b_i^2{c}_i-3c_i)T_2(x_i)+(b_i^3-3d_i)T_3(x_i)\hfill \end{array}$$

(16)

Substituting Equations (14)–(16) into Equation (11), the approximate value of the post-compensation signal $z_i$ is expressed as:

$$\begin{array}{ccc}\hfill z_i& \approx & \frac{b_1^4-3b_1{d}_1}{b_1^5}(b_i{T}_1(x_i)+c_i{T}_2(x_i)+d_i{T}_3(x_i))\hfill \\ & & -\frac{b_1^2{c}_1}{b_1^5}(6b_i{c}_i{T}_1(x_i)+b_i^2{T}_2(x_i)+2b_i{c}_i{T}_3(x_i))\hfill \\ & & -\frac{b_1{d}_1}{b_1^5}[(2b_i^3-3b_i)T_1(x_i)+(12b_i^2{c}_i-3c_i)T_2(x_i)+(b_i^3-3d_i)T_3(x_i)]\hfill \\ & \approx & \frac{b_i}{b_1}{T}_1(x_i)+\frac{b_1^2{c}_i-b_i^2{c}_1}{b_1^3}{T}_2(x_i)+\frac{b_1^3{d}_i-b_i^3{d}_1}{b_1^4}{T}_3(x_i)\hfill \end{array}$$

(17)

Normally, the approximate value of $b_1$ is 1. Thereby, Equation (18) could be represented as:

$$z_i\approx b_i{T}_1(x_i)+\frac{b_1^2{c}_i-b_i^2{c}_1}{b_1^3}{T}_2(x_i)+\frac{b_1^3{d}_i-b_i^3{d}_1}{b_1^4}{T}_3(x_i)$$

(18)

Equation (18) allows for the conclusion that the nonlinear compensation process just suppresses the distortion product without affecting linear components which contain the weak signals.

It can be known from Equations (1) and (17) that the uncompensated amplitudes of the second-order and third-order terms are $c_i$ and $d_i$, respectively in the ith channel; while the post-compensation amplitudes of the second-order and third-order terms are $\frac{b_1^2{c}_i-b_i^2{c}_1}{b_1^3}$ and $\frac{b_1^3{d}_i-b_i^3{d}_1}{b_1^4}$, respectively. Therefore, the improvement degrees of the second-order and third-order nonlinear distortion components in the $i$th channel, $I_{i2}$ and $I_{i3}$ are obtained:

$$\begin{array}{cc}\hfill I_{i2}& =10\lg \frac{{(c_i)}^2}{{(\frac{b_1^2{c}_i-b_i^2{c}_1}{b_1^3})}^2}=20\lg \left|\frac{c_i}{\frac{b_1^2{c}_i-b_i^2{c}_1}{b_1^3}}\right|\hfill \\ & =20\lg \left|\frac{b_1^3{c}_i}{b_1^2{c}_i-b_i^2{c}_1}\right|\hfill \\ \hfill I_{i3}& =10\lg \frac{{(d_i)}^2}{{(\frac{b_1^3{d}_i-b_i^3{d}_1}{b_1^4})}^2}=20\lg \left|\frac{d_i}{\frac{b_1^3{d}_i-b_i^3{d}_1}{b_1^4}}\right|\hfill \\ & =20\lg \left|\frac{b_1^4{d}_i}{b_1^3{d}_i-b_i^3{d}_1}\right|\hfill \end{array}$$

(19)

where, “lg” represents Naperian logarithm in base 10.

$G_i$ is defined as the ith channel’s gain. $P{E}_i$ means the difference between $G_i$ and $G_1$ (the first channel’s gain) and it is represented by:

$$P{E}_i=G_i-G_1$$

(20)

and:

$$\begin{array}{l}{G}_{i1}=20\lg \left|\frac{b_i}{b_1}\right|=P{E}_i\\ G_{i2}=20\lg \left|\frac{c_i}{c_1}\right|\\ G_{i3}=20\lg \left|\frac{d_i}{d_1}\right|\end{array}$$

(21)

where, $G_{i1}$, $G_{i2}$ and $G_{i3}$ are the linear, second-order and third-order components of $P{E}_i$, respectively.

The authors make a group of experiments on the values of $P{E}_i$ under different signal frequencies. The experimental results are from a widely-used 16-channel array receiver which is also applied to the nonlinear blind compensation processing of the array system. Among the values of $P{E}_i$ for each frequency point, for statistical convenience, only the maximum and minimum absolute values are listed in

Table 1. The maximum and minimum absolute values of PE_i versus different frequency points.

Signal Frequency (MHz)	Computational Complexity
	Maximum |PEi| (dB)	Minimum |PEi| (dB)
5.02	PE14 = 1.166	PE5 = 0.143
13.92	PE10 = 0.583	PE3 = 0.008
16.97	PE5 = 0.582	PE6 = 0.022
19	PE10 = 0.697	PE3 = 0.029
21.6	PE10 = 0.318	PE2 = 0.019
28.02	PE5 = 0.976	PE10 = 0.04

where the absolute values of $P{E}_i$ range from approximately 0 to near 1. Consequently, in the real array receiving system, the reasonable estimation that $G_{i1}$ ranges from −1 dB to 1 dB can be obtained. Moreover, as shown in Table 1, most of the measurements are much smaller than 1. Therefore, for computational convenience, $G_{i2}$ and $G_{i3}$ both ranging from −1 dB to 1 dB is rational.

Based on Equations (19) and (21), the authors obtain the calculation results for $I_{i2}$ and $I_{i3}$, as shown in Figure 2, which demonstrates that $I_{i2}$ is approximately 16–32 dB and $I_{i2}$ is about 12–24 dB. Therefore, the results prove that this novel method could significantly compensate each channel’s nonlinear distortion of array receiver.

3. Proposed Mitigation Architecture for Array Receiver

It is the principle of the proposed algorithm that large signals create significant distortion product. The objective function is written as the minimal residual sum of squares between small signals and the nonlinear model of large signals.

This paper proposes, first, the blind compensation architecture in detail followed by the proposed IWI method for calculating nonlinear model parameters. Second, a SRA-TFC algorithm is given to separate time-domain large and small signals. Finally, the authors describe the comparisons to traditional SVD-based algorithms regarding the compensation performance and the computational complexity.

3.1. Nonlinear Blind Mitigation Structure

Figure 3 demonstrates the proposed structure implemented in the blind compensation unit is made up of two basic modules of identification and compensation. The model parameters of the first-channel are extracted adaptively.

To dramatically improve the SFDR of the multi-channel array receiver, the nonlinear behavior model of the first channel RF front-end can be approximately expressed by the Hammerstein model, which is a cascade of static nonlinearity blocks followed by linear dynamical blocks [33]. The relationship between input and output can be expressed as:

$$v(n)={\displaystyle \sum _{d=1}^D{\displaystyle \sum _{r_d=0}^{N_d-1}h(r_d,d){L_1}^d(n-r_d)}}$$

(22)

where $v(n)$ is the output signals of the discrete Hammerstein model; $L_1(n)$ means large signals in the first channel; $D$ expresses the maximum value of the nonlinear order $d(1\le d\le D)$; $N_d$ represents the memory depth of the $d$th order Hammerstein kernel; $h(r_d,d)$ is the model kernel coefficient of the $d$th order with memory depth $r_d$; the total number of the Hammerstein kernel is:

$$P={\displaystyle \sum _{d=1}^D{N}_d}$$

(23)

Furthermore, the major function of the identification module is to generate identification parameters of a nonlinear model according to the objective function minimizing the power of $N$ points error signal $e(n)$ where:

$\mathit{\omega }(n)={[h(0,1)h(1,1)\dots h(N_1-1,1)h(0,2)h(1,2)\dots h(N_D-1,D)]}^T$ are the kernel coefficients of the nonlinear model whose original value is an $P$-dimensional zero column vector.

The nonlinear model kernel $\mathit{u}(n)$ is a $P$-dimensional column vector produced by large signals and denoted as:

$$\begin{array}{cc}\hfill \mathit{u}(n)=& [L_1(n)L_1(n-1)\dots L_1(n-N_1+1)\hfill \\ & {(L_1(n))}^2{(L_1(n-1))}^2\dots {(L_1(n-N_2+1))}^2\hfill \\ & ........\hfill \\ & {(L_1(n))}^D{(L_1(n-1))}^D\dots {(L_1(n-N_D+1))}^D{]}^{\mathrm{T}}\hfill \end{array}$$

(24)

and:

$$J(\mathit{\omega })={[\mathit{e}(n)]}^{\mathrm{T}}\mathit{e}(n)={\displaystyle \sum _{k=1}^N{[e(k,\mathit{\omega })]}^2}$$

(25)

where, $e(k,\mathit{\omega })$ stands for the $k$th point $\mathit{e}(n)$; $\mathit{e}(n)={[e(1,\mathit{\omega }),e(2,\mathit{\omega })\dots e(N,\mathit{\omega })]}^{\mathrm{T}}$ ($1\le k\le N$) is $N\times 1$; $J(\mathit{\omega })$ means the power of $N$ points $\mathit{e}(n)$.

$$e(k,\mathit{\omega })=S_1(n-k+1)-{\mathit{u}}^{\mathrm{T}}(n-k+1)\mathit{\omega }(n)$$

(26)

where, $S_1(n-k+1)$ and $\mathit{u}(n-k+1)$ signify small signals and nonlinear model kernel of the $k$th point in the first channel, respectively.

It can be obtained from Equations (25) and (26):

$$\begin{array}{ccc}\hfill J(\mathit{\omega })& =& {\displaystyle \sum _{k=1}^N{[e(k,\mathit{\omega })]}^2}={\displaystyle \sum _{k=1}^N{[e(k,\mathit{\omega })]}^{\mathrm{T}}e(k,\mathit{\omega })}\hfill \\ & =& {\displaystyle \sum _{k=1}^N{[S_1(n-k+1)-{\mathit{u}}^{\mathrm{T}}(n-k+1)\mathit{\omega }]}^{\mathrm{T}}}[S_1(n-k+1)-{\mathit{u}}^{\mathrm{T}}(n-k+1)\mathit{\omega }(n)]\hfill \\ & =& {\displaystyle \sum _{k=1}^N{S}_1{}^{\mathrm{T}}(n-k+1)S_1(n-k+1)}-{\displaystyle \sum _{k=1}^N{S_1}^{\mathrm{T}}(n-k+1){\mathit{u}}^{\mathrm{T}}(n-k+1)\mathit{\omega }}(n)\hfill \\ & & -{\displaystyle \sum _{k=1}^N{\mathit{\omega }}^{\mathrm{T}}\mathit{u}(n-k+1)S_1(n-k+1)}+{\displaystyle \sum _{k=1}^N{\mathit{\omega }}^{\mathrm{T}}\mathit{u}(n-k+1){\mathit{u}}^{\mathrm{T}}(n-k+1)\mathit{\omega }}(n)\hfill \end{array}$$

(27)

The objective function is to minimize $J(\mathit{\omega })$, therefore, the authors compute the gradient of $J(\mathit{\omega })$ with respect to $\mathit{\omega }(n)$:

$$\begin{array}{l}-{\displaystyle \sum _{k=1}^N{[{S_1}^{\mathrm{T}}(n-k+1){\mathit{u}}^{\mathrm{T}}(n-k+1)]}^{\mathrm{T}}}-{\displaystyle \sum _{k=1}^N\mathit{u}(n-k+1)S_1(n-k+1)}\\ +2[{\displaystyle \sum _{k=1}^N\mathit{u}(n-k+1){\mathit{u}}^{\mathrm{T}}(n-k+1)]{\mathit{\omega }}_{opt}(n)}=0\end{array}$$

(28)

Continuing to earn Equation (29):

$$\mathit{K}(n){\mathit{K}}^{\mathrm{T}}(n){\mathit{\omega }}_{opt}(n)=\mathit{K}(n)\mathit{Y}(n)$$

(29)

where $\mathit{Y}(n)={[S_1(n),S_1(n-1),S_1(n-2),\dots S_1(n-N+1)]}^{\mathrm{T}}$ is $N\times 1$, $\mathit{K}(n)=[\mathit{u}(n),\mathit{u}(n-1),\mathit{u}(n-2),\dots \mathit{u}(n-N+1)]$ is $P\times N$ and the kernel vector of the nonlinear behavior model ${\mathit{\omega }}_{opt}(n)$ means the solution of the objective function.

Simplifying Equation (29) to earn Equation (30):

$$\mathit{K}{\mathit{K}}^{\mathrm{T}}{\mathit{\omega }}_{opt}=\mathit{K}\mathit{Y}$$

(30)

where, $\mathit{K}(n)$, ${\mathit{\omega }}_{opt}(n)$ and $\mathit{Y}(n)$ is abbreviated to $\mathit{K}$, ${\mathit{\omega }}_{opt}$ and $\mathit{Y}$, respectively.

Last, based on the proposed compensation strategy, the authors acquire:

$$\begin{array}{cc}{z}_i=y_i-{{\mathit{U}}_i}^{\mathrm{T}}{\mathit{\omega }}_{opt}& 1\le i\le m\end{array}$$

(31)

where ${\mathit{U}}_i$ (the facilitation form of ${\mathit{U}}_i(n)$) denotes the pipeline model kernel of the $i$th channel, while ${\mathit{U}}_i(n)$ is expressed as:

$$\begin{array}{cc}\hfill {\mathit{U}}_i(n)=& [y_i(n)y_i(n-1)\dots y_i(n-N_1+1)\hfill \\ & {(y_i(n))}^2{(y_i(n-1))}^2\dots {(y_i(n-N_2+1))}^2\hfill \\ & ........\hfill \\ & {(y_i(n))}^D{(y_i(n-1))}^D\dots {(y_i(n-N_D+1))}^D{]}^{\mathrm{T}}\hfill \end{array}$$

(32)

The specific procedures to resolve Equation (30) are discussed as the following:

The orders of the Hammerstein series are set from 1 to 5 in Equation (30). Due to a very great power difference among various order terms, obviously $\mathit{K}{\mathit{K}}^{\mathrm{T}}$ is considered as an ill-conditioned matrix. Consequently, the calculation of the model kernel vector can be converted into that of ill-conditioned linear equations. Here, the IWI method is present to calculate the ill-conditioned equations. Representing a given iteration $l$, the iterative formula is:

$$\begin{array}{cc}\left(\mathit{A}+\lambda \mathit{I}\right){\mathit{X}}^{\left(l+1\right)}=\mathit{b}+\lambda {\mathit{X}}^{\left(l\right)}& (\mathit{A}=\mathit{K}{\mathit{K}}^{\mathrm{T}}\mathit{b}=\mathit{K}\mathit{Y})\end{array}$$

(33)

where $\lambda$ represents a non-zero constant, $l$ expresses the iteration index, ${\mathit{X}}^{\left(0\right)}$ is the initial value of the equations solution, ${\mathit{X}}^{\left(l\right)}$ is the approximate estimate at $l$th iteration.

Figure 4 illustrates the concrete steps of IWI algorithm, where ${\mathit{A}}_2=\mathit{A}+\lambda \mathit{I}$, ${\mathit{r}}^{\left(l\right)}$ and ${\mathit{e}}^{(l)}$ indicate the margin and error at iteration $l$, respectively, $l++$ means $l$ increasing by 1 in every iteration.

3.2. SRA-TFC Method

The key point of the proposed nonlinear blind compensation method is to gain large signals $L_1(n)$ and small signals $S_1(n)$ blindly and precisely. The following are the details regarding the proposed blind signal separation method.

According to Figure 5, the SRA-TFC algorithm can be divided into the following four steps:

• The time domain distorted signal $y_1\left(n\right)$ and pure noise $n_1\left(n\right)$ from the first channel RF front-end needs to be extracted.
• $y_1\left(n\right)$ and $n_1\left(n\right)$ are both translated into the corresponding signal and noise in frequency-domain via the Discrete Fourier transform (DFT) technique and denoted as $Y_s\left(f\right)$ and $Y_n\left(f\right)$, respectively.
• Each point of $Y_s\left(f\right)$ is compared with the setting power spectrum threshold. Next, frequency-domain large signals $L_1\left(f\right)$ are obtained with the method that the points in $Y_s\left(f\right)$, whose values are below the threshold, are replaced by the corresponding points in $Y_n\left(f\right)$, whereas the others remain. Conversely, frequency-domain small signals, $S_1\left(f\right)$ is gained by the replaced points in $Y_s\left(f\right)$ whose values are above the threshold.
• By the IDFT technique, the authors turn $L_1\left(f\right)$ and $S_1\left(f\right)$ to time-domain large signals, $L_1(n)$ and time-domain small signals, $S_1(n)$, separately.

3.3. Comparisons between SRA-TFC Method and Traditional SVD-Based Method

To demonstrate the advantages of the proposed method, the authors compare the SRA-TFC signal separation method to the traditional SVD-based signal separation method in terms of the compensation performance and the computational complexity.

3.3.1. Comparison of Compensation Results

The simulation compensation results are provided to testify the advantages of the SRA-TFC signal separation method. To focus on the comparison of these two signal separation methods, the parameters of any one of the channels are adaptively extracted to alleviate the nonlinearity of this channel. It is noteworthy that the number of the sampling points for power spectrum analysis is 65,536 in Figure 6 and Figure 7.

Figure 6a depicts that the original simulation signals cover the equal-amplitude three-tone strong sinusoidal wave signal (large signals) and the weak 16-QAM (quadrature amplitude modulation) modulation signal (target signal). The frequencies of large signals are 10.3 MHz, 10.6 MHz and 11.8 MHz and the center frequency and bandwidth of the weak modulation signal is 21.5 MHz and 1 MHz, respectively. Moreover, the power spectrum difference between large signals and the target signal is about 75 dB and the sample rate is 70 MHz. A Hammerstein model with D = 5 and N_d = [2, 2, 2, 2, 2] is adopted. As shown in Figure 6b, the signal and carrier frequencies are chosen arbitrarily, the only condition being that some large nonlinear distortion components should be in the target signal region to dramatically reduce the SINR of the target signal.

Regarding Figure 7, the blue spectrum illustrates the compensation performance when the SRA-TFC signal separation method is employed, while the red spectrum indicates the mitigation result in the application of the SVD-based signal separation method. While both signal separation methods perform adequately, the proposed method demonstrates slightly better overall mitigation performance.

Furthermore, as depicted in Figure 6b, the SINR of the weak 16-QAM signal is dramatically reduced due to the large disturbance stemming from the large signals. Moreover, after the processes of frequency selection, filtering and demodulation for the weak 16-QAM signal, its constellation is very messy, as shown in Figure 8a. However, as illustrated in Figure 8b,c, benefiting from the proposed algorithm, the nonlinear distortion in the target signal region is obviously depressed. Finally, a regularly-distributed constellation is obtained. Moreover, compared to Figure 8b, the constellation in Figure 8c is more regular because the novel algorithm has better compensation performance in the target signal region.

This experiment illustrates that the nonlinear blind compensation method can improve the SINR of the weak 16-QAM signal but not change the weak signal itself, which is a crucial capacity of the parameter estimation for target signals.

3.3.2. Comparison between Computational Complexities

The computational complexities of the proposed method and the traditional SVD-based algorithm are discussed in this section. Let the data length $N=M\times M$ where $M$ is applied to be the measurement base of complexity.

Separating large signals and small signals from the one-channel output signal, the required calculations of the traditional SVD-based algorithm include $O\left(M^3\right)$ times addition, $O\left(M^3\right)$ iterations of multiplication, $O\left(M^2\right)$ iterations of division and $O\left(M^2\right)$ times square-root, to the knowledge of Reference [27].

According to the procedures described in Section 3.2, the required computational complexity of the proposed method is displayed as the following:

• $N$ points $y_1\left(n\right)$ and $N$ points $n_1\left(n\right)$ are first extracted and then separately converted into $Y_s\left(f\right)$ and $Y_n\left(f\right)$ via DFT. The $N$ points DFT requires both $N\log N$ times of multiplication and $N\log N$ times of addition. Therefore, the required calculations to acquire $Y_s\left(f\right)$ and $Y_n\left(f\right)$ are $2N\log N$ iterations of multiplication and addition.
• It requires $2N$ times multiplication and $N$ times addition to calculate the power spectral density of $Y_s\left(f\right)$. Moreover, with the application of above threshold detection method, it needs $2N$ iterations of addition to achieve $L_1\left(f\right)$ and $S_1\left(f\right)$.
• The authors separately convert $L_1\left(f\right)$ and $S_1\left(f\right)$ into $L_1(n)$ and $S_1(n)$ by IDFT. This step requires $2N\log N$ iterations of multiplication and addition.

According to Steps 1–3, the proposed signal separation method requires $O\left(N\right)=O\left(M^2\right)$ iterations of multiplication and addition. The authors summarize the computational complexities of the two methods in

Table 2. Comparisons of the computational complexities between two signal separation methods.

Signal Separation Method	Computational Complexity
	Addition	Multiplication	Division	Square-Root
traditional SVD-based method	$O\left(M^3\right)$	$O\left(M^3\right)$	$O\left(M^2\right)$	$O\left(M\right)$
SRA-TFC method	$O\left(M^2\right)$	$O\left(M^2\right)$	without	without

. It can be observed that the computational complexity of the proposed method is $O\left(M\right)$ times less than that of SVD-based method.

Compared to the traditional SVD-based algorithm, in summary, the proposed method has much lower computational loads but achieves a slightly better mitigation performance.

4. Experimental Results and Analysis

An experimental verification platform is built to validate the performance of the proposed algorithm. Moreover, the real-world experimental results show the impact of the nonlinearity mitigation of the array receiver on the source localization accuracy of target signals and evaluate the performance of the proposed algorithm.

4.1. Nonlinearity Mitigation Performance for ASP System

Relative to the other array geometries, uniform circular array (UCA) has several advantages such as its 360° azimuth angle coverage and additional elevation angle information [34,35]. Additionally, it has wide application in 2-D DOA estimation. Therefore, according to Figure 1, the shortwave array signal used for testing is from an actual UCA with sixteen elements. The experimental array receiver is a 16-channel receiver and its synchronous sampling rate is $f_s=70\mathrm{MHz}$.

The model parameters of the first channel are extracted and then to alleviate the nonlinear distortion components of all channels in this experiment. Additionally, the nonlinear inverse behavior model for all channels’ RF front-ends is employed by a Hammerstein model with D = 5 and N_d = [3, 3, 3, 3, 3]. Figure 9a,b show the power spectrum of the original and post-compensation array signal in the first channel. Since the input array signals of all channels are similar and their compensation results are also similar, it is only necessary to show the power spectrum of the signal before and after compensation in other channels such as the second one, which are shown in Figure 9c,d.

Moreover, Figure 9a,b illustrates that the SFDR of the first channel has a 10–20 dB improvement in the whole band by carrying out the presented compensation method. Meanwhile, the same conclusion of the second channel from Figure 9c,d can be drawn. Therefore, it can be deduced that the SFDR of the entire array receiver in in the full band rises by 10–20 dB through the implementation of the proposed compensation method.

Meanwhile, three target signals (two weak ones and a strong one) are chosen to dissect the effect of the nonlinearity mitigation of array signal processing system on the 2-D DOA estimation performance of array signal. Displayed in Figure 10a–c (the local power spectrum graphs of the first channel), there are three array signals of which center frequencies are 16.97 MHz, 19 MHz and 13.92 MHz, respectively. Meanwhile, the powers of the two array signals in Figure 10a,b are both almost 20 dB lower than the power of the array signal in Figure 10c. Therefore, the three targets are named weak target signal-1, weak target signal-2 and strong target signal, respectively.

Figure 10a,b illustrate that the two weak array signals are both seriously interfered by unwanted in-band nonlinearity and that the SINRs of the two signals are below 0 dB. However, after compensation, the SINRs of the two signals are both increased by approximately 10 dB. Therefore, the two array signals are the typical examples the authors research.

Contrarily, as shown in Figure 10c, the power spectrums of the strong array signal before and after compensation are almost similar, owing to a total lack of the nonlinear distortion components. Additionally, the SINRs of the original and post-compensation signals are approximately 20 dB.

4.2. 2-D DOA Performances of Target Signals

The direction-finding performances of the above-mentioned array signals before compensation versus after compensation are discussed. Figure 11 shows the 16-channel array receiver output before compensation is processed through the frequency selection and channelization before their 2-D DOAs (elevation angle $\alpha$ and azimuth angle $\beta$) are estimated as is the compensated array receiver output where, the channelized sampling rate and bandwidth are $f_n$ = 10 KHz and 8 KHz, respectively. Moreover, the authors employ the MUSIC algorithm to estimate their 2-D DOAs of array signal before and after compensation. Ultimately, the authors provide explicit comparisons about these two estimation results.

The array receiver output with 118 s in succession is used to estimate the 2-D DOA parameters based on the MUSIC algorithm. Where one-second array signal data is estimated at a time. Thereupon, one hundred and eighteen independent Monte Carlo trials are performed for the following experiments. Moreover, the root mean square error (RMSE) is used for the performance measure. The RMSE of elevation and azimuth DOA ($RMSE(\alpha )$ and $RMSE(\beta )$) are defined as:

$$RMSE(\alpha )=\sqrt{\frac{1}{118}{\displaystyle \sum _{t=1}^{118}{({\alpha }_t-\tilde{\alpha })}^2}}$$

(34)

$$RMSE(\beta )=\sqrt{\frac{1}{118}{\displaystyle \sum _{t=1}^{118}{({\beta }_t-\tilde{\beta })}^2}}$$

(35)

where ${\alpha }_t$ and ${\beta }_t$ mean the estimation of the elevation angle and azimuth angle at the $t$th second (the $t$th Monte Carlo trial), separately, while $\tilde{\alpha }$ and $\tilde{\beta }$ stand for the corresponding average ones.

4.2.1. 2-D DOA Performance in the Case of Weak Target Signals

These two experiments are to compare the 2-D DOA estimation results of weak target signals before and after compensation.

Figure 12 depicts the DOA estimation performances of weak target signal-1 before and after compensation by 118 trials and demonstrates that the estimation result of the post-compensation signal is much more accurate than that of the original signal. Meanwhile, the same phenomenon also happens on weak target signal-2, as shown in Figure 13. These phenomena are observed more directly from the azimuth-elevation scatter figures, Figure 14 and Figure 15.

Moreover, by carrying out the compensation algorithm, the $RMSE(\alpha )$ of weak target signal-1 is reduced from 15.62° to 3.31° and the $RMSE(\beta )$ is decreased from 65.81° to 10.03°. Concurrently, for weak target signal-2, its $RMSE(\alpha )$ and $RMSE(\beta )$ are reduced from 25.66° to 3.94° and from 34.98° to 2.48° by employing the blind compensation algorithm.

The above two experimental results indicate that the proposed algorithm can dramatically increase the 2-D DOA estimation accuracy of the weak target signals.

4.2.2. 2-D DOA Performance in the Case of Strong Target Signal

Another crucial concern to the authors is the azimuth-elevation direction finding performance analysis of strong array signal before compensation versus after compensation.

It can be seen in Figure 16 that the azimuth-elevation angle estimation results of strong target signal before and after compensation are identical. The experimental result suggests that the proposed method just eliminates the nonlinearity without affecting the strong signal and its parametric estimation accuracy.

5. Conclusions

Rather than developing different parametric estimation methods, the authors investigated a nonlinear blind compensation technique in array signal processing. A novel blind compensation algorithm was presented to improve the SFDR of an array receiver and to obtain remarkable improvements of the parameters estimation performances for weak target signals such as their 2-D DOAs employing the MUSIC algorithm. These purposes have already been verified by the repeated experiments on a real-world array signal from a UCA. The main advantages of the proposed algorithm are summed up in the following points:

• During the blind compensation process, the parameters of the identification module and the compensation module are totally independent of each other, which can improve the efficiency of array signal processing and increase its dependability.
• The suggested algorithm accomplished in pure time-domain saves a great deal of the system hardware scale. In addition, the proposed algorithm is only necessary to set power threshold rather than multi-stopband/multi-passband digital filters with extremely high performance of constant change, which is apparently more flexible and convenient to handle the situation of multiple signals with different power levels or wide ranges of bandwidth.
• The blind compensation strategy for multi-channel RF front-ends designed is that the model parameters of any one of the channels are extracted to mitigate the nonlinear distortion components of all channels synchronously. It has the advantages of reducing a mountain of computational loads and avoiding the inconsistency of the array compensation performance caused by the iterative computation error, especially the phase disturbance of the array signal after compensation.

In conclusion, the proposed algorithm can improve the SFDR of the array receiver, achieving high parameters estimation accuracy of array signal with low computational complexity. It is predicted to have wide and promising applications in array signal processing.