分数阶傅立叶变换的最优阶数论文

1、现代数字信号处理学号： 140808040219学生所在学院：测试与光电工程学院 学生姓名：任 课 教 师 ：李 志 农 教师所在学院：测试与光电工程学院2015年 1 月分数阶傅立叶变换的最优阶数摘 要：传统傅立叶变换描述了信号时域或频域的特性，而不是描述信号时频特性，于是人 们提出了一系列新的时频分析理论和方法来处理非平稳信号，分数阶傅立叶变换为其中的一种。本文主要介绍了它的定义、性质，还有它的离散算法，介绍了求最优阶数的方法，主要 是峰值搜索算法。最后进行仿真验证，利用MATLAB对一个已知的chirp信号求解最优阶数。关键字：傅立叶变换；分数阶傅立叶变换；峰值搜索算法；MATLAB ;

2、最优阶数1引言传统的傅立叶变换在所有的信号处理工具中是应用最广泛、 研究最成熟的数 学工具，作为一种线性算子，传统傅立叶变换可视为在时频平面上， 信号从时间 轴逆时针旋转二到频率轴，而FRFT作为FT的广义形式可理解为对信号旋转2任意角度的线性算子，从而可以得到信号的任意阶次或者任意分数阶傅立叶域上 的FRFT表示，并且在保留了传统的 FT所有性质和优点的基础之上又增添了 新的优势。2 FRFT的定义及其性质1.1 FRFT的定义如图1所示如果把信号的分数阶傅立叶变换看作是从时间 -频率平面旋转的 话，那么傅立叶变换就相当于在时频平面中逆时针旋转了7. 2角度，从时间域变换到频率域。令：二px

3、2，p是一个分数，那么就可以在时频平面内以任意角 度的旋转定义线性算子 R = Rp二2，记作F p，我们就可以把傅立叶变换推广到任 意角度即分数阶傅立叶变换。图1.(t,w)平面旋转角度变成(u,v)平面分数阶傅立叶变换的定义为：Xp(u) =Fpx(u) =Kp(u,t)x(t)dt(2-1)其中是核函数Kp(u,t),1 - /COt?/ + Ua = 2 rmexp /(cot a - Uiesc a), a r nn 2;r26t hexp(-y),(2-2)这里a =兰L,0 c p c2 , p=1或-1时退化成为常规的傅立叶变换和逆变换。2分数阶傅立叶变换和经典傅立叶变换具有以

