分数阶傅里叶变换

分数阶傅里叶变换的MATLAB 仿真计算以及几点讨论

在Haldun M. Ozaktas 和 Orhan Arikan 等人的论文《Digital computation of the fractional Fourier transform 》中给出了一种快速计算分数阶傅里叶变换的算法， 其MATLAB 计算程序可在https://www.360docs.net/doc/e29149289.html,.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 等人的计算程序在附表中给出。

程序如下：

clear

clc

%构造方波?????<=其它

,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<=m

x(k)=1;

else

x(k)=0;

end

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,'-');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 transform

error(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 cases

if (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.5

if (a>2.0), a = a-2; f = flipud(f); end

if (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.5

alpha = a*pi/2;

tana2 = tan(alpha/2);

sina = sin(alpha);

f = [zeros(N-1,1) ; interp(f) ; zeros(N-1,1)];

% chirp premultiplication

chrp = exp(-i*pi/N*tana2/4*(-2*N+2:2*N-2)'.^2);

f = chrp.*f;

% chirp convolution

c = 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 multiplication

Faf = chrp.*Faf;

% normalizing constant

Faf = exp(-i*(1-a)*pi/4)*Faf(N:2:end-N+1);

function xint=interp(x)

% sinc interpolation

N = 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 fft

N = length([x(:);y(:)])-1;

P = 2^nextpow2(N);

z = ifft( fft(x,P) .* fft(y,P));

z = z(1:N);

从图中可见，当旋转角度0→α时，分数阶Fourier 变换将收敛为方波信号)(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》中介绍了离散分数阶傅里叶变换的算法，并给出了计算仿真图形（错误！未找到引用源。）二者吻合得很好。

图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.m

N = 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;

else

if (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); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 f1

clear 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(:);

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

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*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');