Air-Induced Passive Intermodulation in FDD MIMO Systems: Algorithms and Measurements

In this article, we model and develop effective cancellation schemes for passive intermodulation (PIM) distortion induced by external objects in the vicinity of the transmitter—referred to as air-induced PIM—in frequency-division duplex (FDD) multiple-input–multiple-output (MIMO) systems. PIM is interference generated from the transmitted signals undergoing nonlinear transformation, which, in FDD systems, may cause desensitization of the receiver chain. First, we present a general model of the received air-induced PIM signal, with an arbitrary number of dual-carrier TX chains active and an arbitrary number of PIM sources causing interference, on all the intermodulation frequencies. Then, the derived model is used to develop a cancellation scheme based on a complete set of basis functions (BFs) in a rank-2 dual-carrier MIMO system with four active carriers. To alleviate the high complexity of the aforementioned scheme, we then propose a novel alternative cancellation scheme with much reduced complexity, leveraging on the physical modeling of the system, which is capable of handling any number of PIM sources in the system. RF measurement-based experimentations carried out with real-life equipment evidence excellent cancellation capabilities of the complete BF model, which can be retained with much reduced complexity with the proposed alternative technique.

In the lower end of the available spectrum, most frequency range 1 (FR1) bands utilize frequency-division duplex (FDD) as the duplexing scheme, as is the case, for example, with 5G NR band n3 [1], [2]. In order to avoid self-interference (SI) directly from the transmitter (TX) signals, the transmission and reception are carried out on different center frequencies in FDD systems. However, such a system is now possibly vulnerable to interference from PIM distortion, as depicted in Fig. 1. If the receiver (RX) operates on or near the intermodulation frequencies, the PIM products can fall on the RX passband, which can possibly desensitize the RX path by drowning the desired received signals from user equipment (UE). All FDD transceivers incorporate a duplexer between the TX and RX chains, as depicted in Fig. 2, which protects the RX path from intermodulation stemming from the TX path. However, any PIM generated after the duplexer (e.g., in connectors, antennas, or nearby metallic objects) can still desensitize the RX path. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Fig. 2. Simplified block diagram presentation of the considered system, where N t transceivers located in a BS simultaneously transmit two aggregated CCs, per transceiver, and receive nonlinear interference stemming from N P air-PIM sources. Some relevant modeling notations are also shown. The PIM canceller block can be implemented, as described in Section III.
Pragmatic suppression of the PIM distortion can be achieved by backing off the TX power; however, such will greatly deteriorate the system power efficiency and the corresponding network coverage and is, therefore, not a congenial solution.
Since PIM is commonly generated by the internal circuitry after the duplexer of the transceiver, the simplest solution to such conducted-PIM mitigation is to use components that are highly linear. This unfortunately increases the manufacturing costs of the devices, and it is, therefore, not an appealing option [5], [7], [9], [10]. In the air-induced PIM scenarios, the obvious solution would be to manually remove the PIM source from the radiation field of the TX, which, in some cases, could turn out to be impossible, if, for example, the PIM source is part of the built environment. Thus, PIM mitigation techniques have been studied in the literature, and they continue to hold scientific and practical relevance in the 5G era. These are shortly reviewed next.

A. Prior Art
In general, the PIM cancellation methods can be divided into analog and digital domain suppressors. In the analog domain, a preliminary work in [5] introduced a nonlinear imposer network, which was shown to suppress PIM by up to 40 dB, while, in [6], an engineered PIM source was used for mitigating PIM. Furthermore, the work in [22] proposed a framework for the adaptation of the modeling and consecutive suppression of high-order PIM distortion. In [7], simulations were carried out to validate PIM mitigation methods by introducing a compensation signal or an additional passive PIM source accompanied by a phase shifter. The works in [23] and [24], in turn, seek to suppress PIM by microstrip coating. Recently, Chen et al. [8] proposed a tunable nonlinear transmission line strip to cancel PIM.
In the digital domain, in [9] and [10], the conducted PIM of two CC systems was modeled with a simple third-order polynomial model and subsequently canceled. The work in [11] modeled and canceled the second-order PIM products of a diplexer. The potential coexistence of nonlinear power amplifiers (PAs) and PIM distortion is considered in [12] and [18], where reduced complexity models through decoupling of the nonlinearity and memory were introduced. In [16], basic modeling and digital cancellation of air-induced PIM were pursued in a simple two-CC scenario. Finally, our preliminary work in [17] addresses the air-induced PIM cancellation in a multiple-input-multiple-output (MIMO) scenario with four simultaneous carriers.
In addition to PIM modeling and cancellation, the works in [3], [4], [25], [26], [27], [28], [29], and [30] have considered the cancellation of the intermodulation distortion imposed by active components, primarily the TX PA systems. However, these works do not consider PIM, yet it is noted that similar behavioral modeling principles apply when the source of the intermodulation distortion is an active component. Intermodulation distortion can also be mitigated through digital predistortion (DPD), where the inverse of the nonlinear block is applied to the signal digitally, thus linearizing the TX output [31], [32]. In addition, it is mentioned here that inband full-duplex (IBFD) technology [33], [34], [35], where the RX operates simultaneously and on the same frequency as the TX, is dependent on the mitigation of the so-called SI signal, which, without mitigation, desensitizes the RX. In IBFD transceivers, the SI signal is present even if the entire system is completely linear. Therefore, in the IBFD case, one must model the RX signal at the TX frequency, while, in FDD PIM cancellation, the signals are modeled at the intermodulation frequencies.

