分数阶余弦变换

分数阶余弦变换（共7篇）

分数阶余弦变换 篇1

0引言

随着世界通信技术的发展, 既需要在有效利用频带传输多种信息, 更需要保证通信安全, 为了保证信息安全, 有混沌加密, RSA加密等多种加密方式。

而基于分数阶Hilbert变换的单边带通信, 更是将传统Hilbert的解析信号抑制负频谱推广到了分数阶Fourier域。这种方法最大的优点是, 在传输信号的同时可以利用分数阶Hilbert变换阶数作为加密密钥, 保证通信安全。但是加密信号的信号只有单加密密钥, 即密钥空间较小[1]。

基于多抽样速率滤波器可以在频谱分析中将信号分为几个互不重叠的子带, 分别进行编码, 压缩等处理, 可以起到节省计算工作量及存储空间的目的。

本文利用多抽样速率滤波器组, 提出了一种可以将信号进行多路传输, 各路独立加密的对称密钥加密方案。

1分数阶Hilbert变换原理

1.1 分数阶Fourier变换

分数阶Fourier变换的定义为[2]:

式中:undefined为整数, p为分数阶Fourier的变换阶数。

虽然实信号的分数阶Fourier变换不满足共轭对称性, 但是任意信号f (t) 的分数阶Fourier变换可以看成是信号undefined的Fourier变换。如果r (t) 为实信号, 可以用指数信号undefined调制, 调制后为undefined, 调制后信号的分数阶Fourier变换为:

式中:Fπ/2表示传统Fourier变换 (看做是旋转角度为π/2的分数阶Fourier变换) , Aα与u无关, undefined关于u为偶对称, 而实信号r (t) 的Fourier变换满足共轭对称性。因此保留f (t) 的分数阶Fourier正谱Fα[f] (u) (u≥0) , 就能得到实信号r (t) 的正频谱, 从而恢复实信号[3]。

1.2 分数阶Hilbert变换

类似于传统Hilbert变换, 定义信号x (t) 的分数阶Hilbert变换为[4]:

式中:undefined;undefined为传统Hilbert变换的核函数, 即h (t) =1/πt (t≠0) , α为分数阶Fourier变换的旋转角度。

相应的解析信号:

式中:undefined为信号x (t) 的分数阶Hilbert变换。此解析信号undefined只包含信号x (t) 的分数阶Fourier正谱分量[5]。

设信号x (t) 是实信号r (t) 经过指数信号调制后得到的信号, 即:undefined, 则解析信号undefined和r (t) 的关系为:

1.3 基于分数阶Hilbert变换的单边带通信

信号发射前必须经过调制, 根据式 (3) , 式 (4) 生成解析信号, 此解析信号只包含原信号的分数阶Fourier域的正频谱部分[6]。如图1所示。

接收端接收到解析信号后, 利用式 (5) 恢复, 然后解调即可恢复x (t) 。如图2所示。

设正确角度为α, 若接收端不按照α为角度进行恢复, 假设以角度β解调, 可将β代入式 (5) 进行整理, 则生成的实信号r′ (t) 为:

由此可以看出, 信号不但节省了频带, 同时还保持了一个加密密钥, 接收端只有满足角度α=β时, 才能精确重构原信号。

2基于分数阶Hilbert变换的多速率信号处理多路加密方法及性能分析

2.1 基于分数阶Hilbert变换的加密原理

多路分数阶Hilbert变换加密通信框图如图3所示。

在图3中, 待发送信号首先使用多抽样速率滤波器组分为N路。N通道滤波器组如图4[7]所示。

实际应用中, 为了高效计算, 常采用滤波器组的多相结构。设E (z) 和R (z) 分别是H (z) 和G (z) 的多相分解形式, 定义为:

式中:Q=M/N (M为滤波器长度) ;k=0, 1, 2, …, N-1;hl为H0～HN-1中第l个滤波器。

式中:undefined为滤波器长度) ;k=0, 1, 2, …, N-1;gl为G0～GN-1中第l个滤波器[8]。

若满足:

则系统为完全重建系统。式中:c和λ为常数;I为单位矩阵。

经过多抽样速率滤波器组分为N路后, 分别进行各路加密, 各路使用式 (3) , 式 (4) 加密, 分别设定独立的旋转角度。

在接收端, 各路根据式 (5) , 按照双方约定的密钥进行解密。

2.2 基于多抽样速率滤波器组的多加密密钥系统的实现

综上所述, 提出一种可以将单路信号使用多个加密密钥进行加密的系统, 如图5所示。

图5中, x0 (n) ～xN-1 (n) 为图4中原信号x (n) 经过N通道滤波器组经过分析滤波器H0 (z) ～HN-1 (z) 分成不同的子带信号。使用式 (3) , 式 (4) 加密后进行传输。接收端利用式 (5) , 对加密信号undefined进行解密得到xα0 (n) ～xαN-1 (n) , 将这些信号进行插值恢复原采样率和经过综合滤波器组处理, 最后归并各路信号可以恢复原信号。恢复信号和原信号的关系为:

式中:c和n0为常数。

由图5可以看出, 加密密钥由一个拓展为N个, 即密钥空间增大了。

2.3 性能分析

单路加密算法的计算复杂度分别是复杂度为O (n) 的乘法, O (nlog n) 的卷积运算。当n足够大时, 后者支配前者。所以系统总的计算复杂度为O (nlog n) 。其中, n为数据长度[9]。

多加密密钥系统中, 应采用高效算法, 即使用多相结构, 并且将分析滤波器组的多相分量与抽取器、插值器和综合滤波器组的多相分量分别进行等价变换。这样可以使得卷积和乘法运算在低抽样率的一端, 使得卷积和乘法的计算复杂度减少。设采用N个加密算子, 即多抽样速率滤波器组共N路。

根据以上讨论。则原信号x (n) 经过抽取后分别与分析滤波器组的多相分量进行卷积运算, 每路卷积运算的计算复杂度为undefined, N路总的计算复杂度为undefined;加密和解密的计算复杂度占支配地位的是卷积运算。进行各路加密, 此时各路数据长度为n/N, 每路加密和传输后解密的复杂度为undefined, N路总的加密解密计算复杂度为undefined;同理, 最后和综合滤波器组的各多相分量进行卷积运算后插值恢复, 总的计算复杂度也为undefined。所以, 系统总的计算复杂度为undefined。因此, 与单路加密系统相比, 多加密系统不但使得系统总的计算复杂度降低, 并且增大了密钥空间[10]。

3仿真实验

仿真实验系统为2通道分数阶Hilbert变换加密系统, 输入信号为方波, 分析滤波器组的多相分量为, 综合滤波器组的多相分量为, 分别以π/6和π/5作为第1路和第2路的加密密钥, 经传输后, 分别以正确密钥 (π/6和π/5) 和错误密钥 (π/6.5和π/5) 进行解密。仿真实验结果如图6所示。

由图6 (b) 可见, 当密钥正确时, 可以重建原信号;若某一路解密, 设第一路角度变为π/6.5, 第二路保持不变, 则经过综合滤波器后两路信号和合并重建的信号如图6 (c) 所示。因此, 只有知道N个通道独有的加密密钥, 才能解密重建原信号。

4结语

根据分数阶Hilbert变换理论, 可以得到分数阶域信号的解析表示, 并且可以使用分数阶相位参数作为单边带通信的加密密钥。基于多通道滤波器组的多路加密系统不但使计算复杂度更小, 同时还可以将单路信号的加密密钥拓展为N个, 从而增大了密钥空间, 因此通信的安全性更高。该系统可用于存在多发射端, 多接收端的系统, 有一定的实用价值。

参考文献

[1]李雪梅, 陶然, 李炳照.基于FRFT的Hilbert变换[J].兵工学报, 2009, 30 (2) :221-225.

[2]陶然, 邓兵, 王越.分数阶傅里叶变换及其应用[M].北京:清华大学出版社, 2009.

[3]ZAYED A I.On the relationship between the Fourier andfractional Fourier transform[J].IEEE Signal ProcessingLetters, 1996, 3 (12) :310-311.

