分数阶傅里叶变换

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

在Haldun M. Ozaktas 和 Orhan Arikan等人的论文《Digital computation of the

fractional Fourier transform》中给出了一种快速计算分数阶傅里叶变换的算法，

其MATLAB计算程序可在.tr/~haldun/fracF.m 上查到。现在基于该程序，对一方波其它,01,1)(ttx进行计算仿真。

注：网上流传较为广泛的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)(ttx

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

z = ifft( fft(x,P) .* fft(y,P));

z = z(1:N);

从图中可见，当旋转角度0时，分数阶Fourier变换将收敛为方波信号)(tx；当2时，收敛为csin函数。

对于线性调频chirp信号Xk=exp(-j0.01141t2)，k=-32,-31……32，变换后的信号波形图如下

几点讨论

一，目前的分数阶傅里叶变换主要有三种快速算法：

1，B. Santhanamand和 J. H. McClellan的论文《The discrete rotational Fourier

transform》中，先计算离散FRFT的核矩阵，再利用FFT来计算离散FRFT。

2，本文中采用的在Haldun M. Ozaktas 和 Orhan Arikan等人的论文《Digital

computation of the fractional Fourier transform》所述的算法，是将FRFT分解为信号的卷积形式，从而利于FFT计算FRFT。

3，Soo-Chang Pei和 Min-Hung Yeh等人在《Two dimensional discrete fractional

Fourier transform》和《Discrete frac-tional fourier transformbased on orthogonal

projections》中，采用矩阵的特征值和特征向量来计算FRFT。

二，Ozaktas 在《Digital computation of the fractional Fourier transform》所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。在C. Candan和 M.A. Kutay等人的论文《The discrete Fractional

Fourier Transform》中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（错误!未找到引用源。）二者吻合得很好。

图 1 C. Candan和 M.A. Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较

三，在Luis B. Almeida 的论文《The Fractional Fourier Transform and Time

Frequency Representations》中给出了方波的分数阶傅立叶变换图形（图 2）

图 2 Almeida 的论文中给出的方波的分数阶傅立叶变换图形

该图形与讲义中的图形相符。本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论用那种算法程序，怎样调整输入信号，在2时，虚部都不为零，这与Almeida和讲义中的图形并不一致。而在Haldun M. Ozaktas 和 Orhan

Arikan等人的论文《Digital computation of the fractional Fourier transform》中只给出了幅值的绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论

图 3 本文中计算的2时，实部与虚部分布 附：

Haldun M. Ozaktas 和 Orhan Arikan等人的论文《Digital computation of the

fractional Fourier transform》所述的算法程序

%FAST COMPUTATION OF THE FRACTIONAL FOURIER TRANSFORM

%by M. Alper Kutay, September 1996, Ankara

%Copyright 1996 M. Alper Kutay

%This code may be used for scientific and educational purposes

%provided credit is given to the publications below:

%

%Haldun M. Ozaktas, Orhan Arikan, M. Alper Kutay, and Gozde Bozdagi,

%Digital computation of the fractional Fourier transform,

%IEEE Transactions on Signal Processing, 44:2141--2150, 1996.

%Haldun M. Ozaktas, Zeev Zalevsky, and M. Alper Kutay,

%The Fractional Fourier Transform with Applications in Optics and

%Signal Processing, Wiley, 2000, chapter 6, page 298.

%

%The several functions given below should be separately saved

%under the same directory. fracF(fc,a) is the function the user

%should call, where fc is the sample vector of the function whose