CN106301417A - 一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法 - Google Patents

Abstract

本发明一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，属于抗干扰技术领域。本方法包括1在天线收发共用的通信系统中，双工器发射下行信号，再接收上行链路信号；2对1接收的上行链路信号进行稀疏分数阶傅立叶变换，估计幂级数模型的参数信息；3利用2估计得到的参数信息和下行信号重建PIM干扰信号；4判断3重建的PIM干扰信号的幅值是否超过阈值并决定跳至步骤5还是步骤6；5采用自适应滤波算法进一步进行PIM干扰自适应滤波，得出经自适应滤波后的PIM干扰信号；6修正上行链路信号，完成PIM干扰对消。本方法无需考虑参数估计收敛困难问题，利用稀疏分数阶傅立叶变换处理LFM信号，实现了低运算量的PIM干扰抑制。

Description

一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法

技术领域

本发明涉及一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，属于无线通信、卫星通信及地面通信系统的抗干扰技术领域。

背景技术

无源互调(Passive Intermodulation，简称PIM)是指由两个或两个以上频率成分经过无源器件非线性效应产生的一族产物，其频率为基波频率的线性组合。当线性组合产物落入通信接收机的通带内即形成干扰。无源互调广泛存在于通信系统中，是难以完全消除的现象，最早在上世纪70年代就被发现。随着通信需求的提高，高功率和高灵敏度已成为发展方向，通信所需的天线尺寸越来越大，受平台限制被迫采用收发共用技术。在双工器隔离度低的天线收发共用通信系统中，无源互调产物更易进入接收频带形成干扰，且无法通过频域滤波的方式将其滤除。以卫星转发器为例，下行大功率发射信号的无源互调产物会落入上行接收频带内，对通信系统造成影响。

为实现实时PIM干扰对消，爱立信公司、北京理工大学和华为公司均已申请专利，提出了基于数字信号处理的自适应PIM干扰对消方法，但是此类方法复杂度高、收敛速度慢，不能很好地适应PIM信号的时变特性。

其中，专利一为爱立信公司申请的发明专利，专利申请号为：US8855175B2，标题为“Low Complexity All-Digital PIM Compensator”，该专利提出了无需测试模式的自适应PIM干扰对消方法，该方法需要通过测试信号确定PIM参数个数，从而利用参数估计模型得到PIM模型参数，估计出PIM干扰信号。该方法复杂度高，资源消耗较大，且参数估计的收敛速度无法保证。

专利二为北京理工大学申请的发明专利，专利申请号为：201510547503.5，标题为“一种基于导频的无源互调干扰对消方法”，该专利采用导频时隙和数据传输时隙的双时隙机制对PIM干扰信号进行实时估计和抑制，在导频时隙，快速估计PIM模型参数，在数据传输时隙，利用已知发射导频信号及PIM模型参数，重建PIM干扰信号，并从接收信号中除去，完成PIM干扰对消。该方法可显著提高收敛性能，实现实时对消宽带信号引起的PIM干扰，但是对PIM模型参数的估计精度较差，对信道环境要求较高。

此外，北京理工大学、爱立信公司和华为公司申请的其他相关专利都没有解决参数估计收敛困难的问题。本专利致力于提出一种无需考虑收敛问题的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法。选用能覆盖整个发射频带的线性调频信号(LinearFrequency Modulation，简称LFM)作为导频信号，利用稀疏分数阶傅立叶变换，实现对落入接收频带内PIM的对消。

发明内容

本发明的目的在于解决现有无源互调干扰对消技术中参数估计收敛困难的问题，并进一步降低运算量，提高无源互调干扰对消效果和实时性，提出了一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法。

一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，简称本方法，包括以下步骤：

步骤一：在天线收发共用的通信系统中，双工器发射下行信号，再接收上行链路信号；

其中，下行信号包括导频信号和非导频信号，导频信号产生的PIM干扰信号s_PIM(t)可以表示为如下公式(1)：

