% 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 = freturnend
if (a==2), Faf = flipud(f)returnend
if (a==1), Faf(shft,1) = fft(f(shft))/sNreturnend
if (a==3), Faf(shft,1) = ifft(f(shft))*sNreturnend
% reduce to interval 0.5 <a <1.5
if (a>2.0), a = a-2f = flipud(f)end
if (a>1.5), a = a-1f(shft,1) = fft(f(shft))/sNend
if (a<0.5), a = a+1f(shft,1) = ifft(f(shft))*sNend
% 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 = 2^nextpow2(N)
z = ifft( fft(x,P) .* fft(y,P))
z = z(1:N)
首先,是关于分数阶傅里叶变换阶次的表示问题。现在很多研究者都习惯用分数阶傅里叶域相对时域的旋转角度来表示神誉册相对分数阶傅里叶域的阶次,这没有对“分数阶”三字的进行充分认识。所谓分数阶傅里叶域是指该变换域相对时域的旋转角度是90度的分数倍游宏,不同于以往的FFT、IFFT分别为90度的+1、-1倍,因此称该种变换为分数阶傅里叶变换,因此虚槐,使用旋转角度来表示分数阶的阶次极为不妥。其次,是现在很多研究者总喜欢拿分数阶傅里叶变换与短时傅里叶变换、小波变换以及威格纳变换这些时频分析做比较,从此怀疑分数阶傅里叶变换的实用价值。我认为,这是对分数阶傅里叶变化的广义时频分析特点认识不足所致。首先,必须声明的一点是分数阶傅里叶变换并不是传统意义上的时频分析,它只是一种广义的时频分析,它并没有完全解决传统傅里叶变换在时频定位以及分辨率上的局限性。所以,我们要做的工作是在分数阶傅里叶变换的基础上构建新的时频分析体系,例如短时分数阶傅里叶变换、分数阶小波包变换等。最后,就是现在部分研究者总是用D-CHIRP变换替代分数阶傅里叶变换,也就是说忽略最后一个chirp乘积。我认为这并没有将信号变换到分数阶傅里叶域,也就是说只是变换到传统的傅里叶域,并没有真正意义上用chirp信号做载波。在分数阶傅里叶变换可分为三步进行,第一步,乘以时间chirp基;第二步,做FFT运算;第三步,乘以频率chirp基。欢迎分享,转载请注明来源:内存溢出
评论列表(0条)