第5章离散算法




教学视频



分数阶傅里叶变换作为傅里叶变换的广义变换，是处理时变信号的有力工具。同时，FRFT与Wigner分布、小波变换等时频信号分析工具具有密切的关系，广泛应用于多个研究领域，如微分方程求解、光信号处理、扫频滤波器、模式识别和时频信号分析等。
正如快速傅里叶变换(Fast Fourier Transform,FFT)的出现大大推动了傅里叶变换的发展一样，FRFT的快速算法将是FRFT在信号处理中获得成功应用的基础。这使得对离散分数阶傅里叶变换(Discrete Fractional Fourier Transform,DFRFT)及其快速算法的研究显得尤为重要。
由FRFT的定义可以看出，DFRFT的计算将比离散傅里叶变换(Discrete Fourier Transform,DFT)的计算复杂得多。理想的DFRFT应该尽可能具备连续FRFT所具备的所有特性。同时，离散信号的DFRFT可以逼近其所对应连续信号的连续FRFT，是离散变换作为连续变换离散形式的最基本条件。此外，DFRFT的计算复杂度要尽可能得小，即实现离散变换的计算高效性。基于上述分析，理想的DFRFT需满足如下准则。
(1) 近似性： DFRFT与连续FRFT的近似性。
(2) 边界性： 1阶运算退化为DFT，即F1=F，F为DFT算子。
(3) 酉性： FpH=F-p。
(4) 旋转相加性： FpFq=Fp+q。
(5) 计算高效性。
近年来，国内外学者提出了多种DFRFT的定义及快速算法，然而遗憾的是，迄今为止，尚无一种定义能够很好地满足上述所有要求。目前的DFRFT主要分为三类： ①采样型DFRFT； ②特征分解型DFRFT； ③线性加权型DFRFT。因为采样型DFRFT和特征分解型DFRFT是实际应用中广泛使用的两类离散变换，本章将重点对这两类DFRFT进行系统的分析与比较。

5.1采样型离散分数阶傅里叶变换
一个直接并且简单的DFRFT定义方法是直接采样连续FRFT核得到DFRFT核矩阵。这方面的工作主要有： 
(1) Kraniauskas等通过在时域和分数阶傅里叶域直接采样得到DFRFT定义为［12］

Xα(kF)=Aαej12cotα(kF)2∑N-1n=0x(nT)ej12cotα(nT)2-j(2π/N)nk(5.1)

其中，T为时域采样间隔，F=2π/(NTcscα)是分数阶傅里叶域采样间隔。
(2) Ozaktas推导了两种高效并且精确计算FRFT的方法［3］。这种算法把时域原始函数的N个采样点映射为分数阶傅里叶域的N个采样点，并且这种算法的计算复杂度为O(NlogN)。
(3) Pei定义了另一种采样类型的DFRFT，这种方法对连续FRFT在时域和分数阶傅里叶域选择合适的采样间隔，使得DFRFT具有可逆性［4］。更重要的是，该DFRFT具有更低的计算复杂度，在所有近似连续FRFT的DFRFT类型中，这种DFRFT具有最低的计算复杂度。
因为Ozaktas和Pei的采样型算法目前使用最普遍，所以接下来对这两种算法进行详细介绍。




5.1.1Ozaktas采样型离散分数阶傅里叶变换
Ozaktas提出根据连续FRFT的积分定义式，将FRFT的复杂积分变换分解成若干简单的计算步骤,然后经过两步的离散化处理得到一个离散卷积的表达式，这样便可以利用FFT来计算FRFT。这里，称此方法为改进的采样型分数阶傅里叶变换(Improved DFRFT，IPDFRFT)。具体地，在文献［3］中,Ozaktas分别给出了两种不同的分解方法。在两种方法的公式推导中，都运用了技巧cscα-tan(α/2)=cotα。首先，我们介绍一下该算法所采用的一种特殊技巧，即原始信号在作FRFT的数值计算之前先进行量纲归一化处理，它对FRFT离散化的实现起到了重要作用。
5.1.1.1量纲归一化原理
首先谈谈信号的紧致性问题。若一个函数的非零值仅限定在一个有限区间内，则称这个函数是紧致的。从理论上讲,一个信号不可能在时域和频域同时紧致，然而在实际中所研究的信号一般认为是时宽带宽有限的。这种理论和实际的差异并不矛盾，因为只要信号的绝大部分能量集中在时频域的有限区域，就可以近似认为信号在时域和频域同时紧致。
假定原始信号在时间轴和频率轴上都是紧致的，其时域表示限定在区间-Δt2,Δt2，而其频域表示限定在区间-Δf2,Δf2，Δt和Δf分别表示信号的时宽和带宽。信号的时宽带宽积为N=ΔtΔf，信号的有限区间表示等价于假定信号的能量绝大部分集中在-Δt2,Δt2×-Δf2,Δf2内。对于一给定的函数类，这一假定可以通过选择足够大的Δt和Δf予以保证。由于原始信号的时域和频域具有不同的量纲，而且尺度也不统一，这给FRFT的离散化处理带来许多不便。分解型算法采用了一种特殊的量纲归一化技巧，将时域和频域分别转换成无量纲的域，并且将信号的时宽和带宽统一起来。具体方法是引入一个具有时间量纲的尺度因子S，并定义新的尺度化坐标x=t/S，ν=fS。新的坐标系x,ν实现了无量纲化。信号在新坐标系中被限定在区间-Δt/2S,Δt/2S和-ΔfS/2,ΔfS/2内。为使两个区间的长度相等，选择S=Δt/Δf，则两个区间长度都等于无量纲量Δx=ΔtΔf，即两个区间归一化为-Δx/2,Δx/2。归一化以后信号的WignerVille分布限定在以原点为中心、直径Δx的圆内，如图5.1所示。最后，根据采样定理对归一化后的信号进行采样，采样间隔为1/Δx，采样点数为N=Δx2。需要注意的是，今后出现在FRFT的表达式中的信号将是经过量纲归一化处理的信号。


图5.1归一化前后信号的时频支撑区域




5.1.1.2两种实用的量纲归一化方法
在5.1.1.1节的量纲归一化原理中，原始连续信号首先经过量纲归一化处理，然后再对连续函数以1/Δx为间隔采样得到N点样本值，这种处理方法的特点是先量纲归一化再采样。我们在实际应用中发现，这种处理方法只是一种原理性的方法，在实际工程应用中不具有可操作性，因为我们在实际工程中所能得到的往往不是原始连续信号，而是按照一定采样率进行采样后得到的离散信号。要想将分解型算法成功地应用于实际工程计算，就必须解决对这种实际的离散信号进行量纲归一化处理的问题。文献［5］针对此问题给出了两种实用的量纲归一化方法。
1. 离散尺度化法
所谓离散尺度化法是指，直接对离散数据作尺度伸缩变换，使得尺度化后的离散数据正好等价于对原始连续信号作量纲归一化后再采样所得的结果，如图5.2所示。其关键是要选择合适的时宽Δt、带宽Δf、尺度因子S以及归一化宽度Δx。信号的时宽比较容易确定，直接取为观测时间T，即Δt=T，同时以信号的中点作为时间原点，信号的时域表示限定在区间-T/2,T/2。信号的带宽确切值我们并不知道，但是在实际中我们知道信号的采样率fs。根据采样定理，采样率一定大于信号最高频率的2倍。信号带宽Δf的选取并不要求是最小值，只要满足将信号的全部能量包含在其中即可。我们将带宽直接取为采样率是完全合理的，即Δf=fs，信号的频域表示限定在区间-fs/2,fs/2。在确定了信号的时宽和带宽之后，可以得到尺度因子S和归一化宽度Δx分别为


S=Δt/Δf=T/fs(5.2)

Δx=ΔtΔf=Tfs(5.3)


离散数据原来的采样间隔为Ts=1/fs，对离散数据作尺度变换后，采样间隔变为

T′s=1/Tfs=1/Δx(5.4)

而原来的时域区间-T/2,T/2，经尺度变换后变成区间-Δx/2,Δx/2。因此，以采样频率为带宽，以观测时间为时宽，直接对离散数据作尺度伸缩变换，所得结果与原始连续信号作量纲归一化后再采样所得的结果完全相同。经过这样的分析处理，离散数据就可直接进行FRFT数值计算。


图5.2离散尺度化法示意图




2. 数据补零/截取法
离散尺度化法是通过对离散数据的在时间域上的伸缩来实现归一化。信号尺度的伸缩必然会导致原有信号的某些特征发生畸变，例如对一个chirp信号进行尺度伸缩将使它的调频率变大或变小。数据补零/截取法可以使原有信号不发生畸变而又实现量纲归一化，其关键仍然是选择合适的时宽Δt、带宽Δf、尺度因子S以及归一化宽度Δx。首先将时间原点定在数据的中点。为了保证原有信号不发生畸变，尺度因子只能选1，即S=1。先将时宽定为观测时间，即Δt=T，带宽定为采样频率，即Δf=fs。在确定归一化宽度Δx时，分两种情况。第一种情况： 若带宽值大于时宽值(fs＞T)，则Δx直接取两者的大值，Δx=fs。由于原始数据的采样间隔为1/fs，时间区间在-T/2,T/2，而归一化后要求采样间隔仍为1/fs，时间区间增加为-fs/2,fs/2。因此，通过在-fs/2,-T/2和fs/2,T/2区间以同样的采样间隔作数据补零来人为地增加信号的时宽，从而实现了信号的时宽带宽归一化，如图5.3所示，这就是数据补零法实现归一化的原理。第二种情况： 当时宽值大于带宽值(T＞fs)，则Δx取两者的小值，即Δx=fs。由于原始数据的采样间隔为1/fs，时间区间在-T/2,T/2，而归一化后要求采样间隔仍为1/fs，时间区间减小为-fs/2,fs/2。因此需要对原有数据作截取，只取出在区间-fs/2,fs/2内的数据，从而实现了信号的时宽带宽归一化，如图5.4所示，这就是数据截取法实现归一化的原理。


图5.3数据补零法示意图






图5.4数据截取法示意图




5.1.1.3第一种分解方法
为了方便起见,这里重写FRFT的定义式如下： 

Xp(u)=Aα∫+∞-∞expjπu2cotα-2utcscα+t2cotαx(t)dt(5.5)

其中,Aα=exp-jπsgnsinα/4+jα/2sinα1/2,α=pπ2。假定阶次p∈-1,1，将上式分为以下三步运算


g(t)=exp-jπt2tanα/2x(t)(5.6)

g′(u)=Aα∫+∞-∞expjπβu-t2g(t)dt(5.7)

Xp(u)=exp-jπu2tanα/2g′(u)(5.8)


其中，g(t)和g′(u)是两个中间结果，β=cscα，-π/2≤α≤π/2。要实现连续FRFT的数值计算必须对以上每个分解步骤都进行离散化处理，下面是具体的实现过程。
(1) 如式(5.6)所示，信号x(t)被一个线性调频信号exp-jπt2tanα/2所调制。为了对调制信号g(t)进行离散化处理，首先需要确定它的带宽。因为信号的时域支撑区间为-Δx/2,Δx/2，则线性调频信号的最高瞬时频率为tanα/2Δx/2,它的双边带宽为tanα/2Δx。因为信号x(t)与线性调频信号相乘对应于两者在频域的卷积，因此调制信号g(t)的总的双边带宽可确定为1+tanα/2Δx。当角度满足-π/2≤α≤π/2时，chirp调制信号g(t)的带宽最高可达到原信号x(t)带宽的2倍,即2Δx。为了满足采样定理，我们应当对g(t)以1/2Δx为间隔采样。如果x(t)的样本值的采样间隔为1/Δx，那么就需要对这些样本进行二倍插值，然后再与线性调频信号的离散采样值相乘，以得到所希望的g(t)的采样。关于插值的方法可以参考有关文献［6］。
(2) 如式(5.7)所示，信号g(t)与一个线性调频信号expjπβt2作卷积。因为g(t)是带限信号，所以线性调频信号也可以用其带限形式代替而不会有任何影响，也就是

g′(u)=Aα∫+∞-∞expjπβt-u2g(t)dt=Aα∫+∞-∞ht-ug(t)dt(5.9)

其中

h(u)=∫Δx-ΔxHνexpj2πνudν(5.10)

这里

Hν=1βejπ/4exp-jπν2/β(5.11)

是线性调频信号expjπβt2的傅里叶变换。其中，函数h(u)需要利用如下的Fresnel积分来求解

f(z)=∫z0expπz2/2dz(5.12)

于是，式(5.9)的离散形式为

g′m2Δx=Aα∑Nn=-Nhm-n2Δxgn2Δx(5.13)

这一离散卷积可以利用FFT快速计算。
(3) 根据式(5.8)得到FRFTXp(u)的以1/2Δx为采样间隔的样本值Xpm2Δx。由于假定x(t)的所有变换都是带限的，它们位于区间-Δx/2,Δx/2,所以需要对Xpm2Δx进行二倍抽取，以得到离散采样XpmΔx。
归纳起来，上述方法从唯一描述连续信号的x(t)的N个离散采样xnΔx开始，最后得到唯一描述Xp(u)的N个离散采样XpmΔx。如果令x-和X-p分别表示x(t)和Xp(u)的N个离散样本的列向量，则整个过程可以写作

p=FpIx-,FpI=DΛHlp
ΛJ(5.14)

式中，D和J分别是对应内插和抽取运算的矩阵； 矩阵Λ为对角矩阵，它对应为线性调频函数乘法； 矩阵Hlp对应卷积运算。FpI使得我们可以利用原函数的离散采样得到FRFT的离散采样，这是对DFRFT矩阵定义的基本要求。
本计算方法只适用于-1≤p≤1的情况。当阶次p位于该区间之外，则可利用FRFT的基本性质来得到所需结果。
5.1.1.4第二种分解方法
5.1.1.3节介绍的第一种分解算法存在一个问题，就是它必须计算Fresnel积分。本节介绍另一种分解算法则避免了这个问题。FRFT的定义式也可以改写作如下形式：

Xp(u)=Aαexpjπγu2∫+∞-∞exp-j2πβutexpjπγt2f(t)dt(5.15)

其中，γ=cotα，β=cscα，α=pπ/2。由上式可以明显看出，我们还可以将其分解成如下三个步骤来运算： 


g(t)=expjπγt2x(t)(5.16)

g′(u)=∫+∞-∞exp-j2πβutg(t)dt(5.17)

Xp(u)=Aαexpjπγu2g′(u)(5.18)


信号x(t)首先被chirp信号expjπγt2调制。为了对调制信号g(t)进行离散化处理，需要先确定它的带宽。因为x(t)为量纲归一化后的信号，它在所有分数阶傅里叶域(包括时域和频域)上的宽度都限定在区间-Δx/2,Δx/2内，或者说，它的Wigner分布限定在以原点为中心、直径为Δx的圆内。当限定阶次在0.5≤|p|≤1.5范围时，γ≤1,chirp调制信号expjπγt2x(t)的最高频率为0.51+γΔx≤Δx。这样,以1/2Δx为采样间隔并利用香农内插公式可将其表示为

expjπγt2x(t)=∑Nn=-Nexpjπγn2Δx2xn2Δxsinc2Δxt-n2Δx(5.19)

将式(5 19)代入式(5.15)，并交换积分和求和顺序，便得到


Xp(u)=Aαexpjπγu2∑Nn=-Nexpjπγn2Δx2xn2Δx·

∫+∞-∞exp-j2πβutsinc2Δxt-n2Δxdt(5.20)


上式中的积分项可以计算得到

∫+∞-∞exp-j2πβutsinc2Δxt-n2Δxdt=exp-j2πβun2Δx12Δxrectβx2Δx(5.21)

在0.5≤|p|≤1.5的范围，矩形函数rectβx2Δx在变换函数的支撑区|x|≤Δx/2将总是等于1。于是，我们可以写出


Xp(u)=Aα2Δx∑Nn=-Nexpjπγu2exp-j2πβun2Δx·

expjπγn2Δx2xn2Δx
(5.22)

上式中时域变量已经实现了离散化，而分数阶傅里叶域变量仍然保持连续。接下来需要对分数阶傅里叶域变量进行离散化。以1/2Δx为采样间隔，在-Δx/2,Δx/2内对分数阶傅里叶域变量采样，即令x=m/2Δx，代入上式得到


Xpm2Δx=Aα2Δx∑Nn=-Nexpjπγm2Δx2-j2πβmn2Δx2+

jπγn2Δx2xn2Δx,-N≤m≤N

(5.23)

这是一个有限求和，使得可以利用原函数的离散样本值求出FRFT的离散样本值。但是如果直接以上式作计算，计算复杂度为ON2，它的运算量仍然很大。将一个恒等式mn=12m2+n2-m-n2代入上式并经过一些化简后可得到


Fp［x］m2Δx=Aα2Δxexpjπγ-βm2Δx2·

∑Nn=-Nexpjπβm-n2Δx2

expjπγ-βn2Δx2xn2Δx,-N≤m≤N
(5.24)

式中的求和部分为离散卷积形式，该卷积可以用FFT快速计算，其总的计算复杂度为ONlogN。
和第一种方法一样，假定取二倍的插值和抽取，第二种方法也是从唯一表示函数x(t)的N个样本xnΔx出发，最后得到唯一表示FRFTXp(u)的N个样本XpmΔx。如果令x-和p分别表示x(t)和Xp(u)的N个离散样本的列向量，则第二种分解算法的整个过程可以用矩阵表示为

p=FpⅡx-，FpⅡ=DKpJ(5.25)

其中，


Kp(m,n)=Aα2Δxexpjπγm2Δx2-j2πβmn2Δx2+jπγn2Δx2,

|m|≤N,|n|≤N
(5.26)

