alpha=a*pi/2;
tana2=tan（alpha/2）;
sina=sin（alpha）;
f=[zeros（N-1,1）;interp（f）;zeros（N-1,1）];
%chirppremultiplication
chrp=exp（-i*pi/N*tana2/4*（-2*N+2:
2*N-2）'.^2）;
f=chrp.*f;
%chirpconvolution
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）;
%chirppostmultiplication
Faf=chrp.*Faf;
%normalizingconstant
Faf=exp（-i*（1-a）*pi/4）*Faf（N:
2:
end-N+1）;
functionxint=interp（x）
%sincinterpolation
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）;
functionz=fconv（x,y）
%convolutionbyfft
N=length（[x（:
）;y（:
）]）-1;
P=2^nextpow2（N）;
z=ifft（fft（x,P）.*fft（y,P））;
z=z（1:
N）;
从图中可见，当旋转角度
时，分数阶Fourier变换将收敛为方波信号
；当
时，收敛为
函数。
对于线性调频chirp信号Xk=exp（-j0.01141t2），k=-32,-31……32，变换后的信号波形图如下
几点讨论
一，目前的分数阶傅里叶变换主要有三种快速算法：
1，B.Santhanamand和J.H.McClellan的论文《ThediscreterotationalFouriertransform》中，先计算离散FRFT的核矩阵，再利用FFT来计算离散FRFT。
2，本文中采用的在HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》所述的算法，是将FRFT分解为信号的卷积形式，从而利于FFT计算FRFT。
3，Soo-ChangPei和Min-HungYeh等人在《Twodimensionaldiscretefractional
Fouriertransform》和《Discretefrac-tionalfouriertransformbasedonorthogonalprojections》中，采用矩阵的特征值和特征向量来计算FRFT。
二，Ozaktas在《DigitalcomputationofthefractionalFouriertransform》所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。
在C.Candan和M.A.Kutay等人的论文《ThediscreteFractionalFourierTransform》中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（错误！
未找到引用源。
）二者吻合得很好。
图1C.Candan和M.A.Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较
三，在LuisB.Almeida的论文《TheFractionalFourierTransformandTimeFrequencyRepresentations》中给出了方波的分数阶傅立叶变换图形（图2）
图2Almeida的论文中给出的方波的分数阶傅立叶变换图形
该图形与讲义中的图形相符。
本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论用那种算法程序，怎样调整输入信号，在
时，虚部都不为零，这与Almeida和讲义中的图形并不一致。
而在HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》中只给出了幅值的绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论
图3本文中计算的
时，实部与虚部分布
附：
HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》所述的算法程序
%FASTCOMPUTATIONOFTHEFRACTIONALFOURIERTRANSFORM
%byM.AlperKutay,September1996,Ankara
%Copyright1996M.AlperKutay
%Thiscodemaybeusedforscientificandeducationalpurposes
%providedcreditisgiventothepublicationsbelow:
%
%HaldunM.Ozaktas,OrhanArikan,M.AlperKutay,andGozdeBozdagi,
%DigitalcomputationofthefractionalFouriertransform,
%IEEETransactionsonSignalProcessing,44:
2141--2150,1996.
%HaldunM.Ozaktas,ZeevZalevsky,andM.AlperKutay,
%TheFractionalFourierTransformwithApplicationsinOpticsand
%SignalProcessing,Wiley,2000,chapter6,page298.
%
%Theseveralfunctionsgivenbelowshouldbeseparatelysaved
%underthesamedirectory.fracF（fc,a）isthefunctiontheuser
%shouldcall,wherefcisthesamplevectorofthefunctionwhose
%fractionalFouriertransformistobetaken,and`a'isthe
%transformorder.Thefunctionreturnsthesamplesofthea'th
%orderfractionalFouriertransform,undertheassumptionthat
%theWignerdistributionofthefunctionisnegligibleoutsidea
%circlewhosediameteristhesquarerootofthelengthoffc.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[res]=fracF（fc,a）
%Thisfunctionoperatesonthevectorfcwhichisassumedto
%bethesamplesofafunction,obtainedatarate1/deltax
%wheretheWignerdistributionofthefunctionfisconfined
%toacircleofdiameterdeltaxaroundtheorigin.
%（deltax^2isthetime-bandwidthproductofthefunctionf.）
%fcisassumedtohaveanevennumberofelements.
%Thisfunctionmapsfctoavector,whoseelementsarethesamples
%ofthea'thorderfractionalFouriertransformofthefunctionf.
%Thelengthsoftheinputandouputvectorsarethesameifthe
%inputvectorhasanevennumberofelements,asrequired.
%Operatinginterval:
-2<=a<=2
%Thisfunctionusesthe`core'functioncorefrmod2.m
N=length（fc）;
%iffix（N/2）~=N/2
%error（'Lengthoftheinputvectorshouldbeeven'）;
%end;
fc=fc（:
）;
fc=bizinter（fc）;
fc=[zeros（N,1）;fc;zeros（N,1）];
flag=0;
if（a>0）&&（a<0.5）
flag=1;
a=a-1;
end;
if（a>-0.5）&&（a<0）
flag=2;
a=a+1;
end;
if（a>1.5）&&（a<2）
flag=3;
a=a-1;
end;
if（a>-2）&&（a<-1.5）
flag=4;
a=a+1;
end;
res=fc;
if（flag==1）||（flag==3）
res=corefrmod2（fc,1）;
end;
if（flag==2）||（flag==4）
res=corefrmod2（fc,-1）;
end;
if（a==0）
res=fc;
else
if（a==2）||（a==-2）
res=flipud（fc）;
else
res=corefrmod2（res,a）;
end;
end;
res=res（N+1:
3*N）;
res=bizdec（res）;
res
（1）=2*res
（1）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[res]=corefrmod2（fc,a）
%CorefunctionforcomputingthefractionalFouriertransform.
%Validonlywhen0.5<=abs（a）<=1.5
%Decompositionused:
%chirpmutiplication-chirpconvolution-chirpmutiplication
deltax=sqrt（length（fc））;
phi=a*pi/2;
N=fix（length（fc））;
deltax1=deltax;
alpha=1/tan（phi）;
beta=1/sin（phi）;
x=[-ceil（N/2）:
fix（N/2）-1]/deltax1;
fc=fc（:
）;
fc=fc（1:
N）;
f1=exp（-1i*pi*tan（phi/2）*x.*x）;%multiplicationbychirp!
f1=f1（:
）;
fc=fc.*f1;
x=x（:
）;
clearx;
t=[-N+1:
N-1]/deltax1;
hlptc=exp（i*pi*beta*t.*t）;
cleart;
hlptc=hlptc（:
）;
N2=length（hlptc）;
N3=2^（ceil（log（N2+N-1）/log
（2）））;
hlptcz=[hlptc;zeros（N3-N2,1）];
fcz=[fc;zeros（N3-N,1）];
Hcfft=ifft（fft（fcz）.*fft（hlptcz））;%convolutionwithchirp
clearhlptcz;
clearfcz;
Hc=Hcfft（N:
2*N-1）;
clearHcfft;
clearhlptc;
Aphi=exp（-i*（pi*sign（sin（phi））/4-phi/2））/sqrt（abs（sin（phi）））;
xx=[-ceil（N/2）:
fix（N/2）-1]/deltax1;
f1=f1（:
）;
res=（Aphi*f1.*Hc）/deltax1;%multiplicationbychirp!
if（fix（N/2）~=N/2）
res2（1:
N-1）=res（2:
N）;
res2（N）=res
（1）;
res=res2;
end;
res=res（:
）;
clearf1
clearHc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionxint=bizinter（x）
N=length（x）;
im=0;
ifsum（abs（imag（x）））>0
im=1;
imx=imag（x）;
x=real（x）;
end;
x2=x（:
）;
x2=[x2.';zeros（1,N）];
x2=x2（:
）;
xf=fft（x2）;
ifrem（N,2）==1%N=odd
N1=fix（N/2+1）;N2=2*N-fix（N/2）+1;
xint=2*real（ifft（[xf（1:
N1）;zeros（N,1）;xf（N2:
2*N）].'））;
else
xint=2*real（ifft（[xf（1:
N/2）;zeros（N,1）;xf（2*N-N/2+1:
2*N）].'））;
end;
if（im==1）
x2=imx（:
）;
x2=[x2.';zeros（1,N）];
x2=x2（:
）;
xf=fft（x2）;
ifrem（N,2）==1%N=odd
N1=fix（N/2+1）;N2=2*N-fix（N/2）+1;
xmint=2*real（ifft（[xf（1:
N1）;zeros（N,1）;xf（N2:
2*N）].'））;
else
xmint=2*real（ifft（[xf（1:
N/2）;zeros（N,1）;xf（2*N-N/2+1:
2*N）].'））;
end;
xint=xint+i*xmint;
end;
xint=xint（:
）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionxdec=bizdec（x）
k=1:
2:
length（x）;
xdec=x（k）;
xdec=xdec（:
）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionF2D=fracF2D（f2D,ac,ar）
[M,N]=size（f2D）;
F2D=zeros（M,N）;
ifac==0
F2D=f2D;
else
fork=1:
N
F2D（:
k）=fracF（f2D（:
k）,ac）;
end;
end;
F2D=conj（F2D'）;
ifar~=0
fork=1:
M
F2D（:
k）=fracF（F2D（:
k）,ar）;
end;
end;
F2D=conj（F2D'）;
升级会员 