[4]ZAYED A.A convolution and product theorem for the frac-tional Fourier transform[J].IEEE Signal ProcessingLetters, 1998, 5 (4) :101-103.

[5]TAO Ran, LI Xue-Mei, WANG Yue.Generalization of thefractional Hilbert transform[J].IEEE Signal ProcessingLetters, 2008, 15:365-368.

[6]ZAYED A I.Hilbert transform associated with the fraction-al Fourier transform[J].IEEE Signal Processing Letters, 1998, 5 (8) :206-208.

[7]陶然, 张惠云, 王越.多抽样率数字信号处理理论及其应用[M].北京:清华大学出版社, 2007.

[8]胡广书.现代信号处理教程[M].北京:清华大学出版社, 2004.

[9]DASGUPTA Sanjoy.算法概论[M].王沛, 译.北京:清华大学出版社, 2008.

[10]宗孔德.多抽样率数字信号处理[M].北京:清华大学出版社, 1996.

分数阶余弦变换 篇2

分数阶微分算子在分数阶控制系统的求解计算中具有重要作用, 其固有特性决定了理论上的分数阶控制器实际上是无穷维线性滤波器, 这就要求在分数阶控制器的实现过程中使用有限带宽的整数阶控制器对其进行近似[4]。而在计算机控制系统中, 分数阶微分算子的离散化是分数阶控制器数字化实现的关键[5]。一般地, 有两种离散化方法:一种是直接离散化方法, 如Euler算子的直接幂级数展开法 (power series expansion, PSE) [6]、Tustin算子的连分式展开法 (continued fractional expansion, CFE) [7]和Muir递归展开法 (Muir-recursion) [8]等;另一种是间接离散化方法, 一般是先进行连续时间频域拟合, 然后将拟合得到的s函数离散化, 如Oustaloup算法[9]、Charef算法[10]等。

Tustin变换是连续系统离散化的常用方法。研究了基于Tustin变换的分数阶微分算子直接离散化方法并比较了基于PSE、CFE和Muir递归展开法的近似离散化方法。

1 分数阶微积分基础

1.1 分数阶微积分的定义

分数阶微积分是经典整数阶微积分的自然扩展, 它允许微积分的阶次为任意实数甚至可以为复数。常用的两种分数阶微积分定义是Grünwald-Letniko定义和Riemann-Liouville定义[11]。对连续可导函数f (t) 的r阶Grünwald-Letnikov微积分定义如下

式 (1) 中, [·]表示整数部分, 为分数阶微积分算子, r为正表示微分, r为负表示积分, a和t分别为积分上下限。

Riemann-Liouville微积分定义如下

式 (2) 中, n-1<r<n, Г[·]为伽马函数。在实际的物理系统和工程应用中, 上述两个微积分定义是等效的[12]。

拉普拉斯变换常用于解决用传统微分方程描述的工程问题, 同样可以用于分数阶微积分。根据Riemann-Liouville微积分定义, 考虑零初始条件分数阶微积分的拉普拉斯变换为

式 (3) 中, 0<r<1, 表示拉普拉斯变换, s为拉普拉斯算子。

1.2 分数阶动态系统的离散化

分数阶动态系统可用如下分数阶微分方程描述[13]

式 (4) 中, , ak (k=0, …, n) 和bk (k=0, …, m) 为常数, αk (k=0, …, n) 和βk (k=0, …, m) 为任意实数。不失一般性, 可假设αn>αn-1>…>α0及βm>βm-1>…>β0。考虑式 (3) , 则分数阶动态系统的传递函数为

为进一步得到分数阶系统的离散化模型, 首先要对分数阶微积分算子s±r进行离散化近似, 首先将拉普拉斯算子用所谓的生成函数s=ω (z-1) 来表示, 于是可得到动态系统离散时间传递函数的一般表达式为

式 (6) 中, 生成函数ω (z-1) 通常表达成复变量z或移位算子z-1的函数。此生成函数及其展开式决定了分数阶微积分算子逼近的形式和系数。例如, 当使用后向差分进行离散化时, ω (z-1) = (1-z-1) /T, 其中T为采样周期。通过幂级数展开可以得到离散化的分数阶微分算子为FIR形式的数字滤波器[14]。

2 分数阶微分算子的离散化

分数阶微积分算子的离散化方法对分数阶控制系统的具体实现具有十分重要的影响。本文采用Tustin变换为生成函数, 并应用不同的展开方法对式 (7) 进行展开得到分数阶微分算子的离散化近似。

2.1 基于PSE的离散化

分数阶微分算子离散化最简单的方法是对Tustin算子直接幂级数展开。通过幂级数展开可以得到其离散时间传递函数, 这实际上是一个FIR形式的数字滤波器。对其进行适当的截断处理[6], 则可以得到近似的离散时间传递函数, 其表达式如下

式 (8) 中, PSE{u}表示函数u的幂级数展开, n为近似阶数。文献[15]给出了一种Tustin算子直接幂级数展开的IIR形式数字滤波器算法, 其表达式如式 (9) 。

式 (9) 中, Pp和Qq分别是阶数为p和q的多项式。

2.2 基于CFE的离散化

基于PSE的分数阶微分算子离散化得到的是多项式函数, 而有理函数在计算或插值中往往要优于多项式。连分式是对函数或数值的一种有效近似形式, 它比幂级数展开的收敛速度快, 并且可以应用到复数空间[7]。因此, 应用连分式展开来进行分数阶微分算子的近似, 能得到更好的近似效果。一般地, 函数G (z) 可用连分式形式表示为

式 (10) 中, ai (z) 和bi (z) 为z的有理函数或常量。于是通过截断处理, 可以得到近似的有理函数G^ (z) 为

式 (11) 中, CFE{u}表示对函数u的连分式展开, Pp和Qq分别是阶数为p和q的多项式。

2.3 基于Muir-recursion的离散化

如何得到分数阶微分算子离散化的递推算法是分数阶控制数字实现的关键。文献[8]介绍了Muir递推算法, 该算法最初应用于石油探测的数据处理。目前, 常与Tustin变换结合用于分数阶微分算子的离散化。不失一般性, 假设r∈[0, 1], 分数阶微分算子的递推公式为

式 (12) 中, A0 (z-1, r) =1, 且

式 (13) 中, n为奇数时cn=r/n, n为偶数时cn=0。对于已给定的近似阶数n, An (z-1, r) 和An (z-1, -r) 可用MATLAB符号工具箱得到, 从而分数阶微分算子sr的近似离散化传递函数可由式 (14) 计算。

3 仿真与比较

为比较以上三种基于Tustin变换的分数阶微分算子离散化方法 (以下分别简记为PSE+Tustin、CFE+tustin、Mur+Tustin) , 取10以内奇数近似阶次, 采样周期T=0.001 s对分数阶微分算子s0.5进行离散化近似。以下图1~图3分别为PSE+Tustin、CFE+Tustin和Mur+Tustin方法的伯德图。其中, 图例Gf表示s0.5的连续域伯德图, Gp3、Gc3、Gr3分别代表近似阶次为3的s0.5离散化数字滤波器的伯德图, 其他图例含义以此类推。

由伯德图可见, 就以上三种离散化近似方法而言, 在低频段幅频特性和相频特性的误差都较大, 而在高频段则较小, 且频率越低误差越大。就某一种方法而言, 离散化近似的幅频特性结果要好于其相频特性, 且近似阶数越高近似效果越好。然而, 需要特别注意的是, 在频率很高时, 各种方法的近似效果又都急剧变差。

相比较而言, CFE+Tustin方法在较宽频带内对幅频特性和相频特性的近似效果都要好于另外两种方法, 也就是说采用较低阶次可以得到另外两种方法取较高阶次时的近似结果, 这就易于数字滤波器的设计和数字化控制器的实现。但从仿真计算过程来看, 当进一步提高近似阶次以期进一步提高近似精度时, 将耗费过多的计算时间。

本文中的PSE+Tustin方法和Mur+Tustin方法都是递推方法, 为比较增加近似阶次时近似误差的变化, 定义误差指数如公式 (15) 。