B. Contributions and Novelty
In this article, building on our initial work in [17], we focus on the modeling and digital cancellation of air-induced PIM stemming from the interaction of the radiated electromagnetic wave with external materials outside the transceiver system. Specifically, we consider a challenging rank-2 dual-carrier MIMO FDD system where four CCs are transmitted simultaneously. Besides our initial work in [17], the authors are not aware of any works explicitly aiming to model and cancel air-induced PIM in MIMO scenarios although the issue of PIM has been mentioned in MIMO literature, e.g., in [36], [37], and [38]. To this end, the so-called full basis function (BF) canceller solution is first described. Then, to alleviate the involved complexity, we also propose an alternative, novel yet simple cancellation method, based on the physical modeling of the system. This method is first described for one dominating air-PIM source, introduced in [17], while it is then also extended to cover multiple simultaneous sources of air-induced PIM. All the models are validated through RF measurements at 5G NR band n3 with true base-station (BS) hardware, evidencing excellent PIM suppression capabilities.
Compared to our initial work in [17] and the overall available literature in the field, the main novelties and contributions of this article can be summarized as follows.
1) Rigorous mathematical modeling of the system, where N P over-the-air (OTA) PIM sources cause nonlinear interference to the RX path, is presented. 2) The complete BF method, introduced initially in [17], is formulated and presented in a complete manner, covering both the cancellation processing and the parameter estimation; 3) In order to alleviate the modeling and processing complexity, an alternative cancellation scheme leveraging the air-PIM channel coefficients is presented, which, unlike the initial channel coefficient method in [17], can suppress PIM stemming from multiple coexisting air-PIM sources. 4) The complexity of each introduced cancellation scheme is studied and quantitatively assessed in terms of floating-point operations (FLOPs). 5) Extensive RF measurements are provided with true BS hardware, evidencing that the proposed methods can efficiently cancel air-induced PIM toward the thermal noise floor, and also that the simplified method can essentially match the cancellation performance of the much more complex full BF model in various cases. The rest of this article is organized as follows. First, in Section II, we provide the modeling fundamentals of the FDD MIMO transceiver under air-induced PIM. Then, in Section III, the proposed digital cancellation methods are described, including the core cancellation processing and the related parameter estimation methods, while also assessing the involved processing complexity. Section IV provides an extensive set of RF measurement-based results. Finally, conclusions are drawn in Section V, while Appendices I and II provide selected technical details related to the BFs and the gradients utilized in the algorithms. Mathematical Notations: Standard complex baseband signal modeling is adopted throughout this article. Time-dependent variables are denoted with a generic time index t in brackets. The imaginary unit is denoted as j . Boldfaced characters indicate vectors and matrices. The complex conjugate and Hermitian transpose are denoted as (·) * and (·) H , respectively. Finally, denotes the convolution operation.

II. FDD MIMO SYSTEM UNDER AIR-INDUCED PIM
In this section, we introduce the overall system model. We consider a dual-carrier MIMO system with N t TXs, each transmitting two CC signals, as shown in Fig. 2. The CC signals are denoted as x 1,n [t] and x 2,n [t], n ∈ {1, 2, . . . , N t }, such that x 1,n [t] are transmitted at frequency f 1 and x 2,n [t] are transmitted at f 2 . For analytical purposes and notational convenience, we assume that f 1 < f 2 . For notational purposes, these frequencies refer to complex baseband equivalent modeling, while the corresponding physical RF frequencies are denoted as f RF,1 and f RF,2 .
Let us assume that there are N P air-PIM sources in the vicinity of the transceiver such that they are time-aligned at the BB sampling rate of the CCs. The transmit signals propagate through a wireless line-of-sight channel to the PIM sources. Assuming a static frequency flat channel response, the combined signal at the input of the PIM source s can be written as where a s,1,n and a s,2,n are the channel coefficients corresponding to the channel responses of the CCs 1 and 2, for TX n, toward the sth PIM source. Following the generalized memory polynomial (GMP) model [39], the sth PIM source output can be written as where M 1 and M 2 are the global precursor and postcursor memory depths, respectively, P 1 and P 2 are the polynomial orders considered in the model, K 1 is the envelope lead, K 2 is the envelope lag, and g is the model coefficients.
The PIM source outputs y s [t] then propagate through wireless channels back to the RX. The received signalr PIM,n [t] at RX n over the whole band can be modeled as where b s,n is the corresponding channel coefficient and η[t] is the additive white Gaussian noise. The total received wideband signal can be seen in (5), as shown at the bottom of the next page, whereγ are the total lumped coefficients of each summed signal ψ, defined in (1) and (2), containing channel coefficients b, the GMP model coefficients g, and the phase shifts stemming from the delays of the memory and envelope lag/lead. Finally, a channel selection filter w n at the RX limits the bandwidth (BW) of the signalr PIM,n [t] to produce a bandlimited version r PIM,n [t] as where denotes the convolution operation. In Section III, building on these fundamental signal models, we concentrate on the derivation and implementation of the PIM cancellation block in Fig. 2.