与FpⅠ一样，FpⅡ也使得原函数的样本值转换为其FRFT的样本值。
虽然上述推导过程只是在假定0.5≤p≤1.5的条件下得到的，但是我们可以利用FRFT的旋转相加性，很方便地将阶次范围扩展到0≤|p|≤0.5或1.5≤|p|≤2范围时的情况。利用FRFT的旋转相加性

Fp=Fp-1+1=Fp-1F1(5.27)

我们可最后得到如下公式：


Xpm2Δx=Aα′2Δx∑Nn=-Nexpjπγ′m2Δx2-j2πβ′mn2Δx2+
jπγ′n2Δx2X1n2Δx
(5.28)

其中，α′=p-1π/2,γ′=cotα′,β′=cscα′，X1(u)表示x(t)的傅里叶变换。
IPDFRFT对连续FRFT具有很好的近似度，同时具有高效性。然而，IPDFRFT 无法满足酉性与可加性，从而导致当使用IPDFRFT 处理信号恢复的应用时存在一些近似误差。总的来说，IPDFRFT 仅满足理想离散傅里叶变换准则中的(1)，(2)和(5)。
5.1.2Pei采样型离散分数阶傅里叶变换
这种方法由Pei等提出［4］。与5.1.1节方法不同，虽然它也是从连续FRFT定义式出发，但是它不对连续FRFT表达式分解，而直接对输入/输出变量实现采样，然后通过限定输入/输出采样间隔来保持变换的可逆性，实现了具有解析表达式(ClosedForm Expression)的采样型离散分数阶傅里叶变换(CFDFRFT)。以下给出它的推导过程。将连续FRFT的定义式表达为

Xp(u)=1-jcotα2πej12u2cotα∫+∞-∞e-jutcscαej12t2cotαx(t)dt(5.29)

其中，p为变换阶次，α=pπ/2。为了推导DFRFT，我们首先对连续FRFT的输入函数x(t)和输出函数Xp(u)进行采样，采样间隔为Δt和Δu，得到

y(n)=xnΔt，Yp(m)=XpmΔu(5.30)

其中，n=-N,-N+1,…，N，m=-M,-M+1,…，M。这里我们不从t=0和u=0开始采样是因为我们希望使直流成分位于中心。将上式代入连续FRFT定义式，得到


Yp(m)=1-jcotα2πΔt
ej2m2Δu2cotα

∑Nn=-Ne-jmnΔtΔucscα
e
j2n2Δt2cotαy(n)
(5.31)

以上公式可以写为

Yp(m)=∑Nn=-NKp(m,n)y(n)(5.32)

其中

Kp(m,n)=1-jcotα2π
Δtej2
m2Δu2cotαe-j
mnΔtΔucscαej2
n2Δt2cotα(5.33)

为了使式(5.32)可逆，当M≥N时，我们需要使它的逆变换等于Kp(m,n)的Hermitian(共轭转置)矩阵，即

y(n)=∑Mm=-MKp(m,n)Yp(m)(5.34)

联立式(5.32)和式(5.34)得到


y(n)=∑Mm=-M∑Nk=-NKp(m,n)Kpm,ky(k)

=Δt22πsinα∑Mm=-M∑Nk=-Nej2cotαk2-n2Δt2ejm(n-k)ΔtΔucscαy(k)
(5.35)

为了使上式中对m的求和等于δ(n-k)，即

∑Mm=-Mej
m(n-k)ΔtΔucscα=δ(n-k)(5.36)

那么，需要满足

ΔuΔt=
2πSsinα2M+1(5.37)

其中，|S|是与2M+1互为质数的整数。这样，式(5.50)变为

Kp(m,n)=1-jcotα2πΔte
j2
m2Δu2cotαe-j
2πnmS2M+1ej2n2Δt2
cotα(5.38)

这样，得到


∑Mm=-M∑Nk=-NKp(m,n)Kpm,ky(k)

=2M+12πsinαΔt2y(n)

=2M+12πsgnsinαsinαΔt2y(n)
(5.39)

对Kp(m,n)做归一化处理以满足式(5.35)，于是得到变换矩阵Kp(m,n)为

Kp(m,n)=sgnsinαsinα-jcosα2M+1ej2

m2Δu2cotαe-j
2πnmS2M+1ej2
n2Δt2cotα(5.40)

为简便起见，选择S=sgnsinα=±1，式(5.40)改写为

Kp(m,n)=sinα-jsgnsinαcosα2M+1
ej2
m2Δu2cotα
e-j2πnmsgn(sinα)2M+1ej2
n2Δt2cotα(5.41)

于是，我们针对sinα＞0和sinα＜0得到以下两个DFRFT公式。
(1) sinα＞0，即α∈2Dπ+0,π： 

Yp(m)=sinα-jcosα2M+1
ej2
m2Δu2cotα
∑Nn=-Ne-j2πnm2M+1
ej2n2Δt2cotαy(n)(5.42)


(2) sinα＜0，α∈2Dπ+-π,0： 

Yp(m)=-sinα+jcosα2M+1
ej2m2Δu2cotα
∑Nn=-Nej2πnm2M+1
ej2n2Δt2cotαy(n)(5.43)

另外，必须满足限制条件M≥N和

ΔtΔu=2πsinα2M+1(5.44)

可以看到，当M=N且α=π/2时，式(5.42)简化为DFT； 当α=-π/2时，式(5.43)简化为IDFT。我们也可以看到，当α=Dπ，D为整数时，Δt和Δu没有合适的选择。即当α=Dπ时，不能用式(5.42)和式(5.43)定义DFRFT。事实上，这种情况可以用下式来定义


Yp(m)=y(m),α=2Dπ(5.45)

Yp(m)=y-m,α=2D+1π(5.46)


我们注意到，在式(5.44)中，如果sinα很小，Δt和Δu也必须很小，采样点数将增加，这将增加DFRFT的计算量。因为对于连续FRFT来说，有

Xp(u)=Fp-1F1［x］(t)(u)(5.47)

所以，当sinα很小时，我们可以先作x(t)采样信号的DFT，然后再计算阶次为p-1的DFRFT。因此，我们将上面的DFRFT变为


Yα(m)=Ce-j2m2Δu2tanα∑Nr=-N∑Nn=-Nej2π
sgn(cosα)rm2M+1×


e-j2r2Δf2tanα
e-j2πnr2N+1y(n)
(5.48)

其中


ΔuΔf=2π|cosα|2M+1(5.49)

C=|cosα|+jsgn(cosα)sinα(2M+1)(2N+1)(5.50)


因为ΔtΔf=2π2N+1，Δf=2πΔt2N+1，所以在sinα≈0时，我们可以定义修正的DFRFT如下： 


Yp(m)=Ce-j2m2Δu2tanα∑Nr=-N∑Nn=-Ne
j2πrmsgn(cosα)2M+1·

e-j2π2r2tanα
(2N+1)2·Δt2e-j2πnr2N+1y(n)

(5.51)

其中，Δu=2N+1cosαΔt2M+1。另外，应当注意到，在上面DFRFT推导中已经进行了归一化，因此在利用上述的DFRFT计算连续FRFT结果时必须考虑这个归一化因子。
CFDFRFT通过对采样间隔的合理限制，使其具备可逆性。在计算复杂度方面，CFDFRFT 包含两个chirp 乘法和一个FFT 运算。因此，它的总运算量为2P+P2log2P，其中P=2M+1为输出序列的长度。相较于IPDFRFT，它具有更小的计算复杂度。总的来说，CFDFRFT 满足理想DFRFT准则中的(1)，(2)，(3)和(5)。对于旋转相加性来说，虽然它不满足旋转相加性，但通过一定转换能够从一个分数阶傅里叶域得到另一个分数阶傅里叶域的结果，具体内容可参见文献［4］。
总结说来，采样型DFRFT都满足理想DFRFT准则中的(1)，(2)和(5)。其中，CFDFRFT还满足准则(3)。因此，对于采样型DFRFT，当我们只是为了利用离散变换去计算连续FRFT时是非常有用的。这种DFRFT把原始函数的N个采样值映射为FRFT的N个采样值。这种形式的DFRFT具有很好地逼近连续FRFT的精度，并且可以利用FFT获得运算量为O(NlogN)的数值算法。在一些应用中，我们只是希望DFRFT可以很好地逼近连续FRFT，而不利用旋转相加性，这时通过把连续FRFT用这种DFRFT取代，在连续分数阶傅里叶域推导的各种信号处理算法可以直接应用到离散信号处理上。此外，由于这种DFRFT具有闭合形式，所以在一些应用中有利于推导一些性质。由于采样型DFRFT的这些优点，它广泛应用于分数阶傅里叶域非均匀采样和重构、chirp信号检测和参数估计、分数阶傅里叶域滤波、分数阶傅里叶域多采样率理论等。



教学视频



5.1.3稀疏离散分数阶傅里叶变换
与传统的DFT相比，DFRFT多了一个参数自由度，分数阶次未知，这就需要调频率的搜索，也需要大量的去斜(Dechirp)操作以及FFT运算，导致其计算量高出很多，这成为工程实际应用上的难题。为了降低DFRFT的计算复杂度，工程上采用了很多DFRFT的数值计算方法，其中Pei采样型DFRFT是常使用的降低计算复杂度的算法之一。然而，对于大数据量的输入信号，在应用Pei采样型DFRFT算法之后，依然难以解决计算复杂度过高的问题。
针对这一问题，工程上常使用分段式处理的方法［78］。该方法的思路是，将整个运算模块分成一个个微小的计算单元，分布式并行处理，最后将所得的一个个微小单元整合出最后的结果。但是，从整体来看，不仅总的时间资源并没有节约，另外的划分与组合的操作过程反而增加额外计算量，而且这种方法会降低频谱分辨率。针对这一难题，一种名为ZoomFFT的算法［9］被提出，保持了频谱的高分辨率，但该算法会导致频域过窄的观测范围。为此，在文献［10］中寻求利用信号的稀疏特性的优势，在保证良好分辨率的同时，降低了复杂度。即使该算法需要关于信号稀疏性的先验信息，但依靠信号的稀疏特性，能在提高算法效率和保持高分辨率之间取得了良好的平衡。稀疏傅里叶变换在雷达信号处理领域取得了大量的成功应用［1113］，利用信号稀疏性的优势为本章节的研究指明了方向。

对于常见的信号，具备稀疏性是普遍情况，例如语音信号和图像信号在小波变换下都是稀疏的。分数阶傅里叶域估计目标参数能更好地聚集加速目标回波或一般的线性调频信号能量［1415］，因此刘升恒等提出了稀疏分数阶傅里叶变换(Sparse Fractional Fourier Transform,SFRFT)和对应的线性调频信号检测算法［1617］。值得注意的是，文献［18］提出一种新型的简单实用型的稀疏傅里叶变换方法，因其算法的高效和实用得到了广泛的关注，也是很多改进的稀疏傅里叶算法的基础。对于频域稀疏的大数据量信号，该方法可将DFT的计算量降低至O(log2NKNlog2N)。我们可以将这种方法的思路引入Pei采样型算法，设计新型离散SFRFT算法，大幅降低运算的时间成本。另外,根据信号频域大值分量的数目是否已知，稀疏傅里叶变换算法可被分成确定型和随机型算法。在实际工程应用中，分数谱大值分量的数目很难提前获取，因此我们更倾向于设计随机型的稀疏分数谱估计算法。
5.1.3.1随机型稀疏分数谱估计算法
对于频域大值分量未知且分数域存在稀疏性的非平稳信号，基于文献［11］提出的新型稀疏傅里叶变换算法，可优化Pei采样型算法的处理流程，大幅降低计算复杂度。输入信号的长度越长，我们设计的算法在计算效率方面的优势越突出，尤其对于信号长度超过212的数据有明显的优势，而这样的数据长度在分数域滤波、雷达SAR成像、数据加密及压缩等工程应用中是十分常见的。本节先介绍此随机型SFRFT算法的结构与流程，然后介绍分数旋转角的选取方法，接着分析所提的随机型稀疏分数谱估计算法与Pei采样型算法和稀疏傅里叶算法的关系，最后仿真分析算法性能。
1. 算法流程
一般的稀疏傅里叶算法流程大致上分成采样、定位和估值这三个过程，我们依次按照这三个过程的思路，具体化所提出的随机型SFRFT算法，将其分成9个步骤。
(1) 对输入信号做chirp乘积。
假定原始输入信号f(n)是非周期信号，且在分数域稀疏，满足Dirichlet条件。为了减小载波中chirp基的影响，对原始输入信号f(n)和chirp做乘积

x(n)=f(n)ej2cotαn2Δt2,n∈［1,N］(5.52)

其中，α是DFRFT的旋转角； Δt是输入信号的采样间隔。
(2) 频谱重排。
频谱重排有两个作用： ①以极大的概率分离出大值系数［19］； ②附加相位信息。因此可以在每次的随机循环中打乱频谱邻近点间的关联，并分隔相邻的谱系数。但是，我们已知时域信号，而其频谱只有在对其作DFT后才已知，因此需要对时域信号重排以达到频域重排的效果。为了解释这个过程，我们需要介绍一些定义，并研究时频重排与频域的重排之间的关系。
定义5.1(取模运算)： 给定两个正整数a和b，定义式amod b表示的是a欧几里得除法除以b得到的余数。
定义5.2(模逆)： 若

σ-1s.t.(σ×σ-1)modN=1(5.53)

则σ∈［1,N］对于modN可逆，定义σ-1是σ关于模N的模逆。
定义5.3(时域重排)： 假定σ∈［1,N］为重排因子，与N互质，且modN可逆。为了避免多次重排间有关联，σ优先选择随机质数，至少要为随机的奇数。对于时域信号x(n)，定义重排操作为n→σnmodN，则重排后时域信号的数学表达式为

s(n)=x((σn)modN)，n∈［1,N］(5.54)

定理5.1(频谱重排)： 时域信号x(n)使用重排因子σ重排后，相应的频域信号会根据重排因子σ-1重排，则时域重排信号s(n)和时域信号x(n)两者频域的关系为

S(m)=X((σ-1m)modN)，m∈［1,N］(5.55)

证明： 定义符号右下的角坐标都隐含取模运算，如对长度为N的信号x，其包含角坐标的标注xn指x(nmodN)。则可简化DFT的定义： 

Xm=∑Nn=1ωnmxn，m=0,1,…,N-1(5.56)

则对任意m∈［1,N］有

Sm=DFT{sn}=DFT{xσn}=∑Nn=1ωnmxσn=∑Nn～=1ωσ-1n~mxn~=xσ-1m(5.57)

(3) 加窗。
对重排后的时域信号加窗，以平滑地提取部分信号并减少频谱泄漏。定义平坦窗函数g(n),n∈［1,N］，窗长为w，其频域信号G(m)满足


G(m)∈［1-δ,1+δ］，m∈［-ε′Nε′N］

［0,δ］,m［-εNεN］(5.58)

其中，ε′和ε分别表示通带和阻带的截断因子，δ表示波纹振荡程度。对重排后的时域信号加窗后的信号y(n)=g(n)s(n),n∈［1,N］，则y(n)的支撑集满足supp(y)supp(g)=-w2,w2。
(4) 时域混叠和频域的关系。
当sinα＞0时，假定存在一个正整数B能整除N，则可构造信号


Z(m)=FFT{z(n)}=FFT∑w/B-1i=0y(n+iB),n∈0,B-1(5.59)

当sinα＜0时，则将上式的FFT替换成IFFT。根据文献［18］，可得到频域信号Z(m)和Y(m)的关系

Z(m)=Y(mN/B)，m∈［0,B-1］(5.60)

因此，可以得出结论，时域信号的混叠对应其频域信号的子采样。同理，时域信号的子采样对应频域信号的混叠。针对上述的时域重排和频域子采样的过程，我们可以取一个例子以便于理解，如图5.5所示。


图5.5时域重排和频域子采样对应关系的例子




(5) 哈希映射。
定义一个哈希函数

hσ(m)=((σm)mod N)B/N(5.61)

和一个偏移函数

oσ(m)=(σm)modN-hσ(m)N/B(5.62)

(6) 定位循环。
定义一个集合

J =argmaxmZ(m)(5.63)

假定该集合中含有Z(m)的2k个较大坐标值。由哈希函数，输出原像： 

I={m∈［1,N］|hσ(m)∈J}(5.64)

原像集合I中的元素个数为2kN/B。定位循环认为，改变重排因子，根据像到原像的映射关系，进行多次定位循环，如果存在一些原像的位置以很大的概率保持不变，则其对应真实大值的概率就更大。如图5.6所示，我们选择不同的重排因子，画出像到原像的对应关系，反映出这种规律。


图5.6三次重排后定位循环的坐标映射示意图



(a) 参数设置： N=16,B=8,w=12,k=1,σ=3(σ-1=11)；  
(b) 参数设置： N=16,B=8,w=12,k=1,σ=5(σ-1=13)； 

(c) 参数设置： N=16,B=8,w=12,k=1,σ=7(σ-1=7)






(7) 估值循环。
X(m)中k个大值的估计值可由下述表达式计算得到

X^(m)=Z(hσ(m))e-jπoσ(m)w/NG(oσ(m)),m∈I

0,m∈［1,N］∩I-(5.65)

假定l是估值循环次数，ltotal表示总的循环次数且ltotal=lloc+lest，正整数lloc和lest分别是定位循环和估值循环次数。当l≤lloc时，执行上述步骤(2)~(6)； 当lloc≤l≤ltotal时，执行上述步骤(2)~(7)； 当l＞ltotal时，循环结束。
(8) 计算中值。
采用对实部和虚部分别取中值的方法，输出X(m)的估计值： 