式 (15) 中, Bd (jω) 和Bc (jω) 分别表示数字滤波器和连续微分算子的幅频特性或相频特性;ωu和ωl分别为所考查频率的上下界。图4示出了这一变化趋势。可见, 随着近似阶数的提高, 这两种方法的近似误差都快速减小, 并在阶数较高时趋于稳定。但从实际的仿真计算过程来看, PSE+Tustin算法比Mur+Tustin算法的运算速度更快, 在近似阶数超过25阶以后差别尤其显著。

4 结束语

分数阶微分算子的离散化是分数阶控制器数字实现的关键。论文研究了基于Tustin变换的分数阶微分算子直接离散化方法, 并对PSE+Tustin、CFE+Tustin和Mur+Tustin三种方法进行了比较。

分数阶余弦变换 篇3

传统的傅里叶变换(FFT)对平稳信号的处理效果很好,但当信号频率随时间变化时,FFT就显得有些力不从心了。分数阶傅里叶变换(Fractional Fourier Transform,FRFT)可以很好地弥补FFT的不足,特别是处理线性调频信号(LFM)时,能够得到令人满意的结果。

FRFT也称为角度傅里叶变换(AFT) 或者旋转傅里叶变换(RFT),其定义式为:

$X{}_p(u)=R{}^{\alpha }[x(t)]={\displaystyle {\int }_{-\infty }^{\infty }{\rm K}}{}_p(t,u)x(t)\text{d}t(1)$

式中变换核取作:

${\rm K}{}_p(t,u)=\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}(\frac{1}{2}u{}^2\cot \alpha -ut\text{c}\text{s}\text{c}\alpha +\frac{1}{2}t{}^2\cot \alpha )},\alpha \ne n\text{π}\\ \delta (t-u),\alpha =2n\text{π}\\ \delta (t+u),\alpha =(2n+1)\text{π}\end{array}$

则:

$\begin{array}{l}X{}_p(u)=\{F{}^{D}x\}(u)={\displaystyle {\int }_{-\infty }^{\infty }x}(t)k{}_p(t,u)\text{d}t\\ =\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }{\displaystyle {\int }_{-\infty }^{\infty }x}(t)\text{e}{}^{\text{j}\frac{t{}^2}{2}\cos \alpha }\text{e}{}^{-\text{j}ut\text{c}\text{s}\text{c}\alpha }\text{d}t,\alpha \ne n\text{π}\\ x(t),\alpha =2n\text{π}\\ x(-t),\alpha =(2n+1)\text{π}\end{array}(2)\end{array}$

其中n为整数,即n∈Z。α=pπ/2称为分数阶Fourier变换的阶数,并有$R{}^{\alpha }=R{}^{p\frac{\pi }{2}}=F{}^p$。Kp(t,u)称为FRFT的核函数。Xp(u)称为x(t)的p阶Fourier变换。FRFT是一种线性算子,记为Fp,他满足以下性质:

(1) FRFT变换为线性算子;

(2) F0[x(t)]=F4[x(t)]=x(t)(恒等变换);

(3) F1[x(t)]=F5[x(t)]=X(ω)(标准Fourier变换);

(4) 广义Fourier变换算子为加性算子,即有Fp+q=FpFq。

2 采用分解方法计算FRFT的步骤

分数阶Fourier变换可以具体分解为以下三个主要的计算步骤:线性调频信号乘法;线性调频信号卷积;另一个线性调频信号乘法。假定p∈[-1,1],则我们可以将经过量纲归一化的信号f(x)的分数阶Fourier变换式(2)分解为以下三步运算:

$f{}_p(x)=\text{e}{}^{-\text{j}\text{π}x{}^2\tan (\alpha /2)}{g}^{\prime }(x)(3)$

和:

$\begin{array}{l}{g}^{\prime }(x)=A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{\text{j}\text{π}\beta (x-x^{\prime }){}^2}g(x^{\prime })\text{d}{x}^{\prime }(4)\\ g(x)=\text{e}{}^{-\text{j}\text{π}x{}^2\tan (\alpha /2)}f(x)(5)\end{array}$

式中,α=pπ/2,β=csc α,而g(x)和g′(x)表示中间结果,并且:

$A{}_{\alpha }=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}(6)$

即是说,分数阶Fourier变换的数值计算的顺序如下:先计算式(式(5)),再计算式(4),最后计算式(3)。下面是每一步计算的有关细节。

第一步:将函数f(x)与线性调频函数相乘(式(5))。注意,g(x)的频率带宽与时间带宽乘积可以是f(x)的相应带宽乘积的两倍,所以要求g(x)的采样间隔为1/(2Δx)。如果f(x)样本值的采样间隔是1/Δx,那么就需要对这些样本值进行插值,然后再与线性调频函数的离散采样值相乘,以得到所希望的g(x)的采样。

第二步:将g(x)与一线性调频函数作卷积式(式(4))。注意,由于g(x)是带限信号,所以线性调频函数也可以用其带限形式代替而不会有任何影响。也就是说,我们可以取:

$\begin{array}{l}{g}^{\prime }(x)=A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{\text{j}\text{π}\beta (x-x^{\prime }){}^2}g(x^{\prime })\text{d}{x}^{\prime }\\ =A{}_{\alpha }{\displaystyle {\int }_{-\infty }^{\infty }h}(x-x^{\prime })g(x^{\prime })\text{d}{x}^{\prime }(7)\end{array}$

式中:

$h(x)={\displaystyle {\int }_{-\text{Δ}x}^{\text{Δ}x}{\rm H}}(\gamma )\text{e}{}^{\text{j}2\text{π}\gamma x}\text{d}\gamma (8)$

这里:

${\rm H}(\gamma )=\frac{1}{\sqrt{\beta }}\text{e}{}^{\text{j}\text{π}/4}\text{e}{}^{-\text{j}\text{π}\gamma {}^2/\beta }(9)$

是函数ejπβx2的Fourier变换。于是,式(7)的离散形式为:

$g^{\prime }(\frac{m}{2\Delta x})={\displaystyle \sum _{n=-{\rm N}}^{{\rm N}}h}\left(\frac{m-n}{2\text{Δ}x}\right)g\left(\frac{n}{2\text{Δ}x}\right)(10)$

这一离散卷积可以利用快速Fourier变换计算。

第三步:计算式(3)得到f(x)的分数阶Fourier变换fp(x)的采样值$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$。由于假定f(x)的所有变换都是带限的,他们位于区间$\left[-\frac{1}{2}\text{Δ}x,\frac{1}{2}\text{Δ}x\right]$,所以需要用因子2对$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$进行二抽一采样,以得到离散采样$f{}_p\left(\frac{m}{2\text{Δ}x}\right)$。因此,对于不是π/2整数倍的角度,分数阶Fourier变换的计算对应以下步骤:

(1) 原信号与一线性调频函数相乘;

(2) Fourier变换(其变元乘以尺度系数csc α);

(3) 再与一线性调频函数相乘;

(4) 乘以一复幅值因子。

3 对信号的FRFT处理及仿真图

首先要给出一个输入信号x,然后根据分解方法编出FRFT快速算法,根据不同的p值输入信号x会生成不同的曲线,其中p∈(0,4)。找出每个p值时FRFT的输出最大值点,组成一个一维数组m1,再从m1中找出一个最大值,该值所对应的p值就是FRFT变换的最佳角度α=pπ/2。

例1 假设输入信号为方波C=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ],长度N=73,在p=0～1之间进行FRFT变换,仿真结果如图1所示。

由图1(d)中可以看出当p=1时,α=π/2就是普通的傅里叶变换,这也验证了分数阶傅里叶变换的正确性。

例2 输入信号为多项式的线性调频信号fmpoly(N,p),x=exp(j*pi*(n-center)·^p/(N/2)^(p-1)/p),其中N为信号长度,p是信号的多项式幂,当p=2时该信号为chirp信号。取x=fmpoly(63,2)。

由图2可以看出,在图2(e)中p=1.5的时候,即α=3π/4时形成了一个冲击信号,说明了在此角度上信号的能量最好地集聚在一点上,由此可以识别出信号的调频系数,检测出信号的参数,这就是FRFT处理LFM信号的显著作用。

