分数阶傅里叶变换

分数阶傅里叶变换的MATLAB 仿真计算以及几点讨论在Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》中给出了一种快速计算分数阶傅里叶变换的算法，其MATLAB 计算程序可在.tr/~haldun/fracF.m 上查到。

现在基于该程序，对一方波进行计算仿真。

⎪⎩⎪⎨⎧<=其它,01,1)(t t x 注：网上流传较为广泛的FRFT 计算程序更为简洁，据称也是Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》使用的算法。

但是根据Adhemar Bultheel 和 Hector E. Martnez Sulbaran 的论文《Computation of the Fractional Fourier Transform 》中提到，Ozaktas 等人的分数阶傅里叶变换的计算程序仅有上述网站这一处，而两个程序的计算结果基本相符。

本文使用较为简洁的计算程序，Ozaktas 等人的计算程序在附表中给出。

程序如下：clearclc%构造方波⎪⎩⎪⎨⎧<=其它,01,1)(t t x dt=0.05;T=20;t=-T:dt:T;n=length(t);m=1;for k=1:n;% tt=-36+k;tt=-T+k*dt;if tt>=-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]);分数阶傅里叶变换计算函数如下：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 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), a = a-2; f = flipud(f); endif (a>1.5), a = a-1; f(shft,1) = fft(f(shft))/sN; end if (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 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(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([-(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 = 2^nextpow2(N);z = ifft( fft(x,P) .* fft(y,P)); z = z(1:N);从图中可见，当旋转角度时，分数阶Fourier 变换将收敛为方波信号；0→α)(t x 当时，收敛为函数。

2πα→c sin 对于线性调频chirp 信号X k =exp(-j0.01141t 2)，k=-32,-31……32，变换后的信号波形图如下几点讨论一，目前的分数阶傅里叶变换主要有三种快速算法：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 discrete fractional Fourier 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 Fourier Transform》中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（Error! Reference source not found.）二者吻合得很好。

图 1 C. Candan和M.A. Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较三，在Luis B. Almeida 的论文《The Fractional Fourier Transform and Time Frequency Representations》中给出了方波的分数阶傅立叶变换图形（图2）图 2 Almeida 的论文中给出的方波的分数阶傅立叶变换图形该图形与讲义中的图形相符。

本文中的仿真结果大致与该图形也相符合，但是令人困惑的是无论用那种算法程序，怎样调整输入信号，在时，虚部2πα→都不为零，这与Almeida 和讲义中的图形并不一致。

而在Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》中只给出了幅值的绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论图 3 本文中计算的时，实部与虚部分布2πα→附：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 scientific 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 Zalevsky, 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 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 diameter is the square root of the length of fc. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[res]=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 origin.% (deltax^2 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 vectors 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.mN = length(fc);% if fix(N/2) ~= N/2% error('Length of the input vector should be even'); % end;fc = fc(:);fc = bizinter(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 (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);elseres = corefrmod2(res,a);end;end;res = res(N+1:3*N);res = bizdec(res);res(1) = 2*res(1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[res]=corefrmod2(fc,a)% Core function for 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):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 = [fc;zeros(N3-N,1)];Hcfft = ifft(fft(fcz).*fft(hlptcz)); % convolution with chirp clear 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)))>0im = 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 = oddN1=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)].'));elsexint=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 = oddN1=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)].'));elsexmint=2*real(ifft([xf(1:N/2); zeros(N,1) ;xf(2*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 == 0F2D = f2D;elsefor k = 1:NF2D(:,k) = fracF(f2D(:,k),ac);end;end;F2D = conj(F2D');if ar ~= 0for k = 1:MF2D(:,k) = fracF(F2D(:,k),ar);end;end;F2D = conj(F2D');。