分数阶傅里叶变换

分数阶傅里叶变换的MATLAB 仿真计算以及几点讨论

在Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》中给出了一种快速计算分数阶傅里叶变换的算法，其MATLAB 计算程序可在.tr/~haldun/fracF.m 上查到。现在

基于该程序，对一方波进行计算仿真。⎪⎩⎪⎨⎧<=其它

,01,1)(t t x 注：网上流传较为广泛的FRFT 计算程序更为简洁，据称也是Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》使用的算法。但是根据Adhemar Bultheel 和 Hector E. Martnez Sulbaran 的论文《Computation of the Fractional Fourier Transform 》中提到，Ozaktas 等人的分数阶傅里叶变换的计算程序仅有上述网站这一处，而两个程序的计算结果基本相符。本文使用较为简洁的计算程序，Ozaktas 等人的计算程序在附表中给出。

程序如下：

clear

clc

%构造方波⎪⎩⎪⎨⎧<=其它

,01

,1)(t t x dt=0.05;

T=20;

t=-T:dt:T;

n=length(t);

m=1;

for k=1:n;

% tt=-36+k;

tt=-T+k*dt;

if tt>=-m && tt<=m

x(k)=1;

else

x(k)=0;

end

end

%确定α的值

alpha=0.01;

p=2*alpha/pi

%调用计算函数

Fx=frft(x,p);

Fx=Fx';

Fr=real(Fx);

Fi=imag(Fx);

A=abs(Fx);

figure,

subplot(2,2,1);

plot(t,Fr,'-',t,Fi,':');title(' α=0.01时的实部和虚部π'); axis([-4,4,-1.5,2]);

subplot(2,2,2);

plot(t,A,'-');title('α=0.01时的幅值');

axis([-4,4,0,2]);

分数阶傅里叶变换计算函数如下：

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 = 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);