分数阶傅里叶变换

分数阶傅里叶变换的MATLAB仿真计算以及几点讨论在HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》中给出了一种快速计算分数阶傅里叶变换的算法，其MATLAB计算程序可在~haldun/fracF.m上查到。现在基于该程序，对一方波进行计算仿真。其它,01,1)(ttx注：网上流传较为广泛的FRFT计算程序更为简洁，据称也是HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》使用的算法。但是根据AdhemarBultheel和HectorE.MartnezSulbaran的论文《ComputationoftheFractionalFourierTransform》中提到，Ozaktas等人的分数阶傅里叶变换的计算程序仅有上述网站这一处，而两个程序的计算结果基本相符。本文使用较为简洁的计算程序，Ozaktas等人的计算程序在附表中给出。程序如下：clearclc%构造方波其它,01,1)(ttxdt=0.05;T=20;t=-T:dt:T;n=length(t);m=1;fork=1:n;%tt=-36+k;tt=-T+k*dt;iftt=-m&&tt=mx(k)=1;elsex(k)=0;endend%确定α的值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]);分数阶傅里叶变换计算函数如下：functionFaf=frft(f,a)%ThefastFractionalFourierTransform%input:f=samplesofthesignal%a=fractionalpower%output:Faf=fastFractionalFouriertransformerror(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);%dospecialcasesif(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;endif(a==3),Faf(shft,1)=ifft(f(shft))*sN;return;end%reducetointerval0.5a1.5if(a2.0),a=a-2;f=flipud(f);endif(a1.5),a=a-1;f(shft,1)=fft(f(shft))/sN;endif(a0.5),a=a+1;f(shft,1)=ifft(f(shft))*sN;end%thegeneralcasefor0.5a1.5alpha=a*pi/2;tana2=tan(alpha/2);sina=sin(alpha);f=[zeros(N-1,1);interp(f);zeros(N-1,1)];%chirppremultiplicationchrp=exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);f=chrp.*f;%chirpconvolutionc=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);%chirppostmultiplicationFaf=chrp.*Faf;%normalizingconstantFaf=exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);functionxint=interp(x)%sincinterpolationN=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)%convolutionbyfftN=length([x(:);y(:)])-1;P=2^nextpow2(N);z=ifft(fft(x,P).*fft(y,P));z=z(1:N);从图中可见，当旋转角度时，分数阶Fourier变换将收敛为方波信号；0)(tx当时，收敛为函数。2csin对于线性调频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等人在《TwodimensionaldiscretefractionalFouriertransform》和《Discretefrac-tionalfouriertransformbasedonorthogonalprojections》中，采用矩阵的特征值和特征向量来计算FRFT。二，Ozaktas在《DigitalcomputationofthefractionalFouriertransform》所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。在C.Candan和M.A.Kutay等人的论文《ThediscreteFractionalFourierTransform》中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（Error!Referencesourcenotfound.）二者吻合得很好。图1C.Candan和M.A.Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较三，在LuisB.Almeida的论文《TheFractionalFourierTransformandTimeFrequencyRepresentations》中给出了方波的分数阶傅立叶变换图形（图2）图2Almeida的论文中给出的方波的分数阶傅立叶变换图形该图形与讲义中的图形相符。本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论用那种算法程序，怎样调整输入信号，在时，虚部2都不为零，这与Almeida和讲义中的图形并不一致。而在HaldunM.Ozaktas和OrhanArikan等人的论文《DigitalcomputationofthefractionalFouriertransform》中只给出了幅值的绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论图3本文中计算的时，实部与虚部分布2附：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