III. PROPOSED DIGITAL PIM CANCELLATION METHODS
For generality, the modeling of Section II assumed an arbitrary number of transmit chains N t and, hence, signals per CC center frequency, and made no assumptions of the corresponding RX frequency. In this section, in order to develop further the closed-form expressions for the received PIM signal, we consider the practical case of N t = 2. Therefore, there are four physical CC signals present in the system, denoted as x 1,1 [t] and x 1,2 [t], which occupy frequency f 1 , and x 2,1 [t] and x 2,2 [t], which occupy frequency f 2 . In addition, we set the RX frequency to 2 f 1 − f 2 , which is the lower third-order intermodulation frequency of the carriers and which, for example, in the case of 5G NR band n3, can fall on the uplink receive band. Accordingly, we will develop models for the received PIM signal for up to polynomial order P 1 = P 2 = 5. However, for clarity, it is noted that the modeling results in Section II up to (5) can be used to develop the model for any number of transmitted carriers on any RX frequency. The introduced methods can be used in the PIM canceller block illustrated in Fig. 2. In this section, the individual signals are expressed as time series, where the samples are stacked in vectors. The vector-form TX signals are, therefore, defined as , which are the signals at frequency f 1 , and x 2, ]] T ∈ C ×1 , which are the signals at frequency f 2 . Here, corresponds to the signal block length.

A. Full BF Method
For brevity, in this section, we neglect the memory effects in our modeling and utilize only the instantaneous model (m = 0). This is done without loss of generality, as it is essentially trivial to extend the model to incorporate the past and future samples, since the memory from m = 0 affects each signal in the same way. With the considered system parameterization, given above, the bandlimited instantaneous received signal at RX n at time instant t, operating on frequency 2 f 1 − f 2 , can be written as shown in (6), as shown at the bottom of the page, stemming from (5) and (7), where noise has now been omitted for the notational simplicity.
Using then the definitions in (1) and (2), we can write the instantaneous received PIM signal, as shown in (8), at the bottom of the next page, where the signal is shown up to polynomial order p = 3 for brevity. The received PIM signal can then be further defined as expressed in (9), as shown at the bottom of the next page, where the total coefficients of the BFs have been lumped to c and d, such that the coefficients c correspond to BFs, where k = 0, and d to BFs where k = 0. Table IV in Appendix I collects the instantaneous BFs at time instant t for polynomial orders p = 3 and p = 5 for readers' convenience. It is noted that the number of PIM sources in the system does not affect the number of considered BFs. Conversely, depending on the value of the envelope delay k, the system employs a different number of BFs, depending on whether k is zero or nonzero.
Next, the BFs in Table IV can be stacked into matrix ϕ ϕ ϕ[t] ∈ C ×C each BF occupying a single column. Here, C is the cardinality, i.e., the number of instantaneous BFs in the system. The instantaneous BFs are generated from the input vectors , and x 2,2 [t] through elementwise multiplications as per Table IV. Incorporating memory into the model can be carried out at his stage by generating the total Then, the cancellation of the received PIM signal can be written as is the nth RX signal stacked in to a vector, and β β β ∈ C C M×1 contains the estimates of the coefficients c and d of each BF shown in (9). For future reference, following (8) and (9), the coefficients ccorresponding to the instantaneous third-order ( p = 3) BFs of the model-can be explicitly expressed as A block diagram of the full BF scheme is illustrated in Fig. 3. Furthermore, due to the linear-in-parameters nature of the  Table IV, the coefficients β β β can be estimated using least squares (LS) fitting as The required steps to cancel PIM with the full BF method are shown and summarized in Algorithm 1 for a given block of received signal samples. In a static environment, the recalculation of the coefficient vector β β β is not always necessary for each processing block.
Since the complexity of the cancellation engine in (10) and especially that of the coefficients estimation in (12) are dependent on the number of BFs in the system, the complexity of the above proposed scheme is extremely high, especially when memory is considered in the model. Inspired by this, in Sections III-B and III-C, we propose alternative cancellation schemes with much reduced complexity, leveraging on the knowledge of the channel coefficients a s,1,1 , a s,1,2 , a s,2,1 , and a s,2,2 .

