第3章分数阶算子及分数阶变换
作为傅里叶变换的广义形式，分数阶傅里叶变换是一种统一的时频分析方法，且可以理解为角度α的时频平面旋转，因此依据分数阶傅里叶变换可以定义一些有用的分数阶算子和分数阶变换。
3.1分数阶算子
在讨论分数阶变换之前，我们先来解释一下分数阶算子的概念。假设有算子O(·)，

O(g(x))=G(y)(3.1)

那么，它的分数阶算子(记作Oa(·)，其中a是实数)满足下列性质： 
(1) 边界性

O0(g(x))=g(x)，O1(g(x))=G(y)(3.2)

(2) 可加性

Ob(Oa(g(x)))=Oa(Ob(g(x)))=Oa+b(g(x))(3.3)

由于参数a可任意选取，一些由经典算子不能解决的问题，可用更为灵活的分数阶算子解决。
容易得出一些运算的分数阶形式，例如，可以将乘法运算O(g(x))=g(x)f(x)的分数阶算子定义为

Oa(g(x))=g(x)fa(x)(3.4)

不过，对大多数运算来说，它们的分数阶算子并不显而易见，我们通过如下方法来求得。
假设算子O可以被分解为O-11、O2和O1，即

G(y)=O(g(x))=O-11O2O1(g(x))(3.5)

且O2的分数阶算子已知，则由此可得O的分数阶算子为

Ga(y′)=Oa(g(x))=O-11Oa2O1(g(x))(3.6)

如果O2是乘法运算，那么它的分数阶算子的导出可进一步简化。由第2章内容，我们知道对于傅里叶变换

G(ω)=F(g(t))=12π∫+∞-∞e-jωtg(t)dt(3.7)

它的分数阶算子，也就是分数阶傅里叶变换(FRFT)定义为

Gp(u)=OpF(g(t))=1-jcotα2π
ej2u2cotα
∫+∞-∞e-jutcscα
ej2t2cotαg(t)dt(3.8)

其中，α=pπ/2。由上式可以看出分数阶傅里叶变换能够分解为如下三步（图3.1）［1］： 


图3.1分数阶傅里叶变换的分解步骤示意图


(1) 乘以chirp信号

y(t)=(1-jcotα)ejt22cotαg(t)




(2) 傅里叶变换(自变量存在尺度转变)，Gp(u)=Ycscαu，其中，

Y(ω)=F1［y］(ω)=12π∫+∞-∞y(t)e-jωtdt

(3) 乘以chirp信号

Gp(u)=eju22cotαGp(u)

因此，我们可以将传统基于傅里叶变换的对偶算子(O、O～)推广到基于分数阶傅里叶变换的分数阶对偶算子(Op、O～p)，其示意图如图3.2所示，图中O表示时域算子，O～表示O的频域对偶算子，Op表示O的p阶分数阶算子，O～p表示Op的p阶分数阶傅里叶域对偶算子，也就是O～的p阶分数阶算子。依据图3.2，我们可以得到分数阶卷积的对偶形式。本章3.3节还将依据对偶转换的概念得到另外一些分数阶对偶算子。



图3.2分数阶对偶算子构造示意图





教学视频








教学视频




3.1.1分数阶卷积
在传统的傅里叶变换理论中，两个函数f(t)与g(t)的卷积定义为

h(t)=(fg)(t)=∫+∞-∞f(t)g(t-τ)dt(3.9)

或简记为h=fg。其运算特性是： 一个域的卷积对应于另一个域的乘积，即


f(t)*g(t)=h(t)FTH(ω)=2πF(ω)×G(ω)(3.10)

2πf(t)×g(t)=h(t)FTH(ω)=F(ω)G(ω)(3.11)


接下来将这种传统的卷积理论直接推广到分数阶傅里叶域中［2］。设s(t)=f(t)×g(t)，则


Sp(u)=Aαej(u2/2)cotα∫+∞-∞f(t)g(t)ej(t2/2)cotα-jtucscαdt

=|cscα|2πej(u2/2)cotα∫+∞-∞Fp(v)Gu-vcscαe-j(v2/2)cotαdv
(3.12)

其中，Aα=1-jcotα，A-α=1+jcotα，α=pπ/2。这里所用的傅里叶变换为如下对称形式： 

F(ω)=12π∫+∞-∞f(t)e-jωtdt(3.13)

同样，令h(t)=f(t)g(t)，则H(ω)=2πF(ω)G(ω)。我们知道Hp(u)可以看作对H(ω)作旋转角度为α-π/2(即p-1阶)的分数阶傅里叶变换。这样，根据式(3.12)可以得到


Hp(u)=Fp-12πF(ω)G(ω)

=secαe-j(u2/2)tanα∫+∞-∞Fp(v)gu-vsecαej(u2/2)tanαdv
(3.14)

显然式(3.12)、式(3.14)这样的表示并不令人满意，因为它们并没有像经典傅里叶卷积定理那样明确的对偶表示形式。因此，根据图3.2，我们可以得到另一种表达形式。
对于任意函数f(t)，定义函数f′(t)=f(t)ejcotα2t2，f″(t)=f(t)e-jcotα2t2，则两个函数f(t)与g(t)的分数阶卷积可以定义为

h(t)=(fpg)(t)=Aα2πe-jcotα2t2f(t)ejcotα2t2*g(t)ejcotα2t2(3.15)

同样，它们的分数阶乘积定义为


s(t)=f×pg(t)=ejcotα2t22πf″g″

=ejcotα2t22πf(t)e-jcotα2t2g(t)e-jcotα2t2
(3.16)

这样，分数阶卷积的实现就可分解成如图3.3所示的结构［3］。


图3.3分数阶卷积分解结构(Cα=cotα/2)


通过推导，可以得到如下的分数阶卷积理论表示［2］。
若h(t)=(fpg)(t)，则


Hp(u)=e-jCαu2Fp(u)×Gp(u)(3.17)

Fpf(t)g(t)ejCαu2=A-αFp×pGp(u)(3.18)


此时，分数阶傅里叶域的卷积与乘积理论的表达形式就比较明确了： 两个函数在参照域的分数阶卷积对应于它们的分数阶傅里叶变换乘积再乘以一个线性调频信号； 两个函数在参照域的乘积再乘以一个线性调频信号对应于它们的分数阶傅里叶变换的分数阶乘积再做幅值调制。
特别地，当α=0(即p=0)时，两个函数的分数阶卷积和分数阶乘积就退化成两个时间函数的经典卷积和乘积； 当α=π/2(即p=1)时，两个函数的分数阶卷积和分数阶乘积就退化成两个频域函数的经典卷积和乘积。
3.1.2分数阶相关
相关(Correlation)是信号处理特别是信号检测理论中一个十分重要的概念。相关运算主要用来对两个信号进行比较，或者从一个强信号中发现某种微弱信号。经典意义下的相关运算只限于在时域或频域进行，随着分数阶傅里叶变换的出现，有必要将相关运算推广到任意分数阶傅里叶域，由此产生了分数阶相关(Fractional Correlation)的概念。分数阶相关是经典相关的推广，下面将由相关的概念出发，导出分数阶相关的几种定义，分析它的时变性质，并介绍分数阶相关在信号检测中的应用。
3.1.2.1分数阶相关的多种定义
若时域表示的输入信号和参考信号为s(t)和h(t)，则两信号的时域相关和频域相关分别定义为


(sh)(τ)=∫s(t)h*(t-τ)dt(3.19)

(sπ/2h)(ν)=∫s(t)h*(t)e-j2πνtdt(3.20)


其中，表示相关符号。在信号处理领域，为了便于数学上表示各种变换关系，会定义一些算子，如傅里叶变换用算子F表示，分数阶傅里叶变换用算子Fp或Fα(α=pπ/2)表示。为了描述相关运算，文献［4］引入了移位算子的概念。时移算子Tτ和频移算子Oν分别定义为

(Tτs)(t)=s(t-τ)(3.21)

和

(Oνs)(t)=ej2πνts(t)(3.22)

将时域相关用时移算子表示可写成

(s0h)(τ)=∫s(β)h*(β-τ)dβ=〈s,Tτh〉(3.23)

同样地，频域相关也可以用频移算子来表示，即

(sπ/2h)(ν)=∫s(t)h*(t)e-j2πνtdt=〈s,Oνh〉(3.24)

从式(3.23)和式(3.24)可以发现，所谓两个信号的时域相关或频域相关，是指将一个信号做时移或频移再与另一个信号做内积的输出结果。时移算子和频移算子虽然不同，但是根据傅里叶变换的性质可以建立两者之间的等效关系。这种关系可以写成

Oν=F-(π/2)TνF(π/2)(3.25)

因为分数阶傅里叶变换可看作傅里叶变换的一种广义形式，仿照上式中的关系，并且将式中的傅里叶变换算子Fπ/2用分数阶傅里叶变换算子Fα代替，可定义一种新的分数阶移位算子Rαρ，它和时移算子之间的关系可写成［4］

Rαρ=F-αTρFα(3.26)

它可以理解为信号在分数阶傅里叶域产生某种位移。可以证明，分数阶移位算子具备如下两个特性


Rαρ1Rαρ2=Rαρ1+ρ2(3.27) 

(Rαρ)-1=Rα-ρ(3.28)


用分数阶移位算子Rαρ作用一个时域信号s(t)将会使信号产生什么样的变化呢？利用分数阶傅里叶变换的性质不难推导出

(Rαρs)(t)=s(t-ρcosα)ej2π(ρ2/2)cosαsinα+j2πtρsinα(3.29)

可以看出，算子Rαρ对信号的作用效果既有时移分量ρcosα，又有频移分量ρsinα。由信号的时频表示出发，我们也可以从信号时频分布的变化来理解算子Rαρ的作用。如图3.4所示，Rαρ的作用是将信号s(t)的Wigner分布在时频平面上沿α角方向移动距离ρ。


图3.4分数阶移位算子Rαρ对信号的时频分布的作用效果


可以看出，式(3.29)定义的分数阶平移算子Rαρ表示时频平面上沿角度α方向(逆时针测得)的径向平移，这样(Rαρs)(t)就是信号s(t)经分数阶傅里叶域平移后的时域表示。可以看出，当α=0和α=π/2时，分数阶平移算子就分别退化为时间平移算子和频率平移算子，即(R0ρs)(t)=s(t-ρ)和Rπ/2ρs(t)=s(t)ej2πtρ。
用分数阶平移算子定义分数阶傅里叶域的相关概念需要把线性时不变系统(Linear TimeInvariant,LTI)推广为线性分数阶移不变系统，参见图3.5和图3.6。输入信号与系统响应函数的分数阶相关可以解释为输入信号与系统脉冲响应h(t)分数阶平移ρ后的内积


CORRps(t),h(t)(ρ)=〈s,Rαρh〉

=e-j2πρ22cosαsinα∫s(t)h*(t-ρcosα)e-j2πtρsinαdt
(3.30)

其中，α=pπ/2，p为分数阶傅里叶变换的阶数。


图3.5一个简单的LTI系统




图3.6对应的LFRT系统



实际上，与时移或频移算子相比，分数阶移位算子是一种更为一般化的移位算子，时移或频移算子是它的两个特例。当α=0时，分数阶移位算子退化为时移算子； 当α=π/2时，分数阶移位算子退化为频移算子。
此外，我们还可导出另外两个等价的分数阶相关定义式，这两种定义可用于分数阶相关的离散计算。
根据分数阶傅里叶变换的酉特性［5］，有〈s,h〉=〈Fαs,Fαh〉，由此，可以导出分数阶相关的第二种定义式，即


sαh(ρ)=〈s,Rραh〉=〈s,F-αTρFαh〉=〈Fαs,TρFαh〉

=∫SαβHαβ-ρ*dβ=Sα0Hα(ρ)
(3.31)

这表明分数阶相关可以通过输入信号与脉冲响应函数分别进行分数阶傅里叶变换后再进行经典相关计算得到。按照此定义，计算分数阶相关的方法是先对信号求旋转角为α的分数阶傅里叶变换，再在与α相应的分数阶傅里叶域内作经典相关运算。其中相关运算比较耗时，可以利用傅里叶变换的卷积性质将相关运算转换为乘积运算。于是，可导出分数阶相关的第三种定义式，即

sαh(ρ)=Sα0Hα(ρ)=F-π/2Sπ/2+α(u)Hπ/2+α(u)*(ρ)(3.32)

这里利用了p+1阶分数阶傅里叶变换等价于对信号的p阶分数阶傅里叶变换作傅里叶变换。这种计算方法首先对信号求转角为α+π/2的分数阶傅里叶变换，然后求它们的乘积，再求逆傅里叶变换。分数阶傅里叶变换和逆傅里叶变换都有快速算法，所以第三种定义是一种计算分数阶相关的有效方法。以上定义均为针对确定信号的分数阶相关定义，第8章将进一步分析针对随机信号的分数阶相关函数。
3.1.2.2分数阶相关的性质
前面介绍了分数阶相关的几种定义及计算方法，下面将介绍其特性及应用。众所周知，相关运算是移不变的，也就是说，如果输入信号有一个延迟，相关输出也有同样一个延迟，但波形不变。根据相关的移不变特性，在一些应用中，含有目标回波信号的输入信号与参考信号作相关运算(匹配滤波)，不论目标信号在什么位置，只要与参考信号相匹配就能被检测到，而且相关峰值的位置表明了被检测目标的位置。分数阶相关一般不具备上述特性，它的最大特点在于它的移变属性。移变特性是分数阶相关最为显著的特性，也是正确地将分数阶相关运用于信号检测的基础。下面将从分数阶相关第一种定义出发，重点分析其移变特性。
首先对分数阶相关定义式(3.30)的两端取模，可以得到

g(ρ)=sαh(ρ)=∫s(t)h*t-ρcosαe-j2πtρsinαdt(3.33)

假定h0(t)为要检测的目标信号，输入信号s(t)=h0t-tinp，参考信号h(t)=
h0t-tref，对函数g(ρ)求导并令其等于0，可得到分数阶相关峰值的位置为

ρpeak=(tinp-tref)secα(3.34)

在ρ=ρpeak处，相关峰值为

g(ρpeak)==∫|h0(t)|2e-j2πttinp-treftanαdt(3.35)

可以看出，分数阶相关的峰值位置和高度与角度α以及目标相对位置均呈现一定关系，我们将针对各种情况来进行分析。当α=0(经典相关)时，相关峰值位置ρpeak=tinp-tref，相关峰值g(ρpeak)=∫|h0(t)|2dt，这表明经典相关峰值高度与目标信号的位置无关，不管目标在何处，相关峰值都保持相同的高度； 对于分数阶相关，当tinp=tref时，峰值位置在零点，且相关峰值达到最大高度g(ρpeak)=∫|h0(t)|2dt； 当tinp≠tref时，随着tinp-tref的增大，g(ρpeak)中复指数项的振荡越来越快，峰值高度不断减小。复指数项振荡频率取决于两个因素，即tinp-tref和tan，这两项的增大会导致指数项的振荡加快，从而使峰值高度降低。峰值高度随着输入信号和参考信号之间的距离增加而降低，而降低的快慢由分数阶的阶数决定。因此，通过调整分数阶角度α,可控制相关峰值高度对输入信号和参考信号位置差异的敏感度。为了量化这种关系，设Δt表示|h0(t)|2的支撑区间，为避免振荡和退化，指数项必须小于π，于是可以得到如下条件

tinp-trefΔt≤cotα(3.36)

这个公式定义了在相关峰值允许的退化情况下输入信号和参考信号的最大位置差异。由式(3.36)我们发现，当α=0时，公式右端为无穷大，无论输入信号和参考信号位置差异多大，相关都有最大峰值； 若0＜α＜π/2，相关峰值会随着位置差异的增加而减小，而α越大，减小的速率越快。
根据分数阶相关的上述移变性质，我们可将其用于某种特定的目标检测系统，该系统只检测参考信号一定范围内的目标，而不检测此范围以外的目标，范围可由分数阶次来控制。文献［5］已经将其成功用于机器人的视觉识别系统中。由于一个机器人在行走时需要发现并躲避前方的障碍物，于是要求它的视觉系统能发现一定距离内的障碍物，在这种视觉识别系统中，利用分数阶相关可以只检测某个距离以内，而忽略范围以外的目标。
在分析了分数阶相关的移变性质的基础上，可以进一步研究分数阶自相关与模糊函数的关系问题。根据第2章的分析，分数阶傅里叶变换与Wigner分布之间有一个重要的关系，即一个信号的分数阶傅里叶变换的模平方等于它的Wigner分布在分数阶域上的投影。那么分数阶自相关与模糊函数究竟是什么关系呢？首先，我们给出经典相关和模糊函数的一组关系。若用AFs(τ,ν)表示信号x(t)的模糊函数，则有


