分?jǐn)?shù)階傅里葉變換

《分?jǐn)?shù)階傅里葉變換》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫(kù)。

1、分?jǐn)?shù)階傅里葉變換的MATLAB仿真計(jì)算以及幾點(diǎn)討論在HaldunM.Ozaktas和OrhanArikan等人的論文《DigitalcomputationofthefractionalFouriertransform》中給出了一種快速計(jì)算分?jǐn)?shù)階傅里葉變換的算法，其MATLAB計(jì)算程序可在www.ee.bilkent.edu.tr/~haldun/fracF.m上查到?，F(xiàn)在基于該程序，對(duì)一方波進(jìn)行計(jì)算仿真。注：網(wǎng)上流傳較為廣泛的FRFT計(jì)算程序更為簡(jiǎn)潔，據(jù)稱也是HaldunM.Ozaktas和OrhanArikan等人的論文

2、《DigitalcomputationofthefractionalFouriertransform》使用的算法。但是根據(jù)AdhemarBultheel和HectorE.MartnezSulbaran的論文《ComputationoftheFractionalFourierTransform》中提到，Ozaktas等人的分?jǐn)?shù)階傅里葉變換的計(jì)算程序僅有上述網(wǎng)站這一處，而兩個(gè)程序的計(jì)算結(jié)果基本相符。本文使用較為簡(jiǎn)潔的計(jì)算程序，Ozaktas等人的計(jì)算程序在附表中給出。程序如下：clearclc%構(gòu)造方波dt=0.05;T=20

3、;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%調(diào)用計(jì)算函數(shù)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時(shí)的實(shí)部和虛部π');axis([-4

4、,4,-1.5,2]);subplot(2,2,2);plot(t,A,'-');title('α=0.01時(shí)的幅值');axis([-4,4,0,2]);分?jǐn)?shù)階傅里葉變換計(jì)算函數(shù)如下：functionFaf=frft(f,a)%ThefastFractionalFourierTransform%input:f=samplesofthesignal%a=fractionalpower%output:Faf=fastFractionalFouriertransformerror(nargchk(2,2,nargin));f=

5、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.5

6、.5if(a>2.0),a=a-2;f=flipud(f);endif(a>1.5),a=a-1;f(shft,1)=fft(f(shft))/sN;endif(a<0.5),a=a+1;f(shft,1)=ifft(f(shft))*sN;end%thegeneralcasefor0.5

7、rp=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);functio

8、nxint=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+