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

《分数阶傅立叶变换的最优阶数论文》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

现代数字信号处理学号：140808040219学生所在学院：测试与光电工程学院学生姓名：任课教师：李志农教师所在学院：测试与光电工程学院2015年1月 分数阶傅立叶变换的最优阶数摘要：传统傅立叶变换描述了信号时域或频域的特性，而不是描述信号时频特性，于是人们提出了一系列新的时频分析理论和方法来处理非平稳信号，分数阶傅立叶变换为其中的一种。本文主要介绍了它的定义、性质，还有它的离散算法，介绍了求最优阶数的方法，主要是峰值搜索算法。最后进行仿真验证，利用MATLAB对一个已知的chirp信号求解最优阶数。关键字：傅立叶变换；分数阶傅立叶变换；峰值搜索算法；MATLAB；最优阶数1引言传统的傅立叶变换在所有的信号处理工具中是应用最广泛、研究最成熟的数学工具，作为一种线性算子，传统傅立叶变换可视为在时频平面上，信号从时间轴逆时针旋转到频率轴，而FRFT作为FT的广义形式可理解为对信号旋转任意角度的线性算子，从而可以得到信号的任意阶次或者任意分数阶傅立叶域上的FRFT表示，并且在保留了传统的FT所有性质和优点的基础之上又增添了新的优势。2FRFT的定义及其性质1.1FRFT的定义如图1.所示如果把信号的分数阶傅立叶变换看作是从时间-频率平面旋转的话，那么傅立叶变换就相当于在时频平面中逆时针旋转了角度，从时间域变换到频率域。令，是一个分数，那么就可以在时频平面内以任意角度的旋转定义线性算子，记作，我们就可以把傅立叶变换推广到任意角度即分数阶傅立叶变换。图1.平面旋转角度变成平面 分数阶傅立叶变换的定义为：(2-1)其中是核函数，(2-2)这里，p=1或-1时退化成为常规的傅立叶变换和逆变换。分数阶傅立叶变换和经典傅立叶变换具有以下的关系：1.分数阶傅立叶变换是线性算子2.周期性(2-3)(2-4)1.1分数傅立叶变换的性质（1）线性分数阶傅立叶变换为线性变换，满足叠加原理：若和分别是原函数和的阶分数傅立叶变换，则有(2-5)（2）旋转相加行(2-6)（3）逆(2-7)（4）酋性(2-8)（5）交换性(2-9) （1）结合性(2-10)（2）周期性(2-11)（3）特征函数(2-12)（4）卷积、相乘、相关①函数、在阶次p分数阶傅立叶域的卷积记作(2-13)②函数、在阶次p分数阶傅立叶域的乘积记作(2-14)③函数、在阶次p分数阶傅立叶域的相关定义为(2-15)（5）时移特性(2-16)（6）频移特性(2-17)（7）尺度特性(2-18)其中（8）微分特性(2-19)（9）积分特性(2-20)3分数阶傅立叶变换的离散的离散算法 分数阶傅立叶变换的出现引起了各个领域研究人员的重视，在工程上也有十分广阔的应用前景。在数字信号处理的应用中，必须采用离散形式的分数阶傅立叶变换（DFRFT），这使得离散分数阶傅立叶变换及其快速算法的研究显得十分重要。目前，DFRFT的离散化算法主要有四种：1.利用来计算离散FRFT的核矩阵，从而利用FFT来计算DFRFT其中W是离散傅立叶变换核矩阵。这种方法实际上缺乏理论基础，而且其离散FRFT矩阵不满足连续FRFT的旋转相加性，因此不能用相同的方法计算逆FRFT。实际计算所产生的误差比较大，与连续FRFT没有相似的输出结果。其计算复杂度与传统傅立叶变换相同，为。2.分解方法根据FRFT的表达式，将FRFT分解为信号的卷积形式，从而利用FFT来计算FRFT。这种方法思想比较直观，计算出的记过与连续FRFT的输出比较接近。但它要经过一次2倍内插和2倍抽取，而且还要进行坐标的无量纲化，实现起来较为烦琐。其计算复杂度为。3.利用矩阵的特征值和特征向量来计算DFRFT这种方法保持了连续FRFT的特征值-特征函数的关系，克服了第一种方法中特征值和特征向量不匹配的缺点。采用了两种正交映射的办法DFT的Hermite特征向量，由此开发出两种快速方法，即OPA方法和GSA方法，这两种方法都有和连续FRFT相近的输出结果，可逆性好。其计算复杂度均为。4.直接对FRFT进行离散化来计算DFRFT。这种方法采用直接将连续FRFT离散化的方法来获得离散FRFT的核矩阵，避开了烦琐的特征值和特征向量匹配问题以及矩阵的正交归一化运算，与连续FRFT有相似的输出结果，该算法的计算复杂度为。在目前的研究中，采用的最多的是分解方法和矩阵特征值和特征向量两种方法，下面的内容重点介绍这两种方法。3.1分解方法所谓分解方法是指根据FRFT的表达式，将FRFT分解为信号的卷积形式，从而利用FFT来计算FRFT。这种算法由H.M.Ozaktas等人提出，其计算速度几乎与FFT相当，被公认为目前计算速度最快的一种FRFT数字计算方法，非常适合于对信号进行实时的FRFT计算。但这种快速算法的运算机理决定了在进行FRFT数值计算之前必须对原始信号进行量纲归一化处理。 3.1.1量纲归一化原理如果一个信号在时间轴或频率轴的一个子集取非零值，并且取非零值的条件限定在有限区间内，则称该信号在时间轴或频率轴上是紧凑的。从理论上讲，任何信号和它的傅立叶变换不可能是同时紧凑的。然而我们实际中需要处理的信号往往是时限的和带限的。信号的时宽带宽积可以用来确定信号的采样频率和采样点数，用于唯一地从离散化的信号中恢复原始信号。假定原始连续信号在时间轴和频率轴上都是紧凑的，其时域表示限定在区间，而其频域表示限定在区间，和分别表示信号的时宽和带宽。信号的时宽带宽积为，根据不确定性定理，的值应当大于1。由于时域和频域具有不同的量纲，为了FRFT计算处理方便，需要将时域和频域分别转换成无量纲的域。具体方法是引入一个具有时间量纲的尺度因子，并定义新的尺度化坐标,新的坐标系实现了无量纲化。信号在新坐标系中被限定在区间和内。为使两个区间的长度相等，选择，则两个区间长度都等于无量纲量，即两个区间归一化为。归一化以后信号的Wigner-Ville分布限定在以原点为中心，直径的圆内，如图3-2所示。最后根据采样定理对归一化后的信号进行采样，采样间隔为，采样点数为.图3.1归一化后的时频支撑区域3.1.2分解算法 可以把（2-1）式改写为(3-1)式可以具体分解为以下几个步骤：（1）用chirp信号调制信号：(3-2)（2）调制信号与另一个chirp信号卷积：(3-3)（3）用chirp信号调制卷积后的信号：(3-4)要实现FRFT的数值计算必须对以上每个分解步骤都进行离散化处理，具体的实现过程如下：（1）对与chirp信号的乘积进行采样假定分数阶次，chirp信号的调频率，经chirp信号调制后所得的信号时宽带宽积可以是原信号的时宽带宽积的两倍，所以要求的采样间隔为，如果原来的采样间隔是，可通过插值的方法获得样本值，然后再chirp信号的离散采样相乘，以得到的采样值。（2）实现与chirp信号的卷积由于是带限信号，所以chirp信号也可取其带限形式。所以有：(3-5)其中而则是chirp信号傅立叶变换(3-6)于是，式的离散形式为(3-7) 这一离散卷积可利用FFT来计算。（3）计算分数阶傅立叶变换的采样值由于在第一步操作时对信号作了2倍内插操作，所以在最后的结果需要对再进行2倍抽取，以得到离散采样。总之上述方法从连续信号的N个离散采样开始，最后得到由的N个离散样本值得注意的是上述方法只适用于的情况，对位于该区间外的情况，可以利用分数阶傅立叶变换的周期性将阶次变到后再进行计算。3.1.3特征值和特征向量方法分解方法虽然计算复杂度较低，但不严格遵守分数阶傅立叶变换的旋转特性，因此不能从变换后的信号通过其逆变换精确地恢复原始信号。为了克服这种方法的缺点和不足，PeiSoo-Chang等人提出一种新的离散化方法，该方法具有与连续情况相似的变换性质和结果，并可以通过其逆变换恢复原始信号。这种方法就是特征值和特征向量方法，它从连续傅立叶变换的特征函数为Hermite函数出发，通过对Hermite函数的离散化近似和正交投影，得到一组与Hermite函数形状相似的DFT矩阵的正交化离散Hermite特征向量，然后，按照连续分数阶傅立叶变换的核函数谱分解表达式，构造离散分数阶傅立叶变换矩阵。3.1.3.1DFT矩阵特征值和特征向量傅立叶变换的特征函数为Hermite-Gaussian函数，其表达式为：(3-8)其中，为阶Hermite多项式。DFT矩阵的特征值及重复度如下表： DFT矩阵F的特征值是，共有1、-j、-1、j四个值，每个特征值对应的特征向量全体组成一个特征子空间，记为，每个特征值的重复度决定了子空间的秩。矩阵S可用于计算DFT矩阵F的特征向量，S的表达式为：(3-9)可以证明矩阵S和F满足乘法交换性，即SF=FS。因此矩阵S的特征向量也是矩阵F的特征向量，但它们对应不同的特征值。3.1.3.2DFT矩阵Hermite特征向量的计算由于DFT矩阵F的特征向量不唯一，我们希望从中找到与Hermite函数相似的特征向量，这样的向量称DFT矩阵Hermite特征向量，为了求Hermite特征向量，下面给出一系列的重要定理。定理1对DFT矩阵的Hermite特征向量而言，它对应的连续函数的扩展方差应当为是信号的采样间隔，连续Hermite函数采样后得到：(3-10)上式也可以看做是方差为1的Hermite-Gaussian函数，以为采样间隔进行采样得到的序列。定理2若序列是由单位方差Hermite-Gaussian函数经过采样到的，那么，可以证明下列近似等式成立：当N为偶数时 (3-11)当N为奇数时(3-12)定理3将序列按照以下方式平移得到当N为偶数时(3-13)当N为奇数时(3-14)则的离散傅立叶变换近似为，即当N足够大时(3-15)从定理2和定理3可以看出，Hermite函数的采样序列近似为DFT矩阵特征向量。对Hermite函数的采样序列作归一化，以记为：(3-16)通过S矩阵可以得到DFT矩阵F的一组实正交特征向量。因此可以将这些特征向量作为DFT特征子空间的基向量。然后计算向量在DFT特征子空间的投影，从而得到DFT矩阵的Hermite特征向量。计算方法如下：(3-17)3.1.3.3DFRFT核矩阵和二维DFRFT4最优变换阶次的选择本文选择最优变换阶次的方法主要是峰值搜索法，目前主要搜索算法有两种二维搜索算法、拟牛顿搜索算法。尽管拟牛顿算法比二维搜索的计算量有所减但是这两种算法的计算复杂度还是不能满足工程实际的实时处理要求,因此我们一维曲线拟合来代替二维搜索的峰值搜索算法,计算量大大简化。下面主要介绍几种算法。 4.1二维搜索算法传统的二维搜索法若以变换阶数P为变量进行扫描搜索的过程中,当参数估计的估计精度要求较高时,为了满足精度要求,变换阶数P必须选择小的搜索步长,这样就会成倍地增加计算的复杂度,如变换阶数在区间[0,2]上,若以步长0.0001进行搜索,则需要进行20000次的FRFT,计算量很大,通常采用步进搜索算法,即先以大步长再以小步长搜索最大值点,详细步骤为:先以0.01为步长进行阶数从0-2搜索最大值点,再以步长0.0001进行阶数从-0.01至+0.01搜索最大值点。这样一共需要进行400次FRFT运算,计算量有所减小。4.2拟牛顿搜索算法对于拟牛顿搜索算法,首先介绍一下拟牛顿迭代法的基本思想。采用牛顿法解非线性方程就是把非线性线性化的一种近似方法。把在点展为泰勒级数得,(4-1)取方程的线性部分作为非线性方程的近似方程,则有(4-2)假设则方程的解为(4-3)再把在点展开为泰勒级数,同样也取其线性部分作为非线性方程的近似方程,假设,则有(4-4)这样即可得到牛顿法的一个迭代序列：(4-5)(4-6)其中是第和是第n次搜索的结果,为第n次搜索的步长系数,为函数在点的尺度矩阵可以通过迭代的方法求得。牛顿迭代法中每次迭代的主要运算为一次一维搜索和函数 的一阶偏导数的计算,这一拟牛顿算法所需要的计算为分数阶傅里叶域的一次扫描搜索和一次迭代搜索,因为迭代搜索所需的搜索和函数的一阶偏导数的计算量远远小于扫描搜索,若扫描点数为m,信号的样本长度为N,则拟牛顿算法的计算复杂度为O(mNlgN),与其它基于二维时频分布的算法(运算量一般为)相比,拟牛顿算法计算量较小。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最优阶图 图4.3P=0.48，P=1.35，P=1.82和P=0.79分数阶傅立叶变换图 参考文献[1]平先军,陶然,周思永等.一种新的分数阶傅立叶变换快速算法[J].电子学报,2001,29(3):406-408[2]张立浩.分数阶傅立叶变换及其在通信中的应用[D].硕士学位论文,燕山大学,2012.[3]张兆祥.基于分数阶傅立叶变换的数字水印算法研究[D].硕士学位论文,华北电力大学,2007.[4]陈恩庆.基于分数阶fourier变换的ofdm系统关键技术研究[D].硕士学位论文,北京理工大学,2007.[5]王仁明,张春生,贾玉坤.基于分数阶傅里叶域均衡的无线测控系统研究[J].航天电子对抗,2012,28(3):47-50.[6]程乃平,席有猷,郝建华.基于分数阶傅里叶变换的LFM干扰抑制算法[J].装备学院学报,2014,25(1):74-77.[7]陈恩庆,陶然.一种基于分数阶傅里叶变换的OFDM系统及其均衡算法[J].电子学报,2007,35(3):410-414[8]殷敬伟,惠俊英,蔡平等.基于分数阶Fourier变换的水声信道参数估计[J].系统工程与电子技术,2007,29(10):1625-1627. 附录AfunctionFaf=frft(f,a)%ThefastFractionalFourierTransform%input:f=samplesofthesignal%a=fractionalpower%output:Faf=fastFractionalFouriertransformerror(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);%dospecialcasesif(a==0),Faf=f;return;end;if(a==2),Faf=flipud(f);return;end;if(a==1),Faf(shft,1)=fft(f(shft))/sN;return;endif(a==3),Faf(shft,1)=ifft(f(shft))*sN;return;end%reducetointerval0.52.0),a=a-2;f=flipud(f);endif(a>1.5),a=a-1;f(shft,1)=fft(f(shft))/sN;endif(a<0.5),a=a+1;f(shft,1)=ifft(f(shft))*sN;end%thegeneralcasefor0.5