分数阶傅里叶变换的离散算法-Ozaktas

离散分数阶Fourier变换(DFRFT)算法 FRFT
这篇文献发表于：
作者：
• 一、分数阶Fourier变换的定义
• 二、分数阶与其他时频分析工具（ Wigner-
Ville分布）的关系
• 三、离散分数阶傅立叶变换的计算
一、分数阶Fourier变换的定义
二、分数阶傅里叶变换与Wigner-Ville分布
1.第一种分解方法
可以把以上改写为
X e x(t)e e dt 1 jcot
j
u2
cot 2


p
2

j
t
2
cot 2

jtu csc
X e x(t)e e dt 1 jcot
j
u2
cot 2


p
2

j
t
2
cot 2

jtu csc
假定p∈[-1,1],经过量纲归一化的信号x(t）的分数阶傅里叶变换， 可以分解为以下三个步骤：
（1）用chirp信号调制信号f(x):
式1
（2）调制信号与另一个chirirp信号调制卷积后的信号：
式3
1
tan

1
sin

tan(

2
)
• 具体细节：第一步：将函数，f (x) 与线性调频函 数相乘(式1)。
• 注意，g(x)的频率带宽与时间带宽乘积可以是， f(x)的相应带宽乘积的两倍，所以要求g(x)的采 样间隔为1／(2Δ x)。如果，( )样本值的采样间 隔是1／Δ x，那么就需要对这些样本值进行插 值，然后再与线性调频函数的离散采样值相乘， 以得到所希望的g(x)的采样。
• 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); • function z = fconv(x,y) • % convolution by fft • N = length([x(:);y(:)])-1; • P = 2^nextpow2(N); • z = ifft( fft(x,P) .* fft(y,P)); • z = z(1:N);
是B频率截断的的，这时，只要采样间隔 / B
f (t) 按间隔 进行采样就不会损失信息，而且，
利用采样序列 f (n); n Z 可按如下公式构造原信号
f (t) f (n) sin (t n) /
nZ
(t n) /
上式Shannon 插值公式。
也就是说 Xp() 的Wigner-Ville分布， 是由 X (t) 的Wigner-Ville分布旋转а 角得到。
• 所以分数阶Fourier变换有一个重要的 性质，分数阶Fourier变换是角度为α的 时频面旋转. 这个性质建立起分数阶 Fourier变换与时频分布间的直接联系, 并且为分数阶Fourier域理解为一种统 一的时频变换域奠定了理论基础, 同时 也为分数阶Fourier变换在信号处理领 域中的应用提供了有利条件。
ω v
α α
u t
三、离散分数阶傅立叶变换的计算
• 目前DFRFT的四种离散化算法
3
1.利用FFT计离散FRFT的变换核 F p ai ( p)W i j0
2.根据FRFT的定义，将其分解为信号的卷积，再 利用FFT计算。 3.利用矩阵的特征值和特征向量计算。 4.直接离散化计算。
在这篇文献中，第二种，采用分解的方法。
• % chirp premultiplication • chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2); • f = chrp.*f; • % chirp convolution • c = 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 multiplication • Faf = chrp.*Faf; • % normalizing constant • Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);
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);
• 首先，看一下Wigner-Ville分布

WX(t, w)
X
(w

2
)X
*
(w

2
)e
j
t
d

X (w) 是 x(t) 傅里叶变换
经过一系列变换后变为

Wx(t, w) 2e2 juv X ( )X*(2u )e2 jv d
由以上可得，等式的右边是 Xp( ) 的Wigner-Ville分布， 左边是 X (t) 的Wigner-Ville分布
• % 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 = tan(alpha/2); • sina = sin(alpha); • f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];
x( )e n 2x
( 2x )2
sin
c[2x(t

n 2x
)]
n N
简要介绍一下Shannon 插值
Shannon定理 到设信号 f (t) L2 (R) ，如果存在 B 0 ，使
fˆ ( ) 0 ， B ， a.e. R
则称 f (t) L2 (R)
为了简化计算，人们提出更加有效的分解计算方法。
假定x(t)的wigner-ville分布限定在以原点为中心，直径为Δ x的
圆内。若令
，则与chirp信号乘积后的信号
x(t )e j cot 在频域具有带宽Δx。可以用Shannon插值表示
N
j (cot ) n2
x(t)e j cot
• 第二步：将g(x)与一线性调频函数作卷积式(式
• (2))。注意，由于g(x)是带限信号，所以线性调 频函数也可以用其带限形式代替而不会有任何 影响。
2、第二种分解方法
Xp
e x(t)e e dt 1 j cot 2
j
u2
cot 2



j
t
2
cot 2

jtu csc