s0s(τ)=AFs(τ,0)(3.37)

sπ/2s(τ)=AFs(0,ν)(3.38)


由于经典相关是分数阶相关的特例，可以将这种关系推广到分数阶相关，下面的定理描述了分数阶相关和模糊函数的关系。
定理3.1： 一个信号的分数阶自相关等于它的模糊函数沿相应α角的径向切面。用公式表示为

sαs(ρ)=AFsρcosα,ρsinα(3.39)

证明： 模糊函数可以看作信号关于时延τ和频移ν的二维联合相关函数。它的定义式为

AFs(τ,ν)=∫st+τ2s*t-τ2e-j2πνtdt(3.40)

信号s(t)的分数阶自相关可表示为

sαs(ρ)=〈s,Rαρs〉=ej2π(ρ2/2)cosαsinα∫s(t)s*t-ρcosαe-j2πtρsinαdt(3.41)

为了变为对称形式，令t=t′+ρ2cosα，做变量替换可得

sαs(ρ)=∫st′+ρ2cosαs*t′-ρ2cosαe-j2πt′ρsinαdt′(3.42)



图3.7分数阶相关与
模糊函数的关系图



将上式与模糊函数的定义式作比较，可建立分数阶自相关和模糊函数之间的关系

sαs(ρ)=AFsρcosα,ρsinα(3.43)

式(3.43)表明，信号在旋转角为α的分数阶傅里叶域上的自相关等于其模糊函数在α角射线方向的切片，如图3.7所示。

3.1.2.3分数阶相关在信号检测中的应用
根据定理3.1，文献［4］提出了一种基于分数阶相关的线性调频(Linear Frequency Modulation,LFM)信号检测和参数估计的有效算法。这种算法主要是利用了分数阶相关和模糊函数的关系，在只需要估计信号的调频率参数的情况下，应用这种算法比较方便。类似的方法还有文献［6］提出的基于模糊函数和Radon变换的算法。由于LFM信号的模糊函数是在模糊平面上的一条过原点的斜线，而斜率为信号的调频率，为了检测参量未知的LFM信号，通过对模糊函数所有过原点的直线求Radon变换，可构造一个检测统计量

D(m)=∫AFs(τ,mτ)dτ(3.44)

若D(m)在某一点存在一个峰值且超过了预先确定的阈值，则认为检测到LFM信号，并且可由峰值点的位置估计出信号的调频率。
式(3.44)中右端的被积函数是过原点、斜率为m的模糊函数的切片。它可以由角度为α=arctan(m)的分数阶自相关函数代替，于是可导出一个等价的检测统计量

L(α)=∫sαs(ρ)dρ(3.45)

若以调频率为参量，检验统计量为

L(m)=∫sarctan(m)s(ρ)dρ(3.46)

虽然由模糊函数导出的检测统计量和由分数阶自相关导出的检测统计量是等价的，但两者的计算量不同。对于由分数阶自相关导出的检测统计量，具体计算包含以下几个步骤： 首先对信号进行采样得到s(n)=s(nT),n=1,2,…,N，采样周期T必须满足采样定理； 然后根据实际应用确定一组角度αk,k=1,2,…,M，计算在每个角度αk上信号的分数阶自相关并求其积分值。若在某个角度出现峰值，则可检测到信号并估计出调频率参数。这里根据第三种定义来计算信号的分数阶自相关，即

sαs(ρ)=F-π/2Sπ/2+α(u)2(ρ)(3.47)

计算过程包括： ①利用分数阶傅里叶变换的快速算法求角度α+π/2的分数阶傅里叶变换； ②求变换后的模平方； ③求逆傅里叶变换。因为分数阶傅里叶变换快速算法的计算量与快速傅里叶变换(Fast Fourier Transform,FFT)相当，为ONlog2N，所以基于分数阶自相关的LFM信号检测算法总的计算量约为OM2Nlog2N+N，基于模糊函数的算法总的计算量为ON2log2N+MN。两者比较，当MN时，用分数阶自相关作信号检测需要相对较小的计算量。
为了说明这种方法的有效性，下面给出一个仿真实例。假设有如下含离散多分量LFM信号的观测信号

s(k)=∑3i=0ej2πfi+mi(π/4096)kk+n(k),k=1,2,…,2048(3.48)

其中，参数为f0=1/1024，f1=5/1024，f2=3/1024，f3=10/1024，m0=0.124，m1=0.136，m2=0.3，m3=0.5，n(k)为加性高斯白噪声。取调频率样本值ml=0.6/200l,l=0,1,2,…,199，求每个调频率ml的检验统计量值L(ml)，图3.8为不同信噪比情况下的仿真输出结果。


图3.8基于分数阶自相关的LFM信号检测在不同信噪比下的输出结果

虚线： 无噪声，点画线： SNR=-6dB,实线： SNR＝-12dB



3.1.3分数阶酉算子和Hermite算子
在时频分析理论中，酉算子和Hermite算子是比较重要的两类算子，酉性是设计变换算子时经常需要考虑的要素之一，而不同的变换域表示往往能够通过某种Hermite算子联系起来，因此，推导分数阶酉算子和Hermite算子也得到了研究人员的广泛关注。基于时移算子和频移算子的概念(这是两种基本的酉算子)，Akay定义了分数阶位移算子T,τ，即分数阶酉算子［7］，如式(3.49)所示。利用Stones 定理可以得到分数阶Hermite算子，如式(3.50)所示。


T,τ［x］(t)=xt-τcose-jπτ2cossin+j2πtτsin(3.49)

Z［x］(t)=costx(t)+sin-j2π·ddtx(t)(3.50)


分数阶酉算子T,τ对信号的作用既有时移分量τcos，又有频移分量τsin，它是将信号在时频平面上沿着角度为的轴移动径向距离τ，所以称为分数阶位移算子。它与分数阶傅里叶变换的关系如式(3.51)所示，可以看出信号范数保持不变。对比式(3.30)和式(3.49)不难发现，分数阶相关与算子T,τ的关系如式(3.52)所示。




Fp+1T,τ［x］(u)=e-j2πτuFp+1［x］(u)

FpT,τ［x］(u)=Fp［x］(u-τ),=pπ/2(3.51)


CORRp［x,y］(η)=〈x,T,τ［y］(u)〉,=pπ/2(3.52)


式中，符号〈·,·〉表示内积。
3.2分数阶变换
由分数阶算子的定义出发，可引申出分数阶变换的概念。利用分数阶傅里叶变换来构造新的分数阶变换主要有如下两种方式，一种是基于分数阶傅里叶变换是傅里叶变换的广义形式，将以传统傅里叶变换为基础的变换推广为以分数阶傅里叶变换为基础的形式，这样得到的分数阶变换主要有分数阶Hilbert变换、分数阶余弦(Cosine)变换、分数阶正弦(Sine)变换、分数阶Hartley变换等； 另一种是利用分数阶傅里叶变换是一种统一的时频分析方法、是对时频面的旋转，基于该性质可以得到一系列新的基于分数阶傅里叶变换的时频分析工具，如分数阶Gabor展开、分数阶小波包变换、分数阶模糊函数和分数阶Wigner分布等。



教学视频



3.2.1分数阶Hilbert变换
Hilbert变换(Hilbert Transform,HT)是一种重要的信号处理工具，已经在通信调制、图像边缘检测等领域得到了广泛应用。通过将对负谱的抑制由频谱扩展为分数阶傅里叶谱，可以得到分数阶Hilbert变换。在此基础上，SooChang Pei利用分数阶傅里叶变换的特征分解型离散算法给出了分数阶Hilbert变换的一种离散表达形式［8］，并将其应用于数字图像边缘检测。文献［9］进一步探讨了分数阶Hilbert变换器的设计和应用问题，提出了多种FIR、IIR分数阶Hilbert变换器的设计方法，并基于分数阶Hilbert变换对信号分数阶傅里叶变换负谱分量抑制提出了一种单边带(SSB)通信系统，利用分数阶Hilbert变换的变换阶数作为解调密钥来实现安全通信。
实信号f(t)的Hilbert变换定义为［11］

x～(t)=x(t)h(t)=1π∫+∞-∞x(τ)t-τdτ(3.53)

其中，h(t)是Hilbert变换的变换核。x～(t)的傅里叶变换为

X～(ω)=Fπ/2x～(ω)=-jsgn(ω)X(ω)(3.54)

其中，sgn(ω)为符号函数。信号x(t)通过Hilbert变换和它的复信号联系起来，可以构成解析信号［11］

x-(t)=x(t)+jx～(t)(3.54)

解析信号的一个重要性质就是保留了f(t)的正频率部分，剔除f(t)的负频率部分。与解析信号类似，反解析信号保留了f(t)的负频率部分，剔除f(t)的正频率部分，反解析信号定义为

x^(t)=x(t)-jx～(t)(3.55)

3.2.1.1分数阶Hilbert变换定义和仿真
1. 定义
文献［12］给出了两种分数阶Hilbert变换的定义，一种是参数为φ=v·π/2，v∈R的传统Hilbert变换相移版本(定义3.1)； 另一种是用分数阶傅里叶变换替换傅里叶变换在传统Hilbert变换中的角色(定义3.2)。这两种定义方式不是等价的，但是都满足在p=1或v=1情况下退化为经典Hilbert变换。利用文献［12］给出的第一种分数阶Hilbert变换定义，Tseng和Pei提出了信号的解析表示和一种以分数阶相位参数为密钥的单边带调制方法［9］。但由于实信号的分数阶傅里叶变换不满足共轭对称性，因此上述两种定义均无法直接利用分数阶傅里叶正谱来重构原实信号。
定义3.1： 信号x(t)的分数阶Hilbert变换为

x～v(t)=cosφx(t)+sinφX1(t)(3.56)

其中，X1(t)表示x(t)的1阶分数阶傅里叶变换。
定义3.2： 信号x(t)的分数阶Hilbert变换为

x～p(t)=F-pH1Fp［x(t)］(3.57)

其中，H1(ω)=expjπ2s(ω)+exp-jπ2s-ω，s(ω)表示阶跃函数，Fp表示p阶分数阶傅里叶变换算子。
此外，Zayed 还给出了另一种不同的分数阶Hilbert变换定义以及相应的解析信号表达式［13］。
定义3.3： 信号x(t)的分数阶Hilbert变换为

x～p(t)=F-1F1x(t)ejcotα2t2H1e-jcotα2t2,α=
π2p(3.58)

解析信号为

x-p(t)=x(t)+jx～p(t)

Zayed证明了由信号的分数阶傅里叶正谱能够恢复解析信号，但该方法不能由信号的分数阶傅里叶正谱恢复实信号，因此也不适用于处理实信号［13］。SooChang Pei等讨论了因果信号和实信号的分数阶Hilbert变换［14］，通过保留分数阶傅里叶谱的奇偶对称分量，达到节省传输带宽的目的，但这一方法计算过程非常复杂，缺乏实用性。
在实际应用中，信号往往是实信号。由分数阶傅里叶变换的性质可知，这类信号的分数阶傅里叶变换不满足共轭对称性，仅从分数阶傅里叶正谱(或负谱)无法恢复原信号，所以基于分数阶傅里叶变换的Hilbert变换及解析信号用于单边带通信受到一定制约。接下来，我们将基于分数阶傅里叶变换的Hilbert变换理论，从分数阶卷积理论出发，给出一种分数阶Hilbert变换新定义(定义3.4)。与已有定义不同的是，依据定义3.4得到的解析信号的分数阶傅里叶正谱能够被用来恢复实信号。因此，可以用分数阶傅里叶变换的旋转角度作为密钥，将定义3.4用于保密单边带通信系统，提高系统的安全性； 也可以用于图像加密，增强图像的保密效果。
根据文献［15］，任意信号f(t)的分数阶傅里叶变换可以看作信号g(t)=cejcotα2t2f(t)(c=1-jcotα)的傅里叶变换。那么，如果y(t)为实信号，调制后可得到f(t)=y(t)ejmt2，当m=-cotα/2时有


Fαf(u)=Aα∫+∞-∞f(t)ej(u2+t2)2cotαe-jutcscαdt

=Aαeju22cotαFπ/2［y］ucscα，α≠nπ
(3.59)

其中，Fπ/2［y］ucscα为y(t)的傅里叶变换。从式(3.59)可以看出，信号f(t)的分数阶傅里叶谱表现为y(t)的频谱乘以线性调频信号eju22cotα。因为Aα与u无关，eju22cotα关于u为偶对称，而且信号y(t)的傅里叶变换满足共轭对称性，因此保留f(t)的分数阶傅里叶正谱Fα［f］(t)u≥0，就能够得到y(t)的正频谱Fπ/2［y］ucscαu≥0，从而恢复出实信号y(t)。
定义3.4： 根据3.1.1节的分数阶卷积理论，定义信号x(t)的分数阶Hilbert变换为

x～α(t)=Aαe-jcotα2t2x′(t)h′(t)(3.60)

式中，x′(t)=x(t)ejcotα2t2； h′(t)=h(t)ejcotα2t2； α为分数阶傅里叶变换的旋转角度。则相应的解析信号为

x-α(t)=x(t)+jx～α(t)(3.61)

该解析信号x-α(t)只包含x(t)的分数阶傅里叶正谱分量，去掉了负谱分量。下面给出几个重要结论。
定理3.2： 信号x(t)的分数阶Hilbert变换的分数阶傅里叶谱为

Fαx～α(u)=-jsgn(u)Fα［x］(u)(3.62)

证明： 既然h(t)的分数阶傅里叶变换为

Fαh(u)=-jsgn(u)ejcotα2u2(3.63)

根据分数阶卷积理论，有

Fαx～α(u)=e-jcotα2u2Fαh(u)Fα［x］(u)(3.64)

将式(3.63)代入式(3.64)，定理得证。
依据定理3.2，显然可以得到

Fαx-α(u)=2Fα［x］(u)，u＞0

Fα［x］(u)，u=0

0，u＜0(3.65)

式(3.65)表明，解析信号x-α(t)只包含x(t)的正谱成分，去掉了负谱成分。
下面证明能够利用分数阶傅里叶正谱得到解析信号


x-α(t)=2∫+∞-∞x(τ)dτ∫+∞0Kαu,τK-α(u,t)du

=Aαe-jcotα2t2∫+∞-∞x(τ)ejcotα2τ2dτA-α∫+∞-∞1+sgn(u)ejcotα2u2·

e-jcotα2u2ejcscαu(t-τ)du
(3.66)

根据式(3.63)，可以得到


x-α(t)=x(t)+jAαe-jcotα2t2∫+∞-∞x(τ)ejcotα2τ2h(t-τ)ejcotα2(t-τ)2dτ

=x(t)+jx～α(t)
(3.67)

式(3.65)、式(3.67)说明，解析信号x-α(t)去除了信号x(t)的负谱分量，并且能够从正谱分量恢复解析信号x-α(t)。下面的定理给出从x-α(t)得到实信号xreal(t)，以及不同角度α时x-α(t)与解析信号x-(t)和反解析信号x^(t)的关系。
定理3.3： 解析信号x-α(t)与实信号xreal(t)及解析信号x-(t)、反解析信号x^(t)的关系如下


x-α(t)=e-jcotα2t2x-(t),0＜α＜π

e-jcotα2t2x^(t),π＜α＜2π(3.68a)

xreal(t)=Reejt22cotαx-α(t)(3.68b)


证明： 根据式(3.65)，有


x-α(t)=2∫+∞0Fα［x］(u)K-α(u,t)du

=2∫+∞0Aαeju22cotαFπ/2［y］ucscαK-α(u,t)du

=2AαA-αe-jt22cotα∫+∞0Fπ/2［y］ucscαejutcscαdu
(3.69)