4 结 语

分数阶傅里叶变换是近二十年来发展起来的一种全新的信号时频分析工具,在很多方面得到了十分广泛的应用。而其快速算法的研究则对扩展其应用领域有着十分重要的意义。本文提出了一种有效并能准确计算FRFT 的新算法。该算法具有易实现、易理解、精度较高等优点,相信FRFT将会受到更广泛的重视,在信号处理领域会有良好的应用前景。

参考文献

[1]平先军,陶然,周思永,等.一种新的分数阶傅里叶变换快速算法[J].电子学报,2001,29(3):406-408.

[2]Lufs B Almeida.The Fractonal Fourier Transform andTime-Frequency Representations[J].IEEE Transactionson Signal Processing,1994,42(11).

分数阶余弦变换 篇4

Fourier变换 (Fourier Transform, FT) 作为最主要的信号分析工具主要用于处理频率不随时间变化的平稳信号, 在时频平面时间轴与频率轴相互垂直, 即Fourier变换是从时间域旋转π/2到频率域。但是Fourier变换通常无法表述信号的时频域性质, 不能表示某种频率分量发生在哪个时间, 而这种性质恰恰是非平稳信号最关键的性质, 这对非平稳信号十分重要。分数阶Fourier变换 (Fractional Fourier Transform, FrFT) 是Fourier变换的一种推广形式, 揭示了信号从时间域到频率域变化过程中所呈现的特征, 即从时间域和频率域同时表示信号旋转π/2分数倍时的特征, 从而克服了传统Fourier变换不能反映非平稳信号的统计量随时间变化的缺陷[1,2,3,4]。

分数阶Fourier变换是一种线性变换, 与经典的Fourier变换有着天然的联系, 又提供了Fourier变换所不具备的某些特点, 而且它与小波变换、Wigner-Ville分布 (WVD) 都有密切的关系[3], 因此近年来分数阶Fourier变换受到了研究者的广泛关注, 相应的离散变换算法也相继提出[4]。此外, 有关分数阶Fourier变换应用的研究也越来越多, 如[5,6,7,8,9,10]等文献就介绍了分数阶Fourier变换在信号处理和通信系统方面的应用。

线性调频 (Chirp) 信号是一种典型的非平稳信号, 它的瞬时频率随时间呈线性变化, 常见于雷达、声纳和移动通信等系统中。对Chirp信号的研究, 是非平稳随机信号处理理论及方法的基础。研究证明, 分数阶Fourier变换就是一种很适合处理Chirp信号的变换。近年来, 将分数阶Fourier变换用于线性调频信号 (包括时不变、时变幅度线性调频信号) 的检测和参数估计引起了越来越多的关注[11,12]。本文介绍了FrFT的定义、性质和简单应用, 分析了高斯脉冲Chirp信号的FrFT以及最优化分析角度的选取问题。

1 分数阶Fourier变换的基础

分析和处理平稳信号最常用和最主要的方法是Fourier变换 (FT) 。Fourier变换建立了信号整个时域与整个频域的对应关系, 其中:

$\begin{array}{l}X(\omega )=1/\sqrt{2\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }x}(t)\text{e}{}^{-\text{j}\omega t}\text{d}t(1)\\ x(t)=1/\sqrt{2\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }X}(\omega )\text{e}{}^{\text{j}\omega t}\text{d}t(2)\end{array}$

时域和频域构成了分析一个信号的两种表达方式, Fourier变换在整体时域上将信号分解为不同的频率分量, 但它没有将时域和频域组合成一个域, 不能提供时间和频率的联合信息。谱X (ω) 只是显示任一频率ω包含在信号x (t) 内的总强度, 无法表述信号的时频域性质, 不能表示某种频率分量发生在哪个时间, 而这种性质恰恰是非平稳信号中的最根本和最关键的性质。如果x (t) 是由几个非平稳分量组成的, 那么时间上的任何变化都会改变X (ω) 。此时, 传统的Fourier变换就不能满足信号分析的要求了。

时频分析的基本思想就是设计时间和频率的联合函数, 用它描述信号在不同时间和频率上的能量密度或强度, 从而克服传统Fourier变换不能反映非平稳信号的统计量随时间变化的缺陷。分数阶Fourier变换 (FrFT) 是借用时频面的概念, 以时间和频率分别为横轴和纵轴, 旋转一定的角度进行的一种线性变换。

传统的Fourier变换X (ω) 就是x (t) 旋转π/2, 即x (t) 由时间轴t变到频率轴ω的表示形式。令:

$\alpha =p\text{π}/2(3)$

并且定义线性算子:

$R{}^{\alpha }=R{}^{p\text{π}/2}(4)$

记作Fp。F2相当于t轴连续两次逆时针旋转π/2, 得到x (-t) ;F3相当于t轴连续三次逆时针旋转π/2, 得到指向-ω轴的函数;F4表示t轴连续四次逆时针旋转π/2, 得到原函数。因此线性算子Rα有以下数学性质:

(1) 零旋转:

$R{}^0=x(t)(5)$

(2) 与Fourier变换等价:

$R{}^{\text{π}/2}=F[x(t)](6)$

(3) 旋转相加性:

$R{}^{\alpha }R{}^{\beta }=R{}^{\alpha +\beta }(7)$

(4) 2π旋转 (恒等变换) :

$R{}^{2\text{π}}=x(t)(8)$

如果角度α以π/2的非整数倍进行旋转, 则得到函数x (t) 的广义Fourier变换, 记作:

$\{R{}^{\alpha }x\}(u)=\{F{}^{p}x\}(u)=X{}_p(u)(9)$

这就是x (t) 的分数阶Fourier变换:

Xp (u) ={Fpx} (u) =∫∞-∞x (t) Kp (t, u) dt (10)

式中:变换核Kp (t, u) 定义为:

${\rm K}{}_p(t,u)=\{\begin{array}{l}\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{t{}^2+u{}^2}{2}\cot \alpha -\text{j}ut\text{c}\text{s}\text{c}\alpha },\alpha \ne n\text{π}\\ \delta (t-u),\alpha =2n\text{π}\\ \delta (t+u),\alpha =(2n+1)\text{π}\end{array}(11)$

式中:n为整数。

从以上的定义可以看出, 分数阶Fourier变换是经典Fourier变换的广义形式, 它包含了信号的时间域和频率域表示。旋转角度为π/2时, 即阶数为1的分数阶Fourier变换就是传统的Fourier变换;不旋转或旋转的角度为2π的整数倍则为信号本身;当旋转角度不在以上两个位置即p为分数时, 它同时从时间域和频率域给出了信号的特征。在讨论分数阶Fourier变换时, 由于角度以2π为模, 所以只需要考虑$0\le \left|p\right|\le 2$的旋转阶数。α角度的分数阶Fourier反变换对应着-α角度的分数阶Fourier变换, 即:

$x(t)={\displaystyle {\int }_{-\infty }^{\infty }X}{}_p(u){\rm K}{}_{-p}(t,u)\text{d}u(12)$

Fourier变换在线性系统分析、光学系统、信息处理系统等方面起着核心作用, 并应用于众多的工程技术领域。作为广义Fourier变换, FrFT具有比Fourier变换更普遍的特性和更广的应用场合, 尤其是普通Fourier变换技术不能解决问题的场合, FrFT更显其优越性。在FrFT研究的启示下, 许多学者推广了分数阶的概念, 得到分数阶卷积和分数阶相关, 还把分数阶的概念应用到Hadamard变换、Hartley变换等, 得到相应的分数阶变换。而且, 分数阶Fourier变换已应用到图象处理的优化图象恢复方面。此外, 将分数阶Fourier变换算子的分数幂推广至复数幂也成为目前研究的一个热点。

在20世纪90年代中期, 分数阶Fourier变换被引入了信号处理领域。在信号处理中, FrFT有很多应用, 其中两个典型的应用是信号滤波和信号分离[5]。实验已经证实分数阶Fourier变换滤波的效果明显优于Fourier变换;特别是对于线性调频Chirp信号, 分数阶Fourier变换能够获得最佳的能量积聚效果, 是分数阶Fourier变换最合适的应用领域之一。