B. Channel Coefficient Method: Single PIM Source
Next, we consider the special case of having only a single air-PIM source in the vicinity of the transceiver, i.e., N P = 1. The cancellation scheme covered in this subsection has been initially introduced in our early work in [17], while it is also stated here for the sake of completeness and readability. To this end, the computational complexity of the full BF method introduced in Section III-A can be alleviated by summing the signals on the same frequencies, as indicated already by (1), where the relative phases and amplitudes of the superimposed signals are taken into account. Consequently, we define two aggregated signals v 1 where the complex multipliers are defined as q 1 = a 1,1,1 /a 1,1,2 and q 2 = a 1,2,1 /a 1,2,2 , which account for the relative phase and magnitude of the signals on the same frequency. From (8), (9), and (11), and when N P = 1, we can see that the coefficients q can then be obtained from the coefficients c i , i = 1, 2, . . . , 6, as follows: This indicates that we only require a third-order instantaneous model of the received signal to determine the q coefficients. For improved accuracy in the implementation of the algorithm, the coefficients q 1 and q 2 can be defined as arithmetic means of the shown expressions in (15) and (16). We can then construct the BF matrix [t] of (10) using the BFs expressed in Table I  Calculate β β β as per Eq. (12) 8: Calculate r n [t] as per Eq. (10) 10: end while BF coefficients of (11). Then, we can find the LS estimate of the required c coefficients as where θ θ θ [t] ∈ C ×6 contains the instantaneous, zero envelope lag (k = 0) third-order BFs of Table IV, stacked to a matrix, such that each BF occupies a single column. This means that the channel coefficient method requires the coefficient estimation to be carried out twice, rather than once as is the case with the full BF modeling. Despite this, the complexity reduction from this technique in the main path and learning processing is substantial. This will be quantified explicitly in Section III-D. Algorithm 2 compiles and summarizes the steps required by the single PIM source channel coefficient method. As in the full BF method, the recalculation of the LS-based variables β β β and c (and, thus, also q 1 and q 2 ) is not necessary on every iteration in a static environment.

C. Channel Coefficient Method: Multiple PIM Sources
Even though the method introduced in Section III-B offers a great reduction in computational complexity compared to the full BF method, its usability is restricted to cases with a single dominant air-PIM source. To this end, in this section, we introduce a more general cancellation scheme with reduced main path complexity compared to the full BF model, which can handle any number of PIM sources in the vicinity of the transceiver. The basic concept of this cancellation scheme is illustrated as a block diagram in Fig. 4.
From the definition of the coefficients c in (11), it stems that the channel coefficient method of Section III-B can only handle a single PIM source-as already noted-since the above coefficients describing the situation with multiple PIM sources contain summations. On the other hand, we can now build an estimate of the received PIM signal by utilizing an instantaneous third-order model, expressed aŝ (18) where the estimated PIM signalr [t] is dependent on all the channel coefficients a s,1,1 , a s,1,2 , a s,2,1 , and a s which measures how well our model matches the measured data. We can further define a loss function for the system as 2 2 , which we aim to minimize by finding optimal values for the channel coefficients a s,1,1 , a s,1,2 , a s,2,1 , and a s,2,2 , and for the lumped coefficients γ s,0,3,0 . For solving this, we employ the Gauss-Newton iterative method, as the problem is now nonlinear-in-parameters in its nature. To this end, by stacking the unknown parameters to vector α α α[t] ∈ C 5N P ×1 , the parameter update following the damped Gauss-Newton method, also known as the Levenberg-Marquardt algorithm, can be written as [40] (20) where J[t] ∈ C ×5N P is the Jacobian matrix of the function is the diagonal loading factor, and I is the appropriately sized identity matrix. The Jacobian matrix J[t] is defined by the partial derivatives of the functionr [t] w.r. t. all the parameters, namely, a s,1,1 ,  a s,1,2 , a s,2,1 , a s,2,2 , and γ s,0,1,0 , for s = 1, 2, . . . , N P . Due to the definition of the instantaneous estimate of the signal r [t], we can define the total Jacobian J[t] using submatrices J s [t] ∈ C ×5 , as  (10) and (12). The steps of the channel coefficient method with multiple PIM sources are shown in Algorithm 3. Again, the recalculation of vector β β β on every iteration is not necessary for a static environment. Now, the BF matrix [t] consists of N P independent copies of the BFs shown and summarized in Table I, i.e., there are N P third-order BFs, and 2N P or 3N P fifth-order BFs, depending on whether k is zero or nonzero. Therefore, the main path complexity of the proposed channel coefficient method with multiple PIM sources has slightly increased compared to the single PIM source method, but we still have a considerable complexity reduction compared to the full BF method when N P is moderate. This is addressed and quantified more concretely in Section III-D.