$s_{PIM}\left(t\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{s}_p{\left(t\right)}^l---\left(1\right)$

其中，s_p(t)为导频信号，表示对每一项a_ls_p(t)^l进行求和，a_l表示幂级数模型阶次l的幅度值，即：导频信号产生的PIM干扰信号采用幂级数模型表示，且下行信号中的导频信号是线性调频信号，又称为LFM信号；

接收的上行链路信号包括上行信号、PIM干扰信号及噪声三部分，可以用如下公式(2)来表示：

s_R(t)＝s_UP(t)+I_PIM(t)*h(t)+n₀(t) (2)

其中，s_R(t)为接收的上行链路信号，s_UP(t)为上行信号，I_PIM(t)为下行信号产生的PIM干扰信号，包括导频信号产生的PIM干扰信号s_PIM(t)及非导频信号产生的PIM干扰信号，h(t)表示PIM干扰信号由双工器发射到接收的耦合信道冲激响应，I_PIM(t)*h(t)表示双工器接收到的PIM干扰信号；其中，*表示卷积运算，n₀(t)为噪声；

步骤二：对步骤一接收的上行链路信号进行稀疏分数阶傅立叶变换，得到不同分数阶傅立叶旋转角度下的幅度值，对幅度值取峰值即得出导频信号产生的PIM干扰信号在分数域峰值的位置和大小，进而估计导频信号产生的PIM干扰信号携带的参数信息；

其中，导频信号产生的PIM干扰信号S_PIM(t)携带参数信息：幂级数模型阶次l和幂级数模型阶次l的幅度值a_l，记为和

步骤二，具体为：

步骤2.1：对导频信号产生的PIM干扰信号进行量纲归一化，输出量纲归一化后的导频信号产生的PIM干扰信号，简称归一化PIM信号；

首先，离散化表示导频信号产生的PIM干扰信号s_PIM(t)，通过如下公式(3)：

$s_{PIM}\left(n\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{s}_p{\left(n\right)}^l=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{e}^{jl(2{\mathrm{\πf}}_{0}n/f_s+{\mathrm{\π\βn}}^2/{f_s}^2)}---\left(3\right)$

其中，s_PIM(n)为导频信号产生的PIM干扰信号的离散化表示，n代表离散化变量，s_p(n)为导频信号s_p(t)的离散化表示；j代表虚数，f₀表示LFM信号的起始频率，β表示LFM信号的调频率，采样间隔为1/f_s，导频信号产生的PIM干扰信号的频域和时域范围分别为[-f_s/2,f_s/2]和[-T/2,T/2]；

然后，引入一个具有时间量纲的归一化因子定义量纲归一化坐标p＝t/S，q＝f·S，借助量纲归一化坐标，将频域和时域范围两个区间：[-f_s/2,f_s/2]和[-T/2,T/2]都归一化为[-Δx/2,Δx/2]，

归一化PIM信号可用如下公式(4)表示：

${s_{PIM}}^{\′}\left(n\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{e}^{jl(2{\mathrm{\πf}}_{0}n/\mathrm{\Δ}x+{\mathrm{\π\β}}^{\′}{n}^2/{\mathrm{\Δx}}^2)}---\left(4\right)$

其中，s_PIM′(n)表示归一化PIM信号，为归一化PIM信号的调频率，1/Δx为归一化PIM信号的采样间隔；

步骤2.2：对步骤2.1输出的归一化PIM信号进行稀疏分数阶傅立叶变换，具体采用Pei采样型离散算法进行稀疏分数阶傅立叶变换，得到稀疏分数阶傅立叶变换结果；

采用Pei采样型离散算法进行稀疏分数阶傅立叶变换的原因在于：由于导频信号产生的PIM干扰信号为多分量LFM信号，此多分量LFM信号在分数域具有稀疏性，可以用稀疏傅立叶变换(Sparse Fourier Transform，简称SFT)代替Pei采样型离散算法中的傅立叶变换，达到降低运算复杂度的目的；包括如下步骤：

步骤Pei.1：将步骤2.1输出的归一化PIM信号与chirp1信号相乘；

相乘后得到序列x(n)用如下公式(5)表示：

$x\left(n\right)=e^{j\frac{{\mathrm{cot\αn}}^2}{2{\mathrm{\Δx}}^2}}{s_{PIM}}^{\′}\left(n\right)---\left(5\right)$

其中，为chirp1信号，cot为反正切函数，α为稀疏分数阶傅立叶变换的旋转角度；

步骤Pei.2：将步骤Pei.1输出的归一化PIM信号与chirp1信号相乘的结果进行稀疏傅立叶变换，得到归一化PIM信号与chirp1信号相乘结果的稀疏傅立叶变换的估计值，简称稀疏傅立叶变换的估计值，具体为：

步骤Pei.2.A：对x(n)进行重新排序，得到重排后的序列，即定义重排方式为P_σ,τ，重排后序列用如下公式(6)表示：

(P_σ，τx)_i＝x_(σi+τ)\N (6)

其中，(P_σ，τx)_i为重排后序列P_σ,τx中的第i个值，σ是从[1,N]中随机选取的奇数，τ是从[1,N]中随机选取的整数，N为采样点数且N＝(Δx)²，(σi+τ)\N表示对索引值σi+τ的模N操作；

步骤Pei.2.B：将步骤Pei.2.A输出的重排后序列与滤波器相乘，得到相乘后序列，具体为：

所述滤波器采用低通滤波器，此滤波器的传输函数为G(∈,∈′,δ,ω)，此滤波器的傅立叶变换需要满足如下公式(7)：

${\hat{G}}_i=\left\{\begin{array}{cc}\[1-\mathrm{\δ},1+\mathrm{\δ}\]& i\∈\[-{\mathrm{\ϵ}}^{\′}N,{\mathrm{\ϵ}}^{\′}N\]\\ \[0,\mathrm{\δ}\]& i\∈\[-\mathrm{\ϵ}N,\mathrm{\ϵ}N\]\end{array}\right.---\left(7\right)$

其中，表示滤波器傅立叶变换的第i个值，∈′和∈分别表示通带截止频率和阻带截止频率；通带和阻带的最大衰减相同，用δ表示；ω表示滤波器的长度；

相乘后序列可用如下公式(8)表示：

y_i＝G_i×(P_σ,τx)_i (8)

其中，y_i为相乘后序列的第i个值，G_i为滤波器传输函数的第i个值；

步骤Pei.2.C：对步骤Pei.2.B输出的相乘后序列进行时域混叠，得到混叠后序列，并计算此混叠后序列的傅立叶变换；

混叠后序列可用如下公式(9)表示：

$z_i=\underset{jj=0}{\overset{N/B-1}{\Σ}}{y}_{i+Bjj}---\left(9\right)$

其中，z_i为混叠后序列z第i个值，jj表示求和变量，范围是0到N/B-1，B为混叠长度，参考值是(k为稀疏度)；

记混叠后序列的傅立叶变换结果为

定义散列函数用如下公式(10)表示：

h_σ(i)＝round(σ×i×B/N) (10)

其中，h_σ(i)为散列函数的第i个值，round(σ×i×B/N)表示对σ×i×B/N进行四舍五入运算；

定义偏移函数用如下公式(11)表示：

o_σ(i)＝σ×i-h_σ(i)×(N/B) (11)

其中，o_σ(i)为偏移函数的第i个值；

步骤Pei.2.D：取出中幅度最大的前d×k个幅值及其坐标，其中，d为稀疏度增益，参考值是将取出的坐标保存在集合J中；再通过散列函数h_σ(i)∈J折算出坐标i∈[1，N]，将折算出的坐标保存在集合I中；

步骤Pei.2.E：利用步骤Pei.2.D折算出的坐标i计算原始序列x的傅立叶变换的估计值，可以表示为如下公式(12)：

${{\hat{x}}_i}^{\′}={\hat{z}}_{h_{\mathrm{\σ}}\left(i\right)}{\mathrm{\ω}}^{\mathrm{\τ}\×i}/{\hat{G}}_{o_{\mathrm{\σ}}\left(i\right)}---\left(12\right)$

其中，表示原始序列傅立叶变换估计值中的第i个值，为中第h_σ(i)个值，h_σ(i)取集合J中保存的坐标，ω表示为步骤Pei.2.B中滤波器的傅立叶变换的第o_σ(i)个值；

步骤Pei.2.F：记步骤Pei.2.D输出的I为I₁，重复步骤Pei.2.A到步骤Pei.2.E L-1次，对每次循环步骤Pei.2.D得到的集合I顺序编号，第R次循环步骤Pei.2.D得到的集合I记为I_R+1，即：步骤Pei.2.A到步骤Pei.2.E共进行L次运算，第r次运算步骤Pei.2.D得到的集合I为I_r；

其中，L表示运算次数，其参考值是L＝log₂N；R为循环变量，变化范围为1到L-1，r为运算次数变量，变换范围为1到L；

步骤Pei.2.G：统计L次运算步骤Pei.2.D得到的集合I，记为I₀＝I₁∪I₂∪I₃∪...∪I_L，统计每个折算出来的坐标i∈I₀出现的次数，保存在集合s_i中，即s_i＝|{r|i∈I_r}|，取出现次数超过L/2次的坐标，保存在集合I′中，即I′＝{i∈I₀|s_i≥L/2}；

步骤Pei.2.H：对集合I′中每个坐标i，从每次运算中取出在步骤Pei.2.E计算得到的原始序列x的傅立叶变换的估计值，记为并取中位值作为最后的估计值，得到的估计值用如下公式(13)表示：

${\hat{x}}_i^{\′\′}=median\left(\right\{{\hat{x}}_i^r|i\∈I^{\′}\})---\left(13\right)$

其中，表示稀疏傅立叶变换的估计值，表示在i∈I′条件下对取中位值；

步骤Pei.3：将步骤Pei.2计算得到的稀疏傅立叶变换的估计值与chirp2信号相乘，得到稀疏分数阶傅立叶变换结果；

得到稀疏分数阶傅立叶变换的结果用如下公式(14)表示：

$\widehat{F_{\mathrm{\α}}}\left(m\right)={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}\sqrt{\frac{(\sin \mathrm{\α}-j\cos \mathrm{\α})\mathrm{sgn}\left(\sin \mathrm{\α}\right)}{M}}={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}{A}_{\mathrm{\α}}---\left(14\right)$

其中，表示稀疏分数阶傅立叶变换的结果，表示chirp2信号，m为离散化变量，sin为正弦函数，cos为余弦函数，sgn为符号函数，A_α表示M为稀疏分数阶傅立叶变换后信号的采样点数，需要满足M≥N，Δu为稀疏分数阶傅立叶变换后信号的采样间隔，需要满足如下公式(15)：

$\frac{\mathrm{\Δ}u}{\mathrm{\Δ}x}=\frac{2\mathrm{\π}\left|sin\mathrm{\α}\right|}{N}---\left(15\right)$

步骤2.3：寻找步骤2.2得到的稀疏分数阶傅立叶变换结果的峰值位置和大小，估计导频信号产生的PIM干扰信号携带的参数信息；

其中，稀疏分数阶傅立叶变换的结果即步骤Pei.3的输出导频信号产生的PIM干扰信号携带的参数信息即幂级数模型阶次l及其幅度值a_l的估计值：和按照式(16)和式(17)估计：

$\hat{l}=-\frac{{\mathrm{cot\α}}^{\′}}{2{\mathrm{\πk\β}}^{\′}}=-\frac{{\mathrm{cot\α}}^{\′}{f_s}^2}{2{\mathrm{\π\β\Δx}}^2}---\left(16\right)$

$a_{\hat{l}}=\frac{max\left|\widehat{F_{\mathrm{\α}}}\left(m\right)\right|}{(2N+1)\left|A_{{\mathrm{\α}}^{\′}}\right|}---\left(17\right)$

其中，表示幂级数模型阶次l的估计值，表示幂级数模型幅度值a_l的估计值，α′为的峰值所对应的稀疏分数阶傅立叶变换的旋转角度，表示取的最大值，表示对进行绝对值运算，|A_α′|表示对A_α′进行绝对值运算，A_α′为的峰值所对应的A_α的值；

步骤三：利用步骤二估计得到的导频信号产生的PIM干扰信号携带的参数信息和下行信号重建PIM干扰信号；

其中，重建的PIM干扰信号用如下公式(18)表示：

${\hat{I}}_{PIM}\left(t\right)=\underset{\hat{l}}{\Σ}{a}_{\hat{l}}{s}_D{\left(t\right)}^{\hat{l}}---\left(18\right)$

其中，为重建的PIM干扰信号，s_D(t)为下行信号，包括导频信号和非导频信号；

步骤四：判断步骤三重建的PIM干扰信号的幅值是否超过阈值，并决定跳至步骤五还是步骤六，具体为：

步骤4.1：若超过阈值，跳到步骤五采用横向滤波器结构进一步进行PIM干扰自适应滤波；

步骤4.2：若不超过阈值，直接输出步骤三重建的PIM干扰信号并跳至步骤六；

其中，阈值的选取由PIM干扰对消后的误码率来决定；

步骤五：采用自适应滤波算法估计PIM干扰信号由双工器发射到接收的耦合信道冲激响应，得出经自适应滤波算法处理后的PIM干扰信号；

其中，自适应滤波算法可为最小均方算法(LMS)，采用M-1个延迟单元的级联横向滤波器实现，M表示横向滤波器的长度，通过不断调整横向滤波器的系数使得均方误差最小；经过自适应滤波算法处理后的PIM干扰信号表示为

步骤六：修正上行链路信号，完成PIM干扰对消；具体为：

步骤6.1：若进行了自适应滤波算法处理，则用步骤五得到的自适应滤波算法处理后的PIM干扰信号修正上行链路信号，修正后的上行链路信号用如下公式(19)表示：

${s_R}^{\′}\left(t\right)=s_R\left(t\right)-{{\hat{I}}_{PIM}}^{\′}\left(t\right)---\left(19\right)$

其中，s_R′(t)表示修正后的上行链路信号；

步骤6.2：若未进行自适应滤波算法处理，则用步骤三得到的重建PIM干扰信号修正上行链路信号，修正后的上行链路信号用如下公式(20)表示：

${s_R}^{\′\′}\left(t\right)=s_R\left(t\right)-{\hat{I}}_{PIM}\left(t\right)---\left(20\right)$

其中，s_R"(t)表示修正后的上行链路信号；

至此，从步骤一到步骤六，完成了一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法。

有益效果

本发明提出的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，与现有无源互调干扰对消方法相比，具有如下有益效果：

1.本发明所提方法计算复杂度低，资源代价小，易于实现；

2.本发明所提方法通过数字实现，性能稳定，系统适应性强；

3.本发明所提方法无需考虑参数估计收敛困难的问题，实时性较好；

4.本发明所提方法首次将稀疏信号处理理论用于PIM干扰对消，利用稀疏分数阶傅立叶变换的处理多分量LFM信号，实现了低运算量的PIM干扰抑制；

5.本发明所提方法消除无源互调干扰更有效，能够实现更高的干扰抑制增益。

附图说明

图1为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例1的系统框图；

图2为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”中本方法及实施例1的实现流程图；

图3为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例2中参数估计方法的实现流程图；

图4为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例3中自适应横向滤波器的结构示意图；

图5为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例3中采用自适应滤波算法估计耦合信道冲激响应的结构示意图；

图6为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例4中对导频信号产生的的PIM干扰信号进行稀疏分数阶傅立叶变换的仿真结果；

其中，x轴为旋转角度，y轴为采样点，z轴为幅度；

图7为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例4中PIM干扰对消前后时域对比图的仿真结果；

其中，x轴为时间，y轴为幅度；

图8为本发明“一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法”实施例4中在上行信号功率与PIM干扰信号功率比(信干比)为0dB条件下对PIM干扰信号进行抑制，得到误码率随Eb/NO变化的仿真结果；

其中，x轴为Eb/NO，y轴为误码率。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，下面结合附图及实施例，对本发明提出的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法做进一步详细的描述。应当理解，此处所描述的具体实施例仅用于解释本发明，并不用于限定本发明。

实施例1

本实施例阐述了本发明所提出的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法所依托的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消系统，简称本系统。系统框图如图1所示，这也是本发明实现本方法，达到PIM干扰对消的原理框图。

从图1可以看出，本系统由双工器、PIM参数估计模块、PIM干扰重建模块、PIM干扰检测模块及PIM干扰自适应对消模块组成；

双工器完成下行信号的发射和上行链路信号的接收，线性调频信号作为下行信号的导频信号通过双工器发射，下行信号包括导频部分和非导频部分，下行信号产生的PIM干扰信号采用幂级数模型表示，包括导频信号产生的PIM干扰信号及非导频信号产生的PIM干扰信号，上行链路信号包括上行信号、PIM干扰信号和噪声三部分，其中，导频信号产生的PIM干扰信号用来估计幂级数模型的参数；

PIM参数估计模块通过对接收到的上行链路信号进行处理，由此估计幂级数模型的参数，PIM参数估计模块对上行链路信号进行稀疏分数阶傅立叶变换，得到不同分数阶傅立叶旋转角度下的幅度值，其中，导频信号产生的PIM干扰信号携带参数信息，根据导频信号产生的PIM干扰信号在分数域峰值的位置和大小，估计导频信号产生的PIM干扰信号携带的参数信息，包括幂级数模型阶次和幂级数模型阶次的幅度值，最后将幂级数模型参数的估计值送至PIM干扰重建模块；

PIM干扰重建模块利用已知的下行信号和PIM参数估计模块得到的幂级数模型参数估计值重建PIM干扰信号，具体按照幂级数模型进行重建，最后将重建的PIM干扰信号送至PIM干扰检测模块和PIM自适应对消模块；

PIM干扰检测模块判断重建的PIM干扰信号的幅值是否超过阈值即是否需要进行PIM干扰自适应对消，若超过阈值则需要进行PIM干扰自适应对消，反之则不需要，判断的结果送至PIM干扰自适应对消模块；其中，阈值的选取由PIM干扰对消后的误码率来决定；

PIM干扰自适应对消模块根据PIM干扰检测模块判断的结果进行相应的处理，若重建的PIM干扰信号的幅值未超过阈值则直接用重建的PIM干扰信号修正上行链路信号，若重建的PIM干扰信号的幅值超过阈值则需要进一步消除接收信号中的PIM干扰信号，采用自适应滤波算法估计耦合信道的冲激响应，用自适应滤波后的PIM干扰信号修正上行链路信号。

图2为发明内容中本方法的实现流程图及本实施例实现PIM干扰对消的流程图，从图中可以看出，本方法包括以下步骤：

步骤A：双工器发射下行信号，再接收上行链路信号，其中，下行信号的导频信号是线性调频信号；

步骤B：对接收的上行链路信号进行稀疏分数阶傅立叶变换，估计导频信号产生的PIM干扰信号携带的参数信息；

步骤C：利用估计得到的导频信号产生的PIM干扰信号携带的参数信息和下行信号重建PIM干扰信号；

步骤D：检测重建的PIM干扰信号的幅值是否超过阈值，并决定跳至步骤E还是步骤F，具体为：

步骤D.1：若超过阈值，跳到步骤E；

步骤D.2：若不超过阈值，跳至步骤F；

步骤E：采用横向滤波器结构进一步进行PIM干扰自适应滤波，输出自适应滤波后的PIM干扰信号；

步骤F：修正上行链路信号，完成PIM干扰对消，具体为：

步骤F.1：若进行了自适应滤波，则用步骤E输出的自适应滤波后的PIM干扰信号修正上行链路信号；

步骤F.2：若未进行自适应滤波，则用步骤C得到的重建PIM干扰信号修正上行链路信号；

至此，从步骤A到步骤F完成了一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法。

本实施例描述了PIM干扰对消的原理和实现方法，可以看出，该方法通过数字实现，性能稳定，系统适应性强，该方法计算复杂度低，资源代价小，易于实现，该方法无需考虑参数估计收敛困难的问题，实时性较好，该方法利用稀疏分数阶傅立叶变换实现了低运算量的PIM干扰抑制。

实施例2

本实施例对本发明步骤二中导频信号产生的PIM干扰信号携带的参数信息的估计过程进行详细的描述，由于大功率下行信号是完全已知的，利用导频信号估计出幂级数模型的参数即可重建PIM干扰信号，实现PIM干扰对消。请参考图3，本发明实施例提供一种幂级数模型参数的估计方法，可包括：

步骤1)：对导频信号产生的PIM干扰信号进行量纲归一化，输出量纲归一化后的导频信号产生的PIM干扰信号，简称归一化PIM信号；

具体步骤与发明内容中步骤2.1相同，本实施例中，幂级数模型的阶次为7.3和8.9，相应阶次的幅度值分别为0.8和0.7，导频信号产生的PIM干扰信号为：

$s_{PIM}\left(n\right)=\underset{i=1}{\overset{2}{\Σ}}{a}_i{e}^{{\mathrm{jl}}_i(2{\mathrm{\πf}}_{0}n/f_s+{\mathrm{\π\βn}}^2/{f_s}^2)}$

其中，l₁＝7.3，l₂＝8.9，a₁＝0.8，a₂＝0.7，初始频率f₀＝0.5MHz，调频率β＝1.72×10⁴MHz，采样频率f_s＝112MHz，量纲归一化后，信号采样点数为N＝(Δx)²＝8192，归一化PIM信号表示为：

其中，为归一化PIM信号的调频率，1/Δx为归一化PIM信号的采样间隔；步骤2)：对步骤1)输出的归一化PIM信号进行稀疏分数阶傅立叶变换，具体采用Pei采样型离散算法进行稀疏分数阶傅立叶变换，得到稀疏分数阶傅立叶变换结果；傅立叶具体步骤与发明内容中步骤2.2相同，包括如下步骤：

步骤2.1)：将步骤1)输出的归一化PIM信号与chirp1信号相乘；；

