Study on High Efficient Algorithm for Cyclic Correntropy Spectral Analysis

1.   引言

周期调频(Periodic Frequency Modulation, PFM)干扰信号，如线性调频连续波、正弦调频连续波干扰，是全球卫星导航系统(Global Navigation Satellite System, GNSS)接收机面临的最常见的干扰类型之一^[1,2]，该类干扰可能产生于恶意干扰机、雷达或民用电台等，具有平均功率低、峰值功率高、带宽大的特点。鉴于宽带周期调频干扰的频发性以及GNSS服务在国防安全和国民经济领域的重要作用，研究人员从卫星导航信号与干扰在时频域、空域以及空时域的差异着手，研究了大量的干扰抑制方法^[3,4]，主要可以分为基于天线阵的干扰抑制方法和基于单天线的干扰抑制方法。基于天线阵的干扰抑制方法利用多个天线的空间分辨率，对多干扰处理能力较强，且对信号损伤较小，但是天线阵的空间成本和经济成本较高，不适用于空间狭小、预算较少的应用^[5]。因此，适用于单天线接收机的干扰检测与抑制方法仍是GNSS抗干扰领域的研究重点之一^[6]。

单天线接收机的干扰抑制方法主要将接收信号转换到变换域进行干扰检测，再利用滤波器或者消隐技术消除干扰成分^[7]。目前，时频分析类方法是单天线接收机处理宽带调频类干扰的最主要手段，典型的时频分析方法有^[3]：短时傅里叶变换(Short-Time Fourier Transform, STFT)、小波包变换(Wavelet Packet Transform, WPT)、Wigner-Ville分布(Wigner-Ville Distribution, WVD)、分数阶傅里叶变换(FRactional Fourier Transform, FRFT)。STFT无法有效地累积信号能量，且存在因窗口宽度不变带来的分辨率不足的问题；离散WPT容易引起频谱混叠和幅值失真现象；WVD等非线性变换带来的交叉项会影响多分量干扰的参数估计精度；离散FRFT的非正交性会使期望信号产生畸变，且其仅对线性调频干扰效能较好。为了进一步提高单天线接收机的抗干扰能力，文献[8]提出多尺度STFT，增强了周期调频干扰在时频域中的聚集程度，降低了消除干扰时对期望信号的损伤。文献[9]采用多个周期干扰频率拐点的短时R’enyi熵提高干扰的参数估计精度，缓解了干扰频率估计误差导致干扰残余较多的问题。文献[10]将平滑伪WVD获取的时频谱图进行能量重排，在增强干扰时频分布聚集性、提升时频参数估计准确率的同时缓解了消除干扰时对期望信号的损伤。虽然这些方法提高了GNSS接收机对宽带调频干扰的抑制能力，但干扰与期望信号在变换域内仍存在较大的交叠，消除干扰时会对期望信号产生较大的损伤。特别是随着电子技术发展，现有小型干扰机即可产生含有多个调频分量的干扰，进一步增加了GNSS接收机抗干扰的处理难度。为了增加周期调频干扰与卫星导航信号的分离度，文献[11]提出通过调整时间采样序列的顺序将分散在较大带宽上的周期调频干扰能量集中在少数频点内，并利用频域消隐技术消除干扰成分，有效地减小了期望信号的损伤，但是该方法需要较多的采样数据以缓解频域消隐时频谱泄漏的影响。

为了进一步缓解数据较少时频谱泄漏问题，并提高单天线接收机对抗多分量周期调频干扰的能力，本文提出一种基于周期截断数据矩阵奇异值分解(Singular Value Decomposition, SVD)的干扰抑制技术，首先通过周期截断的数据重排方法构建观测矩阵，把分散在较大带宽的能量集中到重排数据中的单个频点，以增强干扰的聚集性及其与卫星导航信号的分离度；进而，利用矩阵奇异值分解将观测矩阵中干扰信号与卫星导航信号映射进不同的投影子空间以消除干扰成分。该方法能够有效消除单/多分量周期调频干扰，并且对卫星导航信号损伤较小，可用于高精度卫星导航接收机。

2.   信号模型

单天线接收机的中频数字信号可以表示为

其中，$t = 1,2, \cdots$表示采样时刻，${s_n}$表示第$n$颗卫星的导航信号，$j$表示周期调频干扰信号，$\eta$表示接收机热噪声。

接收机处的第$n$颗卫星的导航信号表示为^[12]