其中,AαA-α=|cscα|/2π。应用∫+∞0ejcscαutdu=12∫+∞-∞1+sgn(u)ejcscαutdu，于是有

x-α(t)=e-jcotα2t2x-(t),0＜α＜π

e-jcotα2t2x^(t),π＜α＜2π
(3.70)

根据式(3.70)，可得式(3.68b)。
式(3.60)、式(3.68)的结论很容易推广到实值图像信号处理。由于实值图像的二维FFT不具有共轭对称性，应用二维Hilbert变换不能恢复原实值图像。但是应用半平面Hilbert变换理论［16］可以恢复实值图像，因此，将其推广可以得到基于分数阶傅里叶变换的半平面Hilbert变换理论。下面给出基于分数阶傅里叶变换的半平面Hilbert变换的解析信号及实值图像恢复的结果。
若y(t1,t2)为实值图像，则关于t1的调制图像f(t1,t2)=y(t1,t2)e-jcotα2t21，与f(t1,t2)对应的分数阶傅里叶变换解析信号为

f～t1α(t1,t2)=f(t1,t2)+jAαe-jcotα2t21f′(t1,t2)h′(t1,t2)(3.71)

式(3.71)的二维分数阶傅里叶变换为

Fαf～t1α(u1,u2)=2Fαf(u1,u2),u1＞0

Fαf(u1,u2),u1=0

0,u1＜0(3.72)

实值信号y(t1,t2)可以根据下面的公式得到

y(t1,t2)=Reejcotα2t21f～t1α(t1,t2)(3.73)

实值图像关于t2的基于分数阶傅里叶变换解析信号及图像恢复情况与上述处理过程类似，这里不再赘述。接下来给出建立在定义3.4相关理论上的两个仿真结果。
2. 仿真
1) 单边带通信系统
单边带通信系统结构如图3.9所示。在图3.9(a)中，实信号x(t)首先进行载波cosωct调制，得到y(t)，然后对y(t)应用分数阶Hilbert变换得到解析信号y-α(t)。由于载波信号cosωct是实信号，所以y(t)也是实信号，因此这里主要讨论对y(t)的处理。在图3.9(b)所示解调过程中，首先对信号y-α(t)进行旋转角度α的chirp解调，然后对解调后的信号取实部得到实信号y(t)，再进行载波cosωct解调，最后得到信号x(t)。


图3.9单边带通信系统框图



(a) 调制过程； (b) 解调过程


在解调系统中需要事先知道旋转角度α，如果使用的旋转角度不是α而是β，且β≠α，那么不可能从接收信号中恢复信号y(t)。为了考查旋转角度为β时恢复y(t)的情况，在式(3.68b)中代入β并整理，得到的结果为

z(t)=cossinα-β2sinαsinβt2y(t)-sinsinα-β2sinαsinβt2y-(t)(3.74)

可以看出，只有在旋转角度满足关系β=α时，才有z(t)=y(t)，否则输出的结果中存在误差，从而得不到y(t)。以单位幅值矩形实信号y(t)为例，取旋转角度α=0.3π。在系统的解调过程中，令旋转角度β在0~π变化。图3.10给出了从解析信号y-α(t)恢复实信号y(t)的情况。图3.10(a)是在β=α=0.3π时得到的与原信号相同的恢复信号z(t)=y(t)； 图3.10(b)表示在β=0.35π时得到的信号z(t)≠y(t)。从图3.10可以看出，当旋转角度β≠α时，波形存在很大误差。图3.11给出了恢复信号与原信号之间的均方误差结果。图3.11表明，只有在β=0.3π时误差最小。因此，基于定义3.4相关理论的保密单边带通信系统将旋转角度α作为密钥可以提高系统的安全性。


图3.10不同旋转角度时恢复的波形


(a) β=α=0.3π； (b) β=0.35π




图3.11不同旋转角度下恢复信号的误差




2) 图像加密
从式(3.72)、式(3.73)可以看出，通过图像分数阶傅里叶变换的半平面数据能恢复原实值图像，因此将图像分数阶傅里叶变换数据的一半作为密图，将旋转角度作为密钥，可以实现对实值图像的加密。同时，由于分数阶傅里叶变换的结果为复数，将密图数据的虚部提取出来，与实部重新组成一幅实值图像作为密图，就可以得到实值加密图像，并且对密图还可以进行重复加密，提高图像的加密效果。图3.12(a)为原实值图像，取旋转角度α=0.3π，实值加密图像如图3.12(b)所示。

对图像解密时，需要事先知道旋转角度。如果旋转的角度为β≠α，则得不到正确的结果。图3.13(a)为选择旋转角度β=0.38π的解密图像，图3.13(b)为正确的解密图像。从图3.13(a)可以看出，在角度不匹配的情况下，已经无法得到原图像的信息。恢复图像与原图像的均方误差如图3.14所示，表明只有在角度β=α=0.4π时，误差最小，说明该方法具有较好的保密效果。因为实值加密图像为实值图像，因此可以重复加密，进一步提高图像的保密性。



图3.12原图像和加密图像 


(a) 原实值图像； (b) 实值加密图像




图3.13不同旋转角度时恢复的图像


(a) β=0.38π； (b) β=α=0.4π





图3.14不同旋转角度下恢复图像的均方误差



3.2.1.2广义分数阶Hilbert变换
1. 定义
定义3.5： 将定义3.1和定义3.2结合起来就能够得到一种两参数的广义分数阶Hilbert变换［12］

x～p,v(t)=F-pHvFp［x(t)］(3.75)

其中

Hv(ω)=expjφs(ω)+exp-jφs-ω(3.76)

φ=
π2v。式(3.76)也可以写成

Hv(ω)=cosφH0(ω)+sinφH1(ω)(3.77)

该两参数广义分数阶Hilbert变换能够通过改变分数阶傅里叶变换阶次p和分数阶Hilbert变换相位φ来实现图像压缩的边缘增强。以占空比是40/256的矩形函数为例，图3.15~图3.17给出了该矩形函数在固定阶次p下不同相位φ的广义分数阶Hilbert变换结果。可以看出： ①当0＜v＜1时，输入信号的负谱成分得到加强； 而在1＜v＜2时，输入信号的正谱成分得到强化。②当p从1降到0.8时，正谱成分增大得多些； 而降到0.5时，则没有明显的区别了。


定义3.6： 文献［17］提出了一种新的广义分数阶Hilbert变换，即

x～p,v(t)=cosφe-jcotα2x(t)+sinφe-jcotα2
x～(t),
α=π2p,φ=π2v(3.78)

式中，x～(t)为x(t)的传统Hilbert变换。相应的解析信号为

x-p,v(t)=e-jt22cotαx(t)-e-jφx～p,v(t)(3.79)

对式(3.79)两边都作p阶分数阶傅里叶变换，得到

Fpx-p,v(u)=Aαeju22cotαF1［x］ucscα-e-jφFpx～p,v(u)(3.80)

由于x～p,v(t)的分数阶傅里叶变换是

Fpx～p,v(u)=Aαeju22cotαF1［x］ucscαcosφ-jsinφsgn(u)(3.81)

将式(3.81)代入式(3.80)，得到

Fpx-p,v(u)=2Aαjsinφe-jφeju22cotαF1［x］ucscα,u≥0

0,u＜0
(3.82)



图3.15矩形函数的广义分数阶Hilbert变换结果，p=1






图3.16矩形函数的广义分数阶Hilbert变换结果，p=0.8





图3.17矩形函数的广义分数阶Hilbert变换结果，p=0.5



从上式可以看出，解析信号的p阶分数阶傅里叶变换只包含正谱成分。当φ≠0,π且α≠0,π时，解析信号含有必需的半个分数阶傅里叶正谱，只是多了一个复因子2Aαjsinφe-jφeju22cotα。既然chirp信号eju22cotα是对称的，且F1［x］ucscα也是共轭对称的，那么实信号x(t)能够通过x-p,v(t)重构，这样我们就能够抑制实信号的一半分数阶傅里叶谱宽。与定义3.5的广义分数阶傅里叶变换定义一样，定义3.6也同样拥有两个参数，但是它们也存在一些不同： ①前者在变换域定义，后者是在时域定义； ②后者的定义用于单边带通信，以两个参数作为双重密钥，而前者则用在图像处理上。
设g(t)=-jejt22cotα+jφcscφx-p,v(t)，则


g(t)=-jejφcscφx(t)-e-jφcosφx(t)+sinφx～1(t)

=x(t)-jcotφx(t)+jcscφcosφx(t)+sinφx～1(t)(3.83)


可以看出g(t)的实部是x(t)，这样我们便能通过下式来恢复x(t)

x(t)=Re-jejt22cotα+jφcscφx-p,v(t)(3.84)

基于式(3.79)和式(3.84)，一个推广的单边带通信系统如图3.18所示。图中，uc是分数阶傅里叶域的载波频率，α是分数阶傅里叶变换角度，φ是分数阶Hilbert变换相位。


图3.18推广的单边带通信系统



(a) 系统框图； (b) 调制框图； (c) 解调框图


输出信号y(t)如下

y(t)=x-p,v(t)ejuctcscα=x-p,v(t)k(t)ejt22cotα(3.85)

其中,k(t)=e-jt22cotαejuctcscα。根据分数阶卷积定理［2］，有

Fpx-p,v(t)k(t)ejt22cotα(u)=A-αeju22cotαF～px-p,v(u)F～pk(u)(3.86)

其中，F～p=e-ju22cotαFp，“”表示卷积，k(t)的分数阶傅里叶变换为

Fpk(u)=2πAα|cscα|eju22cotαδu-uc(3.87)

将式(3.87) 代入式(3.86)，得到

Fp［y］(u)=e-ju2c2cotα+jucucotαFpx-p,vu-uc(3.88)

考虑噪声时，接收信号为y(t)+n(t)，解调输出信号为x(t)+v(t)，其中

v(t)=cscφRe-je-juctejt22cotα+jφn(t)(3.89)

从上式可以看出，噪声增益主要取决于因子cscφ，这与文献［9］的结果类似。通过在区间π/6,π/2对φ取值，可以将cscφ限制在［1,2］范围内，这样能降低噪声的影响。
图3.18所示的单边带通信系统是对传统单边带通信系统的推广。对于解调来说，参数α和φ必须是已知的，否则很难从接收信号中恢复出x(t)。以β,γ替换式(3.84)中的α,φ，并结合式(3.78)和式(3.79)，则解调出的实信号可表示为


f(t)=Re-jcscγejt22cotβ+jγx-α,φ(t)

=Re-jcscγejt22(cotβ-cotα)+j(γ-φ)ejφx(t)-cosφx(t)-sinφx～(t)
(3.90)

令ξ=t22cotβ-cotα+γ-φ，则


f(t)=Re-jcscγejξejφx(t)-cosφx(t)-sinφx～(t)

=sinφsinγRecosξx(t)-sinξx～(t)+jsinξx(t)+jcosξx～(t)

=sinφsinγcosξx(t)-sinξx～(t)
(3.91)

由于cotβ-cotα=sinα-β/sinαsinβ，用sinα-β2sinαsinβt2+γ-φ替换ξ，可以得到


x(t)=sinφsinγcossinα-β2sinαsinβt2+γ-φx(t)-

sinsinα-β2sinαsinβt2+γ-φx～(t)

(3.92)

可以发现，只有当β=α且γ=φ时才有f(t)=x(t)，否则很难从接收信号中得到x(t)。这样，便能以α,φ为双重密钥来增强单边带通信系统的安全性。
2. 仿真
设待发送实信号为矩形脉冲，α=0.3π，φ=0.4π。我们通过两种不同的方式来进行单边带通信。一种方式是直接利用实信号及其去掉chirp信号e-j0.363t2的分数阶Hilbert变换来构成解析信号； 另一种方式是完整利用式(3.79)定义的解析信号形式。图3.19给出了忽略噪声影响后的解调器输出。从图3.19(a)可以看出，由于发送的解析信号不包含x(t)完整信息而导致解调出的信号严重失真。而图3.19(b)的结果正相反，基于完整的定义3.6发送的信号被相当准确地恢复出来。


图3.19解调器输出信号



(a)  没有chirp信号e-j0.363t2； (b) 具有chirp信号e-j0.363t2


接下来，将利用上述的矩形脉冲进行在不同参数下解调的效果比较。式(3.92)指出，只有在β=α且γ=φ时才能实现原信号的不失真解调。那么当β≠α且γ≠φ时，又会有什么样的结果呢？我们将利用解调信号和原信号幅度的均方根误差来进行分析，计算公式如下


E=∑N-1n=0f～(n)-f(n)2

式中，f～(n)和f(n)分别是解调信号和原信号的采样序列，N是采样点数。β取值区间是0.1π,0.9π，γ取值区间是0.1π,π。图3.20给出了相应仿真结果。
从图3.20可以看出，β和γ都能影响解调效果： ①当γ==0.4π时，误差只取决于参数β，且当β∈0.25π,0.35π时，误差迅速降低为0； ②当β=α=0.3π时，误差只取决于参数γ； 同时，曲线两边的误差值都比较大，且缓慢地过渡到谷底； ③当β≠α且γ≠φ时，误差都比较大； 只有当γ=φ=0.4π和β=α=0.3π时，误差才等于0。这样，参数α,φ就能够用作单边带保密通信的密钥。


图3.20解调的幅度误差均方根



3.2.2分数阶余弦变换、分数阶正弦变换
正弦变换、余弦变换和Hartley变换都属于酉变换，已经在图像压缩和自适应滤波方面得到了广泛应用，利用它们与傅里叶变换的关系，我们可以得到分数阶正弦变换、余弦变换和Hartley变换［10］。需要注意的是： ①与分数阶傅里叶变换周期为4不同，分数阶正弦变换、余弦变换和Hartley变换的周期都是2； ②分数阶正弦变换没有偶特征函数，而分数阶余弦变换(Fractional Cosine Transform,FRCT)没有奇特征函数，因此，最好用分数阶正弦变换(Fractional Sine Transform,FRST)来处理奇函数，而用分数阶余弦变换来处理偶函数。
3.2.2.1定义
首先回顾第2章利用傅里叶变换的特征函数分解来导出分数阶傅里叶变换多样性的有关内容。我们知道傅里叶变换G(ω)=F1g(t)可以分解为
(1) g(t)→{am|m=0,1,2,3,…}(3.93)

其中，am=C-1m∫+∞-∞g(t)Hm(t)e-t2/2dt，

Cm=∫+∞-∞H2m(t)e-t2dt(3.94)

Hm(t)是m阶Hermite多项式。


(2) bm=(-j)mam(3.95)

(3) G(ω)=∑+∞m=0bme-ω2/2Hm(ω)(3.96)


若改变第(2)步，即令

bm=e-jm·α·π/2am(3.97)

并且保留第(1)、(3)步不变，即得到相应的分数阶傅里叶变换。余弦变换、正弦变换与傅里叶变换有如下关系


Ocos(g(t))=G(ω)+G(-ω)2

Osin(g(t))=G(ω)-G(-ω)2
(3.98)

其中，G(ω)=F1(g(t))。当m为偶数时，Hm(t)为偶函数； 当m为奇数时，Hm(t)为奇函数。由此可得如下结论： 
(1) 当m为偶数时，exp(-t2/2)Hm(t)是余弦变换的特征函数； 
(2) 当m为奇数时，exp(-t2/2)Hm(t)是正弦变换的特征函数。
它们相应的特征值记作λC(m)和λS(m)。当m为偶数时，λC(m)=(-j)m，λS(m)=0； 当m为奇数时，λC(m)=0，λS(m)=(-j)m-1。类似于分数阶傅里叶变换，可由如下过程导出分数阶余弦变换和分数阶正弦变换。
(1) g(t)→{am|m=0,1,2,3,…}(3.99)

其中，am=C-1m∫+∞-∞g(t)Hm(t)e-t2/2dt，

Cm=∫+∞-∞H2m(t)e-t2dt(3.100)

(2) dm=λαC,S(m)am(3.101)
(3) GαC,S(s)=∑+∞m=0dme-s2/2Hm(s)(3.102)