相乘后得到序列x(n)表示为：

$x\left(n\right)=e^{\frac{{\mathrm{jcot\αn}}^2}{2{\mathrm{\Δx}}^2}}{s_{PIM}}^{\′}\left(n\right)$

其中，α为稀疏分数阶傅立叶变换的旋转角度；具体到本实施例，α取0到π，精度为0.0005π；

步骤2.2)：将步骤2.1)输出的归一化PIM信号与chirp1信号相乘的结果进行稀疏傅立叶变换，得到归一化PIM信号与chirp1信号相乘结果的稀疏傅立叶变换的估计值，简称稀疏傅立叶变换的估计值；傅立叶

具体到本实施例，稀疏傅立叶变换按照以下步骤实现：

步骤2.2.1)：对x(n)进行重新排序，得到重排后的序列，即定义重排方式为P_σ,τ，重排后序列用如下公式表示为(P_σ,τx)_i＝x_(σi+τ)\N；

其中，(P_σ,τx)_i为重排后序列P_σ,τx中的第i个值，具体到本实施例，σ是从[1,8192]中随机选取的奇数，τ是从[1,8192]中随机选取的整数，(σi+τ)\N表示对索引值σi+τ的模N＝8192操作；

步骤2.2.2)：将步骤2.2.1)输出的重排后序列与滤波器相乘，得到相乘后序列，具体到本实施例，所述滤波器采用选用多尔夫-切比雪夫滤波器，此滤波器的传输函数为G(∈,∈′,δ,ω)，其中，∈′和∈分别表示通带截止频率和阻带截止频率；通带和阻带的最大衰减相同，用δ表示；ω表示滤波器的长度；