4、下的关系：1. 分数阶傅立叶变换是线性算子2周期性F0x(t)HF4x(t)Hx(t)(2-3)F1x(t)HF2x(t)HX()(2-4)1.2 分数傅立叶变换的性质(1) 线性分数阶傅立叶变换为线性变换，满足叠加原理：若FPf(t)和FPg(t)分别是原函数f (t)和g(t)的p阶分数傅立叶变换，则有FP(Gf(t) og(t) FPf(t) qFPg(t)(2-5)(2)旋转相加行FP Fq= FP q(2-6)(3)逆(FP)F 扌(2-7)(4)酋性(FP) (FP)H(2-8)(5)交换性FP1F JFP2FP1(2-9)(6)结合性(FPiFP2)FpFf1(FP2FP3)(7

5、)周期性p4k(u)=Fp4kf(x) = Fpf(x) = fp(u)(8)特征函数F p- | =exp(-jpl 二 /2) - |卷积、相乘、相关(9) 函数f (t)、g(t)在阶次p分数阶傅立叶域的卷积记作Fpf*Fpg = fp(tp)*gp(tp)函数f(t)、 g(t)在阶次p分数阶傅立叶域的乘积记作Fpf FpgH fp(tp)gp(tp)函数f(t)、 g(t)在阶次p分数阶傅立叶域的相关定义为CORRpf(t),g(t)二 CPNVpf(t),g (-t)(10)时移特性pj sin ：.cs：.-ju sinF px(t - ) =Xp(u - cos： )e 2(1

6、1)频移特性1 2ppt_j( V sin :.cos：. JJ sin ：)Fpx(t)e=e 2Xp(uvsin：)(12)尺度特性pr jcot 吟 cot：(1Fpx(ct)24 e 2c _jcotacos2、-2 )sin 屮cosXp(u )csin a其中 护=arctan(c2 tan ： ) = p二 /2(13)微分特性F px(t)二 XpWjcos：juXp(u)sin :-(14)积分特性2 2lju tan： uj-z tan、x( )d 二 sec： e2 Xp(z)e 2 dz aa p(2-10)(2-11)(2-12)(2-13)(2-14)(2-15)(

7、2-16)(2-17)(2-18)(2-19)(2-20)3分数阶傅立叶变换的离散的离散算法分数阶傅立叶变换的出现引起了各个领域研究人员的重视，在工程上也有十分广阔的应用前景。在数字信号处理的应用中，必须采用离散形式的分数阶傅立 叶变换(DFRFT)，这使得离散分数阶傅立叶变换及其快速算法的研究显得十分 重要。目前，DFRFT的离散化算法主要有四种：2g31. 利用： W来计算离散FRFT的核矩阵，从而利用 FFT来i=0计算DFRFT其中W是离散傅立叶变换核矩阵。这种方法实际上缺乏理论基础， 而且其离散FRFT矩阵不满足连续FRFT的旋转相加性，因此不能用相同的方 法计算逆FRFT。实际计算

8、所产生的误差比较大，与连续 FRFT没有相似的输 出结果。其计算复杂度与传统傅立叶变换相同，为O(Nlog2N)。2. 分解方法根据FRFT的表达式，将FRFT分解为信号的卷积形式，从而利用FFT来 计算FRFT。这种方法思想比较直观，计算出的记过与连续 FRFT的输出比较接 近。但它要经过一次2倍内插和2倍抽取，而且还要进行坐标的无量纲化， 实 现起来较为烦琐。其计算复杂度为 O(N log2 N)。3. 利用矩阵的特征值和特征向量来计算 DFRFT这种方法保持了连续FRFT的特征值-特征函数的关系，克服了第一种方法 中特征值和特征向量不匹配的缺点。采用了两种正交映射的办法DFT的Hermi

9、te特征向量，由此开发出两种快速方法，即OPA方法和GSA方法，这两种方法都有和连续 FRFT相近的输出结果，可逆性好。其计算复杂度均为 O(N2)。4. 直接对FRFT进行离散化来计算 DFRFT。这种方法采用直接将连续 FRFT离散化的方法来获得离散 FRFT的核矩 阵，避开了烦琐的特征值和特征向量匹配问题以及矩阵的正交归一化运算， 与连 续FRFT有相似的输出结果，该算法的计算复杂度为 O(N2)。在目前的研究中，采用的最多的是分解方法和矩阵特征值和特征向量两种方 法，下面的内容重点介绍这两种方法。3.1分解方法所谓分解方法是指根据 FRFT的表达式，将 FRFT分解为信号的卷积形 式，

10、从而利用FFT来计算FRFT。这种算法由H.M.Ozaktas等人提出，其计算 速度几乎与FFT相当，被公认为目前计算速度最快的一种 FRFT数字计算方 法，非常适合于对信号进行实时的 FRFT计算。但这种快速算法的运算机理决 定了在进行FRFT数值计算之前必须对原始信号进行量纲归一化处理。3.1.1量纲归一化原理如果一个信号在时间轴或频率轴的一个子集取非零值，并且取非零值的条件限定在有限区间内，则称该信号在时间轴或频率轴上是紧凑的。从理论上讲，任何信号和它的傅立叶变换不可能是同时紧凑的。 然而我们实际中需要处理的信号 往往是时限的和带限的。信号的时宽带宽积可以用来确定信号的采样频率和采样 点

11、数，用于唯一地从离散化的信号中恢复原始信号。假定原始连续信号在时间轴和频率轴上都是紧凑的，其时域表示限定在区间 -：t/2,：t/2，而其频域表示限定在区间f/2,f/2， ：t和 f分别表示 信号的时宽和带宽。信号的时宽带宽积为 n，根据不确定性定理，n的值应当大于1。由于时域和频域具有不同的量纲，为了FRFT计算处理方便，需要将时域和频域分别转换成无量纲的域。具体方法是引入一个具有时间量纲的尺度因子S，并定义新的尺度化坐标x=t/S,v = f/ S新的坐标系(x,v )实现了无量纲 化。信号在 新坐标系中被限定 在区间 t/2St/2S和f S/2 ：f S/2内。 为使两个区间的长度相

12、 等，选择S ：tr ：f，贝U两个区间长度都等于无量纲量W f，即两个区间归一化为J：x/2, .x/2。归一化以后信号的 Wigner-Ville分布限定在以原 点为中心，直径 厶x的圆内，如图3-2所示。最后根据采样定理对归一化后的 信号进行采样，采样间隔为1 / x，采样点数为N -：x2.图3.1归一化后的时频支撑区域3.1.2分解算法可以把(2-1)式改写为11 一 j COtaU2 COta 说t2 COh(3-1)Xp(u)expj x(t)expjexp-jtucsc dt只2a2式可以具体分解为以下几个步骤:(1)用chirp 信号调制信号x(t):g(t)二 x(t)ex

13、p - j二t2 tan(： /2)(3-2)(2)调制信号与另一个chirp信号卷积：: 2g：.(u)二 A：.：exp-j二(u-t) csc g(t)dt(3-3)(3)用chirp信号调制卷积后的信号：Xp(u)= g,(u)-j 二t2ta n(: / 2)(3-4)要实现FRFT的数值计算必须对以上每个分解步骤都进行离散化处理，具体的实现过程如下：(1)对x(t)与chirp信号的乘积进行采样假定分数阶次P，-1,1，chirp信号的调频率tan(:/2)乞1， x(t)经chirp信号调制后所得的信号时宽带宽积可以是原信号的时宽带宽积的两倍，所以要求x(t)的采样间隔为1/2

14、：x，如果原来的采样间隔是1/ ：x，可通过插值的方法获 得样本值，然后再chirp信号的离散采样相乘，以得到g(t)的采样值。(2) 实现g(t)与chirp信号的卷积由于g(t)是带限信号，所以chirp信号也可取其带限形式。所以有：: 2 :(3-5)g：.(u) =A：, ：exp-j：(u -t) esc： g(t)dt = A：.： h(u -t)g(t)dt心2其中 h(u) H(v)exp(j二2vu)du,expj- (esc： )u2)而H (v)则是chirp信号傅立叶变换1j兀H(V)1csc 5?旳j V2csc :(3-6)于是，式的离散形式为(3-7)二 A* h