其中，λαC,S(m)分别是余弦变换和正弦变换分数幂的特征值。当m为偶数时，λαC(m)=
exp-jmαπ/2，λαS(m)=0； 当m为奇数时，λαC(m)=0，λαS(m)=exp-j(m-1)απ/2。
由式(3.101)和式(3.102)可得分数阶余弦变换

GαC(s)=∑+∞m=0e-jmαπa2me-s2/2H2m(s)(3.103)

以及分数阶正弦变换

GαS(s)=∑+∞m=0e-jmαπa2m+1e-s2/2H2m+1(s)(3.104)

由于


Gα(s)+Gα(-s)2=∑+∞m=0e-jmαπ/2ame-s2/2Hm(s)+Hm(-s)2(3.105)

Gα(s)-Gα(-s)2=∑+∞m=0e-jmαπ/2ame-s2/2Hm(s)-Hm(-s)2(3.106)


而且当m是偶数时，Hm(t)是偶函数； 当m是奇数时，Hm(t)是奇函数。因此，分数阶余弦变换、分数阶正弦变换与分数阶傅里叶变换有如下关系


GαC(s)=Gα(s)+Gα(-s)2

GαS(s)=ejαπ/2Gα(s)-Gα(-s)2
(3.107)

因此，我们得到分数阶余弦变换的定义式为


GαC(s)=OαC(g(t))=1-jcot2πej(s2/2)cot·

∫+∞-∞cos(stcsc)ej(t2/2)cotg(t)dt

(3.108)

分数阶正弦变换的定义式为


GαS(s)=OαS(g(t))=1-jcot2πej(-(π/2))ej(s2/2)cot·

∫+∞-∞sin(stcsc)ej(t2/2)cotg(t)dt

(3.109)

值得注意的是，分数阶余弦变换和分数阶正弦变换的周期为2，不同于分数阶傅里叶变换以4为周期。


OαC(g(t))=Oα+2C(g(t))
OαS(g(t))=Oα+2S(g(t))
(3.110)

由于这些分数阶变换具有可加性，因而都是可逆的，它们的逆变换为


O-αC(GαC(s))=g(t)，g(t)为偶函数

O-αS(GαS(s))=g(t)，g(t)为奇函数


(3.111)

需要注意的是，分数阶余弦变换没有奇特征函数，分数阶正弦变换没有偶特征函数。对于原始的余弦变换和正弦变换，经过分数阶余弦变换，输入函数的奇部将丢失； 经过分数阶正弦变换，输入函数的偶部将丢失。因此，最好应用分数阶余弦变换处理偶函数，应用分数阶正弦变换处理奇函数。由此，我们限制分数阶余弦变换的输入函数为偶函数并定义单边分数阶余弦变换为


GαC(s)=OαC(g(t))=2-j2cotπej(s2/2)cot

∫+∞0ej(t2/2)cotcos(stcsc)g(t)dt

(3.112)

限制分数阶正弦变换的输入函数为奇函数并定义单边分数阶正弦变换为


GαS(s)=OαS(g(t))=2-j2cotπej(-(π/2))ej(s2/2)cot

∫+∞0ej(t2/2)cotsin(stcsc)g(t)dt

(3.113)

我们注意到单边分数阶余弦变换和单边分数阶正弦变换对偶函数和奇函数的变换结果与分数阶傅里叶变换的相同，即



GαC(s)=GαF(s)，g(t)为偶函数

GαS(s)=ejGαF(s)，g(t)为奇函数

(3.114)

当应用单边分数阶余弦变换和单边分数阶正弦变换处理偶函数和奇函数时，运算复杂度是应用分数阶傅里叶变换时的一半(因为单边分数阶余弦变换和单边分数阶正弦变换的积分限是［0,+∞))。因此，应用单边分数阶余弦变换和单边分数阶正弦变换处理偶函数和奇函数，比应用分数阶傅里叶变换更有效。
3.2.2.2性质
这里我们将引出正则余弦变换和正则正弦变换，它们是分数阶余弦变换和分数阶正弦变换的进一步推广。为了引出正则余弦变换和正则正弦变换，首先给出线性正则变换(Linear Canonical Transform,LCT)的基本定义(有关线性正则变换的详细内容请参阅第13章)。类似地，线性正则变换又是分数阶傅里叶变换的广义形式，它的定义为
当b≠0时，


G(a,b,c,d)F(s)=O(a,b,c,d)F(g(t))

=1j2πbe(j/2)(d/b)s2∫+∞-∞e-j(s/b)te(j/2)(a/b)t2g(t)dt
(3.115)

当b=0时，

G(a,b,c,d)F(s)=O(a,b,c,d)F(g(t))=de(j/2)cds2g(ds)(3.116)

分数阶傅里叶变换是线性正则变换当{a,b,c,d}={cos,sin,-sin,cos}时的特例： 

OαF(g(t))=ejO(cos,sin,-sin,cos)F(g(t))，=απ/2(3.117)

那么，正则余弦变换(Canonical Cosine Transform,CCT)定义为

当b≠0时，


G(a,b,c,d)C(s)=O(a,b,c,d)C(g(t))

=1j2πbe(j/2)(d/b)s2∫+∞-∞coss/be(j/2)(a/b)t2g(t)dt

(3.118)

当b=0时，

G(a,b,c,d)C(s)=O(a,b,c,d)C(g(t))=de(j/2)cds2g(ds)(3.119)

正则正弦变换(Canonical Sine Transform,CST)定义为

当b≠0时，


G(a,b,c,d)S(s)=O(a,b,c,d)S(g(t))=1j2πbe(j/2)(d/b)s2·

∫+∞-∞-jsinst/be(j/2)·(a/b)t2g(t)dt
(3.120)

当b=0时，

G(a,b,c,d)S(s)=O(a,b,c,d)S(g(t))=de(j/2)cd s2g(ds)(3.121)

分数阶余弦变换和分数阶正弦变换是正则余弦变换和正则正弦变换当{a,b,c,d}={cos,sin,-sin,cos}时的特例，并且仅相差一个常数。


GαC(s)=ejG(cos,sin,-sin,cos)C(s),=απ/2(3.122)

GαS(s)=ejejG(cos,sin,-sin,cos)S(s),=απ/2(3.123)


下面给出正则余弦变换的性质，分数阶余弦变换的性质可由正则余弦变换的性质取{a,b,c,d}={cot,1,-1,0}得到，其中=απ/2。分数阶正弦变换和正则正弦变换的性质与分数阶余弦变换和正则余弦变换类似。
(1) 共轭： 

O(a,b,c,d)C(g(t))=O(a,-b,-c,d)Cg(t)

(2) 时移g(t-η)+g(t+η)：


O(a,b,c,d)Cg(t-τ)+g(t+τ)

=e-jacτ2/2G(a,b,c,d)C(s-aτ)e-jcτs+G(a,b,c,d)C(s+aτ)e-jcτs


(3) cos(ηt)调制：


O(a,b,c,d)Ccos(ηt)g(t)

=e-jbdη2/2［G(a,b,c,d)C(s-bη)ejdηs+G(a,b,c,d)C(s+bη)e-jdηs］/2


(4) sin(ηt)调制：


O(a,b,c,d)Csin(ηt)g(t)

=-je-jbdη2/2［G(a,b,c,d)S(s-bη)ejdηs+G(a,b,c,d)S(s+bη)e-jdηs］/2


(5) 导数：

O(a,b,c,d)Cg′(t)=aG′(a,b,c,d)S(s)-cjsG(a,b,c,d)S(s)

(5a) n阶导数(n为奇数)：


O(a,b,c,d)Cg(n)(t)=a2d2ds2-2acjsdds-ac-c2s2n-12·

aG′(a,b,c,d)S(s)-cjsG(a,b,c,d)S(s)


(5b) n阶导数(n为偶数)：

O(a,b,c,d)Cg(n)(t)=a2d2ds2-2acjsdds-jac-c2s2n2 G(a,b,c,d)C(s)

(6) 乘t：

O(a,b,c,d)Ctg(t)=dsG(a,b,c,d)S(s)+jbG′(a,b,c,d)S(s)

(6a) 乘tn(n为奇数)：


O(a,b,c,d)Ctng(t)=-b2d2ds2+2bdjsdds+jbd+d2s2n-12·

jbG′(a,b,c,d)S(s)+dsG(a,b,c,d)S(s)


(6b) 乘tn(n为偶数)：

O(a,b,c,d)Ctng(t)=-b2d2ds2+2bdjsdds+jbd+d2s2n2G(a,b,c,d)C(s)

(7) 时域翻转：

O(a,b,c,d)Cg(-t)=O(a,b,c,d)C(g(t))=G(a,b,c,d)C(s)

(8) Parseval定理：

∫+∞-∞G(a,b,c,d)C(s)2ds=∫+∞-∞Even(g(t))2dt

(9) 广义Parseval定理： 

∫+∞-∞G(a,b,c,d)C(s)H(a,b,c,d)C(s)ds=∫+∞-∞Even(g(t))Evenh(t)dt

(10) g(t)=exp-j(pt2+qt)时的正则余弦变换结果：

G(a,b,c,d)C(s)=1a-2pbej2·c-2pda-2pb·s2e-jq2b2(a-2pb)cosqsab-2pb2

事实上，g(t)=exp(-jpt2)cos(qt)和g(t)=2exp(-jpt2)cos(qt)u(t)时的变换结果同样是上式。

(11) g(t)=1时的正则余弦变换结果：

G(a,b,c,d)C(s)=1aej2·ca·s2

3.2.2.3离散实现
文献［18］中给出了离散余弦变换(Discrete Cosine Transform,DCT)和离散正弦变换(Discrete Sine Transform,DST)核矩阵。
(1) 四种类型的离散余弦变换核矩阵定义
DCTⅠ

CⅠN+1=2NkmkncosmnπN，m，n=0,1,…,N(3.124)

DCTⅡ

CⅡN=2Nkmcosm(n+1/2)πN，m，n=0，1,…,N-1(3.125)

DCTⅢ

CⅢN=2Nkncos(m+1/2)nπN，m，n=0，1,…,N-1(3.126)

DCTⅣ

CⅣN=2Ncos(m+1/2)(n+1/2)πN，m，n=0，1,…,N-1(3.127)

式中，km和kn定义为

km=1/2，m=0,m=N

1，m≠0,m≠N，kn=1/2，n=0,n=N

1，n≠0,n≠N
(3.128)

(2) 四种类型的离散正弦变换核矩阵定义
DSTⅠ

SⅠN-1=2NsinmnπN，m,n=1,2,…,N-1(3.129)

DSTⅡ

SⅡN=2Nkmsinm(n-1/2)πN，m,n=0,1,…,N(3.130)

DSTⅢ

SⅢN=2Nknsin(m-1/2)nπN，m,n=0,1,…,N(3.131)

DSTⅣ

SⅣN=2Nsin(m-1/2)(n-1/2)πN，m,n=0,1,…,N(3.132)

式中，km和kn的定义与式(3.128)中相同。核DCTⅠ和DSTⅠ有对称的结构并以2为周期。周期性意味着重复应用DCTⅠ和DSTⅠ将得到原始的序列。DCTⅣ与DCTⅠ有相同对称性和周期性，但DCTⅡ和DCTⅢ互为逆算子并不具周期性。DCTⅠ和DSTⅠ将被用于导出离散分数阶余弦变换(Discrete Fractional Cosine Transform,DFRCT)和离散分数阶余弦变换(Discrete Fractional Sine Transform,DFRST)： 


CⅠN=2N-1·

1212…1212

12cosπN-1…cos(N-2)πN-112cos(N-1)πN-1

︙︙︙︙

12cos(N-2)πN-1…cos(N-2)2πN-112cos(N-2)(N-1)πN-1

1212cos(N-1)πN-1…12cos(N-2)(N-1)πN-112cos(N-1)2πN-1
(3.133)

SⅠN=2N+1·

sinπN+1sin2πN+1…sin(N-1)πN+1sinNπN+1

sin2πN+1sin4πN+1…sin2(N-1)πN+1sin2NπN+1

︙︙︙︙

sin(N-1)πN+1sin2(N-1)πN+1…sin(N-1)2πN+1sinN(N-1)πN+1

sinNπN+1sin2NπN+1…sinN(N-1)πN+1sinN2πN+1(3.134)


(3) 离散傅里叶变换(Discrete Fourier Transform,DFT)核矩阵
类似于其连续变换，离散傅里叶变换可看作离散信号的π/2旋转。离散傅里叶变换核矩阵定义如下

FN=1N11…11

1W1N…WN-2NWN-1N

︙︙︙︙

1WN-2N…W(N-2)2NW(N-1)(N-2)N

1WN-1N…W(N-1)(N-2)NW(N-1)2N
(3.135)

其中，WN=e-j(2π/N)。
(4) 离散分数阶傅里叶变换核矩阵
根据离散傅里叶变换核矩阵，N点的离散分数阶傅里叶变换(Discrete Fractional Fourier Transform,DFRFT)核矩阵由下式计算

FN,α=VND2α/πNVTN=VN10

e-jα



0e-j(N-1)αVTN(3.136)

其中，VN=［v0|v1|…|vN-1］，vk是k阶离散傅里叶变换Hermite特征向量，α表示变换在时频平面上旋转的角度。当α=0时，FN,α是单位算子。
DCTⅠ和DSTⅠ的特征向量可由离散傅里叶变换的特征向量得到。
① 如果v=［v0,v1,…,vN-2,vN-1,vN-2,…,v1］T是(2N-2)点离散傅里叶变换核矩阵的偶特征向量，F2N-2v=λv   (λ=1,-1)，那么，

v^=v0,2v1,…,2vN-2,2vN-1T(3.137)

为N点DCTⅠ核矩阵的特征向量，其中λ是相应的特征值。

CⅠNv^=λv^(3.138)

② 如果v=［0,v1,v2,…,vN,0,-vN,-vN-1,…,-v1］T是2(N+1)点离散傅里叶变换核矩阵的奇特征向量，F2N+2v=λv   (λ=j,-j)，那么，

v～=2v0,v1,…,vNT(3.139)

为N点DSTⅠ核矩阵的特征向量，其中jλ是相应的特征值。

SⅠNv～=jλv～(3.140)

所有离散傅里叶变换、离散余弦变换和离散正弦变换变换的核矩阵都具有无限的特征值。通过引入一个新的矩阵S就能非常简便地计算出了一组完整的实数离散傅里叶变换特征向量［1921］

S=210…01

12cos(ω)1…00

012cos(2ω)…00

︙︙︙︙︙

000…2cos［(N-2)ω］1

100…12cos［(N-1)ω］(3.141)

这组特殊的特征向量组成了连续HermiteGaussian函数的离散近似，称为离散傅里叶变换Hermite特征向量。对于离散分数阶余弦变换核矩阵，特征向量v^k的特征值为e-jkα(k为偶数)。这样的赋值将会使离散余弦变换核中α=π/2。类似地，在离散分数阶傅里叶变换中，N点离散分数阶余弦变换核矩阵可定义为


CN,α=V^ND^2α/πNV^TN(3.142)

=V^N10

e-2jα



0e-j2(N-1)αV^TN(3.143)


当V^N=v^0v^1…v^2N-2，v^k是从k阶离散傅里叶变换 Hermite特征向量由式(3.137)得到的DCTⅠ特征向量。当α=π/2，离散分数阶余弦变换将成为常规的DCTⅠ。当α=0，CN,α为单位矩阵。用参数α构造N点离散分数阶余弦变换核矩阵的步骤可总结如下： 
① 计算Mc点离散傅里叶变换 Hermite偶特征向量，Mc=2(N-1)。
② 用式(3.137)从DFTⅠ Hermite 偶特征向量计算DCTⅠ特征向量。
③ 用下式确定离散分数阶余弦变换的核矩阵： 

CN,α=V^ND^2α/πNV^TN(3.144)

其中，V^N=v^0v^1…v^Mc-2。v^k是从k阶离散傅里叶变换Hermite特征向量由式(3.137)得到的DCTⅠ特征向量。
和离散分数阶余弦变换的情况相似，离散分数阶正弦变换的推导也建立在离散分数阶傅里叶变换上。特征向量vk(k为奇数)赋值为特征值e-j(k-1)α。因此，N点离散分数阶正弦变换核矩阵定义为


