分数阶傅立叶变换程序汇总（个人收集自网上）
都是从网上收集来的，由于时间比较久，处处都忘记了，如果是谁的原创请和我联系，我在帖子中标出来的第二楼：the fast fractional Fourier transform
程序参考下文中的算法
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.
第三楼：分数阶fourior变换的源码，主要用来处理线性调频信号
第四楼：离散分数阶傅立叶变换，算法参考
C. Candan, M.A. Kutay, and H.M. Ozaktas.
The discrete Fractional Fourier Transform.
IEEE Trans. Sig. Proc., 48:1329--1337, 2000
第五楼：the discrete fractional Fourier transform，算法参考：
S.-C. Pei, M.-H. Yeh, and C.-C. Tseng.
Discrete fractional Fourier-transform based on orthogonal projections.
IEEE Trans. Sig. Proc., 47(5):1335--1348, 1999.
第六楼：a demonstration programfor the previous three routines
第七楼： the discrete fractional Cosine transform，算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001
第八楼：the discrete fractional Sine transform，算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001
第九楼：test for Disfrct.m，算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001
第十楼：the rescaling preprocessing for the frft routine，详细情况参考
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.
[ 本帖最后由 simon21 于 2007-4-26 20:14 编辑 ] 程序参考下文中的算法
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.
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 = ;
% 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()-1;
P = 2^nextpow2(N);
z = ifft( fft(x,P) .* fft(y,P));
z = z(1:N);
程序来自于http://www.cs.kuleuven.ac.be/cwis/research/nalag/research/software/FRFT/frft.m 分数阶fourior变换的源码，主要用来处理线性调频信号
%DISCRETE FRACTIONAL FOURIER TRANSFORM MATRIX GENERATOR
%by Cagatay Candan <candan@ieee.org>, July 1998, Ankara
%Copyright 1998 Cagatay Candan
%This code may be used for scientific and educational purposes
%provided credit is given to the publications below:
%
%This Matlab function generates the discrete fractional
%Fourier transform matrix originally described in:
%Cagatay Candan, M. Alper Kutay, and Haldun M. Ozaktas,
%The discrete fractional Fourier Transform,
%IEEE Transactions on Signal Processing, 48:1329-1337, May 2000,
%(also in Proc ICASSP'99, pages 1713-1716, IEEE, 1999);
%and further described in:
%Haldun M. Ozaktas, Zeev Zalevsky, and M. Alper Kutay,
%The Fractional Fourier Transform with Applications in Optics and
%Signal Processing, Wiley, 2000, chapter 6.
function F=dFRT(N,a,ord)
%function F=dFRT(N,a,ord) returns the NxN discrete fractional
%Fourier transform matrix with fractional order 'a'.
%The optional argument 'ord' is the order of approximation
%of the S matrix (default 2).
%Note: This Matlab file has some subfunctions for generating S_{2k}
% matrices, eigenvector ordering etc. These functions are not
% visible from the Matlab workspace.
global Evec Eval ordp
if nargin==2, ord=2;end;
if (length(Evec)~=N | ordp~=ord),
=dis_s(N,ord);
ordp=ord;
end;
even=~rem(N,2);
F=Evec*diag(exp(-j*pi/2*a*()))*Evec';
%%%%%%
function M=cconvm(v);
%Generates circular Convm matrix
v=v(:);N=length(v);
M=zeros(N,N);dum=v;
for k=1:N,
M(:,k)=dum;
dum=;
end;
%%%%%%
function S=creates(N,ord)
%Creates S matrix of approximation order ord
%When ord=1, elementary S matrix is returned
ord=ord/2;
dum=;s=0;
for k=1:ord,
s=(-1)^(k-1)*prod(1:(k-1))^2/prod(1:2*k)*2*+s;
dum=conv(dum,);
end;
S=cconvm(s)+diag(real(fft(s)));
%%%%%%
function =dis_s(N,ord)
%function =dis_s(N)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors
if nargin==1, ord=2;end;
%%Construct S Matrix
%S=diag(2*cos(2*pi/N*()))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
%S(1,N)=1;S(N,1)=1;
%%
S=creates(N,ord);
%%%%%%
%Construct P matrix
p=N;
r=floor(p/2);
P=zeros(p,p);
P(1,1)=1;
even=~rem(p,2);
for k=1:r-even,
P(k+1,k+1)=1/sqrt(2);
P(k+1,p-(k+1)+2)=1/sqrt(2);
end;
if (even), P(r+1,r+1)=1; end;
for k=r+1:p-1,
P(k+1,k+1)=-1/sqrt(2);
P(k+1,p-(k+1)+2)=1/sqrt(2);
end;
%%%%%%
CS=P*S*P';C2=CS(floor(1:N/2+1),floor(1:N/2+1));S2=CS(floor(N/2+2):N,floor(N/2+2):N);
=eig(C2);=eig(S2);
qvc=;
SC2=P*qvc; %Even Eigenvector of S
qvs=;
SS2=P*qvs; %Odd Eigenvector of S
es=diag(es);ec=diag(ec);
=sort(-ec);
SC2=SC2(:,in);
ec=ec(in);
=sort(-es);
SS2=SS2(:,in);
es=es(in);
if rem(N,2)==0,
S2C2=zeros(N,N+1);
SS2(:,size(SS2,2)+1)=zeros(N,1);
S2C2(:,+1)=SC2;
S2C2(:,+1)=SS2;
S2C2(:,N)=[];
else
S2C2=zeros(N,N);
S2C2(:,+1)=SC2;
S2C2(:,+1)=SS2;
end;
Evec=S2C2;
Eval=(-j).^[ 0:N-2 (N-1)+even];
%END 离散分数阶傅立叶变换- Separate score step Fourier transforms
算法参考
C. Candan, M.A. Kutay, and H.M. Ozaktas.
The discrete Fractional Fourier Transform.
IEEE Trans. Sig. Proc., 48:1329--1337, 2000
function Fa=DFPei(f,a)
% Compute discrete fractional Fourier transform
% of order a of vector f according to Pei/Yeh/Tseng
N=length(f); f=f(:);
shft = rem((0:N-1)+fix(N/2),N)+1;
global hn_saved p_saved
if (nargin==2), p = 2; end;
p = min(max(2,p),N-1);
if (length(hn_saved) ~= N | p_saved ~= p),
hn = make_hn(N,p);
hn_saved = hn; p_saved = p;
else
hn = hn_saved;
end;
Fa(shft,1)=hn*(exp(-j*pi/2*a*).'.*(hn'*f(shft)));
function hn = make_hn(N,p)
even = rem(N,2)==0;
shft = rem((0:N-1)+fix(N/2),N)+1;
% Gauss-Hermite samples
u = (-N/2:(N/2)-1)'/sqrt(N/2/pi);
ex = exp(-u.^2/2);
hn(:,1) = ex; r = norm(hn(:,1)); hn(:,1) = hn(:,1)/r;
hn(:,2)=2*u.*ex; s = norm(hn(:,2)); hn(:,2) = hn(:,2)/s;
r = s/r;
for k = 3:N+1
hn(:,k)=2*r*u.*hn(:,k-1)-2*(k-2)*hn(:,k-2);
s = norm(hn(:,k)); hn(:,k) = hn(:,k)/s;
r=s/r;
end
if (even), hn(:,N)=[]; else, hn(:,N+1)=[]; end
hn(shft,:) = hn;
% eigenvectors of DFT matrix
E = make_E(N,N/2);
for k = 1:4
if even % N even
switch k
case {1,3}
indx = k:4:N+1;
if (rem(N,4) ~= 0 && k==3) || (rem(N,4) == 0 && k==1),
indx(end) = indx(end)-1;
end
case {2,4}
indx = k:4:N-1;
end
else % N odd
indx = k:4:N;
end
hn(:,k:4:N) = orth(E(:,indx)*E(:,indx)'*hn(:,indx));
end
function E = make_E(N,p)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors
%Construct matrix H, use approx order ord
d2 = ; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
d_p = conv(d2,d_p);
st() = d_p; st(1) = 0;
s = s + (-1)^(k-1)*prod(1:(k-1))^2/prod(1:2*k)*2*st;
end;
% H = circulant + diagonal
col = (0:N-1)'; row = (N:-1:1);
idx = col(:,ones(N,1)) + row(ones(N,1),:);
st = ;
H = st(idx)+diag(real(fft(s)));
%Construct transformation matrix V
r = floor(N/2);
even = ~rem(N,2);
V1 = (eye(N-1)+flipud(eye(N-1)))/sqrt(2);
V1(N-r:end,N-r:end) = -V1(N-r:end,N-r:end);
if (even), V1(r,r)=1; end
V = eye(N); V(2:N,2:N) = V1;
% Compute eigenvectors
VHV = V*H*V';
E = zeros(N);
Ev = VHV(1:r+1,1:r+1); Od = VHV(r+2:N,r+2:N);
=eig(Ev); =eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee)); = sort(diag(eo));
%ve = ve(:,inde'); vo = vo(:,indo');
E(1:r+1,1:r+1) = fliplr(ve); E(r+2:N,r+2:N) = fliplr(vo);
E = V*E;
%shuffle eigenvectors
ind = ; ind = ind(:);
if (even), ind()=[]; else ind(N+1)=[]; end
E = E(:,ind');本程序来自http://www.cs.kuleuven.ac.be/cwi ... ftware/FRFT/cdpei.m
[ 本帖最后由 simon21 于 2007-4-26 19:57 编辑 ] 算法参考
S.-C. Pei, M.-H. Yeh, and C.-C. Tseng.
Discrete fractional Fourier-transform based on orthogonal projections.
IEEE Trans. Sig. Proc., 47(5):1335--1348, 1999.
function y=Disfrft(f,a,p)
% Computes discrete fractional Fourier transform
% of order a of vector x
% p (optional) is order of approximation, default N/2
N = length(f); even = ~rem(N,2);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRFT(N,p);
y(shft,1)=E*(exp(-j*pi/2*a*()).'.*(E'*f(shft)));
function E=dFRFT(N,p)
%function E=dFRFT(N,a,p) returns the NxN eigenvectors of the
%Fourier transform matrix
%The optional argument p is the order of approximation
global E_saved p_saved
if (length(E_saved) ~= N | p_saved ~= p),
E = make_E(N,p);
E_saved = E; p_saved = p;
else
E = E_saved;
end;
function E = make_E(N,p)
%Returns sorted eigenvectors and eigenvalues of corresponding vectors
%Construct matrix H, use approx order ord
d2 = ; d_p = 1; s = 0; st = zeros(1,N);
for k = 1:p/2,
d_p = conv(d2,d_p);
st() = d_p; st(1) = 0;
temp = ; temp=temp(:)'./;
s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
% H = circulant + diagonal
col = (0:N-1)'; row = (N:-1:1);
idx = col(:,ones(N,1)) + row(ones(N,1),:);
st = ;
H = st(idx)+diag(real(fft(s)));
%Construct transformation matrix V
r = floor(N/2);
even = ~rem(N,2);
V1 = (eye(N-1)+flipud(eye(N-1)))/sqrt(2);
V1(N-r:end,N-r:end) = -V1(N-r:end,N-r:end);
if (even), V1(r,r)=1; end
V = eye(N); V(2:N,2:N) = V1;
% Compute eigenvectors
VHV = V*H*V';
E = zeros(N);
Ev = VHV(1:r+1,1:r+1); Od = VHV(r+2:N,r+2:N);
=eig(Ev); =eig(Od);
%malab eig returns sorted eigenvalues
%if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee)); = sort(diag(eo));
%ve = ve(:,inde'); vo = vo(:,indo');
E(1:r+1,1:r+1) = fliplr(ve); E(r+2:N,r+2:N) = fliplr(vo);
E = V*E;
%shuffle eigenvectors
ind = ; ind = ind(:);
if (even), ind()=[]; else ind(N+1)=[]; end
E = E(:,ind');
[ 本帖最后由 simon21 于 2007-4-26 20:08 编辑 ] a demonstration programfor the previous three routines
x=0.0:0.02:2*pi; y =cos(x);
clear p_saved hn_saved E_saved
for a=0:0.05:4
fy=Disfrft(y,a);
fys=cdpei(y,a);
fyss=frft(y,a);
% blue,green,red,cyan
figure(1)
subplot(311);plot(x,real());
title(['a = ',num2str(a)]);
subplot(312);plot(x,imag());
subplot(313);plot(x,abs());
pause(0.7);
end the discrete fractional Cosine
算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001function y = Disfrct(f,a,p)
% Computes discrete fractional cosine transform
% of order a of vector f
% p (optional) is order of approximation, default N/2
% S-C Pei, M-H Yeh, IEEE Tr SP 49(6)2001, pp.1198-1207
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRCT(N,p);
y=E*(exp(-j*pi*a*()).'.*(E'*f));
function E=dFRCT(N,p)
%function E=dFRCT(N,p) returns the NxN eigenvectors of the
%Fourier Cosine transform matrix
global EC_saved pC_saved
if (length(EC_saved) ~= N | pC_saved ~= p),
E = make_EC(N,p);
EC_saved = E; pC_saved = p;
else
E = EC_saved;
end;
function E = make_EC(N,p)
% Returns sorted eigenvectors and eigenvalues of corresponding vectors
% Construct matrix H, use approx order p
N1=2*N-2;
d2 = ; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
d_p = conv(d2,d_p);
st() = d_p; st(1) = 0;
temp = ; temp = temp(:)'./;
s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
H = toeplitz(s(:),s)+diag(real(fft(s)));
% Construct transformation matrix V
V = /sqrt(2);
V = ;
% Compute eigenvectors
Ev = V*H*V';
=eig(Ev);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% = sort(diag(ee));
% ve = ve(:,inde');
E = fliplr(ve);
E(end,:) = E(end,:)/sqrt(2);
[ 本帖最后由 simon21 于 2007-4-26 20:11 编辑 ] the discrete fractional Sine transform
算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001function y = Disfrst(f,a,p)
% Computes discrete fractional sine transform
% of order a of vector f
% p (optional) is order of approximation, default N/2
% S-C Pei, M-H Yeh, IEEE Tr SP 49(6)2001, pp.1198-1207
N = length(f);
shft = rem((0:N-1)+fix(N/2),N)+1;
f = f(:);
if (nargin==2), p = N/2; end;
p = min(max(2,p),N-1);
E = dFRST(N,p);
y=E*(exp(-j*pi*a*()).'.*(E'*f));
function E=dFRST(N,p)
%function E=dFRST(N,p) returns the NxN eigenvectors of the
%Fourier Sine transform matrix
global ES_saved pS_saved
if (length(ES_saved) ~= N | pS_saved ~= p),
E = make_ES(N,p);
ES_saved = E; pS_saved = p;
else
E = ES_saved;
end;
function E = make_ES(N,p)
% Returns sorted eigenvectors and eigenvalues of corresponding vectors
% Construct matrix H, use approx order p
N1=2*N+2;
d2 = ; d_p = 1; s = 0; st = zeros(1,N1);
for k = 1:p/2,
d_p = conv(d2,d_p);
st() = d_p; st(1) = 0;
temp = ; temp = temp(:)'./;
s = s + (-1)^(k-1)*prod(temp)*2*st;
end;
H = toeplitz(s(:),s)+diag(real(fft(s)));
% Construct transformation matrix V
r = N;
V = /sqrt(2);
% Compute eigenvectors
Od = V*H*V';
=eig(Od);
% malab eig returns sorted eigenvalues
% if different routine gives unsorted eigvals, then sort first
% = sort(diag(eo));
% vo = vo(:,inde');
E = flipud(fliplr(vo)); test for Disfrct.m
算法参考：
S.-C. Pei, M.-H. Yeh.
Discrete fractional Cosine and Sine transforms.
IEEE Trans. Sig. Proc., 49(6):1198--1207, 2001
% Disfrft
function testfrct(x)
x=[-35:35];
y = zeros(size(x));
y(fix(length(x)/2)+1)=1;
subplot(321),
plot(x,y);
title('symmetric delta')
yf=Disfrft(y,5/6);
subplot(322),
plot(x,real(yf),'b',x,imag(yf),'r');
axis([-32,32,-0.2,0.2]);
title('FrFT for a = 5/6 of symmetric delta')
% Disfrct
x=;
y1(1:1)=1;
y1(2:36)=zeros(1,35);
subplot(323),
plot(x,y1)
axis()
title('half delta')
yc=Disfrct(y1,5/6);
subplot(324),
plot(x,real(yc),'b',x,imag(yc),'r');
axis();
title('FrCT for a = 5/6 of half delta')
% Disfrcst
x=[-35:0];
y1(1:36)=zeros(1,36);
y1(1:1)=-1;
subplot(325),
plot(x,y1)
axis([-35,0,-1,0])
title('half delta')
ys=Disfrst(y1,5/6);
subplot(326),
plot(x,real(ys),'b',x,imag(ys),'r');
axis([-35,0,-.2,0.2]);
title('FrST for a = 5/6 of half delta') the rescaling preprocessing for the frft routine
算法参考：
H.M. Ozaktas, M.A. Kutay, and G. Bozdagi.
Digital computation of the fractional Fourier transform.
IEEE Trans. Sig. Proc., 44:2141--2150, 1996.
function=rescale(signal,a,delta_x1);
% This routine calculates the parameters to transform the input for the
% fractional Fourier transform when the support of the input
% input is not N^2 with N=sqrt(length(signal))
% Parameters: signal = the signal to be transformed
% a = the order of the transform
% delta_x1 = the length of the support of the input signal.
% To compute the frft with these data use
% Output: signal = possibly flipped signal to avoid infinite factor
% a_new = compute the frft of order a_new
% fact = and elementwise multiply the result with fact
N = length(signal);
delta_x2 = sqrt(N);
k = delta_x2/delta_x1;
a_new =mod(a,4);
if k == 1 || a == 0
a_new = a;
fact = 1;
return
elseif a == 2
signal = flipud(signal)
else
phi = a*pi/2;
psi = acot(cot(phi)/k^2);
c = csc(phi)/(k*csc(psi));
x = linspace(-delta_x1/2,delta_x1/2,N);
u = x;
nu = c*u;
a_new = 2*psi/pi;
A_phi = exp(-i*(pi*sign(sin(phi))/4-phi/2))/sqrt(abs(sin(phi)));
A_psi = exp(-i*(pi*sign(sin(psi))/4-psi/2))/sqrt(abs(sin(psi)));
fact = A_phi/(k*A_psi)*exp(i*pi*(cot(phi)/c^2-cot(psi))*(nu.^2))';
end;
谢谢
:handshake:lol
回复 10楼 的帖子
我看不太懂这个程序，感觉对输入输出参数的理解也很模糊，请清楚的同行解释一下此程序吧。在调用frft.m程序前，需要事先调用此函数吗？ 。。。。。。。。。。好帖，留名。 楼主有没有C的程序啊？急急急啊。 晕死，看不懂，看来对信号这方面俺还是张白纸啊