1、分數(shù)階傅里葉變換的MATLAB仿真計算以及幾點討論在Haldun M. Ozaktas 和 Orhan Arikan等人的論文Digital computation of the fractional Fourier transform中給出了一種快速計算分數(shù)階傅里葉變換的算法，其MATLAB計算程序可在.tr/haldun/fracF.m 上查到。現(xiàn)在基于該程序，對一方波進行計算仿真。注：網(wǎng)上流傳較為廣泛的FRFT計算程序更為簡潔，據(jù)稱也是Haldun M. Ozaktas 和 Orhan Arikan等人的論文Digital computation of

2、 the fractional Fourier transform使用的算法。但是根據(jù)Adhemar Bultheel和 Hector E. Martnez Sulbaran的論文Computation of the Fractional Fourier Transform中提到，Ozaktas等人的分數(shù)階傅里葉變換的計算程序僅有上述網(wǎng)站這一處，而兩個程序的計算結(jié)果基本相符。本文使用較為簡潔的計算程序，Ozaktas等人的計算程序在附表中給出。程序如下：clearclc %構(gòu)造方波dt=0.05;T=20;t=-T:dt:T;n=length(t);m=1;for k=1:n; % tt=-3

3、6+k; tt=-T+k*dt; if tt=-m & tt=m x(k)=1; else x(k)=0; endend%確定的值alpha=0.01;p=2*alpha/pi%調(diào)用計算函數(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時的實部和虛部);axis(-4,4,-1.5,2);subplot(2,2,2);plot(t,A,-);title(=0.01時的幅值);axis(-4,4,0,2);分數(shù)階傅里葉變換計算函

4、數(shù)如下：function Faf = frft(f, a)% The fast Fractional Fourier Transform% input: f = samples of the signal% a = fractional power% output: Faf = fast Fractional Fourier transformerror(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 case

5、sif (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 2.0), a = a-2; f = flipud(f); endif (a1.5), a = a-1; f(shft,1) = fft(f(shft)/sN; endif (a0.5), a

6、= a+1; f(shft,1) = ifft(f(shft)*sN; end % the general case for 0.5 a 1.5alpha = a*pi/2;tana2 = tan(alpha/2);sina = sin(alpha);f = zeros(N-1,1) ; interp(f) ; zeros(N-1,1); % chirp premultiplicationchrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2).2);f = chrp.*f; % chirp convolutionc = pi/N/sina/4;Faf = fconv

7、(exp(i*c*(-(4*N-4):4*N-4).2),f);Faf = Faf(4*N-3:8*N-7)*sqrt(c/pi); % chirp post multiplicationFaf = chrp.*Faf; % normalizing constantFaf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1); function xint=interp(x)% sinc interpolationN = length(x);y = zeros(2*N-1,1);y(1:2:2*N-1) = x;xint = fconv(y(1:2*N-1), sinc(-

8、(2*N-3):(2*N-3)/2);xint = xint(2*N-2:end-2*N+3); function z = fconv(x,y)% convolution by fftN = length(x(:);y(:)-1;P = 2nextpow2(N);z = ifft( fft(x,P) .* fft(y,P);z = z(1:N); 從圖中可見，當旋轉(zhuǎn)角度時，分數(shù)階Fourier變換將收斂為方波信號；當時，收斂為函數(shù)。對于線性調(diào)頻chirp信號Xk=exp(-j0.01141t2)，k=-32,-3132，變換后的信號波形圖如下幾點討論一，目前的分數(shù)階傅里葉變換主要有三種快速算法

9、：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 dis

10、crete fractionalFourier transform和Discrete frac-tional fourier transformbased on orthogonal projections中，采用矩陣的特征值和特征向量來計算FRFT。2， Ozaktas 在Digital computation of the fractional Fourier transform所述的算法，其實不是“離散”分數(shù)階傅里葉變換的算法，而是對連續(xù)分數(shù)階傅里葉變換的數(shù)值計算。在C. Candan和 M.A. Kutay等人的論文The discrete Fractional Fourier Tra

11、nsform中介紹了離散分數(shù)階傅里葉變換的算法，并給出了計算仿真圖形（Error! Reference source not found.）二者吻合得很好。圖 1 C. Candan和 M.A. Kutay等人離散分數(shù)階傅里葉變換的算法與連續(xù)分數(shù)階傅里葉變換的比較三，在Luis B. Almeida 的論文The Fractional Fourier Transform and Time Frequency Representations中給出了方波的分數(shù)階傅立葉變換圖形（圖 2）圖 2 Almeida 的論文中給出的方波的分數(shù)階傅立葉變換圖形該圖形與講義中的圖形相符。本文中的仿真結(jié)果大致與該

12、圖形也相符合，但是令人困惑的是無論用那種算法程序，怎樣調(diào)整輸入信號，在時，虛部都不為零，這與Almeida和講義中的圖形并不一致。而在Haldun M. Ozaktas 和 Orhan Arikan等人的論文Digital computation of the fractional Fourier transform中只給出了幅值的絕對值的圖形，并沒有給出實部與虛部的結(jié)果，因此尚需進一步討論圖 3 本文中計算的時，實部與虛部分布附：Haldun M. Ozaktas 和 Orhan Arikan等人的論文Digital computation of the fractional Fourier

13、 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

14、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, chapt

15、er 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%fractional Fourier transform is to be taken, and a is the%transform order. The function returns

16、 the samples of the ath%order fractional Fourier transform, under the assumption that%the Wigner distribution of the function is negligible outside a%circle whose diameter is the square root of the length of fc. % functionres=fracF(fc,a) % This function operates on the vector fc which is assumed to%

17、 be the samples of a function, obtained at a rate 1/deltax % where the Wigner distribution of the function f is confined% to a circle of diameter deltax around the origin. % (deltax2 is the time-bandwidth product of the function f.)% fc is assumed to have an even number of elements.% This function m

18、aps fc to a vector, whose elements are the samples % of the ath order fractional Fourier transform of the function f. % The lengths of the input and ouput vectors are the same if the % input vector has an even number of elements, as required.% Operating interval: -2 = a 0) & (a-0.5) & (a1.5) & (a-2)

19、 & (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;elseif (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); %

20、functionres=corefrmod2(fc,a) % Core function for computing the fractional Fourier transform.% Valid only when 0.5 = abs(a) 0 im = 1; imx = imag(x); x = real(x);end; x2=x(:);x2=x2.; zeros(1,N);x2=x2(:);xf=fft(x2);if rem(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