SN,α=V～ND～2α/πNV～TN(3.145)

=V～N10

e-2jα



0e-j2(N-1)αV～TN(3.146)


其中，V～N=v～0v～1…v～2N-2。v～k是从k阶离散傅里叶变换Hermite特征向量由式(3.139)得到的DSTⅠ特征向量。
以上所述离散分数阶正弦变换核矩阵，当α=π/2，它将简化为一个DSTⅠ核矩阵，对于α=0，它将变为单位矩阵。用参数α计算N点离散分数阶正弦变换核矩阵的步骤总结如下： 
① 计算MS点离散傅里叶变换 Hermite奇特征向量，MS=2(N+1)。
② 从DFT Hermite奇特征向量用式(3.139)计算DSTⅠ特征值。
③ 用下式确定离散分数阶正弦变换的核矩阵

SN,α=V～ND～2α/πNV～TN(3.147)

其中，V～N=v～0v～1…v～Ms-2。v～k是从k阶DFT Hermite特征向量由式(3.139)得到的DSTⅠ特征向量。
由于用来精确计算离散分数阶傅里叶变换、离散分数阶余弦变换和离散分数阶正弦变换矩阵的积的快速算法尚未被推导出，它们的计算一般需要O(N2)次常规的矩阵乘法。

3.2.3分数阶Hartley变换
3.2.3.1定义
Hartley变换与傅里叶变换有如下关系


OHartley(g(t))=Ocos(g(t))+Osin(g(t))

=1+j2·G(ω)+1-j2·G(-ω)
(3.148)

对于m取任何非负整数，exp(-t2/2)Hm(t)是Hartley变换的特征函数，其相应的特征值记作λH(m)。当m为偶数时，λH(m)=(-j)m； 当m为奇数时，λH(m)=(-j)m-1。类似于分数阶傅里叶变换，可通过如下过程推导出分数阶Hartley变换(Fractional Hartley Transform,FRHT)。
(1) g(t)→{am|m=0,1,2,3,…}
其中，am=C-1m∫+∞-∞g(t)Hm(t)e-t2/2dt

Cm=∫+∞-∞H2m(t)e-t2dt(3.149)

(2) dm=λαH(m)am(3.150)
(3) GαH(s)=∑+∞m=0dme-s2/2Hm(s)(3.151)

其中，λαH(m)是Hartley变换分数幂的特征值。当m为偶数时，λαH(m)=exp(-jmαπ/2)； 当m为奇数时，λαH(m)=exp(-j(m-1)απ/2)。显然，存在如下关系

OαH(f(t))=OαC(f(t))+OαS(f(t))(3.152)

因此，积分形式的分数阶Hartley变换为［22］


GαH(s)=OαH(g(t))=1-jcot2πej(s2/2)cot∫+∞-∞ej(t2/2)cot·


(1-jej)cas(stcsc)+(1+jej)
cas(-stcsc)2g(t)dt
(3.153)

其中，cas(·)=cos(·)+sin(·)。分数阶Hartley变换与分数阶余弦变换、分数阶正弦变换和分数阶傅里叶变换有如下关系


GαH(s)=GαC(s)+GαS(s)=1+e(jαπ/2)2
GαF(s)+1-e(jαπ/2)2GαF(-s)(3.154)


由上式可知，分数阶Hartley变换的周期为2，与分数阶傅里叶变换以4为周期不同。

OαH(g(t))=Oα+2H(g(t))(3.155)

3.2.3.2性质
分数阶Hartley变换的性质与分数阶余弦变换类似，分数阶余弦变换的性质可由正则余弦变换的性质取{a,b,c,d}={cot,1,-1,0}时得到。
3.2.3.3离散实现
离散傅里叶变换矩阵F的元素定义如下

Fnk=1Ncos2πknN-jsin2πknN，0≤n,k≤N-1(3.156)

令ω=2π/N，并令

S=210…01

12cos(ω)1…00

012cos(2ω)…00

︙︙︙︙︙

000…2cos［(N-2)ω］1

100…12cos［(N-1)ω］

(3.157)

则S的特征向量构成F的特征向量集，并有［21］

FS=SF

注意到S是实对称矩阵，则它的特征向量为实数并且彼此正交。
N×N DHT矩阵H的元素定义如下

Hnk=1Ncos2πknN+sin2πknN,0≤n,k≤N-1(3.158)

为了进一步讨论方便，定义


Frnk=1Ncos2πknN

Fink=1Nsin2πknN0≤n,k≤N-1
(3.159)

则矩阵H和F可重写为

H=Fr+Fi

F=Fr-jFi(3.160)

由于矩阵Fr和Fi彼此正交，则Fr具有非零特征值的特征向量是Fi特征值为0的特征向量； 同样，Fi具有非零特征值的特征向量是Fr特征值为0的特征向量。基于这一事实，Fr、Fi和H的特征分解可表示如下



Fr=U1U2U3U4I1000

0-I200

0000

0000U1U2U3U4T(3.161)

Fi=U1U2U3U40000

0000

00I30

000-I4U1U2U3U4T(3.162)

H=U1U2U3U4I1000

0-I200

00I30

000-I4U1U2U3U4T(3.163)



其中，Ii是Ni×Ni的单位矩阵。
N=4m时，N1=m+1，N2=m，N3=m，N4=m-1；
N=4m+1 时，N1=m+1，N2=m，N3=m，N4=m；
N=4m+2 时，N1=m+1，N2=m+1，N3=m，N4=m；
N=4m+3 时，N1=m+1，N2=m+1，N3=m+1，N4=m。
矩阵Ui定义如下： 
(1) U1由矩阵S的特征向量v构造，并满足Frv=v； 
(2) U2由矩阵S的特征向量v构造，并满足Frv=-v； 
(3) U3由矩阵S的特征向量v构造，并满足Fiv=v； 
(4) U4由矩阵S的特征向量v构造，并满足Fiv=-v。
设数据向量为x，则它的分数阶Hartley变换yτ定义为

yτ=Hτx(3.164)

当幂τ为1时，离散分数阶Hartley变换变为Hartley变换。当τ=0时，显然y0=x。因为方程Hτ1+τ2x=Hτ1Hτ2x成立，所以角度可加性可以满足。现在的问题是如何计算矩阵Hτ。矩阵Hτ可以由矩阵H特征值的τ次幂得到

Hτ=U1U2U3U4Iτ1000

0(-I2)τ00

00Iτ30

000(-I4)τU1U2U3U4T
(3.165)

在此分解式中有两个变量而使Hτ的计算结果不唯一。结果分别如下： 
(1) 因为以下两式成立：

1τ=e-j2kπτ

(-1)τ=e-j(2k+1)πτk∈Z

矩阵Iτi(i=1,3)和(-Ii)τ(i=2,4)不是唯一的，需加上限制条件才能去除不定性。
(2) 因为U1是由矩阵S的特征向量v构造，满足Frv=v，U1的任意两列向量可交换，因此存在N1!个矩阵U1。类似地，矩阵Ui(i=2,4)也有同样的问题。为了使Ui唯一，必须设定Ui的列向量的顺序。
下面给出一种简单的方法来排除不定性(1)和(2)。为了排除不定性(1)，可选择Iτi和(-Ii)τ如下


Iτi=e-j0πτ00…00

0e-j2πτ0…00

00e-j4πτ…00

︙︙︙︙︙

000…0e-j2(Ni-1)πτ(3.166)

(-Ii)τ=e-jπτ00…00

0e-j3πτ0…00

00e-j5πτ…00

︙︙︙︙︙

000…0e-j(2Ni-1)πτ(3.167)


另外，我们用以下的方法来安排矩阵Ui的列向量以排除不定性(2)。设uim和uin为矩阵Ui的列向量，则存在λm和λn使


Suim=λmuim

Suin=λnuin


加在uim和uin列向量上的约束条件为

λm＞λn，m＜n(3.168)

由于S的特征值不同，矩阵Ui可由以上限制条件唯一确定。
虽然排除不定性(1)和(2)的方法有很多，但上述方法是最简便的。最后，总结一下分数阶Hartley变换的计算过程： 
(1) 计算矩阵S的特征值和特征向量； 
(2) 用式(3.168)构建矩阵Ui(i=1,2,…,4)； 
(3) 用式(3.166)和式(3.105)计算Iτi(i=1,3)和(-Ii)τ(i=2,4)； 
(4) 用式(3.165)计算矩阵Hτ； 
(5) 用式(3.164)计算yτ=Hτx。
当数据向量x为实数，则它的分数阶Hartley变换yτ是复数，当τ是整数时例外。另外，考虑到矩阵H的特征值是1或-1，易得出


H2mx=x(3.169)

H2m+1x=Hx(3.170)


3.2.4分数阶Gabor展开
3.2.4.1Gabor展开
传统Gabor展开用时间和频率的移位函数来表示信号，已被广泛应用于信号的时频分析。Gabor表示的基函数可通过对一个简单的窗函数平移和正弦调制得到，从而得到一个矩形网格状的时频平面结构(图3.21)。时间连续信号x(t)的Gabor展开为

x(t)=∑+∞m=-∞∑+∞k=-∞am,kgm,k(t)(3.171)

这里基函数为

gm,k(t)=g(t-mT)ejΩkt(3.172)

T表示线性时移参数，Ω是频率采样间隔。综合窗函数g(t)归一化为单位能量。Gabor展开的存在性、唯一性、收敛性和数字稳定性取决于参数T和Ω的选择： 临界采样情况是ΩT=2π； ΩT＜2π称为过采样，会导致Gabor系数冗余； ΩT＞2π称为欠采样，会导致信息丢失。
通常一组时频移位窗函数hm,k(t)可构成L2空间的一组非正交基。此时，由于不能利用内积映射，Gabor系数的计算将变得十分复杂。针对这一问题，有文献提出了一种利用称为双正交窗的辅助函数γ(t)的计算方法，Gabor系数am,k可由此得到

am,k=∫+∞-∞x(t)γm,k(t)dt(3.173)

这里的分析函数是

γm,k(t)=γ(t-mT)efΩkt(3.174)

把式(3.173)代入式(3.171)，就可得到这组基的完整性条件

∑+∞m=-∞∑+∞k=-∞gm,k(t)γm,k(t′)=δ(t-t′)(3.175)



图3.21Gabor展开的正交矩形坐标网格



在以上的完整性关系中，利用泊松求和公式可以在分析基和综合基之间产生一种等价却较为简单的双正交条件

2πΩ∑+∞m=-∞g(t-mT)γt-m+k2πΩTT=δk(3.176)

其中，k=0,±1,±2,±3,…； 因子2π/ΩT为过采样的一个测度。
近年来，非矩形时频网格结构中的Gabor展开引起了广泛的关注。与传统的Gabor展开相比，非矩形时频网格更适于非平稳信号的时频分析。

3.2.4.2分数阶Gabor展开
利用具有线性瞬时频率的基函数代替传统的正弦Gabor核，可定义分数阶Gabor展开［23］。信号x(t)的Gabor展开定义为

x(t)=∑+∞m=-∞∑+∞k=-∞am,k,αgm,k,α(t)(3.177)

其中，综合基函数gm,k,α(t)为

gm,k,α(t)=g(t-mT)Wα,k(t)(3.178)

若定义分数阶核为

Wα,k(t)=expj-12t2+kΩsinα2cotα+kΩt(3.179)

其中，Ω和T分别为频率和时间采样间隔，且0≤α≤2π。利用上述的分数阶核生成的基函数gm,k,α(t)的瞬时频率是线性变化的，这种分数阶基函数就可产生如图3.22所示平行四边形状的时频采样网格。


图3.22分数阶Gabor变换的时频坐标网格



分数阶Gabor系数am,k,α可由下式计算得到

am,k,α=∫+∞-∞x(t)γm,k,α(t)dt(3.180)

其中，分析函数为

γm,k,α(t)=γ(t-mT)Wα,k(t)(3.181)

相对于gm,k,α(t)，它们是双正交的。当α=π/2，式(3.177)就成为式(3.171)中的Gabor展开。这样传统Gabor展开就可看作分数阶Gabor展开的特例。接下来我们研究分数阶Gabor展开的完整性和双正交条件。
1. 分数阶基的完整性
把式(3.179)代入式(3.177)就可得到分数阶Gabor展开的完整性

∑+∞m=-∞∑+∞k=-∞gm,k,α(t)γm,k,α(t′)=δ(t-t′)(3.182)

将式(3.178)和式(3.181)中的分析和综合函数代入上式可得


∑+∞m=-∞∑+∞k=-∞g(t-mT)γ(t′-mT)·

expj12t′2-t2cotα-kΩt-t′=δ(t-t′)

(3.183)

显然，当α=π/2，上述的条件简化为经典Gabor展开的完整性条件。
2. 分数阶双正交条件
分析分数阶分析与综合函数基必须满足的双正交性条件。式(3.183)表示的完整性条件也可写作


∑+∞m=-∞g(t-mT)γ(t′-mT)expj12t′2-t2cotα·

∑+∞k=-∞expjkΩt-t′=δ(t-t′)

(3.184)

应用泊松求和公式对式(3.184)中的参数k求和可得

∑kexpjkΩt-t′=2πΩ∑kδt-t′-k2πΩ(3.185)

将式(3.185)代入式(3.184)得到


2πΩ∑+∞m=-∞g(t-mT)γt′-m+2kπΩTTexpj12t-2kπΩ2-t2cotα

∑+∞k=-∞δt-t′-2kπΩ=δt-t′


从上式可以得到分数阶双正交条件是

2πΩ∑+∞m=-∞g(t-mT)γt′-m+2kπΩTTexpj2kπΩkπΩ-tcotα=δk(3.186)

其中，m,k=0,±1,±2,±3,…。注意上式的指数项是由角度为α的分数阶基产生的，并且当α=π/2时，我们得到的双正交性条件就是式(3.186)表示的Gabor展开的双正交性条件。这表明分数阶Gabor展开是经典Gabor展开在非矩形时频网格中的推广。分析窗γ(t)函数可通过解由式(3.185)得到的线性方程得到。利用分析函数基γm,k,α(t)和式(3.180)即可计算分数阶Gabor系数am,k,α。




教学视频



3.2.5短时分数阶傅里叶变换
短时傅里叶变换(ShortTime Fourier Transform,STFT)通过运用一个合适的窗函数g(t)，能够更好地获得信号x(t)频率成分的时间定位。

STx(t,f)=∫+∞-∞xt+t0gt0exp-j2πt0fdt0(3.187)

当然，对于一个纯正弦信号，我们需要一个较宽的窗，而对于类似于冲激函数的脉冲信号，则需要使用较窄的窗。这一规则也同样适用于分析延伸非常广的信号和非常窄的信号。因此，如果我们已知信号的形状，就可以对窗宽作相应的调整。假定当前信号的最小信号宽度并不对应于时域和频域，如图3.23所示。那么我们可以看到，通过对相位平面的坐标旋转，可以得到信号的最优表示(比如最小信号宽度)。


图3.23时频面上的基准轴既不是时间轴也不是频率轴的信号示意图



为了在新的坐标系中表示信号，我们将用到如下性质： 信号的分数阶傅里叶变换相当于时频面的旋转。根据第2章内容，函数x(t)的分数阶傅里叶变换定义如下

Xα(u)=∫+∞-∞Kα(t,u)x(t)dt(3.188)

其中，核函数Kα(t,u)定义如下

Kα(t,u)=expjα/2jsinαexpjπt2+u2cosα-2tusinα(3.189)

通过旋转关系式

t
f=cosα-sinα

sinαcosαu

v(3.190)

可以进一步得出分数阶傅里叶变换核函数的关系式


Kαt0,u-u0=expj2πu0vexp-jπuv

=K-αu0,t-t0expj2πt0fexp-jπtf
(3.191)

定义信号x(t)的α阶短时分数阶傅里叶变换(ShortTime Fractional Fourier Transform,STFRFT)STαx(u,v)为其分数阶傅里叶变换Xα(u)与窗函数g(u)的卷积的傅里叶变换


STαx(u,v)=STXα(u,v)

=∫+∞-∞Xαu+u0gu0exp-j2πu0vdu0
(3.192)

由式(3.191)可得


exp-jπuv∫+∞-∞Xαu+u0gu0exp-j2πu0vdu0