D. Computational Complexity Analysis
Here, we investigate the computational complexity of the aforementioned cancellation schemes in terms of FLOPs.  The computational complexities are shown in Table II, where M = M 1 + M 2 + 1 denotes the total memory, K = K 1 + K 2 refers to the total envelope lead/lag, and i c denotes the number of iterations it takes the Gauss-Newton algorithm in (20) to converge to a solution. The results are shown per output sample in the main path and per processed data block in the learning such that the complexity is normalized with the block size. To obtain the presented results, it is assumed that a complex multiplication takes eight FLOPs, and a complex addition takes two FLOPs. In addition, it is assumed that only the generation of the instantaneous BFs is necessary per iteration. Furthermore, QR decomposition is assumed to be used in the LS and Gauss-Newton algorithms. The complexity of the QR decomposition can be estimated to be 8n 2 FLOPs per block, where n is the number of columns in the matrices [t] or J[t] [41]. For improved accuracy, the q coefficients in (15) and (16) are calculated as the means of the shown definitions, excluding the terms with square roots, as these are resource-hungry operations. Fig. 5 shows the main path and learning complexities of the full BF method and channel coefficient method of Section III-C, with two different system parameterizations and with various numbers of PIM sources N P . Here, we have set the number of iterations for the Gauss-Newton algorithm to converge as i c = 10. It can be seen from Fig. 5 and Table II that both the main path and learning complexities of the full BF method are invariant w.r.t. the number of PIM sources N P . On the other hand, the complexity of both the main path and learning of the channel coefficient method increases with N P . Since the complexity of the generation of the instantaneous BFs increases with N P , the main path complexity of the channel coefficient method surpasses that of the full BF model when sufficiently many PIM sources are considered; however, this occurs only when N P > 10 for both considered parameterizations. Likewise, the complexity of the learning in the channel coefficient methods increases with N P . With the lower parameterization of M = 3 and K = 0, the full BF model complexity is surpassed when N P > 6, and with M = 5 and K = 2 when N P > 13. These results indicate the superiority of the proposed alternative methods in terms of computational complexity when a moderate number of significant PIM sources is considered in the system. It is noted here that, if many significant sources of PIM are located near the transceiver, relocating the transceiver or the corresponding antenna system is likely the optimal solution. Moreover, by setting N P = 1 in Table II, we can see that the complexity to determine the channel coefficients in the multiple PIM source method from Section III-C is 334i c FLOPs. Therefore, the proposed method of Section III-C is approximately as complex as the single PIM source method from Section III-B with a single PIM source per iteration, while also being a more general solution, capable of handling any number of PIM sources. However, since the multiple PIM source method requires multiple iterations to reach the optimal solution, there is a tradeoff between speed and complexity when only a single dominant PIM source is present in the system.
It is worth noting that changes in signal BW do not affect the relative computational complexities presented in this section of any of the models, as long as the parameterization of the algorithms is kept unaltered, since the provided values are normalized per sample and per processed block. Yet, if a high signal BW would require higher memory depth (M) in the modeling of the nonlinear response of the canceller than what is shown in Fig. 5, the relative complexity savings from the introduced channel coefficient methods would be even more substantial.

IV. RF MEASUREMENTS AND RESULTS
In order to validate the introduced PIM cancellation methods, a set of RF measurements was conducted using real-life cellular BS hardware. We will first describe the setup used for the measurements and then present the PIM cancellation results with one, two, and three air-PIM sources while considering also two transmit signal BWs for the CCs.

A. Measurement Setup
The RF measurement setup is shown in Fig. 6(a). All the measurements were conducted within an anechoic chamber, where a dual-TX/RX BS transmits and receives in FDD operation, transmitting four physical carrier waveforms, as described in Section III. The TX RF center frequencies are set to 5G NR band n3 frequencies of f RF,1 = 1819.0 MHz and f RF,2 = 1866.5 MHz, while the RX operates at the n3 receive band, at 1771.5 MHz, which is the lower third-order intermodulation frequency (2 f RF,1 − f RF,2 ) of the TX frequencies.
A PC outside the anechoic chamber is used to control the BS and collect data for postprocessing. In the measurements, 5G NR standard-compliant cyclic prefix orthogonal frequency division multiplexing (CP-OFDM) waveforms of the chosen carrier BWs of 5 and 20 MHz are utilized. The baseband sampling rates of the carriers and the RX signal are 7.68 and 30.72 MHz, respectively, for the two considered BWs. In the postprocessing, where the PIM cancellation takes place, the TX and RX signals are oversampled by a factor of 16 for the BF generation and synchronization tasks. The parameter estimation and the actual cancellation are carried out at the original baseband sampling rates. In the measurements, the obtained RX signal data vectors, captured over 10 ms, are always divided into two: the parameter estimation is carried out using the first half, and the cancellation is performed and evaluated in the second half.
The directive antennas of the BS with an antenna gain of around 16 dBi were aimed toward the tested PIM source(s) such that the TX signals hit the PIM source(s), undergo a nonlinear transformation, and are reflected toward the BS RX paths. Due to the high transmit power required to generate meaningful or observable levels of air-induced PIM, the TX signals are predistorted with a built-in DPD solution to avoid nonlinear distortion stemming from the PAs. The TX signals then fulfill the spectral emission limits set by the 5G NR standards [2]. It is also noted here that, since the intermodulation distortion considered in this work is generated outside the TX, there are no limits for such emissions in the current 5G NR standards. Specifically, there are additional emission requirements and limits defined in [2] for BS TXs in order to protect the BS's own and/or other colocated RXs. These requirements are, however, defined either for conducted measurements or for OTA measurements in a pure chamber and, hence, do not consider the air-PIM induced by the built environment around the BS antenna system in true deployment sites, which is the air-induced PIM problem studied in this article.
In the measurement arrangement, the PIM signals appear essentially synchronously from the PIM sources, as they are within the chamber. As a concrete example, considering a 5-MHz signal sampled at f s = 7.68 MHz, the distance difference between the PIM sources can be approximately 39 m, measured along the line-of-sight of the transceiver, for them to appear within the same sample, which is far greater than the dimensions of the anechoic chamber where the PIM sources lie. This distance is additionally so large that any potential sources close to the 39 m limit would probably not induce any significant PIM distortion since, in general, the larger the distance to the potential PIM source, the smaller the observed PIM on the RX band will be. Two different types of PIM sources were utilized in the measurements: standard off-the-shelf steel wool and purpose-made bow tie antennas, where the nonlinear response is achieved with a diode soldered to the middle of the antenna. The PIM sources are shown in Fig. 6(b) and (c). In the following, results of only RX path 1 are shown, for presentation simplicity, but the results from RX path 2 are similar.