相乘后序列表示为：y_i＝G_i×(P_σ,τx)_i；其中，y_i为相乘后序列的第i个值，G_i为滤波器传输函数的第i个值；

步骤2.2.3)：对步骤2.2.2)输出的相乘后序列进行时域混叠，得到混叠后序列，并计算此混叠后序列的傅立叶变换；

混叠后序列可用如下公式表示为其中，z_i为混叠后序列z第i个值，jj表示求和变量，本实施例中，范围是0到31，B为混叠长度，本实施例中，B＝256；

记混叠后序列的傅立叶变换结果为

定义散列函数：其中，h_σ(i)为散列函数的第i个值，round(σ×i×B/N)表示对σ×i×B/N进行四舍五入运算；

定义偏移函数：o_σ(i)＝σ×i-h_σ(i)×(N/B)；其中，o_σ(i)为偏移函数的第i个值；

步骤2.2.4)：取出中幅度最大的前2个幅值及其坐标，将取出的坐标保存在集合J中；本实施例中，稀疏度增益d＝2，稀疏度k＝1；再通过散列函数h_σ(i)∈J折算出坐标i∈[1,8192]，将折算出的坐标保存在集合I中；

步骤2.2.5)：利用步骤2.2.4)折算出的坐标i计算原始序列x的傅立叶变换的估计值，可以表示为：其中，表示原始序列傅立叶变换估计值中的第i个值，为中第h_σ(i)个值，h_σ(i)取集合J中保存的坐标，ω表示为步骤2.2.2)中滤波器的傅立叶变换的第o_σ(i)个值；