其中，${P_n}$表示接收机处第$n$颗卫星导航信号的功率，${d_n}( \cdot )$表示导航电文数据序列，${c_n}( \cdot )$表示扩频码序列，${f_{d,n}}$表示含多普勒频移的中频载波频率，${\varPhi _n}$表示载波相位。

周期调频信号数学模型可以记为

其中，${f_{{M_k}}}(t)$是第$k$个调频分量的调频函数，其周期为${P_k}$，$A$和$\varphi$分别为干扰信号的幅度和相位，${f_c}$为干扰的中频载波频率。

3.   周期截断数据矩阵中干扰与卫星导航信号分布特性

将$K$个周期调频干扰的周期记为${P_1},{P_2}, \cdots ,{P_K}$，它们一定具有最小公倍数为$P$，则

其中，${Z_1},{Z_2}, \cdots ,{Z_K}$为正整数。

以$cP$($c$为正整数)为间隔对中频数据进行截断，形成如式(5)的观测数据矩阵

其中，$M{\text{ = }}\left\lfloor {{L_S}/cP} \right\rfloor$，$\lfloor ·\rfloor$表示向下取整数，${L_S}$为采样数据长度，$Q$为最大列数。

其中，$m = 1,2, \cdots ,M$,$q = 1,2, \cdots ,Q$。

则观测数据矩阵可以表示为

由式(3)得，$t + cP$时刻的线性调频信号可以表示为

式(8)表明，时间间隔为$cP$的周期调频干扰信号数据之间只相差一个复数因子${{\rm{e}}^{ - {\rm{i}}2\pi {f_c}cP}}$，因此，第$q$列数据可以表示为

其中，${j_{{\varphi _q}}}$,${\varphi _q}$分别为第$q$列数据对应的原数据及其初始相位。因此，重排后的数据是频率为${f_c}mcP$的单频信号。

周期调频干扰信号第$q$列数据可以表示为

其中，$j(q + (m - 1)cP)$表示第$q$列的第$m$个元素。

由式(10)可知，第$q$列数据中每个元素都可以由第1个元素与1个常数相乘得到，把第1个坐标值分别乘以$- {{\rm{e}}^{ - {\rm{i}}2\pi {f_c}(m - 1)cP}}$加到第$m$个元素，则$J_{q}$为

所以干扰信号观测矩阵${\boldsymbol{J}} {\text{ = }} {\text{[}}{{\boldsymbol{J}}_1} \;\; {{\boldsymbol{J}}_2} \;\; \cdots \;\; {{\boldsymbol{J}}_q} \;\; \cdots \;\; {{\boldsymbol{J}}_Q}]$可经过初等变换化为$1 \times Q$的矩阵，此时矩阵$J$表示为

通过矩阵计算可得矩阵${\boldsymbol{J}} \cdot {{\boldsymbol{J}}^{\rm{H}}}$的特征值矩阵为

又因为矩阵${\boldsymbol{J}}$的奇异值就是矩阵${\boldsymbol{J}} \cdot {{\boldsymbol{J}}^{\rm{H}}}$的特征值的平方根^[13]，则矩阵${\boldsymbol{J}}$的奇异值矩阵为

即对周期截断数据矩阵进行奇异值分解，干扰主要集中在第1个奇异值所对应的子空间中。

卫星导航信号一般使用与高斯随机特性类似的伪随机码作为扩频码，所以接收机处扩频后的期望信号的统计特性类似于白噪声的统计特性^[14]。

因为接收机热噪声$\eta (t)$为白噪声，当

时，$h(t)$仍具有白噪声的统计特性。

当期望信号与热噪声的和信号h(t)映射到干扰J的子空间上时，其可以表示为

H的第m行数据可以表示为

其中，$h(q + (m - 1)cp)$表示第m行的第q个元素。

${\boldsymbol{H}} \cdot {\boldsymbol{H}}^{\rm{H}}$可以表示为

显然，${\boldsymbol{H}} \cdot {\boldsymbol{H}}^{\rm{H}}$是一个对称矩阵，则其一定能进行对角化，且对角化后的元素为其特征值，因此，${\boldsymbol{H}} \cdot {\boldsymbol{H}}^{\rm{H}}$有M个特征值，且为非负实数。

又因为H的第$q$列数据可以表示为

