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

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

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

3、可以在时频平面内以任意角度的旋转定义线性算子，记作，我们就可以把傅立叶变换推广到任意角度即分数阶傅立叶变换。图1.平面旋转角度变成平面分数阶傅立叶变换的定义为： (2-1)其中是核函数， (2-2)这里，p=1或-1时退化成为常规的傅立叶变换和逆变换。分数阶傅立叶变换和经典傅立叶变换具有以下的关系： 1.分数阶傅立叶变换是线性算子 2.周期性 (2-3) (2-4)1.2 分数傅立叶变换的性质（1） 线性分数阶傅立叶变换为线性变换，满足叠加原理：若和分别是原函数和的阶分数傅立叶变换，则有 (2-5)（2） 旋转相加行 (2-6)（3） 逆 (2-7)（4） 酋性 (2-8)（5） 交换性 (2

4、-9)（6） 结合性 (2-10)（7） 周期性 (2-11)（8） 特征函数 (2-12)（9） 卷积、相乘、相关 函数、在阶次 p 分数阶傅立叶域的卷积记作 (2-13)函数、 在阶次 p 分数阶傅立叶域的乘积记作 (2-14)函数、 在阶次 p 分数阶傅立叶域的相关定义为 (2-15)（10） 时移特性 (2-16)（11） 频移特性 (2-17)（12） 尺度特性 (2-18)其中（13） 微分特性 (2-19)（14） 积分特性 (2-20)3分数阶傅立叶变换的离散的离散算法分数阶傅立叶变换的出现引起了各个领域研究人员的重视，在工程上也有十分广阔的应用前景。在数字信号处理的应用中，必

5、须采用离散形式的分数阶傅立叶变换（DFRFT），这使得离散分数阶傅立叶变换及其快速算法的研究显得十分重要。目前，DFRFT 的离散化算法主要有四种：1. 利用来计算离散 FRFT 的核矩阵，从而利用 FFT 来计算DFRFT其中 W 是离散傅立叶变换核矩阵。这种方法实际上缺乏理论基础，而且其离散 FRFT 矩阵不满足连续 FRFT 的旋转相加性，因此不能用相同的方法计算逆 FRFT。实际计算所产生的误差比较大，与连续 FRFT 没有相似的输出结果。其计算复杂度与传统傅立叶变换相同，为。2. 分解方法根据 FRFT 的表达式，将 FRFT 分解为信号的卷积形式，从而利用 FFT 来计算FRFT。

6、这种方法思想比较直观，计算出的记过与连续 FRFT 的输出比较接近。但它要经过一次 2 倍内插和 2 倍抽取，而且还要进行坐标的无量纲化，实现起来较为烦琐。其计算复杂度为。3. 利用矩阵的特征值和特征向量来计算 DFRFT这种方法保持了连续 FRFT 的特征值-特征函数的关系，克服了第一种方法中特征值和特征向量不匹配的缺点。采用了两种正交映射的办法 DFT 的 Hermite 特征向量，由此开发出两种快速方法，即 OPA 方法和 GSA 方法，这两种方法都有和连续FRFT 相近的输出结果，可逆性好。其计算复杂度均为。4. 直接对 FRFT 进行离散化来计算 DFRFT。 这种方法采用直接将连续

7、 FRFT 离散化的方法来获得离散 FRFT 的核矩阵，避开了烦琐的特征值和特征向量匹配问题以及矩阵的正交归一化运算，与连续 FRFT有相似的输出结果，该算法的计算复杂度为。在目前的研究中，采用的最多的是分解方法和矩阵特征值和特征向量两种方法，下面的内容重点介绍这两种方法。3.1 分解方法所谓分解方法是指根据 FRFT 的表达式，将 FRFT 分解为信号的卷积形式，从而利用 FFT 来计算 FRFT。这种算法由 H.M.Ozaktas 等人提出，其计算速度几乎与FFT 相当，被公认为目前计算速度最快的一种 FRFT 数字计算方法，非常适合于对信号进行实时的 FRFT 计算。但这种快速算法的运算