步骤2.2.6)记步骤2.2.4)输出的I为I₁，重复步骤2.2.1)到步骤2.2.5)12次，对每次循环步骤2.2.4)得到的集合I顺序编号，第R次循环步骤2.2.4)得到的集合I记为I_R+1，即：步骤2.2.1)到步骤2.2.5)共进行L＝13次运算，第r次运算步骤2.2.4)得到的集合I为I_r；

其中，R为循环变量，变化范围为1到12，r为运算次数变量，变换范围为1到13；

步骤2.2.7)统计13次运算步骤2.2.4)得到的集合I，记为I₀＝I₁∪I₂∪I₃∪...∪I₁₃，统计每个折算出来的坐标i∈I₀出现的次数，保存在集合s_i中，即s_i＝|{r|i∈I_r}|，取出现次数超过6次的坐标，保存在集合I′中，即I′＝{i∈I₀|s_i≥6}；

步骤2.2.8)对集合I′中每个坐标i，从每次运算中取出在步骤2.2.5)计算得到的原始序列x的傅立叶变换的估计值，记为并取中位值作为最后的估计值，得到的估计值用表示为其中，表示稀疏傅立叶变换的估计值，表示在i∈I′条件下对取中位值；

步骤2.3)将步骤2.2)计算得到的估计值与chirp2信号相乘；

最后，得到稀疏分数阶傅立叶变换的结果为：

$\widehat{F_{\mathrm{\α}}}\left(m\right)={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}\sqrt{\frac{(sin\mathrm{\α}-jcos\mathrm{\α})\mathrm{sgn}\left(sin\mathrm{\α}\right)}{M}}={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}{A}_{\mathrm{\α}}$

其中，表示稀疏分数阶傅立叶变换的结果，表示chirp2信号，m为离散化变量，sin为正弦函数，cos为余弦函数，sgn为符号函数，A_α表示M为稀疏分数阶傅立叶变换后信号的采样点数，需要满足M≥8192，Δu为稀疏分数阶傅立叶变换后信号的采样间隔，需要满足：

步骤3)：搜索峰值的位置和大小，估计PIM行为模型的参数；

用LFM信号作为导频信号并且选择幂级数模型时，导频部分的PIM干扰信号为多分量LFM信号，该信号在分数域会出现几个峰值，设置门限值P_th，假设当峰值大于门限值时，相应阶次的PIM干扰信号会影响通信，峰值的搜索过程按照以下步骤实现：

步骤3.1)：在二维分数域搜索最大的幅度值

步骤3.2)：判断搜索到的最大的幅度值是否超过门限值，即

步骤3.2.1)：若超过门限值，按照式(16)和式(17)通过峰值所对应的稀疏分数分数阶傅立叶变换的旋转角度和幅度值分别计算出幂级数模型参数的估计值，包括幂级数模型阶次及其次a_l的估计值，然后通过窄带滤波降低最大峰值，滤波之后信号表示为：其中，M(m)为窄带滤波器；最后，跳到步骤3.1)搜索下一个峰值；

步骤3.2.2)：若不超过门限值，则完成搜索过程；

具体到本实施例，导频信号产生的PIM干扰信号在分数域会出现两个峰值，需要进行两次搜索，可以分别计算得到PIM模型参数的四个估计值和

本实施例对本发明导频信号产生的PIM干扰信号携带的参数信息的估计过程进行了详细的说明，本方法首次将稀疏信号处理的理论用于PIM干扰抑制，利用稀疏分数阶傅立叶变换处理多分量LFM信号实现PIM干扰信号的对消，实现了低运算量的PIM干扰抑制，该方法通过数字实现，性能稳定，系统适应性强。

实施例3