B. Single PIM Source
First, let us examine the cancellation capabilities of the proposed PIM cancellation schemes with a single air-PIM source, i.e., N P = 1, where the PIM source in question is placed around 2.5 m from the BS unit. In Fig. 7, the achieved PIM cancellation is plotted with various TX carrier powers utilizing 5-MHz carriers with the parameterization of M 1 = M 2 = 2 and K 1 = K 2 = 1. The cancellers adopt the fifth-order modeling introduced in Section III; therefore, the full BF model employs 760 BFs in the matrix [t], while the channel coefficient methods employ only 65 BFs. In addition, the performance of a ninth-order canceller employing the channel coefficient method of Section III-C is shown. Furthermore, Fig. 8 shows the power spectral densities (PSDs) of the involved signals in the +31 dBm TX power case. It can be seen that all the fifth-order cancellers achieve the same level of cancellation with all the TX powers, indicating that the drastic reduction in complexity from the channel coefficient models does not compromise the cancellation capability. All of the fifth-order models mitigate the PIM all the way to the noise floor with TX powers up to +27 dBm. The residual PIM with fifth-order models then starts increasing with increased TX power, reaching a level of up to 2.3 dB at +31-dBm TX power, with the cancellation being 17.3 dB and the maximum being 19.6 dB in this case. These results are well in line with the results from our previous work in [17], where cancellation of up to 20 dB was evidenced in a similar setup.
The near 3 dB difference in cancellation compared to the one reported in our previous work [17] can be explained by the sensitivity of the setup, particularly the exact orientation and Power sweep with 5-MHz carriers, steel wool as PIM source, location of the steel wool PIM source w.r.t. the BS. Therefore, in addition to the fifth-order modeling, Figs. 7 and 8 demonstrate the higher cancellation fidelity obtained by the ninth-order model with the channel coefficient methods. The ninth-order model contains additional memory polynomial BFs (k = 0) up to polynomial order P 1 = 9, with the additional instantaneous BFs being shown in Table III. With the added BFs, the cancellation at +31 dBm TX power is now 19.3 dB, merely 0.6 dB from the noise floor. This added performance comes with the price of only adding seven new instantaneous BFs, which with the considered five taps of memory adds up to 35 new BFs, totaling then 100 BFs in the matrix [t]. In order to reach this level of accuracy in the model with the full BF modeling, BFs in the order of a thousand would be required, which again highlights the significant complexity reduction that the channel coefficient methods are able to provide.
The aforementioned measurements and results are repeated with the 20-MHz carriers in Fig. 9, while Fig. 10 plots the corresponding PSDs of the signals at +31 dBm TX power. In these measurements, a purpose-made diode antenna is used as the PIM source. The parameters utilized for these measurements are M 1 = M 2 = 2 and K 1 = K 2 = 3,  meaning that the full BF model incorporates now 1880 BFs in [t], while the channel coefficient models incorporate only 135 BFs. A similar conclusion can be drawn here as with the earlier 5-MHz case, namely, the fifth-order modeling of the channel coefficient models matches extremely well with the performance of the full BF-based fifth-order model, even with increased carrier BW. Again, with high transmit powers, the residual interference power starts increasing, being around 0.9 dB with +31 dBm TX power, where the achieved cancellation is 10.8 dB, while the maximum is 11.7 dB. Unlike in the 5-MHz case, there is no clear advantage in using the ninth-order model in this scenario, as the performance level is on par with the fifth-order models. This is likely due to the nonlinear characteristics of the utilized diode in the PIM source antenna.
In general, it can be observed from both Figs. 7 and 9 that, in order to avoid generating interference to the RX path, the TX power of each carrier needs to be backed off considerably-if no PIM cancellation is carried out. For example, let us allow an interference level of 2 dB above the noise floor in our system, where a single PIM source generates PIM distortion to the RX chain. It can be seen from both Figs. 7 and 9 that employ the cancellation schemes adhere to this limit, even when transmitting at +31 dBm per carrier. On the other hand, if no PIM cancellation is performed, it can be approximated that, in order to induce interference levels of only 2 dB above the noise floor, we would have to set the TX power of each carrier to +22 dBm in the 5-MHz case. This corresponds to a back-off of 9 dB per carrier compared to if a canceller is utilized. Likewise, we would need to back-off the TX power by approximately 7 dB per carrier in the 20-MHz case to induce an interference level of only 2 dB above the noise floor. Thus, by deploying the proposed PIM cancellation methods, the transceiver can increase the TX power by 9 and 7 dB per carrier in the 5-and 20-MHz cases, respectively, compared to the situation without PIM cancellation. This, in turn, increases both the TX power-efficiency and the downlink coverage of the network. Fig. 11 plots the convergence of the Gauss-Newton algorithm used in the single PM source measurements for both the 5-and 20-MHz cases, utilized by the channel coefficient model of Section III-C. The regularization coefficient λ[t] shown in the Gauss-Newton algorithm in (20) is initialized to 0.2 in the 5-MHz case and to 2 in the 20-MHz case. In order to avoid bias in the steady-state operation, the λ coefficient is updated using the following rule: which will effectively remove the biasing after sufficiently many iterations. It can be seen from Fig. 11 that both the 5-and 20-MHz cases converge to a steady-state operation in about eight iterations.