=exp-jπtf∫+∞-∞xt+t0G-αt0exp-j2πt0fdt0

(3.193)

用初始窗函数g(t)的分数阶傅里叶变换作为新的窗函数，再结合式(3.190)中的坐标系转换，我们能够像计算普通信号x(t)的短时傅里叶变换一样直接计算信号的α阶短时傅里叶变换STαx(u,v)


STαx(u,v)=expjπuv-tf
∫+∞-∞xt+t0G-αu0exp-j2πt0fdt0

(3.194)

其中，u,v和t,f的关系如式(3.190)所示。
接下来考虑利用高斯函数g(t)=exp-πct2作为窗函数。它的分数阶傅里叶变换表示式是

G(u)=expjα/2cosα+jcsinαexp-πcu21+tan2α-jc-c-1tanα1+c2tan2α(3.195)

当c=1时，该高斯函数是分数阶傅里叶变换的特征函数，且利用该函数在分数阶傅里叶域滤波对应于相应短时傅里叶变换表达式的旋转。我们先来看一个简单的线性调频信号

x=expjπpt2+j2πqt(3.196)

假定高斯窗exp-πct2是纯谐波信号expj2πqt的最优滤波窗，为了找到高斯窗函数针对该线性调频信号的最优滤波参数，我们需要在分数阶傅里叶域上来求解。根据文献［24］，有

X～α(u)=Xαu-qsinαexpj2πqcosαu-qsinα/2

其中，x～(t)=x(t)expj2πqt。那么把上式应用到chirp信号expj2πpt2的分数阶傅里叶变换expjα/2cosα1+ptanαexpjπu2p-tanα1+ptanα上，可以得到： 当分数阶傅里叶变换角度α=arctanp时，线性调频信号的α阶分数阶傅里叶变换为一纯谐波信号，而当分数阶傅里叶变换角度α=arctanp+π/2时，在qsinα处会出现一个冲激脉冲。重新回到时域，我们可以得到式(3.196)所示线性调频信号所对应的最优高斯窗函数如下

F-α［g］(u)=Fαg(u)=Aexp-πcu21+p2-jc-c-1p1+c2p2(3.197)

该式是把α=arctanp代入式(3.195)得到的。
实际上，需要处理的信号往往不会是理想的纯线性调频信号，有时还会存在多个分量。然而，只要瞬时频率值在时频面上的某一线段(我们将该线段作为基准轴线)方向上变化缓慢，我们就可以找到信号相对较为集中或较为分散的分数阶傅里叶域。
1. 利用分数阶傅里叶矩估计信号宽度
我们知道，信号的时频域宽度可以由它的二阶中心矩来估计。根据第2章介绍的分数阶傅里叶变换矩的相关内容，可以知道，分数阶傅里叶域的信号宽度也同样可以由其二阶分数阶傅里叶变换中心矩来估计。将二阶分数阶傅里叶变换中心矩pα定义如下

pα=∫+∞-∞Rαx(t)2t-mα2dt=wx-mα2(3.198)

其中，mα=∫+∞-∞Rαx(t)2tdt为一阶矩，wα=∫+∞-∞Rαx(t)2t2dt是二阶矩。在分数阶傅里叶域中，角度α的分数阶傅里叶变换一阶矩mα可由下面的关系式推得

mα=m0cosα+mπ/2sinα(3.199)

其中，m0和mπ/2分别是其时域和频域的一阶矩。同样，角度α的分数阶傅里叶变换二阶矩wα可以由其他3个二阶矩wβ、wγ和wμ表示。需要注意的是，γ、β和μ必须是3个不同的角度，且其中任意两者之差不能为π。我们选择3个不同角度的二阶矩w0、wπ/2和wπ/4，直接利用文献［25］中的结果，有

wα=w0cos2α+wπ/2sin2α+wπ/2-w0+wπ/2/2sin2α(3.200)

综合考虑式(3.198)~式(3.200)可得，只需知道3个分数阶傅里叶变换的功率谱，即可以确定所有角度的二阶中心矩pα，也就得到了相应角度分数阶傅里叶域的信号宽度


pα=w0-m0cos2α+wπ/2-w2π/2sin2α+wπ/4-m0mπ/2-w0+wπ/2/2sin2α

=p0cos2α+pπ/2sin2α+wπ/4-m0mπ/2-w0+wπ/2/2sin2α

(3.201)

我们可以通过pα的导数来研究信号的极限宽度对应的分数阶傅里叶域。通过式(3.201)，易得pα的一阶导数如下

dpαdα=pπ/2-p0sin2α+2wπ/4-m0mπ/2-w0+wπ/2cos2α

令上式等于0，就能求得极限宽度所对应的角度αe

tan2αe=2wπ/4-m0mπ/2-w0+wπ/2p0-pπ/2(3.202)

由于分数阶傅里叶变换是以2π为周期，且满足Fα+π［x(t)］=Fαx-t，因此，信号宽度在α∈0,π内必有一对极大值与极小值。由pα在α=αe时的二阶导数d2pαdα2α=αe=2pπ/2-p0cos2αe可得，若pπ/2-p0与cos2αe同号，此时所得信号的极值宽度为最小宽度，否则为最大宽度。因此，只需要通过3个分数阶傅里叶功率谱值，就可求得使信号有最佳聚焦或最大延伸的最优分数阶傅里叶变换角度。
2. 基于短时分数阶傅里叶变换的时变滤波
时域无重叠的信号分量可以通过选择开关对其进行分离，频域无重叠的信号分量可以通过带通滤波器对其进行分离。对于多分量的非平稳信号，如雷达领域内的多分量LFM回波信号，其在时域和频域都是严重混叠的，且信号能量分散，难以设计合适的滤波器将其分离。但在时频域内，不同分量的信号由于其时频变化特性不同，其时频谱在大部分范围内并不会重叠，且能量集中，这启发我们可以从时频平面将各个信号分离。我们知道，多分量时变信号在时频分布中表现为多个能量脊，而从多分量信号中分离出各分量信号的问题就可以定义为从信号的时频分布中提取出各个分量信号的时频表征，进而通过其时频表征还原出时域信号［26］。这种采用时频表征在时频分布上对信号进行处理，然后通过逆变换重建时域信号的操作定义为时变滤波(TimeVarying Filtering,TVF)。
时变滤波要求所用的时频变化具有可逆性质，常用的有短时傅里叶变换、短时分数傅里叶变换、小波变换谱、同步压缩小波变换谱等。以短时傅里叶变换为例，多分量信号的时变滤波定义为

z～i(n)=1N∑N-1k=0Fz［n,k］Bi［n,k］(3.203)

其中，Fz［n,k］为信号的短时傅里叶变换分布，z～i(n)为需要重建的时域信号，Bi［n,k］表示第i个信号分量的时变滤波器。简单有效的时变滤波器定义为： 在存在信号分量的时频域能量聚集处 Bi［n,k］=1，其他区域Bi［n,k］=0，这样的滤波器对于在时频域内不混叠的多分量信号分离非常有效。而对于在时频平面内混叠的信号分量，就要利用不同信号分量在时频平面内的能量聚焦度来设计滤波器，如

H(n,k)=1,STFTx(n,k)≥Th

0，其他(3.204)

式中，Th为依据经验选取的阈值。该时变滤波器要求待分离的多分量信号在时频平面内的能量聚焦度有显著的差异。
众所周知，窗长是短时傅里叶变换和短时分数阶傅里叶变换中的一个重要参数，不同的窗长可使其呈现出不同的时频聚集特性。宽窗可以提供良好的频率分辨率，但会导致时间分辨率变差； 反之，窄的窗函数可以提供好的时间分辨率，但会降低频率分辨率。因此，为了获得更好的能谱集中，频率慢变信号的分析通常需要一个宽窗，而频率快变信号的分析则需要一个窄窗［27］。由于短时分数阶傅里叶变换结合了短时傅里叶变换和分数阶傅里叶变换的优点，因此它可以很好地描述时变信号。尤其是针对快变的频率调制信号，短时分数阶傅里叶变换可以比短时傅里叶变换提供更好的时频分辨能力［2829］。与短时傅里叶变换时频平面内的矩形支撑不同，短时分数阶傅里叶变换对时频平面的划分为平行四边形，而每个平行四边形网格称为一个时间分数域频率单元，如图3.24所示，其两个边长的乘积的倒数为短时分数阶傅里叶变换的二维分辨率。当且仅当使用高斯窗时，时宽带宽积达到最小，短时分数阶傅里叶变换的二维分辨率达到最佳。因此，使用高斯窗的短时分数阶傅里叶变换也称为最优短时分数阶傅里叶变换或高斯短时分数阶傅里叶变换。此时的高斯窗函数形式如下［30］

g(t)=πTx|sinα|Bx,p-1/4exp-Bx,pt22Tx|sinα|(3.205)

式中，Tx和Bx,p分别是信号的时宽和信号p阶短时分数阶傅里叶域的分数域带宽。高斯窗函数的方差为σ2=Txsinα/Bx,p，根据这些参数可以确定最佳时间窗长度。


图3.24短时傅里叶变换和短时分数阶傅里叶变换的时频分割表示



(a) 短时傅里叶变换； (b) 短时分数阶傅里叶变换


与短时傅里叶变换不同的是，除了所用窗函数的长度，短时分数阶傅里叶变换的变换阶次也会影响信号在时频域的能谱集中度。由于分数阶傅里叶变换的时频旋转变换能力，一个信号经过分数阶傅里叶变换后，其频谱宽度可以被展宽或压缩，从而影响其频谱聚焦度。例如，一个调频率为μ0的LFM信号，当对其进行p0=2arccot(-μ0)/π阶次的分数阶傅里叶变换时，该LFM信号在分数域变成一个冲激信号，此时它的功率谱高度集中，幅度值达到最大，此时的阶次p0称作匹配阶次。即使待分析的时变信号不是严格的LFM信号，可以通过采用短时分数阶傅里叶变换将其在窗函数里的每一段成分近似为LFM，其对应的阶次称为准匹配阶次。通过对信号执行准匹配阶次下的短时分数阶傅里叶变换，可以确保其在短时分数阶傅里叶域的时频表征具有高的能谱集中。相反，所选的阶次与信号的匹配或准匹配阶次的差距越大，信号的分数阶傅里叶谱就展宽得越严重，从而造成其时频表征内的能量严重弥散。基于以上分析，对于一个含有两个成分的信号，可以通过调整变换阶次来改变不同信号成分在时频平面内的能谱聚焦度，这样有助于提高特定信号成分稀疏表征的可能性，有利于时变滤波器的设计。如选定的阶次接近信号成分1的匹配阶次，而远离成分2的匹配阶次时，那么成分1的能谱聚集度就会大于成分2的，这时信号成分1相对于成分2就更容易被稀疏表征。合适的阶次可以通过对信号做不同阶次的分数阶傅里叶变换遍历搜索得到，当某个阶次下对应的不同信号成分的频谱幅度差异最大，该阶次就是所需要的合适阶次。当合适的阶次得到后，可以通过式(3.205)得到短时分数阶傅里叶变换最优窗长的参考值。为了更好地说明以上结论，下面给出一个例子来说明基于短时分数阶傅里叶变换的稀疏表征和时变滤波。
考虑一个含有两个分量的信号，其表达式为

x(t)=ax1(t)+x2(t)=aej2π(8t+2arctan(t-5)2)+ej2π(5t+10sin(3.5t))(3.206)



图3.25给定信号不同阶次的分数阶傅里叶变换结果


(a) 时域； (b) p=1； (c) p=0.99； (d) p=0.1


式中，第一个指数项代表频率慢变的时变信号； 第二个指数项表示频率快变的信号； a是两个成分的幅度比，此处设置为1，即两个成分的幅度相同。以100Hz的采样频率对信号进行采样，采样总点数为1400点。然后对这个信号进行分数阶傅里叶变换，图3.25给出了信号时域和不同阶次下的分数阶傅里叶变换结果，可以看出两个成分在时域严重混叠。图3.25(b)为信号在阶次为1时的分数阶傅里叶变换结果，等同于信号的傅里叶变换，高的峰值是慢变信号的频谱，分布在其两侧较低且宽的谱线是快变信号的频谱。观察图3.25(c)信号0.99阶分数阶傅里叶变换的结果，可以看到慢变信号的频谱强度没有明显的变化，而分布在两侧的快变信号的频谱强度有明显的下降，两个成分的频谱强度差异明显增加。从滤波的角度，两个成分的强度差异越大越有利于强的成分的提取。因此，通过调整阶次改变不同成分的频谱差异有助于特定成分的恢复，而此时的阶次0.99也就是慢变信号的准匹配阶次。同样，从图3.25(d)可以发现慢变信号的频谱展宽，幅度严重下降，而快变信号的频谱出现了一定的聚焦，有8个明显的峰值出现，这是在频域所没有的。而且从后续信号的时频表征可以看到这8个峰值对应快变微多普勒信号的8个周期。尽管如此，快变信号的谱和慢变信号的谱还是严重混叠在一起，不能够将二者分离的。


图3.26给定信号的短时傅里叶变换和短时分数阶傅里叶变换结果



(a) 宽窗短时傅里叶变换； (b) 窄窗短时傅里叶变换； 

(c) 0.99阶宽窗短时分数阶傅里叶变换； (d) 0.1阶窄窗短时分数阶傅里叶变换



然后，对上述信号分别进行宽窗和窄窗的短时傅里叶变换，以及0.99阶的宽窗短时分数阶傅里叶变换和0.1阶的窄窗短时分数阶傅里叶变换。其中，短时傅里叶变换中宽窗和窄窗分别为256点和16点的汉明窗。短时分数阶傅里叶变换的宽窗和窄窗的长度由前面提到的最佳参考窗长求法计算得到，分别为237点和19点。将其折算成2的N次幂便于计算，分别取256和16。图3.26展示了变换后的时频图。从图3.26(a)和图3.26(c)的结果可以看出无论是短时傅里叶变换还是短时分数阶傅里叶变换，当使用较宽的窗函数时，慢变信号的能谱集中度非常高，快变信号的谱几乎不可见，这时相对于频率快变成分可以认为频率慢变成分是稀疏的。且由于受到分数阶变换阶次p的影响，图3.26(c)中信号的稀疏度略高于图3.26(a)中的。相反，当短时傅里叶变换和短时分数阶傅里叶变换的窗函数都是窄窗时，快变信号的能谱集中度明显增加，可以清晰地看到其周期性的调制规律。而且使用短时分数阶傅里叶变换后快变信号的能谱集中度和频率分辨率明显都要优于使用短时傅里叶变换的结果，如图3.26(b)和图3.26(d)所示，这样的结果对快变信号的稀疏表征和分离是有利
实际上，在图3.26(d)中，尽管快变信号的能谱集中度增强了，但其强度与慢变信号的能谱差异并不大，所以此时的快变信号仍不能被很好地稀疏表征，难以设计合适的时变滤波器将其分离出。因此，为了能够更好地稀疏表征快变信号成分，除了要用到合适阶次的窄窗短时分数阶傅里叶变换外，还要利用CLEAN技术，该技术可以消除强的慢变信号对弱的快变信号的影响。用同样的例子说明以上情况，此时的两个信号成分的幅度比值a设置为2。图3.27给出了该信号在使用CLEAN消除强分量干扰前后的0.1阶窄窗短时分数阶傅里叶变换的时频分析结果。从图3.27(a)中可以看到受到强分量的影响，即使采用了合适阶次的窄窗短时分数阶傅里叶变换来增强快变信号的能谱集中度，但其能谱强度还是比强分量的弱很多，这就很难对快变信号进行稀疏表征。反观图3.27(b)，当利用CLEAN消除强分量干扰后，慢变信号的谱几乎不可见，此时的快变信号可以被很好地稀疏表征，通过设计简单的阈值滤波器很容易就可将其提取出来。


图3.27给定信号在幅度比值为2时的0.1阶窄窗短时分数阶傅里叶变换结果



(a) CLEAN前； (b) CLEAN后






教学视频








教学视频