同理，${\boldsymbol{H}}^{\rm{H}}\cdot {\boldsymbol{H}}$是一个Q×Q的对称矩阵，存在Q个特征值，且为非负实数。综上，矩阵${\boldsymbol{H}}\cdot {\boldsymbol{H}}^{\rm{H}}$和${\boldsymbol{H}}^{\rm{H}}\cdot {\boldsymbol{H}}$存在min(Q,M)个非0特征值，即H存在min(Q,M)个奇异值，则期望信号和噪声散布在干扰所在子空间中。

4.   基于周期截断数据矩阵奇异值分解的干扰抑制技术

周期调频干扰信号在周期截断数据矩阵中的各列数据里表现为单频干扰，可以通过频域消隐技术或滤波器来消除，然而，离散傅里叶变换会出现频谱泄露，严重影响干扰抑制效能。考虑到干扰和卫星导航信号在观测数据矩阵的SVD子空间上的分布特性差异，本文提出一种基于周期截断数据矩阵奇异值分解的干扰抑制技术，其框图如图1所示。首先，利用接收信号的多重自相关函数估计周期调频干扰的调制周期^[11]；其次，根据估计的周期$P$构建周期截断数据矩阵；然后，采用SVD去除重排数据中的干扰成分；最后，得到重排后的时域信号$y$，用于后续处理。

图  1  本文所提方法原理框图

下载: 全尺寸图片 幻灯片

4.1   所提算法实施方式

由式(5)可得，根据周期调频干扰周期，构建周期截断数据矩阵${\boldsymbol{X}}$为

对矩阵${\boldsymbol{X}}$进行奇异值分解，即

其中，${\boldsymbol{U}} = \left[ {{U_1}\;{U_2} \;\cdots } \right]$和${\boldsymbol{V}} = \left[ {{V_1}\;{V_2}\; \cdots } \right]$分别代表左奇异矩阵和右奇异矩阵，${\boldsymbol{\varSigma}}$代表奇异值矩阵，为${\boldsymbol{\varSigma}} = {\rm{diag}}\left\{ {{\lambda _1}}\;\;{{\lambda _2}}\;\;\cdots \right\}$，$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle\cdot}$}}{\lambda }$下角标表示主对角线的次序且${\lambda _1} \ge {\lambda _2} \ge \cdots \ge 0$。

由式(10)可知，矩阵中的各列数据中的干扰成分有相同的频率，即为单一频率干扰。文献[15]指出，如果仅仅每一列中的数据相位不同，则对应矩阵的秩为1。如果只有周期调频干扰成分，则${\boldsymbol{X}}$为秩为1的矩阵($\lambda = 0$)，换而言之，奇异值分解后干扰主要集中在第1个奇异值所对应的子空间中。当存在与干扰无关的其他信号成分(期望信号与噪声)时，上述结论仍然成立，并且期望信号和噪声将散布在整个空间中。

因此，将最大的奇异值${\lambda _1}$及其对应的左、右奇异向量剔除，获得新的奇异值矩阵$\widetilde {\boldsymbol{\varSigma}} = {\rm{diag}}\left\{ {{\lambda _2}}\;\;{{\lambda _3}}\;\;\cdots \right\}$，左奇异矩阵$\widetilde {\boldsymbol{U}} = \left[ {{{\boldsymbol{U}}_2}}\;\;{{{\boldsymbol{U}}_3}}\;\;\cdots \right]$，右奇异矩阵$\widetilde {\boldsymbol{V}} = \left[ {{{\boldsymbol{V}}_2}}\;\;{{{\boldsymbol{V}}_3}} \cdots \right]$；进而恢复消除干扰成分的数据矩阵

其中，$\widetilde {\boldsymbol{X}} = {[ {{{\tilde x}_1}}\;\;{{{\tilde x}_2}}\;\; \cdots \;\;{{{\tilde x}_q}}\;\; \cdots \;\;{{{\tilde x}_Q}} ]^{\rm{T}}}$，从而获得用于后续处理所需信号y为

另外，当存在多个不同周期的干扰成分时，由于干扰抑制处理时延与计算资源的限制可能无法以所有干扰的最小公倍数的数据长度进行截断。对此，可以借鉴文献[11]中迭代的方式对多个干扰进行逐个消除。

4.2   所提算法计算复杂度分析