C. Multiple PIM Sources
In this section, we increase the number of PIM sources in the system to two and three, i.e., N P = 2 and N P = 3. In these measurements, we employ diode antennas as the air-PIM sources, which are placed 3 m from the BS unit, with a 1 m separation between the diode antennas. Fig. 12 plots the PSDs of the signals of interest when 5-MHz carriers are utilized at +31 dBm TX power with M 1 = M 2 = 2 and K 1 = K 2 = 1. Here, only fifthorder modeling is employed. First, it is evident that the channel coefficient method from Section III-B cannot handle the coexistence of multiple PIM sources in the system, as was expected based on the modeling of the system. Therefore, the achievable cancellation of the single PIM source channel coefficient model of Section III-B is merely 2 dB. On the other hand, the fifth-order full BF canceller is capable of mitigating the received PIM by some 16.6 dB, out of the maximum of 20.9 dB when referencing the noise floor. Similar performance is evidenced with the fifth-order channel coefficient model of Section III-C, where N P is now set to 2. With the chosen parameterization, the full BF model again employs 760 BFs in the matrix [t], while the multiple PIM source channel coefficient model now employs 130 BFs. Still, with the slightly elevated complexity compared to the single PIM source case, a substantial reduction in complexity is gained by employing the channel coefficient model, without penalty to the cancellation performance. Fig. 13 plots the PSDs of the signals of interest when 20-MHz carriers are transmitted at +31-dBm TX power. Here, the parameters are set to M 1 = M 2 = 2 and K 1 = K 2 = 3. As was the case with the 5-MHz carriers, the cancellation performance of the single PIM source channel coefficient model from Section III-B is poor, at only around 2 dB. Meanwhile, both the full BF model and multiple PIM source channel coefficient model from Section III-C achieve cancellation levels of around 13.9 dB, out of a maximum of around 15.2 dB, indicating that the multiple PIM source channel coefficient method is capable of providing the saving in complexity without loss of performance also with higher BWs. In this case, the full BF model incorporates again 1880 BFs, while the multiple PIM source channel coefficient model incorporates 270 BFs.
Next, three PIM antennas are utilized as PIM sources in Fig. 14 when 5-MHz carriers are transmitted. The cancellers follow the same parameterization as above with the 5-MHz carriers: M 1 = M 2 = 1 and K 1 = K 2 = 1. As before, the single PIM source channel coefficient-based canceller of Section III-B is clearly not capable of handling the PIM stemming from multiple coexisting sources. Meanwhile, the full BF canceller and the multiple PM source channel coefficient-based cancellers are able to suppress the PIM by around 16.5 and 16 dB, respectively, while the maximum cancellation is around  Finally, we observe and show the convergence of the canceled signal powers in the multiple PIM source scenarios. Fig. 15 illustrates the convergence of the canceled signal powers in the two PIM source cases with 5-and 20-MHz BWs. The parameter vector α α α[t] is again initialized to all ones, and the λ[t] coefficient is initialized as 0.1 in the 5-MHz case and to 0.15 in the 20-MHz case. The λ[t] coefficient is updated using the rule from (25). The initial drop in the error signal power occurs again in around eight iterations, but, to reach full convergence, around ten iterations are required in the 20-MHz case, while the 5-MHz case requires around 13 iterations to achieve convergence. Fig. 16 shows the convergence of the canceled signal power in the three PIM source cases with 5-MHz carriers and with two different initial conditions for λ[0], namely, λ[0] = 0.1 and λ[0] = 0.2. It can be seen that the initial convergence until  The convergence plots in Figs. 11, 15, and 16 show that, in general, it is difficult to estimate the exact convergence properties of the Gauss-Newton algorithm with different numbers of PIM sources. This is especially true when the number of PIM sources keeps increasing, as the behavior of the error surface w.r.t. the parameters in the vector α α α[t] becomes challenging to predict. This stems from the convergence being dependent on the exact position and orientation of the PIM sources w.r.t. the TX and RX antennas, which inevitably vary between the measurement cases. In Figs. 11, 15, and 16, it can be seen that, in the measured cases, the single PIM source case converges the fastest, which makes intuitive sense. However, the convergence in the two PIM source cases is slightly slower than in the three PIM source case when the same initial value for λ[0] is applied. Furthermore, the sensitivity to the initial conditions increases with the number of PIM sources, as evidenced by Fig 16. The slightly erratic performance when λ[0] = 0.2 could be amended by decreasing λ more slowly, yet that would delay the convergence even further. We omit further analysis of the convergence here, but it presents an interesting possible future work item.

