基于分数阶傅里叶变换的线性调频干扰消除方法
技术领域
本发明主要涉及到信号处理领域，尤其涉及一种基于分数阶傅里叶变换的线性调频干扰消除方法。
背景技术
战场环境下，接收机可能接收到多个平稳以及非平稳干扰，其中线性调频干扰LFM(Linear Frequency Modulation)是一种普遍的非平稳干扰形式，对于LFM干扰，运用分数阶傅里叶变换进行处理可以得到较好的效果。文献1“王鹏,邱天爽,李景春等.基于高斯加权分数阶傅里叶变换的LFM信号参数估计[J].通信学报,2016,37(4):107-115.”提出了一种高斯加权分数阶傅里叶变换方法，应用于低占空比时的LFM信号参数估计，可以有效提高此种情况下的参数估计精度。文献2“张玉恒,吴启晖,王金龙.基于时频加窗短时傅里叶变换的LFM干扰抑制[J].电子与信息学报,2007,29(6):1361-1364.”提出了一种基于时频加窗短时傅里叶变换(TFW-STFT)的LFM干扰抑制算法，所提出的时频窗对LFM干扰具有较好的频域能量聚集性能，因此加窗对信号的影响要小于无聚集性能的短时傅里叶变换，该算法在信噪比损失和系统误比特率上明显优于基于短时傅里叶变换的算法。但是在实际应用过程中发现，在估计LFM信号的初始频率时，往往会出现偏差。
运用分数阶傅里叶变换可以较准确地估计LFM干扰的参数且运算复杂度比Wigner变换低，在众多文献中得到应用，但是很少有研究频率估计偏差的文献。偏差与误差的性质不同，是由LFM干扰的调频斜率决定的偏差量，不能忽略不计。
分数阶傅里叶变换的定义式如下。设接收信号为x(t)，信号x(t)阶数为p的分数阶傅里叶变换FRFT(Fractional Fourier Transform)定义为
式中Kp(t,u)为分数阶傅里叶变换的核函数，其中p为FRFT的的阶次。Kp(t,u)的定义如下：
式中n为整数，α表示旋转角度，α与阶次P的关系为α＝Pπ/2。上式的意义是：1、对于一般情况，也即α非π的整数倍时，核函数采取(2)式第一行的形式；2、对于特殊情况，也即α为π的整数倍时，核函数采取冲激函数的形式，即如(2)式第二、三行所示。
以上为分数阶傅里叶变换的定义，其运算复杂难以实现，采用文献3即Ozaktas提出的采样型离散分数阶傅里叶变换(DFRFT:Discrete Fractional Fourier Transform)算法“Haldun M.Ozaktas.Digital computation of the fractional Fourier transform,IEEE Trans on Signal Processing,44(9):2141～2150”，可以化繁为简。
发明内容
本发明要解决的技术问题是提高LFM干扰频率参数估算精度，从而使在此基础上构建的干扰子空间提高消除干扰的精度，提出了一种基于分数阶傅里叶变换的线性调频干扰消除方法。
为解决该问题，所采用的技术方案是：
步骤1：对所接收的信号进行量纲归一化；
步骤2：估计LFM干扰的调频斜率；
步骤2.1：将所述量纲归一化后的接收信号序列通过分数阶傅里叶变换映射到以变换域变量u为横轴，阶数p为纵轴的平面上；
所述分数阶傅里叶变换为：
步骤2.1.1：取接收信号x(t)，以1/Δy为采样间隔进行采样，设长度为L，得到接收信号数据序列
步骤2.1.2：对接收信号数据序列
步骤2.1.3：取切普信号exp[-jπt2cot(α)]，以1/(2Δy)为采样间隔进行采样，其长度与
步骤2.1.4：因对于x(t)与exp[-jπt2cot(α)]的采样间隔不同，后者采样间隔是前者一半，所以需要对x(t)采样之后的数据序列进行内插，使用香农内插公式对x(t)采样之后的数据序列
步骤2.1.5：对信号g(n)进行卷积运算；
将信号g(n)与因子
步骤2.1.6：将
其中Aα为常数，表示如下：
步骤2.2：使分数阶p在区间(0,2)之间递增，每递增一步进行一次分数阶傅里叶变换，得到接收信号在分数阶变换后的变阶分数阶傅里叶变换图谱；
步骤2.3：在变换图谱上搜索谱峰，谱峰位置的纵坐标对应LFM干扰的调频斜率，设某个谱峰所在点纵坐标为p1，则对应的旋转角度为α1＝p1·π/2，调频斜率即为c＝-cotα1；
步骤2.4：通过步骤2.3找到所有谱峰，得到所有LFM干扰的调频斜率；
步骤3：估计各LFM干扰的初始频率；
设某个LFM干扰对应谱峰的横坐标为u1，则其初始频率为
上式等号右边没有未知量，其中T取归一化之后的值1；
步骤4：根据所估计得到的各LFM干扰的调频斜率和初始频率，运用子空间投影消除干扰。
进一步地，步骤2.2中，0.5≤|p|≤1.5，当p∈(0,0.5)时，首先对接收信号进行傅里叶反变换，然后对结果运用适用于0.5≤|p|≤1.5时的FRFT处理；当p∈(1.5,2)时，首先对接收信号进行傅里叶变换，然后对变换后的结果运用适用于0.5≤|p|≤1.5时的FRFT处理。
进一步地，运用子空间投影消除干扰的方法是：
步骤5.1：根据估计的LFM干扰参数构建干扰信号序列；
用修正后的频率参数构建的信号为xi(n)＝expjπ(cin2+2fin)，i＝1,2,…,K，其中K表示干扰个数，ci、fi分别表示第i个LFM干扰的调频斜率和初始频率，如ci＝0，则该干扰为单频干扰；
取xi(n)序列的L长度取一段数据作为子空间基矢量：
步骤5.2：构建干扰子空间及其正交子空间。
用时域扩展后的干扰信号矢量构建成一个矩阵，
B＝[X1 X2 … XP]
该组矢量张成了干扰的子空间，其正交子空间表示如下:
I为单位矩阵；
步骤5.3：进行子空间投影消除干扰；
将接收信号矢量投影到干扰子空间的正交子空间，如下式
Y＝PX＝PXs+PV
式中X为接收信号矢量，通过投影过程消除了接收信号中的干扰分量，剩余有用信号和白噪声分量Xs和V。
本发明还提供了一种基于分数阶傅里叶变换的线性调频干扰消除系统，包括处理器，以及与所述处理器连接的存储器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时实现上述基于分数阶傅里叶变换的线性调频干扰消除方法的步骤。
本发明还提供了一种计算机可读介质，其上存储计算机程序，所述计算机程序可被处理器执行以实现上述基于分数阶傅里叶变换的线性调频干扰消除方法的步骤。
与现有技术相比，本发明所取得的有益效果是：
本发明基于分数阶傅里叶变换的线性调频干扰消除方法，在多干扰环境下，通过对信号进行分数阶傅立叶变换，在分数阶傅立叶变换的阶数对应的旋转角度递增时，得到变阶分数阶傅里叶变换图谱，通过扫描图谱中的谱峰，根据谱峰位置估计得到各干扰的调频斜率和初始频率。然后根据估计得到的干扰参数，消除干扰。本发明通过理论及实际采样的运算过程中分析发现，由于实际采样时，信号离散化，为使待处理的信号始终在运算窗口之内，需要在信号的前后补零，将导致卷积结果的前面一段接近于零。假设卷积结果在t0之前接近于零，而实际的信号流在零时刻已经开始了。而在t0时刻，信号x(t)的瞬时频率已不再是初始频率，在分数阶傅里叶变换的运算中，x(t)在此时的频率将被当作初始频率，因此所估计的信号的初始频率有偏差，不是真实的初始频率，从而带来运用子空间投影消除干扰的不准确，本发明在发现偏差量存在的情况后，给出了一个初始频率的估算方法，通过实验数据表明，在考虑了偏移量后，初始频率的估计更加精确，从而使得在多干扰的环境下，运用子空间投影可以更精确地消除干扰。
附图说明
图1为本发明系统流程图；
图2为变阶分数阶傅里叶变换，2(a)为阶数值变动时的分数阶傅里叶变换平面图，2(b)阶数值变动时的分数阶傅里叶变换立体图；
图3为基于分数阶傅里叶变换频率参数估计的LFM、单频干扰抑制效果。
具体实施方式
图1至图3示出了本发明基于分数阶傅里叶变换的线性调频干扰消除方法的一种具体实施例，如图1所示，包括以下步骤：
步骤1：对所接收的信号进行量纲归一化；
采用文献3中，运用Ozaktas的方法进行处理的过程包含了对时间和频率的采样，将两者量纲归一化便于分析与处理。即对于信号的时间t，有t∈[0,T]，因非对称区间难于计算，需将信号左移半个作用区间，即t∈[-T/2,T/2]。设频率范围为f∈[-F/2,F/2]。根据文献4“杨清红.线性调频信号参数的快速估计技术研究[D].华南理工大学,2015”的方法，通过引入一个尺度因子S来实现量纲的归一化。设尺度变换之后的新坐标为：
y＝t/S,v＝F*S
尺度变化之后，时域和频域表达式将被限定在新时宽T/S和新带宽F*S，现令尺度因子
步骤2：估计LFM干扰的调频斜率；
步骤2.1：将所述量纲归一化后的接收信号离散化，对离散后的接收信号序列通过分数阶傅里叶变换映射到以变换域变量u为横轴，阶数P为纵轴的平面上；
所述分数阶傅里叶变换为：
由背景技术中的式2可知，当α为π的整数倍时，分数阶傅里叶变换的核函数转化成为冲激函数，由此可以将定义式简化，无需讨论，所以讨论一般情况，即α取值范围为(0,π)∪(π,2π)，也即p的取值范围为(0,2)∪(2,4)时。在此情况下，(2)式表示的核函数只需要取第一项，在此情况下信号x(t)的p阶分数阶傅里叶变换为
式中Aα为常数，表示如下
本实施例中，对于实际运算过程中因采取离散化采样计算，具体实施过程为：
步骤2.1.1：取接收信号x(t)，以1/Δy为采样间隔进行采样，设长度为L，得到离散化后的接收信号数据序列
步骤2.1.2：对接收信号数据序列
由于实际采样时，信号离散化，为使待处理的信号始终在运算窗口之内，需要在信号的前后补零，将导致卷积结果的前面一段接近于零。假设卷积结果在t0之前接近于零，而实际的信号流在零时刻已经开始了。而在t0时刻，信号x(t)的瞬时频率已不再是初始频率，在分数阶傅里叶变换的运算中，x(t)在此时的频率将被当作初始频率，因此所估计的信号的初始频率有偏差，不是真实的初始频率，从而带来运用子空间投影消除干扰的不准确。
步骤2.1.3：取切普信号exp[-jπt2cot(α)]，以1/(2Δy)为采样间隔进行采样，其长度与
步骤2.1.4：因对于x(t)与exp[-jπt2cot(α)]的采样间隔不同，后者采样间隔是前者一半，对x(t)采样之后的数据序列进行内插，使用香农内插公式对x(t)采样之后的数据序列
对于接收信号x(t)，将接收信号x(t)与一个初始频率为0，调频斜率为t2cot(α)的线性调频信号exp(jπt2cot(α))相乘，如下式
g(t)＝exp(jπt2cot(α))x(t) (6)
为利于离散化计算g(t)，需在频域上将exp(jπt2cot(α))x(t)限制在无量纲带宽F内，旋转角度需限制在区间0.5≤|P|≤1.5内，此时(6)式可以用香农内插公式表示为：
根据文献3的分析，因子exp(jπt2cot(α))的带宽时频积为因子x(t)带宽时频积的两倍，为使之匹配，需要做相应设置。在转换到新的坐标之后，时间和频率共同的归一化采样区间宽度为Δy。如果x(t)的采样间隔为1/Δy，则g(t)的采样间隔应为
步骤2.1.5：对信号序列g(n)进行卷积运算，将信号g(n)与因子
步骤2.1.6：将
其中Aα为常数，表示如下：
α表示旋转角度，α与阶次p的关系为α＝pπ/2；
步骤2.1.5和步骤2.1.6也可以通过直接对连续信号g(t)进行计算，然后再进行离散化，具体如下：
对连续信号g(t)进行积分计算，与exp(-j2πutcscα)相乘之后求积分，积分的结果为自变量u的函数，不妨设该函数为g′(u)，
按照文献4“杨清红.线性调频信号参数的快速估计技术研究[D].华南理工大学,2015”的运算思路，这一积分类似于傅里叶变换，可推导运算如下：
式中rect(·)表示矩形函数，不妨设s为该函数自变量，当矩形函数自变量s满足-1≤s≤1时，有rect(s)≡1。而当0.5≤|P|≤1.5时，因变换域变量u满足|u|≤Δy/2，此时矩形函数的自变量在±1之间，恒等于1。则(9)式可简化为
将g′(u)与切普信号Aαexp(-jπu2cot(α))相乘就得到了FRFT表达式；
其中Aα为常数，表示如下：
α表示旋转角度，α与阶次p的关系为α＝pπ/2；
将g′(u)与切普信号Aαexp(-jπu2cot(α))相乘，如下式
在上式基础上，在u域对其进行采样，使之离散化，可以得到
运用半角公式之中的tan(α/2)＝cscα-cotα，可将上式转化为如下形式
令
步骤2.2：使分数阶p在区间(0,2)之间递增，每递增一步进行一次分数阶傅里叶变换，得到接收信号在分数阶变换后的变阶分数阶傅里叶变换图谱；
本实施例中，线性扫频信号的时频脊线在时频分布图上是关于其中点对称的，绕其中点翻转180度形状不变，所以α取值为(0,π)时已经遍历了所有情况，对应的，分数阶p取值区间即(0,2)。
由步骤2.1看出，在做分数阶傅里叶变换时，阶次p限制在区间0.5≤|p|≤1.5内，对于阶次未涵盖到的开区间(0,0.5)和开区间(1.5,2)，不妨设阶次p1∈(0,0.5)，阶次p2∈(1.5,2)，则有
式中Fp表示阶次为p的分数阶傅里叶变换，而根据分数阶傅里叶变换的性质，
步骤2.3：在变换图谱上搜索谱峰，谱峰位置的纵坐标对应LFM干扰的调频斜率，设某个谱峰所在点纵坐标为p1，则对应的旋转角度为α1＝p1·π/2，调频斜率即为c＝-cotα1；
一般情况下，p是未知的，可利用分数阶傅里叶变换的性质查找出p来。将p在区间(0,2)之间递增，在每一步进行一次分数阶傅里叶变换，当步进值递增到准确p值时，分数阶傅里叶变换会出现尖峰。
设所需要检测的LFM干扰信号表示式为
根据上文所设，p逐渐由0增加到2，在此过程中，LFM干扰对应的p值将会被找到。如图2(a)和2(b)所示阶数值变动时的分数阶傅里叶变换平面图和立体图。在此情况下对LFM干扰做分数阶傅里叶变换，有下式
由上式可知，当且仅当c＝-cotα时，积分式(10)内t的二次项可以被消除。因为根据(6)式的统计特性，只有当二次项被消除时，接收信号x(t)的p阶傅里叶变换才会出现谱峰。也即，当分数阶傅里叶变换的阶数对应的旋转角度递增，递增到与LFM干扰在时频面上的脊线垂直时，二次项将会被消除，此时LFM干扰的调频斜率c＝-cotα。
判断p值是否为1，如p值为1，分数阶傅里叶变换转换为傅里叶变换。将接收信号中的单频干扰映射为变阶分数阶傅里叶变换图谱中的谱峰，谱峰对应的横坐标即单频干扰的归一化频率值。
步骤2.4：通过步骤2.3遍历图谱中所有谱峰，得到所有LFM干扰的调频斜率；
步骤3：估计各LFM干扰的初始频率；
设某个LFM干扰对应谱峰的横坐标为u1，谱峰的纵坐标所对应阶数转换的旋转角为α1时，则其初始频率为
上式等号右边没有未知量，其中T取归一化之后的值1；
本实施例中，在估测了LFM干扰的调频斜率之后，将c＝-cotα代入(16)式，可以得到
但以上分析是在理论情况下，此时积分上下限为无穷，实际中积分区间是有线长度，同前文，设其区间为[0,T]，则有
如前文所述，非对称区间难于计算，需将待处理的信号左移半个作用区间，即t∈[-T/2,T/2]。信号左移半个区间之后，以LFM干扰信号而言，相当于表示式发生了变化，应表述如下
即调频斜率未发生变化，仍然是
由上式可以推出，当
需要指出的是：单频干扰可以看作特殊的LFM干扰，如接收信号中含有单频干扰，则在FRFT的阶数p变化到”1”时，即可出现谱峰，因为此时的分数阶傅里叶变换已经转化为一般的傅里叶变换。
步骤4：根据所估计得到的各LFM干扰的调频斜率和初始频率，运用子空间投影消除干扰。
运用子空间投影消除干扰的基本思想是：设想接收信号中含有的干扰构成一个子空间，找出该子空间最直接的方法就是找出子空间的基矢量。由于各干扰互不相关，所以构建的干扰序列可作为干扰子空间的基矢量，每一个干扰信号对应一个基矢量，以此可以构建干扰子空间。在构建了干扰子空间的基础上，也就可以运用相应的子空间运算达到消除干扰的目的。
步骤4.1：根据估计的LFM干扰参数构建干扰信号序列；
根据步骤2和步骤3所估计的LFM干扰的频率参数之后，构建干扰信号序列
xi(n)＝exp(j2π(fin+ci/2n2)),i＝1,2,…,K (21)式中K表示干扰个数，fi、ci分别表示第i个干扰的初始频率和调频斜率。取xi(n)序列的L′长度数据作为子空间基矢量，以此构建干扰子空间，在此基础上运用子空间投影消除干扰。
取xi(n)序列的L长度段数据作为子空间基矢量：
LFM干扰初始频率的估计结果存在误差，误差尽管微小，但也对子空间投影结果有影响。根据文献5”Zhai Yong-Ping,Zhou Zhu.Linear Frequency ModulationInterference Suppression Based on Error Correction.2019 IEEE9th InternationalConference on Electronics Information and Emergency Communication(ICEIEC),Beijing,China,2019,pp.1-5.”的分析：同样的误差对于子空间投影结果的影响与构建子空间基矢量的长度有关，在同等误差情况下，用长序列作为干扰子空间基矢量引起的干扰抑制效果下降明显，构成干扰子空间基矢量的序列过短也存在频谱泄露的问题。
据上可知，选择合适长度的序列构成干扰子空间基矢量，可以达成良好的抗干扰效果。因为所有的设计均需要经过大量测试实验，大量实验也可以大致确定参数估计的误差范围，也就容易找出最佳的构建子空间的序列长度。这一长度值不需要确切数字。
步骤4.2：构建干扰子空间及其正交子空间。
用时域扩展后的干扰信号矢量构建成一个矩阵，
B＝[X1 X2 … XP](23)
该组矢量张成了干扰的子空间，其正交子空间表示如下
I为单位矩阵；
步骤4.3：进行子空间投影消除干扰；
将接收信号矢量投影到干扰子空间的正交子空间，如下式
Y＝PX＝PXs+PV(25)
式中X为接收信号矢量，通过投影过程消除了接收信号中的干扰分量，剩余有用信号和白噪声分量Xs和V。
抗干扰处理效果分析：
采用场景一，测试本发明的抗干扰效果。将本发明所估计的LFM干扰频率参数用于构建干扰子空间，所估计的LFM干扰的频率参数为：f1＝0.0999、c1＝0.2；f2＝0.2001、c2＝0.2。而单频干扰的频率值估计很精确，两个单频干扰的归一化频率值为：-0.1、0.4。
现对比两种情况：一、用本发明方法估计的频率参数构建干扰子空间，在此基础上通过子空间投影消除干扰；二、用精确的频率参数构建干扰子空间，在此基础上运用子空间投影消除干扰。其中第一种情况即是对本研究最终效果的测试，第二种情况提供对比。因子空间构建与基矢量长度有直接关系，故以子空间长度为横坐标，输出信干噪比为纵坐标，结果如图3所示。
如图3，图中点划线、星划线分别对应于频率参数取精确值、估计值两种情况。由图3可知：1、当用于构建干扰子空间的频率参数取精确值时，序列长度增加对于输出SINR没有影响，始终保持很高的输出SINR；2、用本发明方法估计干扰的频率参数，并据此构建干扰子空间，在此情况下，当用于构建干扰子空间序列较短时，能够获得较高的输出SINR，随着序列长度增加，SINR将会下降，但是序列长度取到32，SINR仍能保持较高，有用信号很容易从噪声中恢复。
由此可知：用本发明方法可以很好的消除单天线接收的数个LFM干扰以及单频干扰。
一种基于分数阶傅里叶变换的线性调频干扰消除系统，包括处理器，以及与所述处理器连接的存储器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时实现上述基于分数阶傅里叶变换的线性调频干扰消除方法的步骤。
一种计算机可读介质，其上存储计算机程序，所述计算机程序可被处理器执行以实现上述基于分数阶傅里叶变换的线性调频干扰消除方法的步骤。
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。
