分数阶傅里叶变换.docx

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

2、ation of 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;

3、 % tt=-36+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(t,A,-);title(=0.01时的幅值);axis(-4,4,0,2);分数阶

4、傅里叶变换计算函数如下：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 spec

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

6、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 premultiplicationchrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2).2);f = chrp.*f; % chirp convolutionc = pi/N/sina/4;Faf

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

8、, 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); 从图中可见，当旋转角度时，分数阶Fourier变换将收敛为方波信号；当时，收敛为函数。对于线性调频chirp信号Xk=exp(-j0.01141t2)，k=-32,-3132，变换后的信号波形图如下几点讨论一，目前的分数阶傅里叶变换主

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 dimensi

10、onal discrete fractionalFourier transform和Discrete frac-tional fourier transformbased on orthogonal projections中，采用矩阵的特征值和特征向量来计算FRFT。二，Ozaktas 在Digital computation of the fractional Fourier transform所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。在C. Candan和 M.A. Kutay等人的论文The discrete Fractional Four

11、ier Transform中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（错误！未找到引用源。）二者吻合得很好。图 1 C. Candan和 M.A. Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较三，在Luis B. Almeida 的论文The Fractional Fourier Transform and Time Frequency Representations中给出了方波的分数阶傅立叶变换图形（图 2）图 2 Almeida 的论文中给出的方波的分数阶傅立叶变换图形该图形与讲义中的图形相符。本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论

12、用那种算法程序，怎样调整输入信号，在时，虚部都不为零，这与Almeida和讲义中的图形并不一致。而在Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computation of the fractional Fourier transform中只给出了幅值的绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论图 3 本文中计算的时，实部与虚部分布附：Haldun M. Ozaktas 和 Orhan Arikan等人的论文Digital computation of the fractional Fourier transform所述的算法程序

13、%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 Kutay, and Gozde

14、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, chapter 6, page 298.%T

15、he 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 the samples of t

16、he 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% be the samples o

17、f 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 maps fc to a vecto

18、r, 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) & (a-1.5) flag =

19、 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); % functionres=coref

20、rmod2(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

21、).);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); if rem(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

22、*N-N/2+1:2*N).); end; xint = xint + i*xmint;end; xint = xint(:); % function xdec=bizdec(x) k = 1:2:length(x);xdec = x(k); xdec = xdec(:); % function F2D=fracF2D(f2D,ac,ar) M,N = size(f2D);F2D = zeros(M,N); if ac = 0 F2D = f2D;else for k = 1:N F2D(:,k) = fracF(f2D(:,k),ac); end;end; F2D = conj(F2D); if ar = 0 for k = 1:M F2D(:,k) = fracF(F2D(:,k),ar); end;end; F2D = conj(F2D);