分数阶傅里叶变换Word文件下载.docx

1、end%确定的值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,=0.01时的幅值axis（-4,4,0,2）;分数阶傅里叶变换计算函数如下：function Faf = frft（f, a）% The fast Fractional Fourier Transform% i

2、nput: 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）;if （a=1）, Faf（sh

3、ft,1） = fft（f（shft）/sN; end if （a=3）, Faf（shft,1） = ifft（f（shft）*sN;% reduce to interval 0.5 a 2.0）, a = a-2; f = flipud（f）;1.5）, a = a-1; f（shft,1） = fft（f（shft）/sN;if （a0.5）, a = a+1; f（shft,1） = ifft（f（shft）*sN;% the general case for 0.5 alpha = a*pi/2;tana2 = tan（alpha/2）;sina = sin（alpha）;f = z

4、eros（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 constantFa

5、f = 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*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（

6、 fft（x,P） .* fft（y,P）;z = z（1:N）;从图中可见，当旋转角度时，分数阶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 和 Or

7、han Arikan等人的论文Digital computation 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 在Digit

8、al computation of the fractional Fourier transform所述的算法，其实不是“离散”分数阶傅里叶变换的算法，而是对连续分数阶傅里叶变换的数值计算。在C. Candan和 M.A. Kutay等人的论文The discrete Fractional Fourier Transform中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（错误！未找到引用源。）二者吻合得很好。图 1 C. Candan和 M.A. Kutay等人离散分数阶傅里叶变换的算法与连续分数阶傅里叶变换的比较三，在Luis B. Almeida 的论文The Fractiona

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

10、绝对值的图形，并没有给出实部与虚部的结果，因此尚需进一步讨论图 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 f

11、or 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

12、 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 th

13、e 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 ath%order fractional Fourier transform, under the assumption that%the Wigner distribution of the function is negligible outside a%circle whose diam

14、eter 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 origin

15、. % （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 ath order fractional Fourier transform of the function f. % The lengths of the input and ouput vectors are the

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

17、,1）;flag = 0;0） & （a0.5） flag = 1; a = a-1;end;-0.5） &0） flag = 2; a = a+1;1.5） &2） flag = 3;-2） &-1.5） flag = 4;res = fc;if （flag=1） | （flag=3） res = corefrmod2（fc,1）;if （flag=2） | （flag=4） res = corefrmod2（fc,-1）;if （a=0） res = fc;elseif （a=2） | （a=-2） res = flipud（fc）; res = corefrmod2（res,a）;res

18、 = res（N+1:3*N）;res = bizdec（res）;res（1） = 2*res（1）;functionres=corefrmod2（fc,a）% Core function for computing the fractional Fourier transform.% Valid only when 0.5 = abs（a） im = 1; imx = imag（x）; x = real（x）;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（

19、N/2）+1; xint=2*real（ifft（xf（1:N1）; zeros（N,1） ;xf（N2:2*N）.）;N/2）;xf（2*N-N/2+1:if （ im = 1） x2=imx（: x2=x2. x2=x2（: xf=fft（x2）; if rem（N,2）=1 %N = odd xmint=2*real（ifft（xf（1: xint = xint + i*xmint;xint = xint（:function xdec=bizdec（x）k = 1: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; for k = 1:N F2D（:,k） = fracF（f2D（:,k）,ac）;F2D = conj（F2Dif ar = 0M,k） = fracF（F2D（:,k）,ar）;