本发明属于雷达及水声信号处理，具体涉及一种用于数字信号处理的分数傅里叶变换高精度快速实现方法、系统及设备。
背景技术：
1、分数傅里叶变换(frft)作为一种流行的信号时频域分析方法在信号处理领域中得到广泛应用，它尤其适合于对线性调频信号做分析。通过对信号做分数傅里叶变换，将该信号投影到一个α-u二维变换域，其中α表示投影角度，它直接匹配信号的调频斜率，u与信号的时延具有一定的关联。通过分数傅里叶变换，可以容易地分析出信号的时频分布，便于信号的检测、识别等后续处理。
2、分数傅里叶变换如果直接按照定义式计算将具有较大的复杂度，不利于实际工程应用。为了有利于分数傅里叶变换的实际应用，目前已经衍生出了该方法的多种快速计算技术及系统。这些快速计算技术主要的原理为：针对每个角度α，将其对应的分数傅里叶变换分解为两个或多个信号的线性卷积形式，然后将这种线性卷积操作用快速傅里叶变换(fft)技术来实现快速计算，从而最终实现分数傅里叶变换的快速计算。上述不同技术的主要区别主要体现在对分数傅里叶变换分解的方式不同，但相同点主要体现为分解为不同信号分量的线性卷积形式。
3、在所述线性卷积的fft快速计算过程中，要用到多次fft及ifft(快速逆傅里叶变换)操作。具体来说，一般要将各个卷积项分别用fft变换到离散傅里叶频域，在离散傅里叶频域中将各个变换项相乘，然后再用ifft将相乘结果反变换回时域，从而获得线性卷积的快速计算结果。虽然fft与ifft技术计算效率较高，但多次应用fft与ifft，相比单次应用来说显然计算量要成倍增加。另一方面，现有frft快速计算技术中，被分解为线性卷积的两个或多个信号(或函数)主要为线性调频形式。针对宽带或超宽带信号处理领域，由于频带宽、以及信号序列较长等特点，被分解出来的信号或函数以及它们之间的运算会必不可免带来不可忽视的计算误差，这种计算误差会对结果有明显影响，导致计算精度损失。
4、公开号为cn 107644004a的中国发明专利申请，公开了一种基于离散分数阶傅里叶变换快速计算方法的数字信号处理方法及装置。但是由于该方法需要计算离散傅里叶变换核矩阵的特征值及特征向量，增加了计算复杂度，损失了计算效率。
5、公开号为cn105783974a的中国发明专利申请，公开了一种线性调频信号的检测、参数估计方法及系统，该方法通过旋转变换、及快速傅里叶变换(fft)来实现分数傅里叶变换的计算。但该方法对应的变换域计算结果分辨精度和变换域范围严格受限于信号的采样点数n和采样频率fs，丧失了设计的灵活性，并且它做的是整个采样频率范围fs上的fft，增加了计算冗余，同样损失了计算效率。
技术实现思路
1、为了克服上述现有技术的缺点，本发明的目的在于提供一种用于数字信号处理的分数傅里叶变换高精度快速实现方法、系统及设备，本方法通过对分数傅里叶变换进行全新形式分解，相比现有技术实现了更快速的计算，具体来说，针对每个角度α，经过一些处理步骤后，仅需使用一次fft即可获得最终的变换结果，较大程度上提高了计算效率。另一方面，本方法在对信号的处理过程中，所组合形成的各信号分量其形式趋于窄带化或单频化，从计算机或处理器实现角度来讲，能以更高的精度实现计算，显著提高了计算精确度。
2、为了实现上述目的，本发明采用的技术方案是：
3、一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，包括以下步骤：
4、s1，根据实际离散化输入信号xs(n)情况及实际需求，确定分数傅里叶变换的角度α的分析范围、离散化输入信号xs(n)的时延分析范围总大小tm、以及时延分析步长δtm；
5、s2，根据步骤s1得到的角度α范围、时延分析范围总大小tm及时延分析步长δtm，计算对应于最大频域范围的扩展分析点数ns以及点数扩展倍数p，并对离散化输入信号xs(n)进行补零处理；
6、s3，计算数字角频率分辨间隔及分数傅里叶变换u轴起始位置的索引值k0；
7、s4，根据步骤s3得到的数字角频率分辨间隔索引值k0，以及步骤s2得到的离散化输入信号xs(n)计算合成信号x's(n)；
8、s5，对步骤s4得到的合成信号x's(n)做抗混叠滤波及下采样处理，得到下采样后信号x'sld(n)；
9、s6，对步骤s5得到的下采样后信号x'sld(n)做快速傅里叶变换(fft)操作，获得分数傅里叶变换的初步结果x'α(kus)；
10、s7，对步骤s6获得的分数傅里叶变换的初步结果x'α(kus)做反转、下采样、截断操作，得到最终的分数傅里叶变换结果xα(kus)。
11、所述步骤s1具体为：
12、针对离散化输入信号xs(n)为瞬时频率随时间线性减小的情况，角度α的分析范围应在区间中选取，针对离散化输入信号xs(n)为瞬时频率随时间线性增长的情况，角度α的分析范围应在区间中选取；设定离散化输入信号xs(n)的时延分析范围为[0,tm]，其中tm表示时延分析范围总大小，时延分析步长为δtm，则可获得时延分析总点数为：
13、
14、其中，int(·)表示四舍五入取整操作；针对时延分析范围起始时刻不为零的情况，对离散化输入信号xs(n)进行超前或延时处理，使其分析起始时刻对齐至零时刻。
15、所述步骤s2具体为：
16、首先求出时延分析步长δtm对应的频率分析步长δfm，以及时延分析范围总大小tm对应的频域分析长度fm，按下式求解：
17、δfm＝δtm|cotα| (1-7)
18、fm＝tm|cotα| (1-8)
19、判断频域分析长度fm是否小于等于最大频域分析长度fs，即采样率，如果是，则继续向下执行，如果否，则返回式(1-7)及式(1-8)重设时延分析范围总大小tm或角度α；
20、按照下式计算对应于最大频域范围的扩展分析点数ns：
21、
22、其中，ns≥nm，nm为时延分析总点数，点数扩倍倍数p表示输入信号总点数n相比于最大频域范围的分析点数ns的倍数，按下式求得：
23、
24、其中，符号表示向上取整，n为离散化输入信号xs(n)的总采样点数；获得点数扩倍倍数p后，得到对应于最大频域范围的总分析点数为pns；
25、然后，根据总分析点数pns，对离散化输入信号xs(n)进行补零处理，即在xs(n)末尾补pns-n个零，使补零后信号的总点数刚好为pns，为简化表述，令补零后的离散化输入信号依然表示为xs(n)，此时时域序列索引n取值为0,1,2,...,pns-1。
26、所述步骤s3具体为：
27、数字角频率分辨间隔直接通过下式得到：
28、
29、其中，pns为总分析点数；
30、当角度α处于时频域第一象限时，索引值k0按下式计算：
31、
32、其中，fbm表示待分析信号xs(n)的中频起始频率，fbm满足关系：fbm＜fs，fs表示最大频域分析长度；
33、当角度α处于时频域第四象限时，索引值k0可按下式计算：
34、
35、其中，p表示点数扩倍倍数，δfm表示频率分析步长。
36、所述步骤s4具体为：
37、当角度α处于第一象限时，x's(n)按下式计算：
38、
39、当角度α处于第四象限时，x's(n)按下式计算：
40、
41、其中，j表示虚数单位，n表示时域序列索引，ts表示采样间隔。
42、所述步骤s5具体为：
43、首先根据下式计算抗混叠滤波器的截止角频率参数ωc：
44、
45、其中，fm表示频域分析长度，fs表示最大频域分析长度；
46、然后根据下式计算下采样倍数d：
47、
48、其中，符号表示向下取整，max(·)表示取一组数的最大值；
49、判断下采样倍数d是否等于1，如果是，则跳过抗混叠滤波及下采样处理，直接进行s6步骤，如果否，则继续执行抗混叠滤波及下采样处理；
50、设抗混叠滤波器长度为nl，nl为奇数，截止角频率为ωc，需要对抗混叠滤波器做超前移位处理，使滤波器的对称中心置于0时刻，即对应序列索引n＝0处，以保证抗混叠滤波器不会对信号产生额外延时，设按上述调整后的抗混叠滤波器为l(n)；
51、接下来对合成信号x's(n)做抗混叠滤波处理，得到滤波后信号x'sl(n)如下式所示：
52、
53、其中，表示循环卷积，对应的循环周期为总分析点数pns，滤波后信号x'sl(n)的长度保持为总分析点数pns不变，时域序列索引n＝0,1,2,...,pns-1；
54、然后对滤波后信号x'sl(n)做下采样处理，下采样后信号x'sld(n)表示为：
55、x'sld(n)＝x'sl(nd) (1-17)
56、其中，d表示下采样倍数，经过下采样处理后，时域序列索引n的取值范围缩短，其取值变为
57、所述步骤s6具体为：
58、按照下式获得分数傅里叶变换的初步结果：
59、
60、其中，j表示虚数单位，us表示u轴离散化步长，ts表示采样间隔，d表示下采样倍数，k'与k分别为修正前、后的变换域序列索引，当角度α处于第一象限时k＝k'-k0，当角度α处于第四象限时k＝k'+k0，其取值为式中的u轴离散化步长us通过下式得到：
61、
62、其中，p表示点数扩倍倍数，δfm表示频率分析步长，注意，如果未经过s5步骤，则(1-18)式中的信号x'sld(n)为x'sl(n)。
63、所述步骤s7具体为：
64、首先根据角度α判断分数傅里叶变换的初步结果x'α(kus)是否需要取反，如果角度α处于第一象限，则无需取反，如果角度α处于第四象限，则需对分数傅里叶变换的初步结果x'α(kus)沿u轴做取反操作，具体按下式进行：
65、
66、其中，pns为总分析点数，d为下采样倍数；
67、然后对x″α(kus)继续做下采样处理，下采样倍数即为点数扩倍倍数p，按下式进行：
68、xα(kus)＝x″α(kpus) (1-21)
69、最后对(1-21)式的结果xα(kus)做截断处理，提取前nm个点，至此，得到最终的分数傅里叶变换结果xα(kus)，其中修正后变换域序列索引k取值为：k＝0,1,2,...,nm-1。
70、一种用于数字信号处理的分数傅里叶变换高精度快速实现系统，其处理原理采用基于所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，包括：
71、合成信号生成模块，根据实际待处理信号情况及实际需求，确定分数傅里叶变换的角度α的分析范围、待分析信号的时延分析范围总大小tm、以及时延分析步长δtm，再根据角度α范围、时延分析范围总大小tm及时延分析步长δtm，计算对应于最大频域范围的扩展分析点数ns以及点数扩展倍数p，并对离散化输入信号xs(n)进行补零处理，然后计算数字角频率分辨间隔及分数傅里叶变换u轴起始位置的索引值k0；最后根据数字角频率分辨间隔索引值k0，以及离散化输入信号xs(n)计算合成信号x's(n)；
72、抗混叠滤波及下采样模块，对合成信号x's(n)做抗混叠滤波及下采样处理，得到下采样后信号x'sld(n)；
73、快速傅里叶变换(fft)模块，对下采样后信号x'sld(n)做快速傅里叶变换(fft)操作，获得分数傅里叶变换的初步结果x'α(kus)；
74、分数傅里叶变换后处理模块，对分数傅里叶变换的初步结果x'α(kus)做反转、下采样、截断操作，得到最终的分数傅里叶变换结果xα(kus)。
75、一种用于数字信号处理的分数傅里叶变换高精度快速实现设备，包括：
76、存储器：用于存储实现所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法的计算机程序；
77、处理器：用于执行所述计算机程序时实现所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法。
78、相对于现有技术，本发明的有益效果在于：
79、1.本发明在计算分数傅里叶变换过程中根据待处理的离散化输入信号xs(n)计算了合成信号x's(n)，相比直接利用xs(n)做处理的方式，实现了输入信号的窄带化或单频化处理，针对宽带或超宽带信号处理领域，该处理方式在数字化计算过程中避免了大数计算，使分数傅里叶变换结果精度更高。
80、2.本发明通过抗混叠滤波、下采样及快速傅里叶变换(fft)的系列操作，实现了在数据信息无损下数据量的显著降低，具有计算更快速的效果。
81、3.本发明采用修正的变换域序列索引k来作为最终的分数傅里叶变换结果xα(kus)的索引，使分数傅里叶变换的起始点索引刚好对应待分析输入信号的时延分析范围起始点，具有信号变换前后索引对应一致的特点，方便了实际使用。
82、4.本发明根据设计需求的输入信号时延分析范围和时延分析步长来计算频域范围的扩展点数及角频率分辨间隔等参数，设计过程及结果严格满足设计需求，而无需额外的后续处理，增加了设计的灵活性。
83、综上所述，本发明具有计算精度高、计算快速、使用方便灵活的特点。
1.一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：包括以下步骤：
2.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s1具体为：
3.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s2具体为：
4.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s3具体为：
5.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s4具体为：
6.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s5具体为：
7.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s6具体为：
8.根据权利要求1所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：所述步骤s7具体为：
9.一种用于数字信号处理的分数傅里叶变换高精度快速实现系统，其处理原理采用基于权利要求1～8任一项所述的一种用于数字信号处理的分数傅里叶变换高精度快速实现方法，其特征在于：包括：
10.一种用于数字信号处理的分数傅里叶变换高精度快速实现设备，其特征在于，包括：