15、 口_N这一离散卷积可利用 FFT来计算。(3) 计算分数阶傅立叶变换Xp(u)的采样值Xp(m/2：x)由于在第一步操作时对信号作了2倍内插操作，所以在最后的结果需要对Xp(m/2x)再进行2倍抽取，以得到离散采样Xp(mAx)。总之上述方法从连续信号x(t)的N个离散采样x(m/cx)开始，最后得到由 Xp(u)的N个离散样本Xp(m/.x)值得注意的是上述方法只适用于 -1乞p1的 情况，对位于该区间外的情况，可以利用分数阶傅立叶变换的周期性将阶次变到 -1乞p叨后再进行计算。3.1.3特征值和特征向量方法分解方法虽然计算复杂度较低，但不严格遵守分数阶傅立叶变换的旋转特 性，因此不能从变

16、换后的信号通过其逆变换精确地恢复原始信号。为了克服这种方法的缺点和不足，Pei Soo-Cha ng等人提出一种新的离散化方法，该方法具有 与连续情况相似的变换性质和结果，并可以通过其逆变换恢复原始信号。这种方法就是特征值和特征向量方法，它从连续傅立叶变换的特征函数为 Hermite函数出发，通过对Hermite函数的离散化近似和正交投影，得到一组与Hermite函数形状相似的 DFT矩阵的正交化离散 Hermite特征向量，然后，按 照连续分数阶傅立叶变换的核函数谱分解表达式，构造离散分数阶傅立叶变换矩 阵。3.1.3.1 DFT矩阵特征值和特征向量傅立叶变换的特征函数 n(u)为Hermi

17、te-Gaussian函数，其表达式为：*n(u)二代Hn(2u)e(3-8)2“4其中An =，n = 0,1，n; Hn(u)为n阶Hermite多项式。DFT矩阵的特征值及重复度如下表：表3-1DFT矩阵特征值及重复度DFT矩阵F的特征值是e对应的特征向量全体组成一个特征子空间，记为 复度决定了子空间的秩。矩阵S-j(i/2)k，共有1、j -1、j四个值，每个特征值Eo,E!,E2,E3，每个特征值的重DFT矩阵F的特征向量，S的表可用于计算N1的重負度的重复度-1的畫整度j的重如度4 mm+ 1mIB J4m+lni+1mmm4m+2m+ 1mm+lmIrti 十 3m+im+lm+