2 高斯脉冲Chirp信号的分数阶Fourier变换

线性调频Chirp信号是一种特殊的非平稳信号, 它的瞬时频率随时间呈线性变化, 广泛地出现在通信、雷达、声纳和地震勘探等系统中。在工程实践中, 高斯调制Chirp信号有着广泛的应用, 而FrFT对于Chirp信号良好的检测和分析效果自然让我们想到能否将此分析用具引入到高斯调制的Chirp信号分析中。本文对此进行了深入研究, 给出了相关解析结果。

高斯信号的标准形式为:

$f(t)=\frac{1}{\sqrt{2\text{π}\delta {}^2}}\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}(13)$

这是一个服从正态分布 (t0, δ) 的函数, 它的Fourier变换和分数阶Fourier变换都具有高斯信号的形式[12]。利用已有的这些性质, 可以研究高斯脉冲Chirp信号的分数阶Fourier变换。

令高斯脉冲Chirp信号为:

$x(t)=\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)}(14)$

将高斯脉冲Chirp信号写成幅度和相位的函数, 得到:

$x(t)=\left|x(t)\right|\text{e}{}^{\text{j}\theta (\omega )}(15)$

其中:

$\begin{array}{l}\left|x(t)\right|=\text{a}\text{b}\text{s}(\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)})=\text{e}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\\ \theta (\omega )=at{}^2+bt+c\end{array}(16)$

式 (14) 所表示的高斯脉冲Chirp信号的分数阶Fourier变换为:

$\begin{array}{l}X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }{\displaystyle \int \text{e}}{}^{-\frac{(t-t{}_0){}^2}{2\delta {}^2}}\text{e}{}^{\text{j}(at{}^2+bt+c)}\text{e}{}^{\text{j}\frac{t{}^2}{2}\cot \alpha }\text{e}{}^{-\text{j}ut\text{c}\text{s}\text{c}\alpha }\text{d}t\\ =\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}{\displaystyle \int \text{e}}{}^{-[(\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2})t{}^2-(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha )t]}\text{d}t\\ ={\rm P}{\displaystyle \int \text{e}}{}^{-({\rm M}{}^{2}t{}^2-2{\rm M}{\rm N}t)}\text{d}t(17)\end{array}$

其中:

$\begin{array}{l}{\rm P}=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2\text{π}}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\\ {\rm M}=\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}\\ {\rm N}=\frac{\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha }{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\end{array}(18)$

所以式 (17) 化为:

$\begin{array}{l}X{}_p(u)={\rm P}{\displaystyle \int \text{e}}{}^{-({\rm M}{}^{2}t{}^2-2{\rm M}{\rm N}t+{\rm N}{}^2-{\rm N}{}^2)}\text{d}t\\ ={\rm P}\text{e}{}^{{\rm N}{}^2}{\displaystyle \int \text{e}}{}^{-({\rm M}t-{\rm N}){}^2}\text{d}t(19)\end{array}$

令Mt-N=Q, 则$t=\frac{Q+{\rm N}}{{\rm M}},\text{d}t=\frac{\text{d}Q}{{\rm M}}$, 因此得到:

$X{}_p(u)=\frac{{\rm P}\text{e}{}^{{\rm N}{}^2}}{{\rm M}}{\displaystyle \int \text{e}}{}^{-Q{}^2}\text{d}Q(20)$

然后取积分限为 (-∞, ∞) , 由于${\displaystyle {\int }_{-\infty }^{\infty }\text{e}}{}^{-Q{}^2}\text{d}Q=\sqrt{\text{π}}$, 所以可得:

$X{}_p(u)=\frac{\sqrt{\text{π}}{\rm P}\text{e}{}^{{\rm N}{}^2}}{{\rm M}}(21)$

再将P, M, N的值代入式 (21) 得:

$\begin{array}{l}X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\cdot \\ \frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{\left(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha \right){}^2}{\frac{2}{\delta {}^2}-\text{j}4a-\text{j}2\cot \alpha }}(22)\end{array}$

将$\frac{(t{}_0/\delta {}^2+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha ){}^2}{2/\delta {}^2-\text{j}4a-\text{j}2\cot \alpha }$的分子分母同时乘以2/δ2+j4a+j2cot α, 则式 (22) 化为:

$X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{-\frac{t{}_0^2}{2\delta {}^2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\cdot$

$\text{e}{}^{\frac{\left(\frac{t{}_0}{\delta {}^2}+\text{j}b-\text{j}u\text{c}\text{s}\text{c}\alpha \right){}^2\left(\frac{2}{\delta {}^2}+\text{j}4a+\text{j}2\cot \alpha \right)}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}(23)$

若t0=0, 可以得到:

$X{}_p(u)=\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\cdot$

$\text{e}{}^{\text{j}c}\text{e}{}^{\text{j}\frac{u{}^2}{2}\cot \alpha }\text{e}{}^{\frac{-\text{j}(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}\text{e}{}^{\frac{-\frac{2}{\delta {}^2}(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}}(24)$

将Chirp信号的分数阶Fourier变换写成幅度和相位的函数, 得到:

$X{}_p(u)=\left|X{}_p(u)\right|\text{e}{}^{\text{j}\rho (\omega )}(25)$

因为$\left|\text{e}{}^{\text{j}A}\right|=1$, 在计算Xp (u) 的幅值时可以忽略所有模为1的部分, 这样就能得到:

$\left|X{}_p(u)\right|=$

$\text{a}\text{b}\text{s}\left(\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{-(u-b\text{s}\text{i}\text{n}\alpha ){}^2}{2\sin {}^2\alpha [\frac{1}{\delta {}^2}+\delta {}^2(2a+\cot \alpha ){}^2]}}\right)(26)$

令:$\sigma {}^2=\sin {}^2\alpha [\frac{1}{\delta {}^2}+\delta {}^2(2a+\cot \alpha ){}^2](27)$

得:

$\begin{array}{l}\left|X{}_p(u)\right|=\text{a}\text{b}\text{s}[\sqrt{\frac{1-\text{j}\text{c}\text{o}\text{t}\alpha }{2}}\\ \frac{1}{\sqrt{\frac{1}{2\delta {}^2}-\text{j}a-\text{j}\frac{\cot \alpha }{2}}}\text{e}{}^{\frac{-(u-b\text{s}\text{i}\text{n}\alpha ){}^2}{2\sigma {}^2}}](28)\end{array}$

由式 (24) 和式 (25) 得:

$\text{e}{}^{\text{j}\rho (\omega )}=\text{e}{}^{\text{j}\left[c+\frac{u{}^2}{2}\cot \alpha -\frac{(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}\right]}(29)$

因此Xp (u) 的相位:

$\rho (\omega )=c+\frac{u{}^2}{2}\cot \alpha -\frac{(4a+2\cot \alpha )(u\text{c}\text{s}\text{c}\alpha -b){}^2}{\frac{4}{\delta {}^4}+(4a+2\cot \alpha ){}^2}$ (30)

由式 (28) 可得$\left|X{}_p(u)\right|$具有高斯函数的形式, 服从正态分布 (bsin α, σ2) 。当方差σ2最小时, 高斯脉冲Chirp信号的分数阶Fourier变换Xp (u) 的能量最大, 因此这时的α就是最优化角度αopt。从式 (27) 可知, σ2最小时, 有2a+cot αopt=0, 即:

$\alpha {}_{\text{o}\text{p}\text{t}}=-\arctan \frac{1}{2a}(31)$

同样, 由式 (3) 可知, 此时Xp (u) 的最优化阶数popt为:

$p{}_{\text{o}\text{p}\text{t}}=-\frac{2}{\text{π}}\arctan \frac{1}{2a}(32)$

由式 (31) 和式 (32) 对比文献[12]可得, Chirp信号和高斯脉冲调制Chirp信号的αopt, popt是一样的, 这说明加上高斯窗的Chirp信号 (即高斯脉冲Chirp信号) 的最优化参数并没有改变。这就为高斯Chirp信号的FrFT最优化分析奠定了坚实的理论基础, 可以充分利用FrFT对于LFM信号最优的能量积聚性质, 将其应用于高斯Chirp信号的分析中, 在低信噪比的环境中这一点更具有实际意义。

3 结 语

分数阶Fourier变换是Fourier变换的广义形式。作为一种新的时频分析工具, 分数阶Fourier变换既与经典的Fourier变换有着天然的联系, 又提供了Fourier变换所不具备的某些特点, 而且它与小波变换、Wigner-Ville分布 (WVD) 都有密切的关系[12,13], 是一种很适合处理线性调频Chirp信号等非平稳信号的变换。

在Chirp信号的检测和参数估计方面, 分数阶Fourier变换具有重要的应用价值。尤其是它所具有的线性特性, 可以方便地用于多分量线性调频信号的检测和参数估计而不受交叉项的干扰, 大量被用于恒定幅度的线性调频信号, 并取得了良好的效果。事实上, 时变幅度线性调频信号在工程实际中有着更为广泛的应用背景, 在工程实践中, 会遇到很多经过幅度调制的线性调频信号, 其中的高斯窗调制就是典型的一种。

从分析结果看, 经过高斯窗调制的Chirp信号保持了原本的特性, 在FrFT这个分析工具下能够获得优异的性能。

分数阶余弦变换 篇5

有关STLFMCW信号检测的问题, 相关文献有诸多论述。其中, 文献[6]总结了Wigner-ville分布、Choi-Williams分布、正交镜像滤波器组、循环平稳分析、Wigner Hough等方法, 提出了一种FRFT与聚类分析相结合的方法。该方法对于单分量信号处理效果较为理想。当信号分量增多时, FRFT检测器的检测能力下降导致聚类很难实施。文献[7]提出了基于极大chirplet变换 ( MCT) 的FMCW信号检测方法, 可以用于处理STLFMCW信号; 但其处理过程较为繁琐、限制条件较多, 不太适于实际工程应用。文献[8]提出了周期Wigner Hough变换 ( PWHT) 的方法, 通过设计LFMCW匹配函数, 实现在非匹配接收条件下, 对LFMCW信号检测的最优处理; 但该方法计算量很大, 目前没有快速算法。同时, 该方法处理STLFMCW信号时, 由于信号自身函数发生变化, 已经不是最优处理。

针对STLFMCW信号的时频分布结构, 提出了FRFT循环处理 ( CFRFT ) 方法。通过分析STLFM-CW信号在FRFT循环域 ( CFRFD) 的分布特征, 推导STLFMCW信号在CFRFD的尖峰值、坐标位置和峰值处的信噪比公式。CFRFT对STLFMCW信号具有比FRFT更强的检测能力, 具有与周期WignerHough变换处理LFMCW信号时类似的处理效果, 计算流程简单, 可实现低信噪比条件下STLFMCW信号的检测。

1 分数阶傅里叶变换

分数傅里叶变换算子Fa通过实变量a将函数x变换为Xa= Fa ( x) , 定义可以表述为整体积分变换

式 ( 1) 中, a为FRFT的阶数, u为分数阶域 ( FR-FD) , 旋转角 α = a2π, 则FRFT变换核Ka ( u, t) 为

式 ( 2) 中,

Xa ( u) 的逆变换为

由式 ( 3) 可知, 信号x ( t) 可以分解为一组系数为Xa ( u) 的正交Chirp基K－ a ( u, t) 的线性组合。随着变换阶数a从0 连续增长到1, 展示出了信号从时域逐步变化到频域的所有特征。

分数阶傅里叶变换具有一些基本性质, 例如

(1) 线性变换

( 2) Parseval关系 ( 能量守恒定律)

(3) 旋转相加性

2 STLFMCW信号模型

STLFMCW信号的每个调制周期包含绝对值相等的正、负调频率两个LFM信号, 可表示为

式中A为幅度, fc为载频, B为调制带宽, T为调制周期, T = 2tm。本文以2 个调制周期的STLFMCW信号为例进行分析, 其时频分布如图1 所示。

在高斯白噪声背景下, 信号模型可以表示为

式中, w ( t) 是均值为0, 方差为 σ2w的复高斯白噪声, 信号的输入信噪比为SNRin=A2/ σ2w。

3 STLFMCW信号FRFT分析

国内外学者提出了很多种快速近似分数阶傅里叶变换算法［9—11］。实际应用中, 需要处理的是一组原始连续信号经采样后得到的离散观测数据。采用文献[11]提出的二相型算法, 精度高、计算速度更快。量纲归一化处理时, 采用文献[12]中的离散尺度化法。信号的FRFT可以理解为, 信号在以角在时频平面逆时针旋转后的分数阶傅里叶域 ( FRFD) 上的投影［13］。

单分量LFM信号s ( t) = s1 ( t) 在最佳旋转角的FRFD具有最好的能量聚集性。根据FRFT的性质2, 信号s ( t) 的能量在FRFD域坐标 ( a1, u1) 处聚集成尖峰, 如图2、图3 所示。参数设置: 信号载频fc= 10 MHz, 线性调频带宽B = 10MHz, 调制周期T = 8 μs, 采样频率fs= 80 MHz。

由文献[14]可知, 信号s ( t) 在FRFD域坐标 ( a1, u1) 处尖峰值为

式 ( 11) 中,

将式 ( 12) 代入式 ( 11) 中, 由于N ~ 1, 则有

S为归一化因子,

检测器的输出信噪比决定信号检测的性能。检测器在尖峰值处的信噪比定义为

当时, 则输出信噪比, 仅比匹配接收机的输出信噪比SNR=NSNRin低3 d B。4

STLFMCW信号, 根据FRFT性质1,

信号分量s1 ( t) 、s2 ( t) 、s3 ( t) 、s4 ( t) 在FRFD域形成4 个等幅度尖峰坐标分别为 ( a1, u1) 、 ( a2, u2) 、 ( a3, u3) 、 ( a4, u4) , 如图4、图5 所示。尖峰坐标存在以下关系:

可见, STLFMCW信号在阶数为频域附近全部重叠, 存在交叉项, 形成伪峰值, 降低了FRFT的检测性能。

4 STLFMCW信号FRFT循环处理

对于任意能量有限信号x ( t) , 0 ≤ t ≤ Td, 其FRFT循环处理 ( CFRFT) 定义式为

式 ( 22) 中, 为循环处理的周期,

假设STLFMCW信号调制周期数为, 根据FRFT性质2, STLFMCW信号在FRFD满足能量守恒。则易证: 当时, 信号的CFRFT在FRFT循环域 ( CFRFD) 坐标 ( as, us) 处形成尖峰, 即取得最大值

因此, 在CFRFD进行二维峰值搜索即可完成STLFMCW信号的检测与参数估计。由于高斯白噪声的相关函数为冲击函数, 其功率谱函数为常数, 导致高斯白噪声在FRFD接近平坦分布。因此, 信号的CFRFT各分量在做循环处理时, 其方差可以认为是线性叠加。则有, 信号尖峰 ( a1, u1) 坐标处的输出信噪比为

当时, 则输出信噪比, 仅比匹配接收机输出信噪比SNR=PNSNRin低3 d B, 可以看出, 随着STLFMCW信号分量的增加, 其FRFT循环处理检测器可以实现信号在FRFT循环域的准脉冲积累, 依然可以保持与FR-FT检测器对于LFM信号相类似的检测能力。STLFMCW信号的CFRFT如图6、图7所示。

可见, 通过FRFT循环处理, STLFMCW信号在CFRFD形成尖峰, 其检测过程与FRFT检测器相类似。

5 仿真验证

仿真参数条件设置: STLFMCW信号包含2 个调制周期, 调制周期T = 8 μs, 信号载频fc= 10MHz, 线性调频带宽B = 10 MHz, 采样频率fs= 80MHz。利用Monte Carlo方法进行1 000 次计算机仿真实验, 比较FRFT检测器与CFRFT检测器的性能, 获得仿真实验结果。

接收机特性曲线 ( ROC) 可以直观反映接收机的性能。两种检测器的ROC曲线如图8 所示。横轴Pfa代表虚警概率, 纵轴Pd代表检测概率。通过图8 可以看出, SNRin= - 15 d B时信号的CFRFT检测器ROC曲线明显优于SNRin= - 12 d B时信号的FRFT检测器。SNRin= - 18 d B时CFRFT检测器ROC曲线略差于SNRin= - 12 d B时信号的FRFT检测器, 此时, 两种检测器的性能都变得较差, 很难实现信号的有效检测。仿真结果与理论推导相接近。

6 结束语

分数阶余弦变换 篇6

关键词：分数阶Fourier变换,多普勒系数,水声通信

0 引言

多普勒效应是移动通信中产生信号畸变的一个重要原因,尤其是在水声信道中,由于声波在水中传播速度远低于电磁波在空气中的传播速度,水声通信中的多普勒效应要远大于空气中无线电通信的多普勒效应。水声通信中的多普勒效应会带来高误码率的出现,甚至是完全不能实现通信任务。分数阶傅里叶变换(Fractional Fourier Transform,简称FRFT)是一种广义形式的傅里叶变换,可以理解为时频平面上的旋转算子。信号变换到分数域上同时包含了时域与频域的信息,使FRFT具有许多传统傅里叶变换不具备的性质,可用来处理类似于LFM信号的非平稳信号。[1]本文对水声通信中多普勒效应的估计进行研究,利用分数阶Fourier变换估计经过信道后的LFM信号调频斜率的变化量值,进而得到多普勒系数的估计。

1 分数阶Fourier变换原理

分数阶Fourier变换有着若干不同的定义[2],不同定义有着不同的物理解释,但是它们之间是相互等价的.本文采用如下定义:

定义在时域的信号f(t),其FRFT为:

式(1)中,p为FRFT的阶数,FP为FRFT算子,Kp(u,t)为FRFT的核函数,定义如下:

式(2)中,α=pπ/2为旋转角度,又称旋转因子;

分数阶Fourier变换具有几条十分重要的性质,在此给出两个需要用到的重要性质:

(1)线性:

(2)阶数可加性(旋转可加性):

对于阶数可加性特殊地有:当P1=-P2=P时,式(4)可化为:

由式(5)可以得出:具有角度α=pπ/2的FRFT变换的逆变换,即为具有角度α=-pπ/2的FRFT变换:

由式(6)所示的变换可以看出,信号f(t)的FRFT变换与Fourier变换原理是相似的,只是函数展开的基底不同。FRFT变换将f(t)表示成由具有线性频率调制的复指数函数集合K-p(t,u)组成的一组基函数,即FRFT的基底是u域上的一组正交的LFM基。因此对于某调频斜率的LFM信号,存在一个分数阶数使LFM信号的能量聚集于一最大值,称之为与此调频斜率相匹配的“最佳”分数阶数。LFM信号调频斜率与“最佳”分数阶数有确定的对应关系[2]:

2 基于FRFT的多普勒系数估计原理

要利用分数阶Fourier变换对线性调频脉冲信号进行处理得到多普勒系数的估计,首先要分析多普勒效应对线性调频脉冲信号的影响。

线性调频脉冲信号是雷达和声纳中常用的信号形式,其复域表达式为:

式(8)中,f0为中心频率,k=B/T为调频斜率,B为信号带宽。在发射端其瞬时相位和频率为:

假设存在多普勒效应,在接收端收到线性调频信号的瞬时相位与频率为[3]:

式(10)中,D=1+△,△为多普勒系数。

对比式(9)和式(10)可得,多普勒效应带给LFM脉冲信号以下三方面的影响:

(1)LFM脉冲信号的时域展宽或压缩;

(2)中心频率由f0变为Df0;

(3)信号调频斜率由k变为D2k。

从第(2)、(3)方面的影响可以看出,只要得到接收LFM脉冲信号中心频率或调频斜率的估计,即可求得多普勒系数的估计。本文通过对调频斜率的估计得到多普勒系数的估计。具体实现方法如下:通过改变分数阶次p,改变旋转因子α=pπ/2,求接收LFM信号的FRFT,得到对应于新调频斜率的最佳阶次。在进行旋转搜索时,可以在原LFM信号对应最佳阶次p附近进行搜索,以减小计算量。由式(7)可得:

另外有接收到LFM信号的调频斜率,代入式(11)中有:

此处因水中载体运动速度的限制,有|△|远小于1,1+△恒为正数,式(12)可化为:

式(13)即为基于分数阶Fourier变换多普勒系数的估计的解析式。

3 仿真研究

本节旨在研究基于FRFT的多普勒系数估计性能,利用水声通信中经常作为同步信号的线性调频信号进行多普勒估计,以期为水声通信系统的设计提供验证性支持。

仿真参数:采样频率fs=48k Hz,LFM信号带宽4~8k Hz,脉宽40ms。

假设水中声速为1500m/s,一般水面舰艇的运动速度为40km左右即±20m/s范围内。为方便分析研究,假设接收船锚定,则有多普勒系数△在±20/1500范围内,根据式(9)、(10)可得接收到的LFM信号调频斜率变化范围为:k'∈[0.9735k 1.0268k],按照LFM信号参数可得其调频斜率为k=B/T。为了后续进行离散化的FRFT计算,需要对k进行尺度归一化[4]有:k=B/fs。接收LFM信号调频斜率变化范围对应的理论最佳阶次范围为:p'∈[1.05151.0543]。因此在变换域进行最佳阶次搜索时直接在此范围内搜索即可:每次搜索从1.0515开始以一定步长向1.0543搜索,求取最佳阶次。

如图一所示,FRFT对给定的LFM信号在变换域存在一聚焦峰,并且有该聚焦峰与LFM信号的最佳阶次相对应。

表一是基于FRFT的多普勒系数估计结果,本节通过Matlab仿真研究了不同多普勒、不同信噪比条件下基于FRFT的多普勒系数估计。采用的噪声为带限高斯白噪声,测得的多普勒系数为多途中最大途径对应的多普勒系数。

图二为基于FRFT的多普勒系数估计偏差以及标准差。由表一、图二分析可得:随着信噪比的增加,多普勒系数的估计结果逐渐逼近真值,基于FRFT多普勒系数的估计偏差与标准差在逐渐减小,在高信噪比下该方法是近似无偏的估计方法,能有效地估计水声通信中的多普勒效应。另外基于FRFT的多普勒系数的估计性能不会随着相对运动速度的增大或者减小而显著变化,在多普勒系数较大且多变的水声通信信道中具有一定的稳健性。

4 结束语

本文首先简要介绍了分数阶傅里叶变换的原理,推导出利用LFM信号的分数阶傅里叶变换进行多普勒系数估计的方法,对该方法进行了仿真研究。通过理论推导与仿真研究得出:在多普勒系数较大且多变的水声通信信道中,基于FRFT的多普勒系数估计方法具有一定的稳健性与有效性。

参考文献

[1]H.M.Ozaktas,Z.Zalevsky,M.A.Kutay.The Fractional Fourier Transform with Application in Optics and Signal Processing.John Wiley&Sons Ltd.2001.

[2][4]陶然,齐林,王越.分数阶Fourier变换原理及应用[M].北京:清华大学出版社,2004.

分数阶余弦变换 篇7

无刷直流电机 (BLDCM) 因其结构简单、运行效率高、调速性能好, 已经越来越广泛地应用在各个领域中。无刷直流电机驱动系统中逆变器的功率半导体管器件是最易发生故障的薄弱环节, 其故障率约占整个系统故障的82.5%[1]。因此, 准确定位逆变器故障, 据此进行容错控制, 是提高无刷直流电机驱动系统运行安全性和可靠性的根本。

逆变器功率管故障可分为短路和开路。功率管短路故障已有成熟的解决方案, 即通过硬件电路检测其D-S压降, 可准确定位故障管[2]。而功率管开路时, 电机往往还能继续运行, 不易发现, 从而导致二次故障, 引发更大的事故。文献[3]对现有的逆变器功率管开路故障诊断方法进行了归类和总结, 可分为基于电流量的诊断方法和基于电压量的诊断方法。文献[4]提出了一种基于傅里叶变换和神经网络的逆变器故障检测与诊断方法, 利用了加窗傅里叶变换提取故障特征。文献[5]对电机三相电流进行小波变换提取故障特征, 并通过BP神经网络识别故障类型。

本文采用分数阶傅里叶变换 (FRFT) 对输出相电流进行故障特征提取, 并与时域、傅里叶变换及小波变换提取的特征进行比较、分析, 寻求最优的故障特征, 最后运用支持向量机进行故障诊断。

1 分数阶傅里叶变换定义及特点

传统的傅里叶变换是分析和处理平稳信号的一种标准和有力的工具, 而对于分析和处理时变的非平稳信号则显得乏力, 这是由于傅里叶变换采用的全局性的基函数所决定的。FRFT方法已经在信号检测、滤波器、图像处理等领域得到广泛的重视和应用[6,7,8,9]。FRFT作为傅里叶变换的一种广义形式, 可以解释为信号在时频平面内坐标轴绕原点逆时针旋转任意角度后构成的分数阶傅里叶域上的表示方法。如图1所示, 其中α=pπ/2, p表示FRFT的分数阶次。

信号的p阶FRFT可表示为:

其中kα (t, u) 为分数阶傅里叶变换的核函数:

常用的核函数有线性核函数、多项式核函数、Sigmoid核函数和径向基核函数。p阶分数阶傅里叶域是在 (t, w) 平面上按逆时针方向旋转的α角度所产生的坐标空间。由于分数阶域的对称性和周期性, p取值一般在0~1之内, 即α取0~π/2。当p=0时, 零旋转对应于信号本身;当p=1时, 信号的FRFT等于经典的傅里叶变换。当p从0变到1时, FRFT能给出信号逐渐从时域变化到频域的所有特性。

分数阶傅里叶变换同小波变换一样, 具有良好的时频局部化特性。可由公式推导得, 将式 (2) 带入式 (1) 得:

其中, 。令y=usecα, 得:

式 (4) 右边的公式与信号x (t) 小波变换为母小波, a为伸缩因子, b为尺度因子) 相比较, 把tan1/2α看成一个伸缩因子, 当tan1/2α增大 (即p增大) 时, 时域窗变宽, 即时域分辨率降低, 而频域窗变窄, 意味着频率分辨率提高, 通过改变阶次p, 可自动调节时域分辨率和频率分辨率。FRFT具有时移、频移和尺度变换等特性, 可应用于信号处理和特征提取。

分数阶傅里叶变化的p值直接影响信号在新的分数阶空间的特征。文献[10]通过穷举法和最后的诊断结果确定最优p值。文献[11]提出一种基于二阶矩极值点的方法确定p值, 具有最佳谱聚集性。文献[12]提出了通过遗传算法和类内类间距离确定最优p值, 获得了较好的诊断结果。

2 实验研究

2.1 仿真实验

采用MATLAB搭建无刷直流电机逆变器模型, 其原理图如图2所示。

T1~T6为功率管, D1~D6为续流二极管, G1~G6为功率管的驱动信号, 由BLDCM的霍尔位置信号判断决定。

仿真参数设置为:电机定子相绕组电阻R=0.6Ω, 定子相绕组自感与互感之差L-M=0.93m H, 转动惯量J=0.00009kg·m2, 额定转速n=3000r/min, 极对数P=4, 48V直流电源供电。无刷直流电机逆变器工作在两两导通方式, 每隔60°电角度换相一次, 每个功率管导通120°。采样频率为5MHz, 逆变器各个功率管分别设置在0.03s时发生故障, 本文只研究单个功率管开路情况。设定电机定子相绕组电阻容差10%, 定子相绕组电感容差5%, 采集A相电流值, 对单个功率管开路及正常工作共7种故障模式各进行50次蒙特卡洛分析, 共采样350组数组。可得到正常工作及单个功率管开路时的A相电流波形如图3所示。

由图3可以看出, T1开路时, A相不再与电源正极相接, A相电流不再为正值;T2开路时, A相电流不再为负值;其它管子开路时, 由于反电动势的影响, A相电流都有明显变化。

2.2 故障特征提取

对于时域分析, 通过采集电流波形信号一个周期的点, 再等间隔采样32个点作为故障特征。对时域采集的特征进行FFT, 取其幅值, 作为频域特征。对于小波分析, 根据电流信号波形和小波函数的相似度, 并且考虑小波的消失矩、正则性、支撑长度等参数, 经比较分析后, 采用了db8小波对电流信号进行6层分解, 将获得的频带能量值作为故障特征[13]。对于分数阶傅里叶域, 以类内类间距离为判据, 对信号进行分数阶傅里叶变换, 寻找最优分数阶傅里叶变换p值, 获得最优故障特征。在欧式距离下可以得到类间离散度Sb和类内离散度Sw公式为:

式中, c为类别数, ni为i类中样本数, Pi是相应类别的先验概率, xk (i) 为i类D维特征向量。mi表示第i类样本集的均值向量, m表示所有各类的样本集总均值向量。一般认为各类类间离散度越大、类内离散度越小, 则类别的可分性越好。所以, 可将Sb/Sw作为可分性判据的标准, 越大说明越有利于分类。当p以一定步长从0变化到1时, 以类内类间距离为判据, 可以找到最优p值。

2.3 基于支持向量机的故障诊断

支持向量机 (SVM) 是由Vapnik和Cortes于1995年提出的。支持向量机是建立在VC维理论和结构风险最小化的基础上的, 因而具有较优的泛化能力, 并巧妙利用内积核函数避免了维数灾难, 通过解一个线性约束的凸二次规划问题得到全局最优解, 因而不存在常规神经网络存在的局部极值等问题[14,15]。支持向量机非常适合解决小样本、非线性及高维模式识别问题。

本文采用了台湾大学林智仁等开发的LIBSVM工具箱进行故障诊断[16]。选用了基于一对一的支持向量机分类器和径向基核函数。径向基函数表达式为:

式中:xi、xj为已知的样本, g为核参数。径向基核函数可以将样本映射到一个更高维的空间, 能够处理当类标签和特征之间的关系是非线性时的样本, 并具有参数少的优点。

将得到的50组数据, 30组作为训练数据, 20组作为测试数据。对数据进行归一化处理, 并通过交叉验证方法得到分类准确率最高惩罚因子c=1和核参数g=16。

3 结果及分析

取p=0:0.001:1, 以类内类间距离为判据, 求得最优p值为0.9130, 得到的故障特征的类内类间距离为0.4473。时域、频域、小波分析及FRFT得到的故障特征的类内类间距离以及诊断结果如表1所示。

由表1可知, 类内类间距离的大小与分类结果基本一一对应, FRFT得到的故障特征有最大的类内类间距离和诊断准确率。小波分析比时域、频域得到的分类结果好, 因为小波分析是一种时频分析方法, 它在时域和频域都具有表征信号局部特征的能力, 特别适宜处理非平稳信号, 而FFT变换只能获得信号在整个频域上的信息, 不适合非平稳信号的分析。小波分解和FRFT都是时频分析方法, 说明在时频域内的特征分析比单纯的在时域或频域内能得到更多的信号细节特征, 有利于提取优质特征数据信息。FRFT与小波分析相比, 通过调整分数阶p值, 将原始信号映射到不同的分数阶域内获得不同的时频分析效果, 通过寻得最优p值, 获得最优的故障特征, 结果表明此方法能提高分类准确率。

4 结语

针对无刷直流电机逆变器功率管开路故障, 本文以类内类间距离为判据, 寻找最优p值, 得到了最优的FRFT故障特征, 并与时域、频域和小波分析得到的故障特征相比较。诊断结果说明此方法可行, 并具有良好的分类性能。

摘要：研究了基于分数阶傅里叶变换 (FRFT) 的无刷直流电机逆变器开路故障诊断技术。首先, 对待测电路的输出电流信号进行分数阶信息分析;其次以类内类间距离为判据, 寻取最优分数阶傅里叶变换p值, 并提取故障特征;最后, 采用基于支持向量机方法进行故障诊断。另外, 对几种常用的信号处理方法得到的结果进行比较, 仿真结果证明, 基于分数阶傅里叶变换的故障诊断技术是可行的。