分数阶傅里叶变换的离散算法-Ozaktas课件
离散分数阶Fourier变换(DFRFT)算法FRFT,这篇文献发表于：,作者：,一、分数阶Fourier变换的定义 二、分数阶与其他时频分析工具（ Wigner-Ville分布）的关系 三、离散分数阶傅立叶变换的计算,一、分数阶Fourier变换的定义,二、分数阶傅里叶变换与Wigner-Ville分布,首先，看一下Wigner-Ville分布,是,傅里叶变换,经过一系列变换后变为,由以上可得，等式的右边是 的Wigner-Ville分布， 左边是 的Wigner-Ville分布,也就是说 的Wigner-Ville分布， 是由 的Wigner-Ville分布旋转角得到。,所以分数阶Fourier变换有一个重要的性质，分数阶Fourier变换是角度为的时频面旋转. 这个性质建立起分数阶Fourier变换与时频分布间的直接联系, 并且为分数阶Fourier域理解为一种统一的时频变换域奠定了理论基础, 同时也为分数阶Fourier变换在信号处理领域中的应用提供了有利条件。,三、离散分数阶傅立叶变换的计算,目前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(2x)。如果，( )样本值的采样间隔是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 2.0), a = a-2; f = flipud(f); end if (a1.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);,% 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);,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 = 2nextpow2(N); z = ifft( fft(x,P) .* fft(y,P); z = z(1:N);,