18、1m达式为:210 1112cos 国1 00a1a2cos2a+0a00 2cos( N -1) jS 二(3-9)SF = FS。因此矩阵S的特征向量可以证明矩阵S和F满足乘法交换性，即 也是矩阵F的特征向量，但它们对应不同的特征值。3.1.3.2 DFT矩阵Hermite特征向量的计算由于DFT矩阵F的特征向量不唯一，我们希望从中找到与 Hermite函数相 似的特征向量，这样的向量称 DFT矩阵Hermite特征向量，为了求 Hermite特 征向量，下面给出一系列的重要定理。定理1对DFT矩阵的Hermite特征向量而言，它对应的连续函数的扩展 方差应当为N/2二)Ts是信号的采样间

19、隔，连续 Hermite函数采样后得到：勺(k) = I h kjexp(k2兀/N)(3-10)J2nn! JN /2Ts&N/2兀丿上式也可以看做是方差为1的Hermite-Gaussian函数，以T - 2 / N为采样间隔进行采样得到的序列。定理2 若序列n(k)是由单位方差Hermite-Gaussian函数经过T 2二/ N 采样到的，那么，可以证明下列近似等式成立：当N为偶数时(N/2) J(3-11)N为奇数时Gj)n n(k):(N/2) A送%(m)ej2加Nm _(N/2)/2(3-12)定理3将序列n(k)按照以下方式平移得到(k)N为偶数时(3-13)ng。*乞 N/