8、机理决定了在进行 FRFT 数值计算之前必须对原始信号进行量纲归一化处理。3.1.1 量纲归一化原理如果一个信号在时间轴或频率轴的一个子集取非零值，并且取非零值的条件限定在有限区间内，则称该信号在时间轴或频率轴上是紧凑的。从理论上讲，任何信号和它的傅立叶变换不可能是同时紧凑的。然而我们实际中需要处理的信号往往是时限的和带限的。信号的时宽带宽积可以用来确定信号的采样频率和采样点数，用于唯一地从离散化的信号中恢复原始信号。假定原始连续信号在时间轴和频率轴上都是紧凑的，其时域表示限定在区间，而其频域表示限定在区间， 和 分别表示信号的时宽和带宽。信号的时宽带宽积为，根据不确定性定理，的值应当大于 1

9、。由于时域和频域具有不同的量纲，为了 FRFT 计算处理方便，需要将时域和频域分别转换成无量纲的域。具体方法是引入一个具有时间量纲的尺度因子，并定义新的尺度化坐标,新的坐标系实现了无量纲化。信号在新坐标系中被限定在区间和内 。 为 使 两 个 区 间 的 长 度 相 等 ， 选 择，则两个区间长度都等于无量纲量，即两个区间归一化为。归一化以后信号的 Wigner-Ville 分布限定在以原点为中心，直径 的圆内，如图 3-2 所示。最后根据采样定理对归一化后的信号进行采样，采样间隔为 ，采样点数为.图3.1 归一化后的时频支撑区域3.1.2分解算法可以把（2-1）式改写为 (3-1)式可以具体

10、分解为以下几个步骤：（1） 用 chirp 信号调制信号： (3-2)（2） 调制信号与另一个 chirp 信号卷积： (3-3)（3） 用 chirp 信号调制卷积后的信号： (3-4)要实现 FRFT 的数值计算必须对以上每个分解步骤都进行离散化处理，具体的实现过程如下：（1）对 与 chirp 信号的乘积进行采样 假定分数阶次 ，chirp 信号的调频率 ，经 chirp 信号调制后所得的信号时宽带宽积可以是原信号的时宽带宽积的两倍，所以要求 的采样间隔为，如果原来的采样间隔是 ，可通过插值的方法获得样本值，然后再 chirp 信号的离散采样相乘，以得到的采样值。（2）实现 与 chir

11、p 信号的卷积由于是带限信号，所以 chirp 信号也可取其带限形式。所以有： (3-5)其中而则是 chirp 信号傅立叶变换 (3-6)于是，式的离散形式为 (3-7)这一离散卷积可利用 FFT 来计算。（3）计算分数阶傅立叶变换的采样值 由于在第一步操作时对信号作了 2 倍内插操作，所以在最后的结果需要对再进行 2 倍抽取，以得到离散采样。总之上述方法从连续信号的 N 个离散采样开始，最后得到由的 N 个离散样本值得注意的是上述方法只适用于的情况，对位于该区间外的情况，可以利用分数阶傅立叶变换的周期性将阶次变到 后再进行计算。3.1.3 特征值和特征向量方法分解方法虽然计算复杂度较低，但

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

13、特征值和特征向量傅立叶变换的特征函数为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.2 DFT 矩阵 Hermite 特征向量的计算

14、由于 DFT 矩阵F的特征向量不唯一，我们希望从中找到与Hermite函数相似的特征向量，这样的向量称DFT矩阵Hermite特征向量，为了求Hermite 特征向量，下面给出一系列的重要定理。定理1 对 DFT 矩阵的 Hermite 特征向量而言，它对应的连续函数的扩展方差应当为是信号的采样间隔，连续 Hermite 函数采样后得到： (3-10)上式也可以看做是方差为 1 的 Hermite-Gaussian 函数，以为采样间隔进行采样得到的序列。定理2 若序列是由单位方差Hermite-Gaussian 函数经过采样到的，那么，可以证明下列近似等式成立：当 N 为偶数时 (3-11)当

15、 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.3