3.2.6分数阶模糊函数、分数阶Wigner分布
早期时频分析的许多发展都是由量子力学中使用的数学方法引导的［31］，其中有些方法的出发点是把一个合适的Hermite算子和一个诸如时间或频率的抽象物理变量联系起来。例如，对一个时域信号s(t)进行运算，那么Hermite时间算子T和Hermite频率算子O可以分别定义为［31］

T(s)(t)=ts(t),O(s)(t)=-j2πddts(t)(3.207)

任何一个Hermite算子都可以构成信号空间上的一个完备正交基［31］。因此，信号在这些特征函数上的投影自然地定义了一种信号变换，弄清楚这一点对于我们要考虑的信号来说是最为基础的。例如，Hermite频率算子O的特征函数是复指数型的，uF(t,f)=ej2πft，而由它们定义的信号变换就是傅里叶变换。

SF(f)=〈s,uF〉=∫s(t)e-j2πftdt(3.208)

其中，〈，〉表示内积，定义为∫g(t)h*(t)dt。对于任意变量的信号变换，比如a，可以类似地通过把信号分解到Hermite算子A和变量a的特征函数uA(t,a)上来定义

SA(a)=〈s,uA〉=∫s(t)u*A(t,a)dt(3.209)

参数化的酉算子也可以一种等价和一致的方式来表示物理变量［32］。例如，时移算子Tτ和频移算子Ov可以分别定义为

Tτ(s)(t)=s(t-τ)，Ov(s)(t)=s(t)ej2πvt(3.210) 

其中，参数τ,v∈R表示时移和频移值的大小。
接下来，我们回顾前述分数阶位移算子的概念，它的定义为

Rρ(s)(t)=s(t-ρcos)e-j2π(ρ2/2)cossin+j2πtρsin(3.211)

其中，Rρ将信号s(t)在时频平面上沿着与时间轴成角的方向做大小为ρ的移位。因此，Rρ将单位时移算子Tτ(=0)和单位频率算子Ov(=π/2)的概念推广到了时频平面的任意方向。也就是说，R0ρ(s)(t)=Tρ(s)(t)，Rπ/2ρ(s)(t)=Oρ(s)(t)。此外，分数阶位移算子Rρ也可以通过时移算子Tτ和频移算子Ov来表示： 

Rρ(s)(t)=e-jπρ2cossinOρsinTρcos(s)(t)(3.212)



图3.28时频图上分数阶变量(r)的表示


因此，分数阶位移算子也可理解为一种特殊的联合时频移位。时移和频移算子是时间变量t和频率变量f的一种表示，类似地，分数阶位移算子可以认为是时间轴和频率轴之间的某条坐标轴相关的“分数阶”变量的表示，该变量可以用r来表示(图3.28)。

根据文献［33］中提出的理论，我们可以导出分数阶变量r的Hermite算子表示为

R=cosT+sinO(3.213)

注意到R也可以简化为Hermite时间算子T和频率算子O的表示，它们分别对应式(3.207)中=0和=π/2的情况。利用Stone定理，酉算子和Hermite算子之间的关系可以表示为［7］

Rρ=e-j2πρR+π/2和R-π/2β=e-j2πβR(3.214)

傅里叶变换是信号的频域表示，与之类似，图3.28所示的分数阶傅里叶域r上的变换可以看作信号在式(3.209)所表示的Hermite算子R的特征函数上的分解。实际上我们可以看到，R的特征函数uR(t,r)形成了分数阶傅里叶变换的基函数［7］。因此，我们可以得到由Hermite分数阶算子R～定义的信号变换如下所示

SR～(r)=〈s,uR〉=Cejπr2cot∫s(t)ejπt2cot-j2πtrcscdt(3.215)

其中，C=1-jcot。当分数阶傅里叶变换旋转角度≠nπ(n为整数)时，式(3.215)与之是完全等价的。
利用上述酉算子和Hermite算子的知识，我们可采用特征函数法推导分数阶模糊函数和分数阶Wigner分布。
假定有我们感兴趣的两个变量a和b，它们的Hermite算子分别表示为A和B，那么我们就可以利用特征函数法来推导它们的联合信号表示。例如，Cohen类时频分布就可以利用式(3.207)中的Hermite时间算子T和频率算子O通过特征函数算子的方法推导出来［31］。最后，Cohen类的双线性、尺度变化时频表示的通用公式可以表示为

(PΨs)(t,f)=Ψτ,vsu+τ2s*u-τ2e-j2π(v(t-u)+τf)dudτdv(3.216)

现在我们考虑对应旋转角度为(图3.28)的时间变量t和分数阶域变量r的联合表示。我们在特征函数法中使用的第一个算子是Hermite时间算子T＝R0，第二个算子是Hermite分数阶算子R，其中0＜||≤π/2。如果=π/2，那么Rπ/2=O，这种时频表征是式(3.216)中Cohen类时频分布的简化。
步骤Ⅰ： 构造一个与函数ej2π(σt+τr)对应的特征函数算子。由于算子T和R的不可交换性(TR≠RT)，函数ej2π(σt+τr)与特征函数算子之间有多种不确定的关系。其中一种著名的特征函数算子称为Wely分类，它可以通过将表达式ej2π(σt+τr)中的变量t和r替换为Hermite算子T和R而得到，如下所示

M(τ,σ)=ej2πσT+j2πτR(3.217)

为了不用研究所有可能的分类情况，Cohen提出了一种简化的方法，称作核函数法。这种方法选定了一种单一的分类，例如式(3.217)所示的Wely分类，将它与一个核函数Ψ(τ,σ)相乘来得到所有可能的特征函数算子。因此，这个“归一化”的特征函数算子可以表示为

MΨ(τ,σ)=Ψ(τ,σ)M(τ,σ)=Ψ(τ,σ)ej2πσT+j2πτR(3.218)

步骤Ⅱ： 将特征函数作为特征函数算子的均值进行积分。与式(3.218)中归一化特征函数算子对应的归一化特征函数可以由下式计算而得［1］


MΨ(s)(τ,σ)＝∫s*(t)［MΨ(τ,σ)s］(t)dt

＝∫s*(t)Ψ(τ,σ)ej2πσT+j2πτRs(t)dt
(3.219)

我们称这种与单位核函数Ψ(τ,σ)＝1相关的特征函数为“分数阶模糊函数”(Fractional Ambiguity Function,FAF)，并用Ys(τ,σ)来表示

Ys(τ,σ)=〈ej2πσO+j2πτR〉(3.220)

式(3.220)中的分数阶模糊函数可以使用算子积分得到

Ys(τ,σ)=∫su+τ2sins*u-τ2sinej2π(σ+τcos)udu(3.221)

令上式的=π/2，则得到经典模糊函数(Ambiguity Function,AF)如下

Yπ/2s(τ,v)=∫st+τ2s*t-τ2ej2πvtdt(3.222)

模糊函数是与Cohen类时频分布核函数相关的特征函数，式(3.218)表示的传统模糊函数与式(3.221)表示的分数阶模糊函数有如下关系

Ys(τ,σ)=Yπ/2s(τsin,σ+τcos)(3.223)

我们可以从式(3.223)看到分数阶模糊函数将传统模糊函数的自变量“连接”起来，自变量使用三角函数进行了加权。
步骤Ⅲ： 最后，通过计算式(3.219)中归一化特征函数的二维傅里叶变换，我们可以得出广义的联合分数阶表示。这种通用的联合分数阶表示可以用式(3.221)中的分数阶模糊函数表示为

(PΨs)(t,r)=Ψ(τ,σ)Ys(τ,σ)e-j2π(σt+τr)dσdτ(3.224)

通过令Ψ(τ,σ)与Cohen类时频表示任一成员的核函数相等，我们可以找到一个对应的分数阶联合分数阶表示。例如，仍然用单位核函数这个特例，令式(3.224)中的Ψ(τ,σ)=1，并进行化简，可以得到“分数阶Wigner分布”


S,WD(τ,r)=Ys(τ,σ)e-j2π(σt+τr)dσdτ

=1sin∫st+τ2s*t-τ2e-j2π(τ/sin)(r-tcos)dτ

(3.225)

对于=π/2，R=O的情况，分数阶Wigner分布退化为Wigner分布。和分数阶模糊函数相似，分数阶Wigner分布与Wigner分布之间具有如下关系

S,WD(t,r)=1sinSπ/2,WDt,rsin-tcot(3.226)




教学视频








教学视频




3.2.7分数阶小波包变换
将分数阶傅里叶变换和小波包变换的概念相结合，可定义分数阶小波包变换。信号f(t)∈L2(LR)的短时傅里叶变换定义为12π∫+∞-∞e-jutg(t-τ)f(t)dt，这里g(t)是窗函数，小波变换中有类似的描述。相对于分解子波ψ，连续小波变换(Continuous Wavelet Transform,CWT)就变为1a∫+∞-∞ψt-baf(t)dt，这里a＞0并且ψ是完整化的，即L2模‖ψ‖=1。
小波包变换(Wavelet Packet Transform,WPT)可以看作短时傅里叶变换与连续小波变换的“混合”，即12πa∫+∞-∞e-jutψt-baf(t)dt。换句话说，小波包变换就是一个加窗信号的傅里叶变换，并且加窗函数是经a展宽和b平移的小波。
如前所述，分数阶傅里叶变换是传统傅里叶变换的推广。如果考虑一个信号在时间轴的表示及其傅里叶变换在频率轴的表示，可以将傅里叶变换看作信号逆时针旋转π/2角度。那么分数阶傅里叶变换就成为旋转角度非π/2整数倍的变换。
一个给定信号f(t)的分数阶小波包变换(Fractional Wavelet Packet Transform,FRWPT)WPTαf(u,a,b)为［31］

WPTαf(u,a,b)=1-jcotα2πeju22cotα∫+∞-∞ψt-baf(t)ejt22cotαejutcscαdt(3.227)

注意WPTα是时间频率比例函数，而且当α=π/2时，分数阶小波包变换与小波变换一致。分数阶小波包变换的计算步骤如下： 
(1) 与小波作乘积； 
(2) 与chirp信号作乘积； 
(3) 傅里叶变换； 
(4) 再与chirp信号作乘积； 
(5) 与复相位因子相乘。
如果定义ψa(t)=1aψta，那么‖ψα‖=‖ψ‖=1。回顾小波包变换的定义，可以把小波包变换看作短时傅里叶变换加窗函数ψα。因而我们可以得到分数阶小波包变换恒等式的分解形式为


∫+∞-∞∫+∞-∞dudbWPTαf(u,a,b)WPTα［g］(u,a,b)

=dudb∫Kα(u,t)ψα(t-b)f(t)dt∫Kα(u,t′)ψα(t′-b)g(t′)dt′

=dbdtψα(t-b)2f(t)g(t)

=‖ψa‖〈f,g〉=〈f,g〉

(3.228)

显然，分数阶小波包变换是一种新的时频分析工具，它与其他时频变换，如Wigner分布、模糊函数以及谱图等的联系及其在信号处理中的应用有待进一步研究。
3.3基于分数阶傅里叶变换的对偶转换
周期函数的傅里叶变换是离散的，离散函数的傅里叶变换是周期的，因此，周期和离散可以看作傅里叶对偶算子(傅里叶共轭)。一个周期函数fper(u)可以定义为任意函数f(u)周期延拓后的结果

fper(u)=∑+∞k=-∞f(u-kΔu)(3.229)

其中，Δu是周期，Δu＞0。一个离散函数fdis(u)可以定义为任意函数f(u)的均匀采样： 

fdis(u)=δu∑+∞k=-∞f(kδu)δ(u-kδu)(3.230)

其中，δu是采样间隔，δu＞0。
按照上述理解思路，不少常用运算对都是傅里叶对偶算子，如坐标乘法与微分、时移和相移、chirp乘法和chirp卷积，这些都是众所周知的例子。尺度算子以参数互为倒数而自对偶。分数阶傅里叶变换算子Fp是傅里叶变换算子F的推广形式。在连续形式上，分数阶傅里叶变换算子介于恒等算子F0=I和傅里叶算子F1=F之间。通过分数阶傅里叶变换可以对上述对偶算子进行转换。当变换阶次从0变到1时，对偶算子中的一个会逐渐转换到另外一个。
周期函数和有限延拓函数是密切相关的，周期函数可看作有限延拓函数的冗余表示； 相应地，有限延拓函数可看作周期函数的紧凑表示。因此，接下来的内容也是与有限延拓函数有关的。
从严格的数学意义上来说，一个函数f(u)和它的傅里叶变换F(μ)=∫+∞-∞f(u)·exp(-j2πμu)du不可能都是有限延拓的。但在实际中，当超出延拓范围之外的信号能量可以忽略时，至少可以假定它们是近似有限延拓的。假定f(u)的延拓范围是Δu，F(μ)的延拓范围是Δμ，并且都是关于原点对称的。那么根据采样理论，采样间隔δu=1/Δμ和δμ=1/Δu分别在u域和μ域都是足够的。换言之，信号在一个域的延拓范围(周期)就对应于另一个域的分辨率。
这意味着任何一个域上的ΔuΔμ个样本就足以完全表征该函数。我们用N来表示这个数，它也称为信号的时宽带宽积或者自由度。在一定意义上，超出Δu和Δμ外的信号能量是可以忽略的。离散傅里叶变换将N个时间样本映射为N个频率样本。两者的精确关系可以通过泊松公式给出［34］

Fper(m/Δu)=1Δμ∑N-1n=0fper(n/Δμ)exp(-j2πmn/N)(3.231)

上式可以看成f(u)的周期延拓和F(μ)的离散傅里叶变换结果。其中，0≤m≤N-1，m∈Z； Δu和Δμ任意取值，Fper(μ)=∑+∞k=-∞F(μ-kΔμ)，fper(u)=∑+∞k=-∞f(u-kΔu)。
3.3.1一般对偶算子及其分数阶版本
A的对偶算子用AD表示，并且满足


AD=F-1AF

A=FADF-1
(3.232)

AD对频域表示F(μ)所起作用与A对时域表示f(u)所起作用相同。本节所讨论的分数阶算子在分数阶域中所起作用与其对应的传统算子在时域中所起作用是一样的。准确地说，就是分数阶算子在p阶分数阶域对信号Fp(up)的作用与原算子在时域对信号f(u)的作用是一样的。用数学形式表示出来如下

Ap=F-pAFp(3.233)

上式是式(3.232)的推广。当p=1，A1=AD时，式(3.233)就变成了式(3.232)。注意： A和Ap是两个不同的算子，它们不是同一算子的不同表示，但是它们分别在时域和第p阶分数阶域的表示却是相同的。为了将第p阶分数阶域中的分数阶算子和算子A的p次幂(Ap)区分开来，我们将前者表示为Ap。可以发现，当p=0时，A0=A，且A0=I； 当p=1时，A1=AD，且A1=A。换言之，Ap是介于算子A和其对偶算子AD之间的，而Ap是介于恒等算子和算子A之间的。式(3.229)对AD也成立

ADp=F-pADFp(3.234)

既然从式(3.233)可以很容易地得出ADp=ApD，因此，将其简化表示为ADp。
接下来，将对一些经典对偶算子及其分数阶推广做个总结性的讨论，有关详细内容可参看文献［35］。
首先，基于坐标乘法算子U和微分算子D在时域中的作用，将这两种算子定义如下


Uf(u)=uf(u)(3.235)

Df(u)=j2π-1f′(u)(3.236)


令A=U，AD=D，可以很容易地看出这两个算子满足式(3.232)，即两者是对偶的。同样，D和-U也是一对对偶算子。可以发现，D在频域所起的作用就是坐标乘法，这与U在时域的作用是一样的。从这个性质可以得到

F(j2π)-1f′(u)(μ)=μF(μ)

这两个算子的分数阶形式定义如下


Up{f}(up)=upf(up)(3.237)

Dp{f}(up)=j2π-1f′(up)(3.238)


从这个定义，可以看出Up和Dp都满足式(3.233)和式(3.234)，且存在下面的几种特殊情况： U0=U，U1=D，U-1=-D，D0=D，D1=-U，D-1=U。
相移算子PH(η)和平移算子SH(ξ)定义为


PH(η)=expj2πηU(3.239)

SH(ξ)=expj2πξD(3.240)


这两个算子分别表示时域和频域上信号的平移，即


PH(η)f(u)=expj2πηuf(u)(3.241)

SH(ξ)f(u)=fu+ξ(3.242)