D. Discussion and Future Considerations
Here, we will shortly discuss selected further details, and potential extensions and future work of the introduced PIM cancellation methods. An important detail of the proposed channel coefficient methods is the assumption that the channel between the transceiver and an individual air-PIM source can be approximated by a single complex tap. As mentioned before, this is a fair assumption in most circumstances and BS deployment scenarios when relatively low BW signals are utilized (such as in the measurements above) since the PIM source requires a line-of-sight channel to the transceiver in order to generate detectable and harmful levels of PIM.
Higher system BW could, in general, introduce frequency selectivity and therefore considerable memory effects to the channel model. However, at higher BWs, the PIM issue is typically less severe due to the reflected PIM power being spread over larger BW (lower PIM spectral density), therefore posing less of a threat to the RX sensitivity.
An additional valid note is that our analysis has assumed the PIM sources to be static, i.e., they do not move, and therefore, the channels between the transceiver and the PIM sources are approximately time-invariant. This is a valid assumption in many real-life situations, such as when the air-induced PIM is generated by nuts and/or bolts in the vicinity of the transceiver. However, in certain situations, the PIM may be generated by nonrigid objects that can move or vibrate, such as a steel fence vibrating in the wind, consequently altering the channel between the transceiver and the PIM source. The full BF method and the single PIM source channel coefficient method of Section III-B can inherently handle these changes within one processing block, as the adaptation or parameter estimation is essentially based on block LS in both models. The performance of the multiple PIM source channel coefficient method of Section III-C, on the other hand, is likely to suffer from moving PIM sources, as the solution can only be obtained iteratively. Though this is strictly speaking a matter of the available computational power.
In this work, we have considered the specific case of two carrier frequencies and two spatially multiplexed signals per carrier while canceling the intermodulation products falling on the frequency 2 f 1 − f 2 . As mentioned before, the generic modeling results in (5) can be used to determine the BFs with any number of signals per frequency, falling on any odd intermodulation frequency. Furthermore, in systems with more than two carrier frequencies, the pairwise modeling and cancellation of the PIM products, say, ψ a [t] and ψ b [t] falling on frequencies of the form 2 f a − f b , can be performed as presented in this article; however, there are then also additional third-order intermodulation frequencies. For example, in the case of three carrier frequencies, the instantaneous third-order BFs falling on the frequency f a + f b − f c take the form of ψ a [t]ψ b [t]ψ * c [t] [4]. In such cases, the presented channel coefficient methods are not directly valid in their current form but would require slight reformulations.
Finally, the PIM sources may also generate harmonic distortion, which is neglected in the system modeling in this work. In general, the harmonic products appear far from the passband in the frequency domain since the second harmonic frequencies are of the form f a + f b . Yet, there can be cases where, say, band n3 and n78 BSs are colocated, and in such cases, harmonic PIM generated from the band n3 BS can

V. CONCLUSION
In this article, we have investigated air-induced PIM and its efficient digital cancellation in a challenging FDD MIMO scenario. As a fundamental basis, extensive system modeling was provided, describing the total PIM signal at a given intermodulation frequency, when an arbitrary number of transceivers operate simultaneously and when there are an arbitrary number of PIM sources causing distortion. Stemming from such a system model, a corresponding digital cancellation scheme based on the full set of polynomial BFs, called the full BF method, was then described. To alleviate the high computational complexity of the full BF approach, an alternative method leveraging the air-PIM channel coefficients, called the channel coefficient method, was formulated. While the basic channel coefficient method allows only for handling the distortion stemming from a single dominant PIM source, we then proposed a novel Gauss-Newton-based iterative method, which can identify the channel coefficients of any number of concurrent PIM sources. The provided computational complexity analysis of all three methods reveals substantial complexity reduction when employing the proposed channel coefficient-based solutions compared to the full BF model-particularly when the number of significant PIM sources is moderate. All the methods were assessed and validated through extensive RF measurements at the 5G NR band n3, utilizing true BS hardware. The measurement results show that the high cancellation fidelity of the full BF model can be retained with much reduced complexity with the proposed channel coefficient-based methods. Our future work topics include addressing the potential time-varying nature of the air-PIM channels in realistic BS deployment sites while also generalizing the PIM cancellation solutions to handle harmonic PIM.

APPENDIX I BASIS FUNCTIONS OF THE FULL BF MODEL
For readers' convenience, the instantaneous BFs of the full BF model are given in Table IV up to the polynomial order of 5. These are stemming from (6).