16、DFRFT 核矩阵和二维 DFRFT4 最优变换阶次的选择本文选择最优变换阶次的方法主要是峰值搜索法，目前主要搜索算法有两种二维搜索算法、拟牛顿搜索算法。尽管拟牛顿算法比二维搜索的计算量有所减但是这两种算法的计算复杂度还是不能满足工程实际的实时处理要求,因此我们一维曲线拟合来代替二维搜索的峰值搜索算法,计算量大大简化。下面主要介绍几种算法。4.1 二维搜索算法传统的二维搜索法若以变换阶数P为变量进行扫描搜索的过程中,当参数估计的估计精度要求较高时,为了满足精度要求,变换阶数P必须选择小的搜索步长,这样就会成倍地增加计算的复杂度,如变换阶数在区间0,2上,若以步长0.0001进行搜索,则需要进行

17、20000次的FRFT,计算量很大,通常采用步进搜索算法,即先以大步长再以小步长搜索最大值点,详细步骤为:先以0.01为步长进行阶数从0-2搜索最大值点,再以步长0.0001进行阶数从-0.01至+ 0.01搜索最大值点。这样一共需要进行400次FRFT运算,计算量有所减小。4.2拟牛顿搜索算法对于拟牛顿搜索算法,首先介绍一下拟牛顿迭代法的基本思想。采用牛顿法解非线性方程就是把非线性线性化的一种近似方法。把在点展为泰勒级数得, (4-1)取方程的线性部分作为非线性方程的近似方程,则有 (4-2)假设则方程的解为 (4-3)再把在点展开为泰勒级数,同样也取其线性部分作为非线性方程的近似方程,假设

18、,则有 (4-4)这样即可得到牛顿法的一个迭代序列： (4-5) (4-6)其中是第和是第n次搜索的结果,为第n次搜索的步长系数,为函数在点的尺度矩阵可以通过迭代的方法求得。牛顿迭代法中每次迭代的主要运算为一次一维搜索和函数的一阶偏导数的计算,这一拟牛顿算法所需要的计算为分数阶傅里叶域的一次扫描搜索和一次迭代搜索,因为迭代搜索所需的搜索和函数的一阶偏导数的计算量远远小于扫描搜索,若扫描点数为m,信号的样本长度为N,则拟牛顿算法的计算复杂度为O(mNlgN),与其它基于二维时频分布的算法(运算量一般为)相比,拟牛顿算法计算量较小。5 仿真验证本文中利用二维搜索法求最优阶次，该信号是一个chirp

19、信号以相加的形式调制简单的高斯信号的混合信号，对该信号进行求解它的最优阶数。编写的分数阶傅立叶变换函数见附录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.3 P=0.48，P=1.35，P=1.82和P=0.79分数阶傅立叶变换图参考文献1 平先军,陶然,周思永等. 一种新的分数阶傅立叶变换快速算法J. 电子学报,2001,29(3):406-408

20、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 陈恩庆, 陶然.一种基于分数阶傅里叶

21、变换的OFDM系统及其均衡算法J. 电子学报, 2007, 35(3):410-4148 殷敬伟,惠俊英,蔡平等. 基于分数阶Fourier变换的水声信道参数估计J. 系统工程与电子技术, 2007, 29(10):1625-1627.附录Afunction Faf = frft(f, a)% The fast Fractional Fourier Transform% input: f = samples of the signal% a = fractional power% output: Faf = fast Fractional Fourier transformerror(narg

22、chk(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 = 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

23、to interval 0.5 < a < 1.5if (a>2.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% the general case for 0.5 < a < 1.5alpha = a*pi/2;tana2 = tan(alpha/2);sina = sin(alpha);f = zeros(N-1,1) ; i

24、nterp(f) ; zeros(N-1,1);% chirp premultiplicationchrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.2);f = chrp.*f;% chirp convolutionc = pi/N/sina/4;Faf = fconv(exp(i*c*(-(4*N-4):4*N-4)'.2),f);Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);% chirp post multiplicationFaf = chrp.*Faf;% normalizing constantFaf = e

25、xp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);function xint=interp(x)% sinc interpolationN = length(x);y = zeros(2*N-1,1);y(1:2:2*N-1) = x;xint = fconv(y(1:2*N-1), sinc(-(2*N-3):(2*N-3)'/2);xint = xint(2*N-2:end-2*N+3);function z = fconv(x,y)% convolution by fftN = length(x(:);y(:)-1;P = 2nextpow2(N);z = ifft( fft(x,P) .* fft(y,P);z = z(1:N);附录Bclc;clear all;close all;fs=30;