20、2-1n(kN),N/2Ek EN -1N为奇数时；(k) =n(k),Ok z(N -1)/2n(k N),(N1)/2 Ek 乞 N 1(3-14)则；(k)的离散傅立叶变换近似为(-j)n(k)，即当N足够大时N/2 (一 dfN:戶严k/N(3-15)DFT矩阵特征(3-16)从定理2和定理3可以看出，Hermite函数的采样序列近似为 向量。对Hermite函数的采样序列作归一化，以un记为： u _(0),%(1),(N 1)Tn |可(0),可，可(N -1)T|通过S矩阵可以得到 DFT矩阵F的一组实正交特征向量。因此可以将这些特 征向量作为DFT特征子空间的基向量。然后计算向

21、量 un在DFT特征子空间 的投影，从而得到 DFT矩阵的Hermite特征向量。计算方法如下：山=嘉汕n,Vk；Vk(3-17)(n -k)mod4 -03.1.3.3 DFRFT 核矩阵和二维 DFRFT4最优变换阶次的选择本文选择最优变换阶次的方法主要是峰值搜索法，目前主要搜索算法有两种二维搜索算法、拟牛顿搜索算法。尽管拟牛顿算法比二维搜索的计算量有所减但 是这两种算法的计算复杂度还是不能满足工程实际的实时处理要求，因此我们一维曲线拟合来代替二维搜索的峰值搜索算法，计算量大大简化。下面主要介绍几 种算法。4.1二维搜索算法传统的二维搜索法若以变换阶数P为变量进行扫描搜索的过程中,当参数估

22、计的估计精度要求较高时，为了满足精度要求，变换阶数P必须选择小的搜索步 长,这样就会成倍地增加计算的复杂度，如变换阶数在区间0,2上,若以步长0.0001进行搜索,则需要进行20000次的FRFT,计算量很大,通常采用步进搜索 算法，即先以大步长再以小步长搜索最大值点，详细步骤为：先以0.01为步长进行 阶数从0-2搜索最大值点Pmax,再以步长0.0001进行阶数从Pmax-0.01至+ 0.01搜 索最大值点。这样一共需要进行400次FRFT运算,计算量有所减小。4.2拟牛顿搜索算法对于拟牛顿搜索算法，首先介绍一下拟牛顿迭代法的基本思想。采用牛顿法 解非线性方程f(x)就是把非线性线性化的

23、一种近似方法。把f(X)= 0在X0点展为 泰勒级数得，f (n) (x )f (x) = f (Xo) (x-Xo) f (x) (x-Xo)n0(4-1)n!取方程的线性部分作为非线性方程f (X)= 0的近似方程，则有f(x0)(x-X0)f(x) =0(4-2)假设f (X)= 0则方程的解为X1- f(X0)/ f(x)(4-3)再把f (X)在X1点展开为泰勒级数，同样也取其线性部分作为非线性方程f (x) = 0的近似方程，假设，则有(4-4)(4-5)(4-6)X2 二为-f(X1)/ f(X1)Pn4f半L/n-Un 十-1 1 n n 丿nXp(u)2Xn 1 二 Xn -

24、 f(Xn)/ f(Xn) 0|Xp(u)|2cuPF这样即可得到牛顿法的一个迭代序列：U=tln其中是第Pn和Un是第n次搜索的结果，n为第n次搜索的步长系数，H.为函2数Xp(u)在(4,山)点的尺度矩阵可以通过迭代的方法求得。 牛顿迭代法中每次 迭代的主要运算为一次一维搜索和函数 Xp(u)2的一阶偏导数的计算，这一拟牛顿算法所需要的计算为分数阶傅里叶域的一次扫描搜索和一次迭代搜索，因为迭代搜索所需的搜索和函数的一阶偏导数的计算量远远小于扫描搜索 ，若扫描点数 为m,信号的样本长度为N,则拟牛顿算法的计算复杂度为 O(mNlgN),与其它基 于二维时频分布的算法(运算量一般为o(N3)相

25、比,拟牛顿算法计算量较小。5仿真验证本文中利用二维搜索法求最优阶次，该信号是一个chirp信号以相加的形式调制简单的高斯信号的混合信号，对该信号进行求解它的最优阶数。编写的分数 阶傅立叶变换函数见附录 A所示，主函数见附录 B所示。根据运行的结果为 P=0.79为所求的最优阶数。在该主程序中还用到了阶数 P=0.48,P=1.35,和P=1.82 和所求的最优阶数进行对比，很明显P=0.79是它的幅值最大。初始信号见下图4.1所示，最优阶图见图4.2所示。图4.1初始信号图图4.2最优阶图冋48分数阶Fourier换團尸0 79分数阶Fcwier变换團151050. JAFJ-1 .比分数阶F

26、ounw变换劉.亚分數阶Fouwr变换舉1510亠- 50列:0505图4.3 P=0.48, P=1.35, P=1.82和P=0.79分数阶傅立叶变换图参考文献1 平先军，陶然，周思永等.一种新的分数阶傅立叶变换快速算法J.电子学 报,2001,29(3):406-4082 张立浩 . 分数阶傅立叶变换及其在通信中的应用 D. 硕士学位论文 ，燕山大学 ，2012.3 张兆祥 . 基于分数阶傅立叶变换的数字水印算法研究 D. 硕士学位论文 ，华北电力大 学，2007.4 陈恩庆 . 基于分数阶 fourier 变换的 ofdm 系统关键技术研究 D . 硕士学位论文 ，北京理工 大学 ，

27、2007.5 王仁明 ， 张春生 ， 贾玉坤 .基于分数阶傅里叶域均衡的无线测控系统研究 J. 航天电子对 抗， 2012， 28(3):47-50.程乃平，席有猷，郝建华.基于分数阶傅里叶变换的LFM干扰抑制算法J.装备学院学报2014， 25(1):74-77.7 陈恩庆 ， 陶然.一种基于分数阶傅里叶变换的 OFDM 系统及其均衡算法 J. 电子学报 ， 2007， 35(3):410-4148 殷敬伟 ，惠俊英，蔡平等 . 基于分数阶 Fourier 变换的水声信道参数估计 J. 系统工程与电 子技术 ， 2007， 29(10):1625-1627.附录 Afunction Faf

28、= frft(f, a)% The fast Fractional Fourier Transform% input: f = samples of the signal% a = fractional power% output: Faf = fast Fractional Fourier transform error(nargchk(2, 2, nargin);f = f(:);N = length(f);shft = rem(0:N-1)+fix(N/2),N)+1;sN = sqrt(N);a = mod(a,4);% do special casesif (a=0), Faf =

29、f; return; end;if (a=2), Faf = flipud(f); return; end;if (a=1), Faf(shft,1) = fft(f(shft)/sN; return; end if (a=3), Faf(shft,1) = ifft(f(shft)*sN; return; end % reduce to interval 0.5 a 2.0), a = a-2; f = flipud(f); endif (a1.5), a = a-1; f(shft,1) = fft(f(shft)/sN; end if (a0.5), a = a+1; f(shft,1)

30、 = ifft(f(shft)*sN; end % the general case for 0.5 a 1.5 alpha = a*pi/2;tana2 = tan(alpha/2);sina = sin(alpha);f = zeros(N-1,1) ; interp(f) ; zeros(N-1,1);% chirp premultiplicationchrp = exp(-i*pi/N*ta na2/4*(-2*N+2:2*N-2)2);f = chrp.*f;% chirp convolutionc = pi/N/sina/4;Faf = fcon v(exp(i*c*(-(4*N-4):4*N-4)42),f);Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);% chirp post multiplicationFaf = chrp.*Faf;% normalizing constantFaf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);function xint=interp(x)% sinc interpolationN = length(x);