精华内容
下载资源
问答
  • 离散分数阶傅里叶变换是功能强大的信号处理工具,在非平稳状态下具有广泛的应用信号。 在本文中,我们提出了一种稀疏离散分数阶傅里叶变换(SDFrFT)算法,以减少处理在分数阶傅里叶域中稀疏表示的大型数据集时的...
  • 目前网络关于分数阶傅里叶变换可靠的代码来源主要有一下几个: (1)来源一: M.A. Kutay. fracF: Fast computation of the fractional Fourier transform, 1996. www.ee.bilkent.edu.tr/~haldun/fracF.m. (2)...
  • 分数阶傅里叶变换作为傅里叶变换的广义形式,广泛地划分科学计算和研究,离散分数阶傅里叶变换是其恢复应用的关键。对一种基于数特征分解的方法进行了改进,并进行计算机仿真。从而提高了离散分数阶Fourier变换的...
  • 基于离散分数阶傅里叶变换的安全语音内容认证算法
  • 稀疏分数阶傅里叶变换及其在雷达运动目标检测中的应用
  • 基于中心类型DFT矩阵特征分解的MA-CDFRFT(Muhiangle Centered Discrete Fractional Fourier Transform)算法在计算一组离散分数阶傅里叶变换DFRFT(Discrete Fractional Fourier Transform)时充分利用FFT运算来...
  • 离散分数阶Fourier变换(DFRFT)算法FRFT 这篇文献发表于: 作者: 一、分数阶Fourier变换的定义 二、分数阶与其他时频分析工具( Wigner-Ville分布)的关系 三、离散分数阶傅立叶变换的计算 一、分数阶Fourier变换的...

    离散分数阶Fourier变换(DFRFT)算法FRFT 这篇文献发表于: 作者: 一、分数阶Fourier变换的定义 二、分数阶与其他时频分析工具( Wigner-Ville分布)的关系 三、离散分数阶傅立叶变换的计算 一、分数阶Fourier变换的定义 二、分数阶傅里叶变换与Wigner-Ville分布 首先,看一下Wigner-Ville分布 是 傅里叶变换 经过一系列变换后变为 由以上可得,等式的右边是 的Wigner-Ville分布, 左边是 的Wigner-Ville分布 也就是说 的Wigner-Ville分布, 是由 的Wigner-Ville分布旋转а角得到。 所以分数阶Fourier变换有一个重要的性质,分数阶Fourier变换是角度为α的时频面旋转. 这个性质建立起分数阶Fourier变换与时频分布间的直接联系, 并且为分数阶Fourier域理解为一种统一的时频变换域奠定了理论基础, 同时也为分数阶Fourier变换在信号处理领域中的应用提供了有利条件。 t ω u v α α 三、离散分数阶傅立叶变换的计算 目前DFRFT的四种离散化算法 在这篇文献中,第二种,采用分解的方法。 1.第一种分解方法 可以把以上改写为 假定p∈[-1,1],经过量纲归一化的信号x(t)的分数阶傅里叶变换, 可以分解为以下三个步骤: (1)用chirp信号调制信号f(x): (2)调制信号与另一个chirp信号卷积: (3)用chirp信号调制卷积后的信号: 式 1 式 2 式 3 具体细节:第一步:将函数, 与线性调频函数相乘(式1)。 注意,g(x)的频率带宽与时间带宽乘积可以是,f(x)的相应带宽乘积的两倍,所以要求g(x)的采样间隔为1/(2Δx)。如果,( )样本值的采样间隔是1/Δ x,那么就需要对这些样本值进行插值,然后再与线性调频函数的离散采样值相乘,以得到所希望的g(x)的采样。 第二步:将g(x)与一线性调频函数作卷积式(式 (2))。注意,由于g(x)是带限信号,所以线性调频函数也可以用其带限形式代替而不会有任何影响。 2、第二种分解方法 为了简化计算,人们提出更加有效的分解计算方法。 假定x(t)的wigner-ville分布限定在以原点为中心,直径为Δx的 圆内。若令 ,则与chirp信号乘积后的信号 在频域具有带宽Δx。可以用Shannon插值表示 简要介绍一下Shannon 插值 Shannon定理 到设信号 ,如果存在 ,使 , , 则称 是B频率截断的的,这时,只要采样间隔 按间隔 进行采样就不会损失信息,而且, 可按如下公式构造原信号 上式Shannon 插值公式。 利用采样序列 3、MATLAB程序 function Faf = 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 cases if (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 to interval 0.5 < a < 1.5 if (a>2.0), a = a-2; f = flipud(f); end if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end % the general case for 0.5 < a < 1.5 alpha = a*pi/2; tana2

    展开全文
  • 基于FPGA离散分数阶傅里叶变换算法的实现.pdf
  • 根据Ozaktas1996年的文章《Digital Computation of the Fractional Fourier Transform》编写的分数阶傅里叶变换程序。
  • frft(离散傅里叶分数变换)的matlab实现代码。利用其对单分量、多分量的LFM信号进行了检测和估计,结果验证了frft的正确性。
  • Face recognition systems suffers from illumination,expression and occlusion and so on.This paper presented a novel discrete fractional Fourier features method based on sparse principal component ...
  • 分数阶傅里叶变换(FRFT)之数值计算入门概要

    万次阅读 多人点赞 2018-05-25 21:48:58
    本文主要简述了作为一名初学者在对分数阶傅里叶变换的摸索过程中的一些浅显的看法,本文仅针对分数阶傅里叶变换离散计算部分,并没有对分数阶傅里叶变换中重要的思想方法和应用进行讨论。 1.关于代码 目前网络...

    分数阶傅里叶变换入门概要

    引言
    本文主要简述了作为一名初学者在对分数阶傅里叶变换的摸索过程中的一些浅显的看法,本文仅针对分数阶傅里叶变换的离散计算部分,并没有对分数阶傅里叶变换中重要的思想方法应用进行讨论。


    1.关于代码
    目前网络关于分数阶傅里叶变换可靠的代码来源主要有一下几个:
    (1)来源一:
    M.A. Kutay. fracF: Fast computation of the fractional Fourier transform, 1996.
    www.ee.bilkent.edu.tr/~haldun/fracF.m.
    (2)来源二:
    J. O’Neill. DiscreteTFDs:a collection of matlab files for time-frequency analysis, 1999.
    ftp.mathworks.com/pub/contrib/v5/signal/DiscreteTFDs/.
    (3)来源三
    Adhemar Bultheel and H´ector E. Mart´ınez Sulbaran.frft:Computation of the Fractional Fourier Transform,2004.
    www.cs.kuleuven.ac.be/~nalag/research/software/FRFT/.
    (4)来源四:
    A. Bultheel.frft2:A two-phase implementation of the fractional Fourier transform,2011.

    这里比较推荐的是来源三及其网站。2004年的这篇文章《Computation of the Fractional Fourier Transform》对前面的(1)和(2)进行了详细的对比,分析了两个程序的优劣并进行了改进,(3)的代码较为简洁也比较实用,其对应网站上还有很多相关的代码。文末附有下载链接,供参考。


    2.关于文献

    • 目前国内陶然教授研究得较为深入,其分别在04年和09年出版了《分数阶Fourier变换的原理与应用》、《分数阶傅里叶变换及其应用》,两本内容相似,参考其一即可。
    • L.B. Almeida. The fractional Fourier transform and time-frequency representation. IEEE Trans. Sig. Proc., 42:3084{3091, 1994.
      这篇文章是谈论分数阶傅里叶变换性质和基本思想的经典文章。文中详细介绍了分数阶傅里叶变化及其在信号处理中的应用。文中主要讨论的是连续的情况,文中讨论了分数阶傅里叶变换的性质,计算步骤,一些常见信号的分数阶傅里变换对,其与常见的WVD,短时傅里叶变换的关系。
    • H. Ozaktas, O. Arıkan, M. Kutay, G. Bozdagi, Digital computation of the fractional Fourier transform, IEEE Trans. Sig. Proc. 44 (9) (1996) 2141–2150.
      这篇文章是分数阶傅里叶变换数值计算的经典文章。文中给出了两种数值计算方法,但注意这里并不是真正的离散分数阶傅里叶变换,是用离散抽样近似连续的情况,然后用常见的 FFT F F T 方法处理。文中值得注意的是讨论了量纲归一化的作用。网上流传的很多matlab程序均基于此或者在此基础上改进。
    • A. Bultheel, H. Mart´ınez Sulbaran, An introduction to the fractional Fourier transform and friends, Cubo Matematica Educacional 7 (2) (2005) 201–221.
      这篇文章写的内容涵盖比较广,分析了分数阶傅里叶变换中的一些基本内容,分析了其与其他变换的关系,另外值得注意的是这里提到了线性正则变换。
    • A. Bultheel, H. Mart´ınez Sulbaran, Computation of the fractional Fourier transform, Appl. Comput. Harmonic Anal. 16 (3) (2004) 182–202.
      这篇文章如前所述,算是对程序分析比较清楚的一篇,详细研究matlab实现可以看。
    • Bultheel A. A two-phase implementation of the fractional Fourier transform[J]. Journal of Non Crystalline Solids, 2011, 352(9-20):VII.
      这篇文章针对特殊结构,在上文之上进行了一些改进。
    • 赵兴浩, 邓兵, 陶然. 分数阶傅里叶变换数值计算中的量纲归一化[J]. 北京理工大学学报, 2005, 25(4):360-364.
      这篇文章给出了量纲归一化的两种实用的实现方法。

    3.算法浅析
    本文主要讨论来源三中给出的算法,首先梳理一下基础知识。
    基本的连续型分数阶傅里叶变换形式如下

    fα(ξ)=Fα(ξ)=Kα(ξ,x)f(x)dx f α ( ξ ) = F α ( ξ ) = ∫ − ∞ ∞ K α ( ξ , x ) f ( x ) d x

    其中
    Kα(ξ,x)=Cαexp{iπ(2xξsinα(x2+ξ2)cotα)} K α ( ξ , x ) = C α exp ⁡ { − i π ( 2 x ξ sin ⁡ α − ( x 2 + ξ 2 ) cot ⁡ α ) }

    Cα=1icotα=exp{i[πsgn(sinα)/4α/2]}|sinα| C α = 1 − i cot ⁡ α = exp ⁡ { − i [ π s g n ( sin ⁡ α ) / 4 − α / 2 ] } | sin ⁡ α |

    为了方便计算,上式经整理可得到如下形式:
    fα(ξ)=CαTt(ξ)+Ts(ξx)[Tt(x)f(x)]dx f α ( ξ ) = C α T t ( ξ ) ∫ − ∞ + ∞ T s ( ξ − x ) [ T t ( x ) f ( x ) ] d x

    其中
    Tt(x)=exp{iπtx2} ,t=tan(α/2) ,s=csc(α) T t ( x ) = exp ⁡ { − i π t x 2 }   , t = tan ⁡ ( α / 2 )   , s = − csc ⁡ ( α )

    通过上式我们可以很明显的看出,要计算这样一个式子应该分为三步:
    1.将 f(x) f ( x ) 与线性调频信号 Tt T t 相乘
    2.将上述结果再与线性调频信号 Ts T s 进行卷积运算
    3.与线性调频信号 Tt T t 相乘

    上述的讨论均是在连续域,这种连续的变换在实际中并无法计算,通常需要对连续信号抽样和插值以进行数值计算。
    利用香农内插公式和数值积分运算可得,FRFT的离散计算公式如下:

    fα(k2Δ)Cα2Δexp[iπtan(α/2)(k2Δ)2]l=NN1exp[iπ(kl2Δ)2cscα]{exp{iπ(l2Δ)2tan(α/2)f(l2Δ)}} f α ( k 2 Δ ) ≈ C α 2 Δ exp ⁡ [ − i π tan ⁡ ( α / 2 ) ( k 2 Δ ) 2 ] ∑ l = − N N − 1 exp ⁡ [ i π ( k − l 2 Δ ) 2 csc ⁡ α ] { exp ⁡ { − i π ( l 2 Δ ) 2 tan ⁡ ( α / 2 ) f ( l 2 Δ ) } }

    上式看起来可能比较繁杂,若按照前面的分步计算方法整理下会清晰许多,如下
    fα(xk)xkCαkEt[ξk]l=NN1Es[xkl]{Et[xl]f[xl]} f α ( x k ) ≈ x k C α k E t [ ξ k ] ∑ l = − N N − 1 E s [ x k − l ] { E t [ x l ] f [ x l ] }

    其中
    Et[ξk]=exp{iπtξ2k} ,t=tan(α/2) ,s=csc(α) ,xk=ξk=k/2Δ  E t [ ξ k ] = exp ⁡ { − i π t ξ k 2 }   , t = tan ⁡ ( α / 2 )   , s = − csc ⁡ ( α )   , x k = ξ k = k / 2 Δ  

    这里的 Δ Δ 是量纲归一化后的时域或频域尺度。
    这里的离散计算流程同上,但区别是由于离散处理中数据是有限长的要考虑到频域的展宽效果,因此需要先插值再抽取。另外为了加快计算,其中的卷积可以用FFT进行计算,可以利用FFRT的一些性质将 α α 转换到[0.5,1.5)区间内进行计算,详细的细节讨论参见来源三的参考文献。


    4.代码参考

    分数阶傅里叶变换的四种可靠代码参考

    展开全文
  • 从网上下载了frft的matlab程序,没有进行归一化处理,在仿真的过程中角度是对的,的,所以斜率是对的,但是初始频率一直是错的,请问是我的方法出现了什么问题吗? 附程序段: ph=0; %正弦信号的初始相位 ...
  • 水下目标散射回波在时域、频域混叠在一起, 而且... 给出了离散分数阶傅里叶变换对声散射分量的分辨能力和计算精度与发射信号带宽和观测时间之间的关系. 实验数据处理表明, 建立的分数阶傅里叶变换域的全方位模型与目标
  • 满足在数字信号处理器D SP (d ig ita l signa l p ro ce sso r) 上进行离散分数阶傅里叶变换 D FR F T (d isc re te f rac t io na l fo u r ie r t ran sfo rm ) 实时计算的要求, 通过对多种D FR F T 计算方法进 行...
  • 适合分数阶傅里叶变换的初学者,学习实现frft的Python编码参考,仅针对分数阶傅里叶变换离散计算部分。很适合处理非平稳信号,尤其是chirp类信号,具有良好的时频域特性。通过对分数域中的图像能量和幅度相位分布...
  • 本文主要是基于Haldun M. Ozaktas, Orhan等人的论文Digital Computation of ...我们还讨论了离散分数傅里叶变换的定义。 引言 略…… 预备知识 分数傅里叶变换 定义{Ff}(x)\{\mathcal{F}f\}(x){Ff}(x)表示对f(x)f(x)f

    本文主要是基于Haldun M. Ozaktas, Orhan等人的论文Digital Computation of the Fractional Fourier Transform 翻译而成。如有错误之处还请指出。

    摘要

    本文给出了一种高效、精确计算分数傅里叶变换的算法。对于时间带宽积为N的信号,该算法计算复杂度为O(N logN)。我们还讨论了离散分数傅里叶变换的定义。

    引言

    略……

    预备知识

    分数傅里叶变换

    定义 { F f } ( x ) \{\mathcal{F}f\}(x) {Ff}(x)表示对 f ( x ) f(x) f(x)做一次傅里叶变换, { F a f } ( x ) \{\mathcal{F}^af\}(x) {Faf}(x)表示连续对 f ( x ) f(x) f(x) a a a次傅里叶变换。
    根据傅里叶变换的定义,有 { F 2 f } ( x ) = f ( − x ) \{\mathcal{F}^2f\}(x)=f(-x) {F2f}(x)=f(x) { F 4 f } ( x ) = f ( x ) \{\mathcal{F}^4f\}(x)=f(x) {F4f}(x)=f(x)(注:使用信号与系统里的傅里叶变换定义公式推导时可能会差 2 π 2\pi 2π倍,需要将 ω \omega ω换成 2 π f 2\pi f 2πf并以 f f f作为积分变量计算),可见傅里叶变换周期为 a = 4 a=4 a=4,于是设定 a a a的定义域为 0 < ∣ a ∣ < 2 0<|a|<2 0<a<2,负数时表示逆傅里叶变换。

    a a a次傅里叶变换用卷积表示:
    在这里插入图片描述
    其中 ≡ \equiv 表示恒等于, B a ( x , x ′ ) B_a(x,x') Ba(x,x)表示卷积核,注意 x ′ x' x并不是表示 x x x的导数,只是一个新的变量,用于区分 x x x
    这个卷积核有前人推导过了,表达式如下:
    在这里插入图片描述
    其中 i i i表示虚数单位。

    a = 0 a=0 a=0时, B a ( x , x ′ ) = δ ( x − x ′ ) B_a(x,x')=\delta(x-x') Ba(x,x)=δ(xx) δ \delta δ表示冲激函数,显然,根据冲激函数定义,当且仅当 x ′ = x x'=x x=x时,积分有非零值,因此 { F 0 f } ( x ) = f ( x ) \{\mathcal{F}^0f\}(x)=f(x) {F0f}(x)=f(x).
    a = ± 2 a=\pm2 a=±2时, B a ( x , x ′ ) = δ ( x + x ′ ) B_a(x,x')=\delta(x+x') Ba(x,x)=δ(x+x),同理,此时 { F ± 2 f } ( x ) = f ( − x ) \{\mathcal{F}^{\pm2}f\}(x)=f(-x) {F±2f}(x)=f(x).
    a = 1 a=1 a=1时, A ϕ = 1 A_\phi=1 Aϕ=1, B a ( x , x ′ ) = e − 2 π i x x ′ B_a(x,x')=e^{-2\pi i x x'} Ba(x,x)=e2πixx,卷积核为傅里叶变换核, { F 1 f } ( x ) = F ( x ) \{\mathcal{F}^1f\}(x)=\mathcal{F}(x) {F1f}(x)=F(x)
    a = − 1 a=-1 a=1时, A ϕ = 1 A_\phi=1 Aϕ=1, B a ( x , x ′ ) = e 2 π i x x ′ B_a(x,x')=e^{2\pi i x x'} Ba(x,x)=e2πixx,卷积核为逆傅里叶变换核, { F − 1 f } ( x ) = F − 1 ( x ) \{\mathcal{F}^{-1}f\}(x)=\mathcal{F}^{-1}(x) {F1f}(x)=F1(x)

    傅里叶变换对如下:
    F ( x ) = F ( f ( x ′ ) ) = ∫ − ∞ + ∞ f ( x ′ ) e − 2 π i x x ′ d x ′ F(x)=\mathcal{F}(f(x'))=\int_{-\infty}^{+\infty}f(x')e^{-2\pi ixx'}dx' F(x)=F(f(x))=+f(x)e2πixxdx
    f ( x ′ ) = F − 1 ( F ( x ) ) = ∫ − ∞ + ∞ F ( x ) e 2 π i x x ′ d x f(x')=\mathcal{F}^{-1}(F(x))=\int_{-\infty}^{+\infty}F(x)e^{2\pi ixx'}dx f(x)=F1(F(x))=+F(x)e2πixxdx
    (与平常课本上学的在形式上不太一样,不过物理意义是相同的,只不过是把角频率换成了频率)

    可以证明有如下性质成立:
    F a 1 F a 2 = F a 1 + a 2 \mathcal{F}^{a_1}\mathcal{F}^{a_2}=\mathcal{F}^{a_1+a_2} Fa1Fa2=Fa1+a2

    傅里叶变换的核函数,其完备集合是Hermite-Gaussian 函数:
    在这里插入图片描述
    (这块不细讲了)
    为方便起见, { F a f } ( x ) \{\mathcal{F}^{a}f\}(x) {Faf}(x)简写为 f a ( x ) f_a(x) fa(x)

    和维纳分布的关系以及分数阶傅里叶域的概念

    如果时域函数为 f ( x ) f(x) f(x),则维纳分布函数为:
    在这里插入图片描述
    粗略来讲,维纳分布函数给出了信号在时域和频域上的能量分布,有如下特点:
    在这里插入图片描述
    分数阶傅里叶变化的出现,丰富了时域和频域的概念,以时域中的时间为x轴,频域中的频率为y轴,可以画出时频域平面,即分数傅里叶域。维纳分布函数也可以通过几何的方法,改写为:
    在这里插入图片描述
    在这里插入图片描述
    分数阶傅里叶变换突破了时域和频域的概念,可以做到任意角度的变换。

    时域、频域和维格纳空间中的紧凑型

    略……

    离散傅里叶变换

    在这里插入图片描述

    计算连续分数阶傅里叶变换的方法

    上述计算连续分数阶傅里叶变换的公式,很难用解析的方法求解,因此我们通常使用数值计算方法,比如用指数或者二次函数近似,但是这样通常情况下需要大量的样本。当 a a a接近 0 0 0 ± 2 \pm2 ±2时,情况更加严重。

    有一种快速计算方法是利用 F a = F 1 F a − 1 \mathcal{F}^a=\mathcal{F}^1\mathcal{F}^{a-1} Fa=F1Fa1,其中 F a − 1 \mathcal{F}^a-1 Fa1可以直接得到。
    还有一种方法用到了核函数的频谱分解。
    尽管这两种方法都期望能得到正确的结果,但是我们不会进一步考虑他们,因为他们的计算复杂度为 O ( N 2 ) O(N^2) O(N2)

    分数阶傅里叶变换的快速算法

    分数阶傅里叶变换是一类更普遍的变换,通常被称为linear canonical transformations 或 quadratic-phase transforms。这些变换通常能分解成一系列简单的操作,比如线性调频卷积、线性调频相乘、缩放和原始傅里叶变换。这里我们主要研究两种分解方法,对应两种独特的方法。

    方法一

    首先,我们将分数阶傅里叶变换分解成一个线性调频相乘跟着一个线性调频卷积跟着另一个线性调频相乘。在这个方法中,假设 a ∈ [ − 1 , 1 ] a\in [-1,1] a[1,1],可以将分数阶傅里叶变换公式变为:
    在这里插入图片描述
    在这里插入图片描述
    这里 g ( x ) g(x) g(x) g ′ ( x ) g'(x) g(x)代表中间结果,其中 ′ ' 不表示导数, β = csc ⁡ ϕ \beta=\csc\phi β=cscϕ

    在第一步(公式16)中,我们把 f ( x ) f(x) f(x)乘以了一个线性调频函数,首先我们要获得 g ( x ) g(x) g(x)的采样,需要内插。接下来是把 g ( x ) g(x) g(x)卷积上一个线性调频函数。注意到 g ( x ) g(x) g(x)带宽有限,因此线性调频函数可以等价为一个带宽有限函数:
    在这里插入图片描述
    注意到式(19)是线性调频函数 exp ⁡ ( i π β x 2 ) \exp(i\pi \beta x^2) exp(iπβx2)的傅里叶变换。 h ( x ) h(x) h(x)可以认为是菲涅尔积分:
    在这里插入图片描述
    g ′ ( x ) g'(x) g(x)被采样后,结果为:
    在这里插入图片描述
    这个卷积可以通过快速傅里叶变换计算。
    目前采样间隔为 1 / 2 Δ x 1/2\Delta x 1/2Δx,需要使用抽删操作,大批删减样本,使样本间隔变为 1 / Δ x 1/\Delta x 1/Δx

    总而言之, f ( x ) f(x) f(x) f a ( x ) f_a(x) fa(x)都是以 1 / Δ x 1/\Delta x 1/Δx为间隔采样的,均为 N N N个样本,,将这些样本以列向量表示,分别为 f \mathbf{f} f f a \mathbf{f}_a fa,则总程序可以表示为:
    在这里插入图片描述
    D \mathbf{D} D J \mathbf{J} J表示抽删和内插操作的矩阵, Λ \mathbf{\Lambda} Λ表示线性调频相乘, H l p \mathbf{H}_{lp} Hlp对应卷积操作。

    如果我们只关心计算和绘制给定的连续函数 f ( x ) f(x) f(x)的分数阶傅里叶变换,那么抽删和内插操作可以省去。
    注意 a a a的定义域是 − [ 1 , 1 ] -[1,1] [1,1],如果 a a a在区域之外,可以利用分数阶傅里叶变化的对称和周期性将 a a a转换到定义域内。

    方法二

    接下来我们关注一个改进的方法,这个方法不需要菲涅尔积分。分数阶傅里叶变换公式可以写为:
    在这里插入图片描述
    其中 α = cot ⁡ ϕ , β = csc ⁡ ϕ \alpha=\cot\phi,\beta=\csc\phi α=cotϕ,β=cscϕ

    在一堆看不懂的假设后, e i π α x ′ 2 f ( x ′ ) e^{i\pi\alpha x'^{2}}f(x') eiπαx2f(x)可以用香农内插公式书写:
    在这里插入图片描述
    其中 N = ( Δ x ) 2 N=(\Delta x)^2 N=(Δx)2,改变求和顺序后可以得到:
    在这里插入图片描述
    里面的积分等于
    在这里插入图片描述
    窗函数可以去掉,公式写为
    在这里插入图片描述
    采样后的变换函数可以写为:
    在这里插入图片描述
    这是一个有限的求和,允许我们依照原函数的采样求得分数傅里叶变换的采样。直接计算需要 O ( N 2 ) O(N^2) O(N2)次乘法,但是接下来我们介绍如何用 O ( N log ⁡ N ) O(N\log N) O(NlogN)的时间复杂度计算。把上式通过代数变换变为:
    在这里插入图片描述
    可以很明显的看到,求和符号里面是关于 e i π β ( n / 2 Δ x ) e^i\pi\beta(n/2\Delta x) eiπβ(n/2Δx)和线性调频调制后的 f ( x ) f(x) f(x)线性卷积。线性卷积可以使用快速傅里叶变换(fft)计算得到,而fft的计算复杂度为 O ( N log ⁡ N ) O(N\log N) O(NlogN)
    得出的结果样本在乘以一个线性调频函数,就是最终结果。
    计算过程用矩阵和向量总结如下:
    在这里插入图片描述
    注意 ∣ m ∣ , ∣ n ∣ < N |m|,|n|<N m,n<N
    在这个方法中假设了 0.5 ≤ ∣ a ∣ ≤ 1.5 0.5\leq|a|\leq1.5 0.5a1.5,不过我们可以使用分数傅里叶变换的可加性将其拓展到任意定义域,比如对于范围 0 ≤ ∣ a ∣ ≤ 0.5 0\leq|a|\leq0.5 0a0.5,可以使用以下变换式:
    在这里插入图片描述
    注意最后一项不是表示乘积,是表示连续做两次变换。

    在这种计算方法中,分数阶傅里叶变换和原傅里叶变换之间的关系为:
    在这里插入图片描述
    可以看出
    在这里插入图片描述
    其中 F \mathbf{F} F表示原傅里叶变换矩阵.


    代码与注释

    论文的原理部分翻译到此结束了,下面是按照原论文的思想,参考了其它博主的代码实现,加了一些注释便于理解。

    function Faf = myfrft(f, a)
    %分数阶傅里叶变换函数
    %输入参数f为原始信号,a为阶数
    %输出结果为原始信号的a阶傅里叶变换
    N = length(f);%总采样点数
    shft = rem((0:N-1)+fix(N/2),N)+1;%此项等同于fftshift(1:N),起到翻转坐标轴的作用
    sN = sqrt(N);%看原文中对离散傅里叶变换的定义,有这个乘积项
    a = mod(a,4);%按周期性将a定义在[0,4]
    
    %特殊情况直接处理
    if (a==0), Faf = f; return; end%自身
    if (a==2), Faf = flipud(f); return; end%f(-x)
    if (a==1), Faf(shft,1) = fft(f(shft))/sN; return; end%f的傅里叶变换
    if (a==3), Faf(shft,1) = ifft(f(shft))*sN; return; end%f的逆傅里叶变换
    
    %利用叠加性将阶数变换到0.5 < a < 1.5
    if (a>2.0), a = a-2; f = flipud(f); end%a=2是反转
    if (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end%a=1是傅里叶变换
    if (a<0.5), a = a+1; f(shft,1) = ifft(f(shft))*sN; end%a=-1是逆傅里叶变换
    
    %开始正式的变换
    alpha = a*pi/2;
    tana2 = tan(alpha/2);
    sina = sin(alpha);
    
    f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];%使用香农插值,拓展为4N
    % 以下操作对应原论文中公式(29% 线性调频预调制
    chrp = exp(-1i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);
    f = chrp.*f;
    % 线性调频卷积
    c = pi/N/sina/4;
    Faf = fconv(exp(1i*c*(-(4*N-4):4*N-4)'.^2),f);
    Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi);
    % 线性调频后调制
    Faf = chrp.*Faf;
    % 乘以最前面的A_Phi项
    Faf = exp(-1i*(1-a)*pi/4)*Faf(N:2:end-N+1);
    
    end
    
    function xint=interp(x)%香农插值
    % sinc interpolation
    N = 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);
    end
    
    function z = fconv(x,y)%利用fft快速计算卷积
    N = length([x(:);y(:)])-1;%计算最大点数
    P = 2^nextpow2(N);%补零
    z = ifft( fft(x,P) .* fft(y,P));%频域相乘,时域卷积
    z = z(1:N);%去零
    end
    

    测试

    使用如下代码进行测试:

    close all
    a=0:0.25:4;%分数阶傅里叶变换阶数
    
    %生成一个窗函数
    fx=zeros(500,1);
    fx(150:250)=1;
    
    for ai=a
        figure
        F=myfrft(fx,ai);
        plot(abs(F))
        title("a="+num2str(ai))
        grid on
        ylim([0,5])
    end
    

    首先生成一个窗函数,长度为500,在150到250点处为1,其余为0。
    学过信号与系统的同学都知道,窗函数的傅里叶变换是一个Sa函数(也称Sinc函数),选取傅里叶变换的阶数为0到4之间,间隔为0.25,观察每一个结果。我将分数阶傅里叶变换的结果做成了一个动图,可以更加直观的理解分数阶傅里叶变换的物理意义和时频域的概念。
    请添加图片描述

    展开全文
  • 分数阶傅里叶变换离散算法,通过将定义式分解成三步实现的采样型算法。
  • 傅里叶变换不论在数学界还是在工程领域都是一种非常重要的分析工具,一位傅里叶变换的表达式为(1-1-1),其中f(t)f(t)f(t)是一个属于能量有限的信号空间的信号,傅里叶变换是一种可逆变换,其逆变换的表达式为(1-1-2...

    一 引言

    ​ 傅里叶变换不论在数学界还是在工程领域都是一种非常重要的分析工具,一位傅里叶变换的表达式为(1-1-1),其中 f ( t ) f(t) f(t)是一个属于能量有限的信号空间的信号,傅里叶变换是一种可逆变换,其逆变换的表达式为(1-1-2),有的文献逆变换前的系数为 1 2 π \frac 1{2\pi} 2π1,本文的定义是为了表现出傅里叶变换的正交性,因此对一个信号做傅里叶变换可以看作是把该信号与复指数信号做内积(投影)的结果.
    F ( w ) = 1 2 π ∫ − ∞ ∞ f ( t ) e − j w t d t (1-1-1) F(w)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty f(t)e^{-jwt} dt \tag{1-1-1} F(w)=2π 1f(t)ejwtdt(1-1-1)

    f ( t ) = 1 2 π ∫ − ∞ ∞ F ( w ) e j w t d w (1-1-2) f(t)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty F(w)e^{jwt} dw \tag{1-1-2} f(t)=2π 1F(w)ejwtdw(1-1-2)

    因此其具有鲜明的物理意义,即看一个信号中都包含了什么频率成分,以及各个成分所占的比重为多少.今天我们来聊一聊分数傅里叶变换,分数傅里叶变换有多种的定义形式,我们来聊聊台湾学者 c . c . s h i h c.c.shih c.c.shih所定义的分数傅立叶变换.对于分数傅里叶变换,读者可以这么理解:如果把信号的频域看作是信号做一次傅里叶变换所得到的结果,把信号的时域看作是信号做零次傅里叶变换的结果,那么对一个信号做一次分数傅里叶变换可以看作是信号由时域到频域过渡的状态,即信号的时频信息.至于为什么要研究分数阶的傅里叶变换,在以后的文章会进行讲解.本文的主要目的是探究 c . c . s h i h c.c.shih c.c.shih所定义的分数傅立叶变换的表达形式以及其离散快速算法的实现


    二 分数傅里叶变换及其离散算法

    1. c . c . s h i h c.c.shih c.c.shih分数傅立叶变换的原理

      为了说清楚 c . c . s h i h c.c.shih c.c.shih所定义的分数傅立叶变换,先引入傅里叶变换的’周期性’,首先定义 f i ( x ) , i = 0 , 1 , 2 , 3 f_i(x),i=0,1,2,3 fi(x),i=0,1,2,3为对原始信号 f ( t ) f(t) f(t)进行了 i i i次的傅里叶变换的结果, x x x代表着自变量要么是时间 t t t,要么是角频率 w w w.我们知道,对一个信号做傅里叶变换的表达式为(2-1-1),现在我们做一次有趣的实验,对一次傅里叶变换所得的信号再做一次傅里叶变换,见式(2-1-2),可以发现结果是时域信号的翻转.
      f 1 ( x ) = 1 2 π ∫ − ∞ ∞ f ( t ) e − j w t d t = F ( w ) (2-1-1) f_1(x)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty f(t)e^{-jwt} dt=F(w) \tag{2-1-1} f1(x)=2π 1f(t)ejwtdt=F(w)(2-1-1)

      f 2 ( x ) = 1 2 π ∫ − ∞ ∞ F ( w ) e − j w t d w = f ( − t ) (2-1-2) f_2(x)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty F(w)e^{-jwt} dw=f(-t) \tag{2-1-2} f2(x)=2π 1F(w)ejwtdw=f(t)(2-1-2)

      那么对这次的结果再做一次傅里叶变换,见式(2-1-3),可以发现结果式频域信号的翻转.如果对这一次的结果再做一次傅里叶变换呢?见式(2-1-4),我们竟然得到了原始信号本身,这意味着对一个信号连续做4次的傅里叶变换,最后竟然会回到信号本身.
      f 3 ( x ) = 1 2 π ∫ − ∞ ∞ f ( − t ) e − j w t d t = F ( − w ) (2-1-3) f_3(x)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty f(-t)e^{-jwt} dt=F(-w) \tag{2-1-3} f3(x)=2π 1f(t)ejwtdt=F(w)(2-1-3)

      f 4 ( x ) = 1 2 π ∫ − ∞ ∞ F ( − w ) e − j w t d w = f ( t ) (2-1-4) f_4(x)=\frac 1{\sqrt[]{2\pi}}\int_{-\infty}^\infty F(-w)e^{-jwt} dw=f(t) \tag{2-1-4} f4(x)=2π 1F(w)ejwtdw=f(t)(2-1-4)

      c . c . s h i h c.c.shih c.c.shih的观点是:如果把 f 0 ( x ) , f 1 ( x ) , f 2 ( x ) , f 3 ( x ) f_0(x),f_1(x),f_2(x),f_3(x) f0(x),f1(x),f2(x),f3(x)当作基本态,那么将分数傅里叶变换定义为这4种状态的线性叠加,就定义了一种新的分数傅里叶变换,见式(2-1-5),叠加的系数是阶数的函数.
      F α ( x ) = A 0 ( α ) f 0 ( x ) + A 1 ( α ) f 1 ( x ) + A 2 ( α ) f 2 ( x ) + A 3 ( α ) f 3 ( x ) (2-1-5) F^{\alpha}(x)=A_0({\alpha})f_0(x)+A_1({\alpha})f_1(x)+A_2({\alpha})f_2(x)+A_3({\alpha})f_3(x) \tag{2-1-5} Fα(x)=A0(α)f0(x)+A1(α)f1(x)+A2(α)f2(x)+A3(α)f3(x(2-1-5)
      定义符号 A j ( α ) A_j(\alpha) Aj(α)为代表着第j个和阶数 α \alpha α有关的权重系数.同时我们也可以得到当阶数是整数的时候,系数的表达式为式(2-1-6),其中 j = 0 , 1 , 2 , 3 , l = 0 , 1 , 2 , 3 j=0,1,2,3,l=0,1,2,3 j=0,1,2,3,l=0,1,2,3.那么新的问题来了,这些和阶数相关的系数应该如何求解呢? c . c . s h i c.c.shi c.c.shi给出了相关定理,通俗地讲,对一个信号先做 α \alpha α阶的分数傅里叶变换再做 β \beta β阶的分数傅里叶变换所得到的结果应该是对信号做 ( α + β ) (\alpha+\beta) (α+β)阶的分数傅里叶变换,同时也是对一个信号先做 β \beta β阶的分数傅里叶变换再做 α \alpha α阶的分数傅里叶变换,见式(2-1-7),这样给人的感觉很直观.对于任意的 α \alpha α β \beta β,分数傅里叶变换 F s α F_s^{\alpha} Fsα具有可加性.
      A j ( l ) = δ ( j − l ) (2-1-6) A_j(l)=\delta(j-l)\tag{2-1-6} Aj(l)=δ(jl)(2-1-6)

      F s α F s β = F s β F s α = F s α + β (2-1-7) F_s^{\alpha}F_s^{\beta}=F_s^{\beta}F_s^{\alpha}=F_s^{\alpha+\beta} \tag{2-1-7} FsαFsβ=FsβFsα=Fsα+β(2-1-7)

      于是在满足上述定理的情况下,于是可以得到函数方程组(2-1-8),该方程组满足式(2-1-6)的边界条件.我们将其写成矩阵的表达形式,见(2-1-9).
      { A 0 ( α + β ) = A 0 ( α ) A 0 ( β ) + A 1 ( α ) A 3 ( β ) + A 2 ( α ) A 2 ( β ) + A 3 ( α ) A 1 ( β ) A 1 ( α + β ) = A 0 ( α ) A 1 ( β ) + A 1 ( α ) A 0 ( β ) + A 2 ( α ) A 3 ( β ) + A 3 ( α ) A 2 ( β ) A 2 ( α + β ) = A 0 ( α ) A 2 ( β ) + A 1 ( α ) A 1 ( β ) + A 2 ( α ) A 0 ( β ) + A 3 ( α ) A 3 ( β ) A 3 ( α + β ) = A 0 ( α ) A 3 ( β ) + A 1 ( α ) A 2 ( β ) + A 2 ( α ) A 1 ( β ) + A 3 ( α ) A 0 ( β ) (2-1-8) \begin{cases} A_0(\alpha+\beta)=A_0(\alpha)A_0(\beta)+A_1(\alpha)A_3(\beta)+A_2(\alpha)A_2(\beta)+A_3(\alpha)A_1(\beta)\\ A_1(\alpha+\beta)=A_0(\alpha)A_1(\beta)+A_1(\alpha)A_0(\beta)+A_2(\alpha)A_3(\beta)+A_3(\alpha)A_2(\beta)\\ A_2(\alpha+\beta)=A_0(\alpha)A_2(\beta)+A_1(\alpha)A_1(\beta)+A_2(\alpha)A_0(\beta)+A_3(\alpha)A_3(\beta)\\ A_3(\alpha+\beta)=A_0(\alpha)A_3(\beta)+A_1(\alpha)A_2(\beta)+A_2(\alpha)A_1(\beta)+A_3(\alpha)A_0(\beta) \end{cases} \tag{2-1-8} A0(α+β)=A0(α)A0(β)+A1(α)A3(β)+A2(α)A2(β)+A3(α)A1(β)A1(α+β)=A0(α)A1(β)+A1(α)A0(β)+A2(α)A3(β)+A3(α)A2(β)A2(α+β)=A0(α)A2(β)+A1(α)A1(β)+A2(α)A0(β)+A3(α)A3(β)A3(α+β)=A0(α)A3(β)+A1(α)A2(β)+A2(α)A1(β)+A3(α)A0(β)(2-1-8)

    [ A 0 ( α + β ) A 1 ( α + β ) A 2 ( α + β ) A 3 ( α + β ) ] = [ A 0 ( β ) A 3 ( β ) A 2 ( β ) A 1 ( β ) A 1 ( β ) A 0 ( β ) A 3 ( β ) A 2 ( β ) A 2 ( β ) A 1 ( β ) A 0 ( β ) A 3 ( β ) A 3 ( β ) A 2 ( β ) A 1 ( β ) A 0 ( β ) ] [ A 0 ( α ) A 1 ( α ) A 2 ( α ) A 3 ( α ) ] (2-1-9) \begin{bmatrix} {A_0(\alpha+\beta)}\\ {A_1(\alpha+\beta)}\\ {A_2(\alpha+\beta)}\\ {A_3(\alpha+\beta)}\\ \end{bmatrix}= \begin{bmatrix} {A_0(\beta)}&{A_3(\beta)}&{A_2(\beta)}&{A_1(\beta)}\\ {A_1(\beta)}&{A_0(\beta)}&{A_3(\beta)}&{A_2(\beta)}\\ {A_2(\beta)}&{A_1(\beta)}&{A_0(\beta)}&{A_3(\beta)}\\ {A_3(\beta)}&{A_2(\beta)}&{A_1(\beta)}&{A_0(\beta)}\\ \end{bmatrix} \begin{bmatrix} {A_0(\alpha)}\\ {A_1(\alpha)}\\ {A_2(\alpha)}\\ {A_3(\alpha)}\\ \end{bmatrix} \tag{2-1-9} A0(α+β)A1(α+β)A2(α+β)A3(α+β)=A0(β)A1(β)A2(β)A3(β)A3(β)A0(β)A1(β)A2(β)A2(β)A3(β)A0(β)A1(β)A1(β)A2(β)A3(β)A0(β)A0(α)A1(α)A2(α)A3(α)(2-1-9)

    ​ 观察由阶数 β \beta β确定的权重所组成的矩阵,可以发现矩阵方程左侧元素下标之和实际上等于矩阵方程右侧对应对应元素下 标之和模4的结果.这不正是两个序列(2-1-10)与(2-1-11)做4点 D F T DFT DFT的结果?
    ( A 0 ( α ) , A 1 ( α ) , A 2 ( α ) , A 3 ( α ) ) (2-1-10) (A_0(\alpha),A_1(\alpha),A_2(\alpha),A_3(\alpha))\tag{2-1-10} (A0(α),A1(α),A2(α),A3(α))(2-1-10)

    ( A 0 ( β ) , A 1 ( β ) , A 2 ( β ) , A 3 ( β ) ) (2-1-11) (A_0(\beta),A_1(\beta),A_2(\beta),A_3(\beta))\tag{2-1-11} (A0(β),A1(β),A2(β),A3(β))(2-1-11)
    ​ 我们还知道,卷积运算在傅里叶变换下都成乘积运算.于是我们定义新的权重序列(见式2-1-12),实际上从序列 A j ( α ) A_j(\alpha) Aj(α)到序 列 B j ( α ) B_j(\alpha) Bj(α)也就是对序列 A j ( α ) A_j(\alpha) Aj(α)做了4点的 D F T DFT DFT,因此式(2-1-9)改写成式(2-1-13),边界条件为式(2-1-14).
    [ B 0 ( α ) B 1 ( α ) B 2 ( α ) B 3 ( α ) ] = 1 N [ 1 1 1 1 1 − i − 1 i 1 − 1 1 − 1 1 i − 1 − i ] [ A 0 ( α ) A 1 ( α ) A 2 ( α ) A 3 ( α ) ] (2-1-12) \begin{bmatrix} {B_0(\alpha)}\\ {B_1(\alpha)}\\ {B_2(\alpha)}\\ {B_3(\alpha)}\\ \end{bmatrix}= \frac{1}{\sqrt[]N} \begin{bmatrix} {1}&{1}&{1}&{1}\\ {1}&{-i}&{-1}&{i}\\ {1}&{-1}&{1}&{-1}\\ {1}&{i}&{-1}&{-i}\\ \end{bmatrix} \begin{bmatrix} {A_0(\alpha)}\\ {A_1(\alpha)}\\ {A_2(\alpha)}\\ {A_3(\alpha)}\\ \end{bmatrix} \tag{2-1-12} B0(α)B1(α)B2(α)B3(α)=N 111111i1i11111i1iA0(α)A1(α)A2(α)A3(α)(2-1-12)

    [ B 0 ( α + β ) B 1 ( α + β ) B 2 ( α + β ) B 3 ( α + β ) ] = [ B 0 ( α ) B 0 ( β ) B 1 ( α ) B 1 ( β ) B 2 ( α ) B 2 ( β ) B 3 ( α ) B 3 ( β ) ] (2-1-13) \begin{bmatrix} {B_0(\alpha+\beta)}\\ {B_1(\alpha+\beta)}\\ {B_2(\alpha+\beta)}\\ {B_3(\alpha+\beta)}\\ \end{bmatrix}= \begin{bmatrix} {B_0(\alpha)}{B_0(\beta)}\\ {B_1(\alpha)}{B_1(\beta)}\\ {B_2(\alpha)}{B_2(\beta)}\\ {B_3(\alpha)}{B_3(\beta)}\\ \end{bmatrix} \tag{2-1-13} B0(α+β)B1(α+β)B2(α+β)B3(α+β)=B0(α)B0(β)B1(α)B1(β)B2(α)B2(β)B3(α)B3(β)(2-1-13)

    B n ( l ) = e − n l j π 2 (2-1-14) B_n(l)=e^{-\frac{nlj\pi}{2}}\tag{2-1-14} Bn(l)=e2nljπ(2-1-14)
    ​ 因为 D F T DFT DFT是正交变换,所以式(2-1-12)的矩阵的逆矩阵也是该矩阵的转置,因此求解出序列 A j ( α ) A_j(\alpha) Aj(α)式(2-1-15).将其整理成 式(2-1-16),整理的过程以 A 0 ( α ) A_0(\alpha) A0(α)为例子(见附),其余请读者自行推导.
    A 0 ( α ) = 1 2 ( 1 + e − j π α 2 + e − j π α + e − 3 j π α 2 ) (2-1-15) A_0(\alpha)=\frac{1}{2}(1+e^{-\frac{j\pi\alpha}{2}}+e^{-j\pi\alpha}+e^{-\frac{3j\pi\alpha}{2}}) \tag{2-1-15} A0(α)=21(1+e2jπα+ejπα+e23jπα)(2-1-15)

    A 0 ( α ) = 2 c o s ( π α 4 ) c o s ( 2 π α 4 ) e − 3 j π α 4 (2-1-16) A_0(\alpha)=2cos(\frac{\pi\alpha}4)cos(\frac{2\pi\alpha}{4})e^{-\frac{3j\pi\alpha}4} \tag{2-1-16} A0(α)=2cos(4πα)cos(42πα)e43jπα(2-1-16)
    ​ 2.离散分数傅里叶变换算法

    ​ 为了说清楚离散的分数傅里叶变换算法,先聊一聊 D F T DFT DFT,对一个信号做N点 D F T DFT DFT写成矩阵的形式为式(2-2-1)
    [ Y ( 0 ) Y ( 1 ) ⋮ Y ( N − 1 ) ] 1 N [ e − j 2 π 0 N e − j 2 π 0 N ⋯ e − j 2 π 0 N e − j 2 π 1 N e − j 2 π 1 N ⋯ e − j 2 π 1 N ⋮ ⋮ ⋱ ⋮ e − j 2 π ( N − 1 ) N e − j 2 π ( N − 1 ) N ⋯ e − j 2 π ( N − 1 ) N ] [ X ( 0 ) X ( 1 ) ⋮ X ( N − 1 ) ] (2-2-1) \begin{bmatrix} {Y(0)}\\ {Y(1)}\\ {\vdots}\\ {Y(N-1)}\\ \end{bmatrix} \frac{1}{\sqrt[]N} \begin{bmatrix} {e^{-\frac{j2\pi0}{N}}}&{e^{-\frac{j2\pi0}{N}}}&{\cdots}&{e^{-\frac{j2\pi0}{N}}}\\ {e^{-\frac{j2\pi1}{N}}}&{e^{-\frac{j2\pi1}{N}}}&{\cdots}&{e^{-\frac{j2\pi1}{N}}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {e^{-\frac{j2\pi(N-1)}{N}}}&{e^{-\frac{j2\pi(N-1)}{N}}}&{\cdots}&{e^{-\frac{j2\pi(N-1)}{N}}}\\ \end{bmatrix} \begin{bmatrix} {X(0)}\\ {X(1)}\\ {\vdots}\\ {X(N-1)}\\ \end{bmatrix} \tag{2-2-1} Y(0)Y(1)Y(N1)N 1eNj2π0eNj2π1eNj2π(N1)eNj2π0eNj2π1eNj2π(N1)eNj2π0eNj2π1eNj2π(N1)X(0)X(1)X(N1)(2-2-1)
    ​ 既然在2.1中我们发现了傅里叶变换的’周期’性那么这种’周期’性是否在 D F T DFT DFT中也有表现呢?对离散信号做一次 D F T DFT DFT可 以表达为 Y = E X Y=EX Y=EX,做两次 D F T DFT DFT相当于 Y = E 2 X Y=E^2X Y=E2X,而 E 2 E^2 E2等于式(2-2-2)
    [ 1 0 ⋯ 0 0 0 ⋯ 1 ⋮ ⋮ ⋱ ⋮ 1 0 ⋯ 0 ] (2-2-2) \begin{bmatrix} {{1}}&{{0}}&{\cdots}&{{0}}\\ {{0}}&{{0}}&{\cdots}&{{1}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {{1}}&{{0}}&{\cdots}&{{0}}\\ \end{bmatrix} \tag{2-2-2} 101000010(2-2-2)
    ​ 所以连续做两次 D F T DFT DFT等价于对时域进行了翻转.以此类推,可以发现, D F T DFT DFT也具有’周期’性根据离散傅里叶变换的矩阵 形式,定义记号(2-2-3)他们的递推关系为(2-2-4)则定义阶数为a的离散分数傅里叶变换为(2-2-5).
    X ( l ) = ( x 0 ( l ) , x 1 ( l ) , . . . , x N − 1 ( l ) ) T (2-2-3) X^{(l)}=(x_0^{(l)},x_1^{(l)},...,x_{N-1}^{(l)})^{T} \tag{2-2-3} X(l)=(x0(l),x1(l),...,xN1(l))T(2-2-3)

    X ( l + 1 ) = E X ( l ) = E l + 1 X ( 0 ) , l = 0 , 1 , 2 , 3 (2-2-4) X^{(l+1)}=EX^{(l)}=E^{l+1}X^{(0)},l=0,1,2,3 \tag{2-2-4} X(l+1)=EX(l)=El+1X(0),l=0,1,2,3(2-2-4)

    X ( α ) = ∑ i = 0 3 A l ( α ) X ( l ) = ( ∑ i = 0 3 A l ( α ) E l ) X ( 0 ) (2-2-5) X(\alpha)=\sum_{i=0}^3A_l(\alpha)X^{(l)}=(\sum_{i=0}^3A_l(\alpha)E^{l})X^{(0)} \tag{2-2-5} X(α)=i=03Al(α)X(l)=(i=03Al(α)El)X(0)(2-2-5)

    ​ 引入矩阵记号(2-2-6),则阶数为 α \alpha α的离散分数傅里叶变换可以写成(2-2-7).
    E α = ∑ i = 0 3 A l ( α ) E l (2-2-6) E^{\alpha}=\sum_{i=0}^3A_l(\alpha)E^{l} \tag{2-2-6} Eα=i=03Al(α)El(2-2-6)

    X ( α ) = E α X 0 (2-2-7) X^{(\alpha)}=E^{\alpha}X^{0} \tag{2-2-7} X(α)=EαX0(2-2-7)

    ​ 这样,对于任何长度为N的数据序列,利用公式(2-2-8)可以得到其 α \alpha α阶的离散分数傅里叶变换
    X ( α ) = ( ( ∑ i = 0 3 2 c o s ( π ( α − i ) 4 ) c o s ( 2 π ( α − i ) 4 e − 3 j π ( α − i ) 4 ) E l ) X 0 (2-2-8) X^{(\alpha)}=((\sum_{i=0}^32cos(\frac{\pi({\alpha}-i)}4)cos(\frac{2\pi({\alpha}-i)}4e^{-\frac{3j\pi({\alpha}-i)}4})E^{l})X^{0} \tag{2-2-8} X(α)=((i=032cos(4π(αi))cos(42π(αi)e43jπ(αi))El)X0(2-2-8)


    三 离散分数傅里叶变换快速算法

    ​ 我们已经知道了对一个长度为N的数据序列如何求其 α \alpha α阶的离散分数傅立叶变换.这一小节中我们将开始探讨以下是否有一种快速算法能够减少直接计算的运算量.回到定义式(2-2-5),我们将其展开,得到(3-1-1).我们可以发现后两项中可以提取出一个E矩阵,然后得到式(3-1-2).
    X ( α ) = A 0 ( α ) X ( 0 ) + A 1 ( α ) E X ( 0 ) + A 2 ( α ) E 2 X ( 0 ) + A 3 ( α ) E 3 X ( 0 ) (3-1-1) X(\alpha)=A_0(\alpha)X^{(0)}+A_1(\alpha)EX^{(0)}+A_2(\alpha)E^2X^{(0)}+A_3(\alpha)E^3X^{(0)} \tag{3-1-1} X(α)=A0(α)X(0)+A1(α)EX(0)+A2(α)E2X(0)+A3(α)E3X(0)(3-1-1)

    X ( α ) = A 0 ( α ) X ( 0 ) + A 2 ( α ) E 2 X ( 0 ) + E ( A 1 ( α ) X ( 0 ) + A 3 ( α ) E 2 X ( 0 ) ) (3-1-2) X(\alpha)=A_0(\alpha)X^{(0)}+A_2(\alpha)E^2X^{(0)}+E(A_1(\alpha)X^{(0)}+A_3(\alpha)E^2X^{(0)}) \tag{3-1-2} X(α)=A0(α)X(0)+A2(α)E2X(0)+E(A1(α)X(0)+A3(α)E2X(0))(3-1-2)

    后两项与前两项结构大致相同,只不过权重系数不同.我们将算法分成两个部分.第一部分计算基本趋势,即对原始数据做 α \alpha α阶离散分数傅里叶变换所得序列 X ( α ) X^{(\alpha)} X(α)的第一个值, x 0 ( α ) x_0^{(\alpha)} x0(α)它应该是式(3-1-3).然后计算其余值,姑且叫做谐波部分吧.这里以求解序列 X ( α ) X^{(\alpha)} X(α)第一个值为例,其余读者可以自行推导.对于序列 X ( α ) X^{(\alpha)} X(α)的第一个值,可以展开成式(3-1-4),
    x 0 ( α ) = ( A 0 ( α ) + A 1 ( α ) ) x 0 ( 0 ) + 1 N ( A 1 ( α ) + A 3 ( α ) ) ( x 0 ( 0 ) + x 1 ( 0 ) + ⋯ + x N − 1 ( 0 ) ) (3-1-3) x_0^{(\alpha)}=(A_0(\alpha)+A_1(\alpha))x_0^{(0)}+\frac{1}{\sqrt[]N}(A_1(\alpha)+A_3(\alpha))(x_0^{(0)}+x_1^{(0)}+\cdots+x_{N-1}^{(0)}) \tag{3-1-3} x0(α)=(A0(α)+A1(α))x0(0)+N 1(A1(α)+A3(α))(x0(0)+x1(0)++xN1(0))(3-1-3)

    x 1 ( α ) = A 0 ( α ) x 1 ( 0 ) + A 2 ( α ) x N − 1 ( 0 ) + E 2 ( A 1 ( α ) x 1 ( 0 ) + A 3 ( α ) x N − 1 ( 0 ) ) (3-1-4) x_1^{(\alpha)}=A_0(\alpha)x_1^{(0)}+A_2(\alpha)x_{N-1}^{(0)}+E^2(A_1(\alpha)x_1^{(0)}+A_3(\alpha)x_{N-1}^{(0)}) \tag{3-1-4} x1(α)=A0(α)x1(0)+A2(α)xN1(0)+E2(A1(α)x1(0)+A3(α)xN1(0))(3-1-4)

    其中 A 1 ( α ) , A 2 ( α ) , A 3 ( α ) , A 4 ( α ) A_1(\alpha),A_2(\alpha),A_3(\alpha),A_4(\alpha) A1(α),A2(α),A3(α),A4(α)可以化简成式(3-1-5).可以发现,他们都含有共同项().将其提取出来.并且定义新矩阵 K ( α ) = L ( α ) c o s ( 2 π α 4 ) e − 3 j π α 4 K(\alpha)=L(\alpha)cos(\frac{2\pi\alpha}{4})e^{-\frac{3j\pi\alpha}4} K(α)=L(α)cos(42πα)e43jπα,其中 L ( α ) L(\alpha) L(α)(式3-1-6)为.于是得到第二部分的计算公式(式3-1-7).读者可以比较一下该算法的运算次数并与 F F T FFT FFT进行对比.
    { A 0 ( α ) = c o s ( 2 π α 4 ) c o s ( π α 4 ) e − 3 j π α 4 A 2 ( α ) = j c o s ( 2 π α 4 ) s i n ( π α 4 ) e − 3 j π α 4 A 1 ( α ) = c o s ( 2 π ( α − 1 ) 4 ) c o s ( π ( α − 1 ) 4 ) e − 3 j π ( α − 1 ) 4 A 3 ( α ) = j c o s ( 2 π ( α − 1 ) 4 ) s i n ( π ( α − 1 ) 4 ) e − 3 j π ( α − 1 ) 4 (3-1-5) \begin{cases} A_0(\alpha)=cos(\frac{2\pi\alpha}{4})cos(\frac{\pi\alpha}{4})e^{-\frac{3j\pi\alpha}4}\\ A_2(\alpha)=jcos(\frac{2\pi\alpha}{4})sin(\frac{\pi\alpha}{4})e^{-\frac{3j\pi\alpha}4}\\ A_1(\alpha)=cos(\frac{2\pi({\alpha}-1)}4)cos(\frac{\pi({\alpha}-1)}4)e^{-\frac{3j\pi({\alpha}-1)}4}\\ A_3(\alpha)=jcos(\frac{2\pi({\alpha}-1)}4)sin(\frac{\pi({\alpha}-1)}4)e^{-\frac{3j\pi({\alpha}-1)}4}\\ \end{cases} \tag{3-1-5} A0(α)=cos(42πα)cos(4πα)e43jπαA2(α)=jcos(42πα)sin(4πα)e43jπαA1(α)=cos(42π(α1))cos(4π(α1))e43jπ(α1)A3(α)=jcos(42π(α1))sin(4π(α1))e43jπ(α1)(3-1-5)

    [ c o s ( π α 4 ) 0 ⋯ 0 j s i n ( π α 4 ) 0 c o s ( π α 4 ) ⋯ j s i n ( π α 4 ) 0 ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 ⋯ c o s ( π α 4 ) 0 j s i n ( π α 4 ) 0 ⋯ 0 c o s ( π α 4 ) ] (3-1-6) \begin{bmatrix} {{cos(\frac{\pi\alpha}{4})}}&{{0}}&{\cdots}&{{0}}&{{jsin(\frac{\pi\alpha}{4})}}\\ {{0}}&{{cos(\frac{\pi\alpha}{4})}}&{\cdots}&{{jsin(\frac{\pi\alpha}{4})}}&{{0}}\\ {\vdots}&{\vdots}&{\ddots}&{{\vdots}}&{\vdots}\\ {{0}}&{{0}}&{\cdots}&{{{cos(\frac{\pi\alpha}{4})}}}&{0}\\ {{jsin(\frac{\pi\alpha}{4})}}&{{0}}&{\cdots}&{{0}}&{{cos(\frac{\pi\alpha}{4})}}\\ \end{bmatrix} \tag{3-1-6} cos(4πα)00jsin(4πα)0cos(4πα)000jsin(4πα)cos(4πα)0jsin(4πα)00cos(4πα)(3-1-6)

    ( x 1 ( l ) , x 2 ( l ) , . . . , x N − 1 ( l ) ) T = K ( α ) ( x 1 ( 0 ) , x 2 ( 0 ) , . . . , x N − 1 ( 0 ) ) T + K ( α − 1 ) ( x 1 ( 1 ) , x 2 ( 1 ) , . . . , x N − 1 ( 1 ) ) T ( x 1 ( 1 ) , x 2 ( 1 ) , . . . , x N − 1 ( 1 ) ) T = E ( x 1 ( 0 ) , x 2 ( 0 ) , . . . , x N − 1 ( 0 ) ) T (3-1-7) (x_1^{(l)},x_2^{(l)},...,x_{N-1}^{(l)})^{T}=K(\alpha)(x_1^{(0)},x_2^{(0)},...,x_{N-1}^{(0)})^{T}+K({\alpha}-1)(x_1^{(1)},x_2^{(1)},...,x_{N-1}^{(1)})^{T} \tag{3-1-7}\\ (x_1^{(1)},x_2^{(1)},...,x_{N-1}^{(1)})^{T}=E(x_1^{(0)},x_2^{(0)},...,x_{N-1}^{(0)})^{T} (x1(l),x2(l),...,xN1(l))T=K(α)(x1(0),x2(0),...,xN1(0))T+K(α1)(x1(1),x2(1),...,xN1(1))T(x1(1),x2(1),...,xN1(1))T=E(x1(0),x2(0),...,xN1(0))T(3-1-7)


    四 总结

    ​ 相比于传统的傅里叶分析,分数傅里叶变换具有先天的时频联合分析的能力.这一点是傅里叶分析所达不到的(具体原因我会在后续与大家一起讨论).像 D F T DFT DFT那样,离散分数傅里叶变换也具有其快速算法,目前在通信系统中的应用也是比较多的,当然目前比较新颖的方向包括了小波分数傅里叶分析,这里并不是说先对信号做分数傅里叶变换,然后对所得结果再做一次的分数傅里叶变换.而是利用 M A L L A T MALLAT MALLAT分解算法以及合成算法的正交矩阵去改变分数傅里叶变换的特征向量,使得能够满足人们对特定场合的需求.(若有错误,请指出,我会虚心请教)


    五 附


    A 0 ( α ) = 1 2 e − j π ( α ) 4 ( e j π ( α ) 4 + e − j π ( α ) 4 ) + 1 2 e − 5 j π ( α ) 4 ( e j π ( α ) 4 + e − j π ( α ) 4 ) = c o s ( j π α 4 ) e − j π ( α ) 4 + c o s ( j π α 4 ) e − 5 j π ( α ) 4 = c o s ( j π α 4 ) e − 3 j π ( α ) 4 ( e 2 j π ( α ) 4 + e − 2 j π ( α ) 4 ) = 2 c o s ( j π α 4 ) c o s ( 2 j π α 4 ) e − 3 j π ( α ) 4 A_0(\alpha)=\frac{1}{2}e^{-\frac{j\pi({\alpha})}4}(e^{\frac{j\pi({\alpha})}4}+e^{-\frac{j\pi({\alpha})}4})+ \frac{1}{2}e^{-\frac{5j\pi({\alpha})}4}(e^{\frac{j\pi({\alpha})}4}+e^{-\frac{j\pi({\alpha})}4})\\ =cos(\frac{j\pi{\alpha}}4)e^{-\frac{j\pi({\alpha})}4}+cos(\frac{j\pi{\alpha}}4)e^{-\frac{5j\pi({\alpha})}4}=cos(\frac{j\pi{\alpha}}4)e^{-\frac{3j\pi({\alpha})}4}\\(e^{\frac{2j\pi({\alpha})}4}+e^{-\frac{2j\pi({\alpha})}4})=2cos(\frac{j\pi{\alpha}}4)cos(\frac{2j\pi{\alpha}}4)e^{-\frac{3j\pi({\alpha})}4} A0(α)=21e4jπ(α)(e4jπ(α)+e4jπ(α))+21e45jπ(α)(e4jπ(α)+e4jπ(α))=cos(4jπα)e4jπ(α)+cos(4jπα)e45jπ(α)=cos(4jπα)e43jπ(α)(e42jπ(α)+e42jπ(α))=2cos(4jπα)cos(42jπα)e43jπ(α)

    展开全文
  • 分数阶傅里叶变换在信号检测与图像处理中的应用研究(李琼) 发展历史 特性 二维分数阶傅里叶变换 基于分数阶傅里叶变换的图像分析 分数阶变换域中图像的能量分布 分数阶傅里叶变换域中的图像的幅度和相位信息 分数...
  • 线性调频信号的分数阶傅里叶变换,完整可用。分数阶傅里叶变换采用Pei采样算法,在傅里叶域呈现宽带特性的线性调频信号,在分数阶域呈现窄带特征。
  • 基于离散多参数分数阶傅里叶变换和混沌映射的双图像加密
  • python中的离散1d和2d分数阶傅里叶变换 用法: frft2d(mat,ax,ay) mat:要转换的数字矩阵 ax,ay:沿x和y轴的变换顺序 返回垫子的frft光谱 偏移(f,a; p) f:待变换离散信号 a:转换顺序 p(可选):近似的...
  • 基于分数阶傅里叶变换统计信息的人脸识别(孙慧静) 文章地址:paper 分数阶傅里叶变换作为新兴的图像处理工具,可以看作信号在时频平面进行任意角度的逆时针旋转。它引入了变换阶次,可以分析信号的时频域信息...
  •  %频域变换最优阶数,从频域中分数阶傅里叶变换和从时域中变换不一样,注意最优阶中的Fs,Kr for i=1:Na xx_f(i,:)=frft(xx_f(i,:),opr); end G2=abs(xx_f); figure plot(G0) grid on figure plot(G1) % 显示原始...
  • 顾名思义,该信号的频率随着时间线性变换,其复数表达形式如下: s(t)=e2jπ(f0t+0.5μt2)s(t)=e^{2j\pi(f_0 t+ 0.5\mu t^2)}s(t)=e2jπ(f0​t+0.5μt2) 根据欧拉公式,其相位项为2π(f0t+0.5μt2))2\pi(f_0 t+ 0.5\...
  • 西 南 交 通 大 学 毕 业 论 文 基于离散采样型的分数阶FOURIER变换算法 研究与实现 年 级: 2011 学 号: 20112627 姓 名: 方威 专 业: 自动化(交通信息工程及控制方向) 指导老师: 汪晓宁 二零一五年六月 院 系 专 业...
  • 分数阶傅里叶变换与图像加密技术,辛怡,陶然,网络多媒体技术的飞速发展使图像内容的安全问题日益严峻。图像加密已成为保护图像内容安全的两种有效手段之一。分数阶傅里叶变换

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 284
精华内容 113
关键字:

离散分数阶傅里叶变换