利用导频信号产生的PIM干扰信号估计出的幂级数模型的参数和下行信号即可重建PIM干扰信号，若重建的PIM干扰信号超过设定的阈值，需要采用横向滤波器结构进一步进行PIM干扰自适应滤波，本实施例在实施例2的基础上对本发明步骤五中采用自适应滤波算法估计PIM干扰信号由双工器发射到接收的耦合信道冲激响应的过程进行详细的描述。

为进一步消除接收信号中的PIM干扰信号，采用自适应滤波算法估计耦合信道冲激响应h(t)，估计过程采用数字信号处理。采用横向滤波器结构，如图4所示，存储用M-1个延迟单元的级联表示，[w₀,w₁,...,w_M-1]表示自适应滤波器的系数，u(n)表示自适应滤波器的输入信号，d(n)表示自适应滤波器的期望信号，e(n)表示误差信号。

图5为采用自适应滤波算法估计耦合信道冲激响应的结构示意图。从图中可以看出，经过自适应滤波器得到的无源互调估计量为：

${{\hat{I}}_{PIM}}^{\′}\left(n\right)=\underset{k=0}{\overset{M-1}{\Σ}}{w}_k^{*}{\hat{I}}_{PIM}(n-k)$

由此得系统输出误差为：

$e\left(n\right)=s_R\left(n\right)-{{\hat{I}}_{PIM}}^{\′}\left(n\right)$

采用LMS自适应滤波器即不断调整横向滤波器系数使得均方误差J最小。

$J=E\left\{e\left(n\right)e{\left(n\right)}^{*}\right\}=E\left\{\left(s_R\right(n)-{{\hat{I}}_{PIM}}^{\′}(n\left)\right){\left(s_R\right(n)-{{\hat{I}}_{PIM}}^{\′}(n\left)\right)}^{*}\right\}$

易知，若求J最小值，需满足▽J＝0；

▽J(n)＝-2P(n)+2R(n)w(n)

其中，

$R\left(n\right)={\hat{I}}_{PIM}\left(n\right){{\hat{I}}_{PIM}}^H\left(n\right)$

$P\left(n\right)={\hat{I}}_{PIM}\left(n\right){s_R}^{*}\left(n\right)$

而耦合信道是不断变化的，需根据估计误差自动调整滤波器参数，即抽头权向量满足其中，w(n)表示n时刻自适应滤波器系数[w₀,w₁,...,w_M-1]组成的向量，μ是步长参数，取正常数。

实施例4

本实施例按照实施例1所述的对消方法、实施例2所述的参数估计方法及实施例3所述的自适应滤波方法进行仿真，具体阐述了对QPSK调制信号执行本发明的对消结果，仿真结果如图6、图7和图8所示。

其中，图6是对导频信号产生的PIM干扰信号进行稀疏分数阶傅立叶变换的结果。x轴为旋转角度，y轴为采样点，z轴为幅度，从图中可以看出，导频信号产生的PIM干扰信号在分数域具有稀疏性，而且稀疏分数阶傅立叶变换在处理多分量LFM信号时不存在交叉性干扰的问题，所以我们可以根据导频信号产生的PIM干扰信号的稀疏分数阶傅立叶变换结果的峰值位置和大小，估计导频信号产生的PIM干扰信号携带的参数信息；

图7是在本实施例环境下PIM干扰对消前后的时域对比图，其中，x轴为时间，y轴为幅度，圆形线表示无PIM干扰的信号，星形线表示有PIM干扰的信号，点形线表示PIM干扰对消后恢复的信号，从图中可以看出，经过PIM干扰对消后恢复的信号几乎与无PIM干扰的信号重合，该算法能够有效消除PIM干扰信号。

图8是在上行信号功率与PIM干扰信号功率比(信干比)为0dB条件下对PIM干扰信号进行抑制，得到误码率随Eb/NO变化的仿真结果，其中，x轴为Eb/NO，单位为dB,y轴为误码率，叉形线表示无PIM干扰的理论误码率曲线，十字形线表示仿真得到的有PIM干扰的误码率曲线，方形线表示仿真得到的无PIM干扰的误码率曲线，圆形线表示PIM干扰对消后的误码率曲线，从图中可以看出，PIM干扰对消后的误码率优于有PIM干扰的误码率，并且，PIM干扰对消后的误码率曲线与无干扰条件下的理论曲线逼近度高，PIM干扰对消后的误码率曲线与无PIM干扰的误码率曲线几乎完全重合，可见该算法能有效抑制无源互调干扰，实现较高的干扰抑制增益，降低误码率。

以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (7)

1.一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，简称本方法，所依托的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消系统，包括双工器、PIM参数估计模块、PIM干扰重建模块、PIM干扰检测模块及PIM干扰自适应对消模块；本方法的特征在于：包括以下步骤：

步骤一：在天线收发共用的通信系统中，双工器发射下行信号，再接收上行链路信号；

步骤二：对步骤一接收的上行链路信号进行稀疏分数阶傅立叶变换，得到不同分数阶傅立叶旋转角度下的幅度值，对幅度值取峰值即得出导频信号产生的PIM干扰信号在分数域峰值的位置和大小，进而估计导频信号产生的PIM干扰信号携带的参数信息；

步骤三：利用步骤二估计得到的导频信号产生的PIM干扰信号携带的参数信息和下行信号重建PIM干扰信号；

步骤四：判断步骤三重建的PIM干扰信号的幅值是否超过阈值，并决定跳至步骤五还是步骤六；

步骤五：采用自适应滤波算法估计PIM干扰信号由双工器发射到接收的耦合信道冲激响应，得出经自适应滤波算法处理后的PIM干扰信号；

步骤六：修正上行链路信号，完成PIM干扰对消；

至此，从步骤一到步骤六，完成了一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法。

2.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：

步骤一中，下行信号包括导频信号和非导频信号，导频信号产生的PIM干扰信号s_PIM(t)可以表示为如下公式(1)：