所提算法主要有自相关运算、矩阵构建、矩阵分解、信号重构等操作。这些操作主要涉及乘法、加法、比较等基础运算，其中单次加法和比较运算消耗资源数以及二者运算次数远少于乘法运算。因此，以算法实施过程中所需乘法次数衡量复杂度，分析如下：(1)考虑到信号参数在短时间内几乎不变，每批次估计1次周期，其需要两次自相关运算，即$2{N_p}(4 + 6{\log_2}{N_p})$次乘法，其中${N_p}$为自相关运算数据长度；(2)构建周期截断数据矩阵和信号重构过程只涉及数据顺序的改变；(3)设干扰周期长度为$Q$，周期个数为$M$，即周期截断矩阵${\boldsymbol{X}}$的行数为$Q$，列数为$M$；奇异值分解第1步需要$MQ + 3{M^2} - 3M$次乘法得到矩阵${\boldsymbol{XX}}^{\rm{H}}$的特征值，再对非零特征值开根号得到矩阵$X$的奇异值，需要${M^3} + {M^2}$次乘法；第2步需要$M(M + 1) + Q(Q + 1)$次乘法得到左、右奇异向量；(4)消除第1个奇异值和其对应的左、右奇异向量后，得到不含干扰成分的数据矩阵$\widetilde {\boldsymbol{X}}$需要${(M - 1)^2} + {(Q - 1)^2}$次乘法。综上，设数据总长度为${N_I}$，迭代次数为${N_j}$，则总共需要乘法运算次数约为$2{N_p}(4 + 6{\log_2}{N_p}){\text{ + }}{N_j} r \left[ \left\lceil {{{{N_I}} /{(MQ)}}} \right\rceil ({\text{3}}MQ + {M^3} + {\text{3}}{M^2} - M + 2{Q^2} - Q + {\text{3}}) \right]$，$\lceil ·\rceil$表示向上取整。

5.   仿真实验

为验证所提方法的性能，采用全球定位系统(GPS) L1频段的C/A码信号，伪码速率为1.023 MHz，下变频后模拟信号中心频率为1.25 MHz，采样频率为5 MHz，信噪比为–15 dB。干扰形式分为单分量场景、双分量干扰1场景和双分量干扰2场景，各干扰场景干扰参数如表1所示。对比方法为：多尺度短时傅里叶变换(MSTFT)法^[8]、自适应小波包系数阈值(WPCT)法^[16]、时域联合分数阶傅里叶变换(FRFT)法^[17]、基于短时R’enyi熵的时频分析方法^[9]和基于数据重排的多分量周期调频干扰抑制方法^[11]。

表  1  各干扰场景干扰参数

干扰场景	调频率(GHz/s)	扫频周期(ms)	带宽(MHz)	起始频率，终止频率(MHz)
单线性调频干扰场景	10	0.2	2	0.1,2.1
双线性调频干扰场景1	10；–2	0.2；0.4	2；0.8	0.1,2.1；2,1.2
双线性调频干扰场景2	10；–2	0.2；0.32	2；0.64	0.1,2.1；2,1.36

下载: 导出CSV 

| 显示表格

对于MSTFT，采用3个尺度，并设置2阶IIR陷波滤波器的极点收缩因子为0.85；对于WPCT，采用“Dmey”母小波函数，小波分解层数为4，并采用软阈值进行干扰检测与消除；对于FRFT，为了搜索干扰信号最佳阶数，搜索的总次数$r$为2000，对于WPCT和FRFT方法，单次处理的数据长度${N_t}{\text{ = }}1024$个采样点。对于文献[9]方法，1阶IIR陷波滤波器的极点收缩因子为0.95。对于文献[11]方法和本文所提方法，估计干扰周期所用数据长度${N_p}$为$2 \times {10^4}$个采样点，干扰抑制性能与单次处理数据长度密切相关，因此设置不同的数据长度(单次处理的干扰周期个数衡量)以分析算法性能，数据长度分别为16, 64, 128个周期。

当输入干噪比(INR)在25～60 dB变化时，图2、图3和图4中的(a), (b)和(c)分别表示在单分量干扰、双分量干扰场景1和双分量干扰场景2下，通过100次蒙特卡洛实验，得到的干扰抑制处理后的GPS信号归一化均方误差(NMSE)，输出信干噪比(SINR)和捕获因子(AF)，捕获因子是指利用相干积分技术对GPS信号进行捕获后的最大相关峰值和第2相关峰值之比，其中相关积分时间为4 ms。

图  2  单分量干扰场景下各方法干扰抑制性能

下载: 全尺寸图片 幻灯片

