分数阶傅里叶变换

1、分数阶傅里叶变换的MATLAB仿真计算以及几点讨论在Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computation of the fractional Fourier transform中给出了一种快速计算分数阶傅里叶变换的算法，其MATLAB计算程序可在.tr/haldun/fracF.m 上查到。现在基于该程序，对一方波进行计算仿真。注：网上流传较为广泛的FRFT计算程序更为简洁，据称也是Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computation of

2、 the fractional Fourier transform使用的算法。但是根据Adhemar Bultheel和 Hector E. Martnez Sulbaran的论文Computation of the Fractional Fourier Transform中提到，Ozaktas等人的分数阶傅里叶变换的计算程序仅有上述网站这一处，而两个程序的计算结果基本相符。本文使用较为简洁的计算程序，Ozaktas等人的计算程序在附表中给出。程序如下：clearclc %构造方波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%调用计算函数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

4、(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 transformerror(nargchk(2, 2, nargin);f = f(:);N = length(f);shft

5、= rem(0:N-1)+fix(N/2),N)+1;sN = sqrt(N);a = mod(a,4); % do special casesif (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.5if (a>2.0)

6、, 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 % 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 premultipl

7、icationchrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.2);f = chrp.*f; % chirp convolutionc = 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 multiplicationFaf = chrp.*Faf; % normalizing constantFaf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1); fun

8、ction 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(-(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);

9、 从图中可见，当旋转角度时，分数阶Fourier变换将收敛为方波信号；当时，收敛为函数。对于线性调频chirp信号Xk=exp(-j0.01141t2)，k=-32,-3132，变换后的信号波形图如下几点讨论一，目前的分数阶傅里叶变换主要有三种快速算法：1，B. Santhanamand和 J. H. McClellan的论文The discrete rotational Fourier transform中，先计算离散FRFT的核矩阵，再利用FFT来计算离散FRFT。2，本文中采用的在Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computati

10、on of the fractional Fourier transform所述的算法，是将FRFT分解为信号的卷积形式，从而利于FFT计算FRFT。3，Soo-Chang Pei和 Min-Hung Yeh等人在Two dimensional discrete fractionalFourier transform和Discrete frac-tional fourier transformbased on orthogonal projections中，采用矩阵的特征值和特征向量来计算FRFT。二，Ozaktas 在Digital computation of the fractional

11、 Fourier transform所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。在C. Candan和 M.A. Kutay等人的论文The discrete Fractional Fourier Transform中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（Error! Reference source not found.）二者吻合得很好。图 1 C. Candan和 M.A. Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较三，在Luis B. Almeida 的论文The Fractional Fourie

12、r Transform and Time Frequency Representations中给出了方波的分数阶傅立叶变换图形（图 2）图 2 Almeida 的论文中给出的方波的分数阶傅立叶变换图形该图形与讲义中的图形相符。本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论用那种算法程序，怎样调整输入信号，在时，虚部都不为零，这与Almeida和讲义中的图形并不一致。而在Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computation of the fractional Fourier transform中只给出了幅值的绝对值的图形，并

13、没有给出实部与虚部的结果，因此尚需进一步讨论图 3 本文中计算的时，实部与虚部分布附：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 scien

14、tific 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 Zalevsk

15、y, 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

16、 vector of the function whose%fractional Fourier transform is to be taken, and a' is the%transform order. The function returns the samples of the a'th%order fractional Fourier transform, under the assumption that%the Wigner distribution of the function is negligible outside a%circle whose di

17、ameter is the square root of the length of fc. % functionres=fracF(fc,a) % This function operates on the vector fc which is assumed to% 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 o

18、rigin. % (deltax2 is the time-bandwidth product of the function f.)% fc is assumed to have an even number of elements.% This function maps fc to a vector, whose elements are the samples % of the a'th order fractional Fourier transform of the function f. % The lengths of the input and ouput vecto

19、rs are the same if the % input vector has an even number of elements, as required.% Operating interval: -2 <= a <= 2% This function uses the core' function corefrmod2.m N = length(fc);% if fix(N/2) = N/2% error('Length of the input vector should be even');% end; fc = fc(:); fc = bi

20、zinter(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 (f

21、lag=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); % functionres=corefrmod2(fc,a) % Core function for

22、computing the fractional Fourier transform.% Valid only when 0.5 <= abs(a) <= 1.5% Decomposition used: % chirp mutiplication - chirp convolution - chirp mutiplication 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)

23、:fix(N/2)-1/deltax1;fc = fc(:);fc = fc(1:N);f1 = exp(-1i*pi*tan(phi/2)*x.*x); %multiplication by chirp!f1 = f1(:);fc = fc.*f1;x = x(:);clear x;t =-N+1:N-1/deltax1;hlptc =exp(i*pi*beta*t.*t);clear t;hlptc = hlptc(:); N2 = length(hlptc);N3 = 2(ceil(log(N2+N-1)/log(2);hlptcz = hlptc;zeros(N3-N2,1);fcz

24、= fc;zeros(N3-N,1);Hcfft = ifft(fft(fcz).*fft(hlptcz); % convolution with chirpclear hlptcz;clear fcz;Hc = Hcfft(N:2*N-1);clear Hcfft;clear hlptc; 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; % multiplication by chirp! if (fix(N/2) =N/2) res2(1:N-1) = res(2:N); res2(N) = res(1); res = res2;end; res = res(:); clear f1clear Hc % function xint=bizinter(x) N=length(x);im = 0;if sum(abs(imag(x)>0 im = 1; imx = imag(x); x = real(x);end; x2=x(:);x2=x2.' zeros(1,N);x2=x2(:)