X～(m)=Medianlest{R{X^(m)}}+jMedianlest{I{X^(m)}(5.66)

(9) 频域调制。
为了将信号从傅里叶域调制到分数阶傅里叶域，我们将式(5.66)乘上指数chirp函数，最终输出结果： 

F^α(m)=X～(m)ejm2Δu22tanα(sinα-jcosα)sgn(sinα)/M(5.67)

式中，Δu表示输出信号的采样间隔，M表示频域信号的长度。对于sinα＜0的情况(当sinα＞0时，FFT代替IFFT)，算法流程如图5.7所示。


图5.7当sinα＜0时SFRFT的算法流程




2. 分数旋转角选取
在一些应用领域中，分数旋转角α的值是已知的，比如在SAR成像和线性调频匹配滤波中。然而在其他的大多数情况，分数旋转角α的值是未知的。最大似然估计(Maximum Likelihood Estimation,MLE)是估计分数旋转角α的常用方法，例如在文献［20］中，利用最大似然估计的离散形式，提出离散调频傅里叶变换的方法，准确估计出旋转角α的值。但是，最大似然函数的二维最大化过程需要巨大的计算量，时间消耗的成本过高，故需要找到可替代的方法。文献［21］中提出离散多项式相位变换的方法，该方法将chirp信号转化成正弦波进而检测频谱大值量的位置，可快速得到分数旋转角α的估计值，避免了MLE法二维函数最大化过程的时间消耗。
下面介绍基于多项式相位变换估计分数旋转角α的具体过程。假设x(n)是一个关于离散实值变量n的复值函数，τ为时延量。定义算子DP1［x(n),τ］和DP2［x(n),τ］的数学形式： 


DP1［x(n),τ］=x(n)(5.68)

DP2［x(n),τ］=x(n)x*(n-τ)(5.69)


算子DP2表示对离散信号x(n)取差分，相当于对其降阶一次［21］。再假定一个算子DPT表示算子DP的DFT，则DPT2表示DP2的DFT。对于一个离散时域信号s(n)=exp{jπμ(nΔt)2}，Δt表示时域采样间隔，μ为调频率，使用算子DPT2后有如下表达式

DPT2［ejπμ(nΔt)2,ω,τ］=DFT{DP2［ejπμ(nΔt)2,τ］}=DFT{ej2πμτnΔt-jπμ(τΔt)2}(5.70)

据此表达式，DPT2［x(n)ejπμ(nΔt)2,ω,τ］的能量聚集在

ω=ω0=2πμτΔt(5.71)

在分数阶傅里叶域估计目标参数能更好地聚集加速目标回波或一般的线性调频信号能量［2021］。另外，由文献［21］，可证明当τ=N2时估计精度最高，故可由μ估计出α。
但是受到信道噪声影响，以此估计出旋转角α的精度受到限制。为了解决这个问题，需要对旋转角α估计结果附近的小范围进行精细搜索，而搜索步长Δα取决于调频率分辨率Δμ的约束。由f=μτΔt可得

Δf=ΔμτΔt(5.72)

取τ=N2，则可得到

Δμ=ΔfτΔt=4T2(5.73)

其中,T为信号时长。结合原始输入信号f(n)的表达式，可计算出搜索步长： 

cotα2=πμ(cotα)′Δα=2πΔμ|Δα|=8πsin2α/T2(5.74)

3. 与Pei采样型算法和稀疏傅里叶算法的关系
根据文献［24］，对于时域信号x(t)，有连续分数阶傅里叶的变换公式： 


{Fαx}(u)=∫+∞-∞Kα(u,t)x(t)dt,0＜|p|＜2,0＜|α|＜π

=1-jcotα2π∫+∞-∞ejt2+u22cotα-jtucscαx(t)dt,α≠Dπ

x(t),α=2Dπ

x(-t),α=(2D±1)π

(5.75)

其中，D为任意整数； Kα(u,t)是CFRFT的核函数，其中u为CFRFT的频率。假设p是CFRFT的阶次，则分数旋转角α=pπ/2。当CFRFT阶次p使得旋转角α=2Dπ+π/2时，根据式(5.75)，该CFRFT的表达式可变成传统的连续傅里叶变换公式。
再回顾一下Pei采样型算法，该算法是在连续FRFT的基础上推导出来的。该算法分别对输入的时域信号和输出的频域信号等间隔采样，采样间隔分别是Δt和Δμ，且需满足约束


ΔtΔμ=2π|sinα|M(5.76)

该约束是为了保证变换的可逆性。令输入信号长度为N，则M必须满足M≥N。当M=N时，DFRFT表达式可写成［4］

{Fαx}(m)=sinα-jcosαMejm2Δu22tanα∑N-1n=0x(n)e
j2n2Δt2cotα-j2πnmM,α∈2Dπ+(0,π)


-sinα+jcosαMejm2Δu22tanα∑N-1n=0x(n)ej2
n2Δt2cotα+j2πnmM,α∈2Dπ+(-π,0)

x(m),α=2Dπ

x(-m),α=(2D+1)π(5.77)

若α≠Dπ，一次FFT运算和输入/输出的两次chirp乘积构成了Pei采样型算法的主要计算量，只考虑复乘次数，则Pei采样型算法的计算复杂度为O2N+N2log2N。
由此可见，即使Pei采样型算法是DFRFT最高效的算法之一，但是其依然有较高的计算复杂度，尤其是输入信号长度N非常大时，FFT运算占据了很大的计算量。幸运的是，信号稀疏是普遍情况，可以采用稀疏傅里叶变换代替一般的FFT运算［18］，进一步提高DFRFT的计算效率。该稀疏傅里叶变换算法是一个滤波、定位和估值的过程。滤波器是由切比雪夫函数和矩形窗函数卷积得到，其作用是将稀疏信号的整个频域分成一个个频率单元，并使得信号在时域和频域都具有良好的聚焦效应。定位和估值的过程与sketching/streaming的算法［25］类似，这种方法的优势在于避免了传统方法所涉及的插值或迭代的操作。因此，我们可以尝试，以Pei采样型算法的流程为基础，融入稀疏傅里叶变换的算法思路，设计改进型的Pei采样型算法。
改进后的Pei采样型算法，即所提的SFRFT算法，其所适用的信号有一定的限制，即要求信号具有一定的稀疏性。在应用所提的SFRFT算法之前，需要判定稀疏度是否符合算法适用要求。判断的方法主要有两种，一种是根据先验信息直接判断； 另一种则是对信号先进行试采样，再由采样信息判断稀疏性是否适合。对于信号数据量大且非平稳，在分数阶傅里叶域稀疏的信号，其大值个数K最好远小于数据点数N，K/N的值越小则稀疏性越好。经验上一般认为稀疏度满足K/N＜1%的信号适合使用所提算法进行分析处理，而此类信号普遍存在于很多领域，如分数域滤波、合成孔径雷达成像，以及GPS定位过程中动态信号快速获取等。
4. 算法性能分析
为了分析所提SFRFT算法的性能，下面分别从分辨能力、计算复杂度和鲁棒性三个方面进行仿真分析。
在第一个仿真实验中，主要分析所提SFRFT算法在多旋转角情况下的分辨能力。在本次仿真实验中，我们设置如下的仿真参数。设置信号的采样频率fs=900Hz，数据长度为N=215，变换域的大值点数为K=5，子采样FFT的长度设为B=1024。假定输入信号有四个频率分量，其初始频率分别为100、200、300和300Hz，相应的chirp率分别为10、11.85、13.85和13.85Hz/s，并设置相应的信噪比为-12、-18、-24和-12dB。噪声服从均值为0的高斯分布。滤波器参数设置为w=22883，δ=10-6，ε′=2×10-4且ε=5×10-4。在估值和定位循环的过程中，设置估计循环参数lloc=4，定位循环参数lest=11。
本仿真实验的所有结果显示在图5.8中。图5.8(a)显示的是输入信号的频谱图。为了研究在多旋转角的情况下，SFRFT与原DFRFT频谱效果的差别，我们设置三个不同的分数旋转角α，每一行匹配不同的分量阶次，得到图5.8(b)~(g)这三行两列的效果图。仿真结果说明，我们提出的稀疏分数阶傅里叶算法在多旋转角的情况下，能够精确估计出稀疏分量的分数域的频率和幅值。另外，如果信号的分数阶傅里叶域的频率分量呈现出不聚焦或是不稀疏的特点，则其谱线的分布有连续性、低幅值的规律。图5.8(f)和图5.8(g)中紧邻谱线的局部放大图表明，即使在稀疏性不是十分理想的情况下，SFRFT依然保证了良好的多分量分辨性能。


图5.8稀疏傅里叶变换多分量分辨性能



(a) 输入信号频域波形； (b) 匹配第一个分量阶次的离散分数阶傅里叶变换； (c) 匹配第一个分量阶次的SFRFT； 

(d) 匹配第二个分量阶次的离散分数阶傅里叶变换； (e) 匹配第二个分量阶次的SFRFT； 

(f) 匹配第三个分量阶次的离散分数阶傅里叶变换； (g) 匹配第三个分量阶次的SFRFT


在第二个仿真实验中，接着对所提SFRFT算法的复杂度仿真分析。以复乘次数为衡量标准，依据算法步骤的过程分析，则所提的SFRFT算法复乘总数为


MSFRFT=2N+(ω+Blog2B/2)×lloc+

(ω+Blog2B/2+2k)×lest+card(I)×ltotal

(5.78)

其中，card(·)表示集合中元素的数目。基于该复乘运算总次数的表达式，可比较SFRFT和Pei采样型DFRFT算法的计算量。在此仿真过程中，我们设置计算大值个数为K=5，循环次数分别设置为lloc=3和lest=8。图5.9显示，与SFRFT算法相比，随着数据长度的增大，DFRFT在计算效率上的优势越大，尤其当数据长度的量级在213及以上时。


图5.9DFRFT和SFRFT的计算复杂度比较




另外，不同于DFRFT算法，SFRFT算法的计算复杂度与输入信号的稀疏度十分相关。为了研究所提SFRFT算法的计算复杂度与稀疏度的关系，给定输入信号的长度为216，以分数阶傅里叶域的非零值个数为坐标变量，其数目依次从10变化到104。相应地，随着信号稀疏度的变化，循环次数也在确保输出结果精度的情况下匹配。从图5.10可以看出，稀疏度越高，SFRFT算法相对于DFRFT算法在计算复杂度方面的优势就更明显。但是，当分数域的非零频点个数大于103后，SFRFT算法对计算效率就没有提升了，反而比DFRFT算法的计算量更大。


图5.10DFRFT和SFRFT的计算复杂度与稀疏度的关系




在第三个仿真实验中，分析所提SFRFT算法的鲁棒性。SFRFT算法在计算效率方面的优势与算法的鲁棒性并不冲突，SFRFT算法在提升计算效率的同时，依然具有良好的鲁棒性［18］。下面以仿真结果证明所提SFRFT算法对输入噪声的鲁棒性。仿真实验的参数设置如下，信号长度N=215，k=3，旋转角α=0.01rad，Monte Carlo仿真次数为20000，输入信噪比SNR变化范围设为［-10dB,30dB］。假设ε表示估计误差，该估计误差是用SFRFT的输出{SFαx}(i)和DFRFT的输出{Fαx}(i)之间的差值来衡量

ε=1k∑i∈(0,N］|{SFαx}(i)-{Fαx}(i){Fαx}(i)|(5.79)

以此估计误差为标准，信号比SNR为横坐标，得出误差变化曲线如图5.11所示。该估计误差的仿真结果证实了在噪声环境下SFRFT算法良好的鲁棒性。


图5.11SFRFT算法鲁棒性与信噪比的关系




5.1.3.2随机型算法的优化
我们提出随机型稀疏分数阶傅里叶算法利用输入信号的稀疏性，使用稀疏傅里叶变换替代FFT过程，有效地降低了计算复杂度。然而，此随机型SFRFT有如下两方面的缺陷。
(1) 当输入信号的长度较短时，所提的随机型SFRFT算法和一般的DFRFT算法相比，在计算量上并没有明显的提升。
(2) 对于噪声干扰较严重的信号，频域大值量的定位估计误差会很大，大幅降低算法的估计精度。
为此，SFRFT算法后来也陆续有优化算法提出［19,21,23,25］，其中有些算法是针对虚警率和检测概率等先验信息的优化版本［22］，并成功应用于微弱雷达目标检测［23］。我们提出优化的稀疏分数阶傅里叶变换(OSFRFT)算法，采取如下的方法应对上述的问题。
(1) 使用NeymanPearson检测处理噪声干扰的信号，实现对噪声干扰后信号的分数频域估计。
(2) 分析噪声对定位循环的影响，通过ParzenRosenblatt窗方法获得相位误差的分布，并提出定位误差的优化策略。
5.1.3.2.1简化SFRFT
在上节稀疏分数阶傅里叶算法的基础上，我们做了简化，将该算法流程分成两个核心模块，即“hashtobins”过程和定位及估计循环过程。
首先介绍“hashtobins”过程，该过程将大值的分数阶傅里叶系数映射到子采样频谱。“hashtobins”过程可以分成三个步骤： 重排、加窗和下采样。子采样因子D=N/B。重排使用三个参数σ,a,b对输入时域信号的重新排序

(Pσ,a,bx)i=xσ(i-a)ωσbi(5.80)

其中，σ∈{1,3,…,N-1}且a,b∈{0,1,…,N-1}。根据文献［27］，上述公式对应的频域次序为

(P^σ,a,bx)σ(i-b)=x^iωσai(5.81)

为了表示方便，令x^=P^σ,a,bx和q=σ(i-b)，其中x^表示信号重排后的频谱。那么，上述表达式可写成

x^q=x^iωσai(5.82)

接着加窗，所加的窗函数与式(5.58)相同，是长度为D的平坦窗函数，则上述信号可写成

y^i=∑i+D/2-1q=i-D/2x^q(5.83)

下一步是下采样，即信号在时域以参数D混叠，则相应的频域为D倍下采样，得到下采样频谱

u^j=y^jD=∑jD+D/2-1q=jD-D/2x^q(5.84)

由上式可知，在下采样后已重排的谱x^q中的系数将映射到u^j，其中q和j的关系为

j=round(q/D)(5.85)

在下采样执行后，继续对其作B点的FFT，所得到的B点频谱是此“hashtobins”过程的输出。值得注意的是，如果两个及以上的大值系数都映射到u^中的同一个位置，则大值系数的定位和估值将会发生冲突，应当尽量避免。
对于时域信号混叠后的B点频域信号，每一个大值系数对应原频域中的D个位置。Hassanieh所提的稀疏傅里叶变换算法［27］，多次改变随机重排参数后执行定位循环，以在D个候选解中决定大值系数在原频谱中的真实位置。对于精确已知的k稀疏信号，Hassanieh的算法可基于文献［18］设计，但是其仅使用两次“hashtobins”过程定位大值系数。第一次“hashtobins”过程中，信号不附带相位偏移，而第二次“hashtobins”过程中，信号被附加一个线性相位偏移项。接着，该算法利用这两次“hashtobins”过程之间的相位差来定位大值系数。该算法由多次迭代组成，在每次迭代中，都有部分的大值量被估计出来，然后将此时的信号减去已估计出来的大值系数所对应的部分，再将其放入下一次循环。这个减法操作在信号的频域完成，以减少每次迭代的复杂度。Hassanieh在文献［27］中提出的另一种算法有类似的思路，但区别在于，其可用于一般信号的定位和估值。该算法设置候选集合，执行多个定位循环，在每次迭代后，候选的可能大值集合都会被缩小。
5.1.3.2.2NeymanPearson检测
假定输入信号是一个线性调频信号，其幅度是A。信号的噪声n服从自适应复高斯噪声n~CN(0,σ2t)。在DFRFT的过程中，当设置合适的旋转角α后，信号的分数阶傅里叶频谱是稀疏的。此时，频域信号x^由复指数信号频谱和噪声频谱n^组成，其中n^~CN(0,σ2f)，σf=Nσt。令n^i(i∈{N})表示噪声频谱的N个采样点。假定il表示信号频域中大值系数的位置，Af是其相应的幅值，则指数信号的频谱可表示成

s^i=Af,i=il

s^i=0,其他
(5.86)

NeymanPearson检测的阈值设为ζ。为了获此阈值，需要分析“hashtobins”过程中的每一步以确定噪声分布。重排过程将原频谱中的大值系数x^il哈希映射到x^ql，映射关系ql=σ(il-b)。加窗后的信号y^i中的每一个位置的数都是x^q中的D点之和。在下采样后，信号x^ql映射到u^jl，并有映射关系jl=round(ql/D)。由于重排后信号的频谱x^ql服从分布x^ql~Af+CN(0,σ2f)，其余的B-1个点都服从均值为0且方差为σ2f的高斯噪声分布，因此可得到u^jl服从的分布

u^jl~Af+CN(0,σ2u)(5.87)

其中，σu=Dσf=DNσt。原信号频谱x^的信噪比是A2f/(2σ2f)，在“hashtobins”过程后该信噪比为原来的1/D。如果没有大值系数落入u^的第j个位置，则u^j的分布是个服从独立同分布的复高斯噪声之和，即

u^j~CN(0,σ2u)(5.88)

因此，当j≠jl时，u^j的幅度服从Rayleigh分布。如果有大值量落入u^的第j个位置，即j=jl，u^j由幅度为Af的信号分量和方差为σ2u的复高斯噪声组成，故u^j的幅度服从Rice分布。由此，可总结出u^的幅度服从以下分布

fρ(ρ)=ρσ2uexp-ρ2+A2f2σ2uI0ρ|Af|σ2u,j=jl

fρ(ρ)=ρσ2uexp-ρ22σ2u,j≠jl(5.89)

其中，ρ是u^j的幅值，I0(z)表示第一类零阶修正贝塞尔函数［28］。
fρ(ρ)的右尾概率可由下式获得

Pr{ρ＞ζ}=Pd=Q1Afσu,ζσu,j=jl

Pr{ρ＞ζ}=Pfa=exp-ζ22σ2u,j≠jl(5.90)

其中，Q1(z)是Marcum Q函数［28］； Pd是检测率； Pfa是虚警率。
可根据期望的Pfa求出阈值ζ，故由上式可得

ζ=σuTn(5.91)

其中，Tn=-2lnPfa是由期望的Pfa决定的名义因子。此时可估计x中的噪声等级，获得噪声方差σu，并接着由上式选择阈值ζ。因此，在使用NeymanPearson检测后，我们提出的SFRFT算法，则不需要精确稀疏性的先验知识。
5.1.3.2.3噪声对定位循环的影响及优化策略
1. 噪声对定位循环的影响
在此部分，我们从两方面对定位循环中的噪声影响进行分析，第一方面是分析相位误差和噪声等级之间的关系，另一方面是分析相位误差对估计结果的影响。
第一部分我们先分析定位循环中相位误差和噪声等级之间的关系。为了建立下采样频谱u^中的相位误差和原谱x^中的噪声之间的关系，我们需要先用x^表示u^。根据重排特性可得
到如下关系： 

x^σ(i-b)=x^iωσai(5.92)

重写上式为

x^q=x^(σ-1q+b)ωa(q+σb)(5.93)

其中，q=σ(i-b)。基于上述表达式，我们分别分析两次“hashtobins”过程中x^和u^的关系。在第一次“hashtobins”过程中，令a=0且ωa(q+σb)=1。基于式(5.93)，则式(5.84)可写成

u^j=∑jD+D/2-1q=jD-D/2x^q=∑jD+D/2-1q=jD-D/2x^(σ-1q+b)(5.94)

当j=jl时，在此求和过程中存在大值系数，则u^jl可进一步表示成

u^jl=∑jlD+D/2-1q=jlD-D/2ν^(σ-1q+b)+Af(5.95)

在第二次“hashtobins”过程中，a=1，因此噪声频谱和指数信号频谱都应当乘上一个相位项，即

u^jl=∑jlD+D/2-1q=jlD-D/2ν^(σ-1q+b)ω(q+σb)+Afω(ql+σb)(5.96)

其中，ql=σ(il-b)。由ω(ql+σb)=ωσil 得

u^jl=ωσil∑jlD+D/2-1q=jlD-D/2ν^(σ-1q+b)ω(q-ql)+Af(5.97)

u^jl和u^jl间的相位差是

(u^jl/u^jl)=2πσilN+∑jlD+D/2-1q=jlD-D/2ν^(σ-1q+b)+Af∑jlD+D/2-1q=jlD-D/2ν^(σ-1q+b)ω(q-ql)+Af(5.98)

其中，q∈{jlD-D/2,jlD-D/2+1,…,jlD+D/2-1}，这表明一旦输入信号x确定，则输入噪声谱的D个样本也确定了。在两次“hashtobins”过程中，σ和b被选择为相同的值。落入u^jl和u^jl的噪声点数一样，差别在于由a导致的相位项。令ν^m=ν^(σ-1q+b)，其中m∈［-D/2,D/2-1］，则上式可变成

(u^jl/u^jl)=2πσilN+∑D/2-1m=-D/2ν^m+Af∑D/2-1m=-D/2ν^mω(m+jlD-ql)+Af(5.99)

其中，jlD-ql∈［-D/2,D/2-1］表示u^在第jl个位置上大值系数ql的偏移量。因此，可得相位差为

err=∑D/2-1m=-D/2ν^m+Af∑D/2-1m=-D/2ν^mω(m+jlD-ql)+Af(5.100)

err可视作ω(m+jlD-ql)的加权平均值。当ω的指数为m时，加权系数为ν^m； 当ω的指数为0时，加权系数为|ν^-(jlD-ql)+Af|。
虽然ω(m+jlD-ql)和Af被固定了，但{ν^-D/2,ν^-D/2+1,…,ν^D/2-1}是D个随机变量，故err也是一个随机变量。相位误差err的概率密度函数(PDF)f(err)需要进一步分析。然而，上式中随机变量的组合是复杂的，因此得到f(err)准确的表达式是比较困难的。
为此，可以采用一种非参数方法估计随机变量的概率密度函数，即用核密度估计法来估计概率密度函数f(err)，称作ParzenRosenblatt窗方法。首先，我们使用Monte Carlo法形成ν^m的样本集，其中ν^m~CN(0,σ2f)。接着，我们使用式(5.100)得到样本集{err1,err2,…,errξ}，其中ξ表示采样的总数目。集合中的所有元素都服从同一个未知的分布f(err)并相互独立。因此，可使用如下表达式估计f(err)

fest(err)=1ξ×h∑ξi=1K(err-erri)(5.101)

其中，K是核函数，h是一个平滑参数。核函数K可选为一个标准窗，即

K(err-erri)=12πexp-(err-erri)22(5.102)

我们在本节的第二部分，分析相位误差对估计结果的影响。对于Hassanieh等在文献［27］中提出的精确“ksparsity”信号的算法，大值系数的索引由(u^jl/u^jl)估计而来，根据下述关系： 

iest=σ-1(round((u^jl/u^jl)N/2π))modN(5.103)

对于无噪声的情况，则有

ν^m=0,m∈［-D/2,D/2-1］(5.104)

此时err=0，式(5.99)可写成

(u^jl/u^jl)=2πσil/N(5.105)

因此，式(5.101)获得的iest和真实的索引il一样。
如果信号被噪声干扰，则有(u^jl/u^jl)=2πσil/N+err。当相位误差在一个确定的范围，即|err|＜(0.5+Γ)2π/N，其中Γ是正整数，σiest也将在一个确定范围，即

σiest=round((u^jl/u^jl)N/2π)∈［σil-Γ,σil+Γ］(5.106)

则有

iest∈{il-σ-1Γ,il-σ-1(Γ-1),…,il+σ-1Γ}(5.107) 

当Γ=0时，|err|＜π/N，这意味着定位循环中没有误差。注意到round((u^jl/u^jl)N/2π)和iest之间的映射是非线性的，则［σil-Γ,σil+Γ］中的指数在映射后将被分离，导致估计误差被放大。如果增加一倍，iest中的误差相应会扩大到σ-1或N-σ-1。若误差存在，iest将会偏离真实的指数il。因此，算法估计性能的精确性取决于误差校正是否有效。
2. 减小定位循环误差
为了解决定位循环中的噪声影响，我们设计了名为“LOCCORR”的定位误差校正架构，其伪代码被概括在算法 1中。该算法由两个阶段组成，在第一阶段中，可由两个定位循环的相位差形成大值系数的候选集合; 在第二阶段中，第三个定位循环被执行来找出候选集合所包含大值系数的真实位置。下面具体介绍“LOCCORR”的两个阶段。






算法1定位误差校正，命名程序为LOCCORR

输入:  u^,u^,x,ξ,J

输出: ω^

1: 由式(5.109)选择Γ;

2: 对于r依次取集合J中的元素执行{

a←u^r/u^r ;

iestr←round(a)N2π;

Icandr= {iestr-σ-1Γ,…,iestr+σ-1Γ};

}

3: Icand=Icand1∪Icand2 ∪…∪IcandR;

4: 对于σ依次取集合{1,3,…，N/2-1}中的元素执行{参数搜索

cnt=1;

再对于i依次取集合Icand中的元素执行{

Ucand(cnt) = roundBN［(σi) mod N］;

cnt=cnt+1;

}

如果Ucand不包括重复的索引，则终止循环。

}


5: u^←hashtobins(x,σ,0,0,B,G); 第三次定位循环

6: J= {m:u^m＞ζ}; 第二次检测阶段

7: ω^←0;

8: 对于u依次取集合Ucand 中的元素执行{

对于m依次取集合J中的元素执行{

如果u=m则{

i = revmap(u);

j = round BN［(σ(i-b)) mod N］;

ω^i=u^j;

}

}

}

9:返回ω^。


首先，介绍“LOCCORR”的第一个阶段。假定u^中有R个系数在阈值以上，我们先由式(5.103)获得每个系数的位置指数，并将其记为{iest1,…,iestr,…,iestR}。对于x中的每个iestr，我们都设置一个集合Icandr，其所包含的元素如下

Icandr={iestr-σ-1Γ,iestr-σ-1(Γ-1),…,iestr+σ-1Γ}(5.108)

正确的坐标指数落入集合Icandr的概率记为PΓ，其中

PΓ=Pr{|err|＜2π(0.5+Γ)/N}(5.109)

PΓ对OSFRFT算法的重构精度有重要作用，故应当合适选取Γ以维持一个确定的概率PΓ。候选集Icand是所有集合的总集，即Icand=Icand1∪Icand2…∪IcandR。
然后，介绍“LOCCORR”的第二个阶段。为了确保Icand中每个位置索引在下采样后落入不同的“bins”，我们搜索第三次定位循环中合适的重排参数σ。u^中的候选集定义为Ucand，故Icand和Ucand之间有一一对应的映射关系。在重排后，索引i∈Icand将被映射到i=(σi)modN。根据式(5.85)，在加窗和下采样过后，该指数将落入u=round(i/D)。而从Ucand到Icand的逆映射关系可由下式定义： 

revmap(u)={i:i∈Icand,round{D［(σi)modN］}=u}(5.110)

执行第三次定位循环获得u^中的大值系数，记作J。Ucand和J的交集是原频谱x^中大值系数的真实位置的索引。最后，u^j的值被分配到x^中的索引i=revmap(u)，其中i和j之间的关系为

j=roundBN［(σ(i-b))modN］(5.111)

5.1.3.2.4算法流程及性能分析
经过频谱重排、NeymanPearson检测和定位误差校正，所提出的OSFRFT的整个算法流程显示在算法2中。





算法2OSFRFT算法


输入: 稀疏信号s;旋转角α

输出: 分数阶傅里叶域频谱F^


1: 算法进程OSFRFT(s,a):

xi=siIphase,i∈{N};  第一阶段

x^=OSFRFTInner(x);  第二阶段

F^m=x^mQphase,m∈{N};  第三阶段

2: 算法进程OSFRFTInner(x):

由噪声等级和稀疏度选择B;

接着选取阈值ζ ;

从奇数集合{1,3,…,N-1}随机均匀选取σ;

从集合{1,2,…,N-1}随机均匀选取b;

选择极小值δ;

选择合适的窗长L; .

由窗函数公式选择合适的参数(0.5D,0.6D,δ,L);

u^←hashtobins(x,σ,0,b,B,G);  第一次定位循环

u^←hashtobins(x,σ,1,b,B,G);  第二次定位循环

J={j:|u^j|＞ζ};  第一次检测阶段

x^←loccorr(J,u^,u^,ζ);

返回x^。

3: 算法进程hashtobins(x,σ,0,b,B,G):

计算y^jD,j=0,…,B-1,其中:

y=G×Pσ,a,b;

u^j=y^jD;

返回u^。

1. 计算复杂度分析
OSFRFT算法的计算复杂度由算法流程中复数乘法运算次数评估，为此可以从OSFRFT算法两个阶段的计算复杂度分别考虑。在OSFRFT算法的第一个阶段中，输入信号需要与Iphase相乘，需要N次复乘。OSFRFT算法的第二阶段包括了三次定位循环，在每次循环中加窗操作和B点FFT操作分别需要L和Blog2B/2次复乘。但是，在x^中仅有k个非零系数，信号与Qphase相乘，因此在第三阶段仅需要k次复乘。结合这三个阶段的复乘次数，OSFRFT算法的总复乘次数为 

MOSFRFT=N+3(L+Blog2B/2)+k(5.112)

2. 精度评估参数
我们使用复原率来量化真实值和OSFRFT算法所得解之间的差距。复原率指的是，OSFRFT算法输出的大值系数复原在正确位置上的概率，并将这个概率即复原率记作Pr。Pr可由如下公式定义

Pr=P2dPΓ(5.113)

其中，Pd和PΓ分别由式(5.90)和式(5.109)给定。
3. B和Γ的选择
为了获取适合的复原率Pr，我们设计B和Γ的选取规则，其伪代码概括在算法 3中。在算法的流程中，根据噪声等级(A和σ)调整B和Γ来维持Pd和PΓ足够高，以便于获取最小可接受的复原率Pr。同时，约束B＞(2Γ+1)×k需要被满足，这是为了保证Icand中的系数可在子采样频谱中分成B个索引指数。






算法3B和Γ的选取规则


输入: N,A,σ2t,k,Pr

输出: B和Γ

1: 对于B依次取集合{2,4,8,…,N/2}中的元素,由式(5.190)计算Pd;

2: B={B: Pd＞Pr};

3: 对B排序,得到排序后的次序以及排序后的集合B;

4: 对于集合B中的每一个元素B执行{

由KDE估计f(err);

对于Γ依次取集合{1,2,…,
B/(2k)}中的元素执行{

如果PΓ＞Pr/P2d,则终止整个循环;

{

如果B＞(2Γ+1)×k,则终止整个循环;

}

5: 返回B和Γ。

5.1.3.2.5仿真分析
1. 相位误差分析
f(err)受参数B和信噪比SNR影响，随后数值仿真分析这两个参数对相位误差的影响程度。设置参数： N=4096，B=256/512，ξ=10000，jlD-ql=0。输入信噪比SNR分别设置成-8.1648dB和-2.1442dB。如图5.12所示，改变参数B或输入信噪比SNR，可画出其对应的概率分布直方图，图中的蓝色柱状为概率密度分布，红色线条是由核密度估计出的PDF。


图5.12相位误差的概率密度分布



(a) B=256，SNR=-8.1648dB； (b) B=256,SNR=-2.1442dB； (c) B=512,SNR=-8.1648dB



2. 调频率估计误差的影响
为了研究调频率估计误差的影响，我们先假定信号由噪声干扰的线性调频信号组成，可写成

si=Aexpj2πf0ifs+0.5μifs2+ni(5.114)

其中，i=-N2,…,N2-1; f0是线性调频信号的中心频率； fs是采样频率； n表示复高斯噪声。文献［18］提出了一种调频率估计算法，即离散多项式相位变换算法，可估计出本信号模型中的调频率，记估计值为μest。为了研究估计的调频率μest与真实调频率μ之间的差距，我们可以记这个差值为μerr。另外，OSFRFT算法的旋转角α=arccot(-μest)。
在OSFRFT算法的第一个阶段，信号s需要乘上二次指数项exp(j(cotα)i2Δt2/2)。这样，我们可得到一个中心频率为f0，调频率为μerr的线性调频信号xi，其可写成

xi=Aexpj2πf0ifs+0.5μerrifs2+ni(5.115)

其中，i=-N2,…,N2-1。对于更大的μerr，频谱幅度Af会减小，这会导致检测概率Pd和定位精度PΓ减小，进而使得复原率Pr降低。
在本节的仿真实验中，OSFRFT算法被用来恢复线性调频分量，并由Monte Carlo仿真估计不同μerr下的复原率。假定信号由一个线性调频分量组成，其中心频率f0=101.5Hz，调频率μ0=-15Hz/s。信号长度为212，采样频率为1000Hz/s，信噪比SNR=0dB。我们测试不同的误差调频率μerr，将其分别取为0.1，0.15，0.2，0.25和0.3Hz/s。
对于每个不同的误差调频率μerr，取不同B和Γ的OSFRFT算法的复原率经由Monte Carlo仿真获取。对于每组的参数设置，我们执行3000次OSFRFT算法循环，使用下式估计复原率Pr： 

P^r=Nlkη(5.116)

其中，k=1是大值系数的数目； η=3000是循环总数； Nl表示频率在101.5Hz处，可正确恢复的大值系数的数目。理论上的复原率Pr可由式(5.113)计算而来。
图5.13中的两幅仿真图分别显示了复原率Pr与μerr、复原率Pr与信噪比SNR之间的关系。图5.13(a)展示了不同情况下的Pr曲线，这与由式(5.113)计算得到的理论值一致。复原率Pr随着μerr的增加而减小，这是由于频谱幅度的减小。可采用增加参数B的数值或采用定位误差校正的方法，补偿Pr减小所造成的影响。从图中可以看出，对于例子中μerr的所有值，当B设置为256且Γ设置为2时，复原率Pr的值在97.5%以上。根据式(5.112)，当B=256，L=N/2且k=1时，OSFRFT算法的复乘次数为MOSFRFT=13313。相比之下，Pei采样型算法需要32768次复乘。由此可证明，在不同μerr值的情况下，OSFRFT算法在保持高复原率的同时，还具有更低的计算复杂度。



图5.13在不同μerr和输入信噪比下复原率Pr的变化曲线



(a) 不同的μerr； (b) 不同的输入信噪比



3. SNR的影响
本次仿真实验研究在不同SNR情况下OSFRFT算法的性能，并证明所得Pr的精度。具体地，OSFRFT算法复原线性调频信号分量，并由Monte Carlo仿真估计不同信噪比下的复原率。信号参数配置： 由一个线性调频分量组成，中心频率为101.5Hz，调频率为-15Hz/s，且信号长度为212。信号的采样频率为1000Hz/s。信号受白高斯噪声干扰，信噪比SNR分别取0，-2，-4，-6和-8。算法的参数设置： B=128256，Γ=1,2，σ=19，b=0，μerr=0Hz/s。另外，参数ξ在不同信号比下调整以保持同样的虚警率Pfa，其中Pfa=10-4。
对于所设置的每一个信噪比SNR，我们都根据式(5.116)估计Pr，估计值记作P^r。理论上的Pr依然由式(5.113)获取。图5.13(b)展示了在不同输入信噪比下Pr的变化曲线。从该图可观察到，理论上的Pr和仿真的P^r在不同的B和Γ下随SNR的变化规律具有良好的一致性，这再次证明了式(5.113)的准确性。随着SNR的降低，复原率Pr减小，为此可增加B的值或是采用定位误差校正的方法来增加复原率Pr。与5.1.3.2.4节的例子一样，设置B=256和Γ=2后，在不同的信噪比下复原率Perr可达到97.5%以上。同样，可计算出OSFRFT算法的复乘次数为MOSFRFT=13313，这低于Pei采样型算法。由此可证明，所提的OSFRFT算法在噪声干扰下保持良好性能的同时，也降低了运算复杂度。
4. 多分量分离和参数估计性能
在本节的仿真实验中，取不同的α代入OSFRFT算法，分离出多个线性调频分量，并将其结果与DFRFT算法比较。信号由4个线性调频分量组成，其中心频率分别为125、225、325和325.1Hz，相应的调频率分别是-8、-10、-12和-12Hz/s。信号的采样频率为900Hz，信号的长度为215。每个线性调频分量对应的幅度都设置为A=0.5。噪声方差取σ2t=4，则相应的信噪比为SNR=-15.0515dB。OSFRFT算法的参数设置： B=2048，σ=19，b=0，ξ=7612.2。相应的Pfa=10-6，使用定位误差校正方法，则有Γ=1。
图5.14表明，当α与一个线性调频分量的调频率一致时，其中心频率可由DFRFT/OSFRFT算法结果的峰值位置获取。由于第三个和第四个线性调频分量的调频率是一样的，第四个线性调频分量的中心频率可由匹配次序后的DFRFT/SFRFT算法的三次循环获得。在图5.14中，DFRFT和OSFRFT算法的仿真结果相互匹配得很好，这表明线性调频分量的中心频率和幅值可由OSFRFT算法精确估计。第三个和第四个线性调频分量间的间隔很小，这两个分量之间的差别只有0.1Hz，这等价于分数阶傅里叶域采样后的四个离散采样间隔。根据图5.14(c)，这两个邻近分量可由OSFRFT算法明显地区分出来，验证了所提OSFRFT算法优秀的频谱分辨率。



图5.14OSFRFT和DFRFT对4个线性调频分量的匹配顺序


图(a)~(c)中μerr的取值分别为-8，-10和-12Hz/s

(a) 第一个线性调频分量的匹配情况； (b) 第二个线性调频分量的匹配情况； 

(c) 第三个和第四个线性调频分量的匹配情况





图5.14（续）


在本节的另一个仿真实验中，我们考虑信号的分数阶傅里叶频谱除了大值系数，还具有大量的小而非零的系数的情况，即称为非精确的“ksparsity”情况。假设信号由三个线性调频分量和一个复指数分量组成，三个线性调频分量的中心频率分别是125、225和325Hz，其相应的幅度为0.5,0.75和1，相应的调频率都设置为-6Hz/s。而所含指数信号的频率设为-325Hz，其幅度为2。假定输入信号无噪声，采样率为900Hz，且输入信号的长度为215。DFRFT/SFRFT算法的旋转角α与线性调频分量的调频率相匹配，而其他参数与之前的实验设置成一样。输入信号的仿真结果，在给定的分数阶傅里叶域展现出了非精确的“ksparsity”情况。如之前的预期一样，在三个线性调频分量匹配的分数阶傅里叶域中，指数信号对应成了宽带分量。与之对比的是，在图5.15(b)~(d)中，三个线性调频分量的频点位置可由OSFRFT算法精确估计出。但由于非零小值的干扰，估计出的大值系数的幅值与真实值稍有不同。


图5.15“ksparsity”情况的数值仿真



(a) 匹配次序后含三个分量的OSFRFT算法和DFRFT算法的分数阶傅里叶频谱； 

(b)~(d) 对应于第一个、第二个、第三个线性调频分量的分数阶傅里叶频谱的局部细节





图5.15（续）


5. 计算复杂度分析
在本节的仿真实验中，我们分别比较了DFRFT算法、SDFRFT算法和OSFRFT算法的计算复杂度。在本节的第一个仿真实验中，我们以复乘次数作为计算复杂度的评估标准。DFRFT算法和SDFRFT算法的复乘次数的计算方法直接从文献［4］和［28］分别获得，而OSFRFT算法的复乘次数由式(5.112)计算得到，其中的参数B=64，参数L=N/2。
从图5.15(a)可以看出，我们所提出的OSFRFT算法的复乘次数明显少于DFRFT算法和SDFRFT算法，大概减少为原来的1/2.5。当信号的长度N≤212时，SDFRFT算法的计算复杂度比DFRFT算法更高，说明当信号长度较小时，SDFRFT算法在计算效率方面就没有了优势。与之对比的是，无论信号的长度如何，所提的OSFRFT算法在计算效率上都具有明显的优势。这是因为SDFRFT算法为了定位大值系数的位置，采用了多次“hashtobins”过程［17］。而所提的OSFRFT算法采用“ksparse”算法的方法［27］，包含两次“hashtobins”过程，另外还需一次“hashtobins”过程校正由噪声干扰导致的定位误差，总计三次“hashtobins”过程。由于“hashtobins”过程占据主要的计算量，而与SDFRFT算法相比，所提的OSFRFT算法具有明显更少的“hashtobins”过程，故OSFRFT算法的计算效率有了显著的提高。
本节的第二个仿真实验比较算法的运行时间，我们设置信号的稀疏度k=5，信噪比SNR=0dB。OSFRFT算法的其他参数设置成B=N/64，σ=67，b=0和L=N/2。本次仿真运行在具有Intel i74600U 2.1GHz GPU和8GB RAM的笔记本电脑上，对于每个N的值执行1000次实验，并取这些实验运行时间的均值。如图5.16(b)中曲线所示，与DFRFT算法和SDFRFT算法相比，所提的OSFRFT算法在运行时间上具有明显的优势，并且这种优势随着信号长度的增加而增大。



图5.16比较计算复杂度



(a) 复乘次数比较； (b) 运行时间比较





图5.17实验场景说明


5.1.3.3应用
本节我们验证OSFRFT算法在连续波雷达信号处理中的应用价值。雷达实验平台为来自Ancortek公司的SDRKIT 580B设备，该设备的模式选为连续波模式，中心频率设为fc=5.8GHz。在本实验中，我们使用一个金属立方体外形的便携电源包，将其在1.25m的高度以初始为0的速度释放。发射天线和接收天线安置在便携电源包的正上方，天线高度相对于地面为1.35m，视轴垂直向下。实验场景的几何说明展示在图5.17中。下落过程花费了0.5s，我们剪切0.09s~0.41s的下落过程用作后处理。信号采样频率fs=25600Hz，采样点数N=213。对于i=-N2,…,N2-1，基带回波信号可表示成si=Aexp(j2πτifc)+vi，其中，τi是天线和目标之间的时延，vi表示复高斯噪声。目标物体的自由落体可看作一个匀加速线性运动： τi=2c×-vmidifs-0.5ζifs2，其中vmid表示在确定的中间时刻即0.25s时目标的瞬时速度，变量ζ表示目标的加速度，c=2.9979×108m/s表示光速，则基带回波信号可写成si=Aexpj2π-2vmidfcc×ifs-ζfccifs2+vi。

目标运动的结果是产生中心频率fmid=2vmidfccHz且调频率μ=-2ζfcc的回波信号。DFRFT/OSFRFT算法可用来分析回波信号并估计目标参数。我们假设目标的加速度等同于目标的自由落体速度g=9.8m/s2，则相应的调频率μest=-2gfcc=-378.9333Hz/s。设置分数阶傅里叶变换的参数α=arccot(-μest)，算法参数B=512，σ=23和b=0。



图5.18基于连续波雷达所收集到数据的算法比较


需估计的平均信号幅度和标准噪声偏移分别设置为=160.6455和σ-t=5.7521，则输入信噪比为25.9105dB。设置Pfa1=10-6和ζ=10947，并使用一个缩短的切比雪夫窗，长度为L=N/2。我们也使用定位误差校正方法，参数Γ设置为1。FFT、DFRFT和OSFRFT算法结果展现在图5.18中。为了归一化，FFT结果图的幅度缩放N倍。与FFT的结果图相比，DFRFT的结果图更为集中。因此，α与调频率相匹配，目标能量集中在分数阶傅里叶域。根据DFRFT的结果，中心的频率为-96.88Hz，因此所估计的速度vest=-96.88×c/(2fc)=-2.5055/m/s。目标在中间时间点的速度vmid为g×0.25=-2.45m/s，这证明了算法良好的速度估计精度。OSFRFT算法在峰值附近的结果与DFRFT算法结果相接近。所有的大值系数都被检测并正确定位，幅度估计误差可忽略。因此，我们认定OSFRFT算法能够精确估计出大值频量的位置和幅度值。根据式(5.112)，k=20，因此MOSFRFT=27532。相比之下，Pei的算法涉及69632个复乘次数。这表明OSFRFT算法复杂度低于Pei的算法的1/2。此实验也证明，对于一个相对较小的N，使用OSFRFT算法，算法复杂度可大幅降低。

5.2特征分解型离散分数阶傅里叶变换
特征分解型DFRFT从连续傅里叶变换的特征函数(Hermite函数)出发，通过对Hermite函数的离散化近似和正交投影，得到一组与Hermite函数形状相似的DFT矩阵的正交化离散Hermite特征向量。然后，仿照连续FRFT的核函数谱分解表达式，构造了离散分数阶傅里叶的变换矩阵。本节将详细介绍这种DFRFT的构造方法及其所依据的相关定理。
5.2.1傅里叶变换的特征值与特征函数
首先给出傅里叶变换的特征值和特征函数。傅里叶变换的特征函数，即满足方程Ff=λf的解，是HermiteGaussian函数。记为ψn(u)，n=0,1,2，…。表达式为

ψn(u)=AnHn2πue-πu2,An=21/42nn!(5.117)

式中，Hn(u)为n阶Hermite多项式。傅里叶变换的特征值为e-inπ/2，即Fψn(u)=e-inπ/2ψn(u)，全部的HermiteGaussian函数组成一个正交集。
5.2.2离散傅里叶变换矩阵的特征值
DFT矩阵F的特征值是e-jπ2k k=0,1,2,…，N-1，它有4个值，即1,-j,-1,j。它的特征值及重复度如表5.1所示。


表5.1DFT矩阵特征值的多样性



N1的重复度-j的重复度-1的重复度j的重复度
4mm+1mmm-1
4m+1m+1mmm
4m+2m+1mm+1m
4m+3m+1m+1m+1m

对4个特征值1,-j,-1,j来说，每个值对应的特征向量全体组成一个特征子空间，记为E0,E1,E2,E3，每个特征值的重复度决定了子空间的秩。对于N维DFT矩阵来说，重复度为M的特征值有M个独立的特征向量。关于DFT矩阵特征向量的计算有如下定理： 
定理5.2： 矩阵S可用于计算DFT矩阵F的特征向量，S的表达式为

S=210…1

12cosω1…0

012cos2ω…0



100…2cosN-1ω(5.118)

可以证明矩阵S和F满足乘法交换性，即SF=FS。因此矩阵S的特征向量也是矩阵F的特征向量，但它们对应不同的特征值。
5.2.3离散傅里叶变换矩阵的Hermite特征向量
从前面的分析可以看到，DFT矩阵F的特征向量不唯一(或特征分解不唯一)。那么，我们希望从中找到与Hermite函数形状相似的特征向量，这样的向量称为DFT矩阵的Hermite特征向量。为了求Hermite特征向量，我们需要给出一系列重要的定理及证明。
定理5.3： 对DFT矩阵的Hermite特征向量而言，它对应的连续函数的扩展方差应当为N/2πTs，Ts是信号的采样间隔。连续Hermite函数采样后得到

n(k)=12nn!N/2TshnkN/2πexp-k2π/N(5.119)

证明： 设n阶连续HermiteGaussian函数写为

Hσd,n(t)=12nn!πσdhntσde-t22σ2d(5.120)

它的傅里叶变换表达式为

Hn(f)=σd2nn!πhnfσde-f2σ2d/2(5.121)

以采样周期Ts对时间采样，式(5.120)改写为

Hσd,n(k)=12nn!πσdhnkTsσde-k2T2s2σ2d,k=1,2，…，N(5.122)

根据采样理论，DFT之后的频率分辨率为2πNTs。以k2πNTs代替式(5.121)中的f：

Hn(k)=σd2nn!πhnk2πNTsσde-k22π2σ2dN2T2s(5.123)

我们希望通过调整方差σd的大小，使得Hermite离散向量经过DFT变换后保持形状不变，即令两式的方差σd相等。式(5.122)的方差为σ2dT2s，式(5.123)的方差为N2T2s4π2σ2d，令两者相等得σd=N2πTs。代入式(5.122)得到

n(k)=12nn!N/2TshnkN/2πexp-k2π/N

上式也可看作方差为1的HermiteGaussian函数以采样间隔2πN进行采样得到的序列。
定理5.4： 若序列n(k)是由单位方差的HermiteGaussian函数以采样间隔T=2πN进行采样得到的，那么可以证明下列近似等式： 
当N为偶数时，

(-j)nn(k)≈1N∑(N/2)-1m=-(N/2)n(m)e-j2πmk/N(5.124)

当N为奇数时,

(-j)nn(k)≈1N∑(N-1)/2m=-(N-1)/2n(m)e-j2πmk/N(5.125)

证明： 这里只证明N为偶数的情况，N为奇数的情况完全类似。因为当方差σd为1时，HermiteGaussian函数n(t)是傅里叶变换的特征函数，即

F1n(t)=e-jnπ/2n(ω)=(-j)nn(ω)

或改写成

(-j)nn(ω)=12π∫+∞-∞n(t)e-jωtdt

将上式的积分区间从(-∞,+∞)截取为(-NT/2,NT/2)，可得到近似等式

(-j)nn(ω)≈12π∫NT/2-NT/2n(t)e-jωtdt(5.126)

之所以可以这样取近似，是因为当N很大时，NT=2πN也很大，并且高斯函数exp(-t2/2)的衰减很快。下面用数值积分来取代连续积分，可得

∫NT/2-NT/2n(t)e-jωtdt≈∑(N/2)-1k=-(N/2)n(kT)e-jωkTT(5.127)

因为当N很大时，T=2π/N很小，所以这是一个合理的近似。将式(5.126)和式(5.127)合并得到

(-j)nn(ω)≈12π∑(N/2)-1k=-(N/2)n(kT)e-jωkTT≈1N∑(N/2)-1k=-(N/2)n(kT)e-jωkT

将式中的ω作离散化处理，令ω=kT并代入上式得

(-j)nn(k)≈1N∑(N/2)-1m=-(N/2)n(m)e-j2πmk/N(5.128)

证毕。
从以上推导过程可以看出，式(5.124)存在两个近似误差，分别是由式(5.125)引起的截取误差和由式(5.127)引起的数值计算误差。当N趋于无穷大时，误差将趋于0。因此，N越大，式(5.126)的近似程度越好。另外，因为Hermite多项式的原因，对于n阶Hermite函数来说，其随时间的衰减率正比于tne-t2/2。因此，Hermite函数的阶次越高，其随时间的衰减也就越慢。
定理5.5： 将序列n(k)按照以下方式平移得到-n(k)。
当N为偶数时，

-n(k)=n(k),0≤k≤N2-1

n(k-N),N2≤k≤N-1(5.129)

当N为奇数时，

-n(k)=n(k)，0≤k≤N-12

n(k-N),N+12≤k≤N-1(5.130)

则-n(k)的DFT近似为(-j)n-n(k)，即当N足够大时，

(-j)n-n(m)≈1N∑N-1k=0-n(k)exp(-j2πkm/N)(5.131)

证明： 这里只证明N为偶数的情况，N为奇数的情况完全类似。

DFT-n(k)=1N∑(N/2)-1k=0n(k)e-j2πkm/N+1N∑N-1k=N/2n(k-N)e-j2πkm/N(5.132)

利用等式关系e-j2πkm/N=e-j2π(k-N)m/N，式(5.132)等号右边第二项可写为

1N∑N-1k=N/2n(k-N)e-j2πkm/N=1N∑-1l=-N/2n(l)e-j2πlm/N(5.133)

将式(5.133)代入式(5.132)并根据定理5.4，有

DFT-n(k)=1N∑(N/2)-1k=-(N/2)n(k)e-j2πkm/N≈(-j)nn(m)(5.134)

其中，m取值范围为(0,N/2-1)。利用等式

e-j2πkm/N=e-j2πk(m-N)/N(5.135)

式(5.134)可改写为

DFT-n(k)=1N∑(N/2)-1k=-(N/2)n(k)e-j2πk(m-N)/N=-jnnm-N(5.136)

其中，m取值范围为(N/2,N-1)。将式(5.135)和式(5.136)两式合并

DFT-n(k)=1N∑N-1k=0-n(k)e-j2πkm/N≈-jn-n(m)(5.137)

证毕。
从定理5.4和定理5.5可见，Hermite函数的采样序列近似为DFT矩阵的特征向量。对Hermite函数的采样序列作归一化，记为


un=-n(0),-n(1),…，-n(N-1)T‖-n(0),-n(1),…,-n(N-1)T‖(5.138)

通过S矩阵可以得到DFT矩阵F的一组实正交特征向量，因此可以将这些特征向量作为DFT特征子空间的基向量，然后计算向量un在DFT特征子空间的投影，从而得到DFT矩阵的Hermite特征向量

u～n=∑(n-k)mod4=0〈un,vk〉vk(5.139)

因为前面已经证明un近似为DFT矩阵的特征向量,因此做投影后误差不会太大,仍然保持Hermite函数的形状。式(5.139)的含义是，先利用S矩阵求一组DFT矩阵的实正交特征向量vk作为特征空间的基，然后求连续Hermite函数的离散样本在DFT特征空间的投影。但是，由式(5.139)所得的向量不是特征空间的正交基。众所周知，若要使DFRFT矩阵满足旋转相加性，离散Hermite特征向量必须是正交的，因此，对式(5.139)所得的向量必须进行正交化处理。容易证明，位于不同特征子空间的向量已经正交，因此只需在每个子空间内部做正交化即可。总之，它的计算流程如下： 
(1) 计算连续Hermite函数的取样向量un； 
(2) 计算矩阵S的特征向量vk； 
(3) 利用式(5.139)计算Hermite特征向量u^n(未正交化)； 
(4) 对u^n进行正交化，得到正交归一化的Hermite特征向量u^k。
按照正交归一化的方法不同，Pei等分别称其为GSA方法和OPA方法。GSA方法是采用GramSchmite进行正交归一化的，而OPA方法则是采用Orthogonal Procrustes方法进行正交归一化的［29］。
5.2.4离散分数阶傅里叶变换核矩阵
连续FRFT的核函数谱分解表达式如下

Kα(t,u)=∑+∞n=0e-jnαHn(t)Hn(u)(5.140)

式中，Hn(t)表示单位方差的n阶Hermite函数。仿照上面连续FRFT核函数谱分解的形式，DFRFT的核矩阵定义为

F2α/π=U^D2α/πU^T=∑N-1k=0ejkαu^ku^Tk，N为偶数

∑N-2k=0ejkαu^ku^Tk+ejkαu^Nu^TN，N为奇数(5.141)

式中的U^，当N为奇数时U^=［u^0
u^1…u^N-1］
，当N为偶数时U^=
［u^0u^1…
u^N-2u^N］
； u^k是归一化的Hermite特征向量； D2α/π是对角阵，如下方式定义

D2απ=
e-j000…0

0e-jα0…0

00e-j2α…0






e-jαN-2

000…e-jα(N-1),N为奇数(5.142)


D2απ=e-j000…0

0e-jα0…0

00e-j2α…0





e-jα(N-2)

000…e-jαN,N为偶数(5.143)

确定变换核F2απ后，信号x(t)的DFRFT可以通过如下公式计算

Xα=F2απx=UD2απUTx(5.144)

与连续FRFT相似，信号x(t)也可以通过逆DFRFT恢复。

x=F-2απXα=UD-2απUTXα(5.145)

需要注意的是，GSA与OPA方法需要计算特征向量烦琐的正交归一化运算，使得计算复杂度增加，实用性能降低。
以上我们给出了特征分解型DFRFT的基本方法。目前，还存在很多特征分解型DFRFT。基本思路总结为，特征分解型DFRFT是由分数阶特征值和相应的离散傅里叶变换矩阵F对应的特征向量所构成。离散傅里叶变换矩阵F的特征值为{1,-j,-1,j}，所对应的特征向量构成4个特征向量空间{Ek,k=1,2,3,4}，其所对应的多样性如表5.1所示。将离散傅里叶变换矩阵F对应的特征值分数化，可构造DFRFT矩阵所对应的特征值。接着，需要构造特征值所对应的特征向量。目前，存在多种构造特征向量的方法，本节我们凝练了特征分解型DFRFT的共同机理。基于此，将众多的特征分解型DFRFT分为两大类： 基于F可交换矩阵的DFRFT和基于采样HermiteGaussian函数的DFRFT。
5.2.5基于F可交换矩阵的离散分数阶傅里叶变换
首先，若存在矩阵A，满足AF=FA，则A的特征函数与F的特征函数是等价的。基于F可交换矩阵的离散分数阶傅里叶变换(FCommuting Matrixbased DFRFT，FCDFRFT)旨在直接利用F可交换矩阵的正交特征向量来构造特征分解型DFRFT的特征向量。
根据FRFT和函数的谱分解表达式，Pei首次提出了基于特征分解型DFRFT［30］。由于可交换的矩阵具有相同的特征向量，因此可通过F可交换矩阵S的特征向量来实现离散傅里叶变换矩阵F的特征向量。最后，对应的特征分解型DFRFT为

Fp=VDαV=∑Nk=0,k≠(N-1+(N)2)e-jkαvkvTk(5.146)

其中，vk是由矩阵S得到的近似连续HermiteGaussian函数的特征向量，V=［v0v1…vN-1］，(N)2=Nmod  2。需要注意的是，当N为偶数时，最后一个特征值分配上存在一个“跳跃”。这个规则和表5.1特征值的多样性是一致的。具体地，Dα是一个对角矩阵，定义如下： 

当N为奇数时，

Dα=10…00

0e-jα…0

…

0……e-j(N-2)α0

00…0e-j(N-1)α(5.147)

当N为偶数时，

Dα=10…00

0e-jα……0



0e-j(N-2)α0

00…0e-jNα(5.148)


此DFRFT为传统的基于F可交换矩阵的离散分数阶傅里叶变换(记为传统FCDFRFT)，其可以很好地近似连续FRFT，同时严格满足酉性、可加性。矩阵S的具体推导见文献［31］。具体地，连续傅里叶变换的特征函数为连续HermiteGaussian函数，同时也是下面的二阶微分方程的特征函数

S{f(t)}=λf(t)(5.149)

其中，S=D2+FD2F-1,D=d/dt 是微分算子，F为连续傅里叶变换。因此，离散的HermiteGaussian函数可以看作连续算子S的离散矩阵S的特征向量。将连续S中的F和D2分别替换为对应的离散傅里叶变换矩阵F和D～2，可以实现离散矩阵S的推导。设Δ =1/N，D～2是利用D2f(t)的如下二阶泰勒展开式求得，连续傅里叶变换的特征函数为连续的HermiteGaussian函数，其同时也是下面的二阶微分方程的特征函数

D2f(t)=d2f(t)dt2=f(t-Δ)-2f(t)+f(t+Δ)Δ2+O(Δ2)(5.150)

最后，S可以实现对连续算子S的O(Δ2)近似。换言之，从S矩阵获得的特征向量被直接看作连续HermiteGaussian函数的近似。这是因为S的特征向量是离散的Mathieu函数,而这个函数会随着N的增加逼近于Hermite函数。根据特征向量中符号的改变数目,我们可以确定所得到的特征向量所对应的连续HermiteGaussian函数,也就是确定相应的HermiteGaussian特征向量的阶数。如果特征向量具有k个符号改变,就认为该特征向量对应第k阶HermiteGaussian函数,其相应的特征值为exp(-jkα)。对于不同的N值,特征值的分配规则归纳于表5.2。


表5.2DFRFT特征值分配规则



NDFRFT的特征值
4mexp(-jkα),k=0,1,2,…,4m-2,4m
4m+1exp(-jkα),k=0,1,2,…,4m-1,4m
4m+2exp(-jkα),k=0,1,2,…,4m,4m+2
4m+3exp(-jkα),k=0,1,2,…,4m+1,4m+2

Pei等构造了F可交换矩阵S，并将S的特征向量直接作为连续HermiteGaussian特征函数的离散近似。在文献［32］ 中，Pei 等构造了新的近似三对角交换矩阵T，T的特征向量比S的特征向量更加逼近连续HermiteGaussian特征函数。此外，指出可利用S+kT的特征向量来构造DFRFT的特征向量，进一步证明利用S+15T 的特征向量构造的DFRFT效果更佳。
为了使得DFRFT更加精确，研究学者提出了多种F可交换矩阵，其目的都是进一步提高式(5.150)二阶微分算子D2的离散精度。相关的工作分为两类： 第一类是用较小的采样间隔对D2进行离散化，即Δ=1/(rN)(r＞1)［32］； 另一类是在式(5.150)中，用高阶近似对D2f(t)进行泰勒展开。
接下来具体分析通过对D2f(t)进行高阶泰勒展开构造的F可交换矩阵，从而实现特征分解型DFRFT。文献［33］中提到，任意的F可交换矩阵都具有如下表示

K=M+FMF-1+F2MF-2+F3MF-3(5.151)

通过特定的M，可以推导出具有O(Δ2k)近似度的矩阵K，其中M可由下式推导：

f″=∑km=1(-1)m-12［(m-1)!］2δ2m(2m)!fkΔ2+O(Δ2k)(5.152)

其中，δ2fk=fk-1-2fk+fk+1,表示二阶中心差分算子。F可交换矩阵S是矩阵K的特殊情况，其具有O(Δ2)近似度。虽然，矩阵K可以实现对于二阶微分算子高的近似度，但其仅限于2k+1≤N成立。
Pei等用一种新的F可交换矩阵结构克服了如上的上界限制，借助K对称型矩阵M(具体定义见文献［34］)，式(5.151)中的F可交换矩阵K可简化为

K=M+FMF-1(5.153)

设置M2=D～2(D～2可由式(5.151)求得)，可构造具有O(Δ2k)近似度的矩阵M2k

M2k=∑km=1(-1)m-12［(m-1)!］2(2m)!(M2)m(5.154)

然而，当k很大时，此过程需要很高的计算复杂度，同时对于任意阶近似都没有闭合形式解。
在文献［35］中，Serbes提出了具有O(Δ∞)近似度的矩阵，且具有闭合形式。在式(5.154)中，当k→∞，无穷阶差分矩阵的闭合形式表达式为


M∞=limk→∞M2k=-arccos2［(M2+2)/2］(5.155)

因此，最终的F可交换矩阵为

K=M∞+FM∞F-1=E+FEF-1(5.156)

其中，E为如下对角矩阵：

diag(E)=-(2πn/N)2,0＜n＜N/2

-(2π(N-n)/N)2,N/2|
＜n＜N-1(5.157)

相比于其他的基于F可交换矩阵的DFRFT，这种方法可得到极好的F可交换矩阵，并且可交换矩阵生成的特征向量最接近连续HermiteGaussian函数。因此我们称此方法为改进FCDFRFT。
上述特征分解型DFRFT具有相同的机理，它们的目的都是构造F可交换矩阵来最大限度地近似连续算子S，其关系见图5.19。接下来，在图5.20中，我们比较了矩形窗函数(当|t|≤2,x(t)=1，否则，x(t)=0)的连续FRFT与矩形窗函数相对应的离散信号的传统FCDFRFT［30］和改进FCDFRFT［35］。连续FRFT所得的结果用实线表示，DFRFT所得的实验值用圆圈表示。实验结果表明，相比较于传统FCDFRFT 的实验结果，改进FCDFRFT的变换结果更加匹配连续FRFT的结果。



图5.19基于F可交换矩阵的DFRFT




图5.20连续FRFT结果与传统FCDFRFT、改进FCDFRFT结果的比较



(a) 矩形窗函数的2απα=0.25阶FRFT与DFRFT； 
(b) 矩形窗函数的2απα=0.75阶FRFT与DFRFT



5.2.6基于采样HermiteGaussian函数的离散分数阶傅里叶变换
上述DFRFT直接利用F可交换矩阵的特征向量作为离散傅里叶变换矩阵F的特征向量。接下来分析的DFRFT将F可交换矩阵的特征向量作为一组完备的正交基，并利用它们生成离散傅里叶变换矩阵F最终的特征向量。相关的工作分为两类： 基于正交化的方法和基于优化思想的方法。
文献［29］提出了一种新的求特征向量的方法。对于连续HermiteGaussian函数进行采样，结果记为{u-m,m=1,2,…,N}。通过证明可得u-m为离散傅里叶变换矩阵F 的近似特征向量，即(-j)mu-m≈Fu-m。然后，对u-m 进行修正得到严格的特征向量

u～m=∑(m-k)mod4=0〈u-m,vk〉vk(5.158)


其中，k=mmod4，u～m满足(-j)mu～m≈Fu～m。式(5.158)中的vk是F可交换矩阵S的特征向量，其作为一组初始标准正交基，生成严格的离散傅里叶变换矩阵F的特征向量。根据表5.2，特征向量构成4个特征子空间U～k,k=1,2,3,4。不同特征子空间的向量相互正交，为了保证DFRFT的酉性与可加性，需要对特征子空间内的向量进行正交化得到正交向量u^m。这个过程需要用到两个方法： GSA 最小化HermiteGaussian函数的采样值与正交特征向量的误差来计算(从低阶到高阶)； OPA 最小化HermiteGaussian函数构成的k的采样值与正交特征向量构成的子空间U^k的误差。最后，特征分解型DFRFT表示为


F2α/π=U^D2α/πU^T(5.159)

此式由式(5.146)中将V替换为U^所得。基于GSA和OPA算法的DFRFT有很好的性能，可以很好地近似连续FRFT。
Hanna提出序列Sequential OPA (SOPA)算法用来产生矩阵F的近似HermiteGaussian特征向量［36］。根据谱展开理论，矩阵F可以表示为

F=∑4k=1λkPk(5.160)

其中，λk为DFT矩阵的特征值； Pk为在F的第k个特征空间上的正交投影矩阵。对特征空间Pk应用奇异值分解技术，有

Pk=VkVHk，k=1,2,…,4(5.161)

其中，Vk为酉矩阵，上标H表示复共轭转置。那么矩阵F的第k个特征空间的标准正交基可以从Vk的列中得到。然后，类似文献［37］中提到的一样，我们先得到DFT矩阵的HermiteGaussian特征向量u～m，然后利用奇异值分解方法对得到的特征向量应用GSA、OPA和SOPA方法得到标准正交特征向量u^k。这里，SOPA 是基于OPA提出的逐次估计的方法，同时SOPA生成的特征向量与GSA所生成的特征向量是一致的。进一步考虑，GSA、OPA 和SOPA 在不同初始特征向量情况下，生成的最后特征向量是不变的。因此，后续的关于生成不同的初始特征向量的研究对于最终的特征分解型DFRFT是没有影响的，详见文献［38］。
下面介绍基于优化方法的特征分解型DFRFT。Pei等在限制Fu^m=(-j)mu^m和u^m1u^m2=0,(m1≠m2)条件下，利用Largrange 乘数法生成特征向量u^m。在计算过程中，这两个限制条件通过实部和虚部分别执行。为了更有效地计算特征向量，文献［39］使用一种直接的分批次的技术评估方法(Direct and Batch Technique Evaluation，DBEOA)。实验证明DBEOA和OPA 的性能是等价的，其中DBEOA的计算复杂度更低。文献［40］致力于对DBEOA方法进行序列估计(Sequential Operation of DBEOA，DSEOA)。可以证明，DSEOA、GSA和SOPA的性能等价。此外，DSEOA相较于GSA和SOPA具有更强的鲁棒性。文献［41］中提出了具有闭合形式解的HermiteGaussian函数型特征向量，其中用到的正交化方法是GSA，其具有IPDFRFT和传统FCDFRFT相似的性能。
总的来说，上述基于采样HermiteGaussian函数的特征分解型DFRFT的关系如图5.21所示。对于任意的初始特征向量,由GSA、OPA、DBEOA和DSEOA生成的最后的特征向量都是相等的。和文献［29］提出的DFRFT相比较，后续DFRFT为严格正交的离散傅里叶变换矩阵F的特征向量的计算提供了不同的方法，在如何更好地使特征向量逼近连续HermiteGaussian函数方面没有明显的提升。尽管最新的方法在计算复杂度方面有一定的减少，但是相比于整个离散变换O(N2)，其效果是甚微的。因此最后的计算复杂度仍然是O(N2)。


图5.21基于采样HermiteGaussian函数的DFRFT的正交方法之间的关系


上述基于采样HermiteGaussian函数的特征分解型DFRFT旨在解决两个问题： ①通过对连续HermiteGaussian函数进行采样得到离散的向量为离散傅里叶变换矩阵F的近似特征向量； ②将①中的近似特征向量转化为严格的特征向量并正交化，保证DFRFT具有酉性和可加性等性质。
至此，我们对主要的特征分解型DFRFT进行了分析。此类DFRFT满足理想DFRFT准则的(1)~(4)。由于特征向量具有正交性，因此准则(3)和(4)显然成立。此类DFRFT满足准则(1)和(2)。但是，特征分解型DFRFT具有较高的计算复杂度O(N2)，在高效性方面有待提升，即不满足准则(5)。
特征分解型DFRFT是目前为止唯一严格意义上的DFRFT定义，这种类型的DFRFT和连续FRFT非常接近，同时还具有称为“分数阶”的独特性质——旋转相加性。当DFRFT满足旋转相加性时，p阶DFRFT的逆变换就是-p阶DFRFT。因此，正变换和逆变换具有同样的表达式，只是相差一个参数，这样通过一个统一的计算机程序就可以容易地把分数阶傅里叶域的信号处理移植到时域。在没有实时性要求的前提下，采样型DFRFT的计算也可以通过特征分解型DFRFT完成。同时，这种特征分解型DFRFT还可以用于需要旋转相加性的地方，比如图像加密。这种类型的DFRFT的缺点是，它不能写成闭合形式，同时还缺少O(NlogN)的计算方法。
5.3线性加权型离散分数阶傅里叶变换
5.3.1基于离散傅里叶变换的线性组合

这种方法的思想十分简单，是研究人员最早探索出的一种FRFT快速算法。该算法利用Taylor级数展开和CaylayHamilton定理，将任意阶次的DFRFT表示为恒等算子、DFT、时间反转算子和IDFT的线性组合，即阶次依次为0，1，2，3的DFRFT的线性组合。该算法具有计算速度快的特点，且符合可逆性和旋转相加性。但是，其计算结果与连续FRFT有较大误差，因此使用率不高。下面简单介绍其原理。
前面讨论过，分数阶傅里叶算子是一个旋转算子，可以理解为传统傅里叶算子的分数阶幂。那么离散的FRFT与传统的离散傅里叶矩阵之间是什么关系呢？该算法直接给出如下式子： 

Fp［x］=W2απ(x)(5.162)

其中，W是传统DFT矩阵，x是信号，Fp表示信号的p阶FRFT，p=2α/π。
从上式可以看出，F4=W4=I=W0=A0，F1=W1，F2=W2。这些性质与连续情况下的性质相似。这几条性质可以将Fp理解为n,k平面上的旋转算子。算子Fp的Taylor级数展开可以写为

Fp=a0(α)I+a1(α)W+a2(α)W2+a3(α)W3(5.163)

我们知道，算子W是一个酉算子，有N个正交归一化的特征向量vi，其特征值分解可以写为

W=∑i∈N1vivHi-∑i∈N2vivHi+j∑i∈N3vivHi-∑i∈N4vivHi(5.164)

其中,N1是属于λ=1的特征向量索引集合，N2、N3、N4是分别属于特征值为λ=-1、λ=j、λ=-j的特征向量索引集合。这里分数阶傅里叶算子Fp是DFT矩阵W的幂，也是酉矩阵，特征值分解可以写为


Fp=∑i∈N1ej4k1αvivHi-∑i∈N2ej(4k2+2)αvivHi+

j∑i∈N3ej(4k3+1)αvi
vHi-∑i∈N4ej(4k4-1)αvivHi

(5.165)

其中，ki∈Z。
这样得到的DFRFT算子具有以下几个性质。
(1) 酉性： FpH=Fp-1=F-p
(2) 角度可加性： FpFq=Fp+q
(3) 周期性： Fp+4k=Fp,k∈Z
那么如何计算DFRFT呢？直接计算需要对DFT矩阵W进行特征值分解来得到分数阶算子Fp，对于较大的阶次来说这是不现实的。我们可以将分数阶傅里叶算子写成如下方程： 

Fp=TΛ2απTH=∑3i=0ai(α)Wi(5.166)

其中，T是DFT矩阵W的特征向量矩阵，从而也是分数阶傅里叶算子Fp的特征向量矩阵。从上式可以得到

Λ2απ=∑3i=0ai(α)ΛWi(5.167)


其中，ΛWi=THWiT是一个对角矩阵。它是含4个未知变量ai(α),i=0,1,2,3的N个线性方程组，其中只有4个是独立的。当N≤4时，会有部分特征值漏掉。对N≥5，DFT矩阵W包含所有4个特征值，并且这些系数独立于阶次N，可以通过上面方程组算出。例如可以取ki=0,i=1,3,4,5，取特征值λ1=1，λ2=-1，λ3=j，λ4=-j，可以唯一确定ai(α),i=0,1,2,3。
当N=5时，由方程组可以得到如下相互独立的线性方程组： 

1111

1-11-1

1j-1-j

1-j-1ja0(α)

a1(α)

a2(α)

a3(α)=expj4k1α

expj2α2k3+1

expjα4k4+1

expjα4k5-1(5.168)

该方程组的解(当ki=0,i=1,3,4,5时)为

a0(α)=121+expjαcosα

a1(α)=121-expjαsinα

a2(α)=-121-expjαcosα

a3(α)=-121+expjαsinα(5.169)

这些系数满足上面讨论的几个性质，而且对于选定的ki来说，它们是唯一的。将系数代入式(5.166)可以得到Fp。众所周知，FFT的计算复杂度为ONlogN，所以这种算法的计算复杂度也为ONlogN。
我们在讨论FRFT的数学定义时指出，FRFT算子可以用传统傅里叶变换算子的特征函数的作用来定义。事实上，可以证明该算法得到的离散分数阶傅里叶算子Fp并不满足连续FRFT算子的这种定义。

Fpφn=exp-jpnπ2φn(5.170)

其中，φn是DFT矩阵的特征向量。有关证明方法这里不详细讨论。这里得到的DFRFT算子Fp只满足当α=kπ2,k∈Z时的特征方程。因此，这种数字计算方法并不具备良好的理论基础，而且得到的结果与连续FRFT还有很大的偏差，因此得不到实际应用。
5.3.2基于离散分数阶傅里叶变换的线性组合
前面的特征分解方法是基于DFT核矩阵的特征分解，并用DFT的Hermite特征向量来构建DFRFT的核矩阵。对于固定点数，DFT的Hermite特征向量和DFRFT核矩阵可以事先计算出来,但是,输入信号与DFRFT核矩阵的乘积仍然需要很大的运算量,因此,特征分解方法的计算量是O(N2)。文献［42］在特征分解方法的基础上为了减小运算量提出了一种线性组合离散算法。它将线性组合的特定阶数从0，1，2，3改为从0开始依次间隔4/N的N个阶数，N为输入的离散样本数。利用FRFT的旋转相加性，该算法可以采用串行的方式实现，即前次变换的结果作为后次变换的输入，那么所需要的特定阶数离散变换的结果数目可以减少为1个，该算法便于超大规模集成电路(VLSI)的实现，但是它需要已知至少一个特定阶数的变换结果，且这些特定阶数依赖于输入信号的样本数目N。
由于本算法是在特征分解方法的基础上得到的，因此这里为了方便重写特征分解方法的DFRFT计算公式。当计算点数N为奇数时，

Xα=∑N-1n=0e-jkανkνTkx(5.171)

当计算点数N为偶数时，

Xα=∑N-2n=0e-jkανkνTkx+e-jNανkνTkx(5.172)

该方法的基本思想是通过一组特定角度(阶次)的DFRFT的加权和来计算任意角度(阶次)的DFRFT。下面给出两个定理。
定理5.6： 设x表示一个N点离散信号(N为奇数)，则信号x的α角DFRFT可以用下式计算

Xα=∑N-1n=0Bn,αXnβ(5.173)

其中，β=2π/N。加权系数Bn,α表示为

Bn,α=IDFTe-jkαk=0,1,2，…，N-1=1N∑N-1n=0e-jkαej2πNnk(5.174)

证明： 因为Bn,α等于e-jkαk=0,1,2，…，N-1的IDFT变换，因此，Bn,α的DFT变换就等于e-jkαk=0,1,2，…，N-1，即

e-jkα=∑N-1n=0Bn,αe-jβnk=∑N-1n=0Bn,αe-j2πNnk(5.175)

因此，特征分解方法的定义式可以转换为


Xα=∑N-1k=0e-jkανkνTkx

=∑N-1k=0∑N-1n=0Bn,αe-j2πNnkνkνTkx

=∑N-1n=0Bn,αXnβ
(5.176)

Bn,α可以通过求e-jkαk=0,1,2，…，N-1的IDFT计算得到。加权系数Bn,α的闭式解为


Bn,α=1N1-e-j(N-1)(α-nβ)1-e-j(α-nβ),α≠nβ

δ(n-k),α=nβ

(5.177)

总之，对于奇数点DFRFT的计算可以通过求一组特定角度的DFRFT的加权和得到，特定角度是(2π/N)的整数倍，加权系数可以通过求IDFT得到。当信号点数为偶数时，DFRFT特征值的分布会有一个跳跃，以上的计算方法不再成立，因此，需要针对偶数情况设计新的计算方法。
定理5.7： 设x表示一个N点离散信号(N为偶数)，则信号x的α角DFRFT可以用下式计算

Xα=∑Nn=0Bn,αXnβ(5.178)

其中，β=2π/N+1，加权系数表示为

Bn,α=IDFTe-jkαk=0,1,2，…，N=1N+1∑Nn=0e-jkαej2πN+1nk(5.179)

证明： 定理5.7和定理5.6相比存在两个不同点： ①定理5.7是N+1项相加，定理5.6是N项相加； ②定理5.6的特定角度是2π/N的整数倍，定理5.7的特定角度是2π/(N+1)的整数倍。采用与定理5.6类似的方法。将序列e-jkαk=0,1,2,…,N表示为Bn,α的DFT,即

e-jkα=∑Nn=0Bn,αe-j2πN+1nk(5.180)

将其代入式(5.106)得


Xα=∑N-2k=0e-jkανkνTkx+e-jNανNνTNx

=∑N-2k=0∑Nn=0Bn,αe-j2πN+1knνkνTkx+∑Nn=0Bn,αe-j2πN+1NnνNνTNx

=∑Nn=0Bn,αXnβ
(5.181)

与定理5.6类似，也可以求得加权系数Bn,α的闭式解

Bn,α=1N+11-e-jN(α-nβ)1-e-j(α-nβ),α≠nβ

δ(n-k),α=nβ
(5.182)

从定理5.6和定理5.7可以得到结论，对于任意角度的DFRFT都可以表示为一组特定角度的DFRFT的加权和，当信号点数N为奇数时，特定角度为2π/N的整数倍； 当信号点数N为偶数时，特定角度为2π/(N+1)的整数倍。不论是奇数还是偶数，加权系数都可以通过IDFT运算获得。但是对于奇数需要计算N点IDFT，对于偶数需要计算N+1点IDFT。可以看出不论是奇数还是偶数，求加权系数都需要计算奇数点的IDFT。因此，通常计算2的幂次方点数的FFT算法不适合，只能使用针对质数长度的Winograd FFT算法。
利用DFRFT的旋转相加性，以上线性组合方法可以通过串行方法实现。以奇数点为例，式(5.167)可以改写成如下形式： 


Xα=BN-1,αF(N-1)2β/πx+BN-2,αFN-22β/πx+…+B1,αF2β/πx+B0,αx

=F2β/π…F2β/πF2β/πBN-1,αXβ+BN-2,αx+BN-3,αx+

BN-4,αx+…+B0,αx
(5.183)

上式表明，任意角度的DFRFT可以仅通过一个特定角度的DFRFT来实现。如果能够快速计算出一个特定角度的DFRFT，就可以利用串行实现方法快速计算出任意角度的DFRFT。另外，由于串行实现方法具有规则的计算结构，因此适合于VLSI实现。
5.4特殊的离散分数阶傅里叶变换
5.4.1ZoomFRFT

以上介绍的FRFT分解型算法是从连续信号f(x)的N个均匀采样出发,最后得到它的FRFTfa(x)在整个分数阶谱区间上的N点等间隔采样值，因此所计算的是在分数阶域上的具有固定分辨率的全局谱。然而在许多实际应用中,人们不仅要了解信号在FRFT全局谱上的分布情况,更对分数阶域上的某段局部谱的细节感兴趣,因此要求FRFT算法应具有对局部谱进行高分辨分析的能力。针对此问题，文献［43］提出一种FRFT高分辨计算方法,这种算法可根据需要灵活选择变换输出的局部谱区域和分辨率,当计算的分辨率很高即采样间隔很小时,可得到局部谱的精细结构。在保持计算量不变的条件下,这种算法的作用就好像摄影中的变焦镜头,能够根据需要方便地调整谱的显示范围和分辨率，因此，文献［43］将这种FRFT的高分辨计算方法称为ZoomFRFT。
下面给出这种高分辨算法的推导过程。由5.1.1节Ozaktas采样型算法可以看出，第二种分解方法包含了两个离散化步骤： 第一步利用香农内插公式对时域变量离散化； 第二步对分数阶傅里叶域变量离散化。高分辨算法在第二步采用了一种更加灵活的离散化方法。假定要计算FRFT在局部谱区间u1,u2上的M点等间隔取样值，u1、u2和M的取值任意，将分数阶傅里叶域变量离散化为u=u0+mΔI，-M/2≤m≤M/2，其中u0=u2-u1/2表示区间中点，ΔI=u2-u1/M-1表示采样间隔。然后将其代入式(5.23)，得到


Xpu0+mΔI=Aα2Δxejπγ(u0+mΔI)2·

∑Nn=-NejπβΔI2Δx(-2mn)e-j2πβu0n2Δxejπγn2Δx2xn2Δx

(5.184)

其中，γ、β、Δx定义见5.1.1节。将恒等式-2mn=m-n2-m2-n2代入上式,最后整理得到


Xpu0+mΔI=Aα2Δxejπγ(u0+mΔI)2e-jπβ12ΔIΔx(mΔI)2
∑Nn=-Nej2πβΔIΔxm-n2Δx2·

ejπ-2βxin2Δx+γ-2βΔIΔxn2Δx2xn2Δx

(5.185)

式(5.185)的计算要求直接确定u0和ΔI两个参数的取值，这在编程实现时很不方便，我们更愿意使用它们的以Δx为基准的相对取值。为此我们引入了两个相对因子，一个是“平移因子”λ=u0/Δx，-0.5≤λ≤0.5，它表示局部谱中心u0在整个谱范围-Δx/2,Δx/2中的相对位置； 另一个是“变焦因子”P=1/2ΔIΔx,它表示局部谱的分辨率ΔI相对于标准分辨率1/2Δx的放大倍数，一般取P为大于1的自然数。将关系式u0=λΔx和ΔI=1/2PΔx代入式(5.185)得到


XpλΔx+m2PΔx=Aα2ΔxejπγΔx2λ2ej2πγλΔxm2PΔxejπγ-Pβm2PΔx2·

∑Nn=-Nejπ1Pβm-n2Δx2ejπγ-1Pβn2Δx2e-j2πβλΔxn2Δxxn2Δx
(5.186)

若记x-［n］=xn2Δx,α,λ,Pm=Fα［x］λΔx+m2PΔx,则式(5.186)可简化为


α,λ,Pm=AαejπγNλ22NejπγλPmejπγ-Pβ4P2Nm2∑Nn=-Nejπβ4PNm-n2·

ejπPγ-β4PNn2e-jπ(βλ)nx-n

(5.187)

上式中的求和部分为离散卷积形式，因此可以利用FFT来实现快速的数值计算。以上是当0.5≤|p|≤1.5时的计算公式，同5.1.1节一样可得当0≤|p|≤0.5或1.5≤|p|≤2时的计算公式。
为了直观展示利用高分辨算法实现FRFT局部谱高分辨分析的效果，我们通过几个仿真实例进行分析。第一个例子我们考虑chirp信号，假定信号含有调频率相同而中心频率不同的三个chirp分量，可写成如下形式： 

x(t)=∑3i=1Aiexpj2πfit+jπμt2rectt/Δx(5.188)

其中，信号幅度分别为A1=2, A2=1, A3=1,中心频率分别为f1=2,f2=-1,f3=-0.8,调频率为μ≈0.3，归一化宽度为Δx=10，以1/Δx=0.1为间隔采样,得到N=100点信号样本，图5.22给出了信号的时域波形。由于chirp信号在与调频率相垂直的分数阶傅里叶域上具有最好的能量聚集性，因此可利用FRFT进行chirp信号的检测和参数估计，对应的FRFT阶次应为p=-0.8145。我们首先利用FRFT标准算法计算出Xp(u)在全局区间［-5,5］内以1/Δx=0.1为间隔采样的N=100点样本输出，如图5.23所示，从中可粗略地观察到信号在整个分数阶傅里叶域上的谱分布情况。

图5.23中chirp信号的能量主要集中在［-2.5,-1.5］和［0.5,1.5］两个窄带区域。由于采样间隔较大，显示的信号波形很粗糙。为了仔细观察局部谱内的波形细节,我们利用窄带高分辨FRFT算法将分辨率提高10倍来显示局部谱的波形，如图5.24和图5.25所示。图中的虚线表示图5.23中的标准分辨率谱线，实线为高分辨率谱线，它的谱分辨率比标准分辨率提高了10倍。


图5.22信号的时域波形




图5.23用FRFT标准算法得到的全局谱




上面的仿真实例充分展示了高分辨算法带来的优异的谱分析性能。通过高分辨分析可以使我们准确地观察到局部谱的每个细节，如主瓣和旁瓣的幅度、位置、宽度和过零点位置等，而标准谱线由于分辨率不够而损失掉很多细节信息。不仅如此，FRFT高分辨算法还可以为chirp信号的参数估计带来好处。在基于FRFT的chirp信号检测与参数估计中，一般通过检测谱线峰值点的位置来估计chirp信号中心频率。由图5.24、图5.25可以看出，标准谱线(虚线)的间隔较宽，很难正好采到连续谱的峰值点，这样,当以低分辨率的谱线峰值位置来估计chirp信号的中心频率参数时，会造成较大的误差； 反之,高分辨谱线(实线)很密,它的峰值点与连续谱的峰值位置误差很小，这样可大大提高chirp信号中心频率的估计精度。


图5.24区间-2.5,-1.5的高分辨局部谱




图5.25区间0.5,1.5的高分辨局部谱



5.4.2单点快速计算
ZoomFRFT算法可实现任意FRFT局部谱的高分辨分析，但在实际应用中仍然存在一些局限，主要体现在以下几方面： ①ZoomFRFT在输入/输出点数相差不大时具有较高的运算效率，当输出点数很少时运算效率低，也就是说，ZoomFRFT算法并不适合输出点数很小的场合； ②虽然ZoomFRFT可以根据需要灵活选择分辨率，但是分数阶谱的采样间隔只能是均匀的，若需要输出分数阶谱的若干非均匀采样值，则ZoomFRFT无法实现； ③在某些FRFT应用中，可能需要计算分数阶谱上的任意一个单采样点的值。为了更加有效地进行FRFT的数值计算，文献［5］又提出了一种FRFT单点快速计算方法。利用FRFT单点快速算法，我们既可以解决上述的非均匀采样点输出的问题，也可以提高少量点输出时的计算效率； 同时，在应用中还可以和ZoomFRFT配合使用，在运用ZoomFRFT计算分数阶谱之后，再运用单点的计算方法对某个感兴趣的点单独计算。因此,这种FRFT单点的计算方法具有较好的实际应用价值。
由于FRFT单点计算中利用了Horner运算以实现快速计算，先对其做一简要介绍。Horner运算是针对幂次多项式的求解而提出的一种基于循环迭代的计算方法，它实际上可以理解为一种幂次多项式从后向前逐次累加运算的过程。设一个幂次多项式表示为

B(z)=b0+b1z+b2z2+…+bMzM(5.189)

则Horner运算的递推过程可以表示为




InitializeS=0

Fori=M,M-1,…,0do

S=bi+Sz


(5.190)

具体的计算过程为
设初始值为S=0
当i=M时，S=bM
当i=M-1时,S=bM-1+bMz
当i=M-2时,S=bM-2+bM-1z+bMz2
……
当i=0时,S=b0+b1z+b2z2+…+bM-1zM-1+bMzM

当循环结束后，S值就是多项式的最终计算结果。可以看出，Horner迭代算法的运算效率较高。因为运用直接计算的方法计算多项式,需要MM-1/2次乘法和M-1次加法运算; 而运用Horner方法计算多项式，只需要M次乘法和M-1次加法运算。Horner算法已经应用于傅里叶变换的单采样点计算中［44］。
下面介绍FRFT单点快速计算的原理。其主要思想是将5.1.4节式(5.22)转化成如式(5.189)所示的幂次多项式形式，然后利用Horner迭代算法进行快速计算。
先考虑一般的任意非零单点计算，再考虑零点计算的特殊情况。若计算FRFT在任意非零点ui≠0处的值，将u=ui代入式(5.22)得到

Xpui=Aα2Δxexpjπγu2i∑Nn=-Nexp-j2πβuin2Δxexpjπγn2Δx2xn2Δx(5.191)

令


bn=expjπγn2Δx2xn2Δx(5.192)

zi=exp-jπβΔxui(5.193)


则式(5.191)的求和部分可表示为

Bzi=∑Nn=-Nbnzni=z-Nib-N+b-N+1zi+…+b0zNi+b1zN+1i+…+bNz2Ni(5.194)

利用Horner迭代运算过程式(5.184)可快速计算出Bzi。最后得到FRFT在分数阶傅里叶域ui点处的结果为

Xpui=Aα2Δxexpjπγu2iBzi=Aα2Δxexpjπγu2iBexpjπβΔxui(5.195)

但是还有一个重要问题需要考虑，即在计算多项式系数bn时需要先算出序列

gn=expjπγn2Δx2，-N≤n≤N(5.196)

若按式(5.196)直接计算序列gn，对每个点都要进行复指数运算，生成序列gn的计算量很大，并且因为gn并非固定序列，它会随着γ=cotα改变，因此无法将gn事先计算好并预存在存储器中。为减小计算量，可以采用递推的方法。因为gn是偶序列，只需计算其在0≤n≤N范围的取值即可。生成序列gn的递推公式推导如下： 若令Dn=expjπγ4Δx22n+1，有递推公式

gn+1=expjπγ4Δx2n+12=gnDn(5.197)

其中，g0=1。再令W=expjπγ2Δx2，有递推公式

Dn+1=expjπγ4Δx22n+1+1=DnW(5.198)

其中，D0=expjπγ4Δx2。这样，只要计算出D0和W的值，就可以先由递推公式(5.192)得到Dn序列，再由式(5.197)得到gn序列。
以上为任意非零单点ui≠0的计算方法，若计算FRFT 在零点的值，则式(5.191)简化为

Xp0=Aα2Δx∑Nn=-Nexpjπγn2Δx2xn2Δx=Aα2Δx∑Nn=-Ngnxn2Δx(5.199)

可见，零点计算不需要Horner迭代运算，只要利用递推方法计算出gn序列，代入式(5.199)即可。
下面通过仿真实例对上述的FRFT单点计算方法进行分析。我们将利用FRFT单点快速算法来计算信号的分数阶傅里叶谱。为了方便与Ozaktas采样型算法所计算的分数阶傅里叶谱比较，我们在仿真中通过对相应频点的重复计算，得到FRFT的N点输出。信号则选取两种典型的信号： chirp信号和矩形信号。分别计算其在不同阶次的变换结果。
例1: chirp信号


从以上仿真可以看出,FRFT的单点快速算法在输出数值较大时，它与分解型算法所得到的结果的计算误差很小，如图5.26(c)和图5.27(c)所示。当输出数值较小时，单点计算的结果与分解型算法所得到的结果有一定误差,如图5.26(b)和图5.27(b)所示。这主要是因为单点快速算法应用了Horner迭代算法，会产生一定的迭代误差积累，当输出数值较小时，积累误差比较明显。因为在实际应用中，我们运用单点快速算法主要是计算一些峰值点的数值，因此在实际应用中,因为迭代所产生的误差可以忽略。另外，利用单点算法进行FRFT的谱计算时,当同样输出N点时，它的计算量与分解型算法的计算量相当。



图5.26chirp信号的FRFT,阶次为p=-0.81



(a) chirp信号 (μ=0.3)； (b) FRFT全局谱； (c) FRFT局部谱； (d) FRFT局部谱

圆圈为单点计算结果,方形为Ozaktas采样型算法的计算结果





图5.26（续）




图5.27矩形信号的FRFT,阶次 p=0.6



(a) 矩形信号； (b) FRFT全局谱［-5,5］； (c) FRFT局部谱［-1,1］； (d) FRFT局部谱［1,5］

圆圈为单点计算结果,方形为Ozaktas采样型算法的计算结果



例2: 矩形信号





教学视频






5.5其他离散分数阶变换
由于DFT和其他酉变换有着密切的关系，比如离散余弦变换、离散Hadamard变换、离散Hartley变换等，因此可以得到这些离散变换的分数阶变换形式。
许多离散变换的“分数阶”形式都是基于前面所讲的特征分解方法，因为这种方法使得所得到的分数阶变换具有旋转相加性，并且当变换阶数等于1时退化为原来的变换。这些离散变换分数阶化的统一方法是： 首先研究离散变换的特征结构，以得到这些变换的特征值和特征向量，然后利用特征分解结构得到相应的分数阶形式。一旦得到了N×N的变换矩阵K的特征值λn和特征向量un，那么分数阶变换矩阵Kp可以表示为

Kp=UΛpUH=∑N-1n=0λpnunuHn(5.200)

其中，矩阵U和Λ由标准正交特征向量un和特征值的分数阶幂λpn组成。注意，这里我们使用共轭转置而不是转置，是因为并不是所有变换的特征向量都是实的。
Pei在文献［45］中引入了离散分数阶Hartley变换(DFRHT)，并讨论了DFRHT的特征结构： 其特征值为{1,-1}，特征向量与F的特征向量一样，都是HermiteGaussian特征向量。
Pei在文献［46］中引入了离散分数阶Hadamard变换。其特征值为{1,-1}，其特征向量通过从最初的两个特征向量利用回归算法精确地计算从阶数2一直到阶数2n的Hadamard特征向量。
Pei和Cariolaro分别从不同角度定义了离散分数阶余弦和正弦变换［4748］。Pei选择DCTⅠ和DSTⅠ来定义离散分数阶余弦Ⅰ变换(DFRCTⅠ)和离散分数阶正弦Ⅰ变换(DFRSTⅠ)［47］。DCTⅠ和DSTⅠ核矩阵的特征值为{1,-1}，并且它们的特征向量可以通过DFT的特征向量的计算获得。对于N点的DCTⅠ特征向量u，其可以写为

u=v0,2v1,…,2vN-2,vN-1T(5.201)


其中，v=v0,v1,…,vN-2,vN-1,vN-2,…,v
1T是(2N-2)点的DFT矩阵的偶特征向量。

N点的DSTⅠ特征向量u可以写为

u=2v1,v2,…,vN-1,vNT(5.202)

其中，v=0,v0,v1,…,vN,0,-vN,-vN-1,…,-v1T是2(N+1)点DFT矩阵的奇特征向量。
不同于文献［47］中定义的DFRCT，Cariolaro从DCTⅡ中得到了实值的DFRCTⅡ。DCTⅡ的特征值对于不同的N是不同的，特征值的星座图是象限对称的［48］。因为DCTⅡ的特征值是不同的，所以它的特征向量是相互正交的，也就可以直接从DCTⅡ矩阵计算其特征向量。DFRCT和DFRST可以用于计算偶信号和奇信号的DFRFT，还可以用于数字水印［47,49］。读者若想对DFRCT、DFRST和分数阶离散广义和偏移DFT，DHT，DCTⅣ和DSTⅣ有深入研究，可以参见文献［4951］。
其他的离散分数阶变换可以参见文献［5254］。在文献［30］中，基于DFRFT的特征分解定义了离散分数阶Hilbert变换。在文献［53］中定义了分数阶随机变换，这个定义同样可以写作式(5.194)的形式，只是矩阵U的特征向量是依赖于一个随机矩阵。在文献［54］中，通过对式(5.194)中的对角矩阵Λ取不同的分数阶幂，DFRFT可以推广到多参数DFRFT。文献［55］中，角度分解方法用到了其他分数阶酉变换中。
参考文献

［1］Kraniauskas P,Cariolaro G,Erseghe T.Method for defining a class of fractional operations［J］.IEEE Trans.Signal Processing,1998,46: 28042807.

［2］Erseghe T,Kraniauskas P,Cariolaro G.Unified fractional Fourier transform and sampling theorem［J］.IEEE Trans.Signal Processing,1999,47: 34193423.

［3］Ozaktas H M,Arikan O,et al.Digital computation of the fractional Fourier transform［J］.IEEE Trans.Signal Processing,1996,44(9): 21412150.

［4］Pei S C,Ding J J.Closedform discrete fractional and affine Fourier transform［J］.IEEE Trans.Signal Processing,2000,48(5): 13381353.

［5］赵兴浩,邓兵,陶然.分数阶Fourier变换数字计算中的量纲归一化［J］.北京理工大学学报,2005,25(4): 360364.

［6］Crochiere R E,Rabiner L R.Interpolation and decimation of digital signalA tutorial review［J］.Proc.IEEE,1981,69(3): 300331.

［7］Bidet E,Castelain D,Joanblanq C.A fast singlechip implementation of 8192 complex point FFT［J］.IEEE J.SolidState Circuits,1995,30(3): 300305.

［8］Pekurovsky D.P3DFFT: A framework for parallel computations of Fourier transforms in three dimensions［J］.SIAM J.Sci.Comput.,2012,34(4): C192C209.

［9］Hoyer E A,Stork R F.The zoom FFT using complex modulation［C］//IEEE Int.Conf. Acoust.Speech Signal Process.(ICASSP),1977: 7881.

［10］Markel J.FFT pruning［J］.IEEE Trans.Audio Electroacoust.,1971,19(4): 305311.

［11］Pang C,Liu S,Han Y.Highspeed target detection algorithm based on sparse Fourier transform［J］.IEEE Access,2018: 11.

［12］Lin J,Feng Y,Liu S.Fast ISAR imaging based on sparse Fourier transform algorithm［C］//2017 Int.Conf.Advanced Infocomm Technology (ICAIT),2017,334339.

［13］Fan T,Shan T,Liu S,et al.A fast pulse compression algorithm based on sparse inverse Fourier transform［C］//2016 CIE Int.Conf.Radar,2016,13621366.

［14］Liu S,Zhang Y D,Shan T.Detection of weak astronomical signals with frequencyhopping interference suppression［J］.Digital Signal Processing,2018,72: 18.

［15］Liu S,Zeng Z,Zhang Y D,et al.Automatic human fall detection in fractional Fourier domain for assisted living［C］//41st IEEE Int.Conf.Acoust.,Speech,Signal Process.(ICASSP),2016,799803.

［16］Liu S,Shan T,Zhang Y D,et al.A fast algorithm for multicomponent LFM signal analysis exploiting segmented DPT and SDFRFT［C］//IEEE International Radar Conference,2015,11391143.

［17］Liu S,Shan T,Tao R,et al.Sparse discrete fractional Fourier transform and its applications［J］.IEEE Trans.Signal Process,2014,62(24): 65826595.

［18］Hassanieh H,Indyk P,Katabi D.Simple and practical algorithm for sparse Fourier transform［C］//23rd Annu.ACMSIAM Symp.Discrete Algorithms.,2012,11831194.

［19］Zhang H,Shan T,Liu S,et al.Parameter optimization of sparse Fourier transform for radar target detection［C］//2020 IEEE Radar Conference,2020,343347.

［20］Xia X G.Discrete ChirpFourier transform and its application to chirp rate estimation［J］.IEEE Trans.Signal Process.,2000,48(11): 31223133.

［21］Peleg S,Friedlander B.Multicomponent signal analysis using the polynomialphase transform［J］.IEEE Trans.Aerospace Electron.Syst.,1996,32(1): 378387.

［22］Zhang H,Shan T,Liu S,et al.Optimized sparse fractional Fourier transform: Principle and performance analysis［J］.Signal Processing,2020,174: 107646.

［23］Liu S,Zhang H,Shan T,et al.Efficient radar detection of weak maneuvering targets using a coarsetofine strategy［J］.IET Radar Sonar Navig.,2021,15(2): 181193.

［24］Almeida L B.The fractional Fourier transform and timefrequency representations［J］.IEEE Trans.Signal Process.,1994,42(11): 30843091.

［25］Nelson J.Sketching and streaming algorithms for processing massive data［J］.XRDS: Crossroads,2012,19(1): 1419.

［26］Zhang H,Shan T,Liu S,et al.Performance evaluation and parameter optimization of sparse Fourier transform［J］.Signal Processing,2021,179: 107823.

［27］Hassanieh H,Indyk P,Katabi D,et al.Nearly optimal sparse Fourier transform［C］//44th Symp.Theory Comput.,2012: 563578.

［28］Simon M K,Alouini M S.Some new results for integrals involving the generalized Marcum 
Qfunction and their application to performance evaluation over fading channels［J］.IEEE Trans.Wireless Commun.,2003,2(4): 611615.

［29］Pei S C,Yeh M H,Tseng C C.Discrete fractional Fourier transform based on orthogonal projections［J］.IEEE Trans.Signal Processing,1999,47(5): 13351347.

［30］Pei S C,Yeh M H.Improved discrete fractional Fourier transform［J］.Optics Letters,1997,22(14): 10471049.

［31］Dickinson B,Steiglitz K.Eigenvectors and functions of the discrete Fourier transform［J］.IEEE Transactions on Acoustics Speech and Signal Processing,1982,30(1): 2531.

［32］Pei S C,Hsue W L,Ding J J.Discrete fractional Fourier transform based on new nearly tridiagonal commuting matrices［J］.IEEE Trans.Signal Processing,2006,54: 38153828.

［33］Candan C·.On higher order approximations for HermiteGaussian functions and discrete fractional Fourier transforms［J］.IEEE Signal Processing Letter,2007,14(10): 699702.

［34］Pei S C,Hsue W L,Ding J J.DFTcommuting matrix with arbitrary or infinite order second derivative approximation［J］.IEEE Transactions on Signal Processing,2009,57(1): 390394.

［35］Serbes A,Durak Ata L.Efficient computation of DFT commuting matrices by a closedform infinite order approximation to the second differentiation matrix［J］.Signal Processing,2011,91(3): 582589.

［36］Hanna M T,Seif N P A,Ahmed W A E M.HermiteGaussianlike eigenvectors of the discrete Fourier transform matrix based on the singular value decomposition of its orthogonal projection matrices［J］.IEEE Trans.Circuits Syst.I,2004,51: 22452254.

［37］Candan C,Kutay M A,Ozaktas H M.The discrete fractional Fourier transform［J］.IEEE Trans.Signal Processing,2000,48: 13291337.

［38］Hanna M T,Seif N P A,Ahmed W A E M.HermiteGaussianlike eigenvectors of the discrete Fourier transform matrix based on the direct utilization of the orthogonal projection matrices on its eigenspaces［J］.IEEE Transactions on Signal Processing,2006,54(7): 28152819.

［39］Hanna M T.Direct batch evaluation of optimal orthonormal eigenvectors of the DFT matrix［J］.IEEE Transactions on Signal Processing,2008,56(5): 21382143.

［40］Hanna M T.Direct sequential evaluation of optimal orthonormal eigenvectors of the discrete Fourier transform matrix by constrained optimization［J］.Digital Signal Processing,2012,22(4): 681689.

［41］Neto J R D O,Lima J B.Discrete fractional Fourier transforms based on closedform HermiteGaussianlike DFT eigenvectors［J］.IEEE Transactions on Signal Processing,2017,65(23): 61716184.

［42］Yeh M H,Pei S C.A method for the discrete fractional Fourier transform computation［J］.IEEE Trans.Signal Processing,2003,51(3): 889891.

［43］赵兴浩,陶然,邓兵,等.分数阶傅里叶变换的快速计算新方法［J］.电子学报.2007,35(6): 10891093.

［44］Orfanidis S J.信号处理导论［M］.清华大学出版社，Prentice Hall,1998.

［45］Pei S C,Tseng C C,et al.Discrete fractional Hartley and Fourier transform［J］.IEEE Trans.Circuits Syst.II,1998,45: 665675.

［46］Pei S C,Yeh M H.Discrete fractional Hadamard transform［C］.Proc.IEEE Int.Symp.,Circuits and Systems,1999,179182.

［47］Pei S C,Yeh M H.The discrete fractional cosine and sine transforms［J］.IEEE Trans.Signal Processing,2001,49: 11981207.

［48］Cariolaro G,Erseghe T,Kraniauskas P.The fractional discrete cosine transform［J］.IEEE Trans.Signal Processing,2002,50: 902911.

［49］Tseng C C.Eigenvalues and eigenvectors of generalized DFT,generalized DHT,DCTIV and DSTIV matrices［J］.IEEE Trans.Signal Processing,2002,50: 866877.

［50］Pei S C,Ding J J.Generalized eigenvectors and fractionalization of offset DFTs and DCTs［J］.IEEE Trans.Signal Processing,2004,52: 20322046.

［51］VargasRubio J G,Santhanam B.On the multiangle centered discrete fractional Fourier transform［J］.IEEE Signal Processing Letters,2005,12: 273276.

［52］Pei S C,Yeh M H.Discrete fractional Hilbert transform［C］//Proc.IEEE Int.Symp.Signal Processing,1998,506509.

［53］Liu Z J,Zhao H F,Liu S T.A discrete fractional random transform［J］.Optics Communications,2005,255: 357365.

［54］Pei S C,Hsue W L.The multipleparameter discrete fractional Fourier transform［J］.IEEE Signal Processing Letters,2006,13: 329332.

［55］Yeh M H.Angular decompositions for the discrete fractional signal transforms［J］.Signal Processing,2005,85: 537547.