图  3  双分量干扰场景1下各方法干扰抑制性能

下载: 全尺寸图片 幻灯片

图  4  双分量干扰场景2下各方法干扰抑制性能

下载: 全尺寸图片 幻灯片

采用MSTFT和WPCT方法进行干扰抑制时，随着INR增大，NMSE增大且SINR和AF减小，是由于这两种时频分析方法的时频分辨率受限，干扰能量越大，其在时频域能量扩散对期望信号影响越大；而基于FRFT的方法在一定程度上提高了周期线性调频干扰的能量聚集性，但数字FRFT存在频谱泄漏的问题，因此干扰抑制性能也随着干扰能量增大而降低。文献[9]的方法利用估计的瞬时频率将干扰信号变换到基带再采用陷波器进行干扰消除，固定参数的陷波器对期望信号的损伤是不变的，但随着干扰能量的增大，其泄露到陷波器外的能量增加，导致干扰抑制性能下降。文献[11]的方法在干扰周期数较少时导致重排后各组数据长度不足，其采用的自适应傅里叶变换截断方法失效，干扰频谱泄漏产生较大影响，且傅里叶变换的分辨率随着数据长度减小而降低，干扰抑制性能下降。本文所提方法的NMSE, SINR和AF总体优于对比算法，且随干扰能量增加恶化不明显，这是因为采用的周期截断重排方法能够有效地将干扰成分聚集在某单一频率，通过奇异值分解将干扰投影在单一子空间中，干扰聚集度较高，且不随干扰能量增大发生扩散；但周期数减少同样导致性能下降，这是因为期望信号成分是分散在所有投影子空间中(干扰子空间中也存在部分期望信号成分)，周期数越少，重排数据矩阵的维数越低，投影子空间个数越少，导致期望信号分散在单一子空间上能量越大，即在干扰子空间上信号成分随维数降低而增大，则在消除干扰时，对期望信号损伤增大。

图3给出了双分量干扰场景1下各算法的干扰抑制性能，表2列举了在输入INR为30 dB和45 dB条件下的评价指标相对于单分量干扰场景的变化情况，“↑”、“↓”分别表示增大和减小。特殊地，在双分量周期调频干扰情况下，两个干扰分量时频分布存在交叠，此时文献[9]提出的方法难以估计干扰分量的调频率，因此其不适用于抑制该类多分量周期调频干扰。相对于单分量干扰，随着干扰分量增多，在STFT, WPCT, FRFT变换后，与期望信号交叠的干扰成分也增加，则抑制干扰时期望信号损失增大，在相同INR条件下干扰抑制性能降低；由于多分量干扰可以通过公周期截断，转换为矩阵中各列数据的单频干扰，因此与单分量干扰场景相比，文献[11]与本文所提方法性能下降不明显。

表  2  单分量干扰场景与双分量干扰场景1性能变化对比

	MSTFT		WPCT		FRFT		文献[11]方法(16个周期)		本文方法(16个周期)
INR (dB)	30	45	30	45	30	45	30	45	30	45
ΔNMSE	↑0.3	↑0.3	↑0.4	↑0.3	↑0.2	↑0.1	0	0	0	0
ΔSINR	↓2.7	↓5.7	↓9.8	↓7.2	↓1.8	↓5.1	↓1.8	↓2.9	0	0
ΔAF	↓2.2	↓1.8	↓3.9	↓2.2	↓10.3	↓1.6	↓1.2	↓0.8	↓0.7	↓0.7

下载: 导出CSV 

| 显示表格

图4给出了双分量干扰场景2下各算法的干扰抑制性能，表3列举了在输入INR为30 dB和45 dB条件下相对于双分量干扰场景1的评价指标变化情况。与双分量干扰场景1下方法性能相比，STFT, WPCT, FRFT方法性能相似，文献[11]与本文所提方法性能明显下降。该场景2中干扰分量的公周期较大，文献[11]与本文方法无法通过单次处理进行干扰抑制，采用两次迭代逐一消除各干扰分量。因此不能将干扰集中在一个频率点上，不同的干扰可能会相互影响，从而导致一些干扰分量无法消除。不同干扰之间的相互作用也会增加GPS信号和干扰之间的重叠，因此干扰抑制效能也会随着干扰能量的增加而降低。但是，在INR<50 dB时所提算法性能仍然优于对比算法。