令A=PH(ξ)，AD=SH(ξ)或者A=PH(ξ)，AD=SH(-ξ)，可以利用式(3.232)将这两个算子联系在一起

［SH(ξ)］D=F-1SH(ξ)F=F-2PH(ξ)F2=PPH(ξ)P=PH(-ξ)

这两个算子的分数阶形式定义如下


PHp(η)=expj2πηUp(3.243)

SHp(ξ)=expj2πηDp(3.244)


可以看出，这两个分数阶算子在p阶分数域的作用与它们的一般形式在时域的作用是一样的，且都满足式(3.233)和式(3.234)。
利用坐标乘和微分算子给出尺度算子M(M)的定义如下

M(M)=exp-jπ(lnM)(UD+DU)(3.245)

其中，M＞0。它在时域里的作用表示为

M(M)f(u)=1/Mf(u/M)(3.246)

尺度算子是自对偶的，在时域的尺度变换对应于频域的逆尺度变换： 

F1/Mf(u/M)=MF(Mμ)

令A=M(M)且AD=M(1/M)时，该尺度算子满足式(3.228)。同理，可将分数阶尺度算子定义如下

Mp(M)=exp-jπ(lnM)(UpDp+DpUp)(3.247)

该算子在p阶分数域的作用与其一般形式在时域的作用是一样的，且都满足式(3.233)和式(3.234)。
Chirp乘算子Q(q)和chirp卷积算子R(r)定义如下


Q(q)=exp(-jπqU2)(3.248)

R(r)=exp(-jπrD2)(3.249)


它们在时域的作用可表示为


Q(q)f(u)=exp(-jπqu2)f(u)(3.250)

R(r)f(u)=exp(-jπ/4)1/rexp(jπu2/r)f(u)(3.251)


令A=Q(q)，AD=R(q)或者A=R(r)，AD=Q(r)，可以看出这两个算子也是傅里叶对偶算子，且满足式(3.232)。将它们用级数展开为


Q(q)=∑+∞k=0(-jπq)kk!U2k=∑+∞k=0(-jπq)kk!(-F-1DF)2k

=∑+∞k=0(-jπq)kk!F-1D2kF=F-1R(q)F

Q(q)=∑+∞k=0-jπqkk!FDF-12k=∑+∞k=0-jπqkk!FD2kF-1=FR(q)F-1


Chirp信号的乘法可以看成以线性调频信号作频率调制。既然Q(q)和R(r)是满足式(3.232)的对偶算子，那么R(r)在时域的作用就等价于Q(q)在频域的作用。因此，R(r)f(u)就等效于对频域函数e-jπrμ2F(μ)作逆傅里叶变换。
定义chirp乘算子Q(q)和chirp卷积算子R(r)的分数阶形式如下


Qp(q)=exp(-jπqU2p)(3.252)

Rp(r)=exp(-jπrD2p)(3.253)


3.3.2离散算子和周期算子以及它们的分数阶形式
根据相移算子和平移算子，可将离散算子DI(Δμ)和周期算子Pε(Δu)定义如下


DI(Δμ)=∑+∞k=-∞PHkΔμ=∑+∞k=-∞exp(j2πkΔμU)(3.254)

Pε(Δu)=∑+∞k=-∞SHkΔu=∑+∞k=-∞exp(j2πkΔuD)(3.255)


其中，Δu＞0和Δμ＞0分别对应于时域和频域的延拓周期。令 δu=1/Δμ，δμ=1/Δu分别表示时域和频域的采样间隔。
根据上面的定义以及相移算子和平移算子的对偶性，可以发现这两个算子也是傅里叶对偶算子


Pε(Δu)=∑+∞k=-∞SHkΔu=∑+∞k=-∞F-1PHkΔuF=F-1DI(Δu)F

Pε(Δu)=∑+∞k=-∞SH-kΔu=∑+∞k=-∞FPHkΔuF-1=FDI(Δu)F-1
(3.256)

通过式(3.254)和式(3.255)，可以得出∑+∞k=-∞SH-kΔu=∑+∞k=-∞SHkΔu，以及对偶关系PH(ξ)D=SH(ξ)，SH(ξ)D=PH(ξ)。
离散算子和周期算子在时域的作用可表示为


DI(Δμ)f(u)=∑+∞k=-∞expj2πku/δuf(u)=δu∑+∞k=-∞δu-kδufkδu(3.257)

Pε(Δu)f(u)=∑+∞k=-∞fu-kΔu(3.258)


其中，我们采用了泊松求和公式的另外一种形式［34］

∑+∞n=-∞δu-nδu=1δu∑+∞n=-∞expj2πnu/δu(3.259)

接下来，定义comb(u)函数如下

comb(u)=∑+∞k=-∞δu-k=∑+∞k=-∞δu+k(3.260)

这样，便能重写式(3.257)和式(3.258)为


DI(Δμ)f(u)=combuΔμf(u)(3.261)

Pε(Δu)f(u)=1ΔucombuΔuf(u)(3.262)


从上面的等式，可以看出DI(Δμ)是将一个时域函数乘以一个δ序列从而得到一个加权脉冲序列； 而Pε(Δu)是对一个时域函数作周期延拓。
既然Δu-1comb(u/Δu)的傅里叶变换就是comb(Δuμ)，那么式(3.261)和式(3.262)右边的傅里叶变换分别为


1ΔμcombμΔμF(μ)=∑+∞k=-∞Fμ-KΔμ(3.263)

combΔuμF(μ)=δμ∑+∞k=-∞δμ-kδμFkδμ(3.264)


从上面的式子可以看到这两个算子的作用在频域是互换的，因为离散算子对应于频域的周期延拓，周期算子对应于一个δ序列的相乘。
离散算子和周期算子的分数阶形式定义如下


DIp(Δμ)=∑+∞k=-∞PHpkΔμ =∑+∞k=-∞expj2πkΔμUp(3.265)

Pεp(Δu)=∑+∞k=-∞SHpkΔu=∑+∞k=-∞expj2πkΔuDp(3.266)


显然上述算子满足式(3.233)


DIp(Δμ)=F-pDI(Δμ)Fp(3.267)

Pεp(Δu)=F-pPε(Δu)Fp(3.268)


基于上面的等式和式(3.256)，可以进一步得到

Pεp(Δu)=F-1DIp(Δu)F=FDIp(Δu)F-1(3.269)

既然PH(η)k=PH(kη)和SH(ξ)k=SH(kξ)，那么将式(3.265)和式(3.266)右边展开可以得到



DIp(Δμ)=∑+∞k=-∞expjπkΔμ2sinαcosαPHkΔμcosαSHkΔμsinα

=∑+∞k=-∞expjπkΔμ2sinαcosαPHΔμcosαkSHΔμsinαk(3.270)

Pεp(Δu)=∑+∞k=-∞expjπkΔu2sinαcosαSHkΔucosαPH-kΔusinα

=∑+∞k=-∞expjπkΔu2sinαcosαSHΔucosαkPH-Δusinαk(3.271)



因为


PHp(η)=exp(jπη2sinαcosα)PH(ηcosα)SH(ηsinα)

=exp(jπη2sinαcosα)exp(j2πηcosαU)SH(ηsinα)

=exp(-jπcotαU2)exp［jπcotα(U+ηsinαI)2］SH(ηsinα)

=exp(-jπcotαU2)SH(ηsinα)exp(jπcotαU2)


所以

PHp(η)=QcotαSHηsinαQ-cotα(3.272)

类似地，我们可以得到

SHp(ξ)=Q-tanαSHξcosαQtanα(3.273)

那么将式(3.272)代入式(3.265)，可以得到


DIp(Δμ)=∑+∞k=-∞Q(cotα)SH(kΔμsinα)Q(-cotα)

= Q(cotα)∑+∞k=-∞SH(kΔμsinα)Q(-cotα)

=Q(cotα)Pε(Δμsinα)Q(-cotα)

(3.274)

同样，将式(3.273)代入式(3.266)，可以得到

Pεp(Δu)=Q-tanαPεΔucosαQtanα(3.275)

既然(AD)p=(Ap)D，基于Q和R，PH和SH的对偶关系以及式(3.256)，可以得到


DIp(Δμ)=QcotαF-1DIΔμsinαFQ-cotα

DIp(Δμ)=F-1RcotαDIΔμsinαR-cotαF

FDIp(Δμ)F-1=RcotαDIΔμsinαR-cotα

Pεp(Δu)=RcotαDIΔusinαR-cotα

(3.276)

类似地，有


Pεp(Δu)=Q-tanαF-1DIΔucosαFQtanα

Pεp(Δu)=F-1R-tanαDIΔucosαRtanαF

FPεp(Δu)F-1=R-tanαDIΔucosαRtanα

DIp(Δμ)=R-tanαDIΔμcosαRtanα
(3.277)

式(3.274)~式(3.277)通过一般的周期和离散算子以及普通chirp信号来表示分数阶周期/离散算子。下面给出分数阶周期算子和分数阶离散算子对函数f(u)的作用表达式


Pεp(Δu)［f(u)］=∑+∞k=-∞exp-jπk2Δu2sinαcosα·

expj2πkΔusinαufu-kΔucosα(3.278)

DIp(Δμ)f(u)=∑+∞k=-∞expjπk2Δμ2sinαcosα·

exp-j2πkΔμcosαufu-kΔμsinα(3.279)


表3.1给出了本节所列出的分数阶算子及其等效表达式。


表3.1分数阶算子及其等效表达式



分数阶算子符号等 效 表 达

分数阶坐标乘UpcosαU+sinαD
分数阶微分Dp-sinαU+cosαD
续表



分数阶算子符号等 效 表 达


分数阶相移PHp(η)exp(-jπη2sinαcosα)PH(ηcosα)SH(ηsinα),

exp(-jπη2sinαcosα)PH(ηcosα)SH(ηsinα)

分数阶平移SHp(ξ)exp(jπξ2sinαcosα)SH(ξcosα)PH(-ξsinα),

exp(-jπξ2sinαcosα)PH(-ξsinα)SH(ξcosα)

分数阶尺度变换Mp(M)F-pM(M)Fp,M-1-pM(1/M)F1+p
分数阶chirp乘Qp(q)R(-tanα)Q(qcos2α)R(tanα),

Q(cotα)R(qsin2α)Q(-cotα)

分数阶chirp卷积Rp(r)R(cotα)Q(rsin2α)R(-cotα)

Q(-tanα)R(rcos2α)Q(tanα)

分数阶离散DIp(Δμ)R(-tanα)DI(Δμcosα)R(tanα)

Q(cotα)PE(Δμsinα)Q(-cotα)

分数阶周期PEp(Δu)R(cotα)DI(Δusinα)R(-cotα)

Q(-tanα)PE(Δucosα)Q(tanα)

参考文献

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

［2］Zayed A I.A convolution and product theorem for the fractional Fourier transform［J］.IEEE Signal Processing Letters,1998,5(4): 101103.

［3］Ozaktas H M,Barshan B.Convolution,filtering,and multiplexing in fractional Fourier domains and their relation to chirp and wavelet transforms［J］.J.Opt.Soc.Am.A,11(2): 547559,1993.

［4］Olcay A.Fractional convolution and correlation via operator methods and an application to detection of linear FM signals［J］.IEEE Trans.Signal Processing,2001,49(5): 979993.

［5］Akay O.Unitary and Hermitian fractional operators and their extension: Fractional Mellin transform,joint fractional representations and fractional correlations［D］.Univ.Rhode Island,Kingston,2000.

［6］Wang M，Chan A K,Chui C K.Linear frequencymodulated signal detection using RadonAmbiguity transform［J］.IEEE Trans.Signal Processing,1998,46(3): 571586.

［7］Akay O,BoudreauxBartels G F.Unitary and hermitian fractional operators and their relation to the fractional Fourier transform［J］.IEEE Signal Processing Letters,1998,5(12): 312314.

［8］Pei S C,Yeh M H.Discrete fractional Hilbert transform［J］.IEEE Trans.Circuits and SystemsII,2000,47(11): 13071311.

［9］Tseng C C,Pei S C.Design and application of discretetime fractional Hilbert transformer［J］.IEEE Trans.Circuits and SystemsII,2000,47(12): 15291533.

［10］Pei S C,Ding J J.Fractional cosine,sine,and Hartley transforms［J］.IEEE Trans.Signal Processing,2002,50(7): 16611680.

［11］Gabor D.Theory of communications［J］.Inst.EE,1946,93(26): 429457.

［12］Lohmann A W,Mendlovic D,Zalevsky Z.Fractional Hilbert transform［J］.Opt.Lett.,1996,21: 281283.

［13］Zayed A I.Hilbert transform associated with the fractional Fourier transform［J］.IEEE Signal Processing Lett.,1998,5 (8): 206208.

［14］Pei S C,Ding J J.Saving the bandwidth in the fractional domain by generalized Hilbert transform pair relations［J］.IEEE Trans.Circuits and SystemsⅡ,2003,4: 8992.

［15］Zayed A I.On the relationship between the Fourier and fractional Fourier transforms［J］.IEEE Signal Processing Lett.,1996,3(12): 310311.

［16］Havlicek J P,Havlicek J W,et al.Skewed 2D Hilbert transforms and computed AMFM models［J］.Proc.IEEE,1998,59: 602606.

［17］Tao R,Li X M,Wang Y.Generalization of the fractional Hilbert transform［J］.IEEE Signal Processing Letters,2008,15: 365368.

［18］Wang Z.Fast algorithm for discrete W transform and for the discrete Fourier transform［J］.IEEE Trans.Acoust.,Speech,Signal Processing,1984,32: 803816.

［19］Pei S C,Yeh M H.Improved discrete fractional Fourier transform［J］.Opt.Lett.,1997,22: 10471049.

［20］Pei S C,Yeh M H,Tseng C C.Discrete fractional Fourier transform based on orthogonal projection［J］.IEEE Trans.Signal Processing,1999,47: 13351348.

［21］Dickinson B W,Steiglitz K.Eigenvectors and functions of the discrete Fourier transform［J］.IEEE Trans.Acoust.,Speech,Signal Processing,1982,30: 2531.

［22］SooChang Pei.Discrete fractional Hartley and Fourier Transforms［J］.IEEE Trans.circuit and systems,1998,45(6): 665675.

［23］Akan A,Shakhmurov V.A fractional Gabor transform［J］.ICASSP01,2001,6: 35293532.

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

［25］Alieva T,Bastiaans M J.On fractional Fourier transform moments［J］.IEEE Signal Processing Letters,2000,7(11): 320323.

［26］Qiao X,Shan T,Tao R,et al.Separation of human microDoppler signals based on shorttime fractional Fourier transform［J］.IEEE Sensors J.,2019,19(24): 1220512216.

［27］张贤达.现代信号处理［M］.北京： 清华大学出版社,1995.

［28］Chen X,Guan J,Bao Z,et al.Detection and extraction of target with micromotion in spiky sea clutter via shorttime fractional Fourier transform［J］.IEEE Trans.Geosci.Remote Sens.,2014,52(2): 10021018.

［29］Pang C,Han Y,Hou H,et al.Microdoppler signal timefrequency algorithm based on STFRFT［J］.Sensors,2016,16(10): 1559.

［30］Tao R,Li Y L,Wang Y.Shorttime fractional Fourier transform and its applications［J］.IEEE Trans.Signal Process.,2010,58(5): 2568 C2580.

［31］Cohen L. Timefrequency analysis［D］.PrenticeHall,Englewood Cliffs,NJ,1995.

［32］Baraniuk R G.Beyond timefrequency analysis: energy densities in one and many dimensions［J］.IEEE Trans.Signal Process,1998,46: 23052314.

［33］Sayeed A M,Jones D L.Integral transforms covariant to unitary operators and their implications for joint signal representations［J］.IEEE Trans.Signal Processing,1996,44: 13651376.

［34］Papoulis A.Signal analysis［M］.New York: McGrawHill,1977.

［35］Sumbl U,Ozaktas H M.Fractional free space,fractional lenses,and fractional imaging systems［J］.J.Opt.Soc.Amer.A,2003,20(11): 20332040.

［36］Huang Y.The fractional wave packet transform［J］.Multidimensional system and Signal processing,1998,9: 399402.