$s_{PIM}\left(t\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{s}_p{\left(t\right)}^l---\left(1\right)$

其中，s_p(t)为导频信号，表示对每一项a_ls_p(t)^l进行求和，a_l表示幂级数模型阶次l的幅度值，即：导频信号产生的PIM干扰信号采用幂级数模型表示，且下行信号中的导频信号是线性调频信号，又称为LFM信号；

接收的上行链路信号包括上行信号、PIM干扰信号及噪声三部分，可以用如下公式(2)来表示：

s_R(t)＝s_UP(t)+I_PIM(t)*h(t)+n₀(t) (2)

其中，s_R(t)为接收的上行链路信号，s_UP(t)为上行信号，I_PIM(t)为下行信号产生的PIM干扰信号，包括导频信号产生的PIM干扰信号s_PIM(t)及非导频信号产生的PIM干扰信号，h(t)表示PIM干扰信号由发射到接收的耦合信道冲激响应，I_PIM(t)*h(t)表示双工器接收到的PIM干扰信号，其中*表示卷积运算，n₀(t)为噪声。

3.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：步骤二中，导频信号产生的PIM干扰信号s_PIM(t)携带参数信息：幂级数模型阶次l和幂级数模型阶次l的幅度值a_l，记为和

步骤二，具体为：

步骤2.1：对导频信号产生的PIM干扰信号进行量纲归一化，输出量纲归一化后的导频信号产生的PIM干扰信号，简称归一化PIM信号；

首先，离散化表示导频信号产生的PIM干扰信号s_PIM(t)，通过如下公式(3)：

$s_{PIM}\left(n\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{s}_p{\left(n\right)}^l=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{e}^{jl(2{\mathrm{\πf}}_{0}n/f_s+{\mathrm{\π\βn}}^2/{f_s}^2)}---\left(3\right)$

其中，s_PIM(n)为导频信号产生的PIM干扰信号的离散化表示，n代表离散化变量，s_p(n)为导频信号s_p(t)的离散化表示；j代表虚数，f₀表示LFM信号的起始频率，β表示LFM信号的调频率，采样间隔为1/f_s，导频信号产生的PIM干扰信号的频域和时域范围分别为[-f_s/2，f_s/2]和[-T/2，T/2]；

然后，引入一个具有时间量纲的归一化因子定义量纲归一化坐标p＝t/S，q＝f·S，借助量纲归一化坐标，将频域和时域范围两个区间：[-f_s/2，f_s/2]和[-T/2，T/2]都归一化为[-Δx/2，Δx/2]，

归一化PIM信号可用如下公式(4)表示：

${s_{PIM}}^{\′}\left(n\right)=\underset{l=1}{\overset{\mathrm{\∞}}{\Σ}}{a}_l{e}^{jl(2{\mathrm{\πf}}_{0}n/\mathrm{\Δ}x+{\mathrm{\π\β}}^{\′}{n}^2/{\mathrm{\Δx}}^2)}---\left(4\right)$

其中，s_PIM′(n)表示归一化PIM信号，为归一化PIM信号的调频率，1/Δx为归一化PIM信号的采样间隔；

步骤2.2：对步骤2.1输出的归一化PIM信号进行稀疏分数阶傅立叶变换，具体采用Pei采样型离散算法进行稀疏分数阶傅立叶变换，得到稀疏分数阶傅立叶变换结果；

采用Pei采样型离散算法进行稀疏分数阶傅立叶变换的原因在于：由于导频信号产生的PIM干扰信号为多分量LFM信号，此多分量LFM信号在分数域具有稀疏性，可以用稀疏傅立叶变换(Sparse Fourier Transform，简称SFT)代替Pei采样型离散算法中的傅立叶变换，达到降低运算复杂度的目的；包括如下步骤：

步骤Pei.1：将步骤2.1输出的归一化PIM信号与chirp1信号相乘；

相乘后得到序列x(n)用如下公式(5)表示：

$x\left(n\right)=e^{j\frac{{\mathrm{cot\αn}}^2}{2{\mathrm{\Δx}}^2}}{s_{PIM}}^{\′}\left(n\right)---\left(5\right)$

其中，为chirp1信号，cot为反正切函数，α为稀疏分数阶傅立叶变换的旋转角度；

步骤Pei.2：将步骤Pei.1输出的归一化PIM信号与chirp1信号相乘的结果进行稀疏傅立叶变换，得到归一化PIM信号与chirp1信号相乘结果的稀疏傅立叶变换的估计值，简称稀疏傅立叶变换的估计值，具体为：

步骤Pei.2.A：对x(n)进行重新排序，得到重排后的序列，即定义重排方式为P_σ，τ，重排后序列用如下公式(6)表示：

(P_σ，τx)_i＝x_(σi+τ)\N (6)

其中，(P_σ，τx)_i为重排后序列P_σ，τx中的第i个值，σ是从[1，N]中随机选取的奇数，τ是从[1，N]中随机选取的整数，N为采样点数且N＝(Δx)²，(σi+τ)\N表示对索引值σi+τ的模N操作；

步骤Pei.2.B：将步骤Pei.2.A输出的重排后序列与滤波器相乘，得到相乘后序列，具体为：

所述滤波器采用低通滤波器，此滤波器的传输函数为G(∈，∈′，δ，ω)，此滤波器的傅立叶变换需要满足如下公式(7)：

${\hat{G}}_i=\left\{\begin{array}{cc}\[1-\mathrm{\δ},1+\mathrm{\δ}\]& i\∈\[-{\mathrm{\ϵ}}^{\′}N,{\mathrm{\ϵ}}^{\′}N\]\\ \[0,\mathrm{\δ}\]& i\∈\[-\mathrm{\ϵ}N,\mathrm{\ϵ}N\]\end{array}\right.---\left(7\right)$

其中，表示滤波器傅立叶变换的第i个值，∈′和∈分别表示通带截止频率和阻带截止频率；通带和阻带的最大衰减相同，用δ表示；ω表示滤波器的长度；

相乘后序列可用如下公式(8)表示：

y_i＝G_i×(P_σ，τx)_i (8)

其中，y_i为相乘后序列的第i个值，G_i为滤波器传输函数的第i个值；

步骤Pei.2.C：对步骤Pei.2.B输出的相乘后序列进行时域混叠，得到混叠后序列，并计算此混叠后序列的傅立叶变换；

混叠后序列可用如下公式(9)表示：

$z_i=\underset{jj=0}{\overset{N/B-1}{\Σ}}{y}_{i+Bjj}---\left(9\right)$

其中，z_i为混叠后序列z第i个值，jj表示求和变量，范围是0到N/B-1，B为混叠长度，参考值是(k为稀疏度)；

记混叠后序列的傅立叶变换结果为

定义散列函数用如下公式(10)表示：

h_σ(i)＝round(σ×i×B/N) (10)

其中，h_σ(i)为散列函数的第i个值，round(σ×i×B/N)表示对σ×i×B/N进行四舍五入运算；

定义偏移函数用如下公式(11)表示：

o_σ(i)＝σ×i-h_σ(i)×(N/B) (11)

其中，o_σ(i)为偏移函数的第i个值；

步骤Pei.2.D：取出中幅度最大的前d×k个幅值及其坐标，其中，d为稀疏度增益，参考值是将取出的坐标保存在集合J中；再通过散列函数h_σ(i)∈J折算出坐标i∈[1，N]，将折算出的坐标保存在集合I中；

步骤Pei.2.E：利用步骤Pei.2.D折算出的坐标i计算原始序列x的傅立叶变换的估计值，可以表示为如下公式(12)：

${{\hat{x}}_i}^{\′}={\hat{z}}_{h_{\mathrm{\σ}}\left(i\right)}{\mathrm{\ω}}^{\mathrm{\τ}\×i}/{\hat{G}}_{o_{\mathrm{\σ}}\left(i\right)}---\left(12\right)$

其中，表示原始序列傅立叶变换估计值中的第i个值，为中第h_σ(i)个值，h_σ(i)取集合J中保存的坐标，ω表示 为步骤Pei.2.B中滤波器的傅立叶变换的第o_σ(i)个值；

步骤Pei.2.F：记步骤Pei.2.D输出的I为I₁，重复步骤Pei.2.A到步骤Pei.2.E L-1次，对每次循环步骤Pei.2.D得到的集合I顺序编号，第R次循环步骤Pei.2.D得到的集合I记为I_R+1，即：步骤Pei.2.A到步骤Pei.2.E共进行L次运算，第r次运算步骤Pei.2.D得到的集合I为I_r；

其中，L表示运算次数，其参考值是L＝log₂N；R为循环变量，变化范围为1到L-1，r为运算次数变量，变换范围为1到L；

步骤Pei.2.G：统计L次运算步骤Pei.2.D得到的集合I，记为I₀＝I₁∪I₂∪I₃∪...∪I_L，统计每个折算出来的坐标i∈I₀出现的次数，保存在集合s_i中，即s_i＝|{r|i∈I_r}|，取出现次数超过L/2次的坐标，保存在集合I′中，即I′＝{i∈I₀|s_i≥L/2}；

步骤Pei.2.H：对集合I′中每个坐标i，从每次运算中取出在步骤Pei.2.E计算得到的原始序列x的傅立叶变换的估计值，记为并取中位值作为最后的估计值，得到的估计值用如下公式(13)表示：

${\hat{x}}_i^{\′\′}=median\left(\right\{{\hat{x}}_i^r|i\∈I^{\′}\})---\left(13\right)$

其中，表示稀疏傅立叶变换的估计值，表示在i∈I′条件下对取中位值；

步骤Pei.3：将步骤Pei.2计算得到的稀疏傅立叶变换的估计值与chirp2信号相乘，得到稀疏分数阶傅立叶变换结果；

得到稀疏分数阶傅立叶变换的结果用如下公式(14)表示：

$\widehat{F_{\mathrm{\α}}}\left(m\right)={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}\sqrt{\frac{(\sin \mathrm{\α}-j\cos \mathrm{\α})\mathrm{sgn}\left(\sin \mathrm{\α}\right)}{M}}={\hat{x}}_i^{\′\′}{e}^{j\frac{1}{2}{\mathrm{cot\αm}}^2{\mathrm{\Δu}}^2}{A}_{\mathrm{\α}}---\left(14\right)$

其中，表示稀疏分数阶傅立叶变换的结果，表示chirp2信号，m为离散化变量，sin为正弦函数，cos为余弦函数，sgn为符号函数，A_α表示M为稀疏分数阶傅立叶变换后信号的采样点数，需要满足M≥N，Δu为稀疏分数阶傅立叶变换后信号的采样间隔，需要满足如下公式(15)：

$\frac{\mathrm{\Δ}u}{\mathrm{\Δ}x}=\frac{2\mathrm{\π}\left|\sin \mathrm{\α}\right|}{N}---\left(15\right)$

步骤2.3：寻找步骤2.2得到的稀疏分数阶傅立叶变换结果的峰值位置和大小，估计导频信号产生的PIM干扰信号携带的参数信息；

其中，稀疏分数阶傅立叶变换的结果即步骤Pei.3的输出导频信号产生的PIM干扰信号携带的参数信息即幂级数模型阶次l及其幅度值a_l的估计值：和按照式(16)和式(17)估计：

$\hat{l}=-\frac{{\mathrm{cot\α}}^{\′}}{2{\mathrm{\πk\β}}^{\′}}=-\frac{{\mathrm{cot\α}}^{\′}{f_s}^2}{2{\mathrm{\π\β\Δx}}^2}---\left(16\right)$

$a_{\hat{l}}=\frac{\max \left|\widehat{F_{\mathrm{\α}}}\left(m\right)\right|}{(2N+1)\left|A_{{\mathrm{\α}}^{\′}}\right|}---\left(17\right)$

其中，表示幂级数模型阶次l的估计值，表示幂级数模型幅度值α_l的估计值，α′为的峰值所对应的稀疏分数阶傅立叶变换的旋转角度，表示取的最大值，表示对进行绝对值运算，|A_α′|表示对A_α′进行绝对值运算，A_α′为的峰值所对应的A_α的值。

4.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：

步骤三中，重建的PIM干扰信号用如下公式(18)表示：

${\hat{I}}_{PIM}\left(t\right)=\underset{\hat{l}}{\Σ}{a}_{\hat{l}}{s}_D{\left(t\right)}^{\hat{l}}---\left(18\right)$

其中，为重建的PIM干扰信号，s_D(t)为下行信号，包括导频信号和非导频信号。

5.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：步骤四，具体为：

步骤4.1：若超过阈值，跳到步骤五采用横向滤波器结构进一步进行PIM干扰自适应滤波；

步骤4.2：若不超过阈值，直接输出步骤三重建的PIM干扰信号并跳至步骤六；

其中，阈值的选取由PIM干扰对消后的误码率来决定。

6.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：步骤五中，自适应滤波算法可为最小均方算法(LMS)，采用M-1个延迟单元的级联横向滤波器实现，M表示横向滤波器的长度，通过不断调整横向滤波器的系数使得均方误差最小；经过自适应滤波算法处理后的PIM干扰信号表示为

7.根据权利要求1所述的一种基于稀疏分数阶傅立叶变换的无源互调干扰对消方法，其特征在于：步骤六，具体为：

步骤6.1：若进行了自适应滤波算法处理，则用步骤五得到的自适应滤波算法处理后的PIM干扰信号修正上行链路信号，修正后的上行链路信号用如下公式(19)表示：

${s_R}^{\′}\left(t\right)=s_R\left(t\right)-{{\hat{I}}_{PIM}}^{\′}\left(t\right)---\left(19\right)$

其中，s_R′(t)表示修正后的上行链路信号；

步骤6.2：若未进行自适应滤波算法处理，则用步骤三得到的重建PIM干扰信号修正上行链路信号，修正后的上行链路信号用如下公式(20)表示：

${s_R}^{\′\′}\left(t\right)=s_R\left(t\right)-{\hat{I}}_{PIM}\left(t\right)---\left(20\right)$

其中，s_R″(t)表示修正后的上行链路信号。