表  3  双分量干扰场景1与双分量干扰场景2性能变化对比

	MSTFT		WPCT		FRFT		文献[11]方法(16个周期)		本文方法(16个周期)
INR(dB)	30	45	30	45	30	45	30	45	30	45
ΔNMSE	0	0	0	↑0.1	0	↑0.1	↑0.1	↑0.1	0	0
ΔSINR	0	↓1.2	↓1.1	↓1.5	0	↓1.5	↓1.5	↓2.5	↓0.4	↓5.4
ΔAF	↓1.1	0	↓1.0	0	↓1.3	↓0.9	↓2.0	↓0.9	↓0.6	↓12.1

下载: 导出CSV 

| 显示表格

为验证所提方法在不同信噪比下的性能，以双分量干扰场景2为例展示各方法性能随信噪比变化情况。具体条件为：干信比为40 dB，信噪比以2 dB的间隔从–20 ～0 dB，通过100次蒙特卡洛实验，得到干扰抑制处理后的GPS信号的NMSE，SINR和AF，由图5可知，本文方法在不同能量的噪声环境中，干扰抑制效果均优于对比算法。

图  5  双分量场景2下干扰抑制效果随信噪比变化情况

下载: 全尺寸图片 幻灯片

以处理相同长度接收信号所需的乘法次数衡量上述实验中各方法的计算复杂度，各方法所需乘法次数的一般表示和在不同场景下处理${N_I}{\text{ = 3}} \times {\text{1}}{{\text{0}}^6}$个采样点的接收信号所需乘法总次数如表4所示。

表  4  计算复杂度对比

方法	计算复杂度	单干扰场景	双干扰场景1	双干扰场景2
本文方法	$2{N_p}(4 + 6{\log _2}{N_p}){\text{ + }}{N_j}\left[ {\left\lceil {{N_I}/(MQ)} \right\rceil ({\text{3}}MQ + {M^3} + {\text{3}}{M^2}{\text{ - }}M + 2{Q^2} - Q + {\text{3}})} \right]$	3.9×108	7.6×108	8.8×108
MSTFT	$2K{N_I}/({N_t} - v){N_t}{\log_2}{N_t} + {N_I}$	9.2×109	9.2×109	9.2×109
WPCT	$8{N_I}{N_t}{\log_2}{N_t}({2^5} - 1) + {({N_I})^2} + {N_I}$	9.0×1012	9.0×1012	9.0×1012
FRFT	${N_j}(r + 2{N_I}/{N_t})(6{N_t}{\log_2}{N_t} + 3{N_t}) + 2{N_I}$	5.1×109	5.1×109	1.0×1010
文献[9]	$2{N_I}{N_t}{\log_2}({N_t}) + {N_I} + {N_t}({N_t} - 2)$	6.1×1010	–	–
文献[11]	$2{N_p}(4 + 6 {\log_2}{N_p}) + {N_j}[8{N_L}(1 + 2 {\log_2}{N_L}) + 8\left\lceil { {N_I}/{N_B} } \right\rceil {N_B} {\log_2}{N_B}]$	9.9×107	9.9×107	1.9×108

下载: 导出CSV 

| 显示表格

其中，文献[11]方法中的${N_L}$, ${N_B}$分别表示DFT的估计长度和频域消隐的数据长度；MSTFT法中的$K$表示多尺度短时傅里叶变换的窗口尺度个数，$v$表示STFT变换窗口重叠点数。由于在短时间内干扰参数几乎不变，对于FRFT的干扰最优阶数估计、文献[11]方法和本文方法中干扰周期估计只运行1次。表4中本文算法和文献[11]方法以16个干扰周期长度为例进行分析。本文方法、文献[11]方法和FRFT方法的计算复杂度与消除干扰所需的迭代次数成正相关；一般地，本文方法计算复杂度高于文献[11]方法，但小于其他对比方法。

6.   结束语

面向提高单天线GNSS接收机对抗周期调频干扰能力的需求，本文提出一种基于周期截断数据矩阵奇异值分解的干扰抑制技术。利用常见调频干扰的广义周期性以及周期截断数据矩阵特点，把在典型时频域内与期望信号交叠严重的干扰集中到与期望信号分离度较大的子空间中，采用奇异值分解的数据空间分解能力将干扰与期望信号映射到不同的子空间以消除干扰成分。仿真结果表明，在单分量和多分量周期调频干扰情况下，干扰抑制效果优于